Identification and Analysis of microRNAs in the SAM and Leaves of Populus tomentosa

Abstract

The shoot apical meristem (SAM) is a crucial tissue located at the tops of plants which can continually grow and differentiate to develop into all aboveground parts. SAM development is controlled by a series of complicated molecular regulation networks, among which microRNAs (miRNAs) and their target genes play key roles. However, little is known about these miRNAs in woody plants. In this study, we used small RNA (sRNA) sequencing to build four libraries derived from shoot tips and mature leaf tissues of Populus tomentosa, and identified 99 known miRNA families. In addition, 193 known miRNAs, including phytohormone-, developmental-, and cellular process-related miRNAs, showed significant differential expression. Interestingly, quantitative real-time reverse transcription polymerase chain reaction (PCR) analysis of miR172, miR164, and miR393 expression showed marked changes in expression patterns during the development of shoot tips. The target genes of these miRNAs were involved in the regulation of hormone responses and stem cell function. In particular, the miR172 target APETALA2 (AP2), involved in the maintenance of stem cells in the shoot apex, was expressed specifically during the initial active stage of development. These findings provide new insights into the regulatory mechanisms of miRNAs involved in SAM development and differentiation in tree species.

1. Introduction

Meristems are stem cell niches in plants that are capable of unlimited self-renewal and generate differentiated descendants [1]. The shoot apical meristem (SAM), located in the shoot apex, gives rise to all aboveground organs, such as the stems, leaves, and flowers, forming the aerial structure of higher plants. In woody plants, the SAM drives aerial growth and organogenesis throughout the life cycle.

Previous studies on the regulatory mechanisms of the shoot apex growth cycle have mainly focused on hormonal regulation [2,3]. Auxin is a vital regulator in the development of the SAM and the formation of secondary structures in woody plants due to its endogenous spatial distribution and polar auxin transport. Inhibition of polar auxin transport prevents leaf formation at the SAM and leads to reduced cell division in the stem cell niches [4]. Cytokinins (CTKs) are involved in meristem maintenance and positively regulate cell division. Reduced CTK content results in a dramatic reduction in the size of the SAM and the leaf-initiation rate [5]. In addition, a newly discovered plant hormone, strigolactone (SL), regulates growth redistribution of shoots by modulating auxin transport in the apical meristem [6]. A series of studies indicated that phytohormones coordinate with each other during the differentiation process of the SAM. In addition to phytohormones, transcription factors also regulate the development of the shoot apex [7,8]. Class III homeodomain-leucine zipper (HD-ZIP III) transcription factors are associated with cell division and CTK response, acting as key transcription factors involved in the regulation of the maintenance of stem cells and governing the initiation of lateral organs from the SAM [9]. Furthermore, maintenance and differentiation of the SAM demands synergistic interactions between hormones and transcription factors. The relationship between class I knotted1-like homeobox (KNOX1) transcription factors and the ratio of CTK to gibberellin (GA) is a vital factor in the switch to organogenesis in the apical meristem [10]. Transcription factors of the auxin response factor (ARF) family alter the spatial distribution of auxin by binding specifically to auxin response elements of genes, and then cooperate with GAs to accelerate cell division and differentiation in the apical meristem [11]. These previous studies were mainly concerned with the physiological and transcriptional changes that occur during the development of the SAM. Few studies have focused on the post-transcriptional regulation of the maintenance and differentiation of stem cell niches in the SAM.

MicroRNAs (miRNAs) regulate gene expression at the post-transcriptional level and play a dominant role in plant stem cell self-renewal and differentiation, normally by suppressing their target mRNAs [12,13]. For instance, miR166 exerts an influence on shoot apical stem cell differentiation by inhibiting expression of the target genes HD-ZIP III and WUS [14,15,16]. In addition, miR156 targets SPL transcription factor genes that have a specific function in the leaf development and apical dominance of higher plants [17]. APETALA2 (AP2) mainly participates in controlling differentiation during meristem development, and is the target gene of miR172 [18]. Moreover, there is an antagonistic relationship between miR156 and miR172 during the development of the apical meristem [19]. The regulatory roles of miRNAs in stem cells have been reported in herbaceous plants such as Arabidopsis (Arabidopsis thaliana), but little is known about their roles in woody plants.

White poplar (Populus tomentosa) is a rapidly growing woody plant that is widely distributed in China as an imporant economic tree species. Therefore, it is important to obtain knowledge of miRNA regulation in the SAM to improve the study of the growth and development in P. tomentosa. Here, we present a deep-sequencing profile on a genome-wide scale to identify miRNAs in the SAM of P. tomentosa. Millions of small RNA (sRNA) reads were obtained, and many known and novel miRNAs in SAM and leaf libraries were established. Several miRNAs that were differentially expressed between the SAM and its surrounding leaves and involved with phytohormones and development were identified and analysed. Furthermore, the dynamic expression patterns of several important function-related miRNAs during shoot growth were also determined.

