Hantao Wang, Ruiting Zhang, Chao Shen, Ximei Li, De Zhu, Zhongxu Lin,*

aState Key Laboratory of Cotton Biology,Institute of Cotton Research of CAAS,Anyang 455000, Henan,China

bNational Key Laboratory of Crop Genetic Improvement, College of Plant Science and Technology, Huazhong Agricultural University, Wuhan 430070,Hubei, China

cCollege of Agronomy and Plant Protection, Qingdao Agricultural University, Shandong Key Laboratory of Dryland Farming Technology,Qingdao 266109,Shandong,China

Keywords:Upland cotton Fiber quality RNA-seq QTL Candidate genes

ABSTRACT With increasing demand for high-quality cotton,it is desirable to identify genes involved in fiber development for molecular improvement of cotton. In this study, 780 differentially expressed genes (DEGs) were identified in developing fibers at 10 days post-anthesis (DPA)in Gossypium hirsutum acc. DH962 and G. hirsutum cv. Jimian 5 using RNA-seq. Of 15 stable QTL for fiber quality identified in the same two parents in previous studies,4,3,6,1,and 1 QTL were associated with fiber length (FL), fiber strength (FS), micronaire (MIC), fiber elongation(FE)and fiber length uniformity ratio(FU),respectively.Integration of DEGs and QTL allowed the identification of 31 genes in 9 QTL regions, of which 25 were highly expressed in fibers based on the transcriptome datasets and 9 were preferentially expressed in different stages of fiber development. Gh_A01G0453 (GhDTX19), Gh_D07G1581 and Gh_D04G0942 were expressed specifically in 5 and 10 DPA fibers, with Gh_D04G0942 showing low expression in other tissues except pistil. Gh_D07G1799 (GhGAUT9),Gh_D11G0326 (GhVPS29), Gh_D11G0333 (GhTCP14), and Gh_D11G0334 (GhNRP2) were preferentially expressed in 5 or 10 DPA fibers; Gh_A01G0397 (GhABCG10) and Gh_D07G1744 were expressed specifically in 20 and 25 DPA fibers. These results suggest candidate genes for molecular improvement of cotton fiber quality.

1. Introduction

Cotton(Gossypium spp.)is the most important natural fiber crop in the world. There are four cultivated species including tetraploids of G. hirsutum (AD1) and G. barbadense (AD2), and diploids of G. herbaceum (A1) and G. arboreum (A2), which are cultivated in ＞80 countries[1].G.hirsutum(Upland cotton)is the main cultivated species in cotton, accounting for over 95% of world cotton fiber production[2].Although fiber quality and yield have been greatly improved in the past several decades by conventional cultivar breeding programs, demand for fiber quality and yield of cotton has increased with the development of the modern textile industry and needs for diverse cotton products.

