Discovery of Conserved and Novel MicroRNAs in Galdieria sulphuraria

1. Background

Red algae represent a large variety of primitive algae that occupy an important position in plant systemic evolution ( 1 ). Throughout the long evolutionary process, a few species have acquired unique biological functions. Galdieria sulphuraria, a hot spring micro-red algal species ( 2 ), often thrive in springs, soils, or volcanic rocks and acidic and thermal environments. This extremophile can perform nearly all biological activities in extreme environments and can tolerate various biotic and abiotic stresses ( 3 ). Based on these characteristics, we speculate that, in addition to being a valuable microalga,  G. sulphuraria is an essential representative organism for exploring the extremophile’s unique biological functions.

Since the establishment of the formal nomenclature and classification of G. sulphuraria ( 3 ), some works have been undertaken to assess this microalga such as investigations of its culture conditions and metabolic regulation ( 4 ), biological functions ( 5 ), and molecular evolution ( 6 ). However, to date, no studies have assessed its noncoding RNA. MiRNAs of approximately 20～24 nt are a class of noncoding sRNAs with essential roles in gene posttranscriptional regulation. Increasing studies have shown that miRNAs play important roles in an organism’ growth and development, metabolic regulation, and response to various stresses ( 7 ). Gene expression level is downregulated by cleavage or translational inhibition mediated by miRNA. Although Cyanidioschyzon merolae  (another extremophilic, hot spring species of the family Cyanidiaceae) miRNAs have been predicted by bioinformatics ( 8 ), miRNAs in G. sulphuraria have not been detected thus far. Thus, identification of the miRNAs in this unique microalga will be very important.

Next-generation sequencing is a fast and reliable method for the identification of unknown miRNAs in different species ( 9 ). Various bioinformatics tools can provide effective supports for the detection of potential miRNAs and their targets in non-model algae. Thus, sRNA-seq combined with bioinformatics will be an effective strategy for identifying miRNAs in G. sulphuraria. Furthermore, the prediction of putative gsu-miRNA targets, as well as their gene function and metabolic pathways, will also be helpful for research on some critical biological processes of this microalga.

2. Objectives

G. sulphuraria occupied an essential position in the red algae evolutionary process has some special biological functions. MiRNA has important functions in plant growth, development, and response to biotic and abiotic stresses. Conserved and novel miRNAs identified in G. sulphuraria will lay a foundation for studying their regulating functions and provide insights into further exploitation of this thermoacidophilic alga.  

3. Materials and Methods

3.1. Algal Cultivation, sRNA Library Construction, and Sequencing

An identified G. sulphuraria strain (No. UTEX2919) was cultured in static culture vessels that contained a modified saline-adenine-glucose medium with a pH value of 2.7 at room temperature. The culture of G. sulphuraria was illuminated with a 500 µmol.m.^-2s^-1 photosynthetic photon flux density. In its maturity stage (60～70 days), the liquid alga with the most living cells was divided into 3 equal groups (the average density of 573×103±14 cells.mL^-1) (Fig. S1). Then each group was enriched and immediately frozen in liquid nitrogen on the sixtieth day. 

The Biozol mixed reagents were added into the quick-frozen 3 algal samples and mixed for 10min. Then chloroform was added into them and centrifuged for 15min at 4 ℃. The isopropanol was added into the extracted upper layers and mixing. The mixtures were incubated at -20 ℃ for 20min and then centrifuged at 4℃ for 10min. Washing and dissolving the final RNA precipitations with 75% ethanol and RNA-free water, respectively. Following previously described methods ( 10 ), 3 high-quality sRNA TruSeq libraries were then constructed based on the reverse transcription of enriched sRNA pools; approximately 12 µg of sRNA for each group was isolated using an Illumina TruSeq Small RNA Sample Prep Kit (Illumina, USA). Quality assessment of the libraries was performed via an Agilent 2100 Bioanalyzer (Agilent, USA) and a Step OnePlus Real-Time PCR System (ABI, USA). Finally, deep sequencing was performed via an Illumina HiSeq 4000 by BGI (Shenzhen, China). Detail flowchart of this study was shown in Figure S2.  