2. Materials and Methods

2.1. Plant Materials

Healthy shoot tips and the surrounding leaves of P. tomentosa were collected from 20-year-old plants growing under natural conditions at Yangzhou University, Yangzhou, China (32.39° N, 119.42° E). To cover the major stages of the reactivating–activation cycle, three developmental stages were chosen as follows: (1) Reactivating stage (Re): The samples were collected in March when no leaves were attached to the shoot tips; (2) initial active stage (IA): The samples were collected in April, with several young light-green leaves near the shoot tips; (3) stable active stage (SA): The samples were collected in June, with some mature dark-green leaves on the stem. These stages are shown in Figure 1. These samples were collected from three independent trees at each time point, and were flash-frozen in liquid nitrogen and stored at −80 °C for further experiments.

2.2. Morphological Observations

Samples corresponding to the three developmental stages were fixed in 2.5% (v/v) glutaraldehyde in 100 mM sodium phosphate buffer (pH 7.2) for 2 h at 4 °C. The samples were then rinsed in 100 mM phosphate buffer three times at 15 min intervals. After successive graded ethanol dehydration (30%, 50%, 70%, 80%, 90%, 95%, and 100%; 15 min each), samples were dried using the Leica automatic critical point drying instrument (Leica EM CPD300, Leica Microsystems GmbH, Wetzlar, Germany). Finally, the shoot tips and leaves were embedded in Spurr’s resin. A Leica microtome was employed to section samples into 2 μm sections, which were then stained with 0.25% (w/v) toluidine blue O (Sigma). The stained semi-thin sections were observed under a Zeiss Axioskop 2 Plus microscope equipped with a computer-assisted digital camera (Carl Zeiss, Jena, Germany).

2.3. sRNA Library Construction and High-Throughput Sequencing

Samples from the IA, which is a rapid growth stage, were used for high-throughput sequencing, as shown in Figure 1. To construct sRNA libraries, total RNAs were extracted from the SAMs and surrounding leaves using the TaKaRa MiniBEST plant RNA extraction kit (Takara, Dalian, Liaoning, China) according to the manufacturer’s instructions. Subsequently, total RNAs from samples were separated into RNAs of different sizes using polyacrylamide gel electrophoresis. The 18–30 nt long sRNAs were excised and recovered, and then ligated to 5′ (AATGATACGGCGACCACCGAGA) and 3′ (CAAGCAGAAGACGGCATACG) RNA adapters using T4 RNA ligase 2. The ligation products were further purified using a Novex 10% TBE-urea gel, reverse transcribed, and then amplified by PCR using SuperScript II reverse transcriptase (Invitrogen, Carlsbad, CA, USA) and Phusion High-Fidelity DNA polymerase (New England BioLabs, Ipswich, MA, USA). The amplification products were purified on a Novex 6% TBE-urea gel. The purified PCR products were used to construct sRNA libraries and were sequenced on the Illumina HiSeq 2500 platform (Illumina, San Diego, CA, USA) at Genepioneer Biotechnologies (Nanjing, China). All raw sequence data have been deposited in the NCBI database under accession number GSE100282.

2.4. Bioinformatic Analysis of Sequencing Data

Raw sequence data were obtained via the RNA high-throughput sequencing process, and then filtered to remove contaminating reads, adapter sequences, poly-A tails, and other artefacts. All high-quality sequences, including those with only a single unique read, were used for further analyses. The GenBank [20] and Rfam [21] databases were employed to annotate the ribosome RNAs (rRNAs), small nuclear RNAs (snRNAs), small nucleolar RNAs (snoRNAs), and transfer RNAs (tRNAs) within the sRNA library [22]. Repeat overlapping sequences were annotated as repeat-associated sRNAs. Matching sRNAs without annotations were stored in an unannotated file. The unannotated reads were mapped to the Populus trichocarpa genome database [23] downloaded by SOAP (short oligonucleotide alignment program) according to the default settings [24]. The mapped reads were used to predict potential novel miRNAs.

To identify known miRNAs from the SAM and leaves of P. tomentosa, the bioinformatics pipeline was used to annotate known miRNAs. In the present study, all miRNAs belonging to families already annotated in the miRBase Registry [25] in the P. trichocarpa database were defined as known [26]. Potentially novel miRNAs and their hairpin precursors were predicted using the MIREAP [27] software with default parameters from the remaining unannotated reads [28]. The novel miRNAs identified were given names with the precursor “novel-m#” (e.g., novel-m0001).

The miRNA frequencies were normalised to one million, based on the total number of clean reads in each library (transcripts per million [TPM] = miRNA count/total count of clean reads × 1,000,000). Differential expression between the SAM and leaf libraries was calculated as the fold change = log₂ (SAM/Leaf). Both |fold change| > 1 and false discovery rate (FDR) < 0.05 were used to define the differentially expressed miRNAs.

