Comparative genome-wide analysis of extracellular small RNAs from the mucormycosis pathogen Rhizopus delemar

Rhizopus delemar is an emerging fungal pathogen causing devastating mucormycosis in immunocompromised individuals. The organism remains understudied and there are urgent needs for new methods of rapid disease diagnosis for timely therapy. Extracellular vesicles with encapsulated RNAs have recently been discovered to have great potential applications for disease diagnoses and treatments. To explore the utilization of ex-RNA in studies of mucormycosis, we have performed RNA-Seq of ex-sRNAs from two clinical strains of R. delemar. Approximately 3.3 and 3.2 million clean reads were obtained from FGSC-9543 and CDC-8219 strains, respectively. The median sequence length of the sRNAs was 22 nts, with a minimum of 18 and a maximum of 30 nts. Further annotation identified 560 and 526 miRNAs from FGSC-9543 and CDC-8219 strains, respectively. miRNA target prediction and analysis of GO and KEGG pathways have revealed that the regulation of metabolism, secondary metabolite biosynthesis, and two-component system signaling are important during growth. We have also validated RNA-Seq by qRT-PCR and Northern blotting analysis of randomly selected miRNAs. Our results show that R. delemar has a rich reservoir of secreted ex-sRNAs and our studies could facilitate the development of improved diagnostic methods as well as elucidating virulence mechanisms for R. delemar infection.

