Genome-wide identification and characterization of microRNAs differentially expressed in fibers in a cotton phytochrome A1 RNAi line

Cotton fiber is an important commodity throughout the world. Fiber property determines fiber quality and commercial values. Previous studies showed that silencing phytochrome A1 gene (PHYA1) by RNA interference in Upland cotton (Gossypium hirsutum L. cv. Coker 312) had generated PHYA1 RNAi lines with simultaneous improvements in fiber quality (longer, stronger and finer fiber) and other key agronomic traits. Characterization of the altered molecular processes in these RNAi genotypes and its wild-type controls is a great interest to better understand the PHYA1 RNAi phenotypes. In this study, a total of 77 conserved miRNAs belonging to 61 families were examined in a PHYA1 RNAi line and its parental Coker 312 genotype by using multiplex sequencing. Of these miRNAs, seven (miR7503, miR7514, miR399c, miR399d, miR160, miR169b, and miR2950) were found to be differentially expressed in PHYA1 RNAi cotton. The target genes of these differentially expressed miRNAs were involved in the metabolism and signaling pathways of phytohormones, which included Gibberellin, Auxin and Abscisic Acid. The expression of several MYB transcription factors was also affected by miRNAs in RNAi cotton. In addition, 35 novel miRNAs (novel miR1-novel miR35) were identified in fibers for the first time in this study. Target genes of vast majority of these novel miRNAs were also predicted. Of these, nine novel miRNAs (novel-miR1, 2, 16, 19, 26, 27, 28, 31 and 32) were targeted to cytochrome P450-like TATA box binding protein (TBP). The qRT-PCR confirmed expression levels of several differentially regulated miRNAs. Expression patterns of four miRNAs-targets pairs were also examined via RNA deep sequencing. Together, the results imply that the regulation of miRNA expression might confer to the phenotype of the PHYA1 RNAi line(s) with improved fiber quality.

Introduction arboretum, AA, 2n = 2x = 26) and a local D-genome progenitor (Gossypium raimondii, DD, 2n = 2x = 26). The genome sequences of two diploid cotton Gossypium raimondii [27] and Gossypium arboretum [28], and the allotetraploid Gossypium hirsutum [29,30] have been determined. Whole genome sequences will definitely accelerate the identification and annotation of cotton miRNAs, especially cotton-specific miRNAs. Cotton miRNAs were first identified in 2007 by two groups using a comparative genomic approach [31,32]. Thirty and thirtyseven miRNAs candidates were identified by the two groups, respectively. Besides homology searches on database, Abdurakhmonov et al. [33] used a direct cloning technique to identify three miRNAs in cotton. Recently high-throughput RNA sequencing have been adopted for identification and expression analysis of miRNAs in cotton. Through small RNA sequencing, 34 conserved miRNAs families were identified in G. hirsutum [34] and 22 conserved miRNA families were expressed in developing cotton ovules [35]. By using the same approach, Pang et al. identified 4 novel miRNA families and proved that miRNAs play a role in cotton fiber development [26]. Recently, 65 conserved miRNA families were identified in cotton leaf and ovule tissues and forty of the identified miRNAs were ovule specific, suggesting that these miRNAs may play important roles in ovule and fiber development [34,36]. Using short fiber mutants, Naoumkina et al. [37] identified 24 conserved and 147 novel miRNA families and revealed that 4 miRNAs were involved in fiber elongation. In comparison to other plant species such as Arabidopsis or rice, less miRNAs had been identified in cotton. With its larger genome size, it is predicted that the cotton genome may encode more than the typical number of miRNAs, including some with unique roles in fiber development.
With the aim to better understand the altered molecular pathways in the PHYA1 RNAi line of tetraploid cotton, and because of potential importance of miRNAs in the complex developmental processes such as fiber quality, miRNA libraries were constructed in this study using fiber RNAs (at different development stages) extracted from the PHYA1 RNAi cotton line and its non-transgenic background Coker 312. The libraries were subjected to multiplex sequencing on the Illumina platform, and conserved and novel miRNAs were then identified and used for differential expression analysis. As results, a total of 77 known miRNAs from 61 different families and 35 novel miRNAs were identified in developing cotton fibers. Of these, 7 miRNAs were differentially expressed in the PHYA1 RNAi line. Our results implied that the regulation of miRNA expression might play an important role in cotton fiber development.

