Skip to main content
Download PDF

Identification and profiling of upland cotton microRNAs at fiber initiation stage under exogenous IAA application

BMC Genomics volume 20, Article number: 421 (2019) Cite this article

Abstract

Background

Cotton is the most essential textile crop worldwide, and phytohormones are critical for cotton fiber development. One example is the role of auxin in fiber initiation, but we know little molecular basis. MicroRNAs (miRNAs) have a significant function in cotton development; nevertheless their role in fiber initiation remains unclear. Here, exogenous IAA was applied to cotton plant before anthesis. Utilizing small RNA sequencing, the mechanism underlying miRNA-mediated regulation of fiber initiation under exogenous IAA treatment was investigated.

Results

With exogenous IAA application, the endogenous IAA and GA contents of IAA treated (IT) ovules were higher than control (CK) ovules at the fiber initiation stage, while endogenous ABA content was lower in IT than CK. Using scanning electron microscopy, we found the fiber number and size were significantly promoted in IT at 0 DPA. Fiber quality analysis showed that fiber length, uniformity, strength, elongation, and micronaire of IT were higher than CK, though not statistically significant, while lint percent was significantly higher in IT. We generated six small RNA libraries using − 3, 0, and 3 DPA ovules of IT and CK, and identified 58 known miRNAs and 83 novel miRNAs together with the target genes. The differential expressed miRNAs number between IT and CK at − 3, 0, 3 DPA was 34, 16 and 24, respectively. Gene ontology and KEGG pathway enrichment analyses for the target genes of the miRNAs expressed in a differential manner showed that they were significantly enriched in 30 terms and 8 pathways. QRT-PCR for those identified miRNAs and the target genes related to phytohormones and fiber development was performed, and results suggested a potential role of these miRNAs in fiber initiation.

Conclusions

The exogenous IAA application affected the relative phytohormone contents in ovule and promoted fiber initiation in cotton. Identification and profiling of miRNAs and their targets at the fiber initiation stage provided insights for miRNAs’ regulation function of fiber initiation. These findings not only shed light on the regulatory network of fiber growth but also offer clues for cotton fiber amelioration strategies in cotton.

Background

Cotton (Gossypium hirsutum L.) has been listed among the foremost economic crops in the world, which offers the greatest quantity of renewable resources for the cottonoracy. Cotton fibers refer to extremely elongated single-cell trichomes initiated from ovule epidermis. Development of cotton fiber has four differing stages with overlapps: initiation, elongation, secondary wall thickening and maturation [1, 2]. Typically, fiber initiation begins at the day of anthesis (0 DPA) and sustains for ~ 5 days. The number of fiber initials determines the number of fibers, and the initiation stage is thus the determining stage contributing to the final yield of cotton [3]. Fibers start elongation immediately after initiation and continue for ~ 25 days. Fibers enter the secondary wall thickening stage at about 20 DPA, and maturation stage at 35 DPA, respectively [4,5,6].

Phytohormones are quite essential for fiber development in cotton. It has been shown that gibberellins (GAs) not only have a positive effect on the fiber initiation and elongation, but also contribute to fiber secondary wall deposition [7, 8]. Ethylene facilitates fiber elongation, probably because of the enhance expression level of key genes in relation with sucrose synthase, tubulin, and expansin. In contrast, fiber growth was found to be inhibited by abscisic acid (ABA) in cotton ovule culture [9]. Indole-3-acetic acid (IAA), is a widespread natural auxin, that is critical for fiber development [10, 11]. Previous studies reported that auxin could promote fiber initiation and fiber units when applied to ovules exogenously [11, 12]. Pioneering in vivo experiments revealed that overexpression of a gene responsible for the auxin biosynthesis could simultaneously improve the yield and quality of cotton fiber [13]. In an previous experiment, with exogenous IAA application under field conditions, fiber number per ovule significantly increased [14]. Despite these findings, the complex mechanism by which auxin regulates fiber development requires further study.

MicroRNAs (miRNAs) are short noncoding RNAs with length ~ 22 nt. They negatively regulate the expression or translation of their target genes [15, 16]. MiRNAs are usually derived from intergenic regions as well as introns of genomes, and they are initially transcribed into primary miRNAs (pri-miRNAs) mainly by RNA polymerase II [17]. Pri-miRNAs are successively processed twice by a protein complex, including Dicer-like 1 (DCL1), into pre-miRNAs (miRNA precursors) with a typical hairpin structure, and then into mature miRNAs [18]. Mature miRNAs are loaded onto RNA-induced silencing complex (RISC) including AGO1 protein, and complementary binding to their target mRNAs for cleavage or translation inhibition [19, 20]. In the RNA-induced silencing process, miRNA has a guidance role through base pairing with the target mRNAs, while AGO proteins undertake an effectors role through recruitment of factors inducing translational repression, mRNA deadenylation and mRNA decay [21]. Extensive research have proved the miRNAs have a basic regulatory function of miRNAs in diver biological processes, including plant organ development, hormone signal transduction, stress tolerance and disease resistance [20, 22,23,24,25,26,27].

Recent studies of small RNAs in cotton have identified and characterized numerous miRNAs. In a study investigating the regulation of cotton plant height, An et al. sequenced small RNAs, identified 226 conserved miRNAs belonged to 32 known miRNA families along with 38 novel miRNAs, and found that several miRNAs were regulated in different height-related mutants [28]. As for cotton fiber, Pang et al. sequenced small RNAs from cotton fiber and other tissues and made identification of 27 conserved, four novel miRNA families [29]. In an investigation of fiber elongation, by sequencing the developing fiber of short fiber mutants Li1 and Li2, 24 conserved and 147 novel miRNA families along their targets were revealed, and four miRNAs were found to be negatively correlated with fiber lengths [30]. In Gossypium arboreum, Liu et al. found 48 known and 16 novel miRNAs, which expressed in elongating fibers in a differential manner [31]. Using a similar small RNA sequencing strategy, the present study focuses on the key role of exogenous IAA in fiber initiation, and miRNA-mediated regulation of fiber initiation. To date, there have not been any systematic studies using miRNA expression profiling and target genes analysis upon application of exogenous IAA during fiber initiation.

Here, we sequenced six libraries from − 3, 0, and 3 DPA ovules of control (CK) and IAA treated (IT) cotton at fiber initiation stage, and bioinformatics analyses were carried out for discovering of known and novel miRNAs. Based upon the expression profiles of the selected miRNAs together with the relevant target genes, we showed potential roles of some miRNAs in fiber initiation, and provide data for a deeper understanding of the molecular mechanism of miRNA in fiber initiation.

Methods

Plant materials