Mucormycosis is a highly invasive and life-threatening fungal infection that occurs in patients with immune deficiency due to hematologic diseases, diabetic ketoacidosis, or organ transplantation. The infection is caused by zygomycetous fungi, including Rhizopus, Mucor, and Rhizomucor 1 . R. delemar (also known as R. oryzae) is the most common etiologic agent of mucormycosis, responsible for ~70% of clinical cases 2,3 . Because of the increasing prevalence of diabetes mellitus, cancer, and organ transplantation, the number of patients at risk is dramatically increasing. Despite potentially disfiguring surgical debridement and adjunctive antifungal therapy, the overall mortality rate is formidably high, ranging from 50-100% in those individuals with disseminated infection or persistent neutropenia 3 . Significantly, multiple incidences of nosocomial mucormycosis outbreaks have been reported in recent years. During 2007-2009, six fatalities among young children occurred in a pediatric hospital in New Orleans, Louisiana. The deceased victims ranged from premature infants to 13 years 4 . More recently, hospital-acquired mucormycosis outbreaks occurred in the Queen Mary Hospital in Hong Kong and two university hospitals in Pittsburgh, Pennsylvania among organ transplant recipients, resulting in multiple fatalities [5][6][7] . Although individual cases occur more frequently, these nosocomial outbreaks highlight the clinical importance of mucormycosis and the urgency of developing rapid and effective diagnosis methods that will improve treatment strategy and patient outcomes.
The signs and symptoms of mucormycosis are non-specific and definitive diagnosis requires direct examination of specimens and culturing of the organism. These processes are time-consuming and limited in sensitivity 8,9 . The resulting delay in diagnosis and implementation of effective treatment during a critical window of time contributes to the high rate of tragic mucormycosis mortalities and underscores the urgent need for a Results EV extraction and exRNA characterization. Previous studies have shown that there are two types of EVs: exosomes and microvesicles. Exosomes originate from the fusion of the multivesicular membrane and are 40-120 nm in size, whereas microvesicles are formed from outward budding of the surface cell membrane and are larger in sizes (50-1000 nm) (reviewed in 11 ). Both vesicles contain miRNA, many non-coding RNA (ncRNA), and mRNA, as well as cytoplasmic and membrane proteins such as receptors and phospholipids. Even though the numbers of exRNAs identified from C. albicans and C. neoformans were relatively low (2,000 and 1,244, respectively) 13,14 , the results still suggested that R. delemar may also produce EVs with encapsulated exRNAs.
In our ongoing research effort to develop an efficient genetic transformation method for R. delemar, we were also interested in the identification of type III RNA polymerase (Pol III) promoter sequences involved in the synthesis of a large variety of small non-coding RNA (ncRNA). Since a significant amount of ncRNA is encapsulated in EVs, systemic characterization of ex-sRNAs may also lead to the identification of such sequences. For this purpose and for exploring exRNA as a biomarker for mucormycosis, we have isolated membrane vesicles from R. delemar FGSC-9543 and CDC-8219 strains through ultrafiltration and ultracentrifugation. For increased representation, EV preparations from three batches were pooled for exRNA extraction and were found to be of high quality (grade A) suitable for RNA-Seq. cDNA was then synthesized and size fractionated through a PAGE gel. Approximately 100 ng cDNA was subjected to deep sequencing using the BGISEQ-500 desktop sequencer (BGI, Shen Zhen, China). sRNA profiling. The R. delemar FGSC-9543 strain is a clinical strain isolated in 1999 from a brain abscess of a diabetic patient who developed fatal rhinocerebral mucormycosis (www.broadinstitute.org). This is a type culture strain also known as RA99-880 or ATCC MYA-4621. The CDC-8219 strain was recovered from a deceased pediatric patient in a New Orleans hospital in 2008 4 . Initial characterizations of the basic traits and genetic information, including DNA sequences of PyrF, PyrG, (RO3G_14508:EIE89797), CNA1, CNA2, CNA3, and CNB1 (RO3G_16532:EIE91821), suggested that CDC-8219 is genetically similar to FGSC-9543 4 (Bruni and Wang, unpublished observations). Approximately 3.3 and 3.2 million clean reads were obtained from FGSC-9543 and CDC8219 strains, respectively. 86.1% (FGSC-9543 strain) and 69.1% (CDC-8219 strain) of these reads were mapped to the genome by Short Oligonucleotide Alignment Program (SOAP) 15 . The median sequence length of the sRNA was 22 nt (45.9% and 49.2% for FGSC-9543 and CDC-8219 strains, respectively), with a minimum of 18 nt (1.6% and 0.92% for FGSC-9543 and CDC-8219 strains, respectively) and a maximum of 30 nt (0.01% and 0.04% for FGSC-9543 and CDC-8219 strains, respectively) (Fig. 1A,B). The sRNA reads were aligned to Genbank and Rfam database using BLAST 16 for identification and removal of ribosomal RNA (rRNA), transfer RNA (tRNA), small nuclear RNA (snRNA), and small nucleolar RNA (snoRNA) associated reads (Supple Fig. 1A-D). Introns and exons were also mapped and removed (Supple Fig. 1E and F). The distribution of unique ex-sRNAs from both strains is summarized in Fig. 2A,B. miRNA annotation and differential expression. miRNA sequences were identified using the miRBase database. Approximately 393 and 362 known miRNAs were identified from FGSC9543 and CDC8219 strains, respectively (ST1). In addition, the reads that were mapped to the antisense exon, intron, or intergenic region of the genome, but not to any other RNAs, were also analyzed using MIREAP (http://sourceforge.net/projects/ mireap/) to identify 167 and 164 novel miRNAs from FGSC9543 and CDC8219, respectively (ST2).
The FGSC9543 and CDC8219 strains are distinct clinical isolates from different sources. Differential expression of known miRNA between the two strains was compiled and presented in a hierarchical clustering heat map (Fig. 3A), and a scatter plot was also generated using the Expdiff method (p value < 0.05) (Fig. 3B). Most of shared miRNA had similar expression profiles, but the expression appeared to be more depressed in CDC-8219 than in FGSC-9543 (Fig. 3B). Certain miRNA exhibited distinctive expression patterns. For example, miR1028a-3p and miR1043-3p were highly expressed in CDC-8219 but not in FGS-9543, while miR1044-3p, miR1047-3p and miR1048-3p were exclusively expressed in FGS-9543 (ST1). For novel miRNA, TargetFinder predicted 1424 targets for 112 miRNAs, whereas psRobot predicted 1231 target genes for 91 miRNAs (ST5) for FGSC-9543. Again, for CDC-8219, TargetFinder predicted 1529 targets for 116 miRNAs, whereas psRobot predicted 1215 target genes for 88 novel miRNAs, respectively (ST6).

GO terms and KEGG pathway annotation of miRNA targets. GO enrichment and function analysis
was performed using the Web Gene Ontology Annotation Plot (WEGO, wego.genomics.org.cn). GO term annotation revealed 37 GO functional categories for known miRNA targets in both FGSC-9543 and CDC-8219 strains. Cellular and metabolic processes and single organism process (biological process), cell and cell component (cellular component), and binding and catalytic activity (molecular function) processes were the most significantly enriched GO terms (Fig. 5A,B). For novel miRNA, the target genes remained similar and they could be grouped into 28 (FGSC-9543) and 27 (CDC-8219) GO terms (Fig. 5C,D).
The KEGG pathway annotation was performed based on scoring and visualization of the pathways collected in the KEGG database. In the analysis of known miRNA targets, 20 pathways were over-represented in both strains, suggesting that these pathways are significantly regulated during growth of the fungus (Fig. 6A,B). Interestingly, most of the pathways have been shown to be involved in the metabolic pathway and biosynthesis of secondary metabolites. The two-component regulatory system that serves as a basic stimulus-response coupling mechanism to allow organisms to sense and respond to changes in many different environmental conditions was found to be the fourth most significantly regulated. Collectively, these pathway analysis results illustrate the potential roles and mechanisms of miRNAs in regulating the general growth and development of R. delemar.

Validation of expression.
In order to validate the expression profile of exRNAs, we randomly selected six of the known miRNAs and performed quantitative real-time PCR (qRT-PCR) and Northern blot analyses. We extracted total RNA from both strains for cDNA synthesis as we opted to use the expression of the PyrF gene, a putative constitutively expressed gene encoding orotate phophoribosyltransferase, as an internal control. Expression of miR1048-3p in FGSC-9543 and miR1028a-3p in CDC-8219 only was consistent with the data from RNA-Seq. Also consistent was miR1023e-3p whose expression could be found in both strains (Fig. 7A). The expression of the remaining miRNAs (miR1043, miR1044, and miR1047) showed slight inconsistency that could be due to variations in fungal growth or during RNA sample preparations (Fig. 7A). For sample preparation, we used a PEG8000/NaCl method to first separate small RNA from total RNA prior to denaturing polyacrylamide gel electrophoresis separation. Similar to the qRT-PCR results, the expression of five miRNAs detected by Northern blot analysis was largely consistent. The only exception was that the expression of 1048-3p was detected in CDC-8219 (Fig. 7B). However, the expression is very low (0.06 in relevance to the expression of the pyrF gene control) in comparison to FGSC-9543 (0.82) (Fig. 7B). Taken together, these findings are mostly agreeable with the data obtained from RNA-Seq.

Discussion
Zygomycetous fungi such as Rhizopus, Mucor, and Rhizomucor are emerging as significant fungal pathogens. Multiple outbreaks of nosocomial mucormycosis, caused by R. delemar, have occurred in recent years, thereby underscoring the importance of improving diagnostic methods as well as furthering our understanding of how these fungi cause disease. However, this group of fungi is poorly studied and therefore lacks many available genetic research tools and resources that have been developed for studying other medically important fungi, including Candida, C. neoformans, and A. fumigatus. In order to improve molecular resources for studying In addition to sRNAs, we have also identified a combined total of 17,884 extracellular mRNAs and extracellular 1,932 lncRNAs from these two strains. We present herein the comparative analysis of secreted sRNAs from two clinical strains of R. delemar. Since this appears to be the first study in which RNA profiling of ex-sRNAs from a zygomycete was extensively performed, it is unknown if these molecular sequences are unique to R. delemar or whether some exRNA molecular signatures might be similarly conserved with other zyogmycetes such as Mucor circinelloides. Further studies are needed to ascertain the uniqueness of R. delemar exRNA signatures among other zygomycetes; however, it should be noted that the optimal treatment of infections caused by each of the zyogmycetes is the same.
Previous studies by others have revealed that exRNA is produced in relatively high abundance in response to environmental stress and during infection of the host. Circulating pathogen-derived exRNAs are readily detectable in bodily fluids, including blood, bronchoalveolar lavage, and urine [17][18][19][20][21] . These circulating, pathogen-derived exRNAs can be easily isolated from an infected host in a non-invasive manner and hold great potential as molecular biomarkers of infection. The development of a sensitive and rapid detection method is especially important for diagnosis of mucormycosis, since the symptoms are often non-specific. Furthermore, traditional methods of identification, either by direct examination or culture of the organism, may delay diagnosis and implementation of effective treatment during a critical window of time. This delay in diagnosis and administration of the appropriate antifungal therapy contributes to the tragically high mortality rate of 50-100% 3 . Therefore, in order to explore the potential application of ex-sRNA as a molecular biomarker of mucormycosis, we sought to profile ex-sRNAs in this preliminary study. We included two clinical strains, the type culture FGSC-9543/RA 99-880 isolated from a brain abscess of a diabetic patient and the CDC-8219 strain isolated from a pediatric patient in New Orleans with fatal mucormycosis. Preliminary characterization of the basic physiological traits revealed that these strains are very similar in growth and morphology (unpublished observations), indicating that our selection of strains is representative of the species. In addition, we were able to obtain several million clean RNA-Seq reads from each strain that also suggests that the reads far exceed the number needed to saturate the transcriptome(s) of the strain grown in monocultures.
Small RNA, including miRNAs, small interfering RNAs (siRNAs) and piwi-interacting RNAs (piRNAs) play a critical role in the regulation of cellular growth and development by regulating the expression of target genes by RNA cleavage or transcriptional silencing 22 . miRNAs have been reported as vital participants in post-transcriptional gene regulation by affecting mRNA stability and translation efficiency in organisms, including humans, fruit flies and nematodes (reviewed in 23,24 ). miRNAs regulate gene expression through targeting the 3′ untranslated region (UTR) of target mRNAs in a seemingly complex pattern, as a single miRNA can regulate hundreds of target genes. Conversely, a single gene can be targeted simultaneously by multiple miRNAs. In contrast to other model organisms, the functions of sRNAs, especially ex-sRNAs, have been relatively unexplored for zygomycetous fungi such as Rhizopus and Mucor spp until now. Here, our studies revealed that the FGSC-9543 strain has 393 known miRNAs (373 targeting 4,810 host genes) and the CDC-8219 strain has 362 miRNAs (338 targeting 4,082 host genes) in the secreted sRNA profiles, respectively. In addition, along with known miRNA and targets, 116 novel miRNAs with 1,848 targets and 118 novel miRNAs with 1,936 targets were also predicted for FGSC-9543 and CDC-8219 strains, respectively. While the miRNA profiles of each strain are very similar, strain-specific miRNAs were also present. Interestingly, statistically significant differential expression patterns were observed for more than 50% of the known miRNAs produced by the two strains. Furthermore, we validated the expression of six known miRNAs by qRT-PCR and five miRNAs by Northern blot analysis and found that the expression levels are very similar to the results obtained from RNA-Seq. Further exploration of miRNA expression and characterization of target genes will expand our understanding of the role individual miRNAs play in the growth, development, and even virulence of R. delemar, and are, thereby, warranted.
There are three ontologies in Gene Ontology (GO): molecular function, cellular component and biological process. Through GO and KEGG pathway analysis of putative targets of miRNA, we found that most of the biological functions are related to carbohydrate utilization and metabolism, and secondary metabolite biosynthesis, pathways that are critical for fungal growth. Interestingly, the two-component system represented one of the major pathways targeted by miRNA in our analysis. This is not surprising since sensing and responding to changes in the environment is central to fungal growth. Based on the GO and KEGG analysis results, it is highly intriguing to consider whether exRNA profile studies involving clinical or in vivo samples would reveal pathways that are important during infection, but not expressed in vitro.
Future studies will be aimed at profiling additional exRNAs including mRNAs and lncRNAs that could uncover novel molecular mechanisms governing growth, virulence mechanisms, and the response to antifungal   4,000 × g for 15 min at 4 °C and then at 15,000 × g for 15 min. The latter step was necessary for removing smaller debris or particles. The liquid was filtered through an ultrafiltration filter with a cutoff of 100-kDa (Millipore). The remaining liquid was precipitated by ultracentrifugation (Beckman) at 100, 000 × g for 1 hour and at 4 °C, and the precipitate containing EVs was re-suspended in phosphate-buffered saline (PBS). The vesicle pellets were washed twice in PBS and the final pellets were suspended in PBS and lyophilized for RNA extraction. EV preparations from three repeats were pooled.
Cloning of ex-sRNAs. RNA was extracted with Trizol (Sigma-Aldrich) and separated on a PAGE gel. The small RNA band at approximately 18-30 nt in size were recovered and cloned with the Small RNA Cloning Kit (TaKaRa). Briefly, sRNA was detected by an Agilent 2100 Bioanalyzer (Agilent Technologies) and 8 μg purified sRNA was ligated with the 3′ and 5′ RNA adaptors included in the kit. cDNA was subsequently synthesized and amplified with the supplemented oligo primers. cDNA was cloned into T-vector pMD20 provided with the kit to generate a sRNA library. The cDNA libraries were loaded onto the flow cell channels of an BGISEQ-500 platform for paired-end 90 bp x 2 sequencing at the Beijing Genomics Institute (BGI), Shenzhen, China.
Identification of miRNAs via bioinformatics analysis. The ligating adapters with unique synthetic oligonucleotide sequences (sequence tags) were subjected to data cleaning analysis to obtain clean tags/reads. The clean tags were annotated into different categories and tRNA, rRNA, and other impurities were removed. The sequences were mapped to the R. delemar genome at ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/149/305/ GCA_000149305.1_RO3/GCA_ 000149305.1_RO3_genomic.fna.gz. Known miRNAs were identified by searching against miRBase (www.mirbase.org). Since the miRNA hairpins are mostly located in intergenic regions, introns or reverse repeat sequence of coding sequence, sRNAs that were mapped to these regions, but not to any other RNA, were predicted to be novel miRNAs. The Rfam database is at http://rfam.janella.org. Expression levels of sRNA were calculated by using Transcript Per Million (TPM) 25 .

miRNA target identification and functional annotation by GO and KEGG analyses.
Once miRNA results were obtained, target prediction was performed using psRobot and TargetFinder software 26,27 . Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were also performed. GO provides a common descriptive framework and functional annotation and classification for analyzing the gene set data (www.geneontology.org), whereas the KEGG pathway database is a recognized and comprehensive database including all kinds of biochemical pathways (www.annotation.jp/KEGG/ko.html) 28,29 .
Quantitative RT-PCR. For quantitative real time RT-PCR (qRT-PCR), 5 μg of total RNA was reverse transcribed into first-strand cDNA using the oligo(dT) or random primers and M-MLV Reverse Transcriptase (Invitrogen). The qRT-PCR reactions were performed following previously established procedures 30,31 . Primer pairs used in this section are listed in Table I. The expression of the PyrF gene that encodes orotate phophoribosyltransferase, a putative constitutively expressed gene, was used as a control. DNA primers are listed in Supple Table 7.
Northern blot Analysis of miRNAs. Northern blot analysis of miRNA was performed following the protocols described previously 32,33 . Briefly, small RNAs were isolated from total RNA by precipitation in 5% PEG8000 and 0.5 M NaCl. Small RNAs were resolved using 14% denaturing polyacrylamide gel electrophoresis, transferred to the nylon membrane (GE healthcare life sciences) using a semi-dry transfer apparatus (Trans-blot semi-dry transfer cell, Bio-Rad), and detected by a non-radioactive biotin-labeling method (Chemiluminescent nuclei acid detection module kit, NEB). Both miRNA and PyrF probes were labeled with biotin at both 5′ terminus and 3′ terminus. Quantitation of the blots was performed using ImageJ and the relative signals were calculated by normalizing to the PyrF control band.