2.5. Prediction of miRNA Targets

Target genes of known and novel miRNAs were predicted using the web-based application psRNATarget [29] and psRobot [30] with the default parameters. Sequences from P. trichocarpa were used as a custom target database. The networks between miRNAs and target genes were constructed using Cytoscape vr3.2.1 [31].

2.6. qRT-PCR Analysis of Differential Gene Expression

All samples from the three developmental stages were used for quantitative real-time reverse-transcription PCR (qRT-PCR), as illustrated in Figure 1. To assess changes in expression levels of miRNAs from the SAM and leaf tissues during the three developmental stages, the expression levels of several differentially expressed function-related miRNAs were determined by poly (A) tailing-based RT-PCR using the miRcute miRNA First-strand cDNA (complementary DNA) Synthesis Kit (Tiangen, Beijing, China) with three replications. The forward primers for the miRNAs were designed using Primer Premier 5.0 (Table S1), and the reverse primers were provided in the miRcute Plus miRNA qPCR Kit (SYBR Green) (Tiangen, Beijing, China). U6 was used as the reference gene of miRNAs [32]. The corresponding target genes were determined by qRT-PCR using the SYBR Premix Ex Taq kit (TaKaRa, Dalian, China) on the Bio-Rad CFX96 real-time PCR detection system (Bio-Rad, Hercules, CA, USA) with three replications. The total RNA was reverse transcribed into cDNA using the PrimeScript RT Reagent Kit (TaKaRa, Dalian, China), and then diluted 10 times with sterile water. The primers of the targets are shown in Table S1, and the reference gene, Actin 2 (ACT2), was selected according to a previous study [33]. qRT-PCR reactions contained 8.5 μL of double-distilled H₂O (ddH₂O), 2 μL mixtures of forward and reverse primer, 2 μL of diluted cDNA, and 12.5 μL of SYBR dye in final volumes of 25 μL. The qRT-PCR cycling program was as follows: Polymerase activation at 94 °C for 30 s, followed by 40 cycles of 94 °C for 5 s, 60 °C for 30 s, and 72 °C for 10 s. The specificity of primer amplicons was tested using the 2^−ΔΔCT relative quantification method to calculate relative changes in the expression levels of miRNAs and corresponding target genes during development [34]. At each stage, the significance of differences between SAM and leaf expression levels was analyzed using Student’s t-test at the probability level of 0.05 (p < 0.05). To test the significant differences of gene expression levels in the same tissue (SAM or leaves) among the three growth stages, we used One-way ANOVA with Duncan’s Multiple Range Tests at p < 0.05. All these data analysis were performed using SPSS 18.0 software for Windows (Table S2) [35].

3. Results

3.1. Morphological and Anatomical Structure Changes of Shoot Tips and Leaves During Growth

The shoot tips and surrounding leaves showed seasonal variation in their morphology and anatomical structures during the Re, IA, and SA stages, as shown in Figure 1. We initially performed morphological observations of the shoot tips and leaves to determine the timing of the reactivating–activation transition. As shown in Figure 2a–c, the shoot tips in the Re were larger than those in other phases, there were no leaves on the stem until IA, and the leaves in the SA were more mature than those in the IA. Furthermore, we observed the micro-structures of the shoot apex and leaves during the various developmental stages. The results showed that the cells of shoot tips were difficult to distinguish from other parts during the Re, as shown in Figure 2d. In the IA, the shoot tips had a distinct tunica-corpus structure, which is one of the characteristics of the SAM [36]. In this developmental stage, apical histological zonation was visible, with several clonally distinct layers of cells, as shown in Figure 2e (arrows). The SAM had no evident histological zonation during the SA, and leaves were derived from the epidermal layer, as shown in Figure 2f. In addition, leaf cells were arranged more compactly in spongy tissues in the SA than in the IA, as shown in Figure 2g,h.

3.2. Sequencing and Analysis of sRNAs in the SAM and Leaves

To identify miRNAs involved in regulatory networks in the SAM and leaves of P. tomentosa, we constructed four sRNA libraries from the SAM and surrounding mature leaf tissue, with two biological replicates, and sequenced them on the Illumina HiSeq 2500 platform. Approximately 12 million raw reads were obtained for each sample. After cleaning up the data, we obtained 6,782,666 (SAM 01), 8,067,503 (SAM 02), 9,195,827 (leaf 01), and 8,775,153 (leaf 02) high-quality clean reads between 18 and 30 nt in length for use in further analyses (Table 1). Among the clean reads, 1,738,201 (leaf 01), 1,653,933 (leaf 02), 1,960,095 (SAM 01), and 2,390,698 (SAM 02) were unique. After removal of rRNA, tRNA, snRNA, and snoRNA sequences and repeat reads, 893,192 (leaf 01), 795,507 (leaf 02), 1,124,667 (SAM 01), and 1,462,646 (SAM 02) sRNA unannotated reads remained (Table 1). Among the unannotated reads, 614,301 (leaf 01), 474,568 (leaf 02), 314,304 (SAM 01), and 462,045 (SAM 02) reads were mapped to the reference genome (Table 1).

