Identification and characterization of miRNAome and target genes in Pseudostellaria heterophylla

miRNAs play a crucial role in the development and growth of plants by inhibiting the function of targeted genes at the post-transcription level. However, no miRNAs in Pseudostellaria heterophylla have been reported and their function in the morphogenesis of organs is still unclear. In this study, a total of 159 conserved miRNAs (belonging to 64 families) and 303 level miRNAs were identified from P. heterophylla. Some of them showed specifically up or down-regulated expression in different tissues and numbers of unigenes involved in Plant-pathogen interaction and MAPK signaling pathway-plant were targeted. The significant negative correlation of expression profiles between 30 miRNAs and their target genes (37 unigenes) was observed, respectively. Further, a large number of genes involved with signal transduction of auxin, zeatin, abscisic acid and, jasmonic acid were targeted. Predicated targets of two miRNAs were validated by 5′RLM-RACE, respectively. A large number of mRNAs from four pathogens were targeted by miRNAs from P. heterophylla and some of them were targeted by miR414. In summary, we reported a population of miRNAs from four different vegetative tissues of P. heterophylla by high throughput sequencing, which was analyzed by combining with the constructed transcriptome. These results may help to explain the function of miRNAs in the morphogenesis of organs and defense of pathogens, and may provide theoretical basis for breeding and genetic improvement of P. heterophylla.


Introduction
MicroRNA (miRNA) is a kind of endogenous single-stranded non-coding small RNA [1], which inhibits the function of targeted genes at the post-transcription level by degrading or inhibiting the translation [2]. Since the first miRNA (lin 4) in C.elegans was identified [3], a large number of miRNAs have been discovered from various species. Recent evidence has indicated that miRNA can not only impact the growth and development of plants, and the differentiation and morphogenesis of organs by directly regulating the target genes [4], but also participate in the signal transduction of plants [5]. miRNAs can also improve the tolerance to environmental stress by regulating the transcription factor in plants [6], such as drought, low temperature, high salt, and low nitrogen. Furthermore, miRNAs may improve the resistance to biological stress of plants by inhibiting the gene expression of pathogenic microorganisms [7,8]. tissues with the base average error rate below 0.02%. Trinity program was used for de novo assembling clean data. 138,888 unigenes and 207,390 transcripts were obtained. The length of unigenes ranged from 201 to 14,793 bp and the N50 length was 1,862. 71,215 unigenes (51.28%) and 121,508 transcripts (58.59%) were annotated based on the public databases.
In the small RNA library, a total of 180 million raw reads and 139 million clean reads were yielded from these 12 libraries (Table 1). These raw data were filtered through several criteria to identify conserved and specific miRNAs. The length distribution of the small RNAs was analyzed and shown in Fig 1. The majority of reads were 21 to 24 nt in length. The 24-nt RNA is the most abundant class in the sRNA library and comprised 35.7% of root bark, 32.7% of root xylem, 31.5% of leaf and 31.8% of stem, respectively. Similar results were reported in potato [24] and cucumber [30], but in contrast to wheat [31].

Identification of novel candidate miRNAs in P. heterophylla
The non-conserved small RNA reads were mapped to unigenes in the transcriptome of P. heterophylla and a total of 303 novel miRNA candidates were identified (S2 Table). A standard stem-loop hairpin secondary structure (SS) was found in all miRNA precursors (not shown). The secondary structure of precursors of novel miRNAs was predicted by using the Mfold web server [32] and a few were shown in Fig 3. These miRNA precursors had folding free energies ranging from -48 to -18 kcal/mol (average, -27.28 kcal/mol). The length of the predicted precursor of these novel miRNAs was 58-101 nt. The MFEI values of all miRNAs were above 0.85 and these sequences were considered as miRNAs [33]. 4 novel miRNAs had high expression level with TPM >1 000 and the highest read abundance (2 662 normalized reads) was detected for Phn-m0190. 73 novel miRNAs with total TPM of less than 10 were found, such as Phn-m0032, Phn-m0045, Phn-m0050, Phn-m0058, and Phn-m0076.

