A global survey of the gene network and key genes for oil accumulation in cultivated tetraploid cottons.

To enrich our knowledge about gene network of fatty acid biosynthesis in cottonseed, we conducted comparative transcriptome to reveal the differences in gene expression between Gossypium hirsutum and Gossypium barbadense during cottonseed development. The prolonged expression period and increased expression abundance of oil-related genes are the main reasons for producing high seed oil content (SOC) in G. barbadense, which manifested as the bias of homeologous gene expression in Dt-subgenome after 25 day postanthesis (DPA). The dynamic expression profile showed that SAD6 and FATA are more important for oil biosynthesis in G. barbadense than that in G. hirsutum. Three key transcription factors, WRI1, NF-YB6 and DPBF2, showed their elite roles in regulating seed oil in cotton. We observed that sequence variations in the promoter region of BCCP2 genes might contribute to its divergence in expression level between the two species. Based on the quantitative trait loci (QTL) information of the seed oil content and utilizing additional G. barbadense introgression lines (ILs), we propose 21 candidate genes on the basis of their differential expression level, of which the GbSWEET and the GbACBP6 showed the potential functional to improve the oil content. Taken together, studying the different expression of oil-related genes and their genetic regulation mechanisms between G. hirsutum and G. barbadense provide new insights to understanding the mechanism of fatty acid biosynthesis network and fatty acid genetic improving breeding in cotton.

Introduction

Worldwide increasing demand for vegetable oil, especially those with high nutritional profile, requires continuous breeding efforts to broaden the sources and increase the yield of oil crop. Cotton (Gossypium), as an important cash crop in the world, provides not only the main natural fibre for the textile industry, but also the 2nd major source of vegetable oil (Sharif et al., 2019). Cottonseed oil contains a large amount of unsaturated fatty acids, of which more than 50% belongs to linoleic acid (C18:2), and also rich in palmitic acid (C16:0) and oleic acid (C18:1) (Konuskan et al., 2017). The higher proportion of polyunsaturated fatty acids in cottonseed oil is beneficial to human health. In addition, the suitable length of carbon chain in cottonseed oil is considered to be an ideal raw material of biofuel (Yesilyurt and Aydin, 2020). Modern upland cotton (G. hirsutum) cultivars are the result of long‐term selection and domestication (Grover et al., 2015; Senchina et al., 2003), which contain great variability of cottonseed oil content (SOC) ranging from 12.60% to 39.33% (Hinze et al., 2015; Zhao et al., 2019). However, low genetic diversity limits the potential for genetic breeding within species to improve SOC (Kaur et al., 2017; Tyagi et al., 2014). The G. barbadense has often been used as an important germplasm donor to improve yield (Li et al., 2016), fibre quality (Jenkins et al., 2018) and disease resistance in G. hirsutum (Zhao et al., 2018a).

In cotton, the interspecific introgression line (IL) population of G. hirsutum × G. barbadense has frequently been used to identify genomic regions associated with interested quantitative traits (Zhang et al., 2014). However, few studies to date have been aimed at using the G. barbadense to improve the SOC in G. hirsutum background. Studying the correlations between cottonseed traits, Wu et al. (2009) and Yu et al. (2012) showed the potential value of G. barbadense in improving SOC in G. hirsutum background. In addition, the survey of seed traits in cotton germplasm resources performed by Hinze et al. (2015) and Shockey et al. (2016) also confirmed that the G. barbadense accessions have the higher SOC characteristic. Hence, studying the biosynthesis in G. barbadense and revealing its genetic regulation mechanism represent an interesting goal in cotton oil breeding.

Lipids and fatty acids (FAs) are a large family in plant, while most edible vegetable oil consists of a few common FAs, including saturated, monounsaturated, and polyunsaturated fatty acids, and stored in seeds in form of triacylglycerol (TAG) on usual (Troncosoponce et al., 2011). The TAG biosynthesis pathway has become one of the hallmarks of lipid biochemistry, which controlled by a series of functional genes (Bates et al., 2013). As next‐generation high‐throughput sequencing has become routine, several studies in oil crop shed new light on understanding the TAG biochemical pathways involved in oilseeds at genomic and transcriptome levels (Gusev et al., 2016). De novo transcriptome sequencing analysis in Arabidopsis (Belmonte et al., 2013), rapeseed (Deng et al., 2015; Wan et al., 2017), soybean (Jones and Vodkin, 2013), camelina (Abdullah et al., 2016), peanut (Liu et al., 2019) and sesame (Wang et al., 2014) not only depicts the expression patterns of lipid genes during the development of oilseeds, but reveals the molecular mechanisms of TAG metabolism and highlights the key enzymes for lipid storage with different oil components (Baud and Lepiniec, 2010). Genetic engineering of those rate‐limiting enzymes involved in TAG has been successful in regulating the levels and composition of FAs in upland cotton. RNA interference (RNAi) of D12 fatty acid desaturase (FAD2), D9 stearoyl‐ACP desaturase (SAD1) and b‐ketoacyl‐acyl carrier protein synthase (KASII) in cotton dramatically increased the accumulation of oleic acids, stearic acids and palmitic acid in cotton seed, respectively (Liu et al., 2002, 2017b). In addition, RNAi down‐regulation of phosphoenolpyruvate carboxyl (PEPC2) involved in tricarboxylic acid (TCA) cycle in cottonseeds up‐regulated most lipid synthesis‐related genes and increased 7.3% of oil content in seed kernel (Zhao et al., 2018b). The global transcriptome profile of development seeds further enriched our understanding of the important role of these functional genes in regulating lipids biosynthesis in cotton (Hovav et al., 2015; Zhao et al., 2018c). In forward genetic analysis, a large number of oil‐related genes also have been identified by genome‐wide association study, such as fatty acyl‐ACP thioesterase B (FATB), acyl carrier protein 5 (ACP5) (Yuan et al., 2018) and KAS III (Du et al., 2018). Additionally, the calcium‐dependent lipid‐binding (CaLB) family protein (Zhao et al., 2019) and the peroxidase gene (PRXR1) (Ma et al., 2019) were also considered to be involved in improving SOC in cotton, which implied that there may be more complex of the genetic regulation of SOC in cotton.