The overall size distributions of the filtered sRNAs from the SAM and leaf libraries are summarized in Figure 3a. The sRNA sizes from the leaf and SAM libraries both exhibited a major peak at 21 nt. The proportions of 21 nt sRNAs were approximately 18.08% (leaf 01), 14.71% (leaf 02), 23.48% (SAM 01), and 20.77% (SAM 02). The leaf libraries showed minor peaks at 22 nt (leaf 01: 16.34%; leaf 02: 14.08%), whereas the SAM libraries showed minor peaks at 24 nt (SAM 01: 16.20%; SAM 02: 19.23%).

3.3. Identification of Known and Novel miRNAs in P. tomentosa

After a sequence similarity search, we identified 349 known miRNAs corresponding to 99 miRNA families. In the SAM, the 328 known miRNAs were grouped into 88 miRNA families. Most of the known miRNA families have several members, and the top five families (miR169, miR171, miR166, miR478, and miR395) possessed 30, 17, 17, 16, and 11 members, respectively, as shown in Figure 3b. In the leaves, the 314 known miRNA sequences were grouped into 89 miRNA families. The numbers of miR478, miR166, miR171, miR169, and miR395 family members in the leaves were 17, 17, 17, 14, and 11, respectively. We also analysed the counts of reads in the known miRNA families, and they indicated large divergences in expression frequency among these miRNA families. The number of reads for conserved miRNAs sequenced from the same or different miRNA families also varied drastically, with miR166, miR159, and miR319 representing the top three families in the SAM and miR166, miR396 and miR159 in the leaf libraries, as shown in Figure 3c. miR166 family members were the most abundant, with 1,422,348.5 and 1,092,946 total reads in the SAM and leaf libraries, respectively. In addition, several miRNAs such as miR1449, miR397, and miR6454 were absent from one of the two libraries (Table S3). This indicated that the miRNA population present in the SAM library differed from that in the leaf library.

The remaining unannotated sRNAs were queried against the reference genome, and the mapped reads were used to identify potential novel miRNAs. A total of 1609 mature miRNA sequences were recognised as potential novel miRNAs (Table S4). The high number of potential novel miRNAs might be contributed by their isomirs (for many miRNAs, the same miRNA has several length and/or sequence variants) [38]. The lengths of these mature miRNAs ranged from 19 to 26 nt, with the most common being 24 nt. The minimal free energy range of these novel miRNA precursors was −66.5 to −18 kcal/mol, with an average of −26.23 kcal/mol.

3.4. Analysis of Differentially Expressed miRNAs (DEMs) and their Targets

We quantitated the miRNAs in each sample and identified those that were significantly differentially expressed between the SAM and leaves (|log₂ fold-change| > 1, FDR < 0.05). The results showed that 105 known miRNAs were differentially expressed between the SAM and leaves (Figure 4a; Table S5). Among them, 47 known miRNAs, corresponding to 12 miRNA families, were significantly up-regulated in the SAM, while 58 known miRNAs, belonging to 21 miRNA families, were down-regulated, as shown in Figure 4b,c. In addition, 88 novel miRNAs were DEMs, of which 26 were up-regulated and 62 were down-regulated in the SAM, as shown in Figure 4b.

To further investigate the functions of the miRNAs, we searched for potential miRNA target genes. In total, 5904 and 5598 genes were predicted to be targets of the known and novel DEMs, respectively (Table S6). Moreover, we performed Gene Ontology (GO) functional analyses to evaluate the putative functions of these potential target genes. Then, GO functional enrichment analysis of the target genes was carried out with the GOseq R package [39]. The results revealed enrichment of important GO subcategories such as “cellular process”, “signalling”, and “cellular component organisation or biogenesis” within the biological process category. Similarly, the “transporter activity” and “molecular transducer activity” subcategories within the molecular function category were enriched (Figure S1a). We further analysed the biological functions of the miRNA-targeted genes using the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database [40]. As shown by the blue arrows in Figure S1b, the most highly represented pathways included “biosynthesis of unsaturated fatty acid”, “zeatin biosynthesis”, and “RNA degradation”.

3.5. DEMs Related to Stem Cell Functions

