Phyton- International Journal of Experimental Botany

DOI: 10.32604/phyton.2021.014437

ARTICLE

Rapid Identification of a Candidate Gene Related to Fiber Strength Using a Superior Chromosome Segment Substitution Line from Gossypium hirsutum × Gossypium barbadense via Bulked Segregant RNA-Sequencing

Qi Zhang1#, Pengtao Li2#, Aiying Liu1, Shaoqi Li1, Quanwei Lu2, Qun Ge1, Junwen Li1, Wankui Gong1, Xiaoying Deng1, Haihong Shang1,3, Yuzhen Shi1,*and Youlu Yuan1,3,*

1State Key Laboratory of Cotton Biology, Key Laboratory of Biological and Genetic Breeding of Cotton, The Ministry of Agriculture, Institute of Cotton Research, Chinese Academy of Agricultural Sciences, Anyang, 455000, China
2School of Biotechnology and Food Engineering, Anyang Institute of Technology, Anyang, 455000, China
3Zhengzhou Research Base, State Key Laboratory of Cotton Biology, Zhengzhou University, Zhengzhou, 450001, China
*Corresponding Authors: Yuzhen Shi. Email: shiyuzhen@caas.cn; Youlu Yuan. Email: yuanyoulu@caas.cn
Received: 26 September 2020; Accepted: 09 December 2020
#These authors made equal contributions

1  Introduction

Cotton is one of the most important commercial crops all over the world, providing the major natural fibers for the textile industry. With the development of modern industry and people’s living standard, higher requirements have been put forward for modern textile mills to achieve better spinning performance [1]. Therefore, obtaining the high-quality fiber has been the vital goal in cotton breeding currently, yet due to genetic-basis complexities of fiber formation, it remains a challenging task for cotton researchers.

There are four cultivated species of Gossypium, namely diploid G. herbaceum and G. arboretum (2n = 2× = 26) and tetraploid G. hirsutum (Upland cotton) and G. barbadense (Sea Island cotton) (2n = 4× = 52) [2,3]. Among them, G. hirsutum and G. barbadense are the most widely cultivated cotton species and can contribute more than 95% fiber production per year, of which the former harbors the characteristics of high yield and moderate fiber performance, and the latter has the advantages of the premium fiber quality yet subjected to the disadvantage of low yield [4,5]. As a new thought in adequate utilization and integration of excellent traits derived from Upland and Sea Island cotton, Chromosome Segment Substitution Lines (CSSLs) lays the theoretical foundation of developing new cotton varieties with high yield, fine quality. The CSSLs are subjected to multi-generation backcrosses with the recurrent parent and subsequently multi-generation self-crosses, combining marker-assisted selection (MAS) to ultimately obtain the Near-isogenic lines (NILs). The CSSLs are the ideal materials for dividing the complex and variable quantitative traits or phenotypes into a single inheritance factor. Therefore, CSSLs have been widely applied to mapping quantitative traits loci (QTL) in tomato (Lycopersicon esculentum) [6], rice (Oryza sativa) [7], wheat (Triticum aestivum) [8], maize (Zea mays) [9], soybean (Glycine max) [10], and cotton [11–13].

Cotton fiber derived from ovule epidermal cells is a single-celled seed trichome. Although all the ovule epidermal cells in cotton have the potential to differentiate into fibers, only 15–25% of them can eventually develop into mature fibers [14].The development of cotton fiber can be divided into five overlapping stages: initiation (−3 to +3 DPA), elongation (3 to 23 DPA), transition, secondary wall synthesis (16 to 40 DPA), and maturation (40 to 50 DPA) [15], which collectively affect the quality traits of fiber quality, including fiber length, strength, micronaire, elongation, and uniformity [16]. With the rapid development of high-through sequencing technology and complete release of Gossypium genomes [17–23], more efficient strategies have been applied to identifying the genes responsible for fiber length and strength in cotton. The study combined re-sequencing technology and genome-wide association study (GWAS) strategy to identify Gh_A10G1256 associated with fiber length, and Gh_A07G1767, Gh_A07G1768, and Gh_A07G1769 related to fiber-strength [24]. Through fine mapping of stable QTL and quantitative real-time PCR (qRT-PCR) verification, a LRR-RLK (leucine-rich repeat protein kinase) family protein (Gh_A07G1749) was identified as a candidate gene for fiber-strength qFS07.1 [25]. Transcriptome analysis was performed on three CSSLs derived from CCRI36 (G. hirsutum) and Hai 1 (G. barbadense) during the secondary wall thickening period, and three differentially expressed genes (DEGs), namely XLOC_036333 (MNS1), XLOC_029945 (FLA8), and XLOC_075372 (snakin-1), were screened to play the significant roles in regulating fiber strength [26]. Therefore, new ways of identifying major QTL for fiber strength in cotton are still needed.

