Genome-wide polyadenylation site mapping datasets in the rice blast fungus Magnaporthe oryzae

Polyadenylation plays an important role in gene regulation, thus affecting a wide variety of biological processes. In the rice blast fungus Magnaporthe oryzae the cleavage factor I protein Rpb35 is required for pre-mRNA polyadenylation and fungal virulence. Here we present the bioinformatic approach and output data related to a global survey of polyadenylation site usage in M. oryzae wild-type and Δrbp35 strains under a variety of nutrient conditions, some of which simulate the conditions experienced by the fungus during part of its infection cycle.

Polyadenylation plays an important role in gene regulation, thus affecting a wide variety of biological processes. In the rice blast fungus Magnaporthe oryzae the cleavage factor I protein Rpb35 is required for pre-mRNA polyadenylation and fungal virulence. Here we present the bioinformatic approach and output data related to a global survey of polyadenylation site usage in M. oryzae wild-type and Δrbp35 strains under a variety of nutrient conditions, some of which simulate the conditions experienced by the fungus during part of its infection cycle.

Background & Summary
Magnaporthe oryzae is an economically relevant fungal pathogen of cereals, and represents an important threat to food security, particularly for cereal-dependent populations 1 . It is responsible for the devastating blast disease in rice and other grasses such as wheat, finger millet and maize 1 . The RNA-binding protein Rbp35 has been characterized in this fungus 2,3 . No clear orthologues of this protein are present in budding yeast. The Δrbp35 mutant shows reduced virulence on leaves and roots of rice. Interestingly, this protein interacts with CFI25, the homologue of human CFIm25 4,5 , indicating that Rbp35 is a component of the fungal polyadenylation machinery. In fact, Rbp35 is required for the alternative polyadenylation (APA) of pathogenicity-associated genes 3 . The pre-mRNA 3 0 end processing occurs in two steps. First, the pre-mRNA is cleaved and, subsequently, a polyadenosine (poly(A)) tail is added 6 . Unfortunately, standard high throughput RNA sequencing does not allow the identification of polyadenylation sites (PASs) since only a small percentage of sequencing reads contain poly(A) tails. To circumvent this problem, several techniques have been developed to enrich the transcript ends of poly(A)+RNAs (reviewed in Marconi et al. 7 ).
Here, we present the datasets derived from a genome-wide PASs profiling within the M. oryzae wild type isolate Guy11 and the Δrbp35 mutant using the 3 0 T-fill method 8 . The datasets include three biological replicates in four nutritional conditions including nitrogen and carbon starvation. These data will help to identify the coordinated RNA networks that regulate M. oryzae growth in starved cells as well as gene networks regulated by Rbp35. These data also provide useful information for enhancing genome annotations and for cross-species comparisons of PASs and PAS usage within the fungal kingdom. The biological analyses of these data have been published 9 .

Methods
Fungal media, growth conditions and infection assays M.oryzae strain Guy11 (WT) and 2D4 (Δrbp35) were grown on CM (complete medium), MM (minimal medium) and DCM (defined complex medium) 10 . Starvation conditions were carried out on MM depleted of nitrogen or carbon sources.

RNA extraction, library preparation and RNA sequencing
Fungal material was harvested from 60 hour-old fungal mycelia grown in liquid CM (24 h) and then transferred for an additional 12 h in a fresh CM, MM, MM-N or MM-C. Three biological replicates were Figure 1. Workflow describing the steps followed to analyze 3 0 T-fill sequencing data. After total RNA extraction, Illumina RNA libraries were created following the specifics of the 3 0 T-Fill protocol 8 . RNA sequencing was then performed on the Illumina HiSeq2000 platform. FastQ reads were processed by removing sequencing adapters and finally aligned onto the reference genome. Alignments were further processed, filtering out low quality mappings, mappings with high A/T content and potential internal priming. Reads were thereafter assigned to overlapping genes up to 400 bp from the annotated gene end. Finally, identified poly(A) sites were clustered together in order to select high-confidence poly(A) sites. independently harvested and extracted for each library preparation. Total fungal RNA was isolated by a standard trizol method 11 and treated with RNase-free DNaseI using Turbo DNA-free kit (Ambion). Library preparation was carried out as described by the EMBL (Heidelberg) sequencing service, and similarly, the sequencing libraries and the 3'T-fill reaction were prepared as previously described 8 . The final libraries were loaded into the cluster station (cBot, Illumina) and the priming buffer was exchanged for the T-fill solution: 101 μl water, 20 μl Phusion buffer HF (5 ×) (NEB), 3 μl dTTPs (10 mM), 0.8 μl genomic DNA Sequencing primer V2 (100 μM) (Illumina), 3 μl non-hot start Phusion polymerase (2 U/μl, NEB) and 2 μl Taq polymerase E (5 U/μl; Genaxxon). The addition of the latter Taq polymerase was crucial for a complete T-fill. After clustering, the samples were sequenced on a HiSeq 2000 (Illumina).  alignments obtained were later filtered removing low-quality mapping (MAPQ o 30, see https:// samtools.github.io/hts-specs/SAMv1.pdf), and mappings with high A/T content ( >80%) and potential internal primings (more than seven consecutive As in the first 5' nucleotides or more than twelve As in the the first sixteen nucleotides, which is the length of the oligodT primers used during the 3 0 T-fill protocol); these alignments were discarded because they are suspected to be potential internal primings with genomic regions rich in As. In fact, a closer analysis showed that they displayed an unusual nucleotide profile ( Fig. 3a and b), confirming that they were not a product of mRNA 3'UTR sequencing. Globally, approximately 2.5% of the alignments were discarded according to the previous reasoning. Alignments that were not discarded during the previous process were further analyzed, and found to mainly align to annotated 3 0 UTRs regions of the genome (Fig. 3c), supporting the idea that these were in fact genuine sequencing from the 3 0 end of mRNA transcripts. On the other hand, the discarded alignments usually derived from non-annotated genomic regions, possibly related with transposable elements or other highly variable nucleotide sequences.

