Novel long small RNAs (lsRNAs) are evolutionary conserved in plants and differentially expressed during epitype inducing conditions in Norway spruce embryos

doi:10.21203/rs.3.rs-3782103/v1

Download PDF

Article

Novel long small RNAs (lsRNAs) are evolutionary conserved in plants and differentially expressed during epitype inducing conditions in Norway spruce embryos

https://doi.org/10.21203/rs.3.rs-3782103/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Small non-coding RNA (sRNA) and the various other RNA species have tissue specific expression profiles. However, some size fractions may be missed or left out in the sequencing library preparation process, sequence data processing and downstream bioinformatic analysis. Here, we performed sRNA profiling in embryonic and vegetative tissues to elucidate a more complete picture of the gymnosperm sRNA populations than previously reported. We detected a novel group of sRNA between 31–34 nt in our Norway spruce sequencing data, with a prevalence in the 32–33 nt fraction, defined as the long small RNAs (lsRNAs). These lsRNAs were predominantly present in embryonic samples of Norway spruce. An in-silico analysis indicated that these lsRNA could originate from transfer RNAs (tRNA) and other non-coding transcripts and thus may target other non-coding RNAs (including tRNAs) or target repetitive elements such as transposons. We identified 18 putative orthologs of genes involved in Piwi-interacting RNA (piRNA) biogenesis but no spruce proteins were true homologs to the animal Piwi-proteins, thus lsRNAs cannot be considered as plant analogs to the animal piRNAs. Among the lsRNAs, tRNA-derived sequences from Asp, Glu and His iso-acceptors were in a majority and these sequences showed 3’ or 5’- bias dependent on the iso-acceptor type putatively targeted.

The lsRNA sized fraction was detected in seeds of all the gymnosperms examined and in Arabidopsis thaliana, suggesting that these comprise a conserved type of sRNAs between gymnosperms and angiosperms. lsRNA levels differed significantly among tissue types and developmental stages, and interestingly their expression was impacted by epitype-inducing temperature conditions. The lsRNAs add to the complexity of the small RNA world and may play a role in epigenetic regulation of gene expression in plants.

Biological sciences/Genetics/Epigenetics

Biological sciences/Genetics/Gene expression

Biological sciences/Genetics/Gene regulation

Biological sciences/Genetics/Plant genetics

Biological sciences/Genetics/Sequencing

gymnosperms

arabidopsis

small RNA

tRNAs

epigenetic memory

sRNA-seq

The largest part of the transcriptomes in plants comprises non-coding RNAs (ncRNAs) (Costa, 2010). Although we are now able to sequence ncRNA of all lengths, the entire ncRNA repertoire and their genomic origins remains to be fully elucidated, particularly for the non-canonical length fractions obtained.

Non-coding RNAs are RNA molecules that are not translated into proteins. The ncRNAs are often functionally important and include transfer RNAs (tRNAs) and ribosomal RNAs (rRNAs), as well as small RNAs such as microRNAs (miRNAs), small interfering RNAs (siRNAs), piwi-interacting RNAs (piRNAs), small nucleolar RNAs (snoRNAs), small nuclear RNA (snRNA), extracellular RNA (exRNA), small Cajal body-specific RNAs (scaRNAs), and the long non-coding RNAs (lncRNAs) (Li et al., 2012). However, non-coding RNAs that regulate protein coding genes can, for simplicity, be divided into just two smaller groupings: small RNAs (sRNA), shorter than 100 nt, and long non-coding RNAs (lncRNAs), longer than 200 nt (Röther and Meister, 2011; Lee, 2012).

Despite the fact that we are now able to detect sRNAs and lncRNAs of various lengths, their exact origins, biogenesis, function in the regulation of gene expression and level of activity remains unknown for at least some of the observed size fractions (Costa, 2010; Röther and Meister, 2011). Multiple pathways generate sRNAs in plants and have been divided into four classes: micro RNA (miRNA), trans-acting small interfering RNA (tasiRNA), natural antisense transcripts RNA (natsiRNA) and repeat-associated sequences RNA (rasiRNA) (Vazquez and Hohn, 2013; Borges and Martienssen, 2015). In a different classification, Axtell (2013) divided all small RNAs in plants into two main groups, those derived from single-stranded precursors with a hairpin structure, referred as hairpin RNAs (hpRNAs) and those derived from double-stranded RNA (dsRNA) precursors, referred to as small interfering RNAs (siRNAs). hpRNAs in plants are further subdivided into two secondary groups, the microRNAs and those that are not microRNAs. The dsRNA is composed of several subgroups such as repeat-associated small interfering RNAs (rasiRNAs) and trans-acting small interfering RNAs (tasiRNAs) (Carthew and Sontheimer, 2009; Liu et al., 2014). The best studied of these subgroups are the miRNAs and tasiRNAs, which direct translational repression and/or degradation of messenger RNAs, playing important roles in development and stress responses (Jones-Rhoades et al., 2006). The rasiRNAs, derived from repetitive DNA elements such as transposons and centromeres, promote heterochromatin formation to control transcripts emerging from repeated sequences and silenced transposons (Klattenhoff et al., 2007). Thus, there is a large diversity in small RNAs described in plants and more are predicted to be discovered in the future (Li et al., 2012; Alves and Nogueira, 2021).

Transcriptomic and genomic sequencing coupled with bioinformatics tools allow for better insight into the predicted genomic origin, content and diversity of the transcriptome. Sequencing of the small RNA transcriptome in Norway spruce embryogenic tissues has revealed longer small RNA peak with lengths of 31–32 nt, which had not been widely described in plants (Yakovlev and Fossdal, 2017). The same length class was described in animal germ line cells as the PIWI-interacting small RNA (piRNAs), not known to be present in plants (Siomi et al., 2011). In plants, a specific class of 23–27 nt long miRNAs was previously reported inArabidopsis thaliana (Gozmanova et al., 2017). Similarly, 27–30 nt sRNAs are present in unicellular eukaryotes (Lee and Collins, 2006), green algae (Li et al., 2023) and pathogenic fungi (Wang et al., 2016; Jia et al., 2017).

The aims of this comparative study, therefore, were: 1- to confirm the presence of larger than canonical (21–24 nt) fraction of sRNAs longer than 30 nt in 9 species as representation of the main gymnosperm groups (Order Pinales: Araucaria cunninghamii, Picea abies, Pinus sylvestris, Sciadopitys verticillata; Order Gnetales: Ephedra fragilis, Gnetum gnemon, Welwitschia mirabilis; Order Cycadales: Cycas revoluta; Order Ginkgogales: Gingko biloba) and 1 angiosperm (Order Brassicales: Arabidopsis thaliana), 2- predict their putative origin and identify their targets in theP. abies genome, 3- examine their expression in somatic embryos ofP. abies epitypes and 4- examine their transcription levels in seeds of the studied species.

Characterization of small RNAs of gymnosperm species

For the nine gymnosperm and an angiosperm species, more than 37.8 million reads of sRNA with lengths between 8–77 nt were obtained from 30 libraries representing 13 different samples (Table S1). Afterin silico filtering of sRNA with length 19–37 nt, we retained over 23.2 mln. reads. Distribution of sRNAs read length (19–37 nt) for the combined 18 libraries of NorwayP. abies somatic embryos is shown in Fig. 1. It is clear visible that the 20–24 ntfraction of the sRNAome was the most abundant, covering 56% of the sRNAs with lengths between 19–37 nt. There were also distinct sRNAs of lengths 31–34 nt, with a majority of reads at 32 to 33 nt (representing 11% of the sRNAs); this fraction was the center of attention in this study. We defined this 32 to 33 nt sRNA fraction in P. abies as long small RNAs (lsRNA). By read percentage, there is a significant difference (p < 0.05) for 32-nt fraction betweenP. abies embryos and buds and seeds (Suppl. Fig.S1).

