Differential miRNA expression in Rehmannia glutinosa plants subjected to continuous cropping

Background The productivity of the medicinally significant perennial herb Rehmannia glutinosa is severely affected after the first year of cropping. While there is some information available describing the physiological and environmental causes of this yield decline, there is as yet no data regarding the changes in gene expression which occur when the species is continuously cropped. Results Using a massively parallel (Solexa) DNA sequencing platform, it was possible to identify and quantify the abundance of a large number of R. glutinosa miRNAs. We contrasted the miRNA content of first year crop plants with that of second year crop ones, and were able to show that of 89 conserved (belonging to 25 families) and six novel miRNAs (six families), 29 of the former and three of the latter were differentially expressed. The three novel miRNAs were predicted to target seven genes, and the 29 conserved ones 308 genes. The potential targets of 32 of these differentially expressed miRNAs involved in the main transcription regulation, plant development and signal transduction. A functional analysis of the differentially expressed miRNAs suggested that several of the proposed targets could be directly or indirectly responsible for the development of the tuberous root. Conclusion We have compared differential miRNAs expression in the first year crop (FP) R. glutinosa plants and second year crop (SP) ones. The outcome identifies some potential leads for understanding the molecular basis of the processes underlying the difficulty of maintaining the productivity of continuously cropped R. glutinosa.


Background
Rehmannia glutinosa L. is a perennial herbaceous species belonging to the Scrophulariaceae family. Its economic importance results from the medicinal activity present in extracts of its tuberous roots [1]. Because of a lack of known undesirable side effects and its relatively low price, the species is extensively used in traditional Chinese clinical practice. Its prime production region is the Huai area of central China, but the climatic and edaphic conditions in Jiaozuo (Henan province) are also conducive for the cultivation of a high quality product. After one season of production, however, disease buildup (and other factors) forces the land to be cultivated with other crops for a period of 15-20 years [2]. Even in the absence of disease pressure, attempts to continuously crop over several seasons have failed to overcome the major decline in productivity, as the tubers are increasingly replaced by fibrous roots, which are unable to develop into tubers [3,4]. Much of the past research aimed at identifying the causative factors for this continuous cropping yield decline has been focused on the physiological activity and autotoxicity of the root exudates [5][6][7]. However, the molecular basis of the species' sensitivity to its own exudate remains unknown. miRNAs (short RNA molecules, on average~21 nucleotides in length) underlie a number of biological phenomena in the animal, plant and virus kingdoms [8], largely at the level of post-transcriptional gene regulation [9][10][11][12]. As their sequences are so highly conserved across the eukaryotes, they are believed to represent an evolutionarily ancient component of gene regulation. They operate via their complementarity to a stretch of mRNA sequence, and affect the level of gene expression by targeting the mRNA molecule for degradation.
Our hypothesis here was that miRNA activity may underlie some at least of the the problems associated with the continuous cropping of R. glutinosa. In order to gain a global picture of the miRNA content of R. glutinosa, we have therefore employed a high throughput parallel sequencing platform (Solexa sequencing) able to generate millions of short (18-30 nt) reads with a high level of accuracy. We have applied this technology to enable the comparative profiling of the miRNA content of plants in their first year of cropping (FP) with those in their second year (SP), with the intention of identifying miRNAs expressed differentially in FP and SP plants.

Results and Discussion
Sequencing and annotation of R. glutinosa miRNAs Solexa sequencing of the FP and SP miRNA libraries yielded, respectively, 17,619,697 and 18,028,647 unfiltered sequence reads. Of these 19.92% (unique 39.37%) were FP-specific, 23.82% (unique 48.17%) were SP-specific and 56.26% (unique 12.46%) were their common respectively. The average number of occurrences of the sequences common to both libraries was 10.5, while that of library-specific reads was not more than 1.2 (Table 1). After discarding the low quality reads, a total of 14,630,881 FP and 15,644,334 SP reads was retained. These sequences represented 6,748,998 and 7,894,661 unique clean reads in FP and SP, respectively. Their size distribution ( Figure 1) showed that~94% of the sequences in both libraries were of length 20-24 nt, with the modal length of 24 nt and the third peak at 21 nt, consistent with the observed length distribution of mature plant miRNAs [33,34].