However, previous studies mainly focused on the upland cotton, the biosynthesis of lipids and their regulatory mechanisms associated with high SOC formation in G. barbadense remain unclear. Elucidation of high‐quality cotton genome (G. hirsutum and G. barbadense) (Wang et al., 2019), and the availability of IL population derived from interspecific crossing allowed the unprecedented opportunity to study the lipids biosynthesis to understand its genetic mechanism of higher SOC in G. barbadense. In our previous research (Zhu et al., 2020), via analysing the genetic effects of interspecific IL population, we found the amazing performance of SOC and identified a lot of related QTL. Therefore, the aim of present study was to explore the mechanism of higher SOC in G. barbadense by analysing the differentially expressed genes during seed development and to propose associated candidate genes using interspecific IL population.

Results

Phenotype variation and transcriptome assembly

Based on the genotype of IL population identified in our previous research (Zhu et al., 2020), we selected the high SOC lines with inconsistent introgression fragments. The phenotype of SOC in these lines showed significantly differences from their recipient parent G. hirsutum cv. ‘Emian22’(P < 0.001), but similar to the donor parent G. barbadense acc. 3‐79 (3‐79) (Figure S1), which indicated that there may be multiple genetic loci controlling high SOC in G. barbadense.

In order to understand the underlying molecular differences of SOC formation between G. hirsutum and G. barbadense, we selected three stages of whole seeds for transcriptome sequencing to obtain a global view of gene expression during seed development. About 9 Gb transcriptome data were generated from each library, and at least 53 million paired‐end reads were obtained (Table S1). Clean reads from Emian22 and 3‐79 were mapped to the G. hirsutum and G. barbadense reference genome, respectively. Nine libraries of gene expression data from different stages were clustered together concordantly in the principal component analysis (PCA), which showed that these data were appropriated for comparative analysis in developing seed (Figure S2a). A total of 26 528 and 28 306 genes were expressed in G. hirsutum and G. barbadense during seed development, respectively, accounting for 37.79% and 39.70 of total genes in each genome (Figure S2b). Notably, the number of expressed genes decreased during seed development in both species, but the time‐specifically expressed genes in G. hirsutum was significantly higher than that in G. barbadense at 30 DPA (1361 vs 789) (Figure S2c).

Estimation of differentially expressed genes (DEGs)

The time‐series differential expression analysis helps us to understand how the different gene expression levels affect embryo development. During seed development, a total of 25 074 and 22 390 genes were differentially expressed in G. hirsutum and G. barbadense, respectively. There was no significant difference in the up‐ or down‐regulated genes between the two species in the comparison of ‘20 DPA vs 10 DPA’. Notably, the huge difference was observed in the ‘30 DPA vs 20 DPA’ comparison, only 673 DEGs were up‐regulated in G. barbadense comparing 3082 DEGs in G. hirsutum, and 1,860 DEGs were down‐regulated in G. barbadense comparing 4885 DEGs in G. hirsutum (Figure 1a). This indicated that a large number of genes may still had higher expression levels from 20 DPA to 30 DPA in G. barbadense compared with G. hirsutum, which is consistent with the result that the number of up‐regulated genes in G. barbadense is more than that in G. hirsutum in the comparison of ‘30 DPA vs 10 DPA’. A previous study reported that the oil content in cottonseed rapidly increased from 20 DPA to 30 DPA (Hovav et al., 2015), so we analysed the specifical DEGs in the comparison ‘30 DPA vs 20 DPA’. We found that only 352 DEGs and 766 DEGs were specifically up‐ and down‐regulated in the G. barbadense from 20 DPA to 30 DPA. However, a total of 3775 DEGs with 1662 DEGs up‐regulated and 2,113 DEGs down‐regulated were specifically detected in G. hirsutum from 20 DPA to 30 DPA (Figure 1b). Down‐regulated expression of a large number of genes during oil accumulation may result in the low SOC in the G. hirsutum.

Figure 1

Differentially expressed genes (DEGs) at different cottonseed developmental stages. (a) Summary of DEGs between the comparison treatments; (b) Venn diagram of specific DEGs in the treatments.

To comprehend the seed development processes, Gene Ontology (GO) categories were performed to classify the functional interpretation and gene product attributes of these DEGs. Under the biological processes, most of the DEGs were concentrated on the similar potential functions between species during seed development. A total of 1104 and 1118 genes were annotated to the lipid metabolic process (GO: 0006629) in the G. hirsutum and G. barbadense reference genome, respectively. However, no significantly enrichment was observed that the DEGs were involved in lipid‐related biological process, suggesting that the high oil content in cottonseed was regulated by a few key genes (Tables S2 and S3). Furthermore, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was used to categorize function annotations for all DEGs and summarized in Table S4 and Table S5. During the time phase from 10 DPA to 20 DPA, 91 DEGs classified in the lipid metabolism were up‐regulated in G. barbadense, and only 61 up‐regulated DEGs were detected in G. hirsutum (Figure S3). Notably, the types of lipid metabolism pathway are not different with the number of expressed genes increased (Table S5). Interestingly, in the group of down‐regulated DEGs in the ‘30 DPA vs 20 DPA’, more lipid‐related DEGs down‐regulated in G. hirsutum than that in G. barbadense (79 vs 16) (Figure S4). Among those, we found four significant enrichment pathways involved in lipid metabolism in G. hirsutum, including fatty acid biosynthesis, fatty acid metabolism, glycerolipid metabolism and biosynthesis of unsaturated fatty acids, while no relevant pathway was observed in G. barbadense (Figure S5). The prolonged expression period of lipid‐related genes may played critical roles in determining the high SOC in G. barbadense.