miRNA regulation plays a role in hormonal responses. We found that 23 predicted targets of 53 miRNAs are ARFs involved in the auxin-activated signalling pathway. Among the 53 miRNAs, 30 miRNAs belonging to the ptc-miR319, ptc-miR159, ptc-miR160, ptc-miR171, ptc-miR390, and ptc-miR6425 families were up-regulated in the SAM, whereas the remaining 23 miRNAs, including ptc-miR396f, ptc-miR473a-3p, ptc-miR397a, and novel-m0531-5p, were down-regulated (Figure S2a). In addition, three putative targets of 12 miRNAs, including five members of the ptc-miR396 family and seven novel miRNAs, were auxin transporter-like proteins. Nine of these miRNAs were down-regulated in the SAM (Figure S2b). Five genes belonging to the transport inhibitor response 1 (TIR1) family were predicted to be the targets of miR393 and miR6425. miR393 displayed lower expression in the SAM, whereas miR6425 showed higher expression (Figure S2c).

Furthermore, we predicted 10 target genes of the ptc-miR159 family (ptc-miR159a, ptc-miR159b), ptc-miR6459a-3p, ptc-miR1444a, ptc-miR473a-5p, ptc-miR171l-5p, and six novel miRNAs involved in zeatin biosynthesis (Figure S3). Of these miRNAs, ptc-miR6459a-3p was up-regulated in the SAM, whereas ptc-miR159a, ptc-miR159b, ptc-miR1444a, ptc-miR473a-5p, and ptc-miR171l-5p were down-regulated.

Nine predicted target genes were AP2-like family proteins. Among them, POPTR_0005s18400, POPTR_0007s10780, POPTR_0008s04490, and POPTR_0016s08530 were targeted by the miR172 family (ptc-miR172a, ptc-miR172b-3p, ptc-miR172c, and ptc-miR172f), and POPTR_0001s03300 and POPTR_0001s19360 were targeted by ptc-miR396g-5p, novel-m1792-3p, and novel-m1796-5p. Moreover, POPTR_0007s14690 and POPTR_0014s01260 were targets of ptc-miR6457b, and POPTR_0006s18330 was the target of novel-m0773-5p and novel-m1179-5p. All these miRNAs were down-regulated in the SAM, as shown in Figure 5a,b. In addition, two target genes were annotated as CUP-SHAPED COTYLEDON 2 (CUC2) and were predicted to be targets of the ptc-miR164 family, ptc-miR1444a, and novel-2205-3p. The expression levels of these miRNAs were also lower in the SAM than in leaves, as shown in Figure 5c,d. POPTR_0003s04860 and POPTR_0001s18930 were annotated as homeobox-leucine zipper protein ATHB-15 and were targets of the ptc-miR172 family, ptc-miR6421-3p, novel-m0791-3p, novel-m2501-5p, and novel-m2517-5p, as shown in Figure 5e,f.

3.6. Expression of miR172 and its Targets

We focused on members of the miR172 family because of their expression patterns in the network, which are shown in Figure 5b,f. miR172 is highly expressed in leaves but significantly down-regulated in the SAM. This miRNA has been shown to target AP2s, including POPTR_0005s18400, POPTR_0007s10780, POPTR_0008s04490, and POPTR_0016s08530. A closer inspection of the sequences revealed two highly conserved AP2/EREBP domains (Figure 6a, red boxes) and a region of approximately 20 nt that appeared to be targeted by miR172 (Figure 6a). Furthermore, we performed qRT-PCR of miR172 and its targets at three different growth stages (the Re, IA, and SA) in the SAM and leaves. The results showed that the miR172 and targets were negatively correlated at specific stages and tissues. For instance, in the SAM, the expression level of ptc-miR172 decreased from the Re to IA (Figure 6b), whereas its targets (POPTR_0007s10780 and POPTR_0005s18400) increased (Figure 6c,d). During the IA, the expression of ptc-miR172 was lower in the SAM than in leaves, but POPTR_0007s10780 and POPTR_0005s18400 showed the opposite pattern (Figure 6b–d). In leaves, miR172 expression increased from the Re to IA, whereas the expression levels of POPTR_0007s10780 and POPTR_0005s18400 decreased (Figure 6b–d).

3.7. qRT-PCR Validation of Other Important miRNAs

To validate the miRNA expression patterns derived from high-throughput sequencing, we performed qRT-PCR to analyse the expression patterns of 10 known miRNAs and three novel miRNAs, using the same samples used for RNA-seq. The results showed that the expression patterns of 10 of the 13 miRNAs agreed with the high-throughput sequencing data (Figure 7a). In addition, we detected the expression of several miRNAs and their targets at three different growth stages. The results showed that ptc-miR164a had a very high expression level in the SAM at the Re, but low expression levels during the IA and SA (Figure 7b). From the Re to IA, ptc-miR319a expression increased in the SAM but decreased in the leaves, and at the IA, its expression was significantly higher in the SAM than leaves (Figure 7c). In the SAM, ptc-miR393c expression significantly decreased from the Re to IA, and was lower in the SAM than in leaves at the IA (Figure 7d). A potential target gene of pct-miR164a, POPTR_0011s11600, displayed high expression levels in the SAM at the IA (Figure 7e). At all the three stages, its expression was negatively correlated with pct-miR164a (Figure 7b) in the SAM and leaves (Figure 7e). POPTR_0003s04860 was a putative target of miR172 and was expressed more highly in the SAM than in leaves at the IA, whereas miR172 showed the opposite pattern (Figure 6b and Figure 7f). As the predicted target of miR393c, POPTR_0001s33030 expression increased from the Re to IA in the SAM. At the IA, its expression was higher in the SAM than in leaves (Figure 7g).