Prediction of miRNA target genes
To explore the biological functions of identified miRNAs, the Target Finder program was used for searching the complimentary mRNA sequences from the transcriptome of P. heterophylla above. Among the conserved targets, numbers of unigenes coding transcription factors were identified (S3 Table), such as GAMYB-like (targeted by miR159), APETALA2 (targeted by    Table). 162 potential target unigenes coding transcription factors were also observed, including MYB (44 unigenes), bHLH (34 unigenes), WRKY (15 unigenes), and MADS (12 unigenes). In addition, 34 potential target genes were founded to be involved in the synthesis and biotransformation of NADPH. Multiple potential regulation patterns were discovered in several unigenes (S5 Table). 8 binding sites of Phn-m051-5p and Phn-m052-5p were identified in the mRNA sequence of ubiquitin-conjugating enzyme E2 24. Transmembrane protein 97 has 11 potential regulation sites targeted by miRNA399 families. Three HD-ZIPIII genes (ATHB8, ATHB15 and REVO-LUTA) also have 9 potential regulation sites targeted by miRNA166 families. Furthermore, several unigenes were only targeted by novel miRNAs, such as multicopper oxidase LPR2, Lascorbate oxidase homolog, ribosomal RNA-processing protein 7 homolog A and pyrrolidone-carboxylate peptidase. In addition, plenty of non-coding unigenes were targeted (S3 and S4 Tables), including 609 non-coding unigenes targeted by 89 conserved miRNAs and 3673 non-coding unigenes targeted by 274 novel miRNAs. These results indicated that miRNAs may also participate in the regulation of non-coding unigenes.

GO and KEGG analysis
GO analysis was performed to explore the potential regulation of predicted targeted unigenes of miRNAs in P. heterophylla. Based on their functional annotations, the targeted genes were divided into 122 categories (S6 Table). In particular, 6 categories have higher enrichment, such as protein kinase activity, protein serine/threonine kinase activity, phosphotransferase activity, alcohol group as acceptor, calcium ion binding, kinase activity, and response to acid chemical (Fig 4). KEGG pathway enrichment analysis showed that the targeted genes of the miRNAs were mainly involved in 8 pathways, including Plant-pathogen interaction, MAPK signaling pathway-plant, plant hormone signal transduction, ubiquitin-mediated proteolysis, phosphatidylinositol signaling system, phenylpropanoid biosynthesis, inositol phosphate metabolism and diterpenoid biosynthesis (S7 Table and Fig 5). Numbers of targeted genes only involving in Plant-pathogen interaction were targeted, including 9 calcium-dependent protein kinases, 4 pto-interacting protein, 3 disease resistance proteins and 2 calcium-binding proteins, 2 cyclic nucleotide-gated ion channels. Several targeted genes were annotated only in MAPK signaling pathway-plant pathways, such as mitogen-activated protein kinase (MPK1, 3,4,8) and reversion-to-ethylene sensitivity (TMEM222). However, partial targeted genes participated in both Plant-pathogen interaction and MAPK signaling pathway-plant pathway, including 6 calmodulins, 6 respiratory burst oxidase homolog proteins, and 3 WRKY transcription factor. These results suggested that miRNAs in P. heterophylla may play an important role in the regulation of cell signal transduction.

Role of miRNAs in terrestrial and underground tissues
For a better understanding of the role of miRNAs in the morphogenesis of organs, tissue-special miRNAs were identified according to their expression profiles (S1 and S2 Tables), and correlation analyses were performed by combining with the transcriptome (S8 Table). 6 conserved miRNAs and 9 novel miRNAs was found to be up-regulated in underground tissues (root bark and xylem) and down-regulated in terrestrial tissues (stem and leaf). 115 potential target genes were identified (S9 Table), including 33 coding unigenes and 82 non-coding unigenes. 16 unigenes targeted by Phn-m0170-3p and 6 genes targeted by Phn-m0173-5p were identified. CRK7, PUB6, and BBX19 showed significantly negative correlation with the expression of Phn-m0173-5p. SODC may be negatively regulated by miR398a-3p. 12 conserved miR-NAs and 13 novel miRNAs showed abundant expression in underground tissues while being restricted to terrestrial tissues. 107 potential target genes were identified, including 37 coding mRNA and 70 non-coding mRNA. Two conserved miRNAs (miRNA169a and miRNA169a) and a conserved miRNA (Phn-m0046-5p) were negatively correlated with NF-YA5 and PP2C08 respectively (Fig 6), which suggested their gene-silencing function. 9 conserved miRNAs and 2 novel miRNAs were significantly down-regulated in leaf (S10 Table). 36 potential target genes were identified, including 13 coding mRNA and 23 non-coding mRNA. TCP showed significantly negative correlation with 3 miRNA319 members, miR319a, miR319a-3p.2-3p, and miR319i (Fig 7). However, 19 conserved miRNAs and 2 novel miRNAs show abundant expression in leaf. 88 potential target genes were identified, including 29 coding unigenes and 59 non-coding unigenes. None of miRNAs showed negative correlation with their target unigenes, while some unigenes have positive correlation with miRNAs, such as GAMYB (targeted by miR159a, miR159c, and miR159d), EAF6 (targeted by miRNA156a) and receptor-like protein kinase At2g23200 (targeted by miR172b).