Plant materials and plant growth
Gossypium hirsutum L. cv. Coker 312 and its phytochrome A1 (PHYA1) RNAi plants were grown in the field at the North Farm R. R. Foil Plant Science Research Center of Mississippi State University. The seeds of these RNAi cotton plants [2] were provided by Uzbekistan collaborator of this study through Technology Transfer Office of USDA-ARS, USA under the USDA-Uzbekistan cooperation programs. Flowers were tagged on the day of anthesis (0 DPA). Cotton bolls (5, 10, 15 DPA) were harvested with three biological replications. Fibers were carefully dissected from cotton bolls, immediately frozen in liquid nitrogen, and then stored at -80˚C for small RNA extraction. used to convert isolated RNAs (500 ng) into barcoded small RNA libraries according to the manufacturer's instructions. To determine the differential expression of miRNAs in fibers between the PHYA1 RNAi line and its non-transgenic background line Coker 312, 18 barcoded small RNA libraries with three biological replicates were constructed by using miRNAs extracted from fibers of the two lines at three different developmental stages of 5-, 10-, and 15-DPA, respectively. Size selection (~150 nt) of the libraries was performed using AMPure XP beads. The concentration and quality of the libraries were determined with a Qubit Fluorimeter (ThermoFisher Scientific) and the Bioanalyzer 2100 (Agilent Technologies). The libraries were pooled and single-end (1 X 50) sequenced on an Illumina HiSeq 2500 system performed by BGI. The sequencing reads were deposited into NCBI Sequence Read Archive (SRA) with accession number SRX2467209-17 (Gossypium hirsutum L. cv. Coker 312) and SRX2467222-30 (RNAi line), respectively. Data files are available at https://www.ncbi.nlm.nih. gov/Traces/study/?acc=SRP096134 (Coker 312) and https://www.ncbi.nlm.nih.gov/Traces/ study/?acc=SRP096136 (RNAi line).