Due to bulked segregant analysis (BAS) method and RNA-seq technology, BSR (bulked segregant RNA-seq) can not only identify the DEGs by comparing the two constructed extreme pools, but also screen the SNPs (single nucleotide polymorphisms) co-isolated with mutant genes by linkage analysis [27]. Due to saving plenty of manpower, material resources, and shortening the breeding cycle, BSR-seq has attracted more attention in recent years, which has been an efficient strategy for rapid QTL mapping and efficient identification of candidate genes in connection with key traits in plants [27–30]. A candidate gene (GRMZM2G055704) responding to waterlogging was identified by BSR-Seq method on waterlogging sensitive and resistant pools in maize [31]. In order to screen the genes responsive to early and late flowering in tree peony, BSR-seq association analysis obtained 291 unigenes, of which 7 DEGs (c42942.graph_c0, c58332.graph_c0, c58361.graph_c0, c57417.graph_c0, c46352.graph_c0, c53143.graph_c0, and c58526.graph_c0) were confirmed to relate with early and late flowering [32]. But it is rare to use CSSL population with BSR method to map fiber strength related genes in cotton.

In the present study, the high-yield and superior-fiber CSSL, MBI9915, was chosen to develop the secondary segregation population BC7F2 with CCRI36. MBI995 was derived from BC5F3:5 CSSL populations constructed by G. barbadense Hai 1 as the donor parent and G. hisutum CCRI36 as the recipient and recurrent parent in the previous study [33]. In light of significant traits of fiber strength on fiber quality, we performed BSR-seq analysis on the fiber samples collected from two parents (MBI9915 and CCRI36) and two extreme pools at 20 DPA, respectively. Plenty of DEGs were identified with comparisons between different samples, DEGs subsequently underwent GO functional enrichment and KEGG pathway analysis. Through SNP-index association analysis method and qRT-PCR technology, candidate genes were finally screened to control fiber strength. These results provide a rapid and efficient method to screen the candidate genes by BSR-seq, which will facilitate the revelation of molecular mechanism of fiber formation in cotton.

2  Material and Methods

2.1 Plant Materials

Two parents, the superior CSSL MBI9915 and its recurrent parent CCRI36, as well as their segregation population BC7F2 were planted at the Anyang experimental farm of Institute of Cotton Research of Chinese Academy of Agricultural Sciences (Anyang, Henan Province) in 2016, where BC7F2:3 populations were planted through plant-to-row method in 2017. The high-yield and wide-adaptability Upland cotton CCRI36 as the recipient parent, and good-quality and disease-resistance Sea Island cotton Hai1 as the donor parent were chosen to conduct the hybridization in 2003, of which the offspring F1 was subsequently subjected to five generations of backcrosses with the recurrent parent CCRI36 and three generations of self-crossings, eventually obtaining the developed BC5F3:5 CSSLs with plant-to-row and plant-to-family method, respectively [33,34]. Based on the multi-year and multi-environment trait surveys, MBI9915 harboring the high-yield and superior-quality fiber was chosen for secondary population construction, and then underwent two-generation backcrossing with CCRI36 and one-generation self-crossing to develop BC7F2 and BC7F2:3 populations, which were planted with row length of 8 m, row spacing of 0.8 m, and plant spacing of 0.25 m.

2.2 Phenotypic Analysis

Naturally opened bolls from the BC7F2 population were harvested per plant in 2016, while 30 naturally opened bolls from BC7F2:3 population per row were harvested in 2017, of which the fiber quality traits (including fiber length, strength, micronaire, elongation, and uniformity) were tested with HFT9000 using HVICC international calibration cotton samples in the Cotton Quality Supervision and Testing Center of the Ministry of Agriculture of China.

2.3 RNA Isolation

Regarding the day tagging flowers as 0 DPA, 1 developing bolls of fiber samples were taken from MBI9915 and CCRI36 at 20 DPA were immediately immersed in ice. Comprehensive consideration on the fiber performances of both BC7F2 and BC7F2:3, 16 ones with high fiber strength and 16 ones with low fiber strength were separately selected from BC7F2:3 lines at 20 DPA to construct high and low extreme pools. After being dissected from the developing bolls, fiber sampled were quickly frozen into the liquid nitrogen, and subsequently stored at −80°C. RNA extractions were conducted on each plant, and the obtained RNA samples were equally mixed to construct the extreme pools. RNA prep Pure Plant Kit (Tiangen, Beijing, China) was used to extract high-quality RNA samples. The quality and quantity of the RNA were verified using 1% agarose gel and spectrophotometer (Nanodrop 2000, Thermo Scientific, USA).