The genetic-standard upland cotton TM-1 was kindly provided by USDA-ARS, College Station, Texas, USA. This material was grown at the farm of Zhejiang University, China. At the day of the day of anthesis (0 DPA), flowers were labeled. The concentration of exogenous IAA was 5*10− 6 mol/L, and then this IAA solution was sprayed onto the cotton plant for three times every other day. After 3 days, bolls at − 3, − 1, 0, 1, 2, 3 and 5 DPA were harvested, and each of them were separated into three biological repeats randomly. Plants without exogenous IAA treatment were used as control group. After careful dissection from every boll, the ovules were frozen in liquid nitrogen immediately and stored at − 80 °C. The following six samples were used for the small RNA sequencing. C_1, C_2 and C_3 were the ovules of − 3, 0 and 3 DPA TM-1 without exogenous IAA. I_1, I_2 and I_3 were the ovules of − 3, 0, 3 DPA TM-1 with exogenous IAA treatment.

Phytohormone determination

Contents of GA, IAA, and ABA were measured using high performance liquid chromatography (HPLC). Standard samples were dissolved in the mobile phase containing methanol/ 0.2% phosphoric acid (1/1, v/v) to a constant volume of 1 ml.

All samples were weighed (1 g) and grounded into powder with liquid nitrogen. Samples were suspended into 4 ml acetonitrile and leached at 4 °C for 24 h. The supernatant was obtained by centrifugation with the speed of 1000 r/m for 10 min, and the sediment was washed thrice with acetonitrile. All the supernatants should be combined and dried by nitrogen at 35 °C. 3 ml volume of phosphoric acid buffer was added, the solution was centrifuged at 1000 r/m for 15 min, and the pH of the filtrate was adjusted to 3. An equal volume of ethyl acetate was used for extraction twice. The obtained organic phase was dried by nitrogen at 35 °C and added to the mobile phase to constant volume of 1 ml. 0.45 syringe filter (Agela, Newark, USA) was used for solution filtering.

Agilent was used for HPLC analysis of phytohormone, and the method was as previous reported by Cheng et al. [32]. Samples were measured in triplicates.

Small RNA libraries construction and sequencing

Extraction of total RNA was conducted out of ovules by utilizing a Total RNA Extraction kit (Aidlab, Beijing, China) in line with the manufacturer’s protocol. 1% agarose gels was performed for inspection of the RNA contamination and degradation. NanoPhotometer® spectrophotometer (IMPLEN, CA, USA) was used to examine the RNA purity. Qubit® RNA Assay Kit on Qubit® 2.0 Fluorometer (Life Technologies, CA, USA) was utilized for detecting the RNA concentration. Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA) to test the RNA integrity.

For every sample, 3 μg total RNA was applied for the construction of Small RNA library. With the manufacturer’s instruction, the NEBNext® Multiplex Small RNA Library Prep Set for Illumina® (NEB, USA.) was employed for generating of the sequencing library. Lastly, Illumina HiSeq 2500/2000 platform was used to sequence the library, thereby generating 50 bp single-end reads.

Analysis of small RNA sequencing data

Fastq format raw reads had been initially processed with application of custom Perl and Python scripts. During such procedure, clean reads would be acquired through removal of the following reads: reads containing poly-N, with 5′ adapter contaminants, without 3′ adapter or the insert tag; reads containing poly A, T, G, or C; and low quality reads from raw data. Meanwhile, calculation of Q20, Q30, and GC-content of the raw data was performed. Clean reads 18 to 30 nt in length were selected. For removal of tags that originated out of protein-coding genes, repeat sequences, rRNA, tRNA, snRNA, and snoRNA, the small RNA tags were mapped to RepeatMasker, Rfam database, or such kinds of data from the certain species itself. The small RNA tags were mapped onto the reference sequence using Bowtie [33] without mismatches for analysis expression and distribution on the reference sequence.

Identification of known miRNAs and prediction of novel miRNAs

To search the known miRNAs among the mapped small RNA sequences, miRBase21.0 was applied as the reference, utilization of mirdeep2 [34] and srna-tools-cli software was performed for making a prediction of the secondary structure of the flanking sequences of small RNAs matching known miRNAs. The first-position base of the identified miRNAs with specific length, bases upon every position of all identified miRNAs, and miRNA counts were calculated by custom scripts. Since novel miRNAs were predictable by the hairpin structure of miRNA precursors, miREvo [35] and mirdeep2 software were utilized together for novel miRNAs prediction based on the dicer cleavage site, the secondary structure, and the minimum free energy of the small RNA sequences, which was unannotated in the known miRNAs identification steps.

Prediction of miRNA target genes

MiRNA target gene prediction was carried out by psRobot_tar in psRobot [36].

Analysis of differentially expressed miRNAs

TPM was utilized for estimating the miRNA expression level using the formula herein under [37]: TPM = (mapped readcount/total reads)*1,000,000. Differential expression (log2 (fold-change) > 1) of miRNAs in the two samples should be assessed by DEGseq (2010) R package. Adjustment of P-value should be done with application of qvalue (< 0.01) [38].

GO and KEGG pathway enrichment analyses

Gene ontology (GO) enrichment analysis was employed for surveying the target genes of miRNAs expressed in a differential manner. GOseq, which was based Wallenius non-central hyper-geometric distribution [39], was carried out to perform a GO enrichment analysis.

KEGG [40] serves as a biological database, which provides the information of genomic, chemical and systemic (http://www.genome.jp/kegg/). It is for investigating biological system functions and utilities [41]. KOBAS software would be employed for testing for statistically significant enrichment of the target genes among the KEGG pathways.

Expression analysis of selected miRNAs and their target genes using real-time qPCR

Extraction of RNA samples was performed by Total RNA Extraction kit (Aidlab, Beijing, China). The reaction of reverse transcription for miRNAs was carried out with Mir-XTM miRNA First-Strand Synthesis kit (TaKaRa, Dalian, China). First strand cDNA synthesis of target genes should be conducted through utilization of PrimeScript™ First-Strand cDNA Synthesis Kit (TaKaRa, Dalian, China). Design of primers was made for qRT-PCR through Primer 5.0 software. All primers have been listed in Additional file 1: Table S1. The qRT-PCR reactions of miRNAs were carried out on the Lightcycler 96 system (Roche) with TB Green Premix Ex Taq (TaKaRa). Each reaction was carried out in triplicate. U6 and UBQ7 were employed to be internal control miRNA and internal control gene, respectively. Calculation of the relative expression level was conducted with utilization of the 2-ΔΔCt method.

Results

Effect of exogenous IAA on the content of endogenous phytohormones

To investigate the effect of exogenous IAA on cotton fiber and phytohormones, we sprayed IAA on cotton plants at − 9, − 7, and − 5 DPA. The analysis revealed significant differences in the endogenous phytohormone contents between CK and IT ovules (Fig. 1). At − 3 and − 1 DPA, there were no significant differences in endogenous IAA content of CK and IT. However, the endogenous IAA contents of IT significantly greater than that of CK at 0, 1, and 3 DPA, along with markedly significant difference at 2 and 5 DPA, indicating that exogenous IAA promoted the secretion of endogenous IAA. For endogenous GA content, the exogenous IAA application generally promoted the endogenous GA content at the fiber initiation stage. The endogenous GA content in IT was significantly greater than in CK at 3 DPA, while the contents in IT were higher than in CK at 0, 1, 2, and 5 DPA, but not significantly. It was suggested that exogenous IAA might have a positive impact on endogenous GA. For endogenous ABA content, it is generally depressed by the exogenous IAA application. The endogenous ABA contents in IT was significantly lower than in CK at − 1 and 0 DPA, but there were no significant differences between IT and CK at other time point. Thus, exogenous IAA might have a negative impact on endogenous ABA at the early stage of initiation.

Fig. 1
figure 1

The contents of endogenous phytohormones including IAA, GA and ABA in control (CK) and IAA treated (IT) ovules. * indicates P < 0.05

Full size image

Effect of exogenous IAA on fiber initiation and fiber quality

To further investigate the impact of exogenous IAA on the cotton fiber, we observed the development of fiber cells by scanning electron microscopy (SEM) of ovule surface at 0 DPA. As shown in Fig. 2, the number of fiber initials on IT ovules was significantly greater than on CK ovule. According to the statistical result, lint fiber number on CK and IT ovules were 4733 and 3867 per mm2, respectively, and lint fiber number on IAA ovules was significantly higher than that on IT ovules. Moreover, the sizes of lint fibers on IT ovules were also significantly larger than those on CK ovules. Above all, it is suggested that exogenous IAA facilitated fiber initiation and early development at the fiber initiation stage.

Fig. 2
figure 2

Differences in the initial development of fibers at 0 DPA between control (CK) and IAA treated (IT) plant. a and b show CK ovules surfaces at 80× and 2000× magnification; (c) and (d) show IT ovules surfaces at 80× and 2000× magnification. b and d are the black frame aeras of (a) and C. It was indicated that the number of lint fiber on ovules of IT was much higher than that of CK, and the size of of lint fiber on ovules of IT was much larger than that of CK

Full size image

Mature fiber quality traits including fiber length (FL), fiber uniformity (FU),, elongation (FE), strength (FS) and micronaire (MIC) were also measured. As shown in Table 1, FL, FU, FE and FS values tended to be higher in IT compared to CK, but the differences were not significant. MIC and lint percent (LP) in CK were 3.63 and 32.41%, respectively, while MIC and LP in IT were 3.90 and 34.07%, which was significantly higher than those in CK. These results indicated that exogenous IAA had an impact on mature fiber quality and positive effect on lint percent.

Table 1 Mature fiber quality trait in control (CK) and IAA treated (IT) plants (LOG 298 bytes)
Full size table

Deep sequencing of small RNA libraries from developing ovules

To identify miRNAs involved in cotton fiber initiation under exogenous IAA treatment, conduction of six small RNA libraries was performed out of IT and CK ovules of − 3, 0 and 3 DPA. A total of 70.98 million raw reads were generated, and after quality control and length selection of the reads, 9,655,792, 10,836,541, 10,112,797, 11,174,398, 10,823,945 and 9,611,924 small RNAs were obtained from C_1, C_2, C_3, I_1, I_2 and I_3, respectively (Table 2). On average, 76.75% of the reads were successfully aligned to the AD genome of G. hirsutum following the removal poor quality reads and adapter sequences. Furthermore, the sRNAs were categorized to be a couple of major classes. Those unique reads were very similar among the six libraries and ranged from 18 to 30 nt after the filtering process, the length of most sRNAs was 20-24 nt, with the 24 nt reads being the most abundant (Fig. 3).

Table 2 Classification statistics of small RNA sequence reads in six small RNA libraries
Full size table
Fig. 3
figure 3

Sequence length distribution of small RNAs from six libraries. Libraries of C_1, C_2, C_3, I_1, I_2, and I_3 were constructed from control and IAA treated ovules at − 3, 0, and 3 DPA. The 24 nt small RNA group is much larger than the others

Full size image

Identification of known miRNAs and prediction of novel miRNAs

For discovering known miRNAs in the six libraries, the mapped small RNA tags were aligned to all mature miRNA and miRNA precursor sequences in miRBase21.0. For CK and IT, 1974 and 2061 specific sequences had been identified as miRNA homologies in CK and IT plants, respectively. Following the removal of repeated sequences, 58 annotated known miRNAs under 50 miRNA families had been identified (Additional file 2: Table S2). Specifically, 53, 53, 52, 51, 52, and 54 of these miRNAs were distributed in C_1, C_2, C_3, I_1, I_2 and I_3, respectively. Three miRNAs, ghr-miR3476-3p, ghr-miR7488, and ghr-miR7497, were found exclusively in IT plants, and ghr-miR7498 was only present in CK plants.

The sequences failing to match known miRNAs were further compared with the cotton genome sequence to identify potential novel miRNA candidates using miREvo and mirdeep2 software [35]. Upon the basic of the characteristic hairpin structure of miRNA precursors, 84 novel miRNAs were identified among all six libraries; the sequences and stem-loop structures are included in Additional file 3: Table S3 and Additional file 4. The lengths of the novel miRNAs ranged from 18 to 24 bp, with 21 nt miRNAs comprising the greatest category, with subsequent by 24 nt miRNAs. First nucleotide base analysis was performed for these novel miRNAs and the results are summarized in Table 3. In the analysis, adenine (A, 36.14%) and uracil (U, 53.01%) were the most dominant nucleotides at the first position.

Table 3 The firstbase and lengths of the 83 novel cotton miRNAs
Full size table

Targets of known and novel miRNAs

Plant miRNAs bind to their targets with perfect or nearly perfect complementarity to direct post-transcriptional repression via mRNA cleavage or translation repressing. Therefore, miRNAs target genes are predictable upon the basic of the principle of high miRNA-target complementarity. psRobot was used to predict putative target genes, and 911 targets were identified for 109 miRNAs, with 437 target genes for 47 known miRNAs and 614 target genes for 63 novel miRNAs (Additional file 5: Table S4). Most miRNAs processed diverse target genes, such as Gh_A01G1274, Gh_A01G1281, and Gh_A03G0632 targeted by ghr-miR156a. Such miRNAs might therefore play a key role in cotton fiber initiation. From another point of view, there were also examples of a single gene targeted by several miRNAs, such as Gh_A03G0276 targeted by ghr-miR166b, ghr-miR1, and ghr-miR3.

In other cases, miRNAs such as ghr-miR7486a and ghr-miR15 had no predicted targets, suggesting that they might regulate gene expression through distinct mechanisms, such as translation inhibition.

Expression of known and novel miRNAs

Correlations among the six libraries were analyzed using Pearson’s correlation, and the results are shown in Fig. 4. The correlation coefficient between adjacent phases was higher than spaced ones both in the CK group (C_1, C_2 and C_3) and IT group (I_1, I_2 and I_3), and within groups as well. On the same day, the correlation coefficient within the same group was higher than that in different group. These findings revealed the dynamic trend of miRNA within the same group and the effect of IAA on miRNA alteration.

Fig. 4
figure 4

Pearson correlation of miRNA expression levels among six samples. Ovule samples of control and IAA treated at − 3, 0, and 3 DPA were named as C_1, C_2, C_3, I_1, I_2, and I_3. Correlation coefficients closer to a value of 1 indicate greater similarity in the expression patterns of two samples

Full size image

For discovering the differential expressed miRNAs between CK and IT, we detected the abundance of miRNA reads, and carried the veen diagram of the number of the differential expressed miRNAs and the heatmaps of the differential expressed miRNAs at − 3, 0 and 3 DPA (Fig. 5). As shown in Fig. 5a, altogether 56 miRNAs were expressed in a differential manner between CK and IT. At − 3 DPA, 34 differentially expressed miRNAs were identified, containing 11 known miRNAs and 23 novel miRNAs, among those 11 miRNAs were up-regulated and 23 miRNAs were down-regulated in IT (Fig. 5b). At 0 DPA, 16 miRNAs expressed in a differential manner, containing 10 known miRNAs and 6 novel miRNAs, and 8 miRNAs were up-regulated and eight miRNAs were down-regulated in IT (Fig. 5c). At 3 DPA, 24 miRNAs expressed in a differential manner, containing 10 known miRNAs and 14 novel miRNAs, and 14 miRNAs were up-regulated and 10 miRNAs were down-regulated in IT (Fig. 5d). Among the differentially expressed miRNAs, 58.11% belong to novel miRNAs, some of which have an essential function in the mechanism of IAA-regulated fiber initiation regulated by IAA.

Fig. 5
figure 5

Differential expressed miRNAs at − 3, 0 and 3 DPA. a showed the venn diagram of the differential expressed miRNAs between control (CK) and IAA treated (IT) at − 3, 0 and 3 DPA (C_1 VS I_1, C_2 VS I_3 and C_3 VS I_3), respectively. Three circles represented three comparative combinations. Sum of the numbers in the circles was the number of the differential expressed miRNAs in the circles corresponding comparative combinations. Sum of the numbers in the overlapping part of circles was the number of the common differential expressed miRNAs in different combinations. b, c and d showed heatmap of different expressed miRNAs between CK and IT at − 3, 0 and 3 DPA. The color indicates the fold-change as log2, from high (red) to low (blue), as indicated in the color scale. On the right are the names of the miRNAs

Full size image

According to the transcript per million (TPM) of differentially expressed miRNAs among six libraries, we merged the differential expressed miRNAs in all compared pairs of CK and IAA at − 3, 0 and 3DPA together to perform cluster analysis. As shown in Fig. 6, the miRNA expression patterns clustered into four classes for both known and novel miRNAs using the hierarchical clustering method. Class A comprised nine miRNAs and the expression of ghr-miR78, ghr-miR1, and ghr-miR5–2 decreased over time both in CK and IT. In class B, 29 miRNAs were clustered. Many miRNAs, including ghr-miR83, ghr-miR55, ghr-miR8, and ghr-miR7512, were down-regulated in IT. In class C, 42 miRNAs were clustered, and the expression levels of many miRNAs increased over time, like ghr-miR31, ghr-miR7505, ghr-miR49, and ghr-miR61. Class D comprised 24 miRNAs. The expression of several miRNAs decreased over time in CK and increased in IT, such as ghr-miR7502 and ghr-miR16. The potential involvement of these miRNAs in fiber initiation will require further research.

Fig. 6
figure 6

Complete hierarchical cluster analysis of differentially expressed miRNAs in cotton ovules. The analysis was performed by comparing the TPM of miRNAs among control and IAA treated (− 3, 0 and 3 DPA) libraries (C_1, C_2, C_3, I_1, I_2, and I_3). The color indicates the log10 (TPM + 1) change from high (red) to low (blue), as indicated in the color scale. MiRNA names are on the right, and the classes (a-d) they belong are shown on the left

Full size image

Functional annotation and classification

For further investigating the possible function of the identified miRNAs in regulating fiber initiation, GO and KEGG pathway analyses was carried out for all the target genes of miRNAs expressed in a differential manner at − 3, 0 and 3 DPA. For the target genes, there was significant enrichment of the following 30 terms, 12 biological processes, 13 molecular functions and five cellular component terms (Fig. 7a). The top three biological processes terms were intracellular signal transduction (GO: 0035556) containing 58 genes, signal transduction (GO: 0007165) containing 61 genes, and signaling (GO: 0023052) containing 61 genes. The most significant cellular component terms were nucleus (GO: 0005634) containing 52 genes, intracellular membrane-bounded organelle (GO: 0043231) containing 64 genes, and membrane-bounded organelle (GO: 0043227) containing 64 genes. The most significant molecular function terms were oxidoreductase activity, acting upon a sulfur group of donors, NAD(P) as acceptor (GO: 0016668), protein-disulfide reductase activity (GO: 0047134), and oxidoreductase activity, acting upon a sulfur group of donors (GO: 0016667), those which contained 51 genes.

Fig. 7
figure 7

GO and KEGG pathway enrichment analyses of total target genes of differentially expressed miRNAs at − 3, 0 and 3 DPA. a presents the GO enrichment analysis result (green: biological process; orange: cellular component; blue: molecular function), and * indicates P < 0.05. b presents the KEGG pathway analysis result, and the color indicates the qvalue change from 0 to 1 as indicated in the color scale and the size of the dot indicates the number of genes enriched in the corresponding pathway

Full size image

The KEGG pathway analysis revealed significant enrichment of the following: Plant-pathogen interaction, mRNA surveillance pathway, Glycerolipid metabolism, Circadian rhythm-plant, Plant hormone signal transduction, Endocytosis, and Spliceosome, Selenocompound metabolism (Fig. 7b). GO and KEGG pathway analyses results indicate that exogenous IAA might regulate fiber development through these particular pathways and processes.

Expression patterns of selected miRNAs and target genes

Selection and validation of several miRNAs and their target genes related to phytohormones and fiber development were performed by qRT-PCR to analyze their expression patterns. As shown in Fig. 8, the expression levels of all the selected miRNAs, ghr-miR393, ghr-miR10, ghr-miR156c, ghr-miR7504a, and ghr-miR36, were negatively correlated with their target genes. Three phytohormone-related genes were selected, TIR1, MYB33 and CRE1. The transcript levels of the auxin receptor gene TIR1 were significantly higher in IT than in CK at − 3, 0, and 3 DPA. MYB33, which encodes a transcription factor that undertakes the positive regulator of ABA, exhibited higher expression levels in IT verses CK at − 3 and 0 DPA and lower expression levels in IT versus CK at 3 DPA. CRE1, a cytokinin receptor gene, had higher expression levels in IT versus CK at − 3 DPA and lower expression level in IT at 0 and 3 DPA. Apart from these phytohormone biosynthesis genes, pre-mRNA splicing factor (PSF) had higher expression levels in IT than CK at − 3, 0, and 3 DPA. The expression patterns of two fiber-related genes, SPL9 and bZIP, were also investigated. SPL9, which can repress trichome formation, had higher expression levels in IT versus CK at − 3 DPA and lower expression levels in IT at 0 and 3 DPA. bZIP, a fiber development-related gene, had higher expression in IT than in CK at − 3, 0, and 3 DPA.

Fig. 8
figure 8

Expression patterns of selected miRNAs and their corresponding target genes. The columns indicated the abundance of miRNAs and target genes for control (CK) and IAA treated (IT) at − 3, 0, and 3 DPA from qRT-PCR results. The data represent the mean values ± SD of the three replicates

Full size image

Discussion

Cotton, as the most important natural fiber source for the textile industry, provides income for roughly 100 million families in 150 countries throughout the world [42]. Recently, it was shown that the major phytohormone IAA was a critical regulator during the development of cotton fiber both in vivo and in vitro [10,11,12,13]. Moreover, pre- and post-anthesis exogenous IAA treatment was found to contribute to fiber initiation [14]. The present study investigated the regulatory mechanism of the exogenous IAA effect on fiber initiation using small RNA sequencing and profiling.

Cotton fibers emerge from the ovule epidermis at 0 DPA [3, 43]. We observed that both the density and size of fibers on the ovule surface at 0 DPA were significantly increased in IT compared to CK. These findings indicated that exogenous IAA promoted fiber initiation, which is consistent with previous research [14]. Mature fiber quality traits, including FL, FU, FE, and FS, tended to be higher in IT than in CK, although the differences were not significant, the findings suggest that exogenous IAA may have a positive effect on fiber quality, which is particularly relevant for the changes in spinning technology [44]. The significantly higher MIC and LP in IT also suggested that exogenous IAA may contribute to cotton fiber development, and this result is consistent with previous studies [14].

It is well established that endogenous phytohormones critically impact fiber development. In our research, the endogenous IAA content of IT was higher than that of CK after exogenous IAA treatment. It was previously reported that IAA promoted fiber initiation [10, 13]. Therefore, the significantly greater number and size of lint fibers in IT versus CK may have been induced by the increased IAA content in IT. The endogenous GA content was greater than in CK from 0 DPA to 5 DPA, which implies that exogenous IAA promoted endogenous GA content. As previously reported, GA may affect fiber initiation and development [7, 14]. Thus, the increased GA content, which was induced by exogenous application, may have also induced the significantly greater number and size of lint fibers in IT. ABA content was lower in IT than in CK, which suggests that exogenous IAA negatively regulate ABA secretion. ABA decreases the positive effect of GA3 and IAA and has been shown to inhibit the growth of the ovule fiber [9, 10]. On the contrary, ABA played a negative role in fiber development. We therefore hypothesize that exogenous IAA might promote fiber initiation by depressing the ABA biosynthesis. Taken together, the present findings suggest that exogenous IAA affects cotton fiber development in part by regulating the phytohormones network in cotton plants, content of several phytohormones, including increasing the endogenous IAA and GA contents and reducing the endogenous ABA content.

Cotton miRNAs were first identified using a comparative genomic approach [45]. Since then, as deep sequencing technology has provided a powerful high-throughput strategy for miRNA identification, l Numarge numbers of cotton miRNAs have been characterized and implicated in plant growth and resistance to biotic and abiotic stresses.Numerous cotton miRNAs are found to be expressed in a differential manner under salinity, drought and heavy metal stress and had negative correlation to their target genes, indicating a function of these miRNAs in abiotic stress tolerance [46,47,48]. Moreover, miRNAs are discovered to take part in the response to biotic stresses, such as Verticillium dahliae and leafroll dwarf polerovirus in cotton [49,50,51]. The mechanism of cotton fiber development is complex and still unclear, involving numerous miRNAs with important roles. Several prior studies have identified hundreds of conserved and novel miRNAs related to fiber development [6, 29, 31, 42, 45, 52,53,54,55,56,57], the first of which were identified and predicted along with their targets by Zhang et al. in 2007. MiR396, miR414, and miR782 have been shown to target genes encoding callous synthase, fiber protein Fb23, and fiber quinine-oxidoreductase, which are closely linked to fiber development [45]. Using small RNA libraries generated from ovules of the fiberless mutant Xu-142w and its wild type control Xu-142 during cotton fiber initiation and development, Kwak et al. identified 22 conserved miRNA families, most of which were significantly differentially expressed such as miR156/157, miR165/166, miR172, miR395, miR397 and miR399. [57]. In 2014, Guan et al. reported that miR828 and miR858 regulated the MYB2 gene involved in cotton fiber development. MYB2 could induce trichomes in a no-trichome Arabidopsis mutant [58]. So far, no study has researched on the effect of exogenous IAA on the miRNA regulation of fiber initiation. In our present work, deep sequencing of small RNAs was performed for IT and CK cotton ovules and millions of small RNA reads were mapped to the cotton genome. In according to preceding reports in cotton [29] and other plants such as rice [59] and maize [60], 21 and 24 nt miRNAs represented the largest fraction of all the miRNAs. 24 nt miRNAs were reported to be loaded onto AGO4 and showed a preference for sRNAs with 5′ terminal adenine [61, 62]. In our results, adenine at the firstbase contributed the most of all the 24 nt miRNAs. Moreover, uridine was the most nucleotides of all the 21 nt miRNAs, which accorded to the earlier report that miRNAs with firstbase of uridine were usually harbored by AGO1 [62]. Totally, we identified 58 known and 83 novel miRNAs under exogenous IAA treatment during fiber development. At the first position of novel miRNAs, A and U were the prominent nucleotides, in accordance with previous research [31]. Furthermore, 34, 16, and 24 miRNAs were found to be differentially expressed at − 3, 0, and 3 DPA ovules, respectively. It is therefore possible that exogenous IAA promoted fiber initiation through the regulatory activity of these miRNAs on their target genes. MiRNAs with decreasing or increasing expression over time, as determined through cluster analysis result, might also be critical for fiber initiation. As a preliminary step to determine whether any of these miRNAs are indeed critical for cotton fiber development, we identified 357 target genes, which were predicted by the miRNAs expressed in a differential manner, and they were subjected to GO and KEGG analyses. A previously analysis by Xie et al. identified 820 genes from miRNAs related to fiber development, which were associated with 1027 GO terms and 78 KEGG pathways [56]. For the target genes identified in the present study, there was a significant enrichment for 31 GO terms and eight pathways, which might provide some clues about miRNA regulation of fiber development.

The process of cotton fiber development is very complicated, which is related to numerous factors. In the present research, several miRNAs with target genes involved in fiber development were selected to detect their expression patterns. TIR1, an auxin receptor, plays a role in mediating the degradation of auxin/IAA and transcription of auxin-regulated [63]. TIR1 transcript levels were higher in IT than in CK at − 3, 0, and 3 DPA, which indicated that exogenous IAA affected the auxin metabolism in the cotton plants. The finding further suggested that exogenous IAA promoted fiber initiation through the regulation of auxin metabolism in plants. MiR393, which regulates TIR1, was identified in the analysis, and its expression was lower in IT than in CK at − 3, 0, and 3 DPA, consistent with its negative regulatory effect on TIR1. As such, miR393 probably played a role in the regulating the development of cotton fiber.MYB33, a transcription gene that positively regulates ABA response, is a target gene of ghr-miR319 [64]. In our research, ghr-miR10 showed high homologies with ghr-miR319, which was also predicted to regulate the MYB33, exhibited lower expression levels in CK than in IT at − 3 and 0 DPA and higher expression levels at 3 DPA. The expression of MYB33 was the opposite. In accordance with a reported research, our results showed that the ABA content was synchronous with the expression of ghr-miR319 [65, 66]. However, the expression pattern of MYB33 was opposite to ABA content. We speculate that there might be a balance among miRna11 cleavage, MYB33 expression and ABA content, similar to that described in a previous report [65] and that these three factors all participate in the regulation of fiber initiation. CRE1, a cytokinin receptor, was identified previously [67]. In our research, CRE1 transcripts, targeted by ghr-miR156c, were more abundant in IT than in CK at − 3 DPA and less abundant at 0 and − 3 DPA, which was opposite to the expression pattern of ghr-miR156c. This result was consistent with a previous report that cytokinins stimulate fiber initiations before flowering and inhibit fiber growth after flowering [68]. The finding also suggests that the positive effect of exogenous IAA on fiber initiation may involve the regulation of cytokinins metabolism. Ghr-miR156, one of the most conserved miRNAs, also targets SPL9, which controls trichome distribution in Arabidopsis [69], a process akin to cotton fiber development. In our analysis, SPL9 transcripts were lower at 0 and 3 DPA in IT, suggesting that the down regulated expression of SPL9 might promote the cotton fiber initiation. Another transcription factor of interest, bZIP, regulates fiber development in Gossypium barbadense L. [70]. Its expression was higher in IT at − 3 and 0 DPA and lower at 3 DPA, indicating that bZIP might have a positive effect on fiber initiation. An opposite expression pattern was found for the corresponding miRNA, ghr-miR36. We also investigated the expression pattern of the pre-mRNA splicing factor gene PSF [53], which is regulated by ghr-miR7504a. Its expression was higher at − 3, 0, and 3 DPA in IT, suggesting that increased mRNA maturation may be required for more biochemical processes under exogenous IAA treatment. Altogether, our analysis identified many candidate miRNAs potentially involved in cotton fiber initiation.

Conclusion

Cotton serves as a foremost fiber plant worldwide, and the development process of cotton fiber is complex. Here, we found that exogenous IAA application on flower buds before anthesis affected the IAA, GA and ABA contents in the ovule and increased fiber initiation and early development at the fiber initiation stage. Using small RNA sequencing, altogether 58 known miRNAs and 83 novel miRNAs with target genes were identified. The differential expression patterns of certain miRNAs and GO and KEGG analyses further illustrated that miRNAs were very important in the regulation of fiber initiation. These results lay the foundation for investigating the complex regulatory network of fiber development, which offer clues for fiber melioration using molecular strategies.

References

  1. Basra AS, Malik CP. Development of the cotton fiber. Int Rev Cytol. 1984;89(6):65–113.

    Article  CAS  Google Scholar 

  2. Qin YM, Zhu YX. How cotton fibers elongate: a tale of linear cell-growth mode. Curr Opin Plant Biol. 2011;14:106–11.

    Article  CAS  Google Scholar 

  3. Taliercio EW, Boykin D. Analysis of gene expression in cotton fiber initials. BMC Plant Biol. 2007;7:22.

    Article  Google Scholar 

  4. Yang Y, Bian S, Yao Y, Liu J. Comparative proteomic analysis provides new insights into the fiber elongating process in cotton. J Proteome Res. 2008;7:4623–37.

    Article  CAS  Google Scholar 

  5. Shi Y, Zhu S, Mao X, Feng J, Qin Y, Zhang L, et al. Transcriptome profiling, molecular biological, and physiological studies reveal a major role for ethylene in cotton fiber cell elongation. Plant Cell. 2006;18:651–64.

    Article  CAS  Google Scholar 

  6. Wang M, Sun R, Li C, Wang Q, Zhang B. MicroRNA expression profiles during cotton (Gossypium hirsutum L) fiber early development. Sci Rep-UK. 2017;7:44454.

    Article  CAS  Google Scholar 

  7. Xiao Y, Li D, Yin M, Li X, Zhang M, Wang Y, et al. Gibberellin 20-oxidase promotes initiation and elongation of cotton fibers by regulating gibberellin synthesis. J Plant Physiol. 2010;167:829–37.

    Article  CAS  Google Scholar 

  8. Bai W, Xiao Y, Zhao J, Song S, Hu L, Zeng J, et al. Gibberellin overproduction promotes sucrose synthase expression and secondary cell wall deposition in cotton fibers. PLoS One. 2014;9:e96537.

    Article  Google Scholar 

  9. Fan L, Meng L, Ni ZY, Hu WR, Wang J. Digital image analysis of expansion growth of cultured cotton ovules with fibers and their responses to ABA. In Vitro Cell Dev-Pl. 2011;47:369–74.

    Article  Google Scholar 

  10. Beasley CA. Hormonal regulation of growth in unfertilized cotton ovules. Science. 1973;179:1003–5.

    Article  CAS  Google Scholar 

  11. Gialvalis S, Seagull RW. Plant hormones alter fiber initiation in unfertilized, cultured ovules of Gossypium hirsutum. J Cotton Sci. 2001;5:252–8.

    CAS  Google Scholar 

  12. Beasley CA, Ting IP. The effects of plant growth substances on in vitro fiber development from fertilized cotton ovules. Am J Bot. 1973;60:130–9.

    Article  CAS  Google Scholar 

  13. Zhang M, Zheng X, Song S, Zeng Q, Hou L, Li D, et al. Spatiotemporal manipulation of auxin biosynthesis in cotton ovule epidermal cells enhances fiber yield and quality. Nat Biotechnol. 2011;29:453.

    Article  CAS  Google Scholar 

  14. Seagull RW, Giavalis S. Pre- and post-anthesis application of exogenous hormones alters fiber production in Gossypium hirsutum L. cultivar Maxxa GTO. J Cotton Sci. 2004;8:105–11.

    CAS  Google Scholar 

  15. Ambros V. The functions of animal microRNAs. Nature. 2004;431:350.

    Article  CAS  Google Scholar 

  16. Zhang B, Wang Q, Pan X. MicroRNAs and their regulatory roles in animals and plants. J Cell Physiol. 2007;210:279–89.

    Article  Google Scholar 

  17. Lee Y, Kim M, Han J, Yeom KH, Lee S, Baek SH, et al. MicroRNA genes are transcribed by RNA polymerase II. EMBO J. 2004;23(20):4051–60.

    Article  CAS  Google Scholar 

  18. Xie Z, Kasschau KD, Carrington JC. Negative feedback regulation of dicer-Like1 in Arabidopsis by microRNA-guided mRNA degradation. Curr Biol. 2003;13:784–9.

    Article  CAS  Google Scholar 

  19. Cesar L, Kasschau KD, Rector MA, Carrington JC. Endogenous and silencing. Plant Cell. 2002;14:1605–19.

    Article  Google Scholar 

  20. Chen X. A MicroRNA as a translational repressor of APETALA2 in Arabidopsis flower development. Science. 2004;303:2022.

    Article  CAS  Google Scholar 

  21. Minju H, V Narry K. Regulation of microRNA biogenesis. Nat Rev Mol Cell Biol. 2014;15:509–24.

    Article  Google Scholar 

  22. Yoon EK, Yang JH, Lim J, Kim SH, Kim SK, Lee WS. Auxin regulation of the microRNA390-dependent transacting small interfering RNA pathway in Arabidopsis lateral root development. Nucleic Acids Res. 2010;38:1382–91.

    Article  CAS  Google Scholar 

  23. Zhang B, Wang Q. MicroRNA-based biotechnology for plant improvement. J Cell Physiol. 2015;230:1–15.

    Article  Google Scholar 

  24. Navarro L, Jones JDG. A plant miRNA contributes to antibacterial resistance by repressing auxin signaling. Science. 2006;312:436.

    Article  CAS  Google Scholar 

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

    Article  CAS  Google Scholar 

  26. Jin W, Wu F. Characterization of miRNAs associated with Botrytis cinerea infection of tomato leaves. BMC Plant Biol. 2015;15(1). https://0-doi-org.brum.beds.ac.uk/10.1186/s12870-014-0410-4.

    Article  Google Scholar 

  27. Zhang H, Jin W, Zhu X, Liu L, He Z, Yang S, et al. Identification and characterization of Salvia miltiorrhizain miRNAs in response to replanting disease. PLoS One. 2016;11:e0159905.

    Article  Google Scholar 

  28. An W, Gong W, He S, Pan Z, Sun J, Du X. MicroRNA and mRNA expression profiling analysis revealed the regulation of plant height in Gossypium hirsutum. BMC Genomics. 2015;16:886.

    Article  Google Scholar 

  29. Pang M, Woodward AW, Agarwal V, Guan X, Ha M, Ramachandran V, et al. Genome-wide analysis reveals rapid and dynamic changes in miRNA and siRNA sequence and expression during ovule and fiber development in allotetraploid cotton (Gossypium hirsutum L.). Genome Biol. 2009;10:R122.

    Article  Google Scholar 

  30. Naoumkina M, Thyssen GN, Fang DD, Hinchliffe DJ, Florane CB, Jenkins JN. Small RNA sequencing and degradome analysis of developing fibers of short fiber mutants Ligon-lintles-1 (Li1) and −2 (Li2) revealed a role for miRNAs and their targets in cotton fiber elongation. BMC Genomics. 2016;17:360.

    Article  Google Scholar 

  31. Wang Y, Ding Y, Liu JY. Identification and profiling of microRNAs expressed in elongating cotton fibers using small RNA deep sequencing. Front Plant Sci. 2016;7:1722.

  32. Li C, Zhao T, Li C, Mei L, Yu E, Dong Y, et al. Determination of gossypol content in cottonseeds by near infrared spectroscopy based on Monte Carlo uninformative variable elimination and nonlinear calibration methods. Food Chem. 2016;221:990–6.

    Article  Google Scholar 

  33. Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10:R25.

    Article  Google Scholar 

  34. Friedländer MR, Mackowiak SD, Li N, Chen W, Rajewsky N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012;40:37–52.

    Article  Google Scholar 

  35. Wen M, Shen Y, Shi S, Tang T. miREvo. An integrative microRNA evolutionary analysis platform for next-generation sequencing experiments. BMC Bioinf. 2012;13:140.

    Article  CAS  Google Scholar 

  36. Wu H, Ma Y, Chen T, Wang M, Wang X. PsRobot: a web-based plant small RNA meta-analysis toolbox. Nucleic Acids Res. 2012;40:22–8.

    Article  CAS  Google Scholar 

  37. Zhou L, Chen J, Li Z, Li X, Hu X, Huang Y, et al. Integrated profiling of microRNAs and mRNAs: microRNAs located on Xq27.3 associate with clear cell renal cell carcinoma. PLoS One. 2010;5:e15224.

    Article  CAS  Google Scholar 

  38. Storey JD. The positive false discovery rate: a bayesian interpretation and the q-value. Ann Stat. 2003;31:2013–35.

    Article  Google Scholar 

  39. Young MD, Wakeeld MJ, Smyth GK, goseq OA. Gene ontology testing for RNA-seq datasets. Genome Biol. 2012;11:R14.

    Article  Google Scholar 

  40. Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, et al. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2008;36:480–4.

    Article  Google Scholar 

  41. Mao X, Cai T, Olyarchuk JG, Wei L. Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics. 2005;21:3787–93.

    Article  CAS  Google Scholar 

  42. Sun R, Li C, Zhang J, Li F, Ma L, Tan Y, et al. Differential expression of microRNAs during fiber development between fuzzless-lintless mutant and its wild-type allotetraploid cotton. Sci Rep-UK. 2017;7:3.

    Article  Google Scholar 

  43. Stewart JMD. Fiber initiation on the cotton ovule (Gossypium hirsutum). Am J Bot. 1975;62:723–30.

    Article  Google Scholar 

  44. Kohel RJ. Cotton germplasm resources and the potential for improved fiber production and quality. In: Basra AS, editor. Cotton fibers. New York: The Haworth Press Inc; 1999. p. p167–82.

    Google Scholar 

  45. Zhang B, Wang Q, Wang K, Pan X, Wang Y, Guo T, et al. Identification of cotton microRNAs and their targets. Gene. 2007;397:26–37.

    Article  CAS  Google Scholar 

  46. He Q, Zhu S, Zhang B. MicroRNA-target gene responses to lead-induced stress in cotton (Gossypium hirsutum L.). Funct Integr Genomic. 2014;14:507.

    Article  CAS  Google Scholar 

  47. Peng Z, He S, Gong W, Sun J, Pan Z, Xu F, et al. Comprehensive analysis of differentially expressed genes and transcriptional regulation induced by salt stress in two contrasting cotton genotypes. BMC Genomics. 2014;15:760.

    Article  Google Scholar 

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

    Article  CAS  Google Scholar 

  49. Zhu Q, Fan L, Liu Y, Xu H, Llewellyn D, Wilson I. miR482 regulation of NBS-LRR defense genes during fungal pathogen infection in cotton. PLoS One. 2013;8:e84390.

    Article  Google Scholar 

  50. Romanel E, Silva TF, Corrêa RL, Farinelli L, Hawkins JS, Schrago CE, et al. Global alteration of microRNAs and transposon-derived small RNAs in cotton (Gossypium hirsutum) during cotton leafroll dwarf polerovirus (CLRDV) infection. Plant Mol Biol. 2012;80:443–60.

    Article  CAS  Google Scholar 

  51. Yin Z, Yan L, 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.

    Article  CAS  Google Scholar 

  52. Miao Q, Deng P, Saha S, Jenkins JN, Hsu CY, Abdurakhmonov IY, et al. Genome-wide identification and characterization of microRNAs differentially expressed in fibers in a cotton phytochrome A1 RNAi line. PLoS One. 2017;12:e0179381.

    Article  Google Scholar 

  53. Wang Z, Xue W, Dong C, Jin L, Bian S, Wang C, et al. A comparative miRNAome analysis reveals seven fiber initiation-related and 36 novel miRNAs in developing cotton ovules. Mol Plant. 2012;5:889–900.

    Article  CAS  Google Scholar 

  54. Xue W, Wang Z, Du M, Liu Y, Liu J. Genome-wide analysis of small RNAs reveals eight fiber elongation-related and 257 novel microRNAs in elongating cotton fiber cells. BMC Genomics. 2013;14:629.

    Article  CAS  Google Scholar 

  55. Liu N, Tu L, Tang W, Gao W, Lindsey K, Zhang X. Small RNA and degradome profiling reveals a role for miRNAs and their targets in the developing fibers of Gossypium barbadense. Plant J. 2015;80:331–44.

    Article  Google Scholar 

  56. Xie F, Jones DC, Wang Q, Sun R, Zhang B. Small RNA sequencing identifies miRNA roles in ovule and fibre development. Plant Biotechnol J. 2015;13:355–69.

    Article  CAS  Google Scholar 

  57. Kwak P, Qin Q, Xu S, Cheng X, Zhi M. Enrichment of a set of microRNAs during the cotton fiber development. BMC Genomics. 2009;10:457.

    Article  Google Scholar 

  58. Guan X, Pang M, Nah G, Shi X, Ye W, Stelly DM, et al. miR828 and miR858 regulate homoeologous MYB2 gene functions in Arabidopsis trichome and cotton fibre development. Nat Commun. 2014;5:3050.

    Article  Google Scholar 

  59. Li T, Li H, Zhang Y, Liu J. Identification and analysis of seven H2O2-responsive miRNAs and 32 new miRNAs in the seedlings of rice (Oryza sativa L. ssp. indica). Nucleic Acids Res. 2011;39:2821–33.

    Article  CAS  Google Scholar 

  60. Kang M, Zhao Q, Zhu D, Yu J. Characterization of microRNAs expression during maize seed development. BMC Genomics. 2012;13:360.

    Article  CAS  Google Scholar 

  61. Wu L, Zhou H, Zhang Q, Zhang J, Ni F, Liu C, et al. DNA methylation mediated by a microRNA pathway. Mol Cell. 2010;38:465–75.

    Article  CAS  Google Scholar 

  62. Mi S, Cai T, Hu Y, Chen Y, Hodges E, Ni F, et al. Sorting of small RNAs into Arabidopsis Argonaute complexes is directed by the 5′ terminal nucleotide. Cell. 2008;133:116–27.

    Article  CAS  Google Scholar 

  63. Nihal D, Sunethra D, Mark E. The F-box protein TIR1 is an auxin receptor. Nature. 2005;435:446–51.

    Article  Google Scholar 

  64. Reyes JL, Chua NH. ABA induction of miR159 controls transcript levels of two MYB factors during Arabidopsis seed germination. Plant J. 2010;49:592–606.

    Article  Google Scholar 

  65. Thiebaut F, Rojas CA, Almeida KL, Grativol C, Domiciano GC, Lamb CR, et al. Regulation of miR319 during cold stress in sugarcane. Plant Cell Environ. 2012;35:502–12.

    Article  CAS  Google Scholar 

  66. Palatnik JF, Wollmann H, Schommer C, Schwab R, Boisbouvier J, Rodriguez R, et al. Sequence and expression differences underlie functional specialization of Arabidopsis microRNAs miR159 and miR319. Dev Cell. 2007;13:115.

    Article  CAS  Google Scholar 

  67. Inoue T, Higuchi M, Hashimoto Y, Seki M, Kobayashi M, Kato T, et al. Identification of CRE1 as a cytokinin receptor from Arabidopsis. Nature. 2001;409:1060–3.

    Article  CAS  Google Scholar 

  68. Chen J, Du X, Zhou X, Zhao H. Levels of cytokinins in the ovules of cotton mutants with altered fiber development. J Plant Growth Regul. 1997;16:181–5.

    Article  CAS  Google Scholar 

  69. Yu N, Cai W, Wang S, Shan C, Wang L, Chen X. Temporal control of trichome distribution by MicroRNA156-targeted SPL genes in Arabidopsis thaliana. Plant Cell. 2010;22(7):2322–35.

    Article  CAS  Google Scholar 

  70. Yuan D, Tu L, Zhang X. Generation, annotation and analysis of first large-scale expressed sequence tags from developing fiber of Gossypium barbadense L. PLoS One. 2011;6:e22758.

    Article  CAS  Google Scholar 

Download references

Acknowledgements

We are grateful to Dr. Zhang (Zhejiang University, China) for technical assistance and Dr. Muhammad Daud Khan (Kohat University of Science and Technology, Pakistan), an English native speaker familiar with the scientific terminology, for the help of correcting the English writing.

Funding

The research work was funded by The National Natural Science Foundation of China (31571715). The funders had no role in the design or execution of the study.

Availability of data and materials

All relevant data are within this article and its additional files.

Author information

Author notes

  1. Tianlun Zhao and Xiaojian Xu contributed equally to this work.

Authors and Affiliations

  1. Department of Agronomy, Zhejiang University, Zhejiang, 310058, Hangzhou, China

    Tianlun Zhao, Xiaojian Xu, Min Wang, Cheng Li, Cong Li, Rubing Zhao, Shuijin Zhu & Jinhong Chen

  2. Zhejiang Province Key Laboratory of Plant Secondary Metabolism and Regulation, Zhejiang Sci-Tech University, Zhejiang, 310018, Hangzhou, China

    Qiuling He

Authors
  1. Tianlun Zhao
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Xiaojian Xu
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Min Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Cheng Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Cong Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Rubing Zhao
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Shuijin Zhu
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Qiuling He
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Jinhong Chen
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

TZ and JC designed the experiments and wrote the manuscript. XX performed the experiments. QH and MW analyzed the data, CL1 (Cheng Li), CL2 (Cong Li), and RZ assisted in editing the article. SZ, QH and JC conducted and supervised the experiments. All authors have read and approved the manuscript.

Corresponding authors

Correspondence to Qiuling He or Jinhong Chen.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Table S1. All primers used in this study. (DOCX 14 kb)

Additional file 2:

Table S2. Names and sequences of the identified known miRNAs in six libraries. (XLSX 10 kb)

Additional file 3:

Table S3. Names and sequences of the identified novel miRNAs in six libraries. (XLSX 11 kb)

Additional file 4:

Structures of all the identified novel miRNAs. (ZIP 1005 kb)

Additional file 5:

Table S4. MiRNAs and the corresponding target genes. (XLSX 30 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhao, T., Xu, X., Wang, M. et al. Identification and profiling of upland cotton microRNAs at fiber initiation stage under exogenous IAA application. BMC Genomics 20, 421 (2019). https://0-doi-org.brum.beds.ac.uk/10.1186/s12864-019-5760-8

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://0-doi-org.brum.beds.ac.uk/10.1186/s12864-019-5760-8

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us