Role of miRNAs in the development of stem
3 conserved miRNAs and 12 novel miRNAs have abundance transcripts in stem (S11 Table). 239 potential unigenes targeted by 15 miRNAs were identified, including 75 coding mRNA and 164 non-coding mRNA. Interestingly, a few miRNAs have several potential target genes, such as 36 coding mRNA targeted by Phn-m0082-3p, 10 coding mRNA targeted by Phn-m0077-5p, and 7 coding mRNA targeted by Phn-m0024-3p. miR156l-5p and Phn-m0224-3p showed low expression in stem and high expression in underground tissues and leaf. The miR156l-5p showed significant negative correlation with SPL5 (Fig 8).

Role of miRNAs in the development of tuberous root
To explore the function of miRNAs in the development of tuberous root, tissue-special miR-NAs in xylem and bark have been identified. 12 novel miRNAs were up-regulated in the root xylem (S12 Table). 96 potential target genes were identified, including 36 coding unigenes and 60 non-coding unigenes. Though 24 coding unigenes targeted by Phn-m0101-5p and 6 coding unigenes targeted by Phn-m0279-3p were observed, none of miRNAs showed negative correlation with their potential target unigenes (Fig 9). However, 16 conserved miRNAs and 5 novel miRNAs showed low expression in the root xylem. 77 potential target genes were identified, including 38 coding unigenes and 39 non-coding unigenes. Three HD-ZIP IIIs showed extreme negative correlation with the miR166 family. It seems that SCL6, BIR2, ARFR, and WDL5, were negatively regulated by the miR171 family, Phn-m0162-5p, miR160d-5p and miR166g-5p in xylem of the tuberous root, respectively. In addition, 7 conserved miRNAs were especially up-regulated in root bark (S13 Table), including 3 members of the miR160 family (miR160f-5p, miR160g and miR160h), 2 members of miR390 family (miR390-3p and miR390d-3p), miR169w and miR171a. However, none of miRNAs showed negative correlation with their potential target genes.

Role of miRNAs in hormone signal transduction
The hormone plays a curial role in the morphogenesis of plant organs. In this study, we analyzed the regulation function of miRNAs on hormone signal transduction by searching their target genes. A large number of unigenes in hormone signal transduction pathways were targeted by identified miRNAs in P. heterophylla, mainly involving signal transduction of auxin, zeatin, abscisic acid, and jasmonic acid (S14 Table). Five gene families in the auxin signal transduction pathway may strongly be regulated by miRNAs from P. heterophylla, including TIR1, Aux/IAA, and ARF. SAUR was targeted by 10 miRNAs, including 2 known miRNAs and 8 novel miRNAs. In CTK signaling transduction, there were several binding sites of miR-NAs in the sequences of AHP, B-ARR and A-ARR. Interestingly, PYL, PP2C and, SnRK2 in the ABA signal transduction pathway may strongly be regulated by miRNAs, mainly composed of novel miRNAs. We also found that a few novel miRNAs may regulate the expression of COI1 and JAZ in jasmonic acid signaling transduction. In addition, gibberellin receptor GID2 in gibberellin signal transduction was targeted by Phn-m0174-5p. BSK in the brassinosteroid signaling transduction pathway was targeted by 3 known miRNAs and 5 novel miRNAs. A few unigenes involving with the ethylene signal transduction were identified, such as MPK6 (targeted by Phn-m0087-5p, Phn-m0154-3p and Phn-m0281-5p), EBF1/2 (targeted by Phn-m0033-5p), and EIN3 (targeted by mir166e-5p and Phn-m0122-3p).

