Genome-Wide Discovery and Analysis of Phased Small Interfering RNAs in Chinese Sacred Lotus

Phased small interfering RNA (phasiRNA) generating loci (briefly as PHAS) in plants are a novel class of genes that are normally regulated by microRNAs (miRNAs). Similar to miRNAs, phasiRNAs encoded by PHAS play important regulatory roles by targeting protein coding transcripts in plant species. We performed a genome-wide discovery of PHAS loci in Chinese sacred lotus and identified a total of 106 PHAS loci. Of these, 47 loci generate 21 nucleotide (nt) phasiRNAs and 59 loci generate 24 nt phasiRNAs, respectively. We have also identified a new putative TAS3 and a putative TAS4 loci in the lotus genome. Our results show that some of the nucleotide-binding, leucine-rich repeat (NB-LRR) disease resistance proteins and MYB transcription factors potentially generate phasiRNAs. Furthermore, our results suggest that some large subunit (LSU) rRNAs can derive putative phasiRNAs, which is potentially resulted from crosstalk between small RNA biogenesis pathways that are employed to process rRNAs and PHAS loci, respectively. Some of the identified phasiRNAs have putative trans-targets with less than 4 mismatches, suggesting that the identified PHAS are involved in many different pathways. Finally, the discovery of 24 nt PHAS in lotus suggests that there are 24 nt PHAS in dicots.

Introduction of 47,573,025 reads) respectively [30]. These two small RNA profiles had been deposited into the NCBI GEO database under the series accession number GSE62217. The genome and cDNA sequences of Chinese sacred lotus (Nelumbo nucifera Gaertn.) were downloaded from NCBI GenBank [28]. The sequences of TAS3, TAS4 loci and their derived tasiRNAs were downloaded from the tasiRNAdb [31].

Computational steps
The unique sequences in the small RNA libraries were mapped to the genome and cDNA sequences of Chinese sacred lotus with SOAP2 [32]. A self-written program was used to scan the genome and cDNA sequences using a window of 210 nt or 240 nt (ten 21 nt or 24 nt) respectively. A two-nucleotide positive offset was used to calculate the positions of siRNAs on the anti-sense strand because the existence of two-nucleotide over-hang at the 39-end of siRNA duplex [11,18,20,21]. Then a P-value was calculated for each of the windows using a modified version of methods in [20], where n was the number of unique 21 nt (or 24 nt) sRNAs mapped within a window, k was the number of phased unique 21 nt (or 24 nt) sRNAs within the window, and m was the number of phases. Similar to previous work [33], m was set to 10 in this study. And a phase score was calculated for each position of the genome and cDNA sequences using the method in [34]. For a window started at a position with more than three phased unique sRNAs, i.e., when k §3, where P i was the number of phased reads at the ith phase from the position, U i was the number of non-phased reads at the ith phase from the position, and m was the number of phases in the window, and k was the number of unique phased siRNAs in the window. m was 10 in this study. The window with a P-value less than 0.05 was extended 100 base pairs at both 59-and 39-ends, then the overlapped windows were merged. The P-values of the merged windows were used to calculate the false positive rates using the method in [35]. The merged windows with a maximal phase scores of larger than predetermined threshold and multiple test corrected P-values of smaller than 0.05 were reported as PHAS loci. The predicted PHAS were named with its chromosome (or scaffold) and a unique serial number for each chromosome. The neighboring PHAS loci were predicted as PHAS clusters if the distances between individual PHAS loci were smaller than 2,000 base pairs. The phased siRNAs of the predicted PHAS loci were reported as phasiRNAs. The phasiRNAs of a PHAS loci were named by adding siR and a serial number to the name of the PHAS loci.
The miRNA binding sites on PHAS and the targets of predicted phasiRNAs were predicted with the HitSensor algorithm [29]. For 21 nt/22 nt miRNAs and phasiRNAs, targets with less than 4 mismatches were kept for analysis. For 24 nt miRNAs and phasiRNAs, targets with less than 6 mismatches were maintained for analysis.
We combined the annotation of genes of Chinese sacred lotus in [28] with alignment results of predicted PHAS sequences to the NCBI Nucleotide Collection (nr/nt) database and the TIGR Repeat database [36].
The phylogenetic trees of the predicted TAS3, TAS4 loci and their derived tasiRNAs were constructed with the Bootstrap Neighboring-Joining algorithm implemented in ClustalX (version 2.1) [37] and visualized with TreeView [38].