3.2. Bioinformatics Analysis of sRNA-Seq Data

The common sRNA reads from three experimental groups were screened as a data set via filtering the redundant data. After removing the reads with 3’ adapter null, insert null, 5’ adapter contaminants, length ≤18nt, and poly-A, the remaining high-quality common reads were retained and then annotated according to the procedure of Axtell and Meyers (2018) ( 11 ). In addition to known miRNAs, the reads identified by sequences alignment analysis in the reported plant miRNAs in miRBase 22.0 ( 12 ), other sRNAs, including rRNAs, tRNAs, snRNAs, and snoRNAs, annotated in Rfam 11.0 and GenBank, and the remaining unannotated sequences with >98% accuracy were retained for prediction of novel miRNAs.

3.3. Identification of miRNAs

In addition to the known miRNAs above, to obtain as many miRNAs as possible in G. sulphuraria, alignment analysis of the ESTs and GSSs of Galdieria against all plant miRNAs and their precursors in miRBase 22.0 was conducted. The above-aligned sequences were screened as the potential candidate conserved gsu-miRNAs. The candidate conserved miRNAs and unannotated sRNAs were further screened and predicted as the additional and novel ones according to the criteria: (i) the sum of mismatches and gaps could be no more than 4 for both the mature miRNA and its precursor; (ii) the length range of the mature miRNA could be 19～25 nt; (iii) the candidate RNA could not be identified as any a non-miRNA; (iv) miRNA/star miRNA duplex contains less than 2 mismatched bases and no large loops; (v) the precursor could form a perfect stem-loop secondary structure via running Mireap with an MFEI value more than 0.85. Furthermore, the potential MIRNA families and their members in G. sulphuraria  were predicted based on the prediction criteria reported by Meyers et al (2008) ( 13 ) by analyzing the alignment results of gsu-miRNAs and their precursors against all plant miRNAs in miRBase 22.0. 

3.4. Conservation Analysis of miRNAs

All miRNAs and fifteen random sRNAs with the size of 20～24 nt ( 11 ) in G. sulphuraria were selected to research their conservation against all known miRNAs in 20 representative plant species stored in miRBase 22.0 via BLASTn search. Sequence alignment analysis was conducted based on the BLASTn-short model with sequences identities percentage=100%, gaps=0, and E-value≤1e-05. The sequence conservation could be evaluated by constructing the heat map using R language with different alignment quality bit-score. Moreover, the maximum likelihood tree among 20 plant species was also constructed by using the Muscle, Gblock, and MEGA 5.0 based on their 18S rRNA sequences downloaded from the SILVA rRNA database project (www.arb-silva.de). 

3.5. Northern Blot and RT-PCR Validation

Four relatively high-expressed miRNAs, miR5813 (396 reads per miRNA), miR8039 (382 reads per miRNA), gsu-miR1a-5p (100 reads per miRNA), and gsu-miR5-5p (67 reads per miRNA) in G. sulphuraria were randomly selected for northern blot validation. Following previous methods ( 14 ), Northern blot detection was conducted based on hybridization with miRNA-complementary DNA oligonucleotides labeled with digoxigenin (Roche, Switzerland). Furthermore, following previously described methods ( 15 ), 20 conserved, and 13 novel gsu-miRNAs were randomly selected for RT-PCR validation. After synthesis of cDNA via a HiScript 1st Strand cDNA Synthesis Kit (Vazyme Biotech, USA), RT-PCR was carried out with TransScript Two-Step RT-PCR SuperMix (Transgen Biotech, China). 

3.6. Prediction of miRNA Targets