Cotton fiber quality traits include fiber length (FL), fiber strength(FS),fiber length uniformity ratio(FU),fiber elongation(FE), micronaire (MIC), and others [3]. They are complex quantitative traits controlled by multiple genes. The rapid development of molecular-marker technology and nextgeneration sequencing provide opportunities for QTL mapping and gene cloning in cotton.In the past two decades,many QTL for fiber quality traits have been reported[3-9].But it is difficult to use these QTL for further research,owing to the use of diverse genetic maps and markers for their discovery.Said et al.[10,11]collected 635 QTL associated with fiber quality from 42 reports and identified 37 QTL hotspots,of which respectively 21 and 27 were associated with fiber quality in intraspecific and interspecific populations.A total of 2091 QTL associated with fiber quality distributed on 26 chromosomes have been collected in the Cotton QTL Database(http://www2.cottonqtldb.org:8081/index).However, owing to the complexity of the tetraploid cotton genome,no major QTL for fiber traits have been cloned to date.

Because of the difficulty of forward genetics research in cotton,RNA-seq has been used to identify candidate genes for fiber development [12-15]. Some novel genes functioning in fiber initiation and elongation have been verified by molecular biology methods: for example, ACTIN1 [16], GhMYB25 [17],GbPDF1 [18], GhCaM7 [19], GhAnn2 [20], GhHOX3 [21], and GbEXPATR [22]. Li et al. [22] identified an α-expansin family protein GbEXPATR using a normalized cDNA library of G.barbadense 3-79 fiber. GbEXPATR overexpression lines showed a lower micronaire value and greater(3.8%-5.8%)fiber strength and fiber length (5.9%-7.7%) than the wild type. QTL mapping could provide candidate regions for the identification of candidate genes, and RNA-seq could reveal transcription information for these genes in different development stages.The integration of candidate regions and transcriptomic information was useful for suggesting candidate genes in cotton [23-26]. Liu et al. [23] identified three candidate genes associated with fiber quality by integrating the results of fine mapping and RNA-seq. Wu et al. [24] identified 212 DEGs mapped in the regions of 24 yield QTL and 11 yield trait QTL hotspots using a microarray-based comparative transcriptome analysis in 10 DPA fibers. Islam et al. [25] found 3 QTL regions controlling three fiber-quality traits and suggested that receptor-like kinase pathway genes might be candidates for superior fiber strength and length in the cotton line MD52ne.

In our previous study,a high-density genetic map containing 1013 loci and 12 QTL for fiber quality was developed using a F2population from a cross between G.hirsutum DH962 and Jimian 5[27].An RIL population developed from the same two parents was phenotyped in eight environments, revealing 64 QTL for fiber quality [8]. Two immortalized backcross populations were then developed using the RIL population crossed with the two parents(DH962 and Jimian 5), respectively, and 76 QTL for fiber quality were detected[3].These previous studies revealed stable QTL for fiber quality [3,8,27]. In the present study, stable QTL intervals from the three previous studies were identified. Transcriptome data from developing fibers at 10 days post-anthesis(DPA)from the two parents, DH962 and Jimian 5, were used to identify differentially expressed genes(DEGs).Molecular information for fiber development wss obtained by integration of DEGs and stable QTL.

2. Material and methods

2.1. Plant materials

DH962,a G.thurberi introgression line,shows high fiber quality[8]. Jimian 5 is a cotton cultivar with high yield [8]. In our previous studies, an F2population, a RIL population and two immortalized backcross populations crossed with the two same parents, DH962 and Jimian 5, were used to identify stable QTL[3,8,27].Cotton bolls of DH962 and Jimian 5 were harvested at 5,10, 15, 20, and 25 DPA. Fibers were separated from ovules in liquid nitrogen. Each sample was divided into three biological replicates and stored at-80 °C until further processing.

2.2. Total RNA isolation, library construction, and sequencing

Fiber samples at 10 DPA obtained from DH962 and Jimian 5 were used for total RNA isolation. Total RNA was extracted using a Sigma Spectrum Plant Total RNA kit (Sigma-Aldrich,St. Louis, MO, USA). Sequencing libraries were constructed using the Illumina TruSeq Stranded RNA Kit (Illumina, San Diego, CA, USA) following the manufacturer's recommendations. The quality and quantity of the libraries were assessed with an Agilent 2100 Bioanalyzer and an ABI StepOnePlus Real-Time PCR System, respectively. DGE profiling data were performed on the Illumina HiSeq 2000 sequencing platform by the Beijing Genomics Institute (Shenzhen, Guangdong,China). Each sample contained two biological replicates.

2.3. Detection of DEGs

Raw sequence data were filtered by removing reads with adaptors, reads containing N ＞10%, and reads with quality value ≤5 base proportion ＞50%. The remaining reads were mapped to the G.hirsutum genome[28]using the Tuxedo suite of tools [29]. Clean reads from each sample were mapped to the genome with Bowtie 2.3.1 (http://bowtie-bio.sourceforge.net/bowtie2/index.shtml) and TopHat 2.1.1 (http://ccb.jhu.edu/software/tophat/index.shtml).Transcripts were predicted using Cufflinks [29] and gene expression levels were calculated as reads per kilo-base of transcript sequence per million base pairs sequenced (FPKM). Genes with P ＜0.05 and an absolute value log2ratio ≥1 were used to select DEGs.

2.4. Functional annotation of DEGs

Gene ontology (GO) enrichment analysis of differentially expressed genes was performed using OmicShare tools (http://www.omicshare.com/tools/Home/Index/index.html) with P ＜0.01. Kyoto Encyclopedia of Genes and Genomes (KEGG)pathway analysis of all DEGs was performed at the Cotton Functional Genomics Database (CottonFGD) (https://cottonfgd.org/)and significant pathways were selected at P ＜0.05.

2.5. Physical locations of QTL regions

Marker sequence data were obtained from the CottonGen database(https://www.cottongen.org/data/download/marker)and our previous study[27].Candidate intervals of stable QTL were determined based on the sequences of the two flanking markers. The physical locations of markers were obtained using BLAST to align the marker sequence to the G. hirsutum genome[28].

2.6. Gene expression pattern and qRT-PCR analysis

The tissue expression levels of DEGs were obtained from previously reported transcriptome data [28]. Heat map and expression clusters were calculated with Hemi 1.0 (http://hemi.biocuckoo.org/).

qRT-PCR experiments were performed using the ABI Prism 7500 (Foster City, California, USA). GhUBQ7 (GenBank accession number DQ116441) was used as the reference gene.Relative expression levels were determined and normalized using three biological replicates. The primers used in this study are listed in Table S1.

3. Results

3.1. Trait values of DH962 and Jimian 5 in seven environments

Fiber trait values for DH962 and Jimian 5 were obtained in seven environments, except for FE, which was recorded in only six environments [8]. FL ranged from 29.7 to 32.98 mm in DH962 and from 25.9 to 27.91 mm in Jimian 5 in seven environments.FS ranged from 29.23 to 35.66 cN tex-1in DH962 and from 24.6 to 26.58 cN tex-1in Jimian 5.FL,FS,and FU were higher in DH962 than in Jimian 5; while, FE and MIC were not significantly different between DH962 and Jimian 5(Fig.1).

3.2. Identification of DEGs between DH962 and Jimian 5

Four libraries in 10-DPA fibers from DH962 and Jimian 5 were used for expression analysis. DH962 yielded 24,730,647 raw reads in two biological replicates and Jimian 5 revealed 26,813,161 raw reads in two biological replicates. After removal of low-quality reads, the clean reads were mapped to the cotton genome. A total of 780 DEGs between DH962 and Jimian 5 (P ＜ 0.05; log2ratio ≥ 1) were identified (Table S2), of which 272 (34.9%) were up-regulated and 508 (65.1%) downregulated in DH962.

3.3. GO and KEGG enrichment analysis of DEGs

To identify the potential biological functions of DEGs between DH962 and Jimian 5, the 780 DEGs were classified into the three main GO categories (biological process, molecular function, and cellular component) and 35 GO terms (Fig. 2; Table S3). In the biological process category, singleorganism metabolic process (106, 13.59%), oxidationreduction process (64, 8.21%), cellular carbohydrate metabolic process (16, 2.05%), generation of precursor metabolites and energy (11, 1.41%), negative regulation of metabolic process (7, 0.90%), energy derivation by oxidation of organic compounds(6, 0.77%), and electron transport chain (6, 0.77%) were the major subcategories (Table S3). In the molecular function category, catalytic activity (263, 33.72%), oxidoreductase activity (67, 8.59%), electron carrier activity (24, 3.08%),peptidase activity, acting on L-amino acid peptides (21,2.69%), and xyloglucan: xyloglucosyl transferase activity (5,0.64%) were the principal subcategories (Table S3). In the cellular component category, extracellular region (9, 1.15%),photosynthetic membrane (7, 0.90%), thylakoid (7, 0.90%),thylakoid part(7,0.90%),apoplast(5,0.64%),and photosystem I (5,0.64%)were the major subcategories (Table S3).

To further comprehend the enriched pathways of these 780 DEGs,KEGG pathway analysis was performed in CottonFGD.All DEGs were assigned to 22 KEGG pathways(P ＜0.05). Metabolic pathways and biosynthesis of secondary metabolites were the two most common pathways (Table 1). Some of the following pathways were associated with fatty acid metabolism: fatty acid degradation, adipocytokine signaling pathway, beta-Oxidation, acyl-CoA synthesis, and fatty acid biosynthesis.Some pathways were associated with carbohydrate metabolism, including the pentose phosphate pathway, trehalose biosynthesis, and D-glucose 1P ≥trehalose (Table 1). Some of the pathways have been reported to be associated with fiber development, including flavonoid biosynthesis, naringenin ≥pelargonidin,and zeatin biosynthesis(Table 1).

3.4.Identification of stable QTL for fiber quality from the same parents

An F2population, a RIL population, and two immortalized BCF1populations have been developed from the same parents, DH962 and Jimian 5 [3,8,27]. The F2population genetic map contained 1013 loci and 50 linkage groups [27],whereas the RIL population genetic map contained 616 loci and 59 linkage groups [8]. Owing to the limited resolution of these maps, some linkage groups from the same chromosomes were not linked and some QTL did not coincide. For example,qFL-c10-1,qFL-c10-2,and qMIC-c10-1 were assigned to different linkage groups. Of 15 stable QTL for fiber quality(Table 2),4,3,6,1,and 1 QTL were associated with FL,FS,MIC,FE, and FU, respectively. Of which, qFL-c10-1 was detected in four populations, and qFE-c22-1 and qFS-c17-2 were identified in three populations. The remaining QTL were identified in two populations or in more than two environments in a single population.

3.5. DEGs identified in fiber quality QTL intervals

Based on the Upland cotton genome [28], candidate regions were identified for 15 QTL (Table 2). A total of 31 DEGs were located in these regions (Table S4). Three, 2, 3, 4, 4, 3, 6, 1, and 5 DEGs were detected in the qMIC-c1/15-2, qFL-c2, qMIC-c9, qMIC-c10-1, qFS-c15, qMIC-c22, qFS-c17-2, qMIC-c17, and qFUc21-1 intervals, respectively, and no DEGs were detected in the qFL-c10-1, qFL-c10-2, qFS-c17-1, qFL-c21-2, qFE-c22-1, and qMICc25 intervals. Of the 31 genes, four were up-regulated and 27 down-regulated in DH962.

3.6. Expression patterns of candidate genes in fiber quality QTL intervals

Fig. 2 - GO enrichment analysis of 780 DEGs.

Table 1-KEGG pathway enrichment analysis of 780 DEGs between DH962 and Jimian 5 in 10 DPA fiber.

For identifying molecular function,the expression patterns of genes located in QTL regions in 10 tissues (root, stem, leaf,petal,stamen,pistil, and fibers at four developmental stages)were obtained from the transcriptome datasets of the Upland cotton genetic standard TM-1[28].Two genes were expressed at very low levels or not at all in the investigated tissues,and 29 genes with FPKM ≥1 in at least one of the 10 investigated tissues were used for further analysis(Fig.S1).

Of the 29 genes,25 were highly expressed in fibers(Fig.S1;Table 3).Among these,11 genes were highly expressed at fiber elongation stages(5 and 10 DPA),five showed high expression in fiber secondary wall biosynthesis stages (20 and 25 DPA),and nine were prominently expressed in all four fiber samples. Nine of the 25 genes were preferentially expressed at different stages of fiber development(Fig.S1).Gh_A01G0453(GhDTX19) and Gh_D07G1581 were expressed specifically in 5 and 10 DPA fibers and showed low or no expression in other tissues;Gh_D04G0942 was also expressed specifically in 5 and 10 DPA fibers and showed low or no expression in other tissues except pistil; Gh_D07G1799 (GhGAUT9), Gh_D11G0326(GhVPS29), Gh_D11G0333 (GhTCP14), and Gh_D11G0334(GhNRP2) were expressed more strongly in 5 or 10 DPA fibers than in other tissues. Gh_A01G0397 (GhABCG10) and Gh_D07G1744 were expressed specifically in 20 and 25 DPA fibers.

Six of the 29 genes were randomly selected for qRT-PCR verification of their expression levels in fiber development stages(5,10,15,20,and 25 DPA).The expression levels of the six genes in the 10 DPA fiber stage were consistent with theDEG results(Fig.3).All six genes showed expression trends in TM-1 similar to those in the transcriptome data(Fig.3,Fig.S1).

Table 2-Stable QTL for fiber quality from different populations with the same parents.

4. Discussion

Upland cotton is an allotetraploid species (2n = 4× = 52, ≈2.5 Gb),and homologous genes are highly similar between At and Dt subgenomes. Some 64.8% of the genome consists of transposable elements [28], increasing the difficulty of mapbased cloning. As an approach to finding genes,RNA-seq has been used to identify candidate genes associated with fiber development [12-15]. In the present study, RNA-seq revealed 780 DEGs between DH962 and Jimian 5 in the 10 DPA fiber libraries. The low DEG numbers may be due to the limited sequencing depth of DEG profiling technology, and some genes may be not identified. However, some important pathways and genes were identified, including acyl-CoA synthesis,fatty acid biosynthesis, flavonoid biosynthesis and zeatin biosynthesis. In addition, 41 transcription factors associated with fiber development were identified, including MYB,bHLH,NAC,TCP,and ERF.

4.1. Fatty acids and flavonoids associated with fiber cell elongation

In previous studies [30-32], fatty acids and flavonoids were associated with fiber development. Of the DEGs, 3 long chain acyl-CoA synthetase (LACS) genes (Gh_D07G1281, Gh_A07G1181, and Gh_A09G0536) were identified (Table S5). Based on the transcriptome datasets of TM-1 [28], the expressions of Gh_D07G1281 and Gh_A07G1181 were upregulated from 5 DPA fiber to 10 DPA fiber, and downregulated in 20 and 25 DPA fibers (Fig. S2). Gh_D07G1281 and Gh_A07G1181 encode long chain acyl-CoA synthetase 1 (GhLACS1), and are homologous to AtLACS1, which is important for the biosynthesis of very-long-chain fatty acids (VLCFAs) [33]. VLCFAs promote cotton fiber cell elongation [31]. Gh_A05G3400 (GhFLS), Gh_D05G1836 (GhDRF), Gh_D06G0041 (GhDRF), Gh_D10G1429 (GhCHS1), Gh_D10G1431 (GhCHS1), Gh_D12G0566 (GhF3H), and Gh_D13G2458 (GhC4H) are key genes in the flavonoid metabolism pathway (Table S5).Six of the seven genes were preferentially expressed in 5 and 10 DPA fibers, with Gh_A05G3400 expressed only in 10 DPA fiber (Fig. S2). These results indicate that genes in the flavonoid biosynthesis pathway are highly expressed in the early stage of fiber development.This result is consistent with those of previous studies [12,13]. The substrate naringenin(NAR) of GhF3H retards fiber development [34], and the flavonoid metabolism pathway can combine with the lignin metabolism pathway to regulate fiber development[35].

4.2. Transcription factor associated with fiber development

Transcription factors (TF) are important during fiber development [17,36-38]. In previous studies, GhMYB109, GhMYB25, and GhMYB25-like were associated with fiber cell elongation or trichome development [17,37,38]. In cotton, upregulating the expression of the NAC gene GhFSN1 increased the secondary wall thickness of fibers but reduced fiber length [39]. TCP TFs influenced fiber development by regulating jasmonic acid biosynthesis and other pathways [40,41]. In the present study, 3 TCP genes, 3 WRKY genes, 4 NAC genes, 10 bHLH genes, 13 MYB genes, and 8 ERF genes were identified (Table S5). Of these, 2 TCP genes (Gh_D01G1783 and Gh_D11G0333), 3 bHLH genes (Gh_A13G0929, Gh_D11G2821,and Gh_D12G0666) were up-regulated from 5 DPA to 10 DPA fibers and down-regulated in 20 and 25 DPA fibers (Fig. S2). In the previous study [41], the TCP gene Gh_D11G0333 (GhTCP14) was preferentially expressed in 6-12 DPA fibers and interacted with some genes (GhSLR1, GhARF6, GhBZR1, GhEIN3, GhGL3, GhMYB23, GhMYB25, GhMYB25-like, and others) associated with fiber development. Gh_D01G1783 (GhTCP8) was coexpressed with Gh_D11G0333 (GhTCP14) and also interacted with GhSLR1, GhARF6, GhGL3, GhMYB23, GhMYB25, and GhTTG1 in yeast cells [41]. Gh_D11G0333 (GhTCP14) also increasedinitiation and elongation of trichomes and root hairs in Arabidopsis [42]. Two bHLH genes, Gh_A13G0929 and Gh_D11G2821,were homologous with AtAIF4,which regulates cell elongation in Arabidopsis [43]. Seven ERF TFs were highly expressed during all stages of fiber development, suggesting that ethylene is important in fiber development.

Table 3-The function annotations of 25 candidate genes related with fiber development.

4.3. Candidate genes located in fiber quality QTL intervals

Based on our previous studies [3,8,27], 15 stable QTL for fiber quality were identified in the present study (Table 2), and 31 DEGs were located in 9 QTL intervals. However, no DEGs were found in the intervals of qFE-c22-1,qFL-c10-1,qFL-c10-2,qFS-c17-1,and qMIC-c25,spanning 309.7,469.4,0.820,971.4,and 557.3 kb,respectively. The difficulty of identifying DEGs in these small intervals is a limitation of DGE profiling technology.Among the 31 DEGs anchored in the QTL regions, 25 genes were highly expressed at four key stages of fiber development(Table 3),and were assigned as candidate genes for fiber QTL. Of the 25 candidate genes, Gh_D07G1799 (GhGAUT9) was involved in pectin biosynthesis, as the qua1 mutant has a pectin-deficient phenotype in Arabidopsis [44]. Gh_D07G1799 (GhGAUT9) was induced by GhGalT1 and GhGalT1 was involved in control of fiber development [45,46]. Gh_D07G1974 (GhRNS2) was homologous to AtRNS2,the downstream target of AtRHS10,and AtRHS10 was a negative controller of root hair growth in Arabidopsis[47].Root hairs of AtRNS2 overexpression lines were much shorter than those of the wild type in Arabidopsis[47].Gh_D11G0323(GhATH1)was homologous to AtATH1, and stem length and cell length were inhibited in AtATH1 overexpression lines[48,49].The TCP gene Gh_D11G0333(GhTCP14)was highly expressed in 5 and 10 DPA fibers and plays critical roles in cotton fiber development,especially in fiber initiation and elongation [41,42]. These results provide candidate gene information for future fine mapping and cloning of genes influencing cotton fiber development.

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

Acknowledgments

This work was supported by the National Natural Science Foundation of China(31871669).