To further categorize the DEGs of these oil‐related genes, a hierarchical clustering analysis was performed on the expression value fragments per kilobase of exon per million mapped reads (FPKM) and visualized in a heatmap. Indeed, functional enrichment revealed that most of genes enriched with fatty acid elongation pathway and fatty acid biosynthesis processes continued to be highly expression from 20 DPA to 30 DPA in G. barbadense (Figure S6). In the meantime, the alpha‐linolenic acid metabolism‐ and linolenic acid metabolism‐related genes were highly activated at late stage in both species, which emphasized the prevalence of the high unsaturated fatty acid synthesis content in cottonseed.

Dynamic changes of expression profiling during oil accumulation

To investigate the potential mechanism underlying the variation in SOC, we compared the oil‐related genes expression during seed development between the two species, with more attention being paid to the TAG biosynthesis. According to the previous report (You et al., 2016), a total of 374 genes involved in fatty acid metabolism were identified in both species (Tables S6 and S7). Yet, only 29.68% and 28.17% of them were highly expressed (FPKM ≥ 10 in at least one stage) during seed development in G. hirsutum and G. barbadense, respectively. And among these genes, most of highly expressed genes are on orthologous positions in the genome between the two species.

Based on the reference oil biosynthesis model of upland cotton (Hovav et al., 2015; Zhao et al., 2018c), the dynamic expression profile was modified and is shown in Figure 2. Overall, the results indicated higher regulated levels of oil‐related genes involved in the entire FA biosynthesis in G. barbadense compared with G. hirsutum. In the process of de novo formation acyl chains in plastids, the activity of the subunit (α‐PHD) of the pyruvate dehydrogenase complex (PDHC) and 2‐oxoacid dehydrogenase acyltransferase (DHLAT) at 30 DPA may provide more acetyl‐CoA precursor for continuous FA synthesis. Intriguingly, as a key regulatory step of the formation of malonyl‐CoA from acetyl‐CoA, the biotin carboxyl carrier protein 2 (BCCP2) displayed a more important role in G. barbadense, while the BCCP1 showed a slightly higher expression level in G. hirsutum, like to what has been reported in the upland cotton with high oil content (Zhao et al., 2018c). The three types of ketoacyl‐AcP synthase (KASI, KASII and KASIII) were poorly expressed in cottonseed except KASI, which was abundant at 10 DPA and 20 DPA. Both enoyl‐ACP reductase (ENR) and polyketide synthase (ER) showed highly expressed in G. barbadense than that in G. hirsutum in the cycling of fatty acid synthase (FAS). As a crucial desaturation enzyme for C18:1 synthesis, the SAD showed a constantly high expression level in developing cottonseed in both species. Surprisingly, the SAD6 may played a decisive role in determining high oil content in G. barbadense with an increased in expression after 20 DPA.

Figure 2

Gene expression network for fatty acid biosynthesis in development cottonseed. The expression abundance (FPKM in log2(FPKM + 1)) is highlighted with gradient colour in 10, 20 and 30 DPA. The upper right and lower left corners represent G. hirsutum (Emian22) and G. barbadense (3‐79), respectively.

The significant differential expression pattern of the TAG biosynthesis was observed in the endoplasmic reticulum (ER), with most enzymes showing low expression levels in G. barbadense in comparison with G. hirsutum. For example, four enzymes involving in the assembly of TAG from glycerol‐3‐phosphate (G3P) and acyl‐CoAs in G. hirsutum, including glycerol‐3‐phosphate acyltransferase (GPAT), lysophosphatidic acid acyltransferase (LPAAT), diacylglycerol acyltransferase (DGAT) and phospholipid diacylglycerol acyltransferase (PDAT), displayed a slightly higher expression level during seed development, excepting PDAT, which was barely expressed in G. barbadense. The decisive enzyme of formation linoleic acid from oleic acid, FAD2 showed extensive expression abundance at 20 PDA and 30 DPA, with higher in G. barbadense. Finally, a series of enzymes related to TAG storage showed high activity at 20 DPA and 30 DPA in both species, which mirrored oil accumulation with temporal in the seed development.

A series of the transcription factors involved in regulation of seed development and maturation have been reported (Santosmendoza et al., 2008). A total of 61 and 59 genes belonging to 14 transcription factors (TFs) with their expression pattern were picked out from G. hirsutum and G. barbadense, respectively (Table S8). The time‐series expression levels of 13 expressed TFs were higher in G. barbadense than in G. hirsutum, except ABI3 and DOF2.1 (Figure 3a). Surprisingly, three TFs showed high and exclusive activity in development seed, especially at later stage in G. barbadense, including WRI1, NF‐YB6 and DPBF2. WRI1, an AP2‐type transcription factor, which was considered to be the most important regulator in seed storage metabolism, could regulate many enzymes in the FA metabolism (Maeo et al., 2009). Gene expression bias is prevalent in the polyploidy genome, and previous study has shown that about 20% of genes showed subgenome‐specific expression in upland cotton (Hovav et al., 2015). The expression pattern observed for WRI1 exhibited homeologous expression bias between the two species, although most of them in genome were barely expressed in general. The expression levels of two WRI1 genes located on chromosomes A10 and D10 in G. barbadense rapidly increased after 10 DPA; however, increased expression in G. hirsutum was only observed for Dt allele and decreased rapidly after 20 DPA (Figure 3b). In addition, the similar expression trend was detected for another two TFs (NF‐YB6 and DPBF2). The prolonged expression period and asymmetric homologous expression of these TFs may contributed to regulate the continuous activity of enzymes involving in FA biosynthesis to improve the oil content in G. barbadense.

Figure 3

