Next Article in Journal
Current Research on Non-Coding Ribonucleic Acid (RNA)
Next Article in Special Issue
Tissue-Specific Transcriptome Analysis Reveals Multiple Responses to Salt Stress in Populus euphratica Seedlings
Previous Article in Journal
Delivery Mode and the Transition of Pioneering Gut-Microbiota Structure, Composition and Predicted Metabolic Function
Previous Article in Special Issue
Comparative Analysis of Soybean Root Proteome Reveals Molecular Basis of Differential Carboxylate Efflux under Low Phosphorus Stress
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Comparative Analysis of Cotton Small RNAs and Their Target Genes in Response to Salt Stress

1
State Key Laboratory of Cotton Biology, Institute of Cotton Research of Chinese Academy of Agricultural Sciences, Anyang 455000, China
2
State Key Laboratory of Crop Biology, College of Agronomy, Shandong Agricultural University, Tai’an 271018, China
*
Authors to whom correspondence should be addressed.
Submission received: 8 October 2017 / Revised: 23 November 2017 / Accepted: 24 November 2017 / Published: 5 December 2017
(This article belongs to the Special Issue Genetic Regulation of Abiotic Stress Responses)

Abstract

:
Small RNAs play an important role in regulating plant responses to abiotic stress. Depending on the method of salt application, whether sudden or gradual, plants may experience either salt shock or salt stress, respectively. In this study, small RNA expression in response to salt shock and long-term salt stress in parallel experiments was described. Cotton small RNA libraries were constructed and sequenced under normal conditions, as well as sudden and gradual salt application. A total of 225 cotton microRNAs (miRNAs) were identified and of these 24 were novel miRNAs. There were 88 and 75 miRNAs with differential expression under the salt shock and long-term salt stress, respectively. Thirty one transcripts were found to be targets of 20 miRNA families. Eight targets showed a negative correlation in expression with their corresponding miRNAs. We also identified two TAS3s with two near-identical 21-nt trans-acting small interfering RNA (tasiRNA)-Auxin Response Factors (ARFs) that coaligned with the phases D7(+) and D8(+) in three Gossypium species. The miR390/tasiRNA-ARFs/ARF4 pathway was identified and showed altered expression under salt stress. The identification of these small RNAs as well as elucidating their functional significance broadens our understanding of post-transcriptional gene regulation in response to salt stress.

1. Introduction

Salt stress is the exposure of plants to high salinity, the main component of which is Na+ ions. It affects various aspects of plant growth and development, such as inhibition of enzymatic activities and reductions in photosynthetic rates [1]. Depending on the method of NaCl application, whether in a single step or gradual, plants experience either salt shock or salt stress mainly with ionic stress component, respectively [2]. Salt shock and salt stress mainly with ionic stress component are two distinct phenomena. Salt shock immediately induces osmotic shock when roots come into contact with solutions containing unfavourably high concentrations of salts in hydroponic systems or in soil [3]. However, the concentration of Na+ reaches a certain toxic level in cell protoplasts of shoots under long-term salt stress [4]. Osmotic stress results in strong, rapid changes in the expression of genes with antioxidant system function, and a few changes in ion-responsive genes [5]. Many genes with altered expression under long-term salt stress maybe involved in responses to toxic cellular concentrations of Na+ ions. There are very few studies in which the effects of salt shock and long-term salt stress are described in parallel experiments [6,7]. A more detailed analysis of the process should broaden our understanding of the molecular events of the mechanisms involved in salt stress response.
The wave of gene activation is a consequence of stepwise transcriptional, post-transcriptional and post-translational regulation [8]. MicroRNAs (miRNAs), an abundant class of short non-coding (21–23 nt) RNAs, are generated from single-strand hairpin RNA precursors [9,10]. They regulate gene expression by binding to targeted mRNA sequences in a perfect or near-perfect complementary site. After that, they lead to targeted sequence cleavage or, in at least a few cases, translational repression in plants [11]. Many processes, such as leaf development, auxin signaling, phase transition, flowering and genome maintenance are regulated in similar ways by different miRNAs [12]. Another type of small RNAs in plants is small interfering RNAs (siRNAs), a class of double-strand RNAmolecules that also act in RNA interference (RNAi)-related pathways [11]. This encompasses heterochromatin-associated siRNAs and trans-acting siRNAs (tasiRNAs), repeat-associated siRNAs (rasiRNAs) and natural antisense transcript siRNAs (nat-siRNAs) [13].
Recently findings have firmly established a role for miRNAs in plant stress responses [14,15]. For example, Arabidopsis miR393 is strongly up-regulated by cold, dehydration, NaCl and abscisic acid (ABA) treatments [16]. It targets the F-box auxin receptors Transport Inhibitor Response Protein 1 (TIR1) and Auxin Signaling F-box Protein 1 (AFB1) [17]. Over-expression of osa-MIR393 leads to more tillers, early flowering and less tolerance to salt and drought [18]. Xie et al. [19] investigated the transcriptional profiles of cotton miRNAs during exposure to drought and salinity stress. miR413 and miR440 were found to be expressed specifically in salt treatments. Some of miRNAs were also found to be expressed differentially in terms of dose dependence and tissue dependence under stress conditions. In cotton roots exposed to salinity, miR399 shows down-regulated expression levels with the increase of salt concentration [20]. The expression levels of miR156, miR157 and miR172 were up-regulated at 0.25% NaCl and then down-regulated at higher salinity concentrations [20]. Many discoveries that stress can regulate miRNA levels, coupled with the identification of stress-associated genes as miRNA targets provided clues about the role of miRNAs in stress responses [21]. For example, PHO2 (phosphate responsive mutant 2), encoding E2 conjugase, has been verified as a target of miR399 in model plant species [22]. In addition, constitutive miRNAs in maize or cotton were shown to be capable of genotype-specific expression in response to salt stress [23,24]. Therefore, a more detailed analysis of the kinetics of miRNAs and their target genes could be helpful in obtaining better insights into the mechanisms and roles of miRNAs in salt-regulatory networks. Cotton (Gossypium hirsutum L.) is considered to be relatively tolerant to salinity [25]. Its growth and productivity are also adversely affected by salt stress. In the past decades, much progress has been made in unraveling the complex stress response mechanisms, and considerable information has become available on the identification of stress responsive genes. However, little is known about the roles of small RNAs in response to different salt stress phases. In this study, the small RNAome, degradome and transcriptome were studied to investigate miRNA expression profiles in response to the effect of sudden and gradual salt application. These data provide new and valuable insights into the expression profiles of miRNAs, other small RNAs and target genes, and their relationships under salt stress.

2. Materials and Methods

2.1. Plant Materials and Growth Conditions

The salt-tolerance cotton cultivar ‘SN-011’ was used in this study. Seeds were sterilized and germinated in vermiculite. Growth conditions were 30/22 °C day/night temperature, 55–70% relative humidity and a 14/10 h light/dark cycle under 450 μmol m−2·s−1 light intensity. At the two-leaf stage, healthy seedlings were placed in pots containing aerated Hoagland nutrient solution. Plants were cultured under non-saline conditions for 10 d to ensure full establishment before starting the salinity treatments. The pH was maintained close to 6.9 by adding H2SO4 or KOH as required. Seedlings were randomly divided into two groups that were treated with salt stress. The remaining seedlings served as control. In all experiments 15 plants were used for each treatment. For the sudden salt application, seedlings were suddenly transferred from normal growth solution into a 150 mM NaCl solution. Leaves were harvested after exposing the seedlings to salt stress for 2, 4 and 8 h. For the gradual salt application, seedlings were exposed to salinity by adding NaCl to the growth medium in 50-mM increments every 12 h, until a final concentration of 150 mM was reached. After exposing the seedlings to salt stress for 1, 3 and 5 d, leaves were harvested directly into liquid nitrogen and stored at −80 °C for subsequent use.

2.2. Measurement of Na+ Contents and Proline Level

Leaves were dried for 72 h at 70 °C, and the dry weight measured. Dried leaves were then digested with 0.1 M HNO3, and Na+ contents were determined by inductively coupled plasma-optical emission spectrometry (ICP-OES) (Optima 2100 DV; Perkin-Elmer, Wellesley, MA, USA) according to the manufacturer’s instructions. Proline level was measured using commercial assay kits (Jiangsu Suzhou Keming Bioengineering Institute, Suzhou, Jiangsu, China) according to the manufacturer’s instructions. Three independent biological replicates were performed.

2.3. Total RNA Isolation andConstruction of Small RNA and Degradome Libraries