4. Discussion

4.1. SAM and Leaf Formation

Leaves in most eudicot species are composed of derivatives from the epidermal layer, the subepidermal layer, and the corpus of the SAM [14,41]. One of the earliest markers of leaf initiation from the shoot apex is the occurrence of observable apical histological zonation in the SAM [11]. This developmental phase coincides with several changes at the biochemical and molecular levels, for instance, in carbohydrates, plant hormones, and miRNAs [42,43,44]. In the present study, morphological and anatomical observations of the SAM and leaves in P. tomentosa provided details of structural changes during their development transitions. Our results showed that the shoot tips had not differentiated into lateral organ primordia in the Re. In the SA, some leaves had derived from the apex region of the SAM. However, in the IA, several clonal layers of cells were visible during the reactivating–activation transition. The SAM could be distinguished from leaf primordia in the samples of the IA. Therefore, the samples corresponding to this developmental phase (IA) were ideal for analysis of sRNAs that were expressed differentially between the SAM and surrounding leaves in P. tomentosa.

4.2. Complexities of sRNA Generation in P. tomentosa

The differentiation processes of shoot tips play a key role in tree growth and development, including maintaining the integrity of the stem cell, cell division, cell enlargement, and the reiterative formation of lateral organs [45]. Increasing evidence has indicated that miRNAs participate in the developmental processes of the SAM [46,47]. Nevertheless, little is known about the non-coding RNA mechanisms involved in the developmental phase between the SAM and its derivatives in woody plants. In the present study, we generated sRNA libraries of the SAM and leaves, with approximately 6–8 million clean reads per sample. From these sRNA sequences, more than 1.9 and 1.6 million unique sequences were obtained for the SAM and leaf libraries, respectively. Notably, the percentage of unannotated unique sRNAs in the SAM library was very high, indicating that many potential miRNAs have not yet been investigated in the SAM of P. tomentosa.

In plants, the lengths of sRNAs generally range from 20 to 40 nt [48]. The length distribution patterns of sRNAs in A. thaliana, Oryza sativa, and Medicago truncatula display a major peak at 24 nt, followed by at 21 nt [19,49,50]. However, a previous study showed that 21 nt sRNAs were the most abundant, followed by 24 nt sRNAs, in some plants such as tomato, grape, and several Populus spp., including P. tomentosa, Populus balsamifera, and Populus euphratica [51,52,53,54,55]. In our study, 21 nt sRNAs were the most abundant class, while 24 nt sRNAs in the SAM and 22 nt sRNAs in the leaves were the second most abundant respectively, suggesting that the length distributions of sRNAs in different tissues in Populus are diverse.

4.3. miRNA Population Involved in SAM Development and Maintenance

Evidence has accumulated indicating that miRNAs might act as a node molecule involved in the coordinated regulation of plant developmental or cell differentiation processes [56,57,58]. Many miRNAs have been investigated in Populus since high-throughput sequencing technology applications have become widespread [56,59,60]. However, little is known about the involvement of miRNAs in the endogenous molecular changes in the SAM. In the present study, several previously identified miRNAs, such as ptc-miR166, ptc-miR159, and ptc-miR319, were over-represented in both libraries. These highly expressed miRNAs were confirmed to be associated with cell division and meristem activity in many previous reports [61,62,63]. In addition, we determined that more than 5000 genes were predicted to be targets of known miRNAs. Among them, numerous potential targets were involved in processes such as hormone responses and cellular processes. This suggests that these high-expression miRNAs may modulate the expression of genes during meristem development and cell division in poplar.

Hormones play key roles in plants, and the hormones auxin and CTK are crucial for positioning and maintaining meristem activity. In the SAM, local auxin accumulation drives organogenesis, and CTKs are essential for cell proliferation [64,65]. In the present study, target genes of miRNAs, including ptc-miR319, ptc-miR160, ptc-miR390, and ptc-miR396, were predicted to be ARFs and auxin transporter-like proteins. Ptc-miR159, ptc-miR1444, and ptc-miR473 target corresponding genes that participate in zeatin biosynthesis. These miRNAs showed significant differential expression between the SAM and leaves, suggesting their different functions in SAM and leaf development. In addition, auxin was reported to exert its functions via the F-box transcription factor family, including four partially redundant genes. In particular, TIR1 is the most extensively researched gene, and is targeted by miR393 [66]. A mutant unable to produce miR393 would lead to a noticeably increased expression level of TIR1 mRNA as well as concomitantly enhanced auxin sensitivity, and would cause pleiotropic effects on plant development including overproduction of lateral organs, developmental abnormalities of leaves, and delayed flowering [67,68]. In the present study, ptc-miR393 targeted several TIR1 genes in our database. According to the qRT-PCR results, the TIR1 gene POPTR_0001s33030 exhibited an opposite expression pattern to that of ptc-miR393 in the SAM and leaves at the Re, IA, and SA. Thus, ptc-miR393 is very likely to regulate meristem development and auxin response by suppressing the expression of the auxin receptor TIR1.