Gene expression trend of oil‐related TFs in cottonseed development. (a) Heatmap showing the expression of oil‐related TFs during seed development (10 DPA, 20 DPA and 30 DPA); (b) homolog expression bias analysis of three key TFs in developing seed.

Genetic variation of oil‐related homologous TF genes

In this study, significant differential expression levels and prolonged expression time of the oil‐related homologous genes between G. hirsutum and G. barbadense were observed, which may be the direct cause of the difference SOC between the two species. In general, driving the expression of a gene is considered to be a function of its promoter, especially the core promoter region closely associated with the transcription start site (TSS) (usually less than 1 kb), which is enriched in cis‐elements/motifs for the binding of regulatory TFs (Ye et al., 2018). Therefore, comparing the variation of the core promoter region of homologous genes can quickly reveal the mechanism of differences in their expression levels. In tetraploid cotton, there are two (Ghir_A05G038350 and Ghir_D04G004710) and three (Gbar_A05G037550, Gbar_A09G001130 and Gbar_D04G004660) BCCP2 genes in G. hirsutum and G. barbadense genome, respectively. Unlike all three BCCP2 genes expressed in high SOC of G. barbadense, only one on A05 expressed with a declining trend during seed development in G. hirsutum (Figure 4a). Via alignment of the relevant promoter sequences of the orthologous BCCP2 genes, a 12 bp deletion located ~38‐bp upstream of the TSS was found in the putative promoter in the corresponding regions of Gbar_A05G037550 (Figure 4b). Interestingly, two predicted TF‐motif were located in this missing sequence, including AP2‐motif (TAGAT) and NF‐YB‐motif (ATTGA), which were closely related to the TFs that regulate the expressions of fatty acid genes shown above. Another BCCP2 gene located on D04 was barely expressed in G. hirsutum, while its orthologous gene in G. barbadense had a higher expression level. Moreover, we found that the basis for the differential expression might result from three adjacent deletions in the putative promoter region of the 3‐79 allele of this gene (Figure 4c). Deletion of the motif in the promoter may cause its declining ability to be regulated by TFs and can explain the altered expression of BCCP2 genes between two species (Figure 4d). The high expression of BCCP2 genes provided abundant precursor for FA biosynthesis, which improved the SOC in G. barbadense. However, no significant structure variation was discovered in the core promoter region of the highly expressed oil‐related genes (Table S9), such as FAD and SAD gene. These results imply that there should be more complicated genetic mechanism of regulating SOC in cotton.

Figure 4

Compression of BCCP2 homolog alleles. (a) Expression levels of these two homolog genes in developing seed; (b) the promoter of the Emian22 allele of Ghir_A05G038350 features a 12‐bp insertion compared with the 3‐79 allele Gbar_A05G037550; (c) the promoter of the Emian22 allele of Ghir_D04G004710 features three insertion compared with the 3‐79 allele Gbar_D04G004660; (d) expression levels of these two homolog genes in developing seed.

Identification of candidate genes for SOC

In previous research, we have identified fifteen QTL for SOC in the G. hirsutum cv. Emian22 × G. barbadense acc. 3‐79 interspecific IL population, which provided a good opportunity to deeply analyse the genetic mechanism of SOC. Conducting gene expression profile in embryo allowed us to identify tissue‐specific genes and to reveal different genetic regulation of lipid metabolism between species. To identify causal genes within the QTL intervals, we examined genes located in the confidence intervals by referring to the TM‐1 genome annotation. A total of 1573 genes were identified in these QTL intervals, ranging from 12 genes in q‐SOCD01‐2 to 491 genes in q‐SOCD01‐1. We then compared the expression pattern of these genes and exacted those that were differentially expressed between the two parental lines. Among these genes, eleven candidate genes associated with 8 QTLs can be proposed on the basis of their gene annotation. However, only four of these candidate genes were associated with the introgression chromosome interval carried by the high SOC ILs. To avoid errors caused by QTL detection, we searched all the genes located in the introgression intervals of high SOC ILs, based on sequence homology, as well as differential expression and potentially lipid metabolism‐related. Finally, another ten genes were picked out (Table 1). Most of the candidate genes showed significant differential expression levels during the synthesis and accumulation of cottonseed oil between G. hirsutum and G. barbadense, especially in the late stage of cottonseed oil accumulation (25 DPA to 35 DPA) (Figure 5).

Table 1

Candidate genes proposed for the SOC accumulation in cotton

Gene ID	QTLs	Ils	Annotation	Expression
				Emian22		3–79
				20 DPA	30 DPA	20 DPA	30 DPA
Ghir_A01G002250	q‐SOCA01‐1	N159	Senescence‐associated gene 29	4.52	6.30	26.66	14.11
Ghir_A01G003150	q‐SOCA01‐1	N159	Pyruvate dehydrogenase kinase	1.36	0.91	5.99	2.79
Ghir_A01G008810	q‐SOCA01‐3	–	Phosphoglycerate kinase	10.77	4.14	43.11	26.53
Gbar_A01G009640	q‐SOCA01‐4	–	Catalase 2	–	–	42.07	35.06
Ghir_A01G010240	q‐SOCA01‐4	–	Gamma carbonic anhydrase 1	7.24	10.23	27.04	22.90
Ghir_A05G025470	q‐SOCA05	–	2‐oxoacid dehydrogenase acyltransferase family protein	3.97	0.73	28.22	25.59
Ghir_A05G025600	q‐SOCA05	–	GDSL‐like Lipase/Acylhydrolase superfamily protein	0	0	34.67	26.21
Ghir_A07G019980	q‐SOCA07	N156	Small ubiquitin‐like modifier 1	0.41	0.43	17.76	10.33
Ghir_A11G033230	–	R120	Acyl‐CoA‐binding protein 6	59.08	35.22	57.39	61.48
Ghir_A12G012590	–	N70	Late embryogenesis abundant (LEA) hydroxyproline‐rich glycoprotein family	0.96	5.34	29.63	32.83
Gbar_A12G012650	–	N70	Pyruvate dehydrogenase E1 alpha	–	–	15.05	11.36
Ghir_A12G013320	q‐SOCA12‐1	N70	Genomes uncoupled 1	0	0.50	9.00	9.40
Ghir_D01G013740	q‐SOCD01‐1	–	Fatty acid desaturase 2	5.07	0.32	12.86	20.54
Ghir_D01G015590	q‐SOCD01‐2	N58	Phosphoglucomutase/phosphomannomutase family protein	0	0	0.65	15.78
Ghir_D03G009400	–	R104/R106	Cytochrome p450 78A9	2.46	0	46.90	23.74
Ghir_D03G011050	–	R104/R106	Enolase 1	10.08	6.23	25.54	25.85
Ghir_D06G002090	–	R104	Glycosyl hydrolase family 35 protein	7.00	5.17	23.29	22.11
Ghir_D07G011430	–	N45	Lecithin:cholesterol acyltransferase family protein	14.69	6.23	67.35	63.36
Ghir_D07G012450	–	N45	NAD(P)‐binding Rossmann‐fold superfamily protein	1.53	0.50	24.81	43.00
Ghir_D07G013130	–	N45	Lipid‐transfer protein/seed storage 2S albumin superfamily protein	91.33	11.10	702.92	131.68
Gbar_A05G017010	–	N45	Peroxidase superfamily protein	–	–	10.70	12.52