Total RNA was isolated using a modified CTAB method with isopropanol instead of lithium chloride for RNA precipitation and purified using RNeasy columns (Qiagen, Hilden, Germany) [26]. RNA quantity and quality were checked by Agilent 2100 Bioanalyzer. The small RNA digitalization analysis based on HiSeq high-throughput sequencing takes the sequencing by synthesis (SBS), which is useful due to its small requirement of sample quantity, high throughput, high accuracy with a simply operated automatic platform. Small RNA libraries were constructed as previously described [27,28]. In brief, total RNA was separated by 15% denaturing polyacrylamide gel electrophoresis to recover the population of small RNAs. The extracted RNAs were then dephosphorylated using alkaline phosphatase and ligated to 5′- and 3′-RNA adaptors by T4 RNA ligase. The adaptor-ligated small RNAs were subsequently transcribed to single-strand cDNA using SuperScript II Reverse Transcriptase. After PCR amplification, PCR products were purified and subjected to SBS sequencing with an Illumina Solexa genome analyzer II at the BGI (Beijing Genomics Institute, Shenzhen, China).
Degradome sequencing provides a comprehensive means of analyzing patterns of RNA degradation. The libraries were constructed as previously described [29]. In brief, poly(A+) RNA molecules were isolated from 150 μg of total RNA using the Oligotex mRNA kit. Subsequently, the 5′-phosphate of the poly(A+) RNA was ligated to a 5′-RNA oligonucleotide adaptor containing a MmeI recognition site by T4 RNA ligase. The products of a reverse transcription reaction were amplified by five PCR cycles and were then digested with MmeI and ligated to a 3′ double-DNA adaptor. The ligation products were purified and amplified by 20 PCR cycles and gel-purified for SBS sequencing.

2.4. Transcriptome Sequencing andde novoAssembly

To maximize the transcript representation in a broad range of biological processes, the total RNAs harvested at 0 h, 4 h and 5 d after salt treatment were pooled in an equal fraction ratio. Poly(A+) mRNAs were isolated using beads with Oligo(dT). Fragmentation buffer was then added to interrupt mRNA into short fragments. The first-strand cDNA was synthesized using random hexamer-primer. The second-strand cDNA was synthesized using buffer, dNTPs, RNaseH and DNA polymerase I, respectively. Short fragments were purified and resolved with Elution buffer for end reparation and tailing A. Sequencing adapters were linked with short fragments. After that, suitable fragments were selected by agarose gel electrophoresis, and were used for the PCR amplification as templates. The library was sequenced using Illumina HiSeq™ 2000 (Illumina, CA, USA). Dirty raw reads which contained adapters, unknown or low quality bases were discarded. Transcriptome de novo assembly was carried out with Trinity, a short reads assembling program [30]. Briefly, Trinity firstly combined reads with certain length of overlap to form contigs. Then the contigs with paired-end reads were connected to form unigene sequences that could not be extended on either end. For function annotation, unigene sequences are aligned by blastx to protein databases of G. hirsutum (e-value < 0.00001), and aligned by blastn to CDS databases of G. hirsutum (e-value < 0.00001), retrieving proteins with the highest sequence similarity with the given Unigenes along with their protein functional annotations.

2.5. Analysis of Small RNA and Degradome Sequencing Data

The very basic figure from small RNA and degradome sequencing is converted into sequence data by the base-calling step. As for the small RNA library, raw sequencing reads were processed into clean full-length reads by the BGI small RNA pipeline [28]. After the raw sequence reads were obtained, adapter–adapter sequences, low quality tags, contaminants and insert sequences beyond 15–30 nt were filtered. The remaining sequences were aligned to Rfam and GenBank databases to annotate rRNA, tRNA, snRNA, snoRNA and those containing polyA tails [31]. Next, sequence tags were searched against the miRNA database, miRBase (Release 20.0, June 2013), in order to identify similar miRNAs. All plant miRNAs were downloaded and were utilized as dataset to perform alignment searches. The reads that were shorter than 24 nt, and had more than one copy in a read library, either shorter/longer or contained up to two mismatches, and no match to known non-coding RNAs, were tested whether they are miRNAs [32]. The small RNA sequences which mapped on unigenes from transcriptome sequencing were used to predicted novel potential miRNAs. The novel miRNAs were identified using the MIREAP program developed by the BGI [33]. For this procedure, the occurrence of both a miRNA and a miRNA* in the deep sequencing data was the basis to predict novel miRNAs. The approximately 23 nt potential miRNA sequence is not located on the terminal loop of the hairpin structure. ThemiRNA has less than six mismatches with the opposite miRNA sequence (miRNA*) on the other arm [34]. The filtered pre-miRNA sequences were folded again using MFOLD software and selected manually [35,36]. MiRNAs expression abundance was analyzed by counting the number of transcripts per million (TPM) clean reads in libraries. We first normalized the read density measurement and then used p-value < 0.01 and the absolute value of |log2Ratio|>1 as a threshold to judge the statistical significance of miRNA expression. From the three data sets, plenty of miRNAs were found to be differentially expressed between libraries.
Degradome sequencing takes the SE50 sequencing strategy and produces 49-nt raw reads. The 3′-adaptor was trimmed before bioinformatics analysis to get real degradome fragments with length of 20–21 nt. Data filtering was carried out to obtain high quality reads as clean reads. The unique tags and total tags were first analyzed. Then, the degradome tags were annotated with rRNA, scRNA, snoRNA, snRNA and tRNA from GenBank and Rfam, and selected out from unannotated tags. The degradome sequences with a single base representing over 70% of the bases were regarded as polyN and tossed. To map every unique degradome tag to only one of these annotations, we followed the priority rule: Rfam > GenBank > PolyN. Unannotated tags were mapped to reference genes (cDNA). The tags mapped to cDNA_sense were used to predict cleavage sites. After a series of annotation of the degradome sequence, the degradation mRNA obtained go through comparison with miRNA to find mRNA-miRNA pairing. Sequences with less than 4 nt mismatches compared with the query miRNA sequences were selected. The G:U pair was counted as 0.5 mismatch. Cleavages induced by different miRNA sequences were regarded as different cleavage events. The reference genes which were predicted to degradome tags were annotated. The T-Plot figure was constructed to show the tag distributions of one gene.

2.6. qRT-PCR Analysis

For cDNA synthesis,1 μg of total RNA was used. First-strand cDNA was synthesized using a PrimeScript RT reagent kit with gDNA eraser (Takara, Dalian, China). Premix Ex TaqTM II(TaKaRa)was used for real-time quantitative PCR (qRT-PCR). Gene-specific primers are shown in Table S1. qRT-PCR was undertaken with a 7500 Fast Real-Time PCR system (Applied Biosystems Inc., Carlsbad, CA, USA). The 25-μL reaction solutions contained 12.5 μL 2XSYBR Premix Ex TaqTM II, 1 μL forward primer, 1 μL reverse primer, 2 μL DNA template and 8.5 μL H2O. For the mature miRNAs, RT-PCR employed a stem-loop primer to detect miRNAs. The primers were based on those described by Chen et al. and Varkonyi–Gasic et al. [37,38]. Their 3′-ends, which were complementary to the six nucleotides at the miRNA 3′-end, were combined with the 44-nt sequence to form a stem-loop structure.The stem-loop reverse transcription reactions were performed using M-MLV Reverse Transcriptase (Invitrogen) according to the supplier’s manual. U6 was used as a reference for small RNA expression validation [39]. The expression levels of target transcripts were normalized using the housekeeping gene UBQ7 as the reference control. Three biological replicates were used for qRT-PCR validation. The 2ΔΔT method was used to calculate relative expression levels of salt-responsive miRNAs and target genes.

2.7. Statistical Analysis

The statistical significance of Na+ contents, proline level, expression profiles of miRNAs and genes was compared by one-way analysis of variance (ANOVA), followed by a Duncan’s multiple range test. All data were analyzed by SPSS10.0 (SPSS Inc., Chicago, IL, USA) software. p values of less than 0.05 were considered to be statistically significant.

3. Results

3.1. Deep-Sequencing of sRNAs under Salt Stress

Cotton seedlings were randomly divided into two groups to experience two types of salt treatment. One group was suddenly transferred from normal growth solution into 150 mM NaCl solution, experienced salt shock for 2 h, 4 h and 8 h. The remaining seedlings were employed gradual applications of salt, by adding NaCl to the growth medium in 50-mM increments every 12 h, until a final concentration of 150 mM was reached. They were subjected to long-term salt stress for 1 d, 3 d and 5 d. We measured Na+ and proline content of seedling leaves under salt stress treatment. The Na+ contents in cotton leaves showed an overall increasing trend with the increase of processing time (Figure 1). The 5-d samples showed 1.70-fold higher Na+ content than the control samples (CK). No significant difference in Na+ content was detected between 8-h samples with sudden NaCl application and 1-d samples with gradual NaCl application. In the 2-h, 4-h and 8-h samples with sudden NaCl application, the proline level was close to those observed in control samples. The content of proline increased after 1 d salt stress treatment under gradual salt application conditions. Based on these results, we selected two time-points (4 h and 5 d) to indicate the effect of sudden and gradual salt application, respectively, for further analysis of salt stress in cotton.
Three separate cDNA libraries of small RNAs were constructed using total RNA obtained from CK and the 4-h and 5-d salt-treatment samples. The small RNA digitalization analysis takes the SBS-sequencing by synthesis. Sequencing raw data were available at Gene Expression Omnibus (GEO: GSE69702). A total of 57,370,554 sequence reads were obtained from the three libraries. After removing the 3p- and 5p-adapter sequences and filtering out low quality ‘n’ sequences, there were 17,026,070, 15,731,477 and 22,590,391 clean reads remaining in the CK, 4-h and 5-d samples, respectively (Table 1). Further analysis identified a total of 12,837,490 unique sequence tags in the three libraries. Approximately 78% (10,023,245) of the unique tags corresponded to singletons. In all three libraries, the majority of the small RNA sequences were 21 nt, followed by 24 nt. The 4-h sample showed higher 21-nt small RNA abundance (48.86%) than CK (40.63%). The frequencies of 21- and 24-nt small RNAs were similar for CK and 5-d samples.