Role of miRNAs in the interaction between plants and pathogens
For a better understanding of the role of miRNAs in the interaction between plants and pathogens, mRNA sequences of 4 microorganisms (Fusarium oxysporum, Leptosphaeria maculans, Microbotryum violaceum, Rhodosporidium toruloides) were downloaded from Genbank and used for searching the target genes of miRNAs from P. heterophylla (Fig 10 and S15 Table). 39

PLOS ONE
miRNAome from Pseudostellaria heterophylla genes from F. oxysporum were targeted by 30 miRNAs of P. heterophylla, including 15 known miRNAs and 15 novel miRNAs. 6 genes were targeted by miR156h, followed by miR414 (4 genes) and miR171a (3 genes). Furthermore, some important genes of F. oxysporum may be regulated by miRNAs from P. heterophylla, such as pectin lyase F (targeted by Phn-m0104-3p), 30S ribosomal protein S17-like protein (targeted by miR171e-3p) and myo-inositol transport protein ITR1 (targeted by miR171a). Interestingly, there were 37, 43, and 27 genes in L. maculans, M. violaceum, and R. toruloides targeted by miR414, respectively. These results suggested that miR414 may play an important role in the interaction between P. heterophylla and its pathogens.

miRNA-induced cleavage of predicted targets
The cleavage site of predicted targets of a conserved miRNA and a novel miRNA in P. heterophylla was further verified by 5 0 RLM-RACE (Fig 11), respectively. HD-ZIP REVOLUTA and delta (24)-sterol reductase were targeted by miR166k and Phn-m0121-5p, respectively. All the fragments of 5 0 RLM-RACE were analyzed on agarose gel, purified, cloned, and sequenced. Delta (24)-sterol reductase (TRINITY_DN86257_c0_g3) was found to be cleaved between 5 and 6. Interestingly, a new cleavage site of HD-ZIP REVOLUTA (TRINITY_DN91558_c0_g3) was found between G: G bases. This result suggested HD-ZIP REVOLUTA may be coordinately regulated by miR166k and Phn-m0140-3p/Phn-m0141-3p.