Targets of conserved and novel gsu-miRNAs were predicted by running psRobot and TargetFinder. As query sequences, gsu-miRNAs were aligned against transcripts and ESTs of hot spring red algal sequences deposited in GenBank. In plants, based on the principle of perfect complementary between miRNAs and their targeted mRNAs, putative targets could be predicted according to the procedure of Xie et al (2007) ( 16 ) .Additionally, the potential miRNA inhibition type could be predicted via running psRNATarget.

3.7. GO, KEGG Pathway, and Cytoscape Network Analyses

The putative genes targeted by gsu-miRNAs were subjected to GO and KEGG analysis. GO terms with a corrected p-value of ≤0.05 are defined as significantly enriched. After which, the directed acyclic graphs (DAGs) were constructed with the Goseq and Topseq packages. Similar to the GO analysis, the enriched KEGG pathways of target genes were predicted with enrichment ratio of all genes number in some pathway to Q-value. Their pathway maps were obtained using KegSketch. The Cytoscape network between gsu-miRNAs and their targeted genes was constructed via running Cytoscape 3.0. 

4. Results

4.1. Characterization of Gsu-sRNAs

The high-quality sRNA TruSeq libraries of G. sulphuraria were constructed (RINs and 28S/18S values were 7.5 and 1.9 on average, respectively) in a laboratory-based environment. Approximately 12.20 M raw sRNA-seq reads were obtained with 10.47 M clean reads (85.80%) retained (Table 1). Table1 Summary of sequencing reads from a G. sulphuraria sRNA library. Type Number of reads Percent (%) Total reads 12,238,136 High quality 12,200,966 100% 3’ adapter null 271,528 2.23% Insert null 557,729 4.57% 5’ adapter contaminants 131,145 1.07% Smaller than 18 nt 771,596 6.32% PolyA 325 0.00% Clean total reads 10,468,643 85.80% Mapping to genome 2,419,809 23.11% Approximately 2.42 M (23.11%) total sRNA reads could be mapped to the reference genome (Table 1). A majority of gsu-sRNA tags were detected between No. 315 and 420 scaffolds of G. sulphuraria (Fig. S3A). The majority of gsu-sRNA reads were distributed between 18 and 24 nt (Fig. S3B). Lengths of 19 and 20 nt were the most frequent sRNA sizes in G. sulphuraria. Approximately 5.05 M total gsu-sRNAs (48.27%) and 0.88 M unique gsu-sRNAs (72.73%) could not be annotated in this study (Fig. S3C, D). 

4.2. Characterization of Gsu-miRNAs

Totals of 120 known or conserved miRNAs and 14 novel miRNAs were identified in G. sulphuraria (Table S1A-C). Moreover, 111 known MIRNA families with 120 members were predicted in G. sulphuraria (Table S1D). Apart from 9 families containing two members, a majority of the MIRNA families, accounting for 91.94% overall, only has one. The nucleotide bias at each position revealed that A and U were the main bases in G. sulphuraria (Fig. 1A). The most biased base at the first position of 5’ nucleotides in novel miRNAs was U (Fig. 1B). The length distribution of conserved and novel gsu-miRNAs ranged from 18 to 25 nt, with most consisting of 20～22 nt (accounting for 62.69%) (Fig. 1C). Partial precursors secondary structures of gsu-miRNAs are presented in Figure S4.

Figure 1. Nucleotide and size characterization of gsu-miRNAs. (A) Nucleotide bias of conserved gsu-miRNAs at each position. (B) Frequency of first-nucleotide bias of novel gsu-miRNAs. (C) Size distribution of gsu-miRNAs.

4.3. Conservative Analysis of miRNAs

Twenty-one out of 149 miRNAs including sRNAs identified in G. sulphuraria were highly conservative against the 6,741 known plant miRNAs with higher sequence similarity evaluation scores (≥32.2) and at least one aligned miRNA sequence (Fig. 2). Except for 3 lowly conservative miRNAs with score values <11.0 (miR164c, miR7-3p, and miR6-5p), we found another two lowly conserved non-miRNAs (snoRNAt-0219861 and snRNA-t745092) in G. sulphuraria across 20 plant species (Table S2A). From the horizontal view of Figure 3, 76 subject sequences can be detected among 18 species based on the query sequence gsu-miR319a. From the vertical view, 36 conserved miRNAs could be detected in Glycine max based on 7 query gsu-miRNAs. Results of average sequences identity percentage demonstrated that miR159a and miR408-3p in G. sulphuraria have the highest similarity (85%) among these species. 