3.2. Small RNA Annotation and miRNA Identification

In order toinvestigate the repertoire of miRNAs in cotton, the reads in the small RNA libraries were first mapped to protein-coding and structural RNA-coding RNAs (rRNAs, tRNAs, snRNAs and snoRNAs) in the Rfam and GenBank databases (Figure S1). Thereafter, the remaining sequences were aligned with all known miRNAs in the central miRNA Registry Database, miRBase database. As a result, a total of 201 conserved miRNAs were identified. They were classified into 88 families on the basis of sequence similarity. Among these conserved miRNA families, miR156 and miR166 contained 11 and 12 members, respectively. miR160, miR171, miR172 and miR396 contained nine members; miR167 and miR399 contained seven members; miR159, miR164 and miR390 contained five members; and miR168, miR169, miR319, miR393 and miR482 contained four members. These miRNA families were evolutionarily conserved and had homologs in more than 14 other plant species. However, 59 miRNA families, such as miR391, miR403 and miR822, comprised only one member and showed a low level of conservation.
miRNA genes mostly originate from independent transcriptional units. Transcriptome sequencing was performed to identify primary transcripts of novel miRNAs. The transcriptome raw data have been deposited in the National Center for Biotechnology Information (NCBI) Sequence Read Archive under SRX2865307. We obtained 52,170,788 reads and assembled them into 159,520 contigs. A total of 81,114 unigenes were identified (Table 1), among which 55,689 were annotated in the NCBI non-redundant protein database and 34,147 in the Swiss-Prot database. The remaining non-encoding sequences were used to predict novel miRNAs based on their secondary structure. As a result, 85 candidate novel miRNAs were detected. For 24 of them, the complementary miRNA* species were cloned in libraries, providing an indication of precise excision from the stem-loop precursor. The lengths of these newly identified miRNAs varied within 21–23 nt (Table 2). The predicted hairpin structures of the pre-miRNAs were 76–262 nt, with the majority (69.23%) 90–150 nt (Figure S2). Each miRNA-precursor had high minimal folding free energy index (MFEI) values (0.64–1.74 with an average of about 1.02). The values were significantly higher than that for tRNAs (0.64), rRNAs (0.59) and mRNAs (0.62–0.66) [40].

3.3. miRNA Expression Patterns under Different Salt Stress Treatment

miRNAs in small RNA libraries had a very broad range of expression. A few conserved miRNA families, such as miR156, miR157, miR166, miR167 and miR168, showed extraordinarily high expression levels in three libraries. miR157 was the most abundant, with 4,217,486, 5,086,032 and 4,576,202 reads in CK, 4-h and 5-d libraries, respectively. A few less-conserved miRNA families, such as miR829, miR858, miR2118 and miR7490, showed relatively lower expression levels, represented by <10 reads. The results were in agreement with previous studies reporting that less-conserved miRNAs were represented at relatively lower levels than conserved miRNAs. Moreover, different members in the same miRNA family also had significantly different expression levels. For instance, the abundance of miR156 members varied within 7–580,406 reads. The varied frequency of miRNAs suggested their distinct regulatory roles in response to salt stress.
miRNA expression abundance in datasets was analyzed by counting the number of transcripts per million (TPM) clean reads in libraries. We first made a comparative analysis of miRNA expression between CK and 4-h libraries (Figure 2). A total of 78 conserved and ten novel miRNAs were differentially expressed between the two samples. Among them, 37 miRNAs were down-regulated and 51 up-regulated. miR780 and miR7498 had the highest down- and up-regulated fold changes, respectively. Next, we compared the normalized expression levels of miRNAs from CK and 5-d libraries, and found 75 miRNAs differentially expressed between the two samples. Among them, 45 conserved and three novel miRNAs, including miR1336, miR1355 and miR1356, were up-regulated in the 5-d compared with CK sample. There were 23 conserved and four novel miRNAs, including miR1340, miR1344, miR1350 and miR1357, down-regulated in the 5-d sample. The most abundant miRNA with differential expression was miR1344, which exhibited a high level (read count 149,649) in CK sample but was less expressed in 5-d library (read count 92,410). Under the two types of salt treatment, a total of 141 miRNAs showed differential expression between all salt treatment and CK samples. Among them, 20 and nine miRNAs showed consistent up- or down-regulation, respectively. Six miRNAs, including miR472, miR771, miR775, miR822, miR823 and miR839, were specifically expressed in the 4-h sample. Three miRNAs, including miR159d, miR166l and miR169d, were specifically expressed in the 5-d sample. These miRNAs had less accumulation of about 2–10 counts. To confirm the expression patterns of the differentially expressed miRNAs, ten miRNAs were selected for a stem-loop RT-PCR experiment (Figure 3). Stem-loop RT-PCR has become the method of choice due to its sensitivity and robust quantification of miRNA levels [41]. Dissociation curve analysis revealed single peaks for the selected primer sets (Figure S3). Real-time quantitative PCR results indicated that seven miRNA showed similar expression patterns to the sequencing results. There was no significant difference between the expression levels of miR482a and miR1355 in the CK and 4-h samples. The remaining three miRNAs tested showed no detectable expression, an observation possibly related to their low accumulation or lower primer set efficiency.

3.4. Retrotransposons and Their Endogenous siRNAs

There are 5578 Gossypium retrotransposon elements deposited in NCBI. In this study, 36 retrotransposons were found in transcriptome sequencing data. Ten retrotransposon sequences were identical to or had some mismatches with transcripts. For example, three genomic survey sequences, DX403101, DX402782 and DX404336 were mapped on unigene CL2334.Contig1 (Figure 4). They were defined as Gorge1 gypsy-like, copia-like and unknown type retrotransposons by Hawkins [42]. Our research showed that some small RNAs had a perfect match with retrotransposon sequences. siRNA is a 21–24 nt long double-strand RNA, each strand of which is 2 nt longer than the other at the 3′-end. According to this structural feature, we aligned small RNAs to each other to find siRNAs meeting this criterion. A total of 184,498 tags were found to form RNA duplexes. These tags were potential siRNA candidates. The majority of them (58.52%) were 24 bp, indicating the function of DCL3 in this processing. Furthermore, we found 3933 siRNAs mapped on the expressed retrotransposons. Interestingly, 2296 siRNAs were mapped on DX404336 and DX403101, and 1259 on DX402782. Their distribution on sequences was uneven (Figure 4). Using qRT-PCR, we examined the expression of the two retrotransposons regions—both were induced in 4-h sample and repressed in 5-d sample (Figure 5).

3.5. Transcriptome-Wide Identification of miRNA Targets in Cotton

To gain insight into the functions of salt-regulated miRNAs, target genes were identified through a degradome sequencing approach. Three libraries, prepared from the CK, 4-h and 5-d samples were constructed. A total of 73,988,644 reads represented by 3,254,054 unique reads from the 5′-ends of uncapped and poly-adenylated RNAs were obtained (GEO: GSE69820). These sequences were compared with transcriptome sequencing data to identify the sliced targets for the known and novel miRNAs.
Based on degradome sequencing, a total of 31 target genes were identified for 20 cotton miRNA families (Table 3). The abundance of transcripts was plotted for each transcript. It is well-established that conserved miRNAs target conserved homologous genes in diverse plant species (Figure 6). In this study, 15 target genes for 10 highly or moderately conserved miRNAs were identified. These included five auxin response factors (ARF) targeted by miR160; two No Apical Meristem domain transcription factor (NAC) targeted by miR164; one HD-Zip protein targeted by miR165; one ARF6 gene targeted by miR167; one NF-YA3 gene targeted by miR169; one AP2 targeted by miR172; one MYB targeted by miR319; one growth-regulating factor (GRF) gene targeted by miR396; and one MYB targeted by miR828. MYB proteins are characterized by their MYB domain. Interestingly, we found that two genes of MYB transcription factor family proteins emerged as the targets of miR319 and miR828. The divergence of miRNA target sites suggested that subgroups of the MYB gene family might be regulated by different miRNAs. To test the biological function of salt-regulated miRNAs, a negative-correlation expression test was undertaken for their target mRNAs using qRT-PCR. Among the 10 targets tested, eight showed a negative correlation of expression with their corresponding miRNAs (Figure 7). Negative correlations between miRNAs and target mRNAs did not occur in all cases. There was no significant difference between the expression levels of ARF 16 (target of miR160) in the 4-h and 5-d samples. The expression of NAD9, a target of miR7505, was significantly down-regulated in 4-hr samples and was up-regulated in 5-d samples.

