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

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.


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.

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 H 2 SO 4 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.

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 HNO 3 , 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.

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.

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.

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 |log 2 Ratio |>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.

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 H 2 O. 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.

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.

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. to the six nucleotides at the miRNA 3′-end, were combined with the 44-nt sequence to form a stemloop 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.

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.

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 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.

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]. 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.

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, DX402782and 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).

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).

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, DX402782and 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).

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. 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 negativecorrelation 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.

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.
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).

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

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

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.

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.