Discussion
In earlier studies, a large numbers of miRNAs were identified from various plants through computational and direct cloning approaches from different tissues of certain periods [19][20][21][22][23], and plenty of them showed specifical expression similar to the coding genes. For example, 89 conserved miRNAs and 259 potato-specific miRNAs were discovered in root, stem, leaf, and tuber of potato [24] and a few of them demonstrated tuberization stage-specific expressions. In this study, 163 conserved miRNAs and 303 level miRNAs were identified and some of them also showed specifically up or down-regulated expression in leaf, stem, xylem and bark of the tuberous root. These results may help to reveal the role of miRNAs in the morphogenesis of different organs.
In fact, the morphogenesis of organs depends on expression pattern of genes in plants [6][7][8]. In this study, the expression patterns of miRNAs in different tissues were analyzed and these roles in organ morphogenesis have been revealed. Two members of miRNA172 were significantly up-regulated in the terrestrial tissues, and might be associated with curl symptoms of leaf and affect the flowering time [18]. A few members of miRNA164 family and miRNA169 family were significantly up-regulated in the underground tissues. miRNA164 may reduce the formation of lateral roots by down-regulating the expression level of NAC1 [34]. miRNA169 may affect the number and size of certain cells in root meristem and the initiation of lateral roots [35], and participate in the response of roots to environmental stresses, such as drought and salt stress [36].
In reality, the formation of tuberous root mainly depends on the formation of secondary xylem from phloem and the expansion of parenchyma cells [37]. Therefore, the development regulation of xylem plays a crucial role in the formation of the tuberous root. In this study, several members of miR171, miR160, miR399, and miR166 families were significantly down-regulated in the xylem of tuberous roots of P. heterophylla. miR171 may affect the growth and development of tuberous roots in P. heterophylla by regulating three target genes of the scarlike family, namely SCL6-II, SCL6-III, and SCL6-IV [38,39], controlling the radial growth of roots [40] and hormone signal [41]. miR166 may affect the number of lateral roots and symbiotic nodules, and induce the ectopic development of vascular bundles in these transgenic roots [42]. In addition, 12 novel miRNAs were significantly up-regulated in the xylem of tuberous root, which suggested its important role in the development of root xylem.
miRNA not only directly regulates the development of plants through its target genes, but also indirectly affects the morphogenesis of organs by regulating the expression of genes involved in hormone signal transduction. Previous studies have shown that ABA, CTK, and JA increased gradually with the development of tuberous roots in P. heterophylla [43], and which involved in the regulation of cell proliferation, cell elongation, cell expansion [44], formation of microtubule microfilaments [45] and accumulation of assimilates in tuberous roots [46]. In this study, several key genes in hormone signal transduction were targeted by miRNAs from P. heterophylla, such as miR393 and miR172. Interestingly, most of unigenes involved in signal transduction of ABA, CTK, and JA were targeted by novel miRNAs, which suggested a different regulation mechanism of miRNAs from P. heterophylla.
Many plants have evolved perfect defense mechanisms, including the miRNA system. miR393 was firstly discovered and regulated by necrotroph pathogens [47]. A variety of new miRNAs could be detected in plants infected by the rice stripe virus [48,49]. In our previous studies, five kinds of pathogens were identified in P. heterophylla, in which M. violaceum, L. maculans, and P. fluorescens are typical harmful bacteria [29]. In this study, a large number of mRNAs from four kinds of pathogens were targeted by miRNAs from P. heterophylla. Several genes were targeted by miR414, while none of unigenes in P. heterophylla were targeted. These results suggested that miR414 may provide a new way to overcome the disease resistance and continuous cropping obstacles of plants. However, miR414 showed extremely low expression in the tissue of P. heterophylla and it is necessary to clarify its function in the defense of pathogens.

Conclusions
In this study, we reported a population of P. heterophylla miRNAs from four different vegetative tissues. A total of 163 conserved miRNAs and 303 level miRNAs were identified and some of them showed specifically up or down-regulated expression in different tissues. Some unigenes involved in Plant-pathogen interaction, MAPK signaling pathway-plant and signal transduction of auxin, zeatin, abscisic acid and jasmonic acid were targeted by miRNAs from P. heterophylla. A large number of mRNAs of four pathogens were targeted by miRNAs from P. heterophylla and several genes were targeted by miR414. Predicted target genes of a conserved and novel miRNA were validated by 5 0 RLM-RACE. These results may help to explain the function of miRNAs in the morphogenesis of organs and defense of pathogens, and provide theoretical basis for breeding and genetic improvement of P. heterophylla.

Materials and methods
Plant material and total RNA isolation P. heterophylla plants (Shitai NO.1) were selected by our group and certified as a national variety by Guizhou Crop Variety Approval Committee in 2016 (Accession number: 2016002). It has become a very popular variety in Guizhou Province because of its excellent quality and high yield. Plants were grown in the experimental field of Shibing County, Guizhou Province, China. The leaves, stem, and tuberous root (cortex and xylem) of P. heterophylla were collected separately from three randomly selected individuals. After cleaning, all samples were immediately frozen in liquid nitrogen. Transzol Plant RNA Extraction Kit (TaKaRa, Tokyo, Japan) was used for the extraction of total RNA according to the manufacturer's instructions. Traces of remaining DNA were removed by using DNase I (Takara, Tokyo, Japan).