2.4 Library Construction and Transcriptome Sequencing

First, high-quality mRNA was enriched from total RNA using Oligo (dT) magnetic beads, and was randomly cleaved into short fragments by fragmentation buffer, then double-chain cDNA was synthesized with mRNA as a template. After purification, the cDNA segments were utilized to construct the cDNA libraries through end-repair, adding the tail of the ploy A, and selecting the size of the fragments. At last, the cDNA libraries were obtained by PCR amplification, and a total of 4 RNA libraries were sequenced using Illumina HiSeqTM 4000 sequence platform (Biomarker Technologies Corporation, Beijing, China).

By high-throughput sequencing, a large amount of raw data was obtained, and stored in FASTQ format file. The raw data was filtered to obtain clean reads. Clean reads were aligned to the reference genome of TM-1 [21], and the position information of reference genome or gene and the sequence characteristic information of the sequencing samples were obtained through STAR software [35].

2.5 Analysis of Differentially Expressed Genes (DEGs)

EBSeq software was used to analyze differentially expressed genes, FPKM (fragments per kb per million of the mapped reads) [36] parameter as an index for reviewing transcripts or gene expression levels. In the detection of differentially expressed genes, Fold Change (FC) ≥ 2 and False Discovery Rate (FDR) ≤ 0.01 are identified as the DEGs screening standard, of which the former represents the ratio of expression between two parents and pools at 20 DPA, and the latter is obtained by correcting the p-value of difference significance. The accepted Benjamini-Hochberg method was used to correct the p-value obtained from the original hypothesis test. In this experiment, p-value ≤ 0.05 was selected as the criterion for screening differentially expressed genes. To understand the functional classification of all the obtained DEGs, BLAST [37] software was used to annotate the coding genes in Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases [38].

2.6 BSR-Associate Analysis