Several studies have established the involvement of miRNAs in the regulation of appropriate stem cell activity. AP2, a second-class A gene, is targeted by miR172. AP2 has been proven to promote stem cell retention, whereas miR172 promotes cell differentiation. In the floral meristem of Arabidopsis, flowers overexpressing miR172-resistant AP2 exhibited numerous stamens flanking an indeterminate floral meristem [18]. In our study, we found that miR172 was down-regulated in the SAM, and its four predicted target genes were AP2 genes. Furthermore, qRT-PCR indicated inverse expression patterns between ptc-miR172 and its two AP2 targets in samples from both the Re and IA. Therefore, ptc-miR172 likely regulates the reactivating–initial activation phase change through translational repression of AP2 genes related to meristem cell retention.

The miR164 family was found to target CUC. Members of the CUC gene family have been verified to control meristem establishment and maintenance [69,70,71]. Overexpression of miR164 resulted in defects in axillary meristem formation, and miR164 mutants developed accessory buds in leaf axils [71,72]. The HD-ZIP transcription factor regulates cell division and meristem activity. Members of the HD-ZIP transcription factor family include PHABULOSA (PHB), PHAVOLUTA (PHV), REVOLUTA (REV), AtHB8, and AtHB15 in Arabidopsis [14]. In our study, ptc-miR164 family members targeted two CUC2 genes, and miR172 targeted ATHB-15. Both miR164 and miR172 were down-regulated in the SAM, and this was validated by qRT-PCR. Moreover, we detected their corresponding targets, which showed opposite expression patterns to those of miR164 and miR172, implying that miR164 and miR172 may play important roles in meristem cell differentiation through repression of their target genes.

5. Conclusions

In this study, sRNA-seq was performed to explore important miRNAs related to SAM development and differentiation in P. tomentosa, resulting in the identification of 349 known miRNAs and 1609 potential novel miRNAs in the SAM and leaves. Among them, 105 known miRNAs and 88 novel miRNAs were differentially expressed, and 5904 and 5598 genes were predicted to be target genes of these DEMs, respectively. Based on the functions of the target DEMs, we constructed networks among the miRNAs and targets associated with hormone responses. Several target genes of ptc-miR393, ptc-miR319, ptc-miR159, and others were ARF and TIR1 genes, and the target genes of ptc-miR159, ptc-miR473, and ptc-miR1444 were involved in zeatin biosynthesis. In addition, after combining the results of the sRNA-seq and qRT-PCR analyses, we identified miR172 and miR164 as two vital regulators of stem cell retention in the SAM of P. tomentosa. These results suggest that a miRNA–mRNA regulatory network is involved in the post-transcriptional regulation of genes associated with SAM development and organogenesis in poplars.

Supplementary Materials

The following are available online at http://www.mdpi.com/1999-4907/10/2/130/s1, Figure S1: Functions of predicted target genes, Figure S2: Heatmaps and networks between differentially expressed miRNAs and their targets involved in IAA signaling pathway, Figure S3: Heatmaps and networks between differentially expressed miRNAs and their targets involved in Zeatin biosynthesis, Table S1: Primers used for qRT-PCR in this study, Table S2: Significant difference analysis of expression levels among three growth stages, Table S3: Known miRNAs belonging to different families expressed in the SAM and leaf of P. tomentosa, Table S4: Potential novel miRNAs identified in P. tomentosa, Table S5: Differentially expressed miRNAs between SAM and leaves, Table S6: miRNAs and predicted targets with annotation.

Author Contributions

B.J. and N.T. conceived and designed the experiments. J.C. and W.L. wrote the paper. J.C., W.L., Z.L., B.Z., L.W., and S.R. performed the experimental work and contributed the bioinformatics analyses. All authors read the script, provided comments, and approved the final manuscript.

Funding

This research was funded by the National Key Research and Development Program (2016YFD0600105) and the Natural Science Foundation of China (grand No. 31670695).

Acknowledgments

We thank Genepioneer Biotechnologies for helping with sRNA-sequencing and technical assistance. We are grateful to Likui Zhu, Kaibiao Ma and Zhichong He for the help of sampling and preparing of shoot tips and leaves in P. tomentosa.