RNA sequencing and PAS detection protocol
Taking all wild-type data sets, 14,593 PASs were reliably assigned to the genomic features listed in the Ensembl Fungi gene annotation for M. oryzae version 27 (http://fungi.ensembl.org/), which contains 13,218 annotated features. Each read was assigned to its overlapping feature (including protein-coding genes and the available ncRNAs), ambiguous cases were assigned to the closest 3' terminal end. To account for inaccurate (usually short) 3'UTR annotations, which are common even in well-annotated genomes 14 , we extended the gene annotation by 400 bp. The selection of this supplemental UTR range was justified by the fact that on average reads were at most 400 bp far from the closest annotated stop codon (Fig. 3d). Reads that could not be assigned to any known feature are not included in this dataset. A PAS was considered to be called with high-confidence if it was detected in at least two of the three replicates according to the following rule: its expression is considered distinct from basal noise if its standard score calculated against the whole gene expression (summing up all PASs expression from that gene) has a confidence level >99% (z-score > 2.58 or z-score o 2.58). The standard score was calculated as: z-score = (value-mean)/std, where "value" is the number of supporting reads for the PAS, "mean" is the mean value of supporting reads for all the PASs in the gene, and "std" is the standard deviation of supporting reads for all the PASs in the gene. PASs with less than five supporting reads in every replicate were discarded. High-confidence PASs that fell within 33 bp of one another (roughly the span of the PAS regulatory region, including the A-rich and the U-rich regions upstream of the cleavage site) were resolved to a single PAS -that being, the one that was most highly expressed in that 33 bp region.

Code availability
A Galaxy workflow describing the procedure followed to reproduce the PAS identification is available at (Magnaporthe oryzae polyadenylation sites for wild-type and Δrbp35 mutant, Data Citation 1), together with additional scripts required to run the workflow together with a summary of the pipeline.

Data Records
All sequencing data have been uploaded to the National Center for Biotechnology Information (NCBI) (Alternative polyadenylation controls pathogenicity-associated mechanisms in the rice blast fungus, Data Citation 2), with an overview of the submission provided in Table 1. These data contain paired-end read sequencing for the 24 samples considered: wild-type and Δrbp35 mutant strains in four different growth conditions (complete medium (CM), minimal medium (MM), nitrogen starvation (MM-N) and carbon starvation (MM-C)), in three replicates. Each archive is stored in FastQ file format.
Illumina deep sequencing read alignments are provided as BAM files (samtools.github.io/hts-specs/ SAMv1.pdf), which are sorted and come with the relative index files (Magnaporthe oryzae polyadenylation sites for wild-type and Δrbp35 mutant, Data Citation 1); Bed Graph files were created in order to ease the visualization of PAS-usage levels and are available in the same deposit, together with Generic Feature Format (GFF) files that represent the high-confidence PASs called by our algorithm. All

Technical Validation
Illumina deep sequencing produced 24 samples providing between between 4,751,592 and 11,517,077 reads per sample (Table 1). FastQC analysis indicated that these reads were generally of high quality (Fig. 2). The correlation between biological replicates, based on calculated gene expression (i.e. number of reads observed for each gene), ranged from 92% to 97%. Correlation matrix analysis revealed a strong similarity between the biological replicates (Fig. 4). Approximately 73% of reads could be aligned to the M. oryzae genome, and about 36% were finally assigned to a high-confidence PAS (as described in methods). Paired-end reads were 44 bp long with a mean mapping quality of 36. Only 1% and 4% of reads contained illumina adapters or poly(A) tail residues, respectively.
Internal priming was detected and filtered out from the sequencing data using the algorithm described in the methods section, where about 2.5% of the total reads were found to result from internal priming.
As described in methods, reads were assigned to genomic features after an additional 400 bp was included beyond the annotated 3'UTR to compensate for incorrect genome annotations. We justify this by noting that mapped reads are commonly found to accumulate as far as 400 bp from protein-coding genes' stop codons. The typical upstream A-rich region which contains the polyadenylation signal, was found within 3 bp from the cleavage site. Consequently, high-confidence PAS detection included the clustering of PASs closer than 33 bp.

Usage Notes
The publicly accessible triplestore (https://dydra.com/markw/polyadenylation-sites-in-magnaportheoryza/) provides a sample SPARQL query representing the query we anticipate to be most informative (i.e. for a given locus, what PASs are used under which conditions). A Web interface has been constructed to enable visualizing the same query results via a Web browser (http://linkeddata.systems/Magnaporthe/ polyA_Sites/). Consistent with the FAIR Principles, machine-readable metadata about this SPARQL endpoint is available in VOiD format by resolving the URL of the endpoint, and machine-readable metadata about the data within the endpoint is available by resolving the named graph URI of the dataset. A track hub was also registered, compatible with genome browsers such as Ensembl 17 , under the name ID "Magnaporthe oryzae poly(A) sites" (Magnaporthe oryzae poly(A) sites, Data Citation 3).