Figure 5

Relative expression levels of selected candidate genes. Quantitative real‐time PCR data were mean ± standard error of three biological replicates compared with internal control GhUBQ7 (SWEET: senescence‐associated gene 29; PGK: phosphoglycerate kinase; GDSL‐like: GDSL‐like Lipase/Acylhydrolase superfamily protein; ACBP6: acyl‐CoA‐binding protein 6; PGM: phosphoglucomutase/phosphomannomutase family protein; PDH‐E1: pyruvate dehydrogenase E1 alpha; PEP: enolase 1; CYP78A9: cytochrome p450 78A9; EMB3003: 2‐oxoacid dehydrogenase acyltransferase family protein; LEA: late embryogenesis abundant (LEA) hydroxyproline‐rich glycoprotein family; MUM2: glycosyl hydrolase family 35 protein; CAT2: catalase 2; SUMO: small ubiquitin‐like modifier 1; GUN1: genomes uncoupled 1; lecithin:cholesterol acyltransferase family protein; NADP: NAD(P)‐binding Rossmann‐fold superfamily protein; PRXR: peroxidase superfamily protein; and GPI: bifunctional inhibitor/lipid‐transfer protein/seed storage 2S albumin superfamily protein) (one‐way ANOVA for two lines; * P < 0.05, ** P < 0.01, *** P < 0.001).

Among these candidate genes, eleven genes were reported to be involved in FA metabolism or central metabolism. On chromosome A01, for instance, the candidate gene pyruvate dehydrogenase kinase (PDK) in the q‐SOCA01‐1 is only carried by IL N159, which has been described in the previous research (Zhu et al., 2020). Another gene, Ghir_A01G002250, encoding a member of the sweet sucrose efflux transporter family protein, has been identified as key player in intercellular transport of sugars (Eom et al., 2015). The significantly higher expression level of this gene may increased the efficiency of sugar transporting to promote seed oil synthesis in G. barbadense during the 20 DPA to 35 DPA (Figure 5a). On the same chromosome, a gene was identified in the q‐SOCA01‐3 with significant differential expression from 25 DPA to 35 DPA between Emian22 and 3‐79 (Figure 5b). This gene encodes a phosphoglycerate kinase (PGK), which is a highly conserved reversible enzyme that participates in both glycolysis and photosynthesis (Li et al., 2019). The production of dihydroxyacetone‐3‐phosphate (DHAP) from the glycolysis pathway is the intermediate carbon donor for FA synthesis, and previous research has confirmed that high activities of PGK can increase the FA content in sunflower seeds (Troncoso‐Ponce et al., 2010). In additional, the candidate gene, Ghir_A05G025600, classified in a GDSL‐motif acyltransferase, was identified in the region of q‐SOCA05. The expression pattern ‐of this gene showed a high degree of coincidence with the period of seed oil accumulation in G. hirsutum (Figure 5c), while several researches reported that overexpression of the GDSL‐type esterase could reduce the total seed fatty acid content and change the seed fatty acid composition in Arabidopsis (Chen et al., 2012; Huang et al., 2015). Higher expression levels of this gene may negatively regulate the fatty acid accumulation during 20 DPA to 30 DPA in G. hirsutum. Another two candidate genes, Ghir_A11G033230 and Ghir_D01G013740, encode acyl‐CoA‐binding protein 6 (ACBP6) and fatty acid desaturase 2 (FAD2) enzymes known to be directly related to FA metabolism (Tables S6 and S7). The expression abundance of ACBP6 in G. barbadense was significantly higher than that in G. hirsutum during the entire seed development process, except the period of 10 DPA (Figure 5d).

During oilseed development, the intermediate products in the central metabolism, particularly in glycolysis pathway, provide the major precursor and energy source for FA and TAG biosynthesis. Three candidate genes participate in this metabolism pathway (Figure 5e, f and g). For example, in the q‐SOCD01‐2, an ortholog of the Arabidopsis phosphoglucomutase (AtPGM3), Ghir_D01G015590, was found. This gene can catalyse the interconversion of glucose 1‐phosphate (G1P) and glucose 6‐phosphate (G6P), and is a key enzyme of central metabolism (Egli et al., 2010). As we know, G6P is the precursor source of G3P synthesis required for FA and TAG biosynthesis (Baud and Lepiniec, 2010). Another candidate gene, Ghir_D03G011050, associated with ILs R104 and R106, encodes a putative plastid‐localized phosphoenolpyruvate (PEP) enolase. As the main compound in the glycolysis pathway, PEP can convert to pyruvate, which in turn forms acetyl‐CoA to provide a carbon source for FA synthesis (Prabhakar et al., 2009). This gene maintained stable expression abundance from 20 DPA to 35 DPA in G. barbadense (Figure 5g), therefore, and caused a continuous supply of carbon source for FA synthesis.