Figure 2. Conservation analysis heat-map between gsu-sRNAs and miRNAs in 20 representative plants with different phylogenetic relationships. The conservation of gsu-miRNAs or other sRNAs in at least one species in a given evolution order is shown using sequences alignment quality bit-score (corresponding miRNAs count) with different shades of blue cells.

Figure 3. Validation of some gsu-miRNAs. (A) RNA gel blot hybridization of digoxigenin-labeled probes for four gsu-miRNAs. (B) Agarose gel of RT-PCR products based on 33 gsu-miRNAs with positive (U6) and negative controls (water).

4.4. Validation of Gsu-miRNAs

Four random miRNAs (two conserved and two novel miRNAs) in G. sulphuraria were subjected to northern blot hybridization with designed oligonucleotide probes (Table S2B). All four mature gsu-miRNAs were successfully detected via triplicate tests (Fig. 3A). RT-PCR with specific stem-loop primer pairs also represents another effective miRNA validation method and was used here. All PCR products of 33 miRNAs in G. sulphuraria were positive based on triplicate tests (Fig. 3B). The RT-PCR primers information is provided in Table S2C. 

4.5. Putative Targets

Altogether, non-redundant 1,505 and 84 putative targets were predicted for the conserved and novel gsu-miRNAs, respectively (Fig. S5). An average of 15 targets per gsu-miRNA could be predicted in this study. Additionally, psRNA Target prediction revealed that cleavage was the main posttranscriptional inhibitory type of conserved and novel gsu-miRNAs with 59.46% and 61.11% inhibition proportion, respectively (Table S3A-D). 

4.6. GO and KEGG Pathway Enrichment

The order of GO classification was BP>CC>MF for conserved gsu-miRNAs (CC>BP>MF for novel gsu-miRNAs) (Fig. 4A, B). The DAGs for the top 5 enriched GOs are illustrated in Figure S6. The three most enriched GOs for conserved gsu-miRNAs in the BP, CC, and MF were related to organic substance metabolic processes, chloroplast thylakoids, and hydrolase activity, respectively (Table S4A-D). Similarly, electron transport chain, chloroplast part and electron carrier activity were the three most enriched GOs for novel gsu-miRNAs in the BP, CC, and MF, respectively (Table S4E-H). 

Top 20 enriched KEGG pathways of gsu-miRNAs were identified (Fig. 4C, D). RNA transport (ko03013) (Table S5A) was considered the most significantly enriched pathway for conserved gsu-miRNAs (Fig. 4C, Fig. S7A). Similarly, the pentose phosphate pathway (PPP) (Table S5B) was considered the most significantly enriched pathway for novel gsu-miRNAs (Fig. 4D, Fig. S7B).

Figure 4. GO and KEGG analysis. (A) GO classification of genes targeted by conserved gsu-miRNAs. (B) GO classification of genes targeted by novel gsu-miRNAs. (C) Top 20 enriched KEGG pathways of genes targeted by conserved gsu-miRNAs. (D) Top 20 enriched KEGG pathways of genes targeted by novel gsu-miRNAs.

4.7. Cytoscape Networks Construction

Cytoscape networks between conserved gsu-miRNAs and their target genes are indicated by purple circles and black lines (Fig. S8A); the most enriched network contained 79 conserved gsu-miRNAs with 1,485 genes. Similarly, the networks between novel gsu-miRNAs and their target genes were also indicated by yellow circles and black lines (Fig. S8B); the most enriched network contained a novel gsu-miR9 that had 23 putative genes.

5. Discussion