Conserved miRNAs
Sequences homologous to non-coding sequences (rRNA, tRNA, small nuclear RNA and small nucleolar RNA) were identified from a search of the GenBank and the Rfam9.1 databases. This resulted in the allocation of 0.93% of the FP and 0.63% of the SP unique miRNAs to this category. When the remaining sequences were queried against known miRNA sequences, the outcome was the identification of 282,063 (unique 300) and 118,011 (unique 251) hits, accounting for, respectively, 1.93% and 0.75% of the FP and SP libraries. A BLASTn search of the genic miRNAs resulted in the identification of 89 sequences, belonging to 25 families. The extent of their conservation across the plant kingdom was shown by an alignment with the whole genome sequences of A. thaliana, soybean, rice, black poplar and grape (Table S1 in Additonal file 1). The most abundant sequences were miR156/157, miR172 and miR165/166; the former accounted for~47% of all conserved miRNAs in the FP library, while the most frequent single conserved sequence in the SP library was miR172 (~39%) ( Figure 2). In both libraries, miR159, miR394 and miR403 were moderately abundant. The five miRNAs miR161, miR397, miR398, miR408 and miR822 were absent from the SP library. It appeared therefore that the miRNA population present in FP plants differed to some extent from that present in SP plants.

Novel miRNAs
A distinguishing feature of miRNAs is the ability of their pre-miRNA sequences to adopt the canonical stem-loop hairpin structure. After removal of the conserved miR-NAs, 13,724,517 FP (6,685,869 unique sequences) and 15,043,261 (7,845,823 unique sequences) SP sequences were aligned with the A. thaliana genome sequence, producing 7,341 (3,158 unique sequences) FP and 7,468 (3,269 unique sequences) SP sequences ( Table 2) whose flanking region (in A. thaliana, at least) was amenable to secondary structure analysis. The application of a set of strict identification criteria for potential miRNA loci [35,36] resulted in the selection of 18 sequences (Additional file 2) across the two libraries which could be considered as likely novel miRNAs (Table S2 in Additional file 1). Except for miR5138, their frequency of occurrence was <40 (Table S2 in Additional file 1), reflecting an expression level considerably lower than that of the majority of the conserved miRNAs. When RT-PCR was applied to these 18 sequences, six were amplifiable from R. glutinosa cDNA template ( Figure 3).

Differentially expressed miRNAs
Evidence for differential expression in FP and SP plants was sought by comparing the frequency of occurrence of the 89 conserved and six novel miRNAs, based on a Poisson distribution approach [37]. The 29 conserved (11 miRNA families) and three novel miRNAs showing the greatest degree of differential expression are listed in Table 3. Of these 32 sequences, 12 showed a greater than two fold difference in expression level between FP and SP plants. Seven of them were more strongly expressed in SP than in FP plants. The expression levels of the most differentially expressed (17 conserved and 3 novel) miRNAs were reanalysed using qRT-PCR. This confirmed that 14 of the former and two of the latter sequences were indeed differentially expressed in FP and SP plants (Figure 4), showing that frequency of occurrence in Solexa runs produces a reasonably accurate prediction for expression level. Expression levels of 4 miRNAs (miR157a, miR167a, miR160a and miR5138) in roots were measured in different times ( Figure 5). miR157a and miR167a were highly expressed in FP, Figure 1 Size distribution of R. glutinosa small RNAs.  while miR160a and miR5138 were quite opposite, with strongly expressing in SP.
Target prediction for the three differentially expressed novel miRNAs The target of most plant miRNAs possesses a single perfect or near perfect complementary site in the coding region [13,15]. Assuming this to be generally the case, the A. thaliana gene space was searched for complementarity with the sequences of the three differentially expressed novel miRNAs. Using a set of rules for predicting novel miRNA potential target genes [14,38], this exercise predicted seven potential targets, with miR5138 and miR5140 both targeting more than one gene (Table S3 in Additional file 1). The targets encoded the following gene products: ICU2 (INCURVATA2), a DNA-directed DNA polymerase, a magnesium transporter CorA-like family protein, an ATP synthase (α chain), a TIR-NBS-LRR protein, a ZIGA4 (ARF GAP-like zinc finger-containing protein ZiGA4) and a DC1 domain-containing protein.