The gene, Ghir_D03G009400, shared in both ILs R104 and R106, encodes a cytochrome p450 monooxygenase (CYP78A9) and was predicted to involve in the seed development. Although overexpression of this gene induced large and seedless fruits in Arabidopsis (Ito and Meyerowitz, 2000), the analysis of genes coexpressed with CYP78A9 showed that most of them were involving in FA metabolism or FA metabolism process (Sotelo‐Silveira et al., 2013). Moreover, examination of the promoter region of this allele gene between parental lines revealed a 21 bp deletion from the −187 bp of the TAG in the Emian22 (Figure 6a), which may explained the differential expression (Figure 5h). Moreover, we also identified three other genes whose expression patterns were consistent with the accumulation of seed oil (Figure 5i, j and k), such as Ghir_A05G025470, Ghir_A12G012590 and Ghir_D06G002090. No evidence has been reported that these genes are directly involved in the oil biosynthesis pathway, but indispensable genes for seed development (Fiume et al., 2016; Mentzen et al., 2008). Furthermore, genes involved in the redox reactions thought to be related FA metabolism have also been identified in this study. In q‐SOCA01‐4, the candidate gene Gbar_A01G009640, as a gene insertion event in the homology region of G. barbadense, coding a peroxisomal catalase 2, may improve the SOC by regulating the metabolism of H2O2 in cottonseed, based on the report that the formation of high oleic was closely related to the H2O2 content in seeds (Liu et al., 2019). However, significantly different expression level was detected at the maturation stage of the cottonseed (30 DPA to 35 DPA) (Figure 5l), which indicated that the biological regulation of cottonseed oil may be more complicated.

Figure 6

Candidate genes in this study. (a) The promoter variation of candidate gene Ghir_D03G009400; (b) structure variation in the coding sequence of the candidate gene Ghir_D07G012450; (c) function analysis of GbSWEET and GbACBP6 in the fatty acid synthesis in yeast.

No annotated biosynthetic gene likely involved in the production of the FA synthesis was found in the q‐SOCA07. Incidentally, an altered expression gene, Ghir_A07G019980, was found in this region (Figure 5m). This gene was described as a small ubiquitin‐like modifier (SUMO) polypeptide, involved in the post‐translational modification process of protein, which can regulate embryo development(Liu et al., 2017a). Previous studies reported that the retrograde signalling of the candidate gene GhGUN1 (Ghir_A12G013320) may be involved in regulating acetyl‐CoA carboxylase 2 in Arabidopsis (Wang et al., 2018), which is essential for the first step of fatty acid synthesis. The significant difference in the expression pattern of this gene between G. hirsutum and G. barbadense indicated that it may played an important role in the formation of the higher SOC in the G. barbadense (Figure 5n). Four candidate genes with differentially expression were selected from the introgression region of IL line N45 (Figure 5o, p, q and r). One significant frameshift mutation occurs in the coding sequence between Ghir_D07G012450 and its allele gene in G. barbadense (Gbar_D07G012560) (Figure 6b). This allele gene encodes a short‐chain dehydrogenase/reductase in G. barbadense, which plays a role in detoxifying reactive carbonyls in plants (Yamauchi et al., 2011). On the contrary, there is no clear function reported of the homolog allele gene in upland cotton. Moving on to the next gene, Ghir_D07G013130, this gene showed extensive expression abundance, especially after 20 DPA (Figure 5r). Therefore, it is speculated to play dual roles as functional inhibitor for lipid‐transfer and seed storage.

Potential functional roles of candidate genes in fatty acid synthesis

A total of 17 candidate genes were selected to construct the overexpression vector pYES2 in the INVScL1 strain of Saccharomyces cerevisiae. Only two genes (GbSWEET and GbACBP6) have the potential to regulate the FA biosynthesis in the yeast (Figure 6c). Overexpression of the GbSWEET in yeast increased about 17.54%, 22.17% and 40.07% of the OC in the three transgenic strains compared with that in the negative control (CK), which indicated that higher expression level of the Gb SWEET in G. barbadense played a crucial role in the enhancing the SOC. Meanwhile, the oil content in another transgenic yeast strains carried GbACBP6 also showed significantly higher than the empty‐vector.

Discussion

Cotton, beyond its leading nature fibre, seed oil is also a valuable source for food and biodiesel. Although an increasing number of works have reported the biochemical and biological functions of oil‐related genes in upland cotton, less attention was paid to identify corresponding genes in sea‐island cotton. Comparative transcriptome analysis of high‐ and low‐SOC cottonseeds provides an effective approach for studying gene differential expression pattern and dissecting candidate genes in oil synthesis. Here, based on deep RNA sequencing, a high‐resolution gene expression network in developing cottonseed was constructed to analyse the genetic regulation of lipid metabolism in G. barbadense. Different from previous study (Zhao et al., 2018c), the present study found that the number of expressed genes decreased during seed development. Compared with the DEGs in G. hirsutum, significantly reduced number of DEGs appeared in the comparison of ‘30 DPA vs 20 DPA’ in G. barbadense with higher SOC, which meat that most of genes still maintain similar expression pattern in 30 DPA as in 20 DPA. Then, GO and KEGG analysis were performed to identify the functional interpretation and pathway enrichment of those DEGs. The pathway enrichment analysis also confirmed that more up‐regulated and less down‐regulated genes in lipid metabolism pathway were detected during the critical period of oil accumulation in G. barbadense than that in G. hirsutum (Tables S4 and S5). Furthermore, several genes related unsaturated fatty acid synthesis showed significantly down‐regulated in G. hirsutum, such as SAD6 and FATA. In short, the high expression level maintained by a large number of genes at 30 DPA is the direct reason for the higher SOC in G. barbadense compared with G. hirsutum.