MicroRNA identification and target prediction
Total raw sequences generated from Illumina instrument were trimmed to remove lowquality reads and adaptor sequences using Trimmomatic [38], and FastQC (http://www. bioinformatics.bbsrc.ac.uk/projects/fastqc/) was used to test data quality. The miRBase database release 21 (http://www.mirbase.org/) and the Gossypium hirsutum L. acc. TM-1 genome sequence [30] were used as references for known conserved miRNA mapping and novel miRNA identification, respectively. SeqMan NGen and ArrayStar (DNASTAR, Inc.) were used for mapping of clean short reads (18-44 nt) to the reference and expression level analysis. The miRNA expression was normalized by the Reads Per Kilobase per Million mapped reads (RPKM) method [39].
Conserved miRNAs were identified by comparing to miRBase with the criterion that tags were similar to their homologues within two mismatches and without gaps [40]. Potential novel miRNAs were firstly BLAST verified against non-coding rRNA, tRNA, small nuclear RNA (snRNA), and small nucleolar RNA (snoRNA) in GenBank, tRNA database (http:// gtrnadb.ucsc.edu), and Rfam database [41]. Sequences with no hits to known miRNAs or other cellular RNAs were mapped to the Gossypium hirsutum L. acc. TM-1 genome. The mireap 0.2 (https://sourceforge.net/projects/mireap/) was used to analyze pre-miRNAs [42]. The psRNATarget was used for target prediction [43]. Sequence folding prediction were made by MFOLD [44]. miRNAs with negative hairpins folding energy from -25.7 to -230.5 kcal mol -1 [45] were selected for further analysis.

Quantitative real-time PCR
Quantitative Real time PCR was performed to validate the expression of miRNAs identified from RNA seq. The isolated fiber small RNAs (500 ng) were first reverse-transcribed to generate cDNAs using a Mir-X miRNA First-Strand Synthesis Kit (Clontech Laboratories). The cDNAs were then used as template for quantitative real-time PCR on LightCycler 480 (Roche Holding AG) by using SYBR Advantage qPCR Premix (Clontech Laboratories). The qRT-PCR was carried out by pre-denaturation at 95˚C for 15 seconds, following by 45 cycles of 3-step PCR with denaturation at 95˚C for 5 seconds, annealing at 55-60˚C for 15 seconds, and final extension at 72˚C for 15 seconds. Gene expression levels were presented as fold-change and calculated using the comparative threshold cycle (C T ) method as described [46] with U6 snRNA as the internal reference. The expression data of four target genes (Gh_Sca142710G01, Gh_Sca006071G01, Gh_A05G2828 and GhD07G0477) were extracted from high throughput RNA-seq for 10 DPA fibers of RNAi line and Coker 312 (unpublished data). All the primers used for miRNA qRT-PCR were listed in Supplementary S1 Table.

Results
High-throughput deep sequencing of 5-, 10-and 15-DPA fiber small RNA libraries from both Coker 312 and PHYA1 RNAi plants The libraries were subject to next generation sequencing on HiSeq 2500. Libraries had a minimum 8,052,792 and a maximum 12,947,430 short reads extracted from raw data ( Table 1). Adaptor sequences and low-quality reads composed of 39.54%~65.63% of raw reads were removed, and the average length of short reads after trim was from 22.5 to 24.2 nt (  Table 2). Most of these reads were clustered into rRNA, tRNA, snRNA, snoRNA and scRNA and removed before miRNA analyses. Among the rest of the annotated reads, 7,004, 7,291, 5,704, 8,698, 6,942, and 5,559 reads were matched to miRBase database ( Table 2). The size distribution of these small RNAs ranged from 15 to 30 nt, and the most abundant length of them was 24 nt in both Coker 312 and RNAi lines (S1 Fig).

Identification of differentially expressed known miRNAs in PHYA1 RNAi cotton fibers
Known and conserved miRNAs were identified by mapping to the miRBase database (release 21), and 77 known miRNAs belonging to 61 miRNA families were found and they are all  Fig 1A).  Among these conserved miRNAs, 39 miRNAs were up-regulated in PHYA1 RNAi compared to Coker 312, and 38 miRNAs were down-regulated in the RNAi line ( Fig 1A). Seven of conserved miRNAs, including miR7503, miR7514, miR399c, miR399d, miR160, miR169b and miR2950, were differentially expressed (fold change ! 2) in the PHYA1 RNAi line compared to Coker 312. It will be interesting to determine whether miRNA-mediated gene regulation could confer to the phenotype of the RNAi line with improved fiber quality.

Identification of novel miRNAs in PHYA1 RNAi cotton fibers
In addition to successful identification of 77 conserved miRNAs in the PHYA1 RNAi and wild type lines, novel miRNAs in fibers were also identified in both RNAi and wild-type cotton lines. After removing known miRNAs and other small RNAs, the stem-loop secondary structure of novel pre-miRNAs was predicted by MFOLD. The stability of an RNA secondary structure is affected by the minimal folding free energy (MFE). The lower the MFE value, the more stable the RNA structure. One of the criteria to annotate a novel miRNA is that it was generated by precise excision from the stem of a stem-loop precursor [47]. A total of 35 novel miR-NAs (named as novel-miR 1 to 35) were identified in the 18 fiber libraries ( Table 3). The sizes of these novel miRNAs ranged from 21 to 24 nt. The 24 nt miRNAs were the most abundant with 27 newly discovered miRNAs. Five, two and one novel miRNAs were found to have the sizes of 21, 22 and 23 nt, respectively. The stem-loop precursors of these novel miRNAs predicted by MFOLD have negative free folding energies from -45.7 to -230.1 kcalmol -1 , which are lower than those of tRNAs (-27.5 kcal mol -1 ) or rRNA (-33 kcal mol -1 ) ( Table 3). All miRNAs Ã of these novel miRNAs have been identified in RNA sequencing. Among these novel miRNAs, none were significantly differentially expressed in PHYA1 RNAi fibers ( 2 fold change) (Fig 1B). The 35 novel miRNAs were identified in cotton fibers for the first time, and their mature miRNA sequences were listed in Table 3.

Prediction and annotation of targets for novel miRNAs
A plant miRNA generally has perfect or near perfect complementarity to its target, and the cleavage happens between the 10 th and 11 th nucleotides from the 5´end of the miRNA. These characters were used to identify the targets of novel miRNAs [48]. By using psRNATarget, a total of 97 targets were predicted for the 35 novel miRNAs (Table 3). Among them, the targets of three novel miRNAs novel-miR14, 22, and 24 were not identified. The top-rated target genes for each of the novel miRNAs were listed in Table 3. The details of the miRNA silencing its target and the predicted cleavage/translation inhibition sites were also presented in Table 3. Among these targets, the gene encoding cytochrome P450-like TATA box binding protein (cytochrome P450-like TBP) has the highest hits, being the target of 9 different novel miRNAs which included novel-miR1, 2, 16, 19, 26, 27, 28, 31 and 32. Plant cytochrome P450s are involved in a wide range of biosynthetic reactions, including fatty acid conjugation, hormones synthesis, and generating defensive compounds. Of these 9 novel miRNAs, only novel-miR27 and novel-miR32 were accumulated slightly higher in the PHYA1RNAi line, and the other 7 miRNAs, novel-miR1, 2, 16, 19, 26, 28, and 31, were all slightly lower in RNAi cotton compared to Coker 312. The novel miRNAs were also predicted to silence an rRNA promoter binding protein and some transcription factors and related proteins, such as a bZIP transcription factor, a global transcription factor group isoform 1 and a WD40-like transcription factor. Genes for ribosomal proteins were also targeted by these novel miRNAs. In addition, many proteins related to translation, such as translation initiation factor and elongation factor 1α were predicted to be the targets of novel-miR2, 10, 26 and 27. The predicted targets of these novel miRNAs also include many enzymes such as ATP synthase, serine threonine-protein   kinase, peroxidase, acetolactate synthase, lipase, endo-β-glucanase, E3 ubiquitin protein ligase, receptor-like protein kinase, CoA reductase, xyloglucan glycosyltransferase, methylase, and some functional proteins, like chaperons, senescence-associated proteins, calcium-binding proteins and the 14-3-3-like protein ( Table 3).

Validation of differentially expressed miRNAs in PHYA1 RNAi line by qRT-PCR
The expression of twelve known miRNAs and seven novel miRNAs in 10-DPA fibers of PHYA1 RNAi and Coker 312 lines was also determined by qRT-PCR (Figs 2 and 3, S1 Table). The abundances of these miRNAs were relatively high in majority of the fiber libraries (Fig 1). The results showed that these 19 miRNAs were successfully detected by qRT-PCR and showed   consistency with the expression profiles analyzed by small RNA sequencing (Figs 2 and 3). Among the known miRNAs, miR162a, 396a/b, 2950 and 160 were up-regulated, and miR172, 399d, 167a/b, 390a/b/c, 164, 166b, 399c, and 169b were down-regulated (Fig 2). For the novel miRNAs, novel-miR22 and novel-miR8 were up-regulated, while novel-miR5, 3, 6, 2, and 4 were down-regulated (Fig 3). The expression level of these novel miRNAs (except for novel-miR4) in the PHYA1 RNAi line was less than two-fold change when compared to Coker 312.

Inverse expression patterns between four pairs of miRNAs and their targets.
Because miRNAs control the degradation of mRNA transcripts, the expression levels of miRNA and its target mRNA will be negatively correlated. To test this, the expression patterns of predicted targets of four novel miRNAs including novel-miR1, 2, 4 and 8 in 10 DPA fiber were analyzed by RNA sequencing and validated by qRT-PCR (Fig 4A). The results showed consistency between these two assays, and the relative expression levels of these miRNA targets were inversely correlated with the accumulation levels of corresponding miRNAs (Fig 4). For example, Gh_Sca142710G01 encoding an rRNA promoter binding protein was up-regulated in the PHYA1 RNAi plant (Fig 4A) compared to Coker 312, and its corresponding miRNA novel-miR1 was down-regulated ( Fig 4B). Novel-miR1 also targeted to cytochrome P450_TBP (The counts was too low to count in RNA Seq), which plays a crucial role in photosynthesis [49]. In addition to novel-miR1, eight of the novel miRNAs identified in upland cotton were predicted to target cytochrome P450 TBP (Table 3). This implied that the better quality fiber of the PHYA1 RNAi line may be partially conferred through the miRNA regulation of cytochrome P450s, which are involved in a wide range of biosynthetic reactions, including fatty acid conjugation, hormones synthesis, and generating defensive compounds [50]. Thus, miRNA regulation on cytochrome P450 might play a positive role during fiber development.

Discussion
Upland cotton (Gossypium hirsutum L.) produces the most commonly used textile fiber in the world. The improvement of cotton fiber quality using different genetics tools has long been a miRNAs in a cotton phytochrome A1 RNAi line key interest for cotton breeders. Using RNAi technology, Abdurakhmonov et al. [2] had developed PHYA1 RNAi cotton plants with improved fiber length, strength and fineness. Recently, the GhMYB25 genes in allotetraploid cotton genome have been successfully targeted for mutagenesis using the CRISPR/Cas9 technology [51]. This finding indicates that the CRISPR/Cas9 system will be an excellent tool for functional analysis of cotton genes in A and D subgenomes and improvement of cotton fiber quality. In this study, a genome-wide miRNAome analysis was performed to identify differentially expressed miRNAs, which might reveal the molecular mechanisms of the fiber quality improvement by PHYA1RNAi in cotton. A total of 77 known miRNAs were identified in cotton fibers. Of these, 7 known miRNAs were differentially expressed (fold change ! 2) in the PHYA1 RNAi plant. Using the published draft genome sequence of the allotetraploid cotton TM-1 line [30] as a reference, 35 novel cotton-specific miRNAs had been also identified and their target prediction was achieved.

Identification of novel miRNAs in cotton
The utilization of allotetraploid TM-1 cotton genome as a reference is very useful for miRNA identification and target prediction. Without the complete genome sequence of Gossypium hirsutum, it is difficult to conduct miRNAome analysis through bioinformatics prediction. First of all, small RNA read libraries generated by previous studies might contain errors [52]. Secondly, it is hard to evaluate the possibility of MIR gene loci duplication in allotetraploid. Thirdly, the isolation/identification of miRNAs in cotton is lagging when compared to other plant species due to the complexity of cotton genome compensation. The major criterion for novel miRNA prediction is the hairpin structure formed by precursor miRNAs. Due to incomplete upstream and downstream sequences, it is difficult to predict the secondary structure of miRNAs [37]. Furthermore, the available transcript libraries in cotton research were mostly from Gossypium raimondii, and thus target prediction had generated extensive redundancy. In this study, the precursors were predicted and identified for all novel miRNAs by using the psRNATarget, and the redundancy of the targets was reduced by the second blasting with TM-1 as the reference, because many of the hits from distinct miRNAs turned out to be the same gene.
Using high-throughput small RNA sequencing, 89 conserved and 8 novel miRNAs have been recently identified in cotton [37]. Based on the annotation of whole genome sequence of the cultivated Gossypium hirsutum TM-1 upland cotton, 602 miRNAs [29] and 301 miRNAs [30] were identified, respectively, in the allotetraploid cotton by two research groups. As a result, further identification and functional validation of cotton miRNAs are still needed, particularly, for cotton specific miRNAs and fiber-development-related miRNAs. In this study, a total of 35 novel miRNAs were identified in the small RNA libraries of both RNAi and control genotypes. None of these novel miRNAs were differentially expressed in RNAi lines or was specific to RNAi genotype although there were subtle difference in expression levels. This is desired objective of any RNAi crop development that positively contributes to its biosafety concerns because if novel non-existed signature is generated by RNAi hairpins incretion into a genome, as a first step alert, RNAi product must be subjected for further risk assessment analyses [3].
Of the novel miRNAs, 16 miRNAs were in the D-subgenome, and another 15 miRNAs were in the A-subgenome. Most of miRNAs identified previously were located in the Dgenome due to limited sequence information of upland cotton genomes [37]. An enlarged pool of sequenced genomes has increased the accuracy and coverage for bioinformatics prediction and analysis of miRNAs. The 24 nt miRNAs were the most abundant among the 35 novel miRNAs, which is consistent with previous studies [53,54]. Unlike the 21 nt miRNAs associated with AGO1, the 24 nt miRNAs are loaded onto AGO4, which has a preference for sRNAs with 5´terminal adenine [55,56]. As expected, most of the 24 nt miRNAs (9) had a 5´terminal adenine nucleotide, which is consistent with the previous study [54]. The predicted hairpin precursors had negative folding free energies from -45.7 to -230.1 kcal mol -1 , which were in agreement with previous studies [55]. By base pairing to mRNAs, microRNAs mediate mRNA degradation or translational repression. Most of the novel miRNAs in this study inhibit target expression through mRNA degradation, and a few of them regulate the expression of the target genes through translational repression.

Cotton fiber elongation-related miRNAs
In this study a total of 77 known miRNAs belonging to 61 miRNA families were identified in elongating fibers. Many of these conserved miRNAs were present at lower levels in RNAi cotton compared to Coker 312, including Gh-miR2950, 169b, 160, 399c and 399d. In contrast, Gh-miR7514, a cotton specific miRNA, was highly enriched in RNAi cotton. The target (ES793451) of Gh-miR2950, also a cotton specific miRNA, was predicted to encode a putative gibberellin 3 hydroxylase, which accumulated at higher levels in fibers [26]. Gibberellin 3 hydroxylase had been shown to control internode elongation in pea [57]. It was reported that Gh-miR2950 might affect fiber cell elongation via GA signaling [26]. The results in this study are consistent with this hypothesis, in that down-regulation of Gh-miR2950 promotes fiber elongation in the PHYA1 RNAi cotton by increasing Gibberellin 3 hydroxylase activity and therefore increasing the levels of the biologically active gibberellin GA1 [58]. miRNA169, 160, 399 are highly conserved miRNAs in plants. It was reported that miR160 and miR169 were significantly expressed at low levels in fibers than in seedlings [59]. miRNA160 had been shown to target three auxin response factors ARF10, 16, and 17 [26]. It is well known that auxin plays an essential role in cotton fiber development [10,11,60]. The results in this study suggested that miR160 can promote fiber development via the auxin signaling transduction pathway through increased expression of ARF10, 16 and 17. miR399 was predicted to target a MYB transcription factor [54]. Many MYBs, such as MYB25 and 109, were reported to promote fiber elongation and their suppression caused shorter fibers. The ABA 8´-hydroxylase catalyzed the first step of inactivation of ABA, and the active levels of ABA might increase by decreasing the expression of ABA 8´-hydroxylase in the RNAi line. In this study, miR7514 was accumulated at higher levels in RNAi cotton, and it targeted the gene for Rho guanylnucleotide exchange factor 7 (RhoGEF7) [61]. The functional analysis of RhoGEF7 in plants has not been documented so far. However, it was reported that ABA could induce the degradation of another Rho guanyl-nucleotide exchange factor RhoGEF1 [62]. Taken together, the results in this study imply that miRNA mediated gene regulation may be involved in fiber improvement of the PHYA1 RNAi cotton. Although the seven known miRNAs were differentially regulated in fiber cells and possibly involved in fiber development, some of them were found to play roles in biotic and abiotic stress responses. In wheat, miR160 was regulated under drought stress [63], miR164 was involved in salt stress response [64], and miR169b was differentially expressed in response to fungal infection [65]. In addition, miR396 in opium poppy was involved in regulation of secondary metabolite synthesis [66]. Another miRNA miRRNVL5 in G. hirsutum was found to play a role in regulation of cotton response to salt stress [67].
Cytochrome P450 TBP has the most abundant hits when the genome was screened with psRNATarget, and it was the target of 9 novel miRNAs identified in this study. Of these 9 miR-NAs, 7 of them had lower levels in RNAi cotton compared to Coker 312. Cytochrome P450 is involved in the biosynthesis of plant hormones (such as ABA, GA, BR) and secondary metabolites (e.g. phenylpropanoids, alkaloids, terpreoids) [68]. Cytochrome P450 enzymes also play critical roles in response to different abiotic and biotic stresses. The cytochrome P450 involved pathways were abundantly present in RNAi cotton, and miR172, a stress related miRNA was reported to target Cytochrome P450 in cotton and corn [69,70]. Taken together, although the expression of these novel miRNAs were not significantly changed in developing fibers, the observation of cytochrome P450 TBP as the target of 9 miRNAs suggests that miRNA mediated gene regulation on cytochrome P450 associated pathways might play a role in improving fiber quality in the PHYA1 RNAi cotton.

The possible functions of cytochrome P450 in cotton improvement
One possible explanation to the highest hits of cytochrome P450 TBP in miRNA target prediction is the high constitution of cytochrome P450 gene (CYP) family in plant genomes. It was reported that CYP genes occupied approximately 1% of the protein coding genes in six sequenced angiosperms, which include grape (Vitis Vinifera), papaya (Carica papaya), poplar (Populus trichocarpa), rice (Oryza sativa), Arabidopsis (Arabidopsis thaliana) and moss (Physcomitrella patens) and were involved in almost every aspects of the plant life [50]. Based on evolutionary analysis, the oldest CYPs are involved in biosynthesis of carotenoid and sterol. The middle branches of CYPs mediate the adaptation to environment, including abiotic stresses and biotic defenses. The newest CYPs are involved in the synthesis of structural components such as lignin, fatty acid, pigments, odorants, and signaling compounds [50]. Besides with improved fiber quality, PHYA1 RNAi plants have vigorous root and vegetative growth, and early flowering due to increased photosynthesis [2,3]. CYPs are light-driven oxidase enzymes in electron transfer chains [49]. Thus this study support the hypothesis that increased photosynthesis of RNAi cotton might be due to miRNA-mediated regulation on the expression of CYPs. In addition, CYPs are also involved in biosynthesis and catabolism of many plant hormones, such as ABA (CYP707), GA (CYP714) and BR (CYP724) [50]. ABA regulates osmotic stress tolerance [71], and GA plays a role in salinity adaptation [72]. As mentioned previously, CYPs also play crucial roles in environmental adaptation, which might explain the better tolerance of RNAi cotton to salinity, drought, and heat stresses.

Conclusion
Eighteen RNA libraries constructed using small RNAs from 5-, 10-, and 15-DPA fibers of PHYA1 RNAi and Coker 312 lines were sequenced using the Illumina HiSeq system. Sixty-one conserved miRNA families and thirty-five novel miRNAs were identified in the upland cotton lines. The targets of 6 conserved miRNAs, which expressed differentially in the RNAi line, were reported to participate in primary cell wall synthesis and phytohormone signaling pathways. The 35 novel miRNAs were identified in cotton for the first time, and their target genes were predicted. Nine novel miRNAs were identified to target cytochrome P450 TBP. Together, the results imply that miRNAs involved in fine-tune gene regulation might confer to the phenotype of the RNAi line with improved fiber quality.