Function of the potential targets of differentially expressed miRNAs
An indication of the genes responsible for the continuous cropping syndrome was sought by an inspection of the 308 potential targets of the 29 differentially expressed miRNAs (Additional file 3) in addition to the seven targets of the novel miRNAs (Table S3 in Additional file 1). Gene ontology categories were assigned to all 315 putative targets according to their cellular component, their molecular function and the biological process(es) they are involved in A. thaliana ( Figure 6). With respect to molecular function, the targets fell largely into nine categories, with the three most overrepresented being nucleic acid binding, metal ion binding and transcription factor activity. Twelve biological processes were identified, with the three most frequent being the regulation of transcription, plant development and signal transduction. The potential targets for the 32 differentially expressed miRNAs mainly involved transcription, plant development and signal transduction. Several of these targets may be directly or indirectly involved in the development of tuberous vs fibrous roots (Figure 7, 8). For example, miR156/157 targets an SPL transcription factor, which when over-expressed in A. thaliana, produces an early flowering phenotype. The over-expression of miR156/157 itself delays flowering [39][40][41]. Thus it is possible that in R. glutinosa, a higher level of expression of miR156/157 (as occurred in FP plants) could prolong root growth and development. The miR160 target is the auxin response factor ARF17, while those of miR167 are ARF6 and ARF8. ARF17 is a negative regulator, while ARF6 and ARF8 are positive regulators of adventitious rooting. These three ARFs share overlapping expression domains, interact genetically and regulate one another's expression at both the transcriptional and post-transcriptional level [42]. Since SP plants express more miR160 and less miR167 than FP plants, it is possible that the balance of ARF protein present is altered by continuous cropping, and hence there is an effect on tuberous root expansion. The target of miR5138 is the gene ICU2, which is negative regulator of floral homeotic genes in A. thaliana. Its overexpression delays flowering, while its knock-out hastens  it [43]. Since this miRNA is more highly expressed in SP than in FP plants, there may be a differential expression of ICU2 and hence an effect on flowering time, with a knock-on effect on tuberous root expansion. These predicted target genes were cloned in R. glutinosa (Table 4 and Additional file 4) and registered as ESTs in NCBI.
Overall, there was a suggestion that the expression of a number of miRNA families may be correlated with the continuous cropping syndrome in R. glutinosa. Whether these miRNAs actually regulate key genes responsible for the syndrome will require experimental demonstration. The identification of these miRNAs has nevertheless succeeded in providing leads for determining the molecular genetic basis of the continuous cropping syndrome in R. glutinosa.

Conclusions
Here we have described the application of a combination of approaches to identify a set of 89 conserved (belonging to 25 families) and six novel R. glutinosa miRNAs, which are differentially, expressed in first and second year crops. We believe that this information could provide initial candidates for the genes responsible for tuberous root expansion, and in particular for the syndrome of continuous cropping yield decline in this medicinally important species.

Methods
Plant material and RNA isolation R. glutinosa cultivar "Wen 85-5" was collected from the Wen Agricultural Institute, Jiaozuo City, Henan Province, China. The first year crop (FP) was grown from   For the measure of differently expressed miRNAs in various development stages of R. glutinosa, FP and SP plants (cultivar "Wen 85-5") were grown in the isolated plots from April 22 to October 22, 2010. Roots of R. glutinosa were collected every month and total RNAs were extracted with TriZOL reagent.

miRNA library construction and sequencing
The small RNAs were ligated sequentially to 5' and 3' RNA/DNA chimeric oligonucleotide adaptors (Illumina), and the resulting ligation products were gel purified by 10% denaturing PAGE, and reverse transcribed. The cDNAs obtained in this way were sequenced on a Genome Analyzer IIx System by Beijing Genomics Institute (BGI) (Shenzhen, China).  a_thaliana/ath1/SEQUENCES/ using SOAP (soap.genomics.org.cn) software [44]. Candidate pre-miRNAs were identified by folding the flanking genome sequence of distinct miRNAs using MIREAP (mireap.sourceforge. net), followed by a prediction of secondary structure by mFold v3.1 [45]. The criteria chosen for stem-loop hairpins were as suggested elsewhere [35,36].
Validation of differential miRNA expression based on qRT-PCR qRT-PCR was performed using an All-in-One™ miRNA Q-PCR detection kit (GeneCopoeia, Inc.) on a BIO-RAD iQ5 real-time PCR detection system (Bio-Rad laboratories, Inc.). Each 20 μl Q-PCR comprised 0.5 μl cDNA, 2 μl 2 μM miRNA forward primer (sequence given in Table S5 in Additional file 1), 2 μl 2 μM reverse primer (Universal Adaptor PCR Primer), 10 μl 2× All-in-One™ miRNA Q-PCR buffer and 5.5 μl nuclease-free water. The reactions were incubated at 95°C for 10 min, then were cycled 36 times through 95°C/10 s, 55°C/20 s and 72°C/10 s. After the reactions had been completed, the threshold was manually set and the threshold cycle (CT) was automatically recorded. All reactions were replicated twice per biological sample. A 4 μl aliquot of each reaction product was subjected to 3% agarose electrophoresis. The relative expression level of the miRNAs was calculated using the 2 -ΔΔCT method [46], and the data were normalized on the basis of 18 s rRNA CT values.