Conflicts of Interest

The authors declare no conflict of interest.

Figures and Table

Figure 1. Samples used for high-throughput sequencing and qRT-PCR analyses of P. tomentosa. (a) Reactivating stage, (b) initial active stage, and (c) stable active stage. The materials used for small RNA sequencing (sRNA-seq) are shown in the blue dashed box and those used for qRT-PCR in the green dashed box. SAM: shoot apical meristem.

Figure 2. Morphological and anatomical structures of the shoot apical meristem (SAM) and leaves of Populus tomentosa. (a–c) Morphological structures of the SAM and leaves at various developmental stages. (d) Semi-thin section of a compound sample of the shoot tip in the reactivating stage. Micro-structures of the SAM (e, red arrow and dashed box) and leaf (g) in the initial active stage. (f) Cross-sectional view of the SAM structures in the stable active stage (red dashed box: Shoot apex; yellow dashed box: Upward growing young leaves). (h) Cross-sectional view of the leaf structures in the stable active stage.

Figure 3. (a) Length distributions of the small RNA (sRNA) sequences. (b) Number of members within the top 20 known miRNA families. (c) Read counts corresponding to the top 20 known microRNA (miRNA) families.

Figure 4. Differentially expressed miRNAs. (a) Scatter diagram of the differential read counts of the known and novel miRNAs. Each dot in the figure represents one miRNA. Red and green dots represent differentially expressed miRNAs. (b) Hierarchical cluster and expression profile analyses of the differentially expressed known and novel miRNAs between the shoot apical meristem (SAM) and leaves. Each column represents a sample, and each row represents a single miRNA. The hierarchical cluster was obtained using MultiExperiment Viewer (MeV). The data in the heatmaps are the values of log10 transcripts per million (TPM). The bar indicates the expression scale for miRNAs in the SAM and leaves; grey represents non-expressed miRNAs. (c) Known miRNA families up- and down-regulated in the SAM compared with leaves.

Figure 5. Heatmaps and networks of the differentially expressed miRNAs and their targets AP2 (a,b), CUC2 (c,d), and ATHB-15 (e,f). In the heatmaps, the coloured blocks represent the log10 transcript per million (TPM) values of miRNAs in the shoot apical meristem (SAM) and leaf. Red indicates highly expressed miRNAs; blue indicates moderately expressed miRNAs; × indicates non-expressed miRNAs.

Figure 6. miR172 and its targets. (a) Conserved regions and miR172 target sites within P. tomentosa AP2 transcripts; red boxes represent two conserved AP2/EREBP domains. (b–d) qRT-PCR analysis of the expression levels of miR172 (b) and its targets POPTR_0007s10780 (c) and POPTR_0005s18400 (d). The columns and error bars indicate means and standard deviations (n = 3), respectively; * indicates significant differences (p < 0.05) between the shoot apical meristem (SAM) and leaves at each stage; Non-overlapping letters (a–c) indicate significant differences (p < 0.05) of SAM (red letters: a, b, c) or leaves (blue letters: a, b, c) among three different stages.

Figure 7. Confirmation of the expression of miRNAs derived from high-throughput sequencing. (a) The expression ratios [log2 (SAM/leaf)] of miRNAs, as determined by deep sequencing (pink bars) and qRT-PCR (blue bars). (b–g) The relative expression levels of ptc-miR164a (b), ptc-miR319a (c), ptc-miR393c (d), CUC2 (e), ATHB15 (f), and TIR1 (g) at three development stages, as determined by qRT-PCR. The columns and error bars indicate the means and standard deviations (n = 3), respectively; * indicates significant differences (p < 0.05) between the shoot apical meristem (SAM) and leaves at each stage; Non-overlapping letters (a–c) indicate significant differences (p < 0.05) of SAM (red letters: a, b, c) or leaves (blue letters: a, b, c) among three different stages.

Table 1

Analysis of sRNAs from leaf and SAM libraries.

	Leaf 01	Leaf 02	SAM 01	SAM 02
raw reads	12,654,992	12,248,285	12,309,512	11,947,552
clean reads	9,195,827	8,775,153	6,782,666	8,067,503
unique reads	1,738,201	1,653,933	1,960,095	2,390,698
rRNA	607,501	659,392	545,748	602,901
snRNA	6966	6521	7805	9105
snoRNA	13,082	13,579	14,493	17,325
tRNA	13,490	13,548	15,334	17,450
repeats	171,272	135,780	216,642	244,514
unannotated reads	893,192	795,507	1,124,667	1,462,646
mapped reads	614,301	474,568	314,304	462,045

Unique reads: Reads remaining after removal of redundant reads. Unannotated reads: Reads remaining after removal of rRNA, snRNA, snoRNA, and tRNA sequences and repeat reads. Repeats: Unique reads that could be mapped to repeat-Repbase [37]. SAM: shoot apical meristem.