3.6. Cotton miR390/tasiRNA-ARF/ARF4 Pathway Analysis

Using a homology search approach, four tasiRNA-ARFs were identified in cotton (Table 4). Based on Arabidopsis tasiRNA-ARF nomenclature, they were named TAS3a D7(+), TAS3a D8(+), TAS3b D7(+), and TAS3b D8(+). They were identical to or had one nucleotide mismatch with their corresponding homologs in Arabidopsis. Transcriptome analysis was further used for the identification of TAS3 loci. TAS3a and TAS3b were mapped on the Unigene17307 and Unigene12729 sequences, respectively.
To investigate the molecular evolution of the two TAS3 loci in Gossypium species, we also identified TAS3 loci in G. arboretum and G. raimondii, the two diploids whose interspecific hybridization led to the tetraploid G. hirsutum. Each diploid Gossypium species had two TAS3 loci, with the same sequences as in G. hirsutum, indicating the highly conserved nature of TAS3 in Gossypium. TAS3 precursors had two recognition sites for miR390. The space between the two sites in Unigene17307 was 261 nt (12 phases); and was 174 nt (eight phases) in Unigene12729. However, our degradome sequencing results showed that only one target site was cleaved in both TAS3 precursor RNAs (Figure 8A). When the expression levels of miR390, TAS3a and TAS3b were monitored in cotton exposed to salt stress (Figure 9A), the miR390 levels reduced as the duration of exposure was extended. In contrast, the level of TAS3b increased. TAS3a level was firstly reduced, and then increased. Their negative correlation of expression supports the cleaving function of miR390 for TAS3, especially TAS3b.
Transcriptome sequencing data was used to predict potential target genes of tasiRNA-ARFs, resulting in detection of six transcripts with one or two recognition sites. The potential target mRNAs showed near perfect complementarity to the 21-nt tasiRNA-ARFs, with <2 nt mismatches. Subsequently, the CDS database of G. arboretum, G. raimondii and G. hirsutum was used to examine the functional similarity of predicted target mRNAs. CotAD_43230, CotAD_76410 and CotAD_58764 were annotated as ARF2 (Table 5). CotAD_00650, CotAD_64060 and CotAD_72124 were annotated as ARF3. CotAD_13694, CotAD_54309 and CotAD_12422 were annotated as ARF4. ARF3 and ARF4 contained two recognition sites. ARF2 contained one single recognition site. All of the tasiRNA-ARF recognition sites were located in coding regions (Figure S4). Of these ARF genes, cleavage sites were identified in ARF4 in our degradome study. We examined the expression levels of ARF4 and tasiRNA-ARF [TAS3b D7(+)] exposed to salt stress. As the duration of exposure increased, the tasiRNA-ARF level reduced and ARF4 level increased (Figure 9B).

4. Discussion

In this study, we analyzed Na+ content and proline levels, and monitored the molecular responses of cotton plants exposed to sudden and gradual salt application. Na+-specific damage is associated with the accumulation of Na+ in the plant cell cytoplasm. The low increase rates of Na+ concentration at the beginning of salt shock suggested that ion stress was at a low level. The Na+ contents reached a high level at five days under long-term salt stress. It suggests that external salt concentration caused an ion imbalance under salt shock and Na+ concentrations accumulated with the increase of processing time under long-term salt stress conditions. The ionic stress component of salinity stress becomes gradually more severe after gradual salt application. Proline is known as a molecular chaperone to protect protein integrity and increase the enzyme activities under salt stress. It is produced from either glutamate or ornithine in plants, whereas in mature plants or during exposure to stress the glutamate pathway usually dominates [43]. In our study, proline levels in cotton leaves significantly increased under long-term salt stress. It seems that many genes involved in the accumulation of proline and other osmolytes should also be induced. Recent advances indicated that salt tolerance in plants is a qualitative and quantitative trait controlled by multiple genes [44,45]. It is essential to pay more attention to the dynamic activity of post-transcriptional and post-translational regulators in this field [7,46].
miRNA expression profiles in response to salt stress have been analyzed in a variety of plant species, including Arabidopsis, tobacco, rice, maize and poplar [7,16,24,47,48]. In this study, we monitored miRNAs and the expression of their target genes under salt shock and long-term salt stress. A total of 141 miRNAs exhibited altered expression. Among them, 20 and nine miRNAs showed consistent up- or down-regulation, respectively. Six miRNAs were specifically expressed in the 4-h sample. Three miRNAs were specifically expressed in the 5-d sample. It has long been apparent that environmental signals can modulate plant responses to detrimental conditions through changes in phytohormone concentrations [49]. The hormonal changes influence not only the normal biomass but also the adaptive response [50]. Patterns of gene expression were suggested to be different in response to salt shock and salt stress [2]. We studied the target genes of miRNAs by degradome sequencing. Most target genes of conserved miRNAs were transcription factors, playing crucial roles in gene regulation and involved in various aspects of plant growth and development. The majority of target genes of novel miRNAs were protein-coding genes involved in metabolism, protein modification and RNA or carbohydrate binding. miR160 and miR167 are well known as targeting members of the ARF family of transcription factors in the plant kingdom [51,52,53]. In our study, six miR160 genes and miR167g showed altered expression under salt stress. Using degradome sequencing, we detected fragments of cotton ARF8, ARF10, ARF16 and ARF17 mRNAs. They were precisely mapped to the end of each fragment of the nucleotide predicted for miR160- and miR167-directed cleavage. This analysis confirmed that the four cotton ARFs were post-transcriptionally regulated by miR160s and miR167s. ARF proteins can either activate or repress transcription, depending on the nature of the middle domain [54]. ARF10, ARF16 and ARF17 activate ARFs with Gln-Leu-Ser-rich middle regions, whereas, ARF6 represses ARFs with Pro-Ser-Thr-rich middle regions [55]. In our study, the expression of ARF10, ARF16 and ARF17 under salt stress was lower than under normal conditions. ARF8 was firstly up-regulated and then down-regulated by salt stress. It seems that the altered expression of ARF6, -10, -16 and -17 by miR160 or miR167 improved the level of auxin, and in turn reset the developmental program and compensated for the biomass damage due to salt stress.
tasiRNAs form a class of plant-specific endogenous small RNAs that function as miRNA-like post-transcriptional negative regulators [13]. In this study, two TAS3 genes, TAS3a and TAS3b, with two near-identical 21-nt tasiRNA-ARFs that coaligned with the phases D7(+) and D8(+), were identified in three Gossypium species by homology search. The same sequences suggested high conservation of TAS3 in Gossypium species. In Arabidopsis, it was recently demonstrated that miR390 cleaves the TAS3 precursor for the production of tasiRNA-ARF, which could cleave the transcripts of ARF3 (ETT) and its close homolog ARF4 [56]. For cotton, cleavage sites for TAS3 precursors and ARF4 were identified in our degradome sequencing results. In addition, cotton ARF4 mRNAs showed a significant negative correlation with the corresponding miR390. These observations suggested that miR390, TAS3 tasiRNAs and their ARF targets defined a stepwise pathway in response to salt stress in cotton (Figure 9C). The degree of conservation observed in land plants suggests that the tasiRNA-ARF pathway plays a fundamental role in plant development. In dicots, theTAS3 tasiRNA-ARFs regulate leaf patterning, developmental timing and rate of lateral root growth [57]. Defects in this pathway could cause phenotypic abnormalities in flower gynoecium formation and leaf heteroblasty in Arabidopsis [58,59]. Adaxial localization of tasiRNA-ARF correlates with the restriction of ARF3/4 localization in the abaxial domain during leaf development [58]. Among the 23 family members of ARFs in Arabidopsis, ARF3 and ARF4 are known to be transcriptional repressors [55]. As abaxial determinants, ARF3/4 have been suggested to be mediators for auxin signaling to partition adaxial and abaxial domains during leaf primordium formation [60]. Our results showed that cotton ARF4 was dramatically up-regulated. Thus, it seems that auxin concentrations were repressed by the regulation of tasiRNA-ARF-dependent ARF4 expression under salt stress. Among the available loss-of-function Arabidopsis mutants (for at least 18 ARF genes), only single mutants of arf2, arf3, arf5, arf7, arf8 and arf19 have phenotypes differing in growth or development [61,62,63,64]. Accordingly, further analysis of cotton tasiRNA-ARFs and their target genes could shed new light on their regulatory roles of auxin in salt stress response.

5. Conclusions