During the de novo biosynthesis of FA, the difference was observed in the expressed gene type between the two species. For example, α‐CT showed high expression levels in G. barbadense, while β‐CT is more active in G. hirsutum. Similar result was observed for the two types of BCCP genes. DNA sequencing of the orthologous promoter regions revealed deletions in the GbBCCP2, moving the regulation binding motif away from the TSS. The WRI1 as an AP2 type TF has been shown to positively regulate the expression of the BCCP2 genes (Maeo et al., 2009). As a carrier of FAs, upland cotton relied on more FATB at the later period. Different from the previous studies in upland cotton, the LACS and PDCT seem no value for the TAG formation in sea‐island cotton, because of its barely expressed level. Via overexpression of LACS and PDCT in sea‐island cotton may be a good idea to improve the SOC in G. barbadense cultivars in the future. The LPAAT2 gene was regarded as a rate‐limiting enzyme in the Kennedy pathway in higher plants, and previous study showed that this gene highly expressed in upland cotton (Wang et al., 2017b). On the contrary, decreasing expression trend was observed for this gene in G. barbadense. The different expression pattern of genes related to FA synthesis implies that different properties or mechanisms of oil biosynthesis exist in the two distantly related cotton species.

TFs have been proved to play an important role in controlling SOC in this study. Most of the oil‐related TFs showed abundant expression level from 20 DPA to 30 DPA in G. barbadense and thus may be the cause of significantly higher expression of FA synthesis genes than G. hirsutum. However, no apparently visible sequence mutation was detected in the core promoter region of these TFs. It is possible either that we cannot detect mutations that exit in these genes or that an unknown regulatory mechanism exists. Therefore, a higher quality reference genome or cloning of gene is needed for future analysis. An example, a reported gene, Gbar_A05G017010, encoding a peroxidase superfamily protein (PRXR), was found in the introgression region of IL line N45. Previous study has confirmed that overexpression of this gene could significant improve the oil content in cottonseed (Ma et al., 2019), while no gene allele was annotated in the homology sequence in the G. hirsutum. Interestingly, this allele can be detected in the recently updated G. hirsutum reference genome data (Huang et al., 2020).

Since upland cotton and sea‐island cotton are allopolyploid species, homolog expression bias of genes is a phenomenon worthy of attention. Via comparing the oil‐related gene expression pattern between published data and this study, we found that slightly higher level of expression bias towards the Dt‐subgenome and those genes high expressed in At‐subgenome were more likely to be regulated (Table S8). We highlighted three TFs with expression patterns significantly related to the oil synthesis in this study, which homolog alleles in At‐subgenome showed lower expression levels in G. hirsutum. This result is similar to previous report for cotton fibre evolution, which showed that domestication has dramatically changed the transcriptome during fibre development (Yoo and Wendel, 2014) and more Dt‐bias domestication selection (Wang et al., 2017a).

To investigate the potential genetic basis of SOC variability between G. hirsutum and G. barbadense, we evaluated the genome variation and gene expression patterns in the candidate intervals of previously identified SOC‐related QTL and the introgression chromosome segments carried by the ILs with high SOC. Genes involved in lipid metabolism, glucose metabolism, redox metabolism and embryo development were identified as candidate genes related to the seed oil accumulation. Upon further analysis by quantitative real‐time PCR between G. hirsutum and G. barbadense, most of the candidate genes were highly consistent with the oil accumulation in cotton. In this study, relying on the fatty acid synthesis system, we identified two genes (SWEET and ACBP6) may influence oil accumulation in cottonseed. These results suggested that the method of selecting the candidate genes used in this study was credible. Correspondingly, in many cases, regarding differentially expressed genes as candidate genes can also be supported by fine‐mapping (Che et al., 2015). Nevertheless, future studies are required for validation of candidate genes identified in this study. The use of omics, especially metabolomics combined with IL resources, to analyse the genetic regulation of lipids is also considered to be a feasible solution (Fernandez‐Moreno et al., 2017; Garbowicz et al., 2018).

Conclusion

The current study provides a first comprehensive view of the genetic regulation of FA metabolism in G. barbadense. The interspecific IL population provided an excellent resource to compare analysis of oil‐related genes differential expression and confirmed in this study. Furthermore, some candidate genes have been reported in association with oil accumulation. This work further enriches our understanding of the gene expression network of FA biosynthesis in cotton and guides us future genetic breeding of cotton oil.

Experimental procedures

Plant materials

The interspecific IL population was derived from the cross between G. hirsutum cv. ‘Emian22’ (low SOC) and G. barbadense acc. 3‐79 (high SOC). The SOC data set presented is based on field‐gown over three years in five environments, which has been described in our previous study (Zhu et al., 2020). The genome sequencing data and genotypic data of the ILs also have been described in that study (Zhu et al., 2020). Of this IL population, 8 with high SOC (N45, N58, N70, N156, N159, R104, R106 and R120) and their parents were grown at an experimental field in Huazhong Agricultural University, Wuhan, China (30.4°N,114.2°E), under normal farming conditions. The SOC traits were phenotyped by measuring the total oil percentage using nuclear magnetic resonance (NMR) analyser following the protocol of Zhao et al. (2019).

RNA‐sequencing dada analysis and oil accumulation gene network construction