nt PHAS loci in Chinese sacred locus
We predicted PHAS loci by using the alignments of small RNA sequencing libraries to the genome and cDNA database respectively. The predicted PHAS loci were combined and merged if necessary. As listed in Table 1, we totally identified 16 and 7 loci corresponding to 21 nt and 24 nt PHAS loci, respectively, when using a phase score threshold of 10 (Pv0:01, as shown in Figure S1a). After relaxing the phase score threshold to 5 (Pv0:03, see Figure S1a), we identified 31 additional 21 nt PHAS loci (in Table S1), three of which are shown in Table 1.
In addition to the annotation from [28], we aligned the predicted PHAS to the NCBI Nucleotide Collection (nr/nt) database and the TIGR Plant Repeat database to refine putative annotation of the predicted PHAS loci (details are given in Table  S1). Furthermore, as miRNAs are critical in generation of phasiRNAs, we predicted miRNA complementary sites on these PHAS loci on both strands. Based on the miRNA complementary sites, three loci are classified as TAS3 and one locus is classified as TAS4, which will be discussed below. As shown in Figure 1a, the largest category, 50%, of 21 nt PHAS loci are rRNA or repeats. There are eight (17%) 21 nt PHAS overlapped to protein coding genes. Four (9%) 21 nt PHAS loci are TAS3 and TAS4 loci, and two (4%) 21 nt PHAS are pre-miRNA159 loci.
Previous studies found that miR1507, miR1509 and miR2118 families trigger several PHAS loci, especially the NB-LRR disease resistance genes that possess one complementary site to miR2118a, a 22 nt miRNA [21]. Our analysis also predicted miR2118 complementary sites on three predicted PHAS loci, scaffold_87_1 (Figure 2a), PHAS_sf_122_1, and scaffold_170_1 (Figure 2b), from which phased siRNAs are originated, suggesting these loci are NB-LRR family members. Furthermore, after aligning their sequences to the NCBI nr/nt database, PHAS_sf_122_1 and scaffold_87_1 are putative disease resistance proteins (with E-values of 8|10 {33 and 0.002, respectively, see Table S1). Another locus, scaffold_8_1 overlapping with a putative disease resistance RPP13-like protein 1 (NNU_004711-RA) also has a miR2118 site in the upstream region ( Figure 2c). However, there is no miR2118 site around scaffold_131_1, which overlaps with a putative disease resistance protein RGA3 (see Table 1). miR828 targeted MYB genes can generate phasiRNAs in apple [26]. Three PHAS loci scaffold_90_1, PHAS_sf_173_1 and PHAS_sf_21_1 possess complementary sites to miR828 in sacred lotus (see Figure 2d, Figure S2a and Figure  S2b). scaffold_90_1 and PHAS_sf_21_1 overlap with two MYB transcription factors (NNU_018790 and NNU_008785 respectively). PHAS_sf_173_1, shown in Figure S2a and S2c, is a potential MYB transcription factor after aligning its sequence to the NCBI nr/nt database (with an E-value of 4|10 {25 , see Table S1).
As listed in Table S1, the 21 nt PHAS loci also form three additional PHAS clusters in addition to scaffold_326_1/2 in Table 1 and Figure 5a to b.
We do not find TAS1 and TAS2 loci in Chinese sacred locus, probably because these two TAS families are less conserved than TAS3 [8].
A substantial number of predicted 21 nt PHAS loci, 20%, are located in regions of unknown genes or intergenic regions, as shown in Figure 1a. They also did not matched to known genes in NCBI nr/nt database and known repeat elements in TIGR Repeat database. Further research is needed to clarify whether these are also PHAS loci or not in sacred lotus.

Discovery of yet another TAS3 locus in Chinese sacred locus
Among the sixteen 21 nt PHAS loci, three are TAS3 loci, as shown in Figure 3a. Two of them were identified in our previous study [30]. The typical 59 and 39 miR390 complementary sites around the phased siRNA regions of these TAS3 loci are given in Figure 3b and c, respectively. TAS3c is unique since (i) its phase starts from two positions (position 10 and 12) of the 39 miR390 complementary site; (ii) TAS3 only encode one conserved tasiARF (tasiRNAs that target ARF family members [14]), as shown in Figure 3a. The conserved tasiRNA originated from position 12 has a two nucleotide shift, named as TAS3c_D6-2(+) (in Figure 3d). The most abundant siRNA from the whole TAS3 gene is also in phase with position 12, as shown in Figure 3e. The targets of the conserved tasiRNAs include 8 ARF family members. One of these ARF genes is targeted by TAS3c siRNAs with much better complementarities, especially by TAS3c_D6-2(+) (as shown in Figure 3g). These results suggest that the two nucleotide shift in the phase of TAS3 may be beneficial to its targeting to the additional ARF member (NNU_003220-RA). TAS3c may be a shorter variant of TAS3 which only encodes one tasiRNA and both the 59 and 39 miR390 sites of this variant are cleavable as mentioned in [8,26]. As shown in Figure 3b, the 59 miR390 site on TAS3c has more matched nucleotides at positions 10, 11 and 19 than TAS3a and TAS3b, suggesting this site is cleavable. We performed conservation analysis for the TAS3 loci and their derived tasiRNAs to show their relations to TAS3 genes in other species as shown in Figure S3. The identified TAS3a and TAS3b have a close relation, but TAS3c is far from the cluster of TAS3a and TAS3b ( Figure S3a), which is in accordance with the unique features of TAS3c discussed above. TAS3c derived tasiRNA is also distant from TAS3a and TAS3b derived tasiRNAs ( Figure S3b and S3c).

Discovery of TAS4 locus
As shown in Figure 4a and b, one 21 nt PHAS locus (sf_39_1) has a typical miR828 complementary site of TAS4 [19] at the 59 side of the phased region. This locus is annotated as a protein of unknown function, NNU_012673-RA [28], with two exons. Together with these results, conservation analysis suggests that this locus is a conserved TAS4 gene (see Figure S4a). The transcription region of this locus is longer than the exon regions of NNU_012673-RA because there are some phased and non-phased siRNAs beyond the 39 end of the second exon of NNU_012673-RA, as shown in Figure 4a and c. One of the siRNAs, sf_39_1_siR4, derived from sf_39_1 is highly conserved to TAS4-siR81(2) reported in Arabidopsis thaliana [19] and other species [31] (see Figure S4b and S4c). Thus, it  PhasiRNAs in Chinese Sacred Lotus is named as TAS4-siR81(2) too. It is interesting that this locus also produces phased 24 nt siRNAs as shown in Figure 4b and d. This suggests that TAS4 transcripts of sacred lotus can be processed by both DCL4 and DCL5 (also named as DCL3b) [8,12,13] to produce both 21 nt and 24 nt phasiRNAs, respectively.

LSU-rRNA derived 21 nt phased siRNAs
Five 21 nt PHAS loci are annotated as LSU-rRNA ( Table 1) and two of them overlap with each other (Figure 5a and b). It was observed that rRNAs could generate small RNAs [39][40][41] which is dependent on two RNA-dependent RNA polymerases, RDR2 and RDR6, and DCL2/3/4 [40]. Since phasiRNA biogenesis also requires RDR6 and DCL4 [8], the pathways for generating small RNAs from PHAS and rRNA loci share some key protein components with each other. Our results suggest that some LSU-rRNAs are processed through the phasiRNA biogenesis pathway to produce 21 nt phasiRNAs.

nt PHAS loci
We identified 7 PHAS loci that generate 24 nt long phasiRNAs using a phase score threshold of 10 (Pv0:001, see Figure S1b), as shown in Table 1 and Figure 6. Fifty two additional 24 nt PHAS loci are predicted when using 5 as the threshold of phase score (Pv0:01, see Figure S1b) as listed in Table S1. The 24 nt PHAS loci form 11 PHAS clusters as shown in Table S1. Similarly, we also aligned the 24 nt PHAS to the NCBI nr/nt database and the TIGR Plant Repeat database to obtain putative annotation of these PHAS loci (details are given in Table S1). As shown in Figure 1b, an overwhelming portion, 88%, of 24 nt PHAS is unknown genes or elements because they have no matches to the known genes and known repeat elements. Compared to 21 nt PHAS, more researches are needed to clarify these PHAS with unclear annotation. Five (8%) 24 nt PHAS are matched to rRNA or repeats. One locus, scaffold_39_1, is TAS4 as mentioned above. Unlike the 21 nt phasiRNAs being initiated by miR828, no significant miRNA complementary sites are predicted to originate 24 nt phasiRNAs at scaffold_39_1, as well as other 24 nt PHAS loci, although existing studies indicate that miR2275 triggers the generation of 24 nt phasiRNAs in rice [12,13]. Only one 24 nt PHAS is aligned to a hypothetical protein (with an E-value of 1|10 {6 ).
24 nt PHAS loci have mainly been reported in monocots [8,12,13]. Together with the 24 nt PHAS loci reported in grapevine [33], our results suggest that dicots may also have 24 nt PHAS loci, although their biogenesis processes are still to be clarified.

nt phasiRNAs and their putative targets
The first sixteen 21 nt PHAS loci (Table 1) totally derived 224 phasiRNAs, as shown in Table S2. The abundance of these phasiRNAs is similar in the flower and leaf small RNA libraries, suggesting that these phasiRNAs have similar functions in these two tissues. The most differently expressed phasiRNA with more than 100   PhasiRNAs in Chinese Sacred Lotus RPTM in both libraries is an siRNA derived from an LSU-rRNA (scaffold_326_3), as shown in Figure 5d. Its abundance has dropped from 2600 RPTM in flower to 1300 RPTM in leaf. It needs further investigation to clarify whether such differential abundance has any significance or not.
Several phasiRNAs have perfect alignments to the transcripts on their antisense strands, which are regarded as cis-targets. Targets from other loci are regarded as trans-targets. Some of the trans-targets with the least numbers of mismatches in their complementary sites are listed in Table 2 except the well-known ARF family members targeted by TAS3 derived siRNAs (which have been shown in Figure 3f). As shown in Table 1 and Figures S4b/S4c, a conserved TAS4 derived phasiRNA (sf_39_1_siR4 or TAS4-siR81 (2)) targets an MYB transcription factor (NNU_018790). And in Table S3, sf_39_1 derived phasiRNAs target five other MYB transcription factors, two of which are targeted by sf_39_1_siR4. In comparison, MYB transcription factors are also targeted by TAS4 derived siRNAs in Arabidopsis [19,26,27,42]. PhasiRNAs derived from putative NB-LRR genes scaffold_170_1 and scaffold_87_1 potentially target as many as 49 and 44 different members of the NB-LRR disease resistance protein family (see Table 2 and Table  S3), respectively. PhasiRNAs derived from another putative NB-LRR gene (NNU_004711-RA), scaffold_8_1, only has one putative trans-target in the NB-LRR disease resistance protein family (NNU_025710-RA). Similarly, phasiRNAs derived from a putative MYB gene (scaffold_90_1) potentially target 9 other MYB family members in trans (see Table S3). These results are consistent with those reported in other species [8,21,26], and similar to siRNAs derived from PPR genes [11,43].
As shown in Table 2, phasiRNAs from the predicted PHAS loci also target many other genes with or without known functions. The biological relevance of these predicted targets still needs additional studies.

nt phasiRNAs and their putative targets
The seven 24 nt PHAS loci in Table 1 totally generate sixty six 24 nt long phasiRNAs (in Table S2). Similar to 21 nt phasiRNAs, these 24 nt phasiRNAs have similar abundance in the flower and leaf sequencing libraries (see Table S2). Some of the predicted trans-targets of these 24 nt phasiRNAs are given in Table 2 and a full list is given in Table S3. As shown in Table 2, sfd_39_1_siR3 potentially targets an Myb-related protein 315 (NNU_003001) with only 2.5 mismatches in a 24 nt complementary site. As mentioned earlier, TAS4a-siR81(2) (i.e., sf_39_1_siR4) potentially targets an MYB transcription factor (NNU_018790), which is a conserved mechanism [19,26,27,42]. sfd_39_1_siR3 is 24nt long and 3 nt downstream of TAS4a-siR81(2), which is potentially beneficial to its  complementarity to NNU_003001 because the 21 nt sequence immediately downstream of TAS4a-siR81(2) has 5.5 mismatches to the same complementary site on NNU_003001. These results suggest that the potential 24 nt phasiRNAs derived from TAS4a (sf_39_1), as shown in Figure 4, might be produced as alternative small RNA guiders to repress a set of MYB transcription factors other than those repressed by the conserved 21 nt TAS4a-siR81 (2). Another 24 nt phasiRNA sfd_88_1_siR7 has 6 targets with less than 3 mismatches (see Table 2), suggesting it is biologically relevant. However, more studies are necessary to verify whether these 24 nt phasiRNAs are really functional and have trans-targets.

Conclusion
Through genome-wide analysis of two small RNA sequencing libraries, we totally identified 23 putative PHAS loci from Chinese sacred lotus, and 83 additional PHAS loci were further identified with a smaller phase score threshold, including a new putative TAS3 and a putative TAS4 loci. Our results show that the predicted TAS4 loci derives both 21 nt and 24 nt phasiRNAs, suggesting that both DCL4 and DCL5 are involved in the biogenesis of phasiRNAs from TAS4. Several PHAS loci are from NB-LRR and MYB loci, which is consistent with existing results in other species. Existing studies reported that small RNAs are derived from PHAS and rRNA through different pathways with shared protein components, RDR6 and DCL4. In this study, we found that some LSU-rRNA could generate phasiRNAs, suggesting that some rRNAs are processed by the PHAS siRNA biogenesis pathway. However their biogenesis still needs further studies. The identification of 24 nt PHAS loci in Chinese sacred lotus, as well as similar results reported in grapevine [33], suggests that dicots may encode 24 nt PHAS.   The phylogenetic tree of TAS3 derived tasiRNAs. (c) The multiple sequence alignment of TAS3 derived tasiRNAs generated with ClustalX (version 2.1) [37]. The sequences of TAS3 loci and derived tasiRNAs were used to construct the phylogenetic trees with the Bootstrap Neighbor-Joining algorithm implemented in ClustalX (version 2.1). Then, the trees were visualized with TreeView [38]. The numbers in the trees are bootstrap values greater than 500 (50%). The lower case letters at the beginnings of the names of TAS3 and tasiRNAs stand for the species, i.e., at (Arabidopsis thaliana), bn (Brassica napus), cl (Cunninghamia lanceolata), cm (Cucumis melo), gm (Glycine max), lj (Lotus japonicus), md (Malus domestica), mt (Medicago truncatula), nn (Nelumbo nucifera (Gaertn)), nt (Nicotiana tabacum), oe (Olea europaea), os (Oryza sativa), ppa (Physcomitrella patens), ppe (Prunus persica), sl (Solanum lycopersicum), ta (Triticum aestivum), vv (Vitis vinifera), and zm (Zea mays). nnTAS3a, nnTAS3b and nnTAS3c are scaf-fold_106_1, scaffold_65_1, and scaffold_10_1, respectively. The tasiRNAs of TAS3a/b/c in Chinese sacred lotus are given in Figure 3d. The phylogenetic tree of TAS4 derived tasiRNAs. (c) The multiple sequence alignment of TAS4 derived tasiRNAs. The legend are the same as those of Figure  S3. The lower case letters at the beginnings of the names of TAS3 and tasiRNAs stand for the species, i.e., at (Arabidopsis thaliana), md (Malus domestica), nn (Nelumbo nucifera (Gaertn)), ppe (Prunus persica), and vv (Vitis vinifera). nnTAS4 is sf_39_1 and nnTAS4-siR81(2) is sf_39_1_siR4. doi:10.1371/journal.pone.0113790.s004 (TIF )  Table S1. The predicted PHAS loci in Chinese sacred lotus. The loci of ID started with scaffold_ and sf_ are predicted with a phase score threshold of 10, and loci of id started with PHAS_ are predicted with a phase score threshold of 5. The loci from the same PHAS cluster are marked with the same cluster ID in the Cluster column. doi:10.1371/journal.pone.0113790.s005 (XLSX)