In this study, global transcriptional profiles of miRNA and other small non-coding RNAs were investigated in responses to salt shock and long-term salt stress. The salt-responsive miRNAs may function under salt stress in changing the normal level of regulated gene transcripts responsible for optimized developmental patterning, and thereby adjusting the durations of phenological phases. The identification of these small RNAs as well as elucidating their functional significance broadens our understanding of post-transcriptional gene regulation in salt stress response.

Supplementary Materials

The following are available online at www.mdpi.com/2073-4425/8/12/369/s1. Table S1: Primer pairs used in real-time quantitative RT-PCR, Figure S1: Distribution of small RNAs among different categories in three libraries, Figure S2: Mature and the predicted fold-back structures of newly identified miRNAs in cotton. Sequences of mature miRNAs are underlined, Figure S3: Dissociation Curve of seven miRNAs for Real-time quantitative PCR, Figure S4: Target sites in ARF2, ARF3 and ARF4 genes of G. arboreum and G. raimondii. Target site was shown as straight lines.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (31201247), Young Elite Scientists Sponsorship Program by CAST (2016QNRC001), and State Key Laboratory of Crop Biology Open Fund (2015KF12).

Author Contributions

Z.Y., Y.L. and W.Y. conceived and designed the experiments; Z.Y., J.W., D.W. and S.W. performed the experiments; Z.Y., X.H. and F.X. analyzed the data; Z.Y. and Y.L. wrote the paper. All authors read and approved the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Xiong, L.; Schumaker, K.S.; Zhu, J.K. Cell signaling during cold, drought, and salt stress. Plant Cell 2002, 14, S165–S183. [Google Scholar] [CrossRef] [PubMed]
  2. Shavrukov, Y. Salt stress or salt shock: Which genes are we studying? J. Exp. Bot. 2013, 64, 119–127. [Google Scholar] [CrossRef] [PubMed]
  3. Munns, R.; Tester, M. Mechanisms of salinity tolerance. Annu. Rev. Plant Biol. 2008, 59, 651–681. [Google Scholar] [CrossRef] [PubMed]
  4. Zhu, J.K. Regulation of ion homeostasis under salt stress. Curr. Opin. Plant Biol. 2003, 6, 441–445. [Google Scholar] [CrossRef]
  5. Munns, R. Genes and salt tolerance: Bringing them together. New Phytol. 2005, 167, 645–663. [Google Scholar] [CrossRef] [PubMed]
  6. Zhao, K.F.; Song, J.; Fan, H.; Zhou, S.; Zhao, M. Growth response to ionic and osmotic stress of NaCl in salt-tolerant and salt-sensitive maize. J. Integr. Plant Biol. 2010, 52, 468–475. [Google Scholar] [CrossRef] [PubMed]
  7. Li, B.; Duan, H.; Li, J.; Deng, X.W.; Yin, W.; Xia, X. Global identification of miRNAs and targets in Populus euphratica under salt stress. Plant Mol. Biol. 2013, 81, 525–539. [Google Scholar] [CrossRef] [PubMed]
  8. Busov, V.B.; Brunner, A.M.; Strauss, S.H. Genes for control of plant stature and form. New Phytol. 2008, 177, 589–607. [Google Scholar] [CrossRef] [PubMed]
  9. Schwab, R.; Palatnik, J.F.; Riester, M.; Schommer, C.; Schmid, M.; Weigel, D. Specific effects of microRNAs on the plant transcriptome. Dev. Cell 2005, 8, 517–527. [Google Scholar] [CrossRef] [PubMed]
  10. Zhang, B.; Pan, X.; Cannon, C.H.; Cobb, G.P.; Anderson, T.A. Conservation and divergence of plant microRNA genes. Plant J. 2006, 46, 243–259. [Google Scholar] [CrossRef] [PubMed]
  11. Axtell, M.J.; Snyder, J.A.; Bartel, D.P. Common functions for diverse small RNAs of land plants. Plant Cell 2007, 19, 1750–1769. [Google Scholar] [CrossRef] [PubMed]
  12. Jones-Rhoades, M.W.; Bartel, D.P.; Bartel, B. microRNAs and their regulatory roles in plants. Annu. Rev. Plant Biol. 2006, 57, 19–53. [Google Scholar] [CrossRef] [PubMed]
  13. Axtell, M.J. Classification and comparison of small RNAs from plants. Annu. Rev. Plant Biol. 2013, 64, 137–159. [Google Scholar] [CrossRef] [PubMed]
  14. Shukla, L.I.; Chinnusamy, V.; Sunkar, R. The role of microRNAs and other endogenous small RNAs in plant stress responses. Biochim. Biophys. Acta 2008, 1779, 743–748. [Google Scholar] [CrossRef] [PubMed]
  15. Nonogaki, H. MicroRNA gene regulation cascades during early stages of plant development. Plant Cell Physiol. 2010, 51, 1840–1846. [Google Scholar] [CrossRef] [PubMed]
  16. Sunkar, R.; Zhu, J.K. Novel and stress-regulated microRNAs and other small RNAs from Arabidopsis. Plant Cell 2004, 16, 2001–2019. [Google Scholar] [CrossRef] [PubMed]
  17. Chen, Z.; Hu, L.; Han, N.; Hu, J.; Yang, Y.; Xiang, T.; Zhang, X.; Wang, L. Overexpression of a miR393-resistant form of transport inhibitor response protein 1 (mTIR1) enhances salt tolerance by increased osmoregulation and Na+ exclusion in Arabidopsis thaliana. Plant Cell Physiol. 2015, 56, 73–83. [Google Scholar] [CrossRef] [PubMed]
  18. Gao, P.; Bai, X.; Yang, L.; Lv, D.; Pan, X.; Li, Y.; Cai, H.; Ji, W.; Chen, Q.; Zhu, Y. Osa-miR393: A salinity- and alkaline stress-related microRNA gene. Mol. Biol. Rep. 2010, 38, 237–242. [Google Scholar] [CrossRef] [PubMed]
  19. Xie, F.; Wang, Q.; Sun, R.; Zhang, B. Deep sequencing reveals important roles of microRNAs in response to drought and salinity stress in cotton. J. Exp. Bot. 2015, 66, 789–804. [Google Scholar] [CrossRef] [PubMed]
  20. Wang, M.; Wang, Q.; Zhang, B. Response of miRNAs and their targets to salt and drought stresses in cotton (Gossypium hirsutum L.). Gene 2013, 530, 26–32. [Google Scholar] [CrossRef] [PubMed]
  21. Sunkar, R.; Chinnusamy, V.; Zhu, J.; Zhu, J.K. Small RNAs as big players in plant abiotic stress responses and nutrient deprivation. Trends Plant Sci. 2007, 12, 301–309. [Google Scholar] [CrossRef] [PubMed]
  22. Bari, R.; Datt Pant, B.; Stitt, M.; Scheible, W.R. PHO2, microRNA399, and PHR1 define a phosphate-signaling pathway in plants. Plant Physiol. 2006, 141, 988–999. [Google Scholar] [CrossRef] [PubMed]
  23. Yin, Z.; Li, Y.; Yu, J.; Liu, Y.; Li, C.; Han, X.; Shen, F. Difference in miRNA expression profiles between two cotton cultivars with distinct salt sensitivity. Mol. Biol. Rep. 2012, 39, 4961–4970. [Google Scholar] [CrossRef] [PubMed]
  24. Ding, D.; Zhang, L.; Wang, H.; Liu, Z.; Zhang, Z.; Zheng, Y. Differential expression of miRNAs in response to salt stress in maize roots. Ann. Bot. 2009, 103, 29–38. [Google Scholar] [CrossRef] [PubMed]
  25. Leidi, E.O.; Saiz, J.F. Is salinity tolerance related to Na+ accumulation in upland cotton (Gossypium hirsutum) seedlings? Plant Soil 1997, 190, 67–75. [Google Scholar] [CrossRef]
  26. Carra, A.; Gambino, G.; Schubert, A. A cetyltrimethylammonium bromide-based method to extract low-molecular-weight RNA from polysaccharide-rich plant tissues. Anal. Biochem. 2007, 360, 318–320. [Google Scholar] [CrossRef] [PubMed]
  27. Yin, Z.; Li, Y.; Han, X.; Shen, F. Genome-wide profiling of miRNAs and other small non-coding RNAs in the Verticillium dahliae-inoculated cotton roots. PLoS ONE 2012, 7, e35765. [Google Scholar] [CrossRef] [PubMed]
  28. Hafner, M.; Landgraf, P.; Ludwig, J.; Rice, A.; Ojo, T.; Lin, C.; Holoch, D.; Lim, C.; Tuschl, T. Identification of microRNAs and other small regulatory RNAs using cDNA library sequencing. Methods 2008, 44, 3–12. [Google Scholar] [CrossRef] [PubMed]
  29. German, M.A.; Pillay, M.; Jeong, D.H.; Hetawal, A.; Luo, S.; Janardhanan, P.; Kannan, V.; Rymarquis, L.A.; Nobuta, K.; German, R.; et al. Global identification of microRNA-target RNA pairs by parallel analysis of RNA ends. Nat. Biotechnol. 2008, 26, 941–946. [Google Scholar] [CrossRef] [PubMed]
  30. Grabherr, M.G.; Haas, B.J.; Yassour, M.; Levin, J.Z.; Thompson, D.A.; Amit, I.; Adiconis, X.; Fan, L.; Raychowdhury, R.; Zeng, Q.; et al. Trinity: Reconstructing a full-length transcriptome without a genome from RNA-Seq data. Nat. Biotechnol. 2011, 29, 644–652. [Google Scholar] [CrossRef] [PubMed]
  31. Rfam. Available online: http://www.sanger.ac.uk/resources/databases/rfam.html (accessed on 9 January 2013).
  32. Kozomara, A.; Griffiths-Jones, S. miRBase: Integrating microRNA annotation and deep-sequencing data. Nucleic Acids Res. 2011, 39, D152–D157. [Google Scholar] [CrossRef] [PubMed]
  33. MIREAP program. Available online: https://sourceforge.net/projects/mireap/ (accessed on 19 December 2012).
  34. Yin, Z.; Li, C.; Han, X.; Shen, F. Identification of conserved microRNAs and their target genes in tomato (Lycopersicon esculentum). Gene 2008, 414, 60–66. [Google Scholar] [CrossRef] [PubMed]
  35. MFOLD software. Available online: http://mfold.rna.albany.edu/ (accessed on 15 June 2013).
  36. Zuker, M. Mfold web server for nucleic acid folding and hybridization prediction. Nucleic Acids Res. 2003, 31, 3406–3415. [Google Scholar] [CrossRef] [PubMed]
  37. Chen, C.; Ridzon, D.A.; Broomer, A.J.; Zhou, Z.; Lee, D.H.; Nguyen, J.T.; Barbisin, M.; Xu, N.L.; Mahuvakar, V.R.; Andersen, M.R.; et al. Real-time quantification of microRNAs by stem-loop RT-PCR. Nucleic Acids Res. 2005, 33, e179. [Google Scholar] [CrossRef] [PubMed]
  38. Varkonyi-Gasic, E.; Wu, R.; Wood, M.; Walton, E.F.; Hellens, R.P. Protocol: A highly sensitive RT-PCR method for detection and quantification of microRNAs. Plant Methods 2007, 3, 12. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  39. Macovei, A.; Tuteja, N. MicroRNAs targeting DEAD-box helicases are involved in salinity stress response in rice (Oryza sativa L.). BMC Plant Biol. 2012, 12, 183. [Google Scholar] [CrossRef] [PubMed]
  40. Zhang, B.H.; Pan, X.P.; Cox, S.B.; Cobb, G.P.; Anderson, T.A. Evidence that miRNAs are different from other RNAs. Cell. Mol. Life Sci. 2006, 63, 246–254. [Google Scholar] [CrossRef] [PubMed]
  41. Unver, T.; Namuth-Covert, D.M.; Budak, H. Review of current methodological approaches for characterizing microRNAs in plants. Int. J. Plant Genom. 2009, 2009, 262463. [Google Scholar] [CrossRef] [PubMed]
  42. Hawkins, J.S.; Kim, H.; Nason, J.D.; Wing, R.A.; Wendel, J.F. Differential lineage-specific amplification of transposable elements is responsible for genome size variation in Gossypium. Genome Res. 2006, 16, 1252–1261. [Google Scholar] [CrossRef] [PubMed]
  43. Kavi Kishor, P.B.; Sreenivasulu, N. Is proline accumulation per se correlated with stress tolerance or is proline homeostasis a more critical issue? Plant Cell Environ. 2014, 37, 300–311. [Google Scholar] [CrossRef] [PubMed]
  44. Parida, A.K.; Das, A.B. Salt tolerance and salinity effects on plants: A review. Ecotoxicol. Environ. Saf. 2005, 60, 324–349. [Google Scholar] [CrossRef] [PubMed]
  45. Munns, R.; James, R.A.; Lauchli, A. Approaches to increasing the salt tolerance of wheat and other cereals. J. Exp. Bot. 2006, 57, 1025–1043. [Google Scholar] [CrossRef] [PubMed]
  46. Xie, F.; Stewart, C.N.; Taki, F.A.; He, Q.; Liu, H.; Zhang, B. High-throughput deep sequencing shows that microRNAs play important roles in switchgrass responses to drought and salinity stress. Plant Biotechnol. J. 2014, 12, 354–366. [Google Scholar] [CrossRef] [PubMed]
  47. Barrera-Figueroa, B.E.; Gao, L.; Wu, Z.; Zhou, X.; Zhu, J.; Jin, H.; Liu, R.; Zhu, J.K. High throughput sequencing reveals novel and abiotic stress-regulated microRNAs in the inflorescences of rice. BMC Plant Biol. 2012, 12, 132. [Google Scholar] [CrossRef] [PubMed]
  48. Frazier, T.P.; Sun, G.; Burklew, C.E.; Zhang, B. Salt and drought stresses induce the aberrant expression of microRNA genes in tobacco. Mol. Biotechnol. 2011, 49, 159–165. [Google Scholar] [CrossRef] [PubMed]
  49. Junghans, U.; Polle, A.; Duchting, P.; Weiler, E.; Kuhlman, B.; Gruber, F.; Teichmann, T. Adaptation to high salinity in poplar involves changes in xylem anatomy and auxin physiology. Plant Cell Environ. 2006, 29, 1519–1531. [Google Scholar] [CrossRef] [PubMed]
  50. Mockaitis, K.; Estelle, M. Auxin receptors and plant development: A new signaling paradigm. Annu. Rev. Cell Dev. Biol. 2008, 24, 55–80. [Google Scholar] [CrossRef] [PubMed]
  51. Mallory, A.C.; Bartel, D.P.; Bartel, B. MicroRNA-directed regulation of ArabidopsisAUXIN RESPONSE FACTOR17 is essential for proper development and modulates expression of early auxin response genes. Plant Cell 2005, 17, 1360–1375. [Google Scholar] [CrossRef] [PubMed]
  52. Wu, M.F.; Tian, Q.; Reed, J.W. ArabidopsismicroRNA167 controls patterns of ARF6 and ARF8 expression, and regulates both female and male reproduction. Development 2006, 133, 4211–4218. [Google Scholar] [CrossRef] [PubMed]
  53. Yang, J.H.; Han, S.J.; Yoon, E.K.; Lee, W.S. Evidence of an auxin signal pathway, microRNA167-ARF8-GH3, and its response to exogenous auxin in cultured rice cells. Nucleic Acids Res. 2006, 34, 1892–1899. [Google Scholar] [CrossRef] [PubMed]
  54. Ulmasov, T.; Hagen, G.; Guilfoyle, T.J. Activation and repression of transcription by auxin-response factors. Proc. Natl. Acad. Sci. USA 1999, 96, 5844–5849. [Google Scholar] [CrossRef] [PubMed]
  55. Tiwari, S.B.; Hagen, G.; Guilfoyle, T. The roles of auxin response factor domains in auxin-responsive transcription. Plant Cell 2003, 15, 533–543. [Google Scholar] [CrossRef] [PubMed]
  56. Fahlgren, N.; Montgomery, T.A.; Howell, M.D.; Allen, E.; Dvorak, S.K.; Alexander, A.L.; Carrington, J.C. Regulation of AUXIN RESPONSE FACTOR3 by TAS3 ta-siRNA affects developmental timing and patterning in Arabidopsis. Curr. Biol. 2006, 16, 939–944. [Google Scholar] [CrossRef] [PubMed]
  57. Yoon, E.K.; Yang, J.H.; Lim, J.; Kim, S.H.; Kim, S.K.; Lee, W.S. Auxin regulation of the microRNA390-dependent transacting small interfering RNA pathway in Arabidopsis lateral root development. Nucleic Acids Res. 2010, 38, 1382–1391. [Google Scholar] [CrossRef] [PubMed]
  58. Hunter, C.; Willmann, M.R.; Wu, G.; Yoshikawa, M.; de la Luz Gutierrez-Nava, M.; Poethig, S.R. Trans-acting siRNA-mediated repression of RTTIN and ARF4 regulates heteroblasty in Arabidopsis. Development 2006, 133, 2973–2981. [Google Scholar] [CrossRef] [PubMed]
  59. Nemhauser, J.L.; Feldman, L.J.; Zambryski, P.C. Auxin and ETTIN in Arabidopsis gynoecium morphogenesis. Development 2000, 127, 3877–3888. [Google Scholar] [PubMed]
  60. Pekker, I.; Alvarez, J.P.; Eshed, Y. Auxin response factors mediate Arabidopsis organ asymmetry via modulation of KANADI activity. Plant Cell 2005, 17, 2899–2910. [Google Scholar] [CrossRef] [PubMed]
  61. Okushima, Y.; Overvoorde, P.J.; Arima, K.; Alonso, J.M.; Chan, A.; Chang, C.; Ecker, J.R.; Hughes, B.; Lui, A.; Nguyen, D.; et al. Functional genomic analysis of the AUXIN RESPONSE FACTOR gene family members in Arabidopsis thaliana: Unique and overlapping functions of ARF7 and ARF19. Plant Cell 2005, 17, 444–463. [Google Scholar] [CrossRef] [PubMed]
  62. Schruff, M.C.; Spielman, M.; Tiwari, S.; Adams, S.; Fenby, N.; Scott, R.J. The AUXIN RESPONSE FACTOR 2 gene of Arabidopsis links auxin signalling, cell division, and the size of seeds and other organs. Development 2006, 133, 251–261. [Google Scholar] [CrossRef] [PubMed]
  63. Nagpal, P.; Ellis, C.M.; Weber, H.; Ploense, S.E.; Barkawi, L.S.; Guilfoyle, T.J.; Hagen, G.; Alonso, J.M.; Cohen, J.D.; Farmer, E.E.; et al. Auxin response factors ARF6 and ARF8 promote jasmonic acid production and flower maturation. Development 2005, 132, 4107–4118. [Google Scholar] [CrossRef] [PubMed]
  64. Harper, R.M.; Stowe-Evans, E.L.; Luesse, D.R.; Muto, H.; Tatematsu, K.; Watahiki, M.K.; Yamamoto, K.; Liscum, E. The NPH4 locus encodes the auxin response factor ARF7, a conditional regulator of differential growth in aerial Arabidopsis tissue. Plant Cell 2000, 12, 757–770. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Measurement of Na+ content and proline level under salt stress. (A) Na+ content; (B) Proline content. The different letters present on the columns indicate significant differences at p < 0.01 among the control samples (CK), 4-h and 5-d salt-treatment samples. Error bars represent standard error of the mean.