The embryos were collected at four developmental stages (10, 20, 30 and 40 day postanthesis (DPA)), representing different stages of seed filling and maturation. The samples were quickly frozen in liquid nitrogen for RNA isolation, and three biological replicates were taken from each time point. High‐quality RNA extraction was performed using the RNAprep Pure Plant Kit (Tiangen, Beijing, China). Total RNA at three stages (10, 20 and 30 DPA) from parents was sequenced with the Illumina HiSeq 2000 system (paired‐end 150 bp). The clean RNA‐seq reads of Emian22 and 3‐79 were mapped to the G. hirsutum and G. barbadense reference genome using HISAT 2.0 software (http://ccb.jhu.edu/software.shtml), respectively. Fragments per kilobase of exon per million mapped reads (FPKM) value were calculated by StringTie for gene expression levels (Pertea et al., 2016), and only criterion of FPKM ≥ 1 in each library was considered an expressed gene. The differentially expressed genes (DEGs) were identified using an R package of DESeq (Anders and Huber, 2010). The genes were defined as differentially expressed following these rules: (1) significance P‐value ≤ 0.05 and false discovery rate (FDR) threshold ≤ 0.05; (2) |log2FoldChange| ≥ 2. Gene Ontology (GO) function annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation were carried out using the online data of the new reference genome project (http://www.cottongen.org), which was updated by our laboratory in previously.

Genes involved in FAs biosynthesis were detected and referenced to Arabidopsis acyl‐lipid metabolism (http://aralip.plantbiology.msu.edu/pathways). The gene network of oil accumulation was constructed and modified by referencing previous reported in upland cotton (Hovav et al., 2015; Zhao et al., 2018c). The expression abundance of genes is counted by their putative functions, independent from their relative expression levels.

Candidate gene selection

The candidate genes for SOC were selected from two parts. In the gene network of lipid biosynthesis, the differential expression gene has any function mutation between G. hirsutum and G. barbadense was considered as a candidate gene, including gene doubling and deletion, structural variation in the promoter region and SNP differences in the coding sequences producing non‐conserved protein changes. In the chromosome regions carried by SOC‐related QTLs and ILs with high SOC, candidate genes selected on the basis of: (i) the corresponding gene annotation being related to the lipids biosynthetic process; (ii) the differential gene expression level between G. hirsutum and G. barbadense being consistent with seed developing stages; (iii) gene function mutation requirements being consistent with above. We named candidate genes in the interest interval that satisfied two or three of as described.

Quantitative PCR

Total RNA was isolated using the RNAprep Pure Plant Kit (Tiangen, Beijing, China) from cottonseeds at 10, 15, 20 25, 30, 35 and 40 DPA. High‐quality RNA was reverse‐transcribed using the SuperScript III reverse transcriptase (Invitrogen, Cat. No. 18080‐093, USA). The ABI Prism 7500 system was used to perform qRT‐PCR, with GhUBQ7 (DQ116441) as the internal control. The real‐time PCR primers are listed in Table S10.

Expression of the candidate genes in yeast fatty acid system

The open reading frame (ORF) of these candidate genes was PCR‐amplified from a G. barbadense acc. 3‐79 ovules cDNA library and inserted into pYES2 (Invitrogen) vector to construct overexpression system in the INVScl1 strain of Saccharomyces cerevisiae (Weidi Biotechnology Co., Shanghai, China). The primers used for vector construction are listed in Table S10. The empty‐vector pYES2 was used as the negative control (CK). The transgenic yeast and CK was induced by the method according to the Ma et al. (2019). The oil content of the transgenic yeast and CK was measured by a tissue triglyceride assay kit (Applygen Technologies Inc., Beijing, China).

Competing interests

The authors declare that they have no competing interests.

Ethics approval and consent to participate

Not applicable.

Authors’ contributions

Zhongxu Lin conceived and designed the project. Xiaojing Li performed the field experiments and investigated phenotypic data. De Zhu, Yu Le and Ruiting Zhang conducted the experiments. De Zhu analysed the data and wrote the manuscript draft. Zhongxu Lin revised the manuscript. All authors discussed the results and approved the final manuscript.

Consent for publication

Not applicable.

Figure S1 SOC of ILs and their parents in 2019. (One ANOVA analysis for two lines and Dunnett’s multiple comparison for multiple lines. *** indicated P < 0.001).

Figure S2 Gene expression cluster and summary of expressed genes at different stages. (a) PCA was performed for all library in each species; (b) summary of expressed gene number in each stage; (c) Venn diagram implies that genes are expressed in specific stages.

Figure S3 KEGG enrichment analysis of up‐regulated DEGs in ‘20 DPA vs 10 DPA’ comparison.

Figure S4 KEGG enrichment analysis of down‐regulated DEGs in ‘30 DPA vs 20 DPA’ comparison.

Figure S5 KEGG enrichment analysis of down‐regulated DEGs in the comparison of ‘30 DPA vs 20 DPA’. The red frames mark the pathways involved in lipid metabolism.

Figure S6 Hierarchical clustering of oil‐related DEGs. Enriched gene ontology categories were shown by each cluster (FDR < 0.1). Value of the colour key refers to the log2 of gene expression (fragments per kilobase of exon per million mapped reads, FPKM).

Click here for additional data file.

Table S1 Summary of RNA_seq data quality in different libraries.

Click here for additional data file.

Table S2 GO enrichment analysis of DEGs in cottonseed developmental stages of G. hirsutum.

Click here for additional data file.

Table S3 GO enrichment analysis of DEGs in cottonseed developmental stages of G. barbadense.

Click here for additional data file.

Table S4 KEGG pathways enrichment analysis of different categories of DEGs in G. hirsutum.

Click here for additional data file.

Table S5 KEGG pathways enrichment analysis of different categories of DEGs in G. barbadense.

Table S6 Expression analysis of genes involved in fatty acid metabolism in developing cottonseed in G. hirsutum.

Table S7 Expression analysis of genes involved in fatty acid metabolism in developing cottonseed in G. barbadense.

Table S8 Oil‐related TFs in cottonseed development between G. hirsutum and G. barbadense.

Table S9 Summary of possible genomic variations of the significant differentially expressed oil‐related genes.

Table S10 Primers used in the research.

Click here for additional data file.