Lots of conserved and novel miRNAs in G. sulphuraria were identified for the first time, which is more than another unique microalga Collomyxa suellipsaidea (118 conserved and six novel miRNAs) ( 15 ). Even so, merely 23.11% gsu-sRNA reads were mapped onto the reference genome with most of them located at No. 315～420 scaffolds. Excluding the sequencing accuracy factor, we speculate that low-quality reference genome data might lead to this ( 11 ). Consistent with the previous research ( 11 ), inappropriate miRNA annotations in miRBase especially for non-model algae might limit more gsu-miRNAs identification (72.73% unannotated gsu-sRNAs ). 

A few highly conserved gsu-miRNAs with a lower percentage (14.09%) than the land plants miRNAs (～30.71%) ( 11 ) are identified against the representative known plant miRNAs in miRBase. We speculate that the miRNAs in G. sulphuraria might be subject to a relatively fast evolutionary rate, which has been confirmed in miRNAs of another microalga Chlamydomonas reinhardtii ( 9 ). As miRNAs evolutionary research on the high plant ( 17 ), further research on algae miRNAs would provide more important information for lower organisms miRNA evolutionary process.

The positive validation rate of gsu-miRNA in this study via Northern blot and RT-PCR was 100%, but the former still be more reliable than the latter. Northern blot analysis has been successfully applied to the construction of a synthetic gene circuit mediated by miRNA in C. reinhardtii  ( 12 ). RT-PCR is also a rapid miRNA validation tool especially for a large number of miRNA samples though their temporal-spatial expression differences may lead to a few negative results. 

Cleavage is the main miRNA inhibition type in a plant ( 11 ), which could be confirmed based on our prediction results (60.28% cleavage ratio on average for gsu-miRNAs). Plant miRNA usually cleaves its target to inhibit target gene expression. However, some studies ( 15 - 17 ) including the gsu-miRNA-target Cytoscape networks showed that this might not be a strict one-to-one. In this study, targets of gsu-miRNAs can sometimes change, which may be beneficial for G. sulphuraria to respond to diverse stresses. Similar speculation has been made in miRNAs of C. reinhardtii under nutrient deprivation ( 9 ). 

Based on GO analysis results, target genes (GO:0071704 and GO:0022900) may be related to the synthesis of some important organic materials in ETC of G. sulphuraria. Target genes (GO:0051540 and GO:0046872) may be related to G. sulphuraria absorption properties for metal ions in wastewater ( 2 ). These genes could provide new clues for environmental governance based on genetic improvement of G. sulphuraria.

The pathway RNA transport (enrichment ratio= 170.73 in conserved gsu-miRNAs) is an essential BP. Different classes of RNAs experience various transport pathways between the nucleus and cytoplasm. Further investigation of the co-regulation of RNAs will elucidate detailed RNA transport processes in G. sulphuraria. Another pathway PPP (enrichment ratio= 529.66 in novel gsu-miRNAs), its metabolites, some key enzymes, is not inactivated under extreme conditions. The detail glycometabolism of G. sulphuraria under extreme environment need to be further clarified.

6. Conclusion

134 miRNAs from 124 MIRNA families were identified in G. sulphuraria. Twenty-one gsu-miRNAs were highly conserved. Some gsu-miRNAs were validated successfully by Northern blot and RT-PCR. Altogether, 1,589 putative miRNA targets were predicted in G. sulphuraria. GO, KEGG pathway, and Cytoscape network results will provide important references for further research and utilization of this microalga.

Acknowledgement

The authors would like to thank Dr. Lindsay L. (WILEY) and for the editorial assistance with the English. All data used in this study were available in the supplementary file. Funding: This work was supported by the National Natural Science Foundation of China (Grant No. 31670208), the Applied Basic Research Programs of Shanxi Province of China (Grant No. 201801D221242) and the 1311 Project of Shanxi.

Conflict of Interest: None.

References