To check for tissue specificity, we also compared the sRNAs length distribution between seeds and terminal buds of both P. abies and P. sylvestris (Suppl. Fig.S1). sRNAs from seeds of both species contained reads with lengths between 31 and 34 nt, peaking at 32 nt and 33 nt. In contrast, even if present, no distinct peaks were observed in the 31–34 nt size range for buds in both species (Suppl. Fig. S1).

In addition, we analyzed the read length distributions of sRNAs in whole seeds of nine gymnosperms and the angiosperm Arabidopsis thaliana, within the 19–37 nt range of sRNAs lengths (Fig. 2). In A. thaliana, there was a peak in reads of 32 nt in length. For most of the gymnosperm beside P. abies, peaks of reads with lengths between 30 and 36 nt were also observed (Fig. 2). However, each species had different prevalent length fractions of lsRNAs. All fourPinales species – P. abies, P. sylvestris, S. verticillata and A. cunninghamii, the most frequently lsRNAs showed two peaks at 32 and 35 nt, while inGnetales species the peak was at 35 nt, in E. fragilis and G. biloba in at 34 nt and in C. revoluta at 30 and 32 nt, and W. mirabilis at 31 nt.

Filtering reads with lengths between 28–37 nt resulted in 6 million reads, representing almost a 25% retention of the total number of small RNAs with lengths between 19–37 nt. From these, the most highly expressed sRNA tags (lsRNAs families) with more than 10 reads in any of the 10 species, were selected. These selected sRNAs included 1608 unique tags represented by over 4 million reads denoted as lsRNAs each marked with a specific ID code (slR#####) and were used in further analysis (Table S1.1). Distribution of lsRNAs lengths is shown in Fig. 3.

Defined lsRNA tags were used for annotation of lsRNAs in all the studied species, which allowed estimates of the extent of their conservation. We found that 221 lsRNA tags were present in transcriptomes of all 10 species, 252 lsRNAs tags were conserved among the nine gymnosperms and approximately five hundred tags were common to specific orders (Table 1). Surprisingly, the number of unique tags, presented uniquely in one single species, was low, less than 10, except forP. abies (50) and the monotypic gymnosperm genus Welwitschia (86). For P. abies, these specific tags were expressed primarily in the embryos. The distribution length of these lsRNAs in the transcriptomes of seeds in species representing the main gymnosperms orders andA. thaliana are shown in Suppl. Fig. S2. There was a prevalence of lsRNAs with lengths of 32 and 33 nt in all species studied, except the unique desert plant W. mirabilis, where the prevailing lsRNAs were 36 nt in length.

Table 1

Analysis of expressed long small RNAs (lsRNAs) in seeds of several gymnosperms’ species and the angiosperm *Arabidopsis thaliana*
Clade	Order	Species	Number of lsRNA tags	Number of common lsRNAs			Number of unique lsRNAs
Clade	Order	Species	Number of lsRNA tags	All species	Gymno-sperms	by orders	Number of unique lsRNAs
Gymnosperms	Pinales (Coniferales)	Picea abies	1292	221	252	573	50
		Pinus silvestris	1131				3
		Araucaria cunninghamii	856				2
		Sciadopitys verticillata	739				0
	Gnetales	Ephedra fragilis	711			485	1
		Gnetum gnemon	513				1
		Welwitschia mirabilis	1091				86
	Cycadales	Cycas revoluta	649			–	9
	Ginkgoales	Gingko biloba	953			–	2
Angiosperms	Brassicales	Arabidopsis thaliana	770		–	–	8

It is important to note the high redundancy of lsRNAs. There were 12 lsRNAs which represented nearly 43% of these transcripts. Among them, the lsRNAs which were most highly expressed in any one of the 10 species were: slR0525 and slR0001 in E. fragilis and A. thaliana ; slR1533 and slR0002 in P. sylvestris and S. verticillata; slR0057 in W. mirabilis ; slR0059 in E. fragilis; slR0010 in E. fragilis and G. biloba ; slR0036 in G. gnemon and slR0657 in A. cunninghamii (Table S6.1).

Analysis of the possible origin of lsRNAs inP. abies embryos

lsRNAs are not plant analogs to animal piRNAs

Animal species express PIWI-interacting RNAs (piRNAs), which are typically 24–32 nt long and are expressed predominantly in the germ line. As the 32 nt lsRNAs have similar length and were nearly absent in the vegetative tissues examined in the present work, we proposed that they might be the plant analogs of piRNAs. We therefore proceeded to obtain evidence for a piRNA-like biogenesis mechanism in P. abies. A detailed similarity search is presented in Supplementary Chap. 1.

We searched forP. abies homologs to the piRNA biogenesis genes in Drosophila. From the 47 gene models in theP. abies draft genome with high similarity to PAZ-domain proteins, we selected 26 full-length transcripts. Six gene models showed high similarity to animal Piwi proteins. However, a detailed amino acid analysis did not reveal the presence of the characteristic residues that distinguish the Piwi and Dicer subfamilies of Argonaute proteins from each other. Thus, it was concluded that theseP. abies proteins were not true orthologs to known animal Piwi proteins. Since Piwi protein coding genes are lacking, the lsRNAs cannot be plant analogs to animal piRNAs

However, in theP. abies genome we identified putative orthologs of genes involved in piRNA biogenesis in Drosophila, such as Tudor proteins (Tud) – 6 gene models, also Zucchini (Zuc) – 8, Shutdown (Shu) – 10, Cutoff (Cuff) – 1, Rhino (chromobox protein) – 1 and many helicases: VASA (DEAD-box helicases + HELICc ) ≥ 30 and UAP56 (ATP-dependent RNA helicase ) ≥ 6. Most of these piRNA biogenesis-like pathway genes were found to have significantly different transcript levels in P. abies embryos in comparison to buds. From 18 tested putative homologues, sixteen were differentially expressed. So, even if there are no Piwi proteins in P. abies, there are numerous piRNA biogenesis protein orthologs, which could potentially be involved in processing of longer fractions of plant sRNAs (such as lsRNA) by an unknown mechanism different to that known for animal piRNAs.

lsRNAs may be derived from non-coding and repetitive element transcript

For analysis of putative sources of lsRNAs, a BLAST search was performed to map all 31–37 nt lsRNAs on databases and transcript libraries. High levels of sequence correspondence, allowing up to only 3 mismatches, gave indications of the putative origins of the lsRNAs. The full list of lsRNA families, indicating putative origins, is presented in Suppl. Table S3.9.

A total of 499 lsRNAs (31%) mapped to 70 gene models defined in the spruce genome project (Table S3.1). As most of these predicted gene models are hypothetical proteins or not annotated (without significant similarities in the NCBI databases), they may be specific for lsRNA production. In parallel, 724 lsRNAs were mapped onto 44 genic contigs, containing predicted pseudogenes (Table S3.2). The mapping revealed that some lsRNAs mapped to both gene models and genic contigs, the majority in the complementary direction.

In addition, 997 lsRNAs mapped onto 248 assembled contigs from previously published transcriptomic data reads with no match to gene models (Table S3.3). With these novel matches, we now interpret these contigs as non-coding transcripts, e.g. putative lncRNAs, tRNA precursors and other forms of ncRNAs. Finally, 92 lsRNAs did not map to any available spruce transcripts or genic sequences, mainly because of the very draft state of the available spruce genome. 47 of them were both highly and differentially expressed in spruce embryos (compared to needle levels).

The lsRNAs were diversely distributed among exonic, intronic, intergenic and repeat sequences of the spruce genome when mapped, and sRNA clusters often overlapped with miRNAs loci (Suppl. Table S3.7). We found 179 lsRNAs matching to precursors of 6 plant miRNAs, and 15 lsRNAs matching the precursors of 3 known animal miRNAs (Suppl. Table S3.6). In some cases, lsRNAs overlapped with the miRNAs and, in others, mapped onto the “reverse” strand of their hairpin. Precursors for these miRNAs mapped onto the spruce genome sequence partially, or in full, at multiple locations. Therefore, such lsRNAs could be companion products of miRNA loci or be specific products from these loci under specific conditions.

lsRNAs may be derived from tRNAs

To better address the likely origin of the highly frequent sequences, we conducted an additional homology search for highly abundant lsRNAs (with more than 300 read counts inP. abies embryos) against a Repbase database. This across-species search enabled us to detect that a large number of reads had predicted homology to Oryza sativa tRNAs, but these highly accumulated lsRNAs showed few matches to known transposable elements and repetitive elements.

Subsequently, the search was focused on tRNA sequences and the 28–37 nt sequences were blasted against the tRNA dataset GtRNAdb. Of 1479 sequences with ≥ 300 reads, almost half (726 sequences) had similarity to tRNA sequences in both A. thaliana andP. abies (Fig. 4A). One third of 31–37 nt lsRNA sequences in seeds of gymnosperm species had high similarity with tRNA fragments and could potentially be derived from tRNAs (Fig. 5). Thus, tRNA-derived sequences may represent a substantial part of the lsRNA populations. In addition, the tDRnamer tool was used to check the similarity and classify some of the lsRNAs as tRNA-derived small RNAs (Table S3.4).

Sequences with similarity to Asp(GUC) tRNA formed a large fraction of the 31–37 nt sRNA population. Asp(GUC) tRNA sequences may be the origin of over 40% of the 31–37 nt sRNA population with ≥ 10 tags (total 5,919,448 read counts) (Fig. 4B). The 31 nt is the most frequent size with a match to Asp (GUC). Additionally, lsRNA sequences matching Gly(CCC), Gly(GCC), Gly(UCC), and His(GUG), with close to 30,000 tag counts were the second largest population. For these sequences, the 32–33 nt length dominated. Indeed, sequences matching tRNA may have explained close to 85% of the 28–37 nt lsRNA population with ≥ 300 tag counts.

To determine whether a similar match to the tRNA population was observed in the angiosperm, we analyzed the 31–37 nt sRNAs population fromA. thaliana seeds. The 31–37 nt population with more than 5 tag counts inA. thaliana showed the highest ratio of Asp(GUC), followed by Gly(CCC), Gly(GCC), Gly(UCC), His(GUG) and Glu(UUC) (Fig. 4C). This trend was similar to that observed in gymnosperm seeds. Thus, the tRNA-derived sequences of the 28–37 nt population are conserved between the angiosperm and the gymnosperms. Since every organism has its own codon usage frequency, we tried to address whether this unequal frequency in slRNA was derived from different frequency of codon usage. A comparison was made of the codon frequency and read counts with tRNA homology using the Codon usage database (Nakamura et al., 1999). Unexpectedly, the read counts did not correlate with codon frequency in protein-coding genes in eitherP. abies or in A. thaliana .

Considering that the typical tRNA is 70 – 90 bp in length, the 28–37 nt size was much smaller than canonical tRNA sequences. To test whether the position bias was observed in the 28–37 nt sRNA, we aligned 28–37 nt tRNA against the rice tRNA sequence (Fig. 4D). The sequences of the 28–37 nt population matched to 5´ or 3´ half position of tRNA and lacked an anti-codon region. Moreover, most tRNA-derived sequences had a 5´ or 3´ bias. For example, the most frequent Asp(GUC) was 5´ biased. The two different His(GUC), HisGUC.1 and HisGUC.2, were both 5´ biased. On the other hand, three different Gly(GCC), GlyGCC.1, and GlyGCC.2 were 3´ biased whereas GlyGCC.3 was not. Hence, the tRNA-derived lsRNA population in seeds was highly position-biased, implying that, potentially, most of these tRNA-derived lsRNAs were specifically generated. Due to the size of the 28–37 nt population, these sRNA may best be classified as tRNA halves.

We find it likely that tRNAs could be the direct source for at least some of the 32 nt lsRNAs produced. 599 lsRNAs (37%) were mapped to 428 tRNA sequences from 4 plant species: spruce (267 tRNAs), poplar (51), arabidopsis (49) and maize (61) (Table S3.5).

lsRNAs expression profiles

Differential expression of lsRNAs under EpI temperatures in P. abies

To confirm the importance of lsRNAs during the development ofP. abies embryos and assess their putative participation in the epigenetic memory formation, comparisons were made of transcript profiles in mature somatic embryos (SE) of two clones A2C and B10W ofP. abies after growing at three different epitype inducing (EpI) temperatures (18°C, 23°C, 28°C) (Fig. 6). For expression analysis, all 1608 defined lsRNA tags were used. Differentially expressed lsRNAs (DES) were identified through the pair-wise comparison of samples by setting the threshold |log₂ RPM ratio| to greater than 1 and p-value < 0.05 and these were examined further.

We found that lsRNA expression in SE was responsive in developing embryos under different temperature conditions. From the 1526 sRNA tags expressed in embryonic tissues, 1130 were DES at different known EpI temperature treatments. All DES were present in all SE samples, and no DES were specific to a single EpI temperature condition.

Based on their transcript profiles at the three different EpI temperatures, 1130 DES slRNA were distributed into 4 groups containing 18 clusters according to their expression trend (Fig. 6):

Group 1: 272 lsRNAs, clusters 1–6; lsRNAs related to high EpI temperature, with highest transcript levels at 28°C and lowest transcript levels at 18°C and 23°C;

Group 2: 562 lsRNAs, clusters 7–13; lsRNAs related to low EpI temperature, with highest transcript levels at 18°C and decreasing transcript levels at 23°C and 28°C;

Group 3: 228 lsRNAs, clusters 14–17, lsRNAs related to intermediate EpI temperature, with highest transcript levels at 23°C and lower transcript levels at 28°C and 23°C;

Group 4: 68 lsRNAs, 1 cluster of sRNA highly expressed at 18°C and 28°C with low expressed at 23°C (opposite to group 3).

Thus, upregulated lsRNAs in either low or high temperature conditions were the majority of DES. The main sRNA transcript profiles are shown in Fig. 6 and detailed description of clusters presented in Suppl. Table S4.1.

Differential lsRNAs levels between tissues and species

To assess the specificity of lsRNA expression in embryonic tissues, we compared the expression profiles in SE, seeds and buds ofP. abies with seeds and buds of P. sylvestris. In addition, we also compared transcript levels of defined lsRNAs in the nine gymnosperm species and A. thaliana. Analysis of DES between different tissues inP. abies showed that over 1500 lsRNAs (more than 90%) differed in levels between embryonic (SE, seeds) and vegetative tissues (buds) (Table S5). 129 lsRNA were highly specific to and expressed only inP. abies embryos (Table S4.2).

Comparison of lsRNA profiles in the seeds of the 9 different gymnosperm species andA. thaliana confirmed the high level of conservancy of lsRNAs among plant species. 241 lsRNAs were common to all gymnosperm species, although their expression profiles differed among the species. The number of expressed lsRNAs varied from 513 in Gnetum to 1292 inP. abies (Suppl. Table S6.1). As indicated earlier, the prevailing lengths of lsRNAs were 32–33 nt, except for G. gnemon (35 nt), W. mirabilis (36 nt) and C. revoluta (31 nt) (Fig. 7, Table S6.2).

In silico prediction of putative lsRNAs targets

For the assignment of functional roles to the entire set of lsRNAs, we tried to predict target gene transcripts forP. abies andA. thaliana (for comparison). Available gene models were applied, using the plant small RNA target analysis server psRNATarget, adjusted for 32 nt sRNAs. Full results for these sets are presented in Supplementary Table S7.

In total for P. abies, 144 lsRNA families targeted 54 protein coding gene models with very high redundancy. Several lsRNAs targeted the same gene model and vice versa, some lsRNA may target several gene models, comprising the 1406 lsRNA–target gene model pairs (Table S7.1). Only 4 target gene models were annotated, a core histone protein, Sodium Bile acid symporter family, EARLY FLOWERING 7-like and a BTB/POZ domain with WD40/YVTN repeat-like protein. Most gene models were unknown and defined as hypothetical proteins. Although theP. abies repetitive sequences (RE) collection is not comprehensive, 2 transposons were detected as putative targets for 8 lsRNA families (Suppl. Table S7.2).

Based on poplar DPTEdb database, 314 lsRNAs could potentially target 56 transposable elements (Table S7.3). Similarly to gene models mentioned above, there was great redundancy in lsRNA – target pairs. Several lsRNAs may target the same TE and vice versa, the one and same lsRNA could target several TEs, totaling 6393 lsRNA–TE pairs. Most putative targets represented retrotransposons (LTR/copia, gypsy, PAO, ERV, MuDR), DNA transposons (Helitron) and unknown REs.

The results were similar for A. thaliana, in which 110 lsRNA families targeted 54 protein coding gene models with very high redundancy (Table S7.4). Among annotated target gene models were bHLH DNA-binding superfamily protein, DNA polymerase subunit, Threonyl-tRNA synthetase and actin filament capping protein.

sRNAs longer than 24 nt found in plant NGS libraries are too often completely disregarded or filtered out during research, presumed to be unimportant degradation products in bioinformatic pipelines. However, sRNAs of non-canonical length may play important roles in gene regulation. In small RNA libraries obtained from somatic embryonic tissues of P. abies, we detected a longer fraction of sRNAs of 31–34 nt, not previously described in plants (Yakovlev and Fossdal, 2017). Our current analysis we report that within the RNAs range the most prevalent size was 24 nt sRNAs (comprising 21%), followed by 21 nt sRNAs (16%), while the 32–33 nt sRNA comprise 11% of the total reads inP. abies somatic embryos. The 32–33 nt sRNA are statistically overrepresented in embryos compared to buds and they were differentially expressed in embryos between different epitype inducing conditions. We here denote these 32–33 nt sRNA as long small RNA (lsRNA). The lsRNAs were present in seeds from all 10 gymnosperms and the angiosperm Arabidopsis thaliana examined, suggesting that these are conserved between all seed plants (gymnosperms and angiosperms) during evolution.

Different gymnosperm orders had lsRNAs of slightly different prevalent lengths.Pinales seed samples had peaks at 31 and 34 nt, inGnetales one at 34 nt, G. biloba at 33 nt and, in C. revoluta, lsRNAs appear shorter with peaks at 29 and 31 nt. It is curious that the lsRNAs prevalent lengths match with the phylogenetic position in gymnosperm orders of the gnepine hypothesis (Wang and Ran, 2014). The most basalCycadales have the shortest lsRNAs, then length and variety increase towards the extantPinales andGnetales sister groups. That this finding fits with the gnepine hypothesis might be a coincidence but may also reflect sRNA related processing machinery adaptations affecting lsRNA production that has occurred during and after the evolutionary divergence of these orders. Among the 9 gymnosperm species analyzed, the most highly expressed lsRNAs differed between species, which resonated with the conclusions of Chen et al. (2018), who noted that lineage-specific sRNAs were likely to perform specific functions in certain plants or groups (Chen et al., 2018).

In the seed of the angiosperm A. thaliana, we found analogous sized sRNA with a peak at 31 nt (putative lsRNA). This finding was a surprise, since theA. thaliana small RNAs was apparently already well described and with just a few reports noting the existence of 30–40 nt sRNAs, only in response to treatment with pathogenic bacteria (Katiyar-Agarwal et al., 2007). These authors denoted them as long siRNAs (lsiRNAs). One of these lsiRNAs, AtlsiRNA-1, represents a natural antisense transcript (natsRNA) specifically induced in response to the bacterium Pseudomonas syringae that carries the effector avrRpt2 (Katiyar-Agarwal et al., 2007). The fact that longer than 30nt sRNAs are typically removed during sRNA library preparations or filtered out during analysis, may be part of the reason that this 32 nt peak has not received attention inA. thaliana seed studies. Despite this omission, we have described for the first time a series of lsRNA tags that are present in gymnosperms and angiosperms that point out the potential importance of lsRNA in these taxa. Taking into account that gymnosperms and angiosperms diverged from a common ancestor ca. 500 million years ago (Vergara-Silva et al., 2000) and the living orders of gymnosperms separated from each other 130 million years ago (Chaw et al., 2000), it is striking that more than 200 specific lsRNA were found in all the 10 species analyzed and more than 250 were shared among the 9 gymnosperms. This finding might reflect the fact that highly conserved sRNAs tend to play key roles in regulation of plant growth and development (Zeng et al., 2019).

The most well studied examples of longer sRNAs are PIWI-interacting RNAs (piRNAs)(Ishizu et al., 2012; Sato and Siomi, 2013) and tRNA-derived sRNAs (Röther and Meister, 2011). piRNAs are usually between 26 and 31 nt long, are mainly derived from transposons, retroelements, and repetitive loci, and only present in germline cells. piRNAs are characterized by the presence of a 5´ uridine and their association with PIWI proteins, the germline specific Argonaut subfamily (Sato and Siomi, 2013). Based mainly on their length, we initially interpreted the 32 nt peak lsRNA as possible plant analogs of the animal piRNAs (Yakovlev and Fossdal, 2017). However, we were unable to identify bona fide PIWI proteins homologs (Nonomura et al., 2007; Zheng et al., 2019) among the defined AGO gene models in the gymnosperm P. abies. Therefore, we should reject our initial hypothesis and now conclude that lsRNAs are not plant analogs to animal piRNAs.

Despite the lack of Piwi proteins in plants we were still able to identify a large set of proteins similar to those involved in piRNA biogenesis in animals (Ishizu et al., 2012). We observed that for 6 of 16 piRNA biogenesis-like the transcripts, levels differed inP. abies embryos (also rich in lsRNAs) compared to needles (with very low lsRNA levels). These gene products could potentially be involved in a yet to be undescribed sRNA biogenesis pathway analogous to, or completely different, to that reported for animal piRNAs.

The exact genomic origin of lsRNAs is enigmatic. Here we predict that lsRNAs derive from varied genomic regions. Direct mapping onto the Norway spruce genome hints at several possible sources that include coding parts of transcripts, introns, longer non-coding RNAs and from specific tRNAs.

Some lsRNA matches are grouped into clusters on the genome and often overlap with predicted siRNAs and miRNAs loci. Therefore, some lsRNAs may be companion products of miRNA and TAS loci or alternatively be specific products from these loci.

At least one third of the lsRNAs have strong similarity to and could originate from tRNAs. Considering that the typical tRNA is 70–90 bp in length, the 28–37 nt size was much smaller than canonical tRNA sequences. It has been reported that the difference between tRF/tRNA-halves and degraded tRNA species shows the position bias, with either 5´ or 3´ of tRNA part observed (Sharma et al., 2016).

tRNA-derived sRNAs include: ~18–28 nt small tRNA-derived fragments (tRFs) and ~ 30–35 nt tRNA-derived halves (tRHs), generated either from pre-tRNAs or mature tRNAs (Martinez, 2018; Zhu et al., 2018). tRDs and tRHs were initially described in animals, including humans (Röther and Meister, 2011). Such sRNA fragments are pronounced in yeast during oxidative stress and similar tRNA fragments were also observed in human cells and plants during oxidative stress using Northern blot analysis (Thompson et al., 2008). The tRHs size fits into the observed length of our lsRNAs. The tRHs can be generated from pre-tRNAs and tRNAs, by cleavage of the anticodon loop, which yields 5′tRHs and 3′tRHs products (Zhu et al., 2018). Zhang et al. (2009) purposely isolated and cloned 30–90 nt RNAs from the phloem sap of pumpkin and discovered several tRH such as 5′tRH-GlyGGA, 5′tRH-MetAUG, 3′tRH-AspGAC, 3′tRH-LysAAG, 3′tRH-GluGAG, 3′tRH-SerAGC, 3′tRH-SerACA, 3′tRH-ProCCT, and 3′tRH-ArgCGA. Similar tRNA-derived small RNAs (tsRNAs) were also identified in rice embryogenic callus, where the 5′tRF from AlaAGC-tRNA was the most abundant, comprising 82% of the total tsRNA pool (Chen et al., 2011). We found similar quantities in this study for both gymnosperms and the angiosperm. InA. thaliana we found the highest ratio of AspGUC, with the second highest being GlyCCC, GlyGCC, GlyUCC and HisGUG and GluUUC tRH-like lsRNAs. InP. abies embryos, the most abundant lsRNAs matched to AspGUC, followed by GlyCCC, GlyGCC, GlyUCC and HisGUG. Asp(GUC) and Gly(GCC) are known as substrates of DNMT2 (Goll et al., 2006). Drosophila Dnmt2 and human DNMT2 are potential RNA methyltransferases (Schaefer et al., 2010; Li et al., 2021).A. thaliana homologs of DNMT2 and NSUN2 are also responsible for tRNA methylation, in particular, Asp-, Gly-, Glu-, and Ser-tRNAs (Burgess et al., 2015). This might suggest that covalent modification in the RNA molecules is involved in the processing of these putative tRNA-derived sRNA or conversely linked to the regulation of RNA methylation. Hence, the tRNA-derived lsRNA population in seeds was highly position-biased, implying that, potentially, most of these tRNA-derived lsRNAs were specifically generated. Due to the size of the 28–37 nt population, these sRNA may best be classified as tRNA halves.

The tRDs and tRHs may have functional roles (Cao, 2023; Wang et al., 2023). Both tRFs purified fromA. thaliana leaves and specific synthetic tRFs deriving from the 5ʹ extremity ofA. thaliana tRNA^Ala and tRNA^Asn efficiently inhibit protein synthesis in vitro. tRFs likely act as general modulation factors of the translation process in plants (Lalande et al., 2020). Recent studies have shown that methylation of RNA by DNMT2 and NSUN2 is involved in RNA transportation from shoot to root inA. thaliana (Yang et al., 2019). Moreover, a rice homolog of Nsun2 appears to be responsible for high-temperature adaptation through translation regulation (Tang et al., 2020). tRHs may be stress- inducible both in animals and plants (Torrent et al., 2018; Hajieghrari and Farrokhi, 2022). The tRNA-derived small RNAs (tsRNAs)-mediated target silencing was further found to be dependent on AGO2 and the fact that 5′ tsRNAs were co-immunoprecipitated with AGO2 suggests that tsRNAs inhibit specific targets via RNA interference (RNAi) (Kumar et al., 2014; Shigematsu and Kirino, 2015; Jehn and Rosenkranz, 2019). Thus, tRNA-derived lsRNA might play a role in adaptive responses to environmental cues and specific stress conditions.

We found that lsRNAs were present naturally in somatic and zygotic embryos of plant seeds, but less abundant in vegetative tissue; however, no functional analysis was made in this study. To obtain insights into their putative role, we studied lsRNAs level changes during embryogenesis, a time when it is known thatP. abies embryogenic tissues are highly sensitive to changing environmental conditions (Kvaalen and Johnsen, 2008). We searched for a possible role for lsRNAs in epigenetic mechanisms by analyzing differentially expressed lsRNAs (DES) patterns in mature somatic embryos cultivated under contrasting epitype-inducing (EpI) temperature conditions. InP. abies epitype embryos, we found 1130 DESs between cooler and warmer EpI conditions grouped into 18 clusters with similar expression profiles. The EpI conditions resulted in a very clear change in lsRNA levels. Such change in expression highlights their dynamic nature, implying direct or indirect involvement in transcriptional and/or translational regulation and, possibly, in epigenetic mechanisms during embryogenesis.

As there are no tools designed specifically for prediction of targets for slRNAs, we tried to use psRNATarget, developed for miRNAs. It does not mean that lsRNAs function the same way as miRNAs, but help to define regions complementary to slRNA sequences, which could be targeted by lsRNAs. The target prediction of lsRNAs showed that only a small number of knownP. abies andA. thaliana gene models could be regulated directly by lsRNAs. Target prediction among repeat and transposable elements was modest. In total, there were 60 targets among the transposable elements. However, the TE databases contain mainly angiosperm species so our search has probably underestimated the TE and repetitive element hits by lsRNAs in gymnosperms. Despite the high degree of conservation of lsRNAs between plants, this finding points to an important role in biological regulation, although their exact role is yet to be determined. In addition, the very high predicted redundancy of lsRNAs means that one single target could potentially be regulated by many lsRNAs. Functional redundancy is thought to introduce a robustness that buffers the deleterious effect of mutations hitting one of the redundant copies (Ambrós et al., 2018).

The possible functional importance of lsRNAs is supported by their presence across different phyla that may have been evolutionary separated for up to 500 million years, namely between angiosperms and gymnosperms. In addition, lsRNAs are differentially expressed in matureP. abies somatic embryos growing under different epitype-inducing temperature conditions. Of course, not all small RNAs are necessarily functionally important and some lsRNAs may be transitory degradation intermediates or simply by-products of aberrant transcription. This being said, we do not rule out the possibility that lsRNAs may have a functional role during propagation, particularly in meristematic cells, but functional studies are needed to confirm or refute such a proposed role.

In conclusion, this is the first detailed in-silico survey of lsRNA, a class of embryo-related 28–37 nt sRNAs, suggesting that at least some of these are conserved between gymnosperm and angiosperm plants. lsRNAs have a prevalent length of 31–32 nt and we recommend that the sRNA fraction of greater than 24 nt in size should be considered potentially functional and not be filtered out or ignored in future RNA analyses. These long small RNAs add to the complexity of the small RNA world and, since they are differentially expressed during epitype inducing conditions during embryogenesis, some may have a direct or indirect role in response to environmental change and in epigenetics.

Plant material and RNA extraction

Total RNA containing all small RNAs was isolated from whole seeds of nine gymnosperm species and Arabidopsis thaliana as a representative of angiosperm plants (Table 2). Additionally, we extracted total RNA from apical buds (with some surrounding needles) of Picea abies and Pinus sylvestris to study the vegetative tissues, giving in total 12 samples (Table 2). No approvals were required for the study, which complied with all relevant regulations.

Table 2

General description of samples, used in current research.
Order	Species	Tissues	Number of libraries
Gymnosperms	Species	Tissues	Number of libraries
Coniferales (Pinales)	Picea abies Pinus silvestris Araucaria cunninghamii Sciadopitys verticillata	embryos, seeds, buds seeds, buds seeds seeds	18, 1, 1 1, 1 1 1
Cycadales	Cycas revoluta	seeds	1
Ginkgoales	Ginko biloba	seeds	1
Gnetales	Ephedra fragil is Gnetum gnemon Welwitschia mirabilis	seeds	1 1 1
Angiosperms	Arabidopsis thaliana	seeds	1

For P. abies, we also used 18 samples of total RNAs containing whole small RNAs fraction extracted from somatic embryos of two epitypes A2C and B10W (genetically identical epitype embryos) generated under one of the three epitype inducing temperature conditions (18°C, 23°C, 28°C) during embryogenesis (Kvaalen, Johnsen, 2008). Nine samples in total were collected from each epitype, representing three maturation stages and three different EpI-temperatures, providing 18 samples in total. Data of epitypes and maturation stages were combined together and considered as 6 biological replicates for each of EpI-temperatures.P. abies epitypes and samples for total RNA extraction are described in detail in our previous works (Yakovlev et al., 2016; Carneros et al., 2017; Yakovlev and Fossdal, 2017).

Each biological replica consisted on around 80 mg of fresh tissue were ground to a fine powder in liquid nitrogen using a mortar and pestle or in 1.5 ml Eppendorf tubes and total RNA extracted using an Epicentre® MasterPure™ Plant RNA Purification Kit (Epicentre, Madison, WI, USA, #MPR09100) following the manufacturer’s instructions. Total RNA preparations were stored at -80°C. The integrity and quantity of the total RNA was assessed using a Agilent 2100 Bioanalyzer with RNA 6000 Nano Kit (Agilent, Santa Clara, CA, USA #5067 − 1511).

Small RNA library construction, sequencing, and bioinformatics

The twelve new sRNA libraries from seeds of 10 species and buds ofP. abies and P. sylvestris, and the 18P. abies embryonic sRNA libraries were constructed in the same way. Each were constructed from 1 µg of total RNA using the Ion Total RNA-Seq Kit v2 for Small RNA Libraries (Catalogue #4476286) and the Ion OneTouch™ 2 System (Thermo Fisher Scientific Inc., Waltham, MA, USA) with Ion PGM™ Template OT2 Kit (Catalogue #4480974) according to the manufacturer’s instructions. Sequencing was performed using an Ion Personal Genome Machine® (PGM™) Sequencer (Thermo Fisher Scientific Inc., Waltham, MA, USA) and each library sequenced individually on a 316 version 2 chip using the Ion PGM™ sequencing 200 kit v2 (#4482006). Data processing was initially completed on an Ion Torrent Server using TorrentSuite software (v.4.0.1). The newly obtained sRNA data and sequencing results from the 18 small RNA libraries, constructed from somatic embryos ofP. abies epitypes, were processed using CLC genomics workbench software (v.23.0) (Qiagen Company, Hilden, Germany) and Excel software.

Identification of 28–37 nt lsRNAs and theirs mapping to theP. abies genome

The whole small RNA pool, combined from twelve libraries, was trimmed to remove single and reads which were too short or too long. The initial filtering criteria were: read frequency greater than 2 and length between 19–37 bp. Second filtering was carried out to select sRNAs with lengths between 28–37 bp and read frequencies higher than 100, to select the most highly expressed larger sRNA tags).

The filtered sRNA tags were named as slR###### and used for annotation and merging of individual library reads. We tried to match lsRNAs onto both coding and non-coding transcripts. To achieve this analysis, we blasted and mapped all defined lsRNAs on database and transcripts libraries. Tentatively named lsRNA tags were mapped to the draftP.abies genome sequence v.1 (Nystedt et al., 2013); PRJEB18220), allowing for 0–2 nucleotides of mismatch using ConGenIE (http://congenie.org/) tools (Nystedt et al., 2013). For non-coding transcripts, multiple databases of repetitive elements were used: Repbase database (Jurka et al., 2005; Bao et al., 2015), RepPop (Zhou and Xu, 2009) and DPTEdb (Li et al., 2016), and also tRNA databases like GtRNAdb: a database of transfer RNA genes detected in genomic sequences (Chan and Lowe, 2008; 2015). Improved nomenclature and tRNA-derived small RNAs naming was performed using the tDRnamer system (Holmes et al., 2022). Codon frequencies in tRNA homologs were compared using the Codon usage database CUTG (codon usage tabulated from GenBank) (Nakamura et al., 1999).

Target prediction

For the assignment of functional roles to the whole set of defined lsRNAs, we tried to predict target gene transcripts in some available gene models and for transposable/repetitive elements using the plant small RNA target analysis server – psRNATarget, adjusted for 32 nt lsRNAs (Dai et al., 2018). As a target for analysis we incorporated into psRNATargetP. abies gene models and repetitive elements (ER) collection (http://congenie.org), MIPS PlantsDB (Nussbaumer et al., 2013) and Dioecious Plants Transposable Elements (TE) Database – DPTEdb for Populus trichocarpa (http://genedenovoweb.ticp.net:81/DPTEdb/index.php) (Li et al., 2016).

DATA accession

Pinus sylvestris small RNA sequences were deposited to the SRA (Short Read Archive, NCBI) with the accession number PRJNA531446 (https://www.ncbi.nlm.nih.gov/sra/?term=PRJNA531446).

Picea abies small RNAs transcripts from embryos were deposited to the SRA with the following accession: submission ID SUB1781210; BioProject ID PRJNA339513, BioSamples: SAMN05592191–SAMN05592208 and accession IDs: SRX2033261 - SRX2033278 (https://www.ncbi.nlm.nih.gov/sra/?term=PRJNA339513).

Gymnosperms andArabidopsis thaliana sRNA transcripts from 24 libraries have the following accession numbers: submission ID: SUB12530980; BioProject ID PRJNA923163; BioSamples: SAMN32691192 – SAMN32691215 and accession IDs: SRR23060401 – SRR23060424 (http://www.ncbi.nlm.nih.gov/bioproject/923163).

sRNAs

small RNAs

lsRNAs

long small RNAs (i.e., 32–33 nt fraction of sRNAs)

repetitive element

transposable element

tRFs

tRNA derived sRNAs

EpI

epitype inducing temperature.

Author Contribution

CGF conceived and designed the research; IY collected the samples; IY and MN performed experiments; IY, MN, MV and YZ analyzed the data; IY, MN and CGF discussed the data and wrote the article; MV prepared the figures; All of the authors have read and approved the final version of this manuscript.

Acknowledgements

The authors thank Adam Vivian-Smith (NIBIO, Norway) for the delivery of gymnosperm seeds and technical help during sequencing. This work was carried out with financial support from the Norwegian Forest and Landscape Institute (GB #135117-1), ProCoGen project (FP7-KBBE-2011-5; # 289841).

Alves, C.S., and Nogueira, F.T.S. (2021). Plant Small RNA World Growing Bigger: tRNA-Derived Fragments, Longstanding Players in Regulatory Processes. Frontiers in Molecular Biosciences 8. doi: 10.3389/fmolb.2021.638911.
Ambrós, S., de la Iglesia, F., Rosario, S.M., Butković, A., and Elena, S.F. (2018). Engineered Functional Redundancy Relaxes Selective Constraints upon Endogenous Genes in Viral RNA Genomes. Genome Biology and Evolution 10(7), 1823-1836. doi: 10.1093/gbe/evy141.
Axtell, M.J. (2013). Classification and comparison of small RNAs from plants. Annual Review of Plant Biology 64(1), 137-159. doi: doi:10.1146/annurev-arplant-050312-120043.
Bao, W., Kojima, K.K., and Kohany, O. (2015). Repbase Update, a database of repetitive elements in eukaryotic genomes. Mob DNA 6, 11. doi: 10.1186/s13100-015-0041-9.
Borges, F., and Martienssen, R.A. (2015). The expanding world of small RNAs in plants. Nat Rev Mol Cell Biol advance online publication. doi: 10.1038/nrm4085.
Burgess, A.L., David, R., and Searle, I.R. (2015). Conservation of tRNA and rRNA 5-methylcytosine in the kingdom Plantae. BMC Plant Biology 15(1), 199. doi: 10.1186/s12870-015-0580-8.
Cao D (2023) Stress-induced tRNA fragments take action in alternative splicing in Arabidopsis. Plant Physiology 193: 877-879
Carneros E, Yakovlev I, Viejo M, Olsen JE, Fossdal CG (2017) The epigenetic memory of temperature during embryogenesis modifies the expression of bud burst-related genes in Norway spruce epitypes. Planta 246: 553-566
Carthew, R.W., and Sontheimer, E.J. (2009). Origins and Mechanisms of miRNAs and siRNAs. Cell 136(4), 642-655.
Chan, P.P., and Lowe, T.M. (2008). GtRNAdb: a database of transfer RNA genes detected in genomic sequence. Nucleic Acids Research 37(suppl_1), D93-D97. doi: 10.1093/nar/gkn787.
Chan, P.P., and Lowe, T.M. (2015). GtRNAdb 2.0: an expanded database of transfer RNA genes identified in complete and draft genomes. Nucleic Acids Research 44(D1), D184-D189. doi: 10.1093/nar/gkv1309.
Chaw, S.-M., Parkinson, C.L., Cheng, Y., Vincent, T.M., and Palmer, J.D. (2000). Seed plant phylogeny inferred from all three plant genomes: Monophyly of extant gymnosperms and origin of Gnetales from conifers. Proceedings of the National Academy of Sciences of the United States of America 97(8), 4086-4091.
Chen, C., Zeng, Z., Liu, Z., and Xia, R. (2018). Small RNAs, emerging regulators critical for the development of horticultural traits. Horticulture Research 5(1), 63. doi: 10.1038/s41438-018-0072-8.
Chen, C.-J., liu, Q., Zhang, Y.-C., Qu, L.-H., Chen, Y.-Q., and Gautheret, D. (2011). Genome-wide discovery and analysis of microRNAs and other small RNAs from rice embryogenic callus. RNA Biology 8(3), 538-547. doi: 10.4161/rna.8.3.15199.
Costa, F.F. (2010). Non-coding RNAs: Meet thy masters. Bioessays 32(7), 599-608. doi: 10.1002/bies.200900112.
Dai, X., Zhuang, Z., and Zhao, P.X. (2018). psRNATarget: a plant small RNA target analysis server (2017 release). Nucleic Acids Research 46(W1), W49-W54. doi: 10.1093/nar/gky316.
Goll, M.G., Kirpekar, F., Maggert, K.A., Yoder, J.A., Hsieh, C.L., Zhang, X., et al. (2006). Methylation of tRNAAsp by the DNA methyltransferase homolog Dnmt2. Science 311(5759), 395-398. doi: 10.1126/science.1120976.
Gozmanova, M., Baev, V., Apostolova, E., Sablok, G., and Yahubyan, G. (2017). "Growing Diversity of Plant MicroRNAs and MIR-Derived Small RNAs," in Plant Epigenetics, eds. N. Rajewsky, S. Jurga & J. Barciszewski. (Cham: Springer International Publishing), 49-67.
Hajieghrari, B., and Farrokhi, N. (2022). Plant RNA-mediated gene regulatory network. Genomics 114(1), 409-442. doi: https://doi.org/10.1016/j.ygeno.2021.12.020.
Holmes, A.D., Chan, P.P., Chen, Q., Ivanov, P., Drouard, L., Polacek, N., et al. (2022). A Standardized Ontology for Naming tRNA-derived RNAs Based on Molecular Origin. bioRxiv, 2022.2005.2006.490965. doi: 10.1101/2022.05.06.490965.
Ishizu, H., Siomi, H., and Siomi, M.C. (2012). Biology of PIWI-interacting RNAs: new insights into biogenesis and function inside and outside of germlines. Genes & Development 26(21), 2361-2373. doi: 10.1101/gad.203786.112.
Jehn, J., and Rosenkranz, D. (2019). tRNA-Derived Small RNAs: The Good, the Bad and the Ugly. Med One 4(4), e190015. doi: 10.20900/mo.20190015.
Jia, J., Lu, W., Zhong, C., Zhou, R., Xu, J., Liu, W., et al. (2017). The 25–26 nt Small RNAs in Phytophthora parasitica Are Associated with Efficient Silencing of Homologous Endogenous Genes. Frontiers in Microbiology 8. doi: 10.3389/fmicb.2017.00773.
Jones-Rhoades, M.W., Bartel, D.P., and Bartel, B. (2006). MicroRNAs and their regulatory roles in plants. Annual Review of Plant Biology 57(1), 19-53. doi: 10.1146/annurev.arplant.57.032905.105218.
Jurka, J., Kapitonov, V.V., Pavlicek, A., Klonowski, P., Kohany, O., and Walichiewicz, J. (2005). Repbase Update, a database of eukaryotic repetitive elements. Cytogenetic and Genome Research 110(1-4), 462-467.
Katiyar-Agarwal, S., Gao, S., Vivian-Smith, A., and Jin, H. (2007). A novel class of bacteria-induced small RNAs in Arabidopsis. Genes & Development 21(23), 3123-3134. doi: 10.1101/gad.1595107.
Klattenhoff, C., Bratu, D.P., McGinnis-Schultz, N., Koppetsch, B.S., Cook, H.A., and Theurkauf, W.E. (2007). Drosophila rasiRNA Pathway Mutations Disrupt Embryonic Axis Specification through Activation of an ATR/Chk2 DNA Damage Response. Developmental Cell 12(1), 45-55. doi: https://doi.org/10.1016/j.devcel.2006.12.001.
Kumar, P., Anaya, J., Mudunuri, S., and Dutta, A. (2014). Meta-analysis of tRNA derived RNA fragments reveals that they are evolutionarily conserved and associate with AGO proteins to recognize specific RNA targets. BMC biology 12, 78. doi: 10.1186/PREACCEPT-5867533061403216.
Kvaalen, H., and Johnsen, O. (2008). Timing of bud set in Picea abies is regulated by a memory of temperature during zygotic and somatic embryogenesis. New Phytologist 177(1), 49-59. doi: doi:10.1111/j.1469-8137.2007.02222.x.
Lalande S, Merret R, Salinas-Giegé T, Drouard L (2020) Arabidopsis tRNA-derived fragments as potential modulators of translation. RNA Biology 17: 1137-1148
Lee, J.T. (2012). Epigenetic regulation by long noncoding RNAs. Science 338(6113), 1435-1439. doi: 10.1126/science.1231776.
Lee, S.R., and Collins, K. (2006). Two classes of endogenous small RNAs in Tetrahymena thermophila. Genes & Development 20(1), 28-33. doi: 10.1101/gad.1377006.
Li, H., Zhu, D., Wu, J., Ma, Y., Cai, C., Chen, Y., et al. (2021). New substrates and determinants for tRNA recognition of RNA methyltransferase DNMT2/TRDMT1. RNA Biol 18(12), 2531-2545. doi: 10.1080/15476286.2021.1930756.
Li, S.F., Zhang, G.J., Zhang, X.J., Yuan, J.H., Deng, C.L., Gu, L.F., et al. (2016). DPTEdb, an integrative database of transposable elements in dioecious plants. Database (Oxford) 2016. doi: 10.1093/database/baw078.
Li, Z., Ender, C., Meister, G., Moore, P.S., Chang, Y., and John, B. (2012). Extensive terminal and asymmetric processing of small RNAs from rRNAs, snoRNAs, snRNAs, and tRNAs. Nucleic Acids Research 40(14), 6787-6799. doi: 10.1093/nar/gks307.
Li, Y., Kim, E.-J., Voshall, A., Moriyama, E.N., and Cerutti, H. (2023). Small RNAs >26 nt in length associate with AGO1 and are upregulated by nutrient deprivation in the alga Chlamydomonas. The Plant Cell. doi: 10.1093/plcell/koad093.
Liu, Y.-X., Wang, M., and Wang, X.-J. (2014). Endogenous Small RNA Clusters in Plants. Genomics, Proteomics & Bioinformatics 12(2), 64-71. doi: https://doi.org/10.1016/j.gpb.2014.04.003.
Martinez, G. (2018). tRNA-derived small RNAs: New players in genome protection against retrotransposons. RNA biology 15(2), 170-175. doi: 10.1080/15476286.2017.1403000.
Nakamura, Y., Gojobori, T., and Ikemura, T. (1999). Codon usage tabulated from the international DNA sequence databases; its status 1999. Nucleic Acids Research 27(1), 292-292. doi: 10.1093/nar/27.1.292.
Nonomura, K., Morohoshi, A., Nakano, M., Eiguchi, M., Miyao, A., Hirochika, H., et al. (2007). A germ cell specific gene of the ARGONAUTE family is essential for the progression of premeiotic mitosis and meiosis during sporogenesis in rice. Plant Cell 19(8), 2583-2594. doi: 10.1105/tpc.107.053199.
Nussbaumer, T., Martis, M.M., Roessner, S.K., Pfeifer, M., Bader, K.C., Sharma, S., et al. (2013). MIPS PlantsDB: a database framework for comparative plant genome research. Nucleic Acids Research 41(D1), D1144-D1151. doi: 10.1093/nar/gks1153.
Nystedt, B., Street, N.R., Wetterbom, A., Zuccolo, A., Lin, Y.-C., Scofield, D.G., et al. (2013). The Norway spruce genome sequence and conifer genome evolution. Nature 497, 579–584. doi: 10.1038/nature12211.
Röther, S., and Meister, G. (2011). Small RNAs derived from longer non-coding RNAs. Biochimie 93(11), 1905-1915. doi: https://doi.org/10.1016/j.biochi.2011.07.032.
Sato, K., and Siomi, M.C. (2013). Piwi-interacting RNAs: biological functions and biogenesis. Essays In Biochemistry 54, 39-52. doi: 10.1042/bse0540039.
Schaefer, M., Pollex, T., Hanna, K., Tuorto, F., Meusburger, M., Helm, M., et al. (2010). RNA methylation by Dnmt2 protects transfer RNAs against stress-induced cleavage. Genes Dev 24(15), 1590-1595. doi: 10.1101/gad.586710.
Sharma, U., Conine, C.C., Shea, J.M., Boskovic, A., Derr, A.G., Bing, X.Y., et al. (2016). Biogenesis and function of tRNA fragments during sperm maturation and fertilization in mammals. Science 351(6271), 391-396. doi: doi:10.1126/science.aad6780.
Shigematsu, M., and Kirino, Y. (2015). tRNA-Derived Short Non-coding RNA as Interacting Partners of Argonaute Proteins. Gene Regulation and Systems Biology 9, GRSB.S29411. doi: 10.4137/grsb.S29411.
Siomi, M.C., Sato, K., Pezic, D., and Aravin, A.A. (2011). PIWI-interacting small RNAs: the vanguard of genome defence. Nat Rev Mol Cell Biol 12(4), 246-258.
Tang, Y., Gao, C.C., Gao, Y., Yang, Y., Shi, B., Yu, J.L., et al. (2020). OsNSUN2-Mediated 5-Methylcytosine mRNA Modification Enhances Rice Adaptation to High Temperature. Dev Cell 53(3), 272-286.e277. doi: 10.1016/j.devcel.2020.03.009.
Thompson, D.M., Lu, C., Green, P.J., and Parker, R. (2008). tRNA cleavage is a conserved response to oxidative stress in eukaryotes. RNA (New York, N.Y.) 14(10), 2095-2103. doi: 10.1261/rna.1232808.
Torrent, M., Chalancon, G., de Groot, N.S., Wuster, A., and Madan Babu, M. (2018). Cells alter their tRNA abundance to selectively regulate protein synthesis during stress conditions. Sci Signal 11(546). doi: 10.1126/scisignal.aat6409.
Vazquez, F., and Hohn, T. (2013). Biogenesis and Biological Activity of Secondary siRNAs in Plants. Scientifica 2013, 12. doi: 10.1155/2013/783253.
Vergara-Silva, F., Martínez-Castilla, L., and Alvarez-Buylla, E.R. (2000). MADS-box genes: development and evolution of plant body plans. Journal of Phycology 36(5), 803-812. doi: 10.1046/j.1529-8817.2000.03654.x.
Wang C, Chen W, Aili M, Zhu L, Chen Y (2023) tRNA-derived small RNAs in plant response to biotic and abiotic stresses. Front Plant Sci 14: 1131977
Wang, S., Li, P., Zhang, J., Qiu, D., and Guo, L. (2016). Generation of a high resolution map of sRNAs from Fusarium graminearum and analysis of responses to viral infection. Scientific Reports 6, 26151. doi: 10.1038/srep26151.
Wang, X.-Q., and Ran, J.-H. (2014). Evolution and biogeography of gymnosperms. Molecular Phylogenetics and Evolution 75, 24 - 40. doi: http://dx.doi.org/10.1016/j.ympev.2014.02.005.
Yakovlev, I.A., Carneros, E., Lee, Y., Olsen, J.E., and Fossdal, C.G. (2016). Transcriptional profiling of epigenetic regulators in somatic embryos during temperature induced formation of an epigenetic memory in Norway spruce. Planta 243(5), 1237-1249. doi: 10.1007/s00425-016-2484-8.
Yakovlev, I.A., and Fossdal, C.G. (2017). In Silico Analysis of Small RNAs Suggest Roles for Novel and Conserved miRNAs in the Formation of Epigenetic Memory in Somatic Embryos of Norway Spruce. Frontiers in Physiology 8(674), 1-17. doi: 10.3389/fphys.2017.00674.
Yang, L., Perrera, V., Saplaoura, E., Apelt, F., Bahin, M., Kramdi, A., et al. (2019). m5C Methylation Guides Systemic Transport of Messenger RNA over Graft Junctions in Plants. Current Biology 29(15), 2465-2476.e2465. doi: https://doi.org/10.1016/j.cub.2019.06.042.
Zeng, J., Gupta, V.K., Jiang, Y., Yang, B., Gong, L., and Zhu, H. (2019). Cross-Kingdom Small RNAs Among Animals, Plants and Microbes. Cells 8(4). doi: 10.3390/cells8040371.
Zhang, S., Sun, L., and Kragler, F. (2009). The Phloem-Delivered RNA Pool Contains Small Noncoding RNAs and Interferes with Translation. Plant Physiology 150(1), 378-387. doi: 10.1104/pp.108.134767.
Zheng, S., Li, J., Ma, L., Wang, H., Zhou, H., Ni, E., et al. (2019). OsAGO2 controls ROS production and the initiation of tapetal PCD by epigenetically regulating OsHXK1 expression in rice anthers. Proceedings of the National Academy of Sciences 116(15), 7549-7558. doi: doi:10.1073/pnas.1817675116.
Zhou, F., and Xu, Y. (2009). RepPop: a database for repetitive elements in Populus trichocarpa. BMC Genomics 10, 14. doi: 10.1186/1471-2164-10-14.
Zhu, L., Ow, D.W., and Dong, Z. (2018). Transfer RNA-derived small RNAs in plants. Science China Life Sciences 61(2), 155-161. doi: 10.1007/s11427-017-9167-5.

No competing interests reported.

Download PDF

Version 1

posted

You are reading this latest preprint version

Novel long small RNAs (lsRNAs) are evolutionary conserved in plants and differentially expressed during epitype inducing conditions in Norway spruce embryos

Status:

Version 1

Abstract

Figures

Introduction

Results

Characterization of small RNAs of gymnosperm species

Analysis of the possible origin of lsRNAs inP. abies embryos

lsRNAs are not plant analogs to animal piRNAs

lsRNAs may be derived from non-coding and repetitive element transcript

lsRNAs may be derived from tRNAs

lsRNAs expression profiles

Differential expression of lsRNAs under EpI temperatures in P. abies

Differential lsRNAs levels between tissues and species

Discussion

Materials and Methods

Plant material and RNA extraction

Small RNA library construction, sequencing, and bioinformatics

Identification of 28–37 nt lsRNAs and theirs mapping to theP. abies genome

Target prediction

DATA accession

Abbreviations

Declarations

Author Contribution

Acknowledgements

References

Additional Declarations

Supplementary Files

Status:

Version 1