Figures
Abstract
The tuberous root of Brassica rapa L. (turnip) is an important modified organ for nutrition storage. A better understanding of the molecular mechanisms involved in the process of tuberous root development is of great value in both economic and biological context. In this study, we analyzed the expression profiles of both mRNAs and miRNAs in tuberous roots at an early stage before cortex splitting (ES), cortex splitting stage (CSS), and secondary root thickening stage (RTS) in turnip based on high-throughput sequencing technology. A large number of differentially expressed genes (DEGs) and several differentially expressed miRNAs (DEMs) were identified. Based on the DEG analysis, we propose that metabolism is the dominant pathway in both tuberous root initiation and secondary thickening process. The plant hormone signal transduction pathway may play a predominant role in regulating tuberous root initiation, while the starch and sucrose metabolism may be more important for the secondary thickening process. These hypotheses were partially supported by sequential DEM analyses. Of all DEMs, miR156a, miR157a, and miR172a exhibited relatively high expression levels, and were differentially expressed in both tuberous root initiation and the secondary thickening process with the expression profiles negatively correlated with those of their target genes. Our results suggest that these miRNAs play important roles in tuberous root development in turnips.
Citation: Li J, Ding Q, Wang F, Zhang Y, Li H, Gao J (2015) Integrative Analysis of mRNA and miRNA Expression Profiles of the Tuberous Root Development at Seedling Stages in Turnips. PLoS ONE 10(9): e0137983. https://doi.org/10.1371/journal.pone.0137983
Editor: Yun Zheng, Kunming University of Science and Technology, CHINA
Received: January 2, 2015; Accepted: August 24, 2015; Published: September 14, 2015
Copyright: © 2015 Li et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited
Data Availability: All relevant data are within the paper and its Supporting Information files.
Funding: This study was supported by Modern Agriculture Industrial Technology System Funding of Shandong Province, China (SDAIT-02-022-04, http://cyjstx.sdny.gov.cn/). The role of funder: study design, decision to publish, and preparation of the manuscript. National Natural Science Foundation of China (31401869, http://www.nsfc.gov.cn/). The role of funder: study design, data collection and analysis, and preparation of the manuscript. China Postdoctoral Science Foundation Funded Project (2013M541948, http://jj.chinapostdoctor.org.cn/V1/Program1/Default.aspx). The role of funder: study design, data collection and analysis, and preparation of the manuscript. Shandong Postdoctoral Science Foundation Funded Project, China (201303030, http://www.sdhrss.gov.cn/cm/). The role of funder: study design, data collection and analysis, and preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
The modified tuberous root in many vegetables, such as radishes, carrots, turnips, and sweet potatoes provides a great food source for human consumption. Tuberous roots are rich in nutrition and dietary fiber [1], and convenient for storage. Some are also of medicinal importance [2–3]. Thus, it is important, in the context of economic value, to study the molecular mechanisms underlying the tuberous root formation and development. In additions, since tuberous roots are important modified organs, studies on tuberous root formation and development is also of high biological value.
Tuberous root development is a complex biological process that is controlled by many internal and external factors. Phytohormones, such as auxin [4,5], cytokinin (CTK) [4–6], abscisic acid (ABA) [6], and Jasmonic acid (JA) [7], are key regulators for tuberous formation and development. Auxin and CTK in combination are crucial for tuberous initiation and formation [4–5], while ABA and CTK are likely the main regulation hormones for the secondary thickening growth of tuberous roots [6]. MADS box gene IbMADS1 was proposed to be an important integrator in the JA and CTK-mediated initiation of tuberization [8], while another MADS box gene SRD1 may be involved in the auxin-mediated initial thickening growth of storage roots [4]. KNOXI genes (class I Knotted1-like homeobox genes) have a possible function in tuberous root development by regulating CTK levels in tuberous roots of sweet potatoes [9]. IbEXP1 (an expansin gene) inhibits the initial thickening growth of tuberous roots by suppressing the proliferation of metaxylem and cambium cells [10]. Sucrose is an important factor regulating tuberous root formation and development [11–13]. Sucrose synthase (SuSy) was reported to play an important role in regulating tuberous development [12, 14]. Many environmental factors such as photoperied [15], drought [16], calcium [17], and hypoxia [18] also impact tuberous root development.
Expression studies have been carried out to explore the underlying molecular mechanisms of tuberous root development, based on a variety of techniques, such as simplified differential display analysis [19], suppression subtractive hybridization technology [3], microarray technology [20], and next generation sequencing and large-scale transcriptome analysis [14, 16, 20–23]. According to the expression studies, thousands of genes belonging to hundreds of biological pathways have been identified in tuberous vegetables, such as sweet potatoes [16, 21–23] and radishes [14], indicating that complex molecular networks and pathways are involved in the tuberous root formation and development. Of these pathways, metabolism pathways are likely the dominant ones [14, 21, 23]. Although remarkable progress has been made in the past, the molecular mechanisms are yet to be understood.
MicroRNAs (miRNAs) are a class of 21–24 nucleotide (nt) non-coding RNAs that regulate gene expression at the post-transcriptional level in many important biological processes, including root development [24–28]. It was reported that miR156 and miR172 are involved in the tuberization process of stem in potatoes (Solanum tuberosum ssp. andigena) [29–30]. Whether miRNAs are involved in the regulation of tuberous root development is still unclear.
Turnip is an important subspecies of Brassica rapa with the A genome (2n = 2x = 20, AA), just like Chinese cabbage [31]. Turnip has a modified tuberous root, while Chinese cabbage has a modified leafy head. Both of them are a good system for studying modified organ formation and development. In addition, the genome sequencing of Brassica rapa has been completed and the genome size is relatively small and simple, offering an excellent model to study the molecular mechanisms on the whole genome scale. Large-scale transcriptome analyses have been performed on the modified leafy head development of Chinese cabbage at both gene [32] and miRNA levels [33]. However, such work has not been done on the tuberous root of turnips. In this study, we analyzed both mRNAs and miRNAs expression profiles of tuberous roots at different seedling stages in turnips using high-throughput sequencing technology. We identified differentially expressed genes (DEGs), and miRNAs (DEMs) among various developing stages of tuberous roots. The study will lead to a better understanding of the underlying molecular mechanisms that regulate tuberous formation and development in the Brassica species.
Materials and Methods
Plant materials and growth conditions
The turnip cultivar “Chang Huang Man Jing” was used in this study. The seeds were carefully selected and sown in soil in square pots (8cm in width, length and height, 2 plants/pot) in a greenhouse at 20 ± 2°C under a 16 h light/8 h dark cycle. The plants were watered once a week with the Hoagland solution. All the plants were grown in the same room until collected for sample preparation. Samples consisting of hypocotyl and main root tissues were collected on day 18 (the early stage before cortex splitting, ES), day 28 (the stage of cortex splitting, CSS) and day 64 (the stage of root thickening, RTS) after sowing, respectively. The samples collected from at least 10 seedlings at each stage were pooled together and stored in liquid nitrogen for RNA extraction. All samples were harvested at the same time during the day. To reduce error rates, two independent biological replicates were used for each stage. The samples were labeled as ES1, CSS1, RTS1 for the first biological replicate, and ES2, CSS2, RTS2 for the second biological replicate, respectively.
RNA isolation
Total RNA was extracted from each sample using Trizol reagent (Invitrogen, Carlsbad CA, USA). The quality of RNA was verified using a NanoDrop ND-1000 Spectrophotometer (Nano-Drop, Wilmington, DE, USA) and a 2100 Bioanalyzer RNA Nanochip (Agilent, Santa Clara, CA, USA), and the quantity of total RNA was determined using the NanoDrop ND-1000 Spectrophotometer. RNA samples were selected for further cDNA library and small RNA (sRNA) library construction according to the following selection criteria: OD260/280≥1.8, OD260/230≥1.8, RIN≥9.0, 28S/18S≥1.6, a concentration of ≥200 ng/μl, and total mass of ≥30 μg. For mRNA sequencing, the total RNA was cleaned up using the RNeasy Plant Mini Kit (74904; Qiagen, Valencia, CA, USA). For sRNA sequencing, RNA fragments of 18–30 nt were isolated and purified from total RNA samples extracted using Trizol reagent after electrophoresing on a denaturing 15% polyacrylamide gel.
Preparation, sequencing and data processing of mRNA libraries
A total of six mRNA libraries were constructed in this study. The preparation and sequencing of the libraries were performed at Beijing Genomics Institute (BGI, Shenzhen, China) according to the manufacturer's protocol (Illumina). The sample libraries were single end sequenced via an IlluminaHiSeqTM 2000 (Illumina Inc., San Diego, CA, USA). Raw sequence reads were generated and filtered using the Illumina pipeline. Adapters, empty reads, and low quality reads (reads in which the percentage of unknown bases (N) greater than 10%, and reads in which the percentage of the low quality bases with a quality value ≤ 5 greater than 50%) were removed. The high-quality, clean sequences were mapped against the B. rapa (Chiifu-401) reference genome (http://brassicadb.org/brad/) using SOAPaligner/soap2 [34]. No more than two mismatches were allowed in the sequence alignment. The clean tags mapped to multiple reference genes were filtered, and the remaining tags were designed as unambiguous clean tags. Quality assessment of reads, statistics of alignment, sequencing saturation analysis and randomness assessments were carried out subsequently to assess the quality of sequencing.
We deposited all the sequence data of the six mRNA libraries into NCBI Sequence Read Archive (http://www.ncbi.nlm.nih.gov/sra/), and the accession numbers are SRX850781, SRX852388, SRX852389, SRX852516, SRX852517, and SRX852518 for ES1, ES2, CSS1, CSS2, RTS1, and RTS2 libraries, respectively.
Identification of DEGs among different seedling stages in the tuberous root of turnips
The gene expression level was calculated using the Reads Per kb per Million reads (RPKM) method [35]. After the gene expression levels of all samples were normalized to RPKM, genes that were differentially expressed among samples were screened according to a strict algorithm [36]. The threshold with a FDR (False Discovery Rate) of ≤0.001 and the absolute value of log2Ratio ≥ 1 was used to assess the significant difference of expression levels of a given gene between two seedling stages. To reduce the false positive rates, only DEGs detected in both biological replicates were considered. To understand the biological functions of DEGs involved in tuberous root development of turnips, Gene Ontology (GO) enrichment analysis using the GO database (http://www.geneontology.org/), and pathway enrichment analysis using the KEGG database (http://www.genome.jp/kegg/) were conducted subsequently.
Preparation, sequencing and data processing of sRNA libraries
A total of six sRNA libraries were constructed and sequenced at Beijing Genomics Institute (BGI, Shenzhen, China). The sRNA Libraries were prepared following the procedure by Wang et. al (2013) [33], and were sequenced on an Illumina HiSeqTM 2000. Raw samples of the 49 nt RNA tags obtained from the HiSeq sequencing were cleaned by removing the low quality tags (reads with more than one base with SQ value ≤10 or more than two bases with SQ values ≤13), 5' adaptor contaminants, oversized insertions, no insert tags, poly A tags and small tags (reads shorter than 18 nt). Length distribution of the remaining clean tags was analyzed to determine the compositions of sRNA samples. Subsequently, we mapped the final clean reads against the Chiifu Chinese cabbage genome (http://brassicadb.org/brad/) using SOAP v1.11 to determine gene expression levels and distributions in the genome [34]. The annotation of the clean sRNA tags and prediction of novel miRNAs were performed following the procedure by Wang et. al (2013) [33]. Putative target genes of the miRNAs were predicted using the rules developed by Allen et al. (2005) [37] and Schwab et al. (2005) [38].
The reads annotated as rRNA, scRNA, snoRNA, snRNA, and tRNA from GenBank and Rfam were removed, while those annotated as conserved miRNAs from the miRBase database and those predicted as novel miRNAs by Mireap (http://sourceforge.net/projects/mireap/) were used for further analysis in their differential expression profiles among different seedling stages in tuberous roots of turnips.
All the sequence data of the six sRNA libraries were deposited into NCBI Sequence Read Archive (http://www.ncbi.nlm.nih.gov/sra/) under the accession numbers SRX857356, SRX856589, SRX856653, SRX856655, SRX856663 and SRX856672 for ES1, ES2, CSS1, CSS2, RTS1, and RTS2 libraries, respectively.
Identification of DEMs among different seedling stages in the tuberous root of turnips
The expression level of miRNA in all samples was normalized as transcripts per million (TPM, TPM = actual miRNA count / total count of clean reads × 106). The fold-change and P-value from the normalized expression level were calculated among different samples. The fold change was calculated as log2 Ratio, and the P-value was calculated using the formula described in Wang et. al (2013) [33]. The threshold with the P-value ≤ 0.05 and an absolute value of log2Ratio ≥ 1 was used to evaluate the significance of the difference of miRNA expression levels. To reduce the false positive rates, only DEMs were detected in both biological replicates were considered.
Quantitative RT-PCR analysis
To confirm the identified DEGs and DEMs, DEMs and their differentially expressed target genes were subjected to the quantitative RT-PCR (qRT-PCR) analysis. For DEM validation, the first strand cDNA was synthesized using the miRcutemiRNAcDNA Synthesis Kit (DP501, Tiangen, Beijing, China), and qRT-PCR analyses were performed using the miRcutemiRNAqPCR Detection Kit (SYBR Green) (FP401, Tiangen, Beijing, China). Spliceosomal RNA U6 was used as a constitutive expression control. PCR reactions were performed in a 96-well plate at 94°C for 2 min, followed by 45 cycles of reaction (94°C for 20 s, followed by 60°C for 34 s). For DEGs validation, the first strand cDNA was synthesized using M-MLV Reverse transcriptase (M1701, Promega, Madison, USA) and qRT-PCR analyses were performed using a PrimeScript™ RT reagent Kit (Perfect Real Time) (DRR037A, Takara, Dalian, China). PCR reactions were performed in a 96-well plate at 95°C for 10 min, followed by 40 cycles of reaction (95°C for 10 s, followed by 60°C for 60 s). Actin gene was used as a constitutive expression control. All qRT-PCR analyses were performed using an IQ5 Real-Time PCR System (BIO-RAD, Hercules, CA, USA). The specific primers for qRT-PCR are described in S1 Table.
Results
mRNA Library construction, sequencing and mapping against the B. rapa (Chiifu-401) reference genome
To study the global transcriptome during the tuberous root development of turnips, a total of six mRNA libraries (three seedling stages during tuberous root development, two biological replicates for each stage) were constructed and sequenced using an IlluminaHiSeqTM 2000. Approximately 12 million raw tags were generated for each library, and more than 99.0% of the raw tags were identified as clean tags (Table 1). We then mapped all the clean tags against the B. rapa reference genome (http://brassicadb.org/brad/), and found that over 71% of the clean tags for each library could be mapped onto the reference genes with ≤2 bp mismatch. Less than 29% of the clean tags for each library could not be aligned to any reference genes due to incomplete sequences, and these tags were designated as unknown tags. More than 67.3% (about 8 million) of the clean tags in each library were mapped to a single gene, while less than 4.7% (about 0.5 million) were mapped to multiple reference genes (Table 1). Unknown tags and tags mapped to multi-genes were filtered out, and the unique clean tags mapped to a single gene were reserved for further DEG analysis.
The saturation for the six libraries was analyzed, and the results showed that the number of detected genes increased with the number of reads until the number of reads reached up to 2 million (S1 Fig), which indicates that the sequencing depth is sufficient for the transcriptome coverage in this study. We also assessed the randomness of RNA fragmentation in the six libraries to estimate whether it was sufficient for a subsequent bioinformatics analysis of gene expression, and the results showed that the reads in each position of the reference gene were distributed evenly and demonstrated highly similar tendencies in all libraries (S2 Fig).
Differentially regulated genes in tuberous root development of turnips
Based on deep sequencing of the six mRNA libraries, approximately thirty thousand genes (about 70% of reference genes in Chinese cabbage) were detected in each library (Table 1, S2 Table). To compare gene expression differences in different stages of tuberous root development in turnips, the gene expression level was calculated using the RPKM method [35]. Gene expression correlations based on RPKM between the two biological replicates were analyzed, and the correlation coefficients (R2) were high (0.97, 0.95, and 0.96 for ES, CSS, and RTS, respectively) (Fig 1A). Genes expressed in both replicates in the same stage were screened. 28,798, 28,315, and 27,288 genes were expressed in both replicates in ES, CSS, and RTS, respectively. Most of them (25,729) were expressed in all three stages, while only a small portion (460–1,365) was expressed at two stages or specifically at one stage (Fig 1B).
(A) Correlation coefficients of two biological duplicates; (B) Venn diagram analyses of stage specific genes in tuberous root of turnips; (C) Venn diagram analyses of stage specific conserved miRNA in tuberous root of turnips; (D) Venn diagram analyses of stage specific novel miRNA in tuberous root of turnips; (E) Differentially expressed genes across all libraries; (F) Differentially expressed conserved miRNAs across all libraries; (G) Differentially expressed conserved novel miRNAs across all libraries.
We compared the gene expression profiles between every two stages in tuberous root development in both biological replicates (ES1 vs. CSS1, ES1 vs. RTS1, CSS1 vs. RTS1, ES2 vs. CSS2, ES2 vs. RTS2, and CSS2 vs. RTS2, where the former is the control and the latter is the experimental group). The results are presented in S2 and S3 Tables. DEGs detected in both biological replicates were screened for subsequent analysis. A total of 2693 (1176 up- and 1517 down-regulated), 4655 (1094 up- and 3561 down-regulated), and 2133 (436 up- and 1677 down-regulated) genes were differentially expressed in both replicates in the stage pairs, ES vs. CSS, ES vs. RTS, and CSS vs. RTS, respectively (Fig 1E). More genes, which were down-regulated, were detected than the up-regulated ones in all pairs (Fig 1E).
Gene Ontology (GO) enrichment analysis of the DEGs was performed to gain insight into the possible mechanisms involved in tuberous root development. All the well-annotated DEGs were functionally clustered into three main groups (cellular component, molecular function, and biological process). For cellular component ontology, a large number of the DEGs in the stage pair, ES vs. CSS, involved in cell, cell part, and organelle, accounted for 63.73%, 63.73%, and 44.15% of all the well-annotated DEGs, respectively. Under the molecular function category, catalytic activity (50.26%) and binding (48.62%) were the main groups, followed by nucleic acid binding transcription factor activity (10.54%) and transporter activity (6.88%) in ES vs. CSS. For biological process category, the main groups of the DEGs in ES vs. CSS were involved in metabolic process (53.74%), cellular process (51.42%), response to stimulus (39.89%), and single-organism process (32.31%). Similar results were obtained for GO enrichment analysis of the DEGs in ES vs. RTS, and CSS vs. RTS (Fig 2).
(A) Gene classification based on GO for DEGs. Top six classes are shown for cellular component, molecular function and biological process, respectively; (B) Gene classification based on GO for target genes of DEMs. Top six classes are shown for cellular component, molecular function and biological process, respectively.
We performed pathway enrichment analysis by mapping all the annotated genes to terms in the KEGG database to identify significantly altered biological pathways associated with the tuberous root development in turnips. The results are shown in Table 2. A total of 17, 20, and 21 pathways were significantly enriched with Q value ≤0.05 in ES vs. CSS, ES vs. RTS, and CSS vs. RTS, respectively. Of these, 11 pathways were significantly enriched in all comparison pairs, indicating that these pathways may be involved in both tuberous root initiation and the secondary thickening process. Four pathways, including plant hormone signal transduction, amino sugar and nucleotide sugar metabolism, other glycan degradation, and cutin, suberine and wax biosynthesis, were only enriched in ES vs. CSS, implying the great importance of these pathways in tuberous root initiation. The plant hormone signal transduction pathways were significantly enriched in ES vs. CSS. Among these seven phytohormones, namely auxin, CTK, brassinosteroid (BR), gibberellin (GA), JA, ethylene, and ABA, signal transduction pathway of auxin had the most DEGs (S4 Table). Metabolic pathways, biosynthesis of secondary metabolites, and plant hormone signal transduction were the top three enriched pathways in ES vs. CSS, while metabolic pathways, biosynthesis of secondary metabolites, and starch and sucrose metabolism were the top three enriched pathways in both ES vs. RTS and CSS vs. RTS (Table 2).
sRNA Library construction, sequencing and mapping against the B. rapa reference genome
sRNA libraries from the six samples of turnips were sequenced using an Illumina HiSeqTM 2000 to reveal the global expression profiles of miRNAs involved in tuberous root development of turnips. Over ten million raw reads were generated from each library. More than 99.4% of the raw reads were identified as clean reads in each library. A total of 3,819,650, 4,555,962, 4,164,927, 4,090,180, 4,665,954, and 4,130,002 unique sRNAs were obtained after collapsing redundancies from the clean reads of the ES1, CSS1, RTS1, ES2, CSS2, and RTS2 libraries, respectively (Table 3). The length distribution of the clean reads in each library was summarized, and the results showed that 24 nt sRNAs were the most abundant in all six sRNA libraries, accounting for 40.53%-51.54% of all sRNAs. The 21 nt class was the second most abundant in the ES1, RTS1, and ES2 libraries, while the 23 nt class was the second most abundant in the CSS1, CSS2, and RTS2 libraries (S3 Fig). We mapped the final clean sRNAs against the reference Chinese cabbage genome (http://brassicadb.org/brad/), and found that 66.72%- 68.96% of the total clean sRNAs and 50.27%- 51.72% of the unique sRNAs of the six libraries could be matched with the reference genome (Table 3). The percentages of mapped total clean sRNAs and unique sRNAs were similar in all of the six sRNA libraries.
Identification of conserved and novel miRNAs in the tuberous root of turnips
The clean sRNAs were annotated following the procedure described in Wang et al. (2013) [33], and the distribution of sRNAs among different categories is shown in S5 Table. A total of 17,465 (0.46%), 21,972 (0.48%), 14,890 (0.36%), 13,422 (0.33%), 16,983 (0.36%), and 14,769 (0.36%) candidate miRNAs were identified from the ES1, CSS1, RTS1, ES2, CSS2, and RTS2 libraries, respectively (S5 Table). We aligned all these miRNA candidates to the miRNA precursor/mature miRNA of all plants deposited in the miRBase database to identify conserved miRNAs. Only perfectly matched or closely related (≤ 2 mismatches) sequences and precursors that could form hairpins were defined as conserved. A total of 173, 187, 162, 165, 179, and 157 conserved miRNAs were identified from the ES1, CSS1, RTS1, ES2, CSS2, and RTS2 libraries, respectively. In all, 279 conserved miRNAs were identified in this study, among which, 98 were expressed in every library (S6 Table). Conserved miRNAs expressed in both replicates in the same stage were screened for further analysis. It was found that 139, 149, and 125 were expressed in both replicates in ES, CSS, and RTS, respectively (Fig 1C). One hundred of the total detected conserved miRNAs were expressed in all three stages, while 2–18 were only expressed at two stages or specifically at one stage (Fig 1C).
Currently, 43 miRNAs (35 miRNAs, 4 miRNA-5p and 4 miRNA-3p) in B. rapa have been deposited in the miRBase database (miRNA Registry, Release 20.0, 2013). In this study, 38, 36, 31, 39, 35, and 31 bra-miRNAs were identified from the ES1, CSS1, RTS1, ES2, CSS2, and RTS2 libraries, respectively (S7 Table). Collectively, 42 of these 43 bra-miRNAs deposited in the miRBase database, except bra-miR5723, were identified in at least one of the six libraries; 27 bra-miRNAs were found in each of the six libraries. However, over half (23–26) of these 43 bra-miRNAs (e.g. bra-miR5720, bra-miR5717, and bra-miR5724) showed a very low level (less than 100 reads) or no expression in all libraries, while only 3–6 (e.g. bra-miR157a and bra-miR164a) were expressed dramatically high (more than 5000 reads) in all or most of the libraries (S7 Table). This finding indicates that the expression of conserved miRNAs varies significantly in the tuberous root of turnips. Bra-miRNAs expressed in both replicates in the same stage were screened, and 36, 34, and 29 were expressed in both replicates in ES, CSS, and RTS, respectively (S7 Table).
Novel miRNAs were predicted by exploring the secondary structure, the Dicer cleavage site and the minimum free energy of the unannotated sRNA tags which could be mapped to the genome using the Mireap program developed by BGI (http://sourceforge.net/projects/mireap/). Accordingly, 292, 292, 279, 266, 287, and 266 novel miRNAs were identified in the ES1, CSS1, RTS1, ES2, CSS2, and RTS2 libraries, respectively (S6 Table). Most of the novel miRNAs (more than 85%) showed very low expression level (less than 100 reads), and only three (novel_mir_14, novel_mir_176, and novel_mir_26) showed an expression level higher than 1000 reads in all six libraries. A total of 912 novel miRNAs were identified, but only 48 were expressed in all six libraries (S6 Table), indicating significant variation of novel miRNA expression in different samples. A total of 119, 128, and 121 were expressed in both replicates in ES, CSS, and RTS, respectively (Fig 1D). Venn diagram analyses of stage specific novel miRNAs in tuberous root of turnips are shown in Fig 1D.
Differentially regulated miRNAs in tuberous root development of turnips
miRNA expression correlations between the two biological replicates were analyzed according to TPM. The correlation coefficients (R2) were 0.98, 0.98, and 0.98 for ES, CSS, and RTS, respectively (Fig 1A). Based on the criteria described in the Methods section, DEMs between different stages in each replicate were identified, and the detailed results are listed in S6 Table. DEMs detected in both the biological replicates were screened for subsequent analysis. In all, 46 (29 up- and 17 down-regulated), 50 (22 up- and 28 down-regulated) and 36 (11 up- and 25 down-regulated) conserved miRNAs (Fig 1F), and 53 (26 up- and 27 down-regulated), 62 (26 up- and 36 down-regulated) and 40 (12 up- and 28 down-regulated) novel miRNAs (Fig 1G) were differentially expressed in both replicates in ES vs. CSS, ES vs. RTS and CSS vs. RTS, respectively. Eight miRNAs (miR156a, miR157a, miR157a-3p, miR5718, novel_mir_4, novel_mir_85, novel_mir_142, and novel_mir_236) were differentially expressed in both replicates of all comparison pairs, indicating that these miRNAs may be involved in both tuberous root initiation and the secondary thickening process. Twenty miRNAs were differentially expressed in ES vs. CSS, implying a potential role of these miRNAs in tuberous root initiation (S6 Table).
We found 128 (73 for 20 conserved DEMs and 67 for 22 novel DEMs), 163 (97 for 26 conserved DEMs and 66 for 23 novel DEMs), and 97 (53 for 12 conserved DEMs and 44 for 13 novel DEMs) putative target genes for DEMs in ES vs. CSS, ES vs. RTS, and CSS vs. RTS, respectively (Table 4, S8 Table). GO enrichment analyses of the target genes of DEMs (both conserved and novel DEMs) detected in both replicates were conducted subsequently. For cellular component ontology, cell, cell part, and organelle were the main groups of the target genes; under the molecular function category, binding, nucleic acid binding transcription factor activity, and catalytic activity were the main groups; for the biological process category, cellular process, metabolic process, and single-organism process were the top three groups (Fig 2B). The results of GO analysis in the target genes of DEMs were similar to those of DEGs in this study (Fig 2).
Pathway enrichment analyses of the target genes of DEMs detected in both replicates were also performed using the KEGG pathway database. No pathway was significantly enriched with Q value of ≤0.05 in all the three comparisons, which may be because of the small numbers of the target genes. More genes were enriched in plant hormone signal transduction pathway in ES vs. CSS and ES vs. RTS, while more were enriched in metabolic pathways in CSS vs. RTS (S9 Table).
The expression profiles of the target genes of DEMs were analyzed according to the high-throughput sequencing of mRNA. The results showed that a number of DEMs had target gene(s) differentially expressed in the DEG analysis. Some differentially expressed target genes had their expression profiles matched up with those of their corresponding DEMs (S8 Table). These DEGs and DEMs are the most important candidates for further study of their functions in tuberous root development of turnips.
Validation of DEGs and DEMs by qRT-PCR
The expression profiles of six DEMs between different stages and their seven differentially expressed target genes were verified by qRT-PCR. The results showed that all seven genes and four miRNAs (miR156a, miR157a, miR172a, and novel-mir-208) exhibited the same expression tendency as the original results (Fig 3, S2 and S6 Tables). The expression profile of miR395b did not match with that of the original results (Fig 3B, S6 Table), probably due to its low expression level. Similarly the expression profile of novel-mir-108 did not match either. Of the four miRNAs exhibiting the same expression tendency as the original results, miR156a, miR157a, and miR172a exhibited relatively high expression level in all samples, while novel-mir-208 expressed extremely low in all samples in the sRNA sequencing.
(A) qRT-PCR validation of DEGs from high-throughput RNA sequencing analyses. Relative expression levels were calculated using actin2 as the reference gene using the formula 2-△△Ct. The values indicate means of three biological replicates ± standard error. (B) qRT-PCR validation of DEMs from high-throughput sRNA analyses. Relative expression levels were calculated using spliceosomal RNA U6 as the reference gene using the formula2-△△Ct. The values indicate means of three biological replicates ± standard error.
Discussion
Quality assessment of sequencing
Tuberous root development is a complex biological process controlled by interactions of genetic, physiological and environmental factors. High-throughput sequencing, owing to its advantages of high coverage and low cost, has been widely used to study complex biological processes of plants in recent years. Here, in order to gain insight into the underlying molecular mechanisms regulating tuberous root development, large-scale analyses of both mRNA and miRNA expression profiles of tuberous roots at different seedling stages in turnips were performed by high-throughput sequencing technology. Good repeatability was achieved, according to the high correlation coefficients (R2>0.95, Fig 1A) between the two replicates on both mRNA and miRNA expression levels. The number of raw reads, the percentages of clean tags, and the percentages of clean tags mapped onto the reference genome in both mRNA and miRNA expression profile analyses were greater than, or similar to, those in previous studies in the brassica species [33, 39–44]. According to the results of sequencing saturation analysis and randomness assessment, the clean tags obtained in this study were sufficient for subsequent bioinformatics analyses of gene expression. The size distribution of sRNAs in each sRNA library was also similar with those in previous studies in Chinese cabbage [33, 42–47]. However, there were some small differences between different libraries. The 21nt sRNAs were the second most abundant class in CSS libraries, while the 23 nt sRNAs were the second most abundant class in RTS libraries [S3 Fig], implying the stage specific functions of these sRNAs in the process of tuberous root development in turnips. In ES libraries, the second major classes of sRNAs were not the same in the two replicates [S3 Fig], which is probably due to sequencing errors.
Genes and biological pathways involved in tuberous root development of turnips
Genome-wide gene expression profiles analysis, based on high-throughput RNA sequencing, could facilitate the revelation of quantitative changes in transcript abundance and the identification of candidate genes. We identified a large number of DEGs among different developing stages of tuberous roots of turnips (Fig 1E, S2 Table), implying a high complexity of molecular regulations in tuberous root development.
Tuberous root is a nutrition storage organ. Sucrose and starch, proteins and amino acid, dietary fiber, vitamin, and other relevant molecules are accumulated in the tuberization process [1, 11, 13]. Many metabolic pathways may be actively involved in this process. Based on our genome-wide DEGs analysis, metabolic and biosynthesis of secondary metabolite pathways were the top two pathways with a large number of DEGs significantly enriched, which is consistent with previous studies on the tuberous root of sweet potatoes [16, 21–23], radishes [14], and cassavas [20]. Of these metabolism pathways, starch and sucrose metabolism pathway and phenylpropanoid biosynthesis pathway are the top two evident pathways, which may be active in both tuberous root initiation and the secondary thickening process. In cassavas, glycolysis/gluconeogenesis was identified to be the most evident pathway involved in tuberous root development [20]. In turnips, this pathway was probably not so actively involved, because we did not see a significant enrichment.
Phytohormones play important roles in the regulation of tuberous formation and development [4–7]. Auxin [4–5], CTK [4–5, 7], and JA [7] were reported to be crucial regulators for tuberous initiation, while auxin [4], CTK [6], and ABA [6] probably regulate the secondary thickening growth of tuberous roots. In our study, plant hormone signal transduction pathways were significantly enriched in ES vs. CSS, but not in ES vs. RTS or CSS vs. RTS (Table 2). This suggests that DEGs related to plant hormone signal transduction may play a more important role in tuberous root initiation than in the secondary thickening growth of tuberous roots in turnips. Besides DEGs involved in auxin, CTK and JA signal transduction pathways, DEGs involved in BR, GA, and ethylene signal transduction pathways were also identified in ES vs. CSS (S4 Table), indicating potential roles of these plant hormones in regulating tuberous root initiation in turnips. Furthermore, over half of the DEGs involved in plant hormone signal transduction pathways in ES vs. CSS belong to auxin and CTK pathways, implying the predominant roles of these two phytohormones in regulating tuberous root initiation in turnips.
Several genes related to tuberous development have been identified in previous studies. MADS box genes, such as IbMADS1, SRD1, and IbAGL17, were identified to be preferentially expressed in tuberous root of sweet potatoes [4, 8, 48]. In this study, four homologous of MADS box genes (Bra017340, Bra004928, Bra006051, and Bra009913) were identified to be possibly involved in tuberous root development in turnips. Bra017340 was more predominantly expressed in CSS libraries than in ES and RTS libraries, indicating a possible positive role of Bra017340 in tuberous initiation process. Bra004928 and Bra006051 were up-regulated in ES vs. RTS, which are likely to be positively involved in the secondary thickness process of tuberous root development in turnips. Bra009913 were down-regulated in both ES vs. RTS and CSS vs. RTS, implying its negative role in regulating tuberous root development. Expansin (EXP1) has been reported to play a negative role in the formation of storage root in sweet potatoes [10]. In our study, the homologous gene EXP1 (Bra007864) and five other Expansin genes (Bra012218, Bra036225, Bra028878, Bra025372, and Bra036499) were identified to be down-regulated in both ES vs. CSS and ES vs. RTS. Previous studies reported SuSy plays a key role in regulating tuberous root development [12], and homologous of the enzyme are up-regulated in tuberous roots of radishes [14]. In this study, one homologous of SuSy 3 was up-regulated, and two homologous of SuSy 6 were down-regulated in ES vs. RTS. No SuSy homologous genes were identified differentially expressed in ES vs. CSS. This indicates that SuSy genes probably are involved in tuberous root thickness, but not in tuberous root initiation. New members of MADS box genes, expansin genes and SuSy were identified for tuberous root development, indicating that high-throughput RNA sequencing is effective in identifying DEGs.
miRNAs involved in tuberous root development of turnips
miRNAs are important regulators of many important biological processes, such as development, hormone signaling, abiotic stress responses, and pathogen responses [49–51]. They negatively regulate their target gene expression at the post-transcriptional level. Large-scale miRNA expression profiles were analyzed on the modified leafy head development of Chinese cabbage, providing a relatively complete view of the miRNAs involved in the modified leafy head development of B. rapa [33]. Whether miRNAs are involved in the modified tuberous root development in B. rapa is still unclear. According to our results, many miRNAs were found differentially expressed in different development stages of tuberous roots in turnips [S6 Table]. In previous studies, it was reported that miR156 and miR172 are involved in tuberous stem development process in potatoes. Overexpression of miR156 reduces the yield of tuberous stem [29], while miR172 overexpression accelerates the tuberization in potatoes [30]. miR156 down-regulates the expression of miR172 by targeting one of the squamosa promoter binding protein-like (SPL) genes, StSPL9, a transranscription factor of miR172 [29]. In this study, we found miR156a was down-regulated in the root tuberization process, while miR172a was up-regulated in CSS and RTS compared with ES, which is coincide with that in stem tuberization process in potatoes. miR157 families are predicted to target the SPL genes at the same binding sites as miR156 families [52]. In our study, similar to miR156a, the expression of miR157a was down-regulated during the tuberization process. They are likely to play the same roles in regulating tuberous root development in turnips. Furthermore, miR156a, miR157a and miR172a exhibited a relatively high expression level in all samples, and their expression patterns were negatively correlated with all or some of their target DEGs. It indicated that these three miRNAs may play important roles in tuberous root development of turnips. Further studies need to be done to fully understand their functions in regulating tuberous root development.
Gene Ontology analyses of the target genes of the DEMs were similar with those of the DEGs (Fig 2). Although no KEGG pathway was significantly enriched due to the small numbers of the target genes, more genes were clustered in the plant hormone signal transduction pathway in ES vs. CSS and ES vs. RTS, and more were enriched in metabolic pathways in CSS vs. RTS (S9 Table). It supports our hypothesis that plant hormone signal transduction pathway plays more important roles in tuberous root initiation.
Since miRNAs negatively regulate gene expression at the post-transcriptional level, the expression profiles of DEMs are expected to be negatively correlated with the target gene expression profiles. Previous studies have shown that some miRNAs were co-expressed with their target genes [53–55]. We analyzed the expression profiles of the putative target genes of DEMs based on the DEGs sequencing, and found that many miRNA expression profiles were not correlated with their target gene expressions. Similar results were obtained in Chinese cabbage for the modified leafy head development [33]. The results indicate that target gene expression is not regulated by miRNAs alone, but many other factors.
qRT-PCR validation of the six DEMs showed that four DEMs (miR156a, miR157a, miR172a, and novel-mir-208) exhibited the same expression tendency as the original results, while two (miR395b, and novel-mir-108) did not match with those of the original results from the high-throughput sRNA sequencing (S6 Table, Fig 3B). This is probably due to the extremely low expression levels (less than 300 reads were detected in each sample) of the two miRNAs detected in the sRNA sequencing.
Conclusions
We analyzed both mRNA and miRNA expression profiles of the tuberous root of turnips at different stages based on high-throughput sequencing technology. Our results indicate that a large number of DEGs and DEMs are involved in tuberous root development. According to the DEG analysis, we propose that metabolism is the dominant pathway in both tuberous root initiation and secondary thickening process; the plant hormone signal transduction pathway may play a more important role in regulating tuberous root initiation; and the starch and sucrose metabolism may be more important for the secondary thickening process. Our analyses for DEMs partially supported these hypotheses. Based on the DEM analyses, miR156a, miR157a, and miR172a may play important roles in tuberous root development in turnips. These findings will contribute to our further understanding of the molecular mechanisms of tuberous root development in turnips.
Supporting Information
S1 Fig. Saturation analysis of the six mRNA libraries.
A, B, C, D, E, and F were saturation analysis of ES1, ES2, CSS1, CSS2, RTS1, and RTS2 library, respectively.
https://doi.org/10.1371/journal.pone.0137983.s001
(TIFF)
S2 Fig. Randomness assessments of the six mRNA libraries.
A, B, C, D, E, and F were randomness assessments of ES1, ES2, CSS1, CSS2, RTS1, and RTS2 library, respectively.
https://doi.org/10.1371/journal.pone.0137983.s002
(TIFF)
S1 Table. qRT-PCR primers for seven DEGs and six DEMs.
https://doi.org/10.1371/journal.pone.0137983.s004
(XLSX)
S2 Table. Genes expressed in each mRNA library.
The genes with Log2 ratio were DEGs in the corresponding comparison.
https://doi.org/10.1371/journal.pone.0137983.s005
(XLSX)
S3 Table. Statistics of differentially expressed genes and miRNAs across all libraries in each replication.
https://doi.org/10.1371/journal.pone.0137983.s006
(XLSX)
S4 Table. Plant hormone signal transduction pathways detected in the study.
https://doi.org/10.1371/journal.pone.0137983.s007
(XLSX)
S5 Table. Distribution of sRNAs among different categories in the six libraries.
https://doi.org/10.1371/journal.pone.0137983.s008
(XLSX)
S6 Table. miRNAs identified in each sRNA libraries.
The genes with Log2 ratio were DEGs in the corresponding comparison. Novel miRNAs shown in bold are ones having miRNA* identified in the study.
https://doi.org/10.1371/journal.pone.0137983.s009
(XLSX)
S7 Table. Bra-miRNAs identified in each sRNA libraries.
https://doi.org/10.1371/journal.pone.0137983.s010
(XLSX)
S8 Table. DEMs and their target gene.
Target genes shown in bold red are DEGs with their expression profiles matched up with their corresponding DEMs, while those shown in bold blue are DEGs with their expression profiles not matched up with their corresponding DEMs.
https://doi.org/10.1371/journal.pone.0137983.s011
(XLSX)
S9 Table. Pathway enrichment analysis of target genes of DEMs detected in both replicates (top three pathways are shown for each comparison).
https://doi.org/10.1371/journal.pone.0137983.s012
(XLSX)
Acknowledgments
We thank Professor Cuihua Wang for providing the plant materials, and for advice on manuscript preparation.
Author Contributions
Conceived and designed the experiments: JJL JWG. Performed the experiments: JJL QD. Analyzed the data: JJL QD FDW. Contributed reagents/materials/analysis tools: YHZ HYL JWG. Wrote the paper: JJL FDW JWG.
References
- 1. Dini I, Tenore GC, Dini A. Saponins in Ipomoeabatatas tubers: Isolation, characterization, quantification and antioxidant properties. Food Chemistry 2009; 113: 411–9.
- 2. Curtis IS. The noble radish: past, present and future. Trends in Plant Science 2003; 8: 305–307. pmid:12878009
- 3. Sun P, Guo Y, Qi J, Zhou L, Li X. Isolation and expression analysis of tuberous root development related genes in Rehmannia glutinosa. Molecular Biology Reports 2010; 37: 1069–1079. pmid:19774491
- 4. Noh SA, Lee HS, Huh EJ, Huh GH, Paek KH, Shin JS, et al. SRD1 is involved in the auxin-mediated initial thickening growth of storage root by enhancing proliferation of metaxylem and cambium cells in sweetpotato (Ipomoea batatas). Journal of Experimental Botany 2010; 61: 1337–1349. pmid:20150515
- 5. Fan M, Liu Z, Zhou L, Lin T, Liu Y, Luo L. Effects of plant growth regulators and saccharide on in vitro plant and tuberous root regeneration of Cassava (Manihot esculenta Crantz). Journal of Plant Growth Regulation 2011; 30: 11–19.
- 6. Wang Q, Zhang L, Guan Y, Wang Z. Endogenous hormone concentration in developing tuberous roots of different sweetpotato genotypes. Agricultural Science in China 2006; 5: 919–927.
- 7. Yasunori K, Keiko OT, Yoshio K. Stimulation of root thickening and inhibition of bolting by Jasmon Acid in beet plants. Plant Production Science 2001; 4: 131–135.
- 8. Ku AT, Huang YS, Wang YS, Ma D, Yeh KW. IbMADS1 (Ipomoea batatas MADS-box 1 gene) is involved in tuberous root initiation in sweetpotato (Ipomoea batatas). Annals of Botany 2008; 102: 57–67. pmid:18463111
- 9. Tanaka M, Kato N, Nakayama H, Nakatani M, Takahata Y. Expression of class 1 Knotted1-like homeobox genes in the storage roots of sweetpotato (Ipomoea batatas). Journal of Plant Physiology 2008; 165: 1726–1735. pmid:18242774
- 10. Noh SA, Lee HS, Kim YS, Paek KH, Shin JS, Bae JM. Down-regulation of the IbEXP1 gene enhanced storage root development in sweetpotato. Journal of Experimental Botany 2013; 64: 129–142. pmid:22945944
- 11. Gupta AK, Singh J, Kaur N. Sink development, sucrose metabolising enzymes and carbohydrate status in turnip (Brassica rapa L.). Acta Physiologiae Plantarum 2001; 23: 31–36.
- 12. Rouhier H, Usuda H. Spatial and temporal distribution of sucrose synthase in the radish hypocotyl in relation to thickening growth. Plant and Cell Physiology 2001; 42: 583–593. pmid:11427677
- 13. Luo X, Huang Q. Relationships between leaf and stem soluble sugar content and tuberous root starch accumulation in cassava. Journal of Agricultural Science 2011; 3: 64–72.
- 14. Wang S, Wang X, He Q, Liu X, Xu W, Li L, et al. Transcriptome analysis of the roots at early and late seedling stages using Illumina paired-end sequencing and development of EST-SSR markers in radish. Plant Cell Reports 2012; 31: 1437–1447. pmid:22476438
- 15. Legnani G, Miller WB. Short photoperiods induce fructan accumulation and tuberous root development in Dahlia seedlings. New Phytologist 2001; 149: 449–454.
- 16. Solis J, Villordon A, Baisakh N, LaBonte D, Firon N. Effect of drought on storage root development and gene expression profile of sweetpotato under greenhouse and field conditions. Journal of the American Society for Horticultural Science 2014; 139: 317–324.
- 17. Yao Y, Min Y, Geng MT, Wu XH, Hu XW, Fu SP, et al. The effects of calcium on the in vitro cassava storage root formation. Advanced Materials Research 2013; 726–731: 4529–4533.
- 18. Eguchi T, Yoshida S. Tuberous root thickening: a rapid response to hypoxia in sweetpotato (Ipomoea batatas (L.) Lam.). Environment Control in Biology 2011; 49: 47–50.
- 19. Tanaka M, Takahata Y, Nakatani M. Analysis of genes developmentally regulated during storage root formation of sweetpotato. Journal of Plant Physiology 2005; 162: 91–102. pmid:15700424
- 20. Yang J, An D, Zhang P. Expression profiling of cassava storage roots reveals an active process of glycolysis/gluconeogenesis. Journal of Integrative Plant Biology 2011; 53: 193–211. pmid:21205184
- 21. Wang Z, Fang B, Chen J, Zhang X, Luo Z, Huang L, et al. De novo assembly and characterization of root transcriptome using Illumina paired-end sequencing and development of cSSR markers in sweetpotato (Ipomoea batatas). BMC Genomics 2010; 11: 726. pmid:21182800
- 22. Xie F, Burklew CE, Yang Y, Liu M, Xiao P, Zhang B, et al. De novo sequencing and a comprehensive analysis of purple sweetpotato (Impomoea batatas L.). transcriptome. Planta 2012; 236: 101–113. pmid:22270559
- 23. Firon N, LaBonte D, Villordon A, Kfir Y, Solis J, Lapis E, et al. Transcriptional profiling of sweetpotato (Ipomoea batatas) roots indicates down-regulation of lignin biosynthesis and up-regulation of starch biosynthesis at an early stage of storage root formation. BMC Genomics 2013; 14: 460. pmid:23834507
- 24. Khan GA, Declerck M, Sorin C, Hartmann C, Crespi M, Lelandais-Brière C. MicroRNAs as regulators of root development and architecture. Plant Molecular Biology 2011; 77:47–58. pmid:21607657
- 25. Sorin C, Declerck M, Christ A, Blein T, Ma L, Lelandais-Brière C, et al. A miR169 isoform regulates specific NF-YA targets and root architecture in Arabidopsis. New Phytologist 2014; 202:1197–211. pmid:24533947
- 26. Guo HS, Xie Q, Fei JF, Chua NH. MicroRNA164 directs NAC1 mRNA cleavage to downregulate auxin signals for lateral root development. The Plant Cell 2005; 17:1376–1386. pmid:15829603
- 27. Wang JW, Wang LJ, Mao YB, Cai WJ, Xue HW, Chen XY. Control of root cap formation by microRNA-targeted auxin response factors in Arabidopsis. The Plant Cell 2005; 17:2204–2216. pmid:16006581
- 28. Yoon EK, Yang JH, Lim J, Kim SH, Kim SK, Lee WS. Auxin regulation of the microRNA390-dependent transacting small interfering RNA pathway in Arabidopsis lateral root development. Nucleic Acids Research 2010; 38:1382–1391. pmid:19969544
- 29. Bhogale S, Mahajan AS, Natarajan B, Rajabhoj M, Thulasiram HV, Banerjee AK. MicroRNA156: A potential graft-transmissible microRNA that modulates plant architecture and tuberization in Solanum tuberosum ssp. Andigena. Plant Physiology 2014; 164:1011–1027. pmid:24351688
- 30. Martin A, Adam H, Díaz-Mendoza M, Zurczak M, González-Schain ND, Suárez-López P. Graft-transmissible induction of potato tuberization by the microRNA miR172. Development 2009; 136:2873–2881. pmid:19666819
- 31. Choi SR, Teakle GR, Plaha P, Kim JH, Allender CJ, Beynon E, et al. The reference genetic linkage map for the multinational Brassica rapa genome sequencing project. Theoretical and Applied Genetics 2007; 115: 777–792. pmid:17646962
- 32. Wang F, Li L, Li H, Liu L, Zhang Y, Gao J, et al. Transcriptome analysis of rosette and folding leaves in Chinese cabbage using high-throughput RNA sequencing. Genomics 2012; 99: 299–307. pmid:22387604
- 33. Wang F, Li H, Zhang Y, Li J, Li L, Liu L, et al. MicroRNA expression analysis of rosette and folding leaves in Chinese cabbage using high-throughput Solexa sequencing. Gene 2013; 532: 222–229. pmid:24055726
- 34. Li R, Yu C, Li Y, Lam TW, Yiu SM, Kristiansen K, et al. SOAP2: An improved ultrafast tool for short read alignment. Bioinformatics 2009; 25: 1966–1967. pmid:19497933
- 35. Mortazavi A, Williams BA, McCue K, Schaeffer L, Word B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nature Methods 2008; 5: 621–628. pmid:18516045
- 36. Audic S, Claverie JM. The significance of digital gene expression profiles. Genome Reseach 1997; 7(10): 986–95.
- 37. Allen E, Xie Z, Gustafson AM, Carrington JC. microRNA-directed phasing during trans-acting siRNA biogenesis in plants. Cell 2005; 121: 207–221. pmid:15851028
- 38. Schwab R, Palatnik JF, Riester M, Schommer C, Schmid M, Weigel D. Specific effects of microRNAs on the plant transcriptome. Developmental Cell 2005; 8: 517–527. pmid:15809034
- 39. Jiang J, Shao Y, Du K, Ran L, Fang X, Wang Y. Use of digital gene expression to discriminate gene expression differences in early generations of resynthesized Brassica napus and its diploid progenitors. BMC Genomics 2013; 14: 72. pmid:23369045
- 40. Yan X, Dong C, Yu J, Liu W, Jiang C, Liu J, et al. Transcriptome profile analysis of young floral buds of fertile and sterile plants from the self-pollinated offspring of the hybrid between novel restorer line NR1 and Nsa CMS line in Brassica napus. BMC Genomics 2013; 14: 26. pmid:23324545
- 41. Song X, Li Y, Liu T, Duan W, Huang Z, Wang L, et al. Genes associated with agronomic traits in non-heading Chinese cabbage identified by expression profiling. BMC Plant Biology 2014; 14:71. pmid:24655567
- 42. Kim B, Yu HJ, Park SG, Shin JY, Oh M, Kim N, et al. Identification and profiling of novel microRNAs in the Brassica rapa genome based on small RNA deep sequencing. BMC Plant Biology 2012; 12: 218. pmid:23163954
- 43. Wang F, Li L, Liu L, Li H, Zhang Y, Yao Y, et al. High-throughput sequencing discovery of conserved and novel microRNAs in Chinese cabbage (Brassica rapa L. ssp. pekinensis). Molecular Genetics and Genomics 2012; 287: 555–563. pmid:22643909
- 44. Yu X, Wang H, Lu Y, de Ruiter M, Cariaso M, Prins M, et al. Identification of conserved and novel microRNAs that are responsive to heat stress in Brassica rapa. Journal of Experimental Botany 2012; 63: 1025–1038. pmid:22025521
- 45. Rajagopalan R, Vaucheret H, Trejo J, Bartel DP. A diverse and evolutionarily fluid set of microRNAs in Arabidopsis thaliana. Genes and Development 2006; 20: 3407–3425. pmid:17182867
- 46. Morin RD, Aksay G, Dolgosheina E, Ebhardt HA, Magrini V, Mardis ER, et al. Comparative analysis of the small RNA transcriptomes of Pinus contorta and Oryza sativa. Genome Research 2008; 18: 571–584. pmid:18323537
- 47. Szittya G, Moxon S, Santos DM, Jing R, Fevereiro MP, Moulton V, et al. High-throughput sequencing of Medicago truncatula short RNAs identifies eight new miRNA families. BMC Genomics 2008; 9: 593. pmid:19068109
- 48. Kim S, Hamada T, Otani M, Shimada T. Isolation and characterization of MADS box genes possibly related to root development in sweetpotato (Ipomoea batatas L. Lam.). Journal of Plant Biology 2005; 48: 387–393.
- 49. Mallory AC, Vaucheret H. Functions of microRNAs and related small RNAs in plants. Nature Genetics 2006; 38: S31–S36. pmid:16736022
- 50. Covarrubias AA, Reyes JL. Post-transcriptional gene regulation of salinity and drought responses by plant microRNAs. Plant Cell and Environment 2010; 33: 481–489.
- 51. Krol J, Loedige I, Filipowicz W. The widespread regulation of microRNA biogenesis, function and decay. Nature Reviews Genetics 2010; 11: 597–610. pmid:20661255
- 52. Bari A, Orazova S, Ivashchenko A. miR156- and miR171-binding sites in the protein-coding sequences of several plant genes. BioMed Research International 2013; Article ID 307145, 7 pages
- 53. Tsang J, Zhu J, Oudenaarden AV. MicroRNA-mediated feedback and feedforward loops are recurrent network motifs in mammals. Molecular Cell 2007; 26: 753–767. pmid:17560377
- 54. Liu H, Kohane IS. Tissue and process specific microRNA–mRNA co-expression in mammalian development and malignancy. PLoS One 2009; 4: e5436. pmid:19415117
- 55. Wang J, Hou X, Yang X. Identification of conserved microRNAs and their targets in Chinese cabbage (Brassica rapa subsp. pekinensis). Genome 2011; 54: 1029–1040. pmid:22111519