Single nucleotide polymorphism (SNP) callings were processed with GATK [39] software toolkit (http://www.broadinstitute.org/gatk/guide/best-practices.php) with strict filtering of SNPs: SNP cluster filtering, near Indel SNP filtering and neighboring InDel filtering [40]. Before doing correlation analysis, SNPs were filtered again. The criteria were as follows: Firstly, SNP loci with multiple genotypes were filtered out, and then SNP loci with read support less than 4 were filtered out, re-filter SNP loci with identical genotypes between mixing pools and SNP loci with recessive mixing pool genes not from recessive parents. Finally, high-quality SNP loci were obtained for SNP-index analysis. SNP-index method was for marker association analysis based on genotype frequency differences between high fiber strength pool and low fiber strength pool [41,42], which was mainly used to find significant differences in genotype frequencies between high fiber strength pool and low fiber strength pool using Δ(SNP-index). Δ(SNP-index) was calculated by subtracting the SNP-index of high fiber strength pool from low fiber strength pool.

2.7 Quantitative Real-Time PCR (qRT-PCR)

Five genes on chromosome 7 were chosen to conduct quantitative real-time PCR for verifying BSR-seq results. Two parents (CCRI36 and MBI9915) at 20 DPA with three technical replicates were used for qRT-PCR. The cDNAs synthesized by using a HiScript® II Q RT SuperMixfor qPCR (+gDNA wiper) kit (Vazyme Biotech, Nanjing, China). ChamQTM Universal SYBR® qPCR Master Mix (Vazyme Biotech, Nanjing, China) kit for real-time PCR was used for qRT-PCR amplification, which were performed on the ABI 7500 fast Real-Time PCR System (Applied Biosystems, USA), following the protocol of the kit. β-Actin gene (Primer sequence; F:5′-ATCCTCCGTCTTGACCTTG-3′, and R: 5′-TGTCCGTCAGGCAACTCAT-3′) was used as an internal control to normalize the genes expression. The primer sequences of five target genes are enlisted in Tab. S1. The qRT-PCR master mix was prepared as 20 μL, the process was employed according to the following steps: One cycle of 95°C for 30 s; 40 cycle of 95°C for 10 s, 60°C for 30 s; one cycle of 95°C 15s, 60°C for 60 s, 95°C for 15s. 2−ΔΔCt method was followed for obtaining the relative gene expression [43].

2.8 Gene Clone, Phylogenetic Analysis, Gene Structure Prediction and Co-Expression Analysis of Gh_A07G0837

RNAs of CSSL MBI9915 and G. barbadense Hai1 were extracted according to the RNA Isolation method. NCBI Primer-blast Web page was used to design the candidate gene’s primers (F: ATGTTCATTTCTTCTTGTGC, R: TCATAAAACTAAAGGCTCCA). TransScript All-in-one First-Strand cDNA Synthesis SuperMix for PCR kit (Transgen Biotech, Beijing, China) was used to synthesis cDNA and TransStartFastPfu Fly DNA Polymerase kit (Transgen Biotech, Beijing, China) was used for amplifying the coding sequence (CDS) of the candidate gene. PCR conditions were one cycle of 95°C for 2 min; 40 cycles of 95°C for 20 s, 60°C for 20 s, 72°C for 40 s; one cycle of 72°C for 5 min. The product was cloned into pEASY-Blunt Zero vector (Transgen Biotech, Beijing, China) for sequencing.

Alignment of Gh_A07G0837 protein sequences in CCRI36 and Hai1, and in four cultivated species were performed by DNAMAN program and NCBI protein-blast, respectively. The phylogenetic tree of proteins in four cultivated species of Gossypium, CCRI36 and Hai1 was constructed with the MEGA program (version 7.0) by the neighbor joining (NJ) method with 1000 bootstraps [44].

Using the protparam tool in the website of ExPASY, the components of protein encoded by Gh_A07G0837 in two materials were analyzed. The location of proteins are predicted by protcomp 9.0. PredictProtein and Pheyre2 were used to predict the secondary and tertiary structure of proteins, respectively. R (3.3.0) software package (https://CRAN.R-project.org/package=pheatmap) was used to conduct co-expression analysis.

3  Results

3.1 Phenotypic Data

On behalf of taking extreme-pool samples with more accurate phenotypic traits, the chosen individuals from BC7F2 population in 2016 were re-planted with plant-to-row to obtain BC7F2:3 populations in 2017, of which 30 naturally opened bolls were subjected to fiber-quality testing (Tab. 1). The average fiber-strength value of 16 individuals from high extreme pools was 39.6 in 2016, while the corresponding value of 16 individuals from low extreme pools in 2016 was 30.0; the mean fiber-strength values of 16 lines from high and low extreme pools in 2017 were 36.9 and 30.6, respectively. Referring to the fiber-strength comparisons between 2016 and 2017 materials, both higher consistency and significant difference were identified, which directly verified the reliability of sampling the individuals in 2016 for BSR-seq analysis.

Table 1: Phenotypic value of fiber strength for two extreme pools in two years

3.2 Sequence Assembling and Analysis

For rapid identification of candidate genes related to fiber strength, BSR-seq was separately conducted on CCRI36, MBI99915, high and low extreme pools, and the obtained raw reads were then subjected to filtering low-quality reads, totally generating 186.622 million (approximately 27.9G data) clean reads from two parents and 369.247 million (approximately 55.2G data) clean reads from two extreme pools (Tab. 2). Subsequently, all the clean data mapped to the reference genome [21], of which the alignment percentage of mapped reads was over 83.97%, while no less than 68.69% were uniquely mapped reads. Meanwhile, the average percentages of Q30 and GC content were 93.83% and 44.56%, respectively, which laterally proved the reliability of our BSR-seq results.

Table 2: Summary of the BSR-Seq outcomes of 4 separately samples

3.3 Analysis of Differentially Expressed Genes (DEGs)

Based on the reference genome TM-1, a large number of differentially expressed genes (Fold change ≥ 2 and FDR < 0.01) were significantly identified by EBSeq [45], to be special 3742 DEGs (793 up-regulated and 2949 down-regulated ones) between the two parents and 3252 DEGs (2913 up-regulated and 339 down-regulated ones) between the extreme pools at 20 DPA (Tab. 3), with 536 genes in common (Fig. 1A). Obviously, there were more DEGs in CCRI45-vs.-MBI9915 group (Fig. 1B) than those in High pool-vs-Low pool group (Fig. 1C), of which the former group identified over 8 times as many down-regulated DEGs, while the latter group identified more than 4 times as many up-regulated DEGs, implying the wide variations between the two groups. Having performing the hierarchical clustering analysis on DEGs from two pairwise-comparison groups, most genes showed diverse expression levels either in two parent groups (Fig. 2A) or in two extreme pools (Fig. 2B), while similar phenomenon occurred on the slight DEGs with high-ratio expression qualities between two groups, which might be closely related to the mechanism and regulation of fiber strength.

GO enrichment analysis was subsequently conducted on the whole DEGs from two parent groups (Fig. 3), two extreme pools (Fig. 3), and common DEGs (Fig. 3), respectively, aiming at identifying the putative functional genes or significant signal pathways during fiber development. Top 20 GO terms were separately enriched into three main categories, and as for the DEGs from both two parents and two extreme pool groups or the common DEGs, cell part and cell were the predominant subcategories in cellular component, catalytic activity and binding were the top subcategories in molecular functional category, and metabolic process, cellular process, and single-organism process were the most enriched subcategories in biological process. In the process of fiber strength enhancement, a lot of secondary wall were come into being, resulting in good fiber strength. DEGs were annotated into these cellular component, molecular functional category and biological process, proved that.

Similarly, the whole DEGs were subjected to KEGG pathway analysis and were mapped onto five top-level subcategories of cellular processing, environmental information processing, genetic information processing, metabolism and organismal systems. At 20 DPA, DEGs in two parent groups (Fig. 4A) were mostly annotated into plant hormone signal transduction (ko04075), carbon metabolism (ko01200), while plant hormone signal transduction (ko04075) and starch and sucrose metabolism (ko00500) were the top two pathway subcategories for DEGs in two pool groups (Fig. 4B). The common DEGs both in two groups were mostly enriched to plant hormone signal transduction (ko04075) pathway subcategories (Fig. 4C). Auxin [36,37], gibberellins [38], brassinosteroid [39], ethylene [40,41], abscisic acid [42,43] and cytokinin [21,44] play important roles in the development of fiber. The results indicated that these DEGs may play important roles in secondary wall thickening and result in fiber enhancement.

Figure 1: Overlapped DEGs in CCRI36-vs.-MBI9915 and High pool-vs.-Low pool. Volcano plot analysis of differentially expressed genes (DEGs) for the pairwise comparisons CCRI36-vs.-MBI9915 (a) and (b) High pool-vs.-Low pool

Figure 2: Hierarchical cluster analysis of differentially expressed genes (DEGs) for the pairwise comparisons CCRI36-vs.-MBI9915 (A) and High pool-vs.-Low pool (B). Note: The abscissa represents the name of the sample and the clustering results of the sample, while the ordinate represents the clustering results of the genes. Different columns in the figure represent different samples, and different rows represent different genes. The color represents the level of gene expression (log2FPKM + 1) in the samples

Table 3: The number of differentially expressed genes in 4 samples at 20 DPA

3.4 qRT-PCR Verification of BSR-Seq Results

To confirm the reliability of BSR-seq results, qRT-PCR was performed on the five DEGs from two parents (CCRI36 and MBI9915) at 20 DPA, namely as Gh_A07G0831, Gh_A07G0837, Gh_A07G0838, Gh_A07G0858 and Gh_A07G0973. According to the BSR-seq data, all the five DGEs showed significantly lower expression levels in MBI9915 than those in CCRI36, which presented the similar difference-pattern between two parents in qRT-PCR results Fig. 5. The relative expression of Gh_A07G0831, Gh_A07G0837, Gh_A07G0838, Gh_A07G0858 and Gh_A07G0973 in CCRI36 were significantly higher than that in MBI9915. The qRT-PCR data were consistent with the BSR-seq results, fully demonstrating the reliability of this study.

3.5 Chromosome Region Confirmation and Candidate Gene Identification

To map candidate genes related to fiber strength, BSR-seq analyses were conducted on two extreme pools at 20 DPA, and 16 samples with high fiber strength and 16 samples with low fiber strength were collected and equally mixed for RNA sequencing, respectively. After being trimmed, the clean reads were aligned to the TM-1 reference genome for SNP identifications between high and low fiber-strength pools, resulting in a total of 831 high-probability SNPs via bulked segregant analysis. All the identified SNPs were quantified with SNP-index, of which the information were synthetically compared between high and low fiber-strength pools to calculate and plot the Δ(SNP-index) against the TM-1 reference genome positions (Supplementary Fig. 1). In this study, associated chromosomal regions were identified with |Δ(SNP-index)| > 0.6365 as the threshold value.

Totally, 18 associated chromosomal regions located in 15 chromosomes (8 ones in A sub-genome and 7 ones in D sub-genome) were identified, and there were two different chromosomal regions in chromosome A05, A07, and A09, respectively (Tab. S2). Referred to the TM-1 genome information, a total of 1981 genes were annotated in 18 associated regions, including 141 ones with non-synonymous mutation, which were generally regarded the potential candidate genes involved in fiber strength formation. Meanwhile, GO enrichment and KEGG pathway analyses were conducted on the whole annotated genes, and the most genes were enriched to cell part and cell, catalytic activity and binding, and cellular process and single-organism process in cellular component, molecular function, and biological process based on GO categories (Fig. 6A), while ribosome (ko03010), carbon metabolism (ko01200), biosynthesis of amino acids (ko01230) and protein processing in endoplasmic reticulum (ko04141) were top 4 pathways with the greatest number of genes (Fig. 6B).

Figure 3: Go database annotation of differentially expressed genes (DEGs) between two parents, pools and overlapped at 20 DPA

Figure 4: KEGG pathway classification of differentially expressed genes (DEGs) between two parents and pools in 20 DPA. Note: The χ-axis indicates the number of genes annotated to this pathway and the number of genes accounted for the total number of annotated genes. The γ-axis is the name of KEGG metabolic pathway

Figure 5: The relative expression of five selected genes in two parents. The χ-axis represents two parents in 20 DPA, T05 (CCRI36) and T06 (MBI9915). The γ-axis indicates relative expression and the value are expressed as the means of three replicates ± SD. All genes were select based on RNA-Seq results, and genes’ relative expression was validated by qRT-PCR

Figure 6: Go database annotation (A) and KEGG pathway classification (B) of genes in the interval through association analysis between two pools at 20 DPA. Note: The χ-axis indicates the number of genes annotated to this pathway and the number of genes accounted for the total number of annotated genes. The γ-axis is the name of KEGG metabolic pathway

Among all the 1981 genes, 69 genes and 93 genes were differentially expressed between two parents and between two extreme pools, respectively, while only 12 common DEGs were identified between two groups. According to annotation information of function genes shown in Tab. 4, 2 genes might play the significant roles during fiber development, namely as Gh_A05G1089 and Gh_A07G0837, which were separately annotated as abscisic acid transport in biological process and Golgi apparatus in cell component. Abscisic acid has been reported as a signal of secondary wall thickening to make great contribution to fiber development, and Golgi apparatus were similarly of importance for plant growth [46–48]. Interestingly, the region of chromosome 7 associated with the correlation analysis was also consistent with the previous QTL localization region, of which segregation populations were constructed by MBI9915 and CCRI36. A associated region with 6.74 Mb on chromosome A07 (chr7: 13.7–20.4 Mb) was overlapped (chrA07:14-16 Mb) [12], therefore the only overlapped DEG Gh_A07G0837 was deemed as a candidate gene related to fiber strength.

Table 4: Functional annotation of 12 important genes

3.6 Clone, Phylogenetic Analysis and Structural Prediction on Gh_A07G0837

For further confirming the contribution of the candidate gene during fiber development, gene clone was conducted on the Gh_A07G0837 in CCRI36 and Hai1, respectively. The coding sequence (CDS) length was 807 bp both in CCRI36 and Hai1, while 2 SNPs were identified to separately locate in the 30th base and the 486th base (Fig. 7A), of which the former resulted in non-synonymous mutation, specifically encoding the aspartic acid in CCRI36 but the glutamic acid in Hai1, while the latter caused synonymous mutation. The CDS length of the candidate gene was only 627 bp in TM-1 reference genome, and the short 181 bp might be the second intron sequence based on annotation information, indicating the probable occurrence of alternative splicing (AS) event in CCRI36 and Hai1. Subsequently, the phylogenetic relationships of Gh_A07G0837 protein sequences were analyzed among the CCRI36, Hai1 and other cotton species, and BLAST results showed the deduced amino acid sequence of Gh_A07G0837 in CCRI36 were highly similar to G. barbadense (99% identity and 95% similarity), while 81% identity and 97%, 95% and 95% similarity to G. hirsutum, G. raimondii and G. arboretum. However, for the deduced amino acid sequence of Gh_A07G0837 in Hai1, G. barbadense (98% identity and 95% similarity) also had high similarity, G. hirsutum, G. raimondii and G. arboretum was 80% identity and 97%, 95%, and 95% similarities (Fig. 7B).

A non-synonymous mutation of the gene in CCRI36 and MBI9915 resulted in the transformation of aspartate at the 10th position of amino acid sequence into glusate. Instability index of protein in CCRI36 was lower than that in MBI9915. The others were listed in Tab. S3. Two proteins were all most likely in intracellular part and were likely in peroxisomal and endoplasmic reticulum. The secondary prediction of protein in CCRI36 contained Alpha helix (30.36%), Beta folding (25.00%) and irregular curl (44.64%), while protein in MBI9915 contained Alpha helix (30.36%), Beta folding (23.66%) and irregular curl (45.98%). The tertiary structure of proteins in CCRI36 and MBI9915 were predicted as oxidoreductase, whose imagines were shown in Supplementary Fig. 2.

Figure 7: Analysis of the Gh_A07G0837 sequence in CCRI36 and Hai1. A. Alignment of the Coding sequences of Gh_A07G0837 in CCRI36 and Hai 1 by DNAMAN software. B: The phylogenetic relationship between the protein sequence in CCRI36, Hai1 and other proteins in some cotton species. Sequences were aligned by ClustalX and analyzed using the MEGA 6.0. Phylogenetic tree was constructed by the Neighbor-joining method with 1000 bootstrap replicates

3.7 Co-Expression of Gh_A07G0837

R software was adopted to conduct a co-expression analysis on Gh_A07G0837. Totally, there were 685 co-expression genes of Gh_A07G0837, of which 144 genes were also differentially expressed between two parents. Most genes were annotated into “plant hormone signal transduction (ko04075)”, and “Biosynthesis of secondary metabolites (ko01110)” pathway in KEGG database (Fig. 8). Gh_A10G0249 and Gh_A11G3069 were annotated to “protein processing in endoplasmic reticulum” pathway (ko04141). Gh_A08G1031 was annotated into “starch and sucrose metabolism” pathway (ko00500). Gh_D02G0430 and Gh_D07G2383 were annotated into “plant hormone signal transduction” pathway(ko04075). Gh_D05G0154 was annotated into “phenylpropanoid biosynthesis” pathway (ko00940).

Figure 8: KEGG functional annotation of 144 co-expression genes of Gh_A07G0837

4  Discussion

4.1 CSSLs Were an Ideal Material to BSA and RNA-Seq

Cotton (G. hirsutum L.), one of the most important cash crops, produces the vast majority of natural fiber, which is widely cultivated all over the world. With the development of modern textile industry and the improvement of human living standard, high yield has not been the only goal in the cotton breeding, while superior fiber quality has increasingly attracted attentions. However, either Upland cotton or Sea Island cotton, as the two most widely cultivated varieties, could not simultaneously meet the above-mentioned demands, of which the former subjects to the ordinary fiber performance and narrow genotypic milieu, and the latter has the disadvantages of poor yield and low adaptability [49,50]. An effective and economical method has been put forward to develop chromosome segment substitution lines (CSSLs) in cotton production, utilizing the conventional breeding (hybridization, backcross and self-crossing) and maker-assisted selection (MAS) technologies to cultivate novel cotton varieties harboring high yield and superior fiber quality. As to fiber strength, one of the critical traits determining cotton quality, its formation principally depends on the secondary wall thickening period, and in-depth studies on candidate gene identification will greatly contribute to the molecular-mechanism revelation of fiber development.

Based on the multiple-year and multiple-environment trait surveys of developed CSSLs BC5F3:5 constructed by G. hirsutum CCRI36 and G. barbadense Hai1, superior MBI9915 was chosen to develop the segregation population F2 and F2:3 together with its receptor and recurrent parent CCRI36, and 41 polymorphic SSR markers screened from the high-density linkage map of CCRI36 × Hai1 [5] were utilized to conduct QTL mapping, resulting in 24 QTLs related with cotton yield and fiber quality [12]. Meanwhile, transcriptome sequencing was performed on superior MBI9915 and ordinary CCRI36 during fiber development (at 10, 15, 20, 25, and 28 DPA), and 1801 identified DEGs were enriched into cell wall organization and response to oxidative stress and auxin pathways [33]. The above-mentioned studies have laid the solid foundations for candidate gene identification related to fiber strength.

4.2 Candidate Genes Related to Fiber Strength

Combining bulked segregant analysis (BAS) method and RNA-seq technologies, BSR-seq could analyze the DEGs located in associated chromosomal regions by comparing two extreme pools, finally rapid identification of the candidate genes relevant to key quantitative traits. In this study, the CSSL MBI9915 and its recurrent parent CCRI36 were used to develop BC7F2 and BC7F2:3 populations, which were subjected to fiber quality tests. From BC7F2 population, 16 samples with high fiber strength and 16 samples with low fiber strength were collected at 20 DPA to construct two extreme pool libraries, and then were separately sequenced together with transcriptome libraries of their parents (MBI9915 and CCRI36).

Based on BSA analysis, 18 regions related to fiber strength were identified. One of regions were same to previous studies. Importantly, another 17 regions were first time to be identified related to fiber strength. Among 18 regions, 12 genes were differently expressed in both two parents and pools. The results indicated that 12 DEGs were likely to be related to fiber strength. Two candidate genes were finally screened based on functional annotation of reference genome, namely as namely as Gh_A05G1089 (abscisic acid transport) and Gh_A07G0837 (Golgi apparatus), of which the former has been reported to make great contribution to fiber development as a signal of secondary wall thickening [47,48], while the latter was involved in fiber development [46]. What is interesting, the region of chromosome 7 associated with the correlation analysis was also consistent with the previous QTL localization region, and the associated region with 6.74 Mb on chromosome A07 (chr7: 13.7–20.4Mb) was overlapped (chrA07:14-16Mb) [12], indicating the only overlapped DEGs Gh_A07G0837 was deemed as a candidate gene related to fiber strength. Subsequently, gene clone and sequence analyses were conducted on the candidate gene, and CCRI36 and Hai1 had two SNPs difference and one amino acid difference that induced the phenotypic difference between two parents. Using NCBI protein-blast, protein sequence in CCRI36 and Hai1 was highly similar with a hypothetical protein GOBAR AA28214 in G. barbadense (98% identity and 95% similarity). For nonsynonymous mutation, whether primary structure or secondary structure, it is different from CCRI36 and MBI9915. It may result in Changes of protein phenotype. In Arabidopsis thaliana (AT4G13010.1), it was one of zinc-binding dehydrogenase family protein, named Chloroplast Envelope Quinone Oxidoreductase Homolog (CEQORH) [51]. CeQORH was associated with the inner membrane of the chloroplast envelope and imported into the organelle by an alternative import pathway [52,53]. Under oxidative stress conditions, CeQORH exhibited the activity of a NADPH-dependent to eliminate highly reactive a, b-unsaturated carbonyls [51]. During the development of fiber strength, cells need to breathe to provide energy, and a series of redox reactions occurred in cell respiration. Therefore, this candidate gene may play quite an important role in the development of cotton fiber strength.

4.3 Co-Expression Analysis of Gh_A07G0837

Having performing the co-expression analysis on candidate genes, many genes were annotated into “plant hormone signal transduction”, “Biosynthesis of secondary metabolites”, which is vital to the development of fiber strength [33]. Gh_A07G0837 cooperates with those genes to regulate the development of fiber strength. Gh_A10G0249 (E3 ubiquitin-protein ligase RMA1H1) is one of Gh_A07G0837 co-expression genes, which was reported that RMA1 cooperates with Gp78 to participate in ER-associated degradation of the cystic fibrosis transmembrane conductance regulator [54]. Surprisingly, Gh_A10G0249 was identified as a down-regulated gene between two parents, just as the same expression pattern as Gh_A07G0837. We guess that it is because this gene is expressed low in MBI9915 to form good fiber strength. Another co-expressed DEG, Gh_A08G1031, also showed was differentially expressed between two parents with down-regulated tend, which was annotated as sucrose synthase 2. Sucrose synthase and invertase are the major enzymes that catalyze sucrose conversion [55], which are highly expressed during the formation of fiber length. The DEG Gh_D02G0430, namely as ethylene-responsive transcription factor 1B, is a down-regulated co-expressed gene between two parents. Ethylene is an important factor participating in plant hormone signal transduction for fiber length, and their low expression in MBI9915 might indicate that inhibiting the expression of these genes is of great significance for the formation of materials with good fiber strength. The down-regulated DEG Gh_D07G2383 with the annotation information of a Protein phosphatase 2C 77, is inhibited by ABA. Gh_D05G0154, a Peroxidase, is expressed highly in fast-elongating cotton fiber cells [56]. Our RNA-Seq showed that the expression level of Gh_D05G0154 is relatively low, especially in MBI9915, and our co-expression analysis on Gh_A07G0837 indicated that those interactional genes might collaborate in the control of regulatory fiber strength function.

5  Conclusion

In this study, we used Chromosome Segment Substitution Lines through BSR-seq technique to map gene related to fiber strength. Association analysis was conducted in this study to obtain some regions related to fiber strength, and 12 DEGs were identified with combination of RNA-seq data. Referred to the relevant QTL results, one candidate gene was screened and analyzed to participated in fiber strength, whose function will be further verified with gene editing and genetic transformation experiments. In this region, Gh_A07G0837 is differentially expressed not only in two parents but also in two pools, which is considered as a candidate gene. The sequence and structure of Gh_A07G0837 in CCRI36 is different from in MBI9915. In previous study, CeQORH(Gh_A07G0837) participated in redox reaction. We deduce that CeQORH is likely through redox reaction to cause the difference of fiber strength in CCRI36 and MBI9915.

Author Contribution: YLY, YZS and QZ conceived and designed the study. QZ and PTL drafted the manuscript. AYL, QWL, SQL, QG, JWL, and WKG participated in segregation population construction and fiber sample collections. XYD, HHS participated in BSR-seq and data analysis. All the authors read and approved the final manuscript.

Funding Statement: This study was supported by the National Natural Science Foundation of China (31801404 and 31621005) and Joint Funds of the National Natural Science Foundation (U1804103), Science and Technology development Project of Henan Province (192102110127 and 182102410041), and the Project of Director (1610162020040904).

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

Appendix

Figure S1: Putative SNPs associated with fiber strength. The polymorphic SNPs were identified by comparing RNA sequence of two pools in 20 DPA. χ-axis represents the physical position of each SNP marker and γ-axis is the linkage probability. Red dotted line indicates the threshold of linkage probability

Figure S2: Prediction of the tertiary structure of Gh_A07G0837 encoding protein

Table S1: Primers are used for qRT-PCR verification

Table S2: Information on 18 linkage domains

Table S3: Analysis of physical and chemical properties of Gh_A07G0837 in CCRI36 and MBI9915

This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.