Figure 1. Measurement of Na+ content and proline level under salt stress. (A) Na+ content; (B) Proline content. The different letters present on the columns indicate significant differences at p < 0.01 among the control samples (CK), 4-h and 5-d salt-treatment samples. Error bars represent standard error of the mean.
Genes 08 00369 g001
Figure 2. Differentially expressed miRNAs between libraries.
Figure 2. Differentially expressed miRNAs between libraries.
Genes 08 00369 g002
Figure 3. Real-time quantitative PCR validated expressions of seven miRNAs. The amount of expression was normalized to the level of U6. The normalized miRNA levels in control were arbitrarily set to 1. The different letters of each miRNAs present on the columns indicate significant differences at p < 0.05 among the CK, 4-h and 5-d samples. Error bars represent standard error of the mean.
Figure 3. Real-time quantitative PCR validated expressions of seven miRNAs. The amount of expression was normalized to the level of U6. The normalized miRNA levels in control were arbitrarily set to 1. The different letters of each miRNAs present on the columns indicate significant differences at p < 0.05 among the CK, 4-h and 5-d samples. Error bars represent standard error of the mean.
Genes 08 00369 g003
Figure 4. Distribution of small RNAs on three expressed retrotransposons.
Figure 4. Distribution of small RNAs on three expressed retrotransposons.
Genes 08 00369 g004
Figure 5. Relative expression levels of two retrotransposons under salt stress. The different letters of each retrotransposons present on the columns indicate significant differences at p < 0.05 among the CK, 4-h and 5-d samples. Error bars represent standard error of the mean.
Figure 5. Relative expression levels of two retrotransposons under salt stress. The different letters of each retrotransposons present on the columns indicate significant differences at p < 0.05 among the CK, 4-h and 5-d samples. Error bars represent standard error of the mean.
Genes 08 00369 g005
Figure 6. Target plot (t-plot) for target genes of less-conserved miRNAs. Arrows indicate the signatures corresponding to the miRNA cleavage site. Partial mRNA sequences of target genes aligned with the miRNAs show perfect matches (straight lines), G-U wobbles (circles).
Figure 6. Target plot (t-plot) for target genes of less-conserved miRNAs. Arrows indicate the signatures corresponding to the miRNA cleavage site. Partial mRNA sequences of target genes aligned with the miRNAs show perfect matches (straight lines), G-U wobbles (circles).
Genes 08 00369 g006
Figure 7. Using qRT-PCR, the expression profiles of ten target genes were examined under salt stress. There were three biological replicates and three technical replicates. The different letters of each target genes present on the columns indicate significant differences at p < 0.05 among the CK, 4-h and 5-d samples. Error bars represent standard error of the mean.
Figure 7. Using qRT-PCR, the expression profiles of ten target genes were examined under salt stress. There were three biological replicates and three technical replicates. The different letters of each target genes present on the columns indicate significant differences at p < 0.05 among the CK, 4-h and 5-d samples. Error bars represent standard error of the mean.
Genes 08 00369 g007
Figure 8. Target plot (t-plot) for cotton TAS3b (A) and ARF4 (B). Arrows indicate the signatures corresponding to the cleavage site. Partial mRNA sequences of target genes aligned with small RNA show perfect matches (straight lines), G-U wobbles (circles).
Figure 8. Target plot (t-plot) for cotton TAS3b (A) and ARF4 (B). Arrows indicate the signatures corresponding to the cleavage site. Partial mRNA sequences of target genes aligned with small RNA show perfect matches (straight lines), G-U wobbles (circles).
Genes 08 00369 g008
Figure 9. miR390/tasiRNA-ARF/ARF4 pathway analysis. Error bars represent standard error of the mean. (A) Relative expression levels of miR390 and TAS3 under salt stress; (B) Relative expression levels of tasiRNA-ARFs and ARF4 under salt stress; (C) A model for the role of the miR390/tasiRNA-ARF/ARF4 module under salt stress.
Figure 9. miR390/tasiRNA-ARF/ARF4 pathway analysis. Error bars represent standard error of the mean. (A) Relative expression levels of miR390 and TAS3 under salt stress; (B) Relative expression levels of tasiRNA-ARFs and ARF4 under salt stress; (C) A model for the role of the miR390/tasiRNA-ARF/ARF4 module under salt stress.
Genes 08 00369 g009
Table 1. Dataset summary of sequencing of small RNA, degradome and transcriptome libraries.
Table 1. Dataset summary of sequencing of small RNA, degradome and transcriptome libraries.
CategorySmall RNADegradomeTranscriptome
CK4 h5 dCK4 h5 d
Total reads17,533,90316,853,76922,982,88222,766,83930,703,66220,880,67957,194,024
High quality17,389,91616,704,86822,900,93122,744,20930,684,33120,837,65554,261,421
Clean reads17,026,07015,731,47722,590,39122,653,37230,649,47220,685,80052,170,788
Unigene sequences4,331,0283,703,5944,802,868864,7201,217,7571,171,57781,114
CK: Control Samples.
Table 2. List of the 24 novel miRNAs in cotton.
Table 2. List of the 24 novel miRNAs in cotton.
NameMature miRNA SequenceGene IDLM (nt)ArmLP (nt)G+C(%)miRNA* SequenceΔG kcal/molMFEIs
miR1327TGAATATTGTTAAAGTAGAAAUnigene15286213′12227.87TCTACTTTAACAATATTCATA−43.301.27
miR1335TGAAGCTGCCACATGATCTAGUnigene16386215′10346.60AGATCATGCTGGCAGCTTCAAC−50.201.05
miR1336TCTTTCCTACTCCTCCCATTCCCL1064.Contig1223′10539.05AAGTGGGATGGGTGGAAAGATT−44.501.09
miR1337AATGGAGGAGTTGGAAAGATTUnigene26231215′9528.42TTTCCAATTCCTCCCATTCCAC−36.701.36
miR1338AATGAGTTAGGCGAGAGGTTCCL5909.Contig2213′16931.36CTCTGGTTTAACTCATTTGTA−51.100.96
miR1339AAAAGCAATGAAGAACGACCTCL5831.Contig1213′10238.24GCTGCGTTCTTCATTGCTTAACA−41.301.06
miR1340ACCGGCCGGGGGACGGACTGGGCL6985.Contig2223′11965.55CAGTCCCGAACCCGTCAGCTGC−62.600.80
miR1341CTTTTTATAGGATAGGGACTGUnigene30026213′10329.13GTTCCTATCTTATTAAAAAGAA−37.101.24
miR1342TGAAGCTGCAAGATGATCTGACL8098.Contig4215′26234.35AGATCATGCTGGCAGCTTCAAC−57.700.64
miR1343CAGAGATCGTCAAAGCATATCCL1118.Contig5215′9437.23TATGCTTTGACGATCTCTGAT−60.901.74
miR1344TTTGCATGAACTAGGAGACGTUnigene10644213′10554.29CGGCTGCTAGTTCATGGATGCC−45.600.80
miR1345AAGCGGAAGCTGAATTAGTTGUnigene14674213′7643.42GAAACTAATCGAGCTCCGTTTGA−26.000.79
miR1346GCTGCCATCTCATGCATTCGGUnigene20900215′16544.85CGATGCATGGCATGGGAGCACCA−59.100.80
miR1347TACAGCTTTAGAAATCATCCCTUnigene13628225′10929.36GGATGATTTCTAAAGCTCTAGA−48.401.51
miR1348AGTGTCTGGGTGGTGTAGTTGGTUnigene4349235′8350.60CCATTCGAACACGGGCTCAGACATTT−23.900.65
miR1349TAACTTGTCTTCGCCCTTCTCUnigene8023213′20134.83GAAGTGCATGGCAAGTTAGA−48.700.70
miR1350TGGCACGGCTCAATCAAATTAUnigene7762215′8938.20ATTTGATTGAGCCATTCCAAC−34.601.02
miR1351ACTCATAATTTAGCAAAGTCGUnigene8111213′12928.68TTGCTGGATTATGAGTCTAAT−61.401.66
miR1352AATGCTTGAGGTGATAGGTTCAUnigene14827225′13750.37AACCATTGCCTCAAGCACTTG−58.800.85
miR1353AAGGCAAAGGAAGAAAAGAGTGACL9674.Contig1235′10538.10ACTCTTTTTTTCCTGCCTTGC−45.901.15
miR1354AAAGTGGATGAAATTTTTAGCUnigene26528213′15936.48TAAAAATTTCATTCATTTCTA−66.101.14
miR1355ATGCACTGCCTCTTCCCTGGCUnigene10226213′10253.92AACAGGCTGAGCATGGATGGA−43.400.79
miR1356AACTGTGAAGCTATAAGGTATUnigene19434215′11037.27ACTTTTGTTTCACATTAC−23.900.65
miR1357AAGCTGTTGATGGCCGGCATGACL9715.Contig1225′9146.15ATGCCTATCATCAGGAGACTCT−30.300.72
LM: length of mature miRNAs; LP: length of precursor; MFEIs: minimal folding free energy indexes.
Table 3. Target genes of miRNAs identified in cotton.
Table 3. Target genes of miRNAs identified in cotton.
miRNA FamilyTarget GeneAlignment ScoreAlignment RangeCleavage SiteCategoryAnnotationAbundance
CK4 h5 d
miR160Unigene5547 (CotAD_09605)4378–3983892ARF 10361254207
CL11606.Contig2 (CotAD_08070)4198–2182082ARF1636154207
CL11606.Contig1 (CotAD_76501)2.5345–3663572ARF1627942332249
CL1202.Contig4 (CotAD_55881)1497–5185092ARF1742083394
CL1202.Contig2 (CotAD_36520)22–22123ARF1747683307
miR164CL1985.Contig3 (CotAD_58864)4.519–39302NAC0610
CL1985.Contig5 (CotAD_43532)4.519–39302NAC0610
miR165CL795.Contig6 (CotAD_15207)493–1131033Class III HD-Zip protein 4111063
miR167Unigene18486 (CotAD_40344)4.53–23152ARF6450125
miR169CL599.Contig2 (CotAD_06772)21095–111511053Nuclear transcription factor Y subunit A-3151688
miR172CL799.Contig8 (CotAD_29570)0761–7827732AP2-EREBP189104277
miR319CL11426.Contig1 (CotAD_73356)2361–3813732MYB84295398
miR396Unigene30805 (CotAD_25614)22–22122GRF0545
miR828Unigene27368 (CotAD_00308)11–22130MYB5101311060
miR894CL861.Contig1 (CotAD_07539)42507–252825214Leucine Rich Repeat Kinase (LRRK)0335
miR1327CL4730.Contig13.5782–8037953Simple Sequence Repeat Marker, mRNA sequence (SSR)0572
miR1344CL10222.Contig14429–4504384Uncharacterized protein47246
miR1355CL9524.Contig2 (CotAD_06422)3.555–76623Copper binding protein 6 (BCP6)58014
miR8155CL11661.Contig1 (CotAD_51244)4735–7537414AUX/IAA105261
miR5083CL11334.Contig1 (CotAD_47266)3.5334–3563412Unknown protein51243
CL5215.Contig1 (CotAD_47976)3.51524–154615363Nucleosome assembly protein73113
miR7484CL11292.Contig4 (CotAD_04736)42–24172Predicted protein0216
miR7502Unigene26285 (CotAD_22251)3.511–32234Trihelix010
CL503.Contig1 (CotAD_21495)4.518–39304Allene Oxide Cyclase (AOC)0120
miR7505Unigene19150 (CotAD_63955)4533–5545453NADH-ubiquinone oxidoreductase subunit 9 (NAD9)7322951
miR7508Unigene9225 (CotAD_26989)4.51–21122Pectin Methylesterase (PME2)3138650
CL11591.Contig1 (CotAD_08283)4.51034–105410454Translocon-associated protein subunit alpha72446
CL1882.Contig3 (CotAD_18999)416–36284Zinc finger family protein,063
miR7512CL6.Contig3 (CotAD_04930)4.513–35264Uncharacterized protein020
CL6.Contig4 (CotAD_04930)4.513–35264Uncharacterized protein020
Unigene19759412–31234Alpha/beta-Hydrolases superfamily protein isoform010
Table 4. Abundance of four tasiRNA-ARFs and genomic position of two TAS3 genes in three Gossypium species.
Table 4. Abundance of four tasiRNA-ARFs and genomic position of two TAS3 genes in three Gossypium species.
NameSequencesAbundanceGenomic Position
CK4 h5 dG. raimondiiG. arboreumG. hirsutum
TAS3a D8(+)TTCTTGACCTTGTAAGGCCTT221Chr9: 41,462,977–41,462,936scaffold13100:41–82At_chr3: 23,630,524–23,630,565
TAS3a D7(+)TTCTTGACCTTGTAAGACCCC22812291Dt_chr8: 37,518,786–37,518,827
TAS3b D8(+)TTCTTGACCTTGTAAGACCTT816897541scaffold493: 67,555–67,595Ca8:60,478,108–60,478,068Dt_chr11: 54,832,692–54,832,651
TAS3b D7(+)TTCTTGACCTTGTAAGACCCA11,71811,0607524At_chr11: 5,189,802–5,189,761
Table 5. List of the ARF2, ARF3 and ARF4 in three Gossypium species and the target sites for tasiRNA-ARFs.
Table 5. List of the ARF2, ARF3 and ARF4 in three Gossypium species and the target sites for tasiRNA-ARFs.
NameG. raimondiiG. arboreumG. hirsutum
Gene identifierTarget SitesGene IdentifierTarget SitesGene IdentifierTarget Sites
ARF2Cotton_D_10022860
Cotton_D_10033064
1328–1351
1331–1352
Cotton_A_03644
Cotton_A_01955
1325–1346
1331–1352
CotAD_43230
CotAD_76410
CotAD_58764
1307–1329
1331–1353
1331–1353
ARF3Cotton_D_10031398
Cotton_D_10033714
1254–1275
1269–1290
1464–1485
1479–1500
Cotton_A_11311
Cotton_A_03933
1272–1293
1320–1341
1482–1503
1530–1551
CotAD_00650
CotAD_64060
CotAD_72124
1254–1275
1281–1303
1104–1125
1464–1485
1491–1512
1314–1335
ARF4Cotton_D_10015852
Cotton_D_10019010
1368–1389
1350–1371
1575–1596
1554–1575
Cotton_A_01738
Cotton_A_11048
1197–1218
1359–1380
1404–1425
1563–1584
CotAD_13694
CotAD_54309
CotAD_12422
1368–1389
1350–1371
1359–1380
1575–1596
1554–1575
1563–1584

Share and Cite

MDPI and ACS Style

Yin, Z.; Han, X.; Li, Y.; Wang, J.; Wang, D.; Wang, S.; Fu, X.; Ye, W. Comparative Analysis of Cotton Small RNAs and Their Target Genes in Response to Salt Stress. Genes 2017, 8, 369. https://0-doi-org.brum.beds.ac.uk/10.3390/genes8120369

AMA Style

Yin Z, Han X, Li Y, Wang J, Wang D, Wang S, Fu X, Ye W. Comparative Analysis of Cotton Small RNAs and Their Target Genes in Response to Salt Stress. Genes. 2017; 8(12):369. https://0-doi-org.brum.beds.ac.uk/10.3390/genes8120369

Chicago/Turabian Style

Yin, Zujun, Xiulan Han, Yan Li, Junjuan Wang, Delong Wang, Shuai Wang, Xiaoqiong Fu, and Wuwei Ye. 2017. "Comparative Analysis of Cotton Small RNAs and Their Target Genes in Response to Salt Stress" Genes 8, no. 12: 369. https://0-doi-org.brum.beds.ac.uk/10.3390/genes8120369

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop