A comparative analysis of small RNAs between two Upland cotton backcross inbred lines with different fiber length:Expression and distribution

Guoyun Liu,Mn Wu,Wenfeng Pei,Xihu Li,Nuohn Wng,Jinjing M,Xinshn Zng,Shuxun Yu,*,Jinf Zhng*,Jiwen Yu,*

a State Key Laboratory of Cotton Biology,Cotton Institute of the Chinese Academy of Agricultural Sciences,Key Laboratory of Cotton Genetic Improvement,Ministry of Agriculture,Anyang 455000,Henan,China

b Department of Plant and Environmental Sciences,New Mexico State University,Las Cruces,NM 88003,USA

Keywords:Allopolyploid Gossypium Fiber elongation miRNA QTL hotspots

A B S T R A C T The cotton fiber is the most important raw material for the textile industry and an ideal model system for studying cell elongation.However,the genetic variation of fiber elongation in relation to miRNA is poorly understood.A high-throughput comparative RNA-seq of two lines differing in fiber length(FL)from a backcross inbred line(BIL)population of G.hirsutum×G.barbadense revealed differentially expressed(DE)miRNAs and their targets in rapidly elongating fibers.A real-time quantitative PCR analysis was further performed to validate the results.A total of 463(including 47 DE)miRNAs were identified,and seven DE miRNAs were co-localized with seven FL quantitative trait loci(QTL)identified in the G.hirsutum×G.barbadense population.Of 82(including 21 DE)targets identified,nine(including one DE)were also co-localized with the seven FL QTL.The relationship between the allopolyploid and its diploid ancestral species with respect to miRNAs and their targets was also characterized.These results will facilitate the understanding of the molecular genetic mechanism of fiber elongation with regards to miRNAs in cotton.

1.Introduction

The cotton(Gossypium spp.)fiber is the most important raw material for the textile industry.The cotton fiber is an ideal model system for studying cell elongation,the cell wall,and cellulose synthesis.Cotton fiber morphogenesis can be divided into four stages:initiation,elongation,secondary wall thickening,and maturation[1].Fiber length(FL)is one of the most important fiber characteristics and is determined primarily in the elongation stage from the day of anthesis to 21-26 days post-anthesis(dpa)[2].Many studies have identified important genes that are involved in the development of cotton fibers[3].For example,the expression levels of genes encoding sucrose and K+transporters are high during the elongation stage[4],and reactive oxygen species(ROS)-regulated Ca2+signaling also affects fiber elongation[5].Another gene,AnxGb6,encoding annexin,affects fiber development by interacting with actin 1,and silencing AnxGb6 expression leads to shorter fibers[6].Plant hormones such as auxin play an important role during the fiber initiation and elongation stages[7].A study comparing wild-type fibers and fiberless(fl)mutant ovules at 10 dpa found that＞100 genes are involved in fiber elongation,including a putative vacuolar(H+)-ATPase catalytic subunit and several arabinogalactanlike proteins[8].

Protein-coding genes may be negatively regulated by microRNAs(miRNAs),a type of non-coding small RNA,which play an important role in regulating gene expression at the post-transcriptional levelin plants.In cotton,111 miRNAs and 2 novelmiRNAs were identified in 0-10 dpa ovules in Xuzhou 142 and its fuzzless/lintless mutant[9],and the study also demonstrated the importance of miRNAs in regulating the development of fiber cells.Similarly,27 conserved miRNA families and 4 novel miRNA families were reported[10]in a study of leaves and ovules in the Upland cotton(G.hirsutum)genetic standard TM-1 and its non-fiber mutant,and 223 target genes were found to be expressed.Other recent studies[11,12]revealed findings in the same genotypes.Additional studies[13,14]used 7-dpa fibers of Xuzhou 142 and 5-20-dpa fibers in CRI 35 to demonstrate that miRNAs regulate the elongation of fiber cells.Most recently[15],miR828 and miR858 were reported to regulate the MYB2 gene expression in Arabidopsis trichomes and cotton fiberinitiation.Overexpression ofmiR157 resulted in smaller floral organs,fewer ovules and decreased seed production by suppressing SQUAMOSA promoter-binding protein-like(SPL)genes[16].However,there remains a lack of information about natural genetic and expression differences of miRNAs in elongating fibers of normal cotton genotypes differing in fiber quality traits.

G.barbadense has much longer,stronger,and finer fibers than G.hirsutum,and introgression lines or backcrossed inbred lines(BILs)with better fiber quality traits in G.hirsutum have been developed by advanced backcrossing[17].Researchers have studied differential gene expression of miRNAs,pre-miRNAs and their putative targets between the two species[18-20].However,their genetic relationship with fiber quality traits such as FL is unclear.

The objectives of this study were to identify differentially expressed miRNAs between two introgressed BILs differing in FL and FL QTL alleles.The two BILs were selected based on our prior QTL mapping results[17]for fiber length in a BIL population to minimize the confounding effects caused by different genetic or species backgrounds.The comparative small RNA transcriptome analysis was followed by a bioinformatic,a degradome,and a qRT-PCR analysis.

2.Materials and methods

2.1.Plant materials and tissue collection

Two BILs,a “Long”(SRR1630807)and a “Short”(SRR1630808)BIL with a statistically significant difference in FL were chosen to identify FL-associated miRNAs.The BILs were developed from a cross between G.hirsutum SG 747 as the recurrent parent and G.barbadense Giza 75,followed by two generations of backcrosses and then four generations of self-pollination at New Mexico State University,Las Cruces,NM,USA.They were grown side by side on the experimental farm of the Institute of Cotton Research,Chinese Agricultural Academy of Science,Anyang,Henan,China,in 2014.Each line was planted in three replications with 20 plants in each replication.Developing ovules at 0 dpa and fibers at 10,15,20,and 25 dpa were collected from each replication in the early morning after flowers were tagged at 0 dpa.The excised ovules and fibers were immediately frozen in liquid nitrogen and stored at-80°C until use.

2.2.Small RNA sequencing and library construction

Total RNA was extracted from each replication using the PureLink Plant RNA Reagent kit(Invitrogen,Carlsbad,CA,USA)according to the manufacturer's instructions.All RNA samples were quantified using a Nanodrop ND 2000 spectrophotometer(NanoDrop,Wilmington,DE,USA).Three RNA samples from each 10-dpa fiber sample of the two BILs were equalized for analysis.For degradome sequencing,total RNAs from 10 dpa fibers of “Long”and “Short”were mixed in an equal molar ratio as one sample.Poly(A+)RNA molecules were isolated,and a 5′RNA adaptor was ligated to the poly(A+)RNA by T4 RNA ligase.Finally,Illumina sequencing technology(Illumina,San Diego,CA,USA)was employed to perform RNA sequencing by BGI(Shenzhen,China).

2.3.Analysis of sequencing data

Small RNA reads,degradome reads,and transcriptome reads were produced using a HiSeq 2000(Illumina,San Diego,CA,USA)at BGI.The small RNA reads were processed into clean,full-length reads by removal of low-quality reads and reads with poly-A tails.Adapter sequences were trimmed from the remaining small RNA sequences,and reads longer than 30 nt or shorter than 18 nt were discarded.The remaining highquality reads were used for further analysis,and small RNA reads were deposited in NCBI with accession number SRX743622.

For the degradome reads,low-quality sequences and adapters were removed,and 20-21-nt sequences were collected into clean reads for a further analysis.The degradome sequencing reads were deposited in NCBI with accession number SRX743623.

The RNAs from 10 dpa fibers were sequenced,with approximately 4.6 Gb of transcriptome data generated for each genotype,using a HiSeq 2000 at BGI.A total of 49,508 and 49,448 genes were annotated in “Long”and “Short”,respectively,by alignments to the sequenced TM-1 genome(http://cgp.genomics.org.cn/page/species/index.jsp) using SOAP2[21].The transcriptome sequence reads were deposited in NCBI with accession number SRP039385.

2.4.Identification of known miRNAs and prediction of novel miRNAs

To identify known miRNAs,the small RNAs were aligned to known miRNAs of all plant species in miRBase(http://www.mirbase.org/index.shtml,release 21.0).Only sequences with no＞2 nucleotide(nt)mismatches were assigned as known miRNAs.

The small RNA sequences were aligned to cotton genome survey sequences(GSS)from the G.hirsutum,G.raimondii,G.arboreum,and G.barbadense genomes(http://cgp.genomics.org.cn/page/species/index.jsp, http://database.chgc.sh.cn/cotton/index.html)to identify potentially novel miRNAs,with fewer than three nt mismatches.Potential pre-miRNAs and secondary structures from the genomic sequences were examined with mireap_0.2(http://sourceforge.net/projects/mireap).Novel miRNAs were identified based on the following standards[22]:(1)the novel miRNA could be cut from the stem of the stem-loop structure,(2)there were fewer than four nt mismatches in the miRNA/miRNA*duplex,and(3)asymmetric bulges were minimal in size(≤2)and frequency(≤1),especially within the miRNA/miRNA*duplex.

MiRNAs with＞10 transcripts per million(TPM)reads in at least one library were selected for a cluster analysis.Based on the Fisher's test performed to compare the two BILs,miRNAs showing at least two fold changes in expression level were assigned as DE miRNAs[11,23].Heat maps were constructed to visualize the DE miRNAs.

Co-localization of miRNAs and QTL hotspots was performed using the Cotton QTL database (http://www.cottonqtldb.org/,release 1.0)based on markers anchoring the regions of QTL hotspots.A BLAST search was performed to locate the above DE miRNA genes and previously reported FL QTL[17]and FL QTL hotspots[24,25].

2.5.Small RNA target prediction

The psRNATarget web server(http://plantgrn.noble.org/psRNATarget/)was used for target prediction.The parameters used were based on previously published papers[26,27].When the alignment score was no＞3 nt,the transcripts were assigned as targets.The most abundant sequence from each miRNA family served as the query and the G.hirsutum,G.raimondii,and G.arboreum genomes(http://cgp.genomics.org.cn/page/species/index.jsp)served as the database.MiRNAs showing at least 1.5 fold changes in expression level were assigned as DE.A BLASTx search was then performed against the NCBI non-redundant protein database of flowering plants(taxonomic ID:3398)to annotate the potential targets.The top three hits(P＜0.001)were chosen as the annotation of the potential function of each predicted target.Finally,Gene Ontology (GO)(http://www.geneontology.org/)categories were used to evaluate the potential functions of the targets using Blast2GO[28].

Unique sequence signatures from the degradome data were aligned to the 10-dpa cotton fiber transcriptome database(NCBI accession number SRP039385)using Perl scripts developed by BGI.These signatures were used to identify the miRNA targets with mismatches of no＞3 nt[29].

2.6.A RNA ligase-mediated rapid amplification of cDNA ends(RLM-RACE)

The GeneRacer kit(Invitrogen,Waltham,MA,USA)was used to perform RLM-RACE.Total RNA(5 μg)from equal mixtures of fibers of“Long”and“Short”were ligated to RNAadapterwithout calf intestinal phosphatase treatment.The cDNA were transcribed using the GeneRacer Oligo dT primer.PCR was performed with 5′adapter primers and 3′gene-specific primers according to the manufacturer's instructions(Table S1).

2.7.The qRT-PCR of miRNAs and their targets

Two micrograms of total RNA from 8 tissues each with 3 replications(0 dpa ovules,and 10,15,20,and 25-dpa fibers)were used for qRT-PCR analysis of miRNAs based on stem-loop reverse transcription PCR[30].TaqMan MicroRNA Assays(Applied Biosystems,850 Lincoln Center Drive,Foster City,CA,USA)and One Step SYBRPrimer Script PLUS RT-PCR kit(Bio Inc.,Nojihigashi,Kusatsu,Shiga,Japan)was used to conduct the qRT-PCR of the miRNA and targets,respectively.The reactions were incubated at 95°C for 10 min,followed by 40 cycles of 5 s at 95 °C and 40 s at 58 °C.Cotton U6 snRNA and His3(histone 3)were used as internal controls for miRNA and mRNA,respectively.Relative expression levels were calculated by the 2-ΔΔCtmethod.All primers are listed in Table S1.

3.Results

3.1.Overview of high-throughput sequencing data

The fiber elongation rate at 10 dpa was highest,with a significant difference in FL between the two BILs[31].To characterize the role of active miRNAs during the fiber elongation stage,deep-sequencing libraries of small RNAs were generated using total RNAs extracted from 10 dpa fibers in the two BILs.In total,17,941,998 and 24,935,739 sequence reads of 18 to 30 nucleotides in length were obtained in the“Long”and “Short”libraries,respectively,using the Illumina sequencing platform.After the low-quality or contaminated reads were removed,15,036,558 and 20,796,489 clean reads remained in the “Long”and “Short”libraries,respectively.The most abundant small RNAs clustered at 24 nt,representing 37.7-37.9%of the reads.The 21-nt reads were the second most abundant class,followed by 22-nt and 23-nt reads(Fig.S1).This result was consistent with findings in rice[32],peanut[33],soybean[34],Medicago truncatula[35],and cotton in other studies[10,36].

Ofthe unique reads in“Long”and“Short”,respectively 74.1%and 73.7%were matched to the sequenced allopolyploid TM-1 genome[37,38].As a comparison,only 24.7%and 25.1%were matched to the sequenced diploid cotton D5genome[39],while only 25.4%and 25.9%were matched to the sequenced diploid cotton A2genome[40].This finding indicates a difference of small RNAs between the At and Dt sub-genomes of G.hirsutum and its ancestral A2or D5genome.The mapping result also shows some similarities in small RNAs between the A2and D5genomes.After the small RNAs were clustered based on annotations to GenBank and Rfam(version 12.0)[41],they were divided into the following classes:rRNA,scRNA,snoRNA,snRNA,tRNA,and miRNA(Table S2).

Degradome libraries for “Long”and “Short”were constructed and sequenced with an Illumina Genome Analyzer[29,42].In total,20,544,138 sequence tags were assayed,containing 99.7%clean tags.After removal of rRNA,snRNA,snoRNA,and tRNA based on annotations of the tags using Rfam and GenBank,the remaining sequences were matched to the annotated genes of the sequenced TM-1 genome[37,38](Table S3).Finally,2,272,664(74.5%)unique reads were mapped to the reference database for detecting candidate targets of miRNAs.

3.2.Known miRNAs in elongating fiber cells of the two BILs

To identify conserved miRNAs within Gossypium,the small RNAs were compared with the miRNA database(miRBase version 21.0,released on July 9,2014)[43].In total,334 conserved cotton miRNA belonging to 65 miRNA families comprising 3,837,205 individual candidate miRNA reads were identified in the small RNA libraries of the two BILs.The miRNA read abundances were measured by normalizing the results to tags per million(TPM)reads.The two miRNA families with the greatest abundance were miR167_1 and miR156,accounting for respectively 48.4%and 22.2%of all the known miRNA reads in “Long”,and 51.4%and 26.9%in“Short”.The miRNA families MIR166,MIR390,MIR164,and MIR172 also showed average abundances with＞500 TPM in both libraries.The expression levels of other miRNA families were much lower.The abundances of many miRNAs differed between the two BILs.For example,MIR167_1 and MIR156 were expressed at a higher level in “Short”than in “Long”,whereas miR166,MIR390,MIR164,and MIR172 were expressed at higher levels in “Long”(Fig.1).

The remaining small RNAs were subjected to a BLASTn search for known miRNA sequences from other plant species deposited in miRBase,leading to the identification of 277 additional small RNAs.For this analysis,these sequences were first mapped to the sequenced TM-1 genome and then the secondary structures of their precursor sequences were predicted to search for hairpin structures.A total of 46 miRNAs containing 3,835,736 reads belonging to 22 miRNA families were identified.Thus,together with the above 334 conserved miRNAs in cotton,380 known miRNAs were identified.

Finally,the expression levels of all the known miRNAs in the two BILs were compared,resulting in the detection of 27 known miRNAs that showed significant expression differences between the two BILs in 10-dpa fibers.Of these miRNAs,81%(22/27)were expressed at higher levels in “Long”than in“Short”(Fig.2-a).

3.3.Identification and verification of novel miRNAs

To identify possible novel miRNAs,the remaining small RNAs were mapped to the sequenced TM-1 genome and subjected to a rigorous secondary structural analysis of their precursors.The minimum free energy(MFE)was calculated.The MFE ranged from-18.1 to-169 kcal mol-1,a range consistent with previous studies[10,14,19].In total,83 novel miRNAs were identified,including 66 in “Long”and 65 in “Short”.

To ensure the accuracy of the identified novel miRNAs,the length and nucleotide bias were further analyzed(Fig.S2).The 21-nt miRNAs were the most dominant class expressed among all the novel miRNAs,followed by 22-nt miRNAs.The first nucleotide of all these 83 novel miRNAs(5′terminal)was most commonly U,in agreement with previous reports in cotton[10,23].In addition to the hairpin structure prediction,detection of miRNA*s(located on the opposite stem-arms of the precursor secondary structure)is another way to confirm novel miRNAs[22].Of the 83 novel miRNAs,39 contained miRNA*sequences.

Among the 83 novel miRNAs,20 were significantly DE between“Long”and “Short”in 10-dpa fibers(Fig.2-b).Relative to “Short”,15(75%)novel miRNAs were upregulated and five downregulated in “Long”,a result similar to that for the known miRNAs.These results suggest that most DE miRNAs(both known and novel)promote fiber elongation.

3.4.Differentially expressed miRNA genes co-localize with fiber-length quantitative trait loci(QTL)and hotspots

To investigate whether the DE miRNAs affect FL,the 4.6%(108/2249)SSR markers that were polymorphic between the two BILs,four FL QTL mapped in the BIL population from which the “Long”and “Short”BILs were derived[17],and 13 FL QTL hotspots in an updated CottonQTLdb,Release 2.0(http://www.cottonqtldb.org/)[24,25],were mapped on the sequenced TM-1 genome.The 47 DE miRNA gene sequences were placed at 86 locations in the TM-1 genome,15 of which on unmapped scaffolds.Only seven DE miRNA genes colocalized with seven FL QTL and QTL hotspots,on chromosomes D05,A07,D08,A11,D11,A12,and D12(Fig.3).Two(miR477b and ghr-miR479)of these colocalized with two FL QTL on chromosomes D11 and A12 previously[17]identified in the BILs.Both of them showed higher expression in “Long”than in “Short”.

Fig.1-Expression of known miRNA families in “Long” and “Short” BILs.

Fig.2-Differentially expressed(DE)known miRNAs and novel miRNAs during the fiber elongation stage at 10 dpa.Differential expression was considered as an at least 2-fold change.(a)Hierarchical cluster analysis of DE known miRNAs during the fiber elongation stage.(b)Hierarchical cluster analysis of DE novel miRNAs during the fiber elongation stage.

3.5.Target prediction and degradome sequencing validation

To characterize the functions of DE miRNAs in cotton fiber elongation,their targets were predicted using all annotated genes of the sequenced TM-1 genome and their expression levels were compared using the 10-dpa fiber transcriptome(NCBIaccession number SRP039385)ofthe same two genotypes.In total,82 genes were identified as candidate target genes of the seven miRNAs co-localizing with QTL(Table S4).The targets for the two BILs were similar,in that“Long”and “Short”carried 79 and 80 target genes,respectively.However,the expression levels of 21 target genes differed significantly between the two BILs.These included genes responsive to auxin,which were targeted by miR160a-5p;several ubiquitin proteins and ACC oxidases targeted by gra-miR477;several genes in the cation channel pathway targeted by novel_miR1206,and several kinases targeted by novel_miR262,among others.Interestingly,a large number of genes associated with biological and environmental stresses were predicted to be targets,including the RPM1-interacting protein(RIN)family and protein kinase superfamily protein(PBS1),suggesting that these classes of genes participate in the development of cotton fiber cells.Of the 21 DE targetgenes,10 were downregulated in“Long”relative to “Short”(Table S4).All 82 target genes were further mapped onto the sequenced TM-1 genome,including nine located in the seven FL QTL and hotspot regions(Fig.3).These nine targets were annotated as genes coding for ARF10,a root hair defective 3-like protein,and an ARF-GAP domain,among others.Only one(CotAD_61076)of these nine targets showed differential expression between the two BILs,but had no annotation.

To validate the accuracy of the predicted results,a degradome analysis was performed.The degraded reads were mapped onto the predicted target sequences with no mismatch.The result showed that among all the 82 predicted target genes,25 genes were verified in the degradome(Table S4).This relatively moderate level of congruence may have been observed because the expression level of many of the targets or the abundance of the cleaved targets was too low to detect.Furthermore,miRNAs can regulate their targets at the translational rather than the transcriptional level[44].A RNA ligase-mediated rapid amplification of cDNA ends(RLMRACE)was further performed to verify the cleavage sites of six selected targets.The results showed that 3 target genes,CotAD_42057,CotAD_26228,and CotAD_17759,were contained in the results of degradome,indicating a relative consistency between the two approaches.Cleavage sites for some gene sequences were not unique(Fig.S3),consistent with a previous report[45].

3.6.Genome distribution of miRNAs and targets in tetraploid cotton and diploid G.arboreum and G.raimondii

The modern descendants,i.e.,G.arboreum(A2)and G.raimondii(D5)of the two ancestral diploid cotton species that hybridized to create the five allopolyploid cotton species including two cultivated(G.hirsutum and G.barbadense),differ twofold in genome size[40].The homoeologous genes(those at the same loci on homoeologous chromosomes in A2and D5)of the predicted targets were identified from the sequenced genomes of G.arboreum and G.raimondii(Table S5).Interestingly,in contrast to the case in the tetraploid species,nearly one third of homoelogous target genes in the A2or D5genomes were not the predicted targets of the DE miRNAs identified in the present study in tetraploid Upland cotton(Table S5),suggesting that more target genes may be regulated in tetraploid cotton than in its ancestral diploid species.For example,three target genes could be regulated by miRNAs in A2but not D5,while nine target genes could be regulated by miRNAs in D5but not A2.Seven targets in the tetraploid had no homologous genes in either D5or A2,while 14 targets in the tetraploid had homologous genes in D5and/or A2.Among the 14 targets,three could be regulated only in A2,nine only in D5,one in both D5and A2,and one in neither D5nor A2.

Fig.3-Chromosome locations of DE miRNAs,FL QTL,and FL QTL hotspots in the genome of Gossypium hirsutum.“FL_QTL_Hotspots(10)” refers to an FL QTL hotspot and the number in parentheses refers to the number of QTL within the hotspot.miRNAs in red were located in QTL or QTL hotspots.

Fig.4-Expression profiling of selected miRNAs and their targets during cotton fiber development.The error bars denote the standard deviations of multiple measurements.

3.7.Expression profiling of selected miRNAs and their targets

The expression patterns of miRNAs and their putative targets may vary during cotton fiber development.To further investigate the expression patterns of miRNAs,their target genes,and the relationship between them,seven DE miRNAs located in seven FL QTL and hotspots,and their targets,were chosen for a qRT-PCR analysis in developing cotton ovules at 0 dpa and in developing fibers at 10,15,20,and 25 dpa.The targets were selected from the predicted targets listed in Table S4.As shown in Fig.4,gra-miR477,miR157d-3p,miR160a-5p,novel_miR1206,miR477b,and ghr-miR479 were expressed at higher levels in “Long”than in “Short”during the fiber elongation stage.The expression level of novel_miR262 became lower with fiber elongation in “Long”,whereas it was highly expressed at 15 dpa in “Short”.The targets showed an expression trend almost opposite to that of the miRNAs(Fig.4).The only exception was CotAD_42057,the target of miR157d-3p,which showed high expression at 10 dpa in“Long”.Considering that other pathways may regulate target gene expression,CotAD_42057 was still expressed at a lower level in “Long”than in “Short”.Thus,almost all the results of qRT-PCR showed that miRNAs were expressed at higher levels in “Long”than in “Short”.

4.Discussion

4.1.miRNAs and their targets in tetraploid cotton

In plants,miRNAs regulate nearly all growth and development,including cotton fiber development.The most recent research has focused on fiber initiation,but understanding of fiber elongation is relatively poor.To our knowledge,previous studies focused only on the expression profiling of miRNAs in a single normal genotype of cotton and its comparison with its fiberless and/or fuzzless mutant.This may be because there were no suitable genotypes with the same genetic background but differing in FL to be used to study fiber elongation without interference from variation in other traits.The present study is the first to employ two BILs from a BIL population with the same genetic background to study the activity of miRNAs during fiber elongation.The two BILs can be considered near-isogenic lines(NILs)based on their high genetic similarity(95.4%based on 2349 markers),as per many other studies using early generations of backcrossing to develop NILs[46,47].

In this study,by the use of G.hirsutum genome sequence information,380 known miRNAs and 83 novel miRNAs were identified in 10-dpa developing fibers.Of these,27 known and 20 novel miRNAs were DE between the two BILs differing in FL.Among all the DE miRNAs,about 78%of the miRNAs were expressed at higher levels in “Long”than in “Short”.Given that miRNAs usually downregulate the expression of their targets,we accordingly hypothesize that the target genes of these DE miRNAs play a negative role in the development of cotton fibers,resulting in longer fibers.But further experiments are required to address this hypothesis.Among the four FL QTL identified previously in the BIL population from which the two BILs were selected,and 13 FL QTL hotspots(containing three of the four FL QTL identified in the BILs),seven DE miRNAs were co-localized with seven FL QTL and hotspots on the genome of upland cotton.The approach of integrating information from miRNA sequencing with QTL linkage mapping and genome sequencing has shed light on the genetic mechanism of fiber elongation.However,it should be pointed out that,although other miRNAs as genetic candidates for FL QTL were ruled out owing to their failure to co-localize with FL QTL,the seven miRNAs co-localizing with FL-QTL should not be considered the only candidate genes for the seven FL-QTL,because these QTL regions are large and may contain several hundred genes.To match a miRNA with a single QTL,one or more large mapping populations would be needed to fine-map each FL QTL to a much narrower region.In the future,a small number of genes coding for miRNAs,if any,and other RNAs or proteins in the regions will be further analyzed for their candidacy as the genetic determining factors for FL.

4.2.Differential expression of miRNAs and their targets among tetraploid cotton and its ancestral diploid cottons

Most of the families of known DE miRNAs identified in the two BILs in this study are conserved in plants.However,there were also 20 DE novel miRNAs enriched in one or the other BIL.After the reads identified by searching other plant miRNAs from miRBase were mapped to the genome of G.hirsutum,it was observed that several known miRNAs in other plants were species-specific,in that only a few miRNAs(46 of 277 in other species)could be matched to G.hirsutum with a hairpin structure.This observation suggested that numerous miRNAs identified in other plants may not be miRNAs in G.hirsutum.These reads may be of other types of small RNA produced in different ways in different species.Testing this speculation would require further experiments.

A recent study[48]compared the conserved miRNAs in diploid cotton species,indicating that most conserved miRNA families are ancient and stable over vast evolutionary timescales.Thus,their functions are conserved across species in Gossypium.While the evolution of miRNAs is not only manifested at the expression level but also their targets,tetraploid cotton also requires specific miRNAs to regulate redundant or novel genes that have gained new functions after the merging of its two ancestral diploid genomes.For example,miR828 regulates MYB2D via ta-siRNAs in tetraploid cotton,although a small amount of MYB2A is cleaved[15].Because tetraploid cotton contains both the At and Dt subgenomes,it is possible that miRNAs from one diploid cotton genome could also regulate target genes originally belonging to the other diploid genome,or vice versa.To further study the relationship of miRNAs and their targets in tetraploid cotton and its diploid cotton ancestors,we investigated whether the seven co-located DE miRNAs regulate the same homologous targets in G.arboreum and G.raimondii.Somewhat unexpectedly,we found that nearly one third of the homologous genes in the diploids could not be targeted by the initial miRNAs identified in the tetraploid cotton,although the miRNAs and their targeting sequences are conserved between homologous genes.In our study,for novel_miR1206 identified in both tetraploid cotton and the two diploid cottons,one of its targets,CotAD_23046 in tetraploid cotton,has homologous genes in its diploid cotton ancestors,but neither of them was a target of novel_miR1206.In comparison with the potential miRNA-mediated gene regulation model in tetraploid cotton proposed recently[49],our study suggests another model in which several miRNAs did not regulate the same target genes in diploid G.arboreum or G.raimondii until tetraploid cotton evolved.However,more experimental data are needed to support this model.

4.3.A potential functional network of miRNAs associated with fiber elongation

The present study identified seven DE miRNAs,co-localizing with FL QTL,that may be associated with the elongation of cotton fiber cells in the two BILs differing in FL.Among known miRNAs,few if any have been reported to be associated with fiber elongation.Previous studies[50,51]have suggested that CaMBP regulates the elongation of cells in Nicotiana tabacum and the cotton fiber.In our study,novel_miR262 was predicted to regulate a calmodulin-binding protein.The expression level of novel_miR262 gradually decreased in “Long”during fiber elongation,whereas it was expressed at a low level during the early stage of fiber elongation in “Short”.As expected,the expression level of its target was lower in “Long”than in“Short” during the early elongation stage. DELLA(Cot_AD17273),a target gene of miR477b,may reduce fiber elongation in concert with HOX3[52].This target was expressed at a higher level in “Short”than in “Long”during fiber elongation,a trend opposite to that of miR477b.It is known[53]that DELLA protein can be degraded in response to gibberellin,and the present study showed that miR477b could also regulate the expression level of DELLA protein at the level of transcription.

miR157d-3p,miR160a-5p and novel_miR1206 were expressed at higher levels in “Long”than in “Short”BIL.Their annotated targets were clustered in disease responses and hormone-mediated signaling pathways.A target of miR157d-3p was predicted to be the gene coding for a potassium channel protein,whereas the gene for auxin response factor(ARF)was predicted to be the target of miR160a-5p and novel_miR1206.ARF is a member of the auxin metabolic pathway,which is known to play an important role in elongating fibers[54].Furthermore,both targets were regulated in G.arboreum by the miRNAs,in agreement with a previous study[49]showing that homologous genes in the Dgenome can contribute to cellelongation.miR157d-3p was also highly expressed during the fiber elongation stage in both G.hirsutum CRI 35[14]and G.barbadense 3-79[45].gra-miR477 was expressed ata higher level during 10-15 dpa in “Long”than in“Short”,whereas its target showed the opposite trend in “Long”and “Short”.The result is consistent with those of previous findings that gra-miR477 was highly expressed in 10-20-dpa fibers in G.hirsutum CRI 35[14]and in 7-25-dpa fibers in G.barbadense 3-79 [45].Its target is annotated as a gene coding for oligopeptidase,which was reported to be associated with the growth of adopted cells[55].In our study,several genes associated with abiotic and biotic stresses such as NBS-LRRwere identified by target prediction and degradome analyses,suggesting that the development of cotton fiber cells may have similarities to stress responses.

5.Conclusions

This study co-localized seven DE miRNAs and their nine targets with seven FL QTL.This result will provide a framework for further studies of the molecular and genomic genetics of fiber elongation when FL genes are introduced from G.barbadense to G.hirsutum for improvement of fiber quality in cotton.

Acknowledgments

The research was supported by grants from the National Natural Science Foundation of China(31621005),the National Key Research and Development Program of China(2016YFD0101400),the National Research and Development Project of Transgenic Crops of China(2016ZX08005005),and the New Mexico Agricultural Experiment Station.

Appendix A.Supplementary data

Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2018.08.004.