cDNA library preparation and small RNA library construction
For a better understanding of miRNAs and mRNAs networks in P. heterophylla, the cDNA library of these tissues was also constructed. The total RNA of four tissues was purified by using Oligo(dT) magnetic beads, respectively. The mRNA was disrupted into small fragments and used for second-strand cDNA synthesis. These cDNA fragments were ligated with the Illumina paired-end sequencing adaptors and sequenced by using the Illumina Hiseq4000 platform. The expression level of each gene in different tissues was normalized by using the reads per kb per million reads (RPKM) and the RPKM of each unigene refers to the sum of RPKM from all isoforms of the same gene.
Small RNA libraries for each sample were prepared independently using Truseq TM Small RNA sample prep Kit (Illumina, U.S.A). Total RNAs of P. heterophylla were firstly size-fractioned on 15% denaturing polyacrylamide gel and small RNAs (18-35nt) fraction was collected. Next, T4 RNA ligase was used to ligate the small RNAs and adapters (5 0 and3 0 ). Small RNAs with adapters were reverse-transcribed using SuperScript II Reverse Transcriptase (Invitrogen, U.S.A) and PCR amplified with PCR. The amplification products were excised from the Novex 6% TBE-Urea Gel. The cDNA samples were sequenced using Illumina Hi-seq2000 at Shanghai Majorbio Bio-pharm Technology Co., Ltd. (Shanghai, China). To obtain the clean reads, reads with low-quality reads, adaptor sequences, and the reads shorter than 18nt and longer than 30nt were removed.

Identification of miRNAs and differential expression analysis of miRNAs
These high-quality sRNA sequences were mapped to the mRNA transcriptome database. These perfectly mapped sRNAs were used for further analysis, while rRNA, tRNA, small nuclear (snRNA), and small nucleolar RNA (snoRNA) were removed. The remaining sRNAs were aligned with the mature miRNA or miRNA precursor from all plant species in the miR-Base Sequence Database (http://www.mirbase.org/, release 22.1, October 2018) with the allowance of two mismatches. These perfectly matched sRNAs were collected and considered to be conserved miRNAs. Potentially novel miRNAs were predicted using MIREAP (http:// sourceforge.net/projects/mireap/), and the secondary structure of the sRNA precursor was calculated by using the RNAfold web server [32] (http://rna.tbi.univie.ac.at/cgi-bin/RNAfold.cgi) with default parameters. Those sRNAs fitting all the published criteria and with no homologous sequence in the public databases were considered to be a putative novel miRNA precursor. The expression level of miRNAs was normalized to one million by total clean reads of miRNAs in each sample sequence. For the novel miRNAs, the names from Phn-miR1 to Phn -miR303 were given.
TPM (Transcripts per million) was used for identifying differentially expressed miRNAs in different tissues. When log FC (log fold expression change) was greater than 1 or less than -1 (FDR<0.001, P-value<0.005), miRNAs were considered to be differentially expressed.

Prediction of miRNA targets
The prediction of miRNA targets was performed by using the Target Finder (https://github. com/carringtonlab/TargetFinder) against the transcriptome sequence data of P. heterophylla in previous study. Maximum expectation and complementarity scoring were used for evaluating the miRNA-target mRNA interactions according to the method described by Allen et al [19].

GO and KEGG analysis
For a better understanding of the function of the target genes and their corresponding metabolic network regulated by these miRNAs, Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) were used. All dates were corrected by Bonferroni. When the corrected p-value (FDR) is less than 0.05, it is considered that the GO or KEGG function is significantly enriched.

Co-expression analysis of miRNAs and their target genes
To analyze the regulation of miRNAs and their target genes, the expression level of Spearman correlation coefficients between the miRNAs and their target genes were performed.

Target validation of miRNAs by 5 0 RLM-RACE
The cleavage site of predicted targets was determined by using a 5 0 RLM-RACE (Takara, Tokyo, Japan) according to the instructions of the manufacturer. Briefly, the RNA adapter was ligated to the mRNA using the T4 RNA ligase and reverse transcribed to cDNA. Gene-specific reverse primers and 5eadapter primers were used to amplify the cleaved transcripts. The primers used to amplify cleavage products of miRNA target genes in P. heterophylla are listed in the S16 Table. The amplifying fragments were cloned to pMD19-T vectors (Takara, Tokyo, Japan) and sequenced by M13 primers.
Supporting information S1  Table. Up-regulated and down-regulated miRNAs and their target genes in leaf. (XLSX) S11 Table. Up-regulated miRNAs and their target genes in stem. (XLSX) S12 Table. Up-regulated and down-regulated miRNAs and their target genes in xylem of the tuberous root. (XLSX) S13 Table. Up-regulated miRNAs and their target genes in bark of the tuberous root. (XLSX) S14 Table.