1. Thangaraj B, Jolley CC, Sarrou I, Bultema JB, Greyslak J, Whitelegge JP, et al. Efficient light harvesting in a dark, hot, acidic environment: the structure and function of PSI-LHCI from Galdieria sulphuraria. Biophys J. 2011; 100(1):135-143. DOI
2. Yoon HS, Müller KM, Sheath RG, Ott FF, Bhattacharya D. Defining the major lineages of red algae Rhodophyta. J Phycol. 2006; 42:482-492. DOI
3. Selvaratnam T, Pegallapati AK, Montelya F, Rodriguez G, Nirmalakhandan N, vanVoorhies W, et al. Evaluation of a thermo-tolerant acidophilic alga, Galdieria sulphuraria, for nutrient removal from urban waste waters. Bioresour Technol. 2014; 156:395-399. DOI
4. Sakurai T, Aoki M, Ju X, Ueda T, Nakamura Y, Fujiwara S, et al. Profiling of lipid and glycogen accumulations under different growth conditions in the sulfothermophilic red alga Galdieria sulphuraria. Bioresour. Technol. 2016; 200:861-866. DOI
5. Jain K, Krause K, Grewe F, Nelson GF, Weber AP, Christensen AC, et al. Extreme features of the Galdieria sulphuraria organellar genomes: a consequence of polyextremophily?. Genome Biol Evol. 2014; 7(1):367-380. DOI
6. Schönknecht G, Chen WH, Ternes CM, Barbier GG, Shrestha RP, Stanke M, et al. Gene transfer from bacteria and archaea facilitated evolution of an extremophilic eukaryote. Science. 2013; 339(6124):1207-1210. DOI
7. D’Ario M, Griffiths-Jones S, Kim M. Small RNAs: big impact on plant development. Trends Plant Sci. 2017; 22(12):1056-1068. DOI
8. Huang A, Wu X, Wang G, Jia Z, He L. Computational prediction of microRNAs and their targets from three unicellular algae species with complete genome sequences. Can J Microbiol. 2011; 57(12):1052-1061. DOI
9. Voshall A, Kim EJ, Ma X, Yamasaki T, Moriyama EN, Cerutti H. miRNAs in the alga Chlamydomonas reinhardtii are not phylogenetically conserved and play a limited role in responses to nutrient deprivation. Sci Rep. 2017; 7(1):5462. DOI
10. Hafner M, Landgraf P, Ludwig J, Rice A, Ojo T, Lin C, et al. Tuschl. Identification of microRNAs and other small regulatory RNAs using cDNA library sequencing. Methods. 2008; 44(1):3-12. DOI
11. Axtell MJ, Meyers BC. Revisiting criteria for plant microRNA annotation in the era of big data. Plant Cell. 2018; 30(2):272-284. DOI
12. Navarro FJ, Baulcombe DC. miRNA-mediated regulation of synthetic gene circuits in the green alga Chlamydomonas reinhardtii. ACS Synth Biol. 2019; 8(2):358-370. DOI
13. Meyers BC, Axtell MJ, Bartel B, Bartel DP, Baulcombe D, Bowman JL, et al. Criteria for annotation of plant microRNAs. Plant Cell. 2008; 20(12):3186-3190. DOI
14. Zhang X, Liu X, Liu C, Wei J, Yu H, Dong B. Identification and characterization of microRNAs involved in ascidian larval metamorphosis. BMC Genomics. 2018; 19(1):168. DOI
15. Yang R, Chen G, Peng H, Wei D. Identification and characterization of miRNAs in Coccomyxa subellipsoidea C-169. Int J Mol Sci. 2019; 20(14):pii E3448. DOI
16. Xie FL, Huang SQ, Guo K, Xiang A, Zhu YY, Nie L, et al. Computational identification of novel microRNAs and targets in Brassica napus. FEBS Lett. 2007; 581(7):1464-1474. DOI
17. Liu T, Fang C, Ma Y, Shen Y, Li C, Li Q, et al. Global investigation of the co-evolution of MIRNA genes and microRNA targets during soybean domestication. Plant J. 2016; 85(3): 396-409. DOI