Impact of water deficiency on leaf cuticle lipids and gene expression networks in cotton (Gossypium hirsutum L.)
BMC Plant Biology volume 22, Article number: 404 (2022)
Water deficit (WD) has serious effect on the productivity of crops. Formation of cuticular layer with increased content of wax and cutin on leaf surfaces is closely related to drought tolerance. Identification of drought tolerance associated wax components and cutin monomers and the genes responsible for their biosynthesis is essential for understanding the physiological and genetic mechanisms underlying drought tolerance and improving crop drought resistance.
In this study, we conducted comparative phenotypic and transcriptomic analyses of two Gossypium hirsutum varieties that are tolerant (XL22) or sensitive (XL17) to drought stress. XL17 consumed more water than XL22, particularly under the WD conditions. WD significantly induced accumulation of most major wax components (C29 and C31 alkanes) and cutin monomers (palmitic acid and stearic acid) in leaves of both XL22 and XL17, although accumulation of the major cutin monomers, i.e., polyunsaturated linolenic acid (C18:3n-3) and linoleic acid (C18:2n-6), were significantly repressed by WD in both XL22 and XL17. According to the results of transcriptome analysis, although many genes and their related pathways were commonly induced or repressed by WD in both XL22 and XL17, WD-induced differentially expressed genes specific to XL22 or XL17 were also evident. Among the genes that were commonly induced by WD were the GhCER1 genes involved in biosynthesis of alkanes, consistent with the observation of enhanced accumulation of alkanes in cotton leaves under the WD conditions. Interestingly, under the WD conditions, several GhCYP86 genes, which encode enzymes catalyzing the omega-hydroxylation of fatty acids and were identified to be the hub genes of one of the co-expression gene modules, showed a different expression pattern between XL22 and XL17 that was in agreement with the WD-induced changes of the content of hydroxyacids or fatty alcohols in these two varieties.
The results contribute to our comprehending the physiological and genetic mechanisms underlying drought tolerance and provide possible solutions for the difference of drought resistance of different cotton varieties.
Cuticular waxes and cutin are the major structural constituents in the cuticle, which is a hydrophobic lipid covering the epidermal cells layer . As an important strategy of plants to adapt to environmental stresses, cuticular waxes and cutin play a key role in drought tolerance by hindering the appearance of cellular dehydration during drought stress [2, 3]. Cuticle wax is complicated mixtures of long chain fatty acids and their derivatives and its components generally include alkanes, aldehydes, fatty alcohols, ketones, wax esters, etc. [4, 5]. Cutin is embedded in or covered with wax and is comprised of C16 and C18 fatty acids and corresponding oxygenated derivatives, e.g. hydroxy-fatty acids and dicarboxylic acids [6,7,8,9].
Recent studies have established that plants could respond to WD by increasing cuticle wax deposition [1, 10, 11]. It has been reported that tobacco, sesame and soybean gained more wax amount on leaves after short period water shortage [12,13,14]. WD also resulted in significant increase in the amount of cutin monomers, largely manifested as a high increase in C16 and C18 dioic acids . Some genes involved in wax and cutin biosynthesis pathway were significantly induced upon drought stress. ECERIFERUM 1 (CER1) encodes aldehyde decarboxylase, which is in charge of converting aldehydes to alkanes [15, 16]. Under low humidity conditions, Arabidopsis cer1 mutant showed male sterility, while plants with CER1 overexpression exhibited the decreased cuticle permeability and slower response to soil WD . Arabidopsis CER3 (WAX2) is also involved in wax biosynthesis  and cer3 mutant showed alterations in wax synthesis . Similarly, the mutants defective in one or more cutin monomers, such as fdh, lacs2, hth, lcr, and att1, showed significant alteration in cutin deposition and cuticle permeability [8, 9, 19,20,21,22], revealing important roles of cutin in the formation of water permeability barrier of the cuticle membrane . However, so far, the function of definite cutin monomers to permeability of plant cuticle is still unclear and the role of wax components and cutin monomers and drought responsive genes networks still needs to be further studied.
Cotton is not only the most important textile fiber crop, but also an important oil crop, which has important economic value [23, 24]. WD and drought stress seriously affect cotton yield and quality . It has been reported that WD reduced stomatal conductance, photosynthetic rate, and transpiration rate, and could restrain accumulation of dry matter of cotton plants up to 50% . Regarding the molecular mechanism of drought resistance of cotton, the mitogen-activated protein kinase (MAPK) pathway and calcium signaling pathway have been implicated in drought stress response . Two MAPK genes, GhMPK2 and GhMPK16 [27, 28], a bZIP transcription factor GhABF2, and a R2R3-type MYB transcription factor GbMYB5 [29, 30] have been demonstrated to be involved in drought stress response. However, few reports have touched on the relationship between the wax components and cutin monomers of cotton leaves and drought stress responses, and little is known about the regulatory metabolic networks involved in drought stress response in cotton.
In this study, we conducted phenotypic and physiological observations on two cotton varieties with different level of drought tolerance under well-watered (WW) and WD conditions, and performed comparative transcriptome analysis using leaf samples collected from WW or WD cotton plants at two developmental stages. The major aims of the study were to uncover the profiles of wax components and cutin monomers of cotton varieties tolerant or sensitive to drought stress under the WW and WD conditions, to identify the major genes and their associated gene networks responding to drought stress, and to understand the potential molecular mechanisms underlying drought tolerance in cotton.
Materials and methods
Plant materials and water treatments
Two Upland cotton (G. hirsutum L.) varieties (XL17 and XL22, susceptible and tolerant to drought stress, respectively) were legally obtained from Cotton Research Institute, Shihezi University and were used as experimental materials. XL22 (drought tolerant) and XL17 (drought sensitive) were screened through drought tolerance experiment in natural population in our previous study. Compared to that of their respective control plants under the WW conditions, the yield of XL22 decreased by 34.2%, while the yield of XL17 decreased by 41.3% (Table S1). To maintain treatment consistency and comparability, the experiments were carried out in an artificial climate chamber at 28–30 °C, a 16L:8D photoperiod, and 40% relative humidity. The pot experiment was carried out in greenhouse, and each pot (16 cm diameter and 13.7 cm height) was filled with cultivation substrate purchased from the Floragard Vertriebs (Germany).
Two treatments with soil water content of 70% (representing WW) or 10% (representing WD) were used to investigate drought response of cotton plants. Soil moisture was controlled by the weight measurement method according to Hanson . The treatment was adopted after seedling emergence and lasted from seedling stage to bud stage, and each treatment was set up with three replications.
Leaf area, the relative water content, and stomatal density of leaves
Leaf area of cotton plants was measured according to the method reported by Drake et al.  using the inverse fourth leaf of each cotton plant using LI 3100C Area Meter (Li-Cor Inc., USA). The relative water content (RWC) of leaves was determined according to Lu et al. . The stomatal density was measured according to Dunn et al.  using a 1 mm2 leaf disc sampled from the fourth leaf. The stomatal density was expressed as the number of pores per unit area (mm2). The stomatal aperture (µm) was calculated using ImageJ (v1.49) from images collected at × 150 magnification. The average of 12 visual fields of each sample was taken as the result of the sample and every experiment was repeated five times.
Cuticular wax and cutin analysis
The epicuticular wax of leaf was examined using a scanning electron microscope (SEM) as described by Djanaguiraman et al. . The total epicuticular wax content was measured using chloroform extract from cotton leaves according to a previously reported method . The composition and amount of wax were analyzed using a published protocol . The constituent analyses were performed using GCMS-QP2020 with a DB-1 column of 30 m × 0.32 mm and film thickness of 0.1 μm. GC–MS analyses were performed as described by Liu et al. . Each compound was quantified against the internal standard by automatic integration of the peak areas.
The analyses of total cutin content and monomer were performed according to Franke et al.  with some modifications. The extraction was performed at room temperature in glass tubes with teflon-lined screw caps. Soluble lipids were removed by dipping the cotton leaf in 10 mL of a methanol:chloroform (1:1, v/v) mixture with continuous agitation for 7 days and the solvent was changed daily. The samples were dried under a gentle stream of nitrogen gas and then the leaf lipid polyesters were depolymerizated in 6 ml methanol/sulfuric acid/chloroform (10:0.5:1/, v/v/v) mixture at 80 °C with occasional agitation. The resulting cutin monomer fraction was derivatized with BFTSA/pyridine (1:1) for 60 min at 70 °C. Samples were analyzed using GCMS-QP2020 according to Liu et al.  and each compound was quantified on the basis of its total ion current as described by Li-Beisson et al. .
RNA isolation, library preparation and sequencing
The extraction, purification and integrity identification of total RNA were carried out according to our previously published methods [24, 37]. The libraries were sequenced on an Illumina Hi-seq platform. Raw reads were firstly processed through in-house perl scripts . And then, clean reads were obtained by removing low quality reads and reads containing adapter or ploy-N. Using Hisat2 v2.0.5, paired-end clean reads were aligned to the TM-1 reference genome . FeatureCounts v1.5.0-p3 was used to count the read numbers mapped to each gene. Fragments per kilobase of transcript per million mapped reads (FPKM) of each gene was calculated based on the length of the gene and the number of reads mapped to the gene.
Differential expression, gene function annotation, enrichment analysis and co-expression networks of differentially expressed genes (DEGs)
Differential expression analysis was performed using the DESeq R package (1.20.0). The resulting p-value was adjusted for controlling the false discovery rate according to Benjamini and Hochberg . The adjusted p-value < 0.05 and |log2 (fold change) |≥ 2 were used as the criteria for identification of DEGs. Gene Ontology (GO) enrichment analysis of DEGs was implemented by the clusterProfiler R package, in which gene length bias was corrected. GO terms with a corrected p-value less than 0.05 were considered to be significantly enriched. The clusterProfiler R package was used to test the statistical enrichment of DEGs in KEGG pathways. The co-expression networks among DEGs was constructed by weighted correlation network analysis (WGCNA) software package according to Ma et al.  and Cheng et al. . The hub gene was screened and correlation networks were drawn according to Ma et al.  and Cheng et al. .
Quantitative RT-PCR validation of DEGs
Gene specific primers were designed using cDNA sequences of the target genes with Primer Premier program (Table S2). qRT-PCR was carried out as described by Cheng et al. . Three independent biological experiments were performed for each sample of each time point. The relative expression levels were calculated with the cotton ubiquitin gene (GhUBI, XM_012634824) as the reference according to Cheng et al. .
Effects of WD on physiological and morphological characteristics of cotton varieties with different level of drought tolerance
Under the WW conditions, i.e., when the soil moisture was kept at 70% of its water holding capacity, the drought tolerant variety XL22 consumed 190 g water per day from the seedling stage (30 days after emergence) to the bud stage (50 days after emergence), and the drought sensitive variety XL17 consumed 206.7 g water per day during the same time period. Under the WD conditions, i.e., when the soil water content was kept at 10% of its water holding capacity, from the seedling stage to the bud stage, XL22 and XL17 consumed 60 g and 90 g water per day respectively (Fig. S1). It thus seems that the drought tolerant variety (XL22) consumed less water than the drought sensitive variety XL17 under both the WW and WD conditions.
Compared to their respective control plants under the WW conditions, both XL17 and XL22 decreased significantly in plant height under the WD conditions at the seedling and the bud stages, particularly at the bud stage, due to a more significant effect of WD on the growth of the drought sensitive variety (XL17) than that of the drought tolerant variety (XL22) (Fig. 1A). From the seedling stage to the bud stage, the plant height of XL22 increased by 57.7 and 12.3% under the WW and WD conditions, respectively. During the same period, the plant height of XL17 increased by 95.3 and 27.2% under the WW and WD conditions, respectively. The results suggested that the effect of WD on plant height growth of XL22 was 45.4%, and that of XL17 was 68.1%. WD reduced the leaf area of XL22 by 43.6 and 38.2% at the seedling stage and the bud stage, respectively. The leaf area of XL17 seemed to be more inhibited by WD, with a decrease of 53.7% and 51.1% observed at the seedling and bud stages, respectively (Fig. 1B). In both XL22 and XL17, the RWC (relative water content) of leaves was always higher under the WW conditions than under the WD conditions. After the WD treatment, the RWC of XL22 leaves was reduced by 15.51% from the seedling stage to the bud stage, while that of XL17 leaves was decreased by 20.8% (Fig. 1C). These results showed that XL17 was more sensitive to WD than XL22, likely due to its weaker ability in maintaining the RWC of plants.
Leaf stomatal density is a primary determinant of water use efficiency. The adaxial side of leaves had higher stomatal density than that of abaxial side for XL22 and XL17 (Fig. 2A). No significant difference in leaf stomatal density per unit leaf area was observed between XL22 and XL17 under the WW conditions. At the seedling stage, the stomatal density of XL17 and XL22 on the abaxial side increased by 22.55 and 43.18% under the WD conditions, respectively, compared to that of their respective control plants under the WW conditions. At the bud stage, the abaxial stomatal density of XL22 increased by 36.73% under the WD conditions, compared to that of XL22 under the WW conditions. However, WD increased only 2.81% (statistically insignificant) of the abaxial stomatal density in XL17 at the bud stage. Under the WD conditions, the stomatal density on the adaxial side also increased, although the relative proportion of increase was not as much as that of on the abaxial side. The highest stomatal density was observed on adaxial side in the bud stage of XL22 under the WD conditions, with an average density of 184 stomata per mm2 (Fig. 2A). WD-induced increase in stomatal density on both leaf surfaces was accompanied by decrease of stomatal aperture in both XL17 and XL22 of the two developmental stages, particularly at the bud stage. For instance, compared to the WW plants, the WD plants of XL22 and XL17 showed a 45.54 and 57.27% decrease of the stomatal aperture, respectively (Fig. 2B). WD-induced changes in stomatal aperture were even more apparent on the adaxial surface than on the abaxial surface (Fig. 2C).
Effect of WD on cuticle lipids in leaves
Cuticular wax formation on the surface of plant leaves is thought to be associated with drought stress tolerance. We thus compared changes of cuticular wax and cutin in the leaves of XL22 and XL17 under the WW and WD conditions to explore the stoma-independent drought tolerant mechanisms. Under the WW conditions, no difference in leaf wax content was evident between XL22 and XL17. After the WD treatment, the leaf surface of both XL22 and XL17 appeared, based on visual inspection, to be covered by more wax at both the seedling and bud stages, particularly at the bud stage. Further observation by SEM, there was no significant difference in waxy crystals between XL22 and XL17 (Fig. S2). After WD treatment, the surface of cotton leaves was covered with more wax, and waxy crystals on leaves front side was more apparent (Fig. S2). The observation was also supported by the result of significant increase of total wax amount per unit leaf area in both varieties. Under the WW conditions, the total wax content of XL22 leaves was 40.35 µg/cm2 and 40.59 µg/cm2 at the seedling stage and bud stage, respectively. The corresponding values in XL17 were 40.16 µg/cm2 and 40.38 µg/cm2. Under the WD conditions, the total wax content of XL22 leaves was 64.06 µg/cm2 and 64.75 µg/cm2 at the seedling stage and bud stage, respectively. The corresponding values in XL17 were 63.02 µg/cm2 and 63.22 µg/cm2. There was no significant difference in the total wax content between XL22 and XL17 leaves under the WD conditions (Fig S3). Leaf wax constituents were mainly divided into alkanes, fatty acids, and primary alcohols (Fig. 3). Alkanes with a carbon chain length ranging from 27 to 33 were the major components of leaf wax constituents and accounted for over 90% increase of the total wax observed in the WD plants, consistent with the result reported in Arabidopsis . Of the six alkanes, C29 and C31 ones had a higher basal level than other four and were also the most significantly induced by WD. The contents of fatty acids and primary alcohols were relatively low in leaf wax, but significant decrease was observed for most of their components in leaves of the WD plants (Fig. 3). The major primary alcohols are octacosanol (C28-OH) and triacontanol (C30-OH) in the total wax content of cotton leaves. Under the WW conditions, the total content of primary alcohols was respectively 4.18 µg/cm2 and 4.04 µg/cm2 at the seedling stage of XL22 and XL17. Under the WD conditions, its content significantly decreased by 62.1% for XL22 and 66.5% for XL17, respectively. The reduction of primary alcohols at the bud stage seemed to be more significant in XL17 than in XL22 under the WD conditions. Under the WD conditions, the total content of primary alcohols was respectively 1.73 µg/cm2 and 1.02 µg/cm2 at the bud stage of XL22 and XL17 (Fig. 3). Generally, XL22 and XL17 had a very similar profile of wax components under both the WW and WD conditions, but reduction of primary alcohols at the bud stage seemed to be more significant in XL17 than in XL22 under the WD conditions.
According to the GC–MS results (Fig. 4), the main cutin monomers identified in cotton leaves included linolenic acid (C18:3n-3), palmitic acid (C16:0), linoleic acid (C18:2n-6), stearic acid (C18:0), oleic acid (C18:1n-9), hexadecane-1,16-dioic acid (C16:0 DCA), and octadecanedioic acid (C18:0 DCA). Hexacosanol (C26-OH) and octacosanol (C28-OH) were also detected in cotton leaves (Fig. 4). The most abundant cutin monomer in cotton leaves was C18:3n-3, accounting for more than 20% of the total cutin content. Similar to the wax profile, the cutin profile was generally similar between XL22 and XL17 at both developmental stages, for instance, the content of polyunsaturated fatty acids, C18:2n-6 and C18:3n-3, decreased significantly in both varieties in the WD plants (Fig. 4). But the WD-induced increase of C16:0 DCA and C18:0 DCA was significant only in drought sensitive XL17 but not in drought tolerant XL22. In addition, although WD-induced reduction of hexacosanol (C26-OH) was observed in both XL22 and XL17 at both developmental stages, a significant reduction of C26-OH was observed only in the drought sensitive XL17 at the bud stage (Fig. 4).
To explore the molecular mechanism underlying drought stress response and the difference of drought tolerance between XL22 and XL17, we compared leaf transcriptomes of XL22 and XL17 at the seedling stage and the bud stage. In total, 24 libraries (three biological repetitions for the two varieties under each treatment) were sequenced and a total of 1,116,540,102 raw reads were generated. The average Q30 of the reads was 94.48% and the read GC content was 44.20%. Raw reads were filtered to remove low quality ones and a total of 1,092,594,914 clean reads were finally used in alignment. Approximately 97.86% of the clean reads could be aligned to the TM-1 reference genome with ~ 93.08% of them being uniquely aligned (Table S3). These results suggested that the RNA-seq data were high quality and suitable for further analyses.
Using the criteria mentioned in Methods, we identified a total of 22,241 genes to be differentially expressed between the WW and WD plants of XL22 or XL17. Of those differentially expressed genes (DEGs), 11,196 (50.34%) were up-regulated (WD vs WW) and 11,045 (49.66%) were down-regulated (Padj < 0.05; Fig. 5). In both XL22 and XL17, there were more DEGs at the seedling stage than at the bud stage. For XL22, compared with the WW treatment, the WD treatment resulted in a slightly higher number of down-regulated genes than that of up-regulated genes at both developmental stages, and for XL17, this was observed only at the seedling stage but not at the bud stage (Fig. 5A). Of the 22,241 DEGs, 15,303 were non-redundant ones, with 4,046 and 3,310 unique to the seedling and bud stage of XL22, respectively, and 2,076 and 705 unique to the seedling and bud stage of XL17, respectively. A total of 363 DEGs were common between XL22 and XL17 at both developmental stages (Fig. 5B). Of the up-regulated genes in both varieties, 2,437 (21.77%) had a 2–3 folds change in expression while 1,079 (9.64%) had an over 5 folds change in gene expression. Among the down-regulated genes, 2,351 (21.29%) had a 2–3 folds expression difference while 903 (8.18%) had a > 5 folds expression difference (Fig. 5C).
DEGs were confirmed by qRT-PCR analysis
To verify the results of the RNA-seq and identify DEGs, the expression levels and expression trend of 9 genes were analyzed by qRT-PCR, including 5 genes related to wax and cutin biosynthesis (Fig. S4). The results of qRT-PCR of the selected genes were highly identical with those of RNA-seq, showing that the results of RNA-seq were highly reliable (Fig. S4).
Gene ontology analysis of DEGs
Gene Ontology (GO) analysis was performed using the identified DEGs between the WW and WD plants. Using the criterion of corrected PValue ≤ 0.05, we found that the 15,303 non-redundant DEGs were enriched for 205 GO terms (Table S4). Of the 205 GO terms, 177 (including 80 unique ones) and 151 (including 55 unique ones) were enriched at the seedling and the bud stage of XL22, respectively, and 130 (including 37 unique ones) and 148 (including 82 unique ones) were enriched at the seedling and the bud stage of XL17, respectively (Fig. 6A). Twenty one GO terms were enriched at both the seedling and bud stages of XL22, and 9 GO terms were enriched at both the seedling and bud stages of XL17. Forty one GO terms were enriched at the seedling stage of both XL22 and XL17 and 30 GO terms were enriched at the bud stage of both XL22 and XL17 (Fig. 6A). A total of 9 GO terms were commonly enriched at the seedling and bud stages of both XL22 and XL17 (Fig. 6B). The top enriched GO terms were quite different between XL22 and XL17. The top two biological process terms enriched in XL22 were related to multi-organism process and DNA replication (Table S5), whereas those enriched in XL17 were related to fatty acid beta-oxidation and fatty acid catabolic process (Fig. 6C), suggesting that different biological processes could be responsible for the variable drought response of XL22 and XL17, although certain drought induced biological processes are shared by the two varieties.
KEGG pathways of DEGs
We further performed Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis according to Kanehisa et al.  using the 15,303 non-redundant DEGs and found that the DEGs were associated with 119 pathways (Table S6). Many of the significantly changed pathways were related to alpha-linolenic acid metabolism, cutin and suberine and wax biosynthesis, fatty acid metabolism, starch and sucrose metabolism, peroxisome, steroid biosynthesis, and biosynthesis of unsaturated fatty acids (Table S6). The KEGG pathways enriched at the seedling stage of XL22 included those involved in metabolism of alpha-linolenic acid, linoleic acid, galactose, starch and sucrose, biosynthesis of unsaturated fatty acids, and biosynthesis of steroid, cutin, suberine and wax (Fig. S5A). At the seedling stage of XL17, pathways related to metabolism of starch, sucrose, fatty acids, alpha-linolenic acid, porphyrin, and chlorophyll were enriched (Fig. S5B). GH_D06G1740 encoding transcription factor bHLH91, GH_D09G1061 encoding alcohol-forming fatty acyl-CoA reductase, GH_D11G3608 encoding delta 12-fatty-acid desaturase (FAD2), and GH_A11G3588 involved in the biosynthesis of unsaturated fatty acids were among the DEGs involved in the pathways related to biosynthesis of cutin, suberine and wax (Fig. S6A-Fig. S6C, Table S6). In addition to the aforementioned GH_A11G3588 and GH_D11G3608, GH_D01G2159 encoding alcohol dehydrogenase class-P, and two genes, GH_A05G0687 and GH_D05G0687, encoding linoleate 13S-lipoxygenase 2–1 were found in the enriched pathways (Fig. S6D-Fig. S6E, Table S6). The pathways enriched at the bud stage of XL22 were related to metabolism of linoleic acid and biosynthesis of steroid, and the related DEGs included GH_A05G0687 and GH_D05G0687 that encode linoleate 13S-lipoxygenase 2–1. These two genes seemed to be also involved in the pathways enriched at the bud stage of XL17 that were related to jasmonate/oxylipin biosynthesis.
Expression profiles of the DEGs involved in wax and cutin biosynthesis
Wax content was induced by drought stress in both XL22 and XL17, we therefore further investigated the DEGs involved in wax biosynthesis. CER1encodes an aldehyde decarbonylase that catalyzes biosynthesis of alkanes, a major step of wax production . The cotton genome contains multiple CER1 homologs. In the seedling stage, four GhCER1 genes (GH_A05G3496, GH_ D04G0747, GH_ A07G1203, and GH_ D07G1182) were highly expressed in leaves. Of these genes, GH_D07G1182 showed the highest expression level and was significantly induced by WD in both XL22 (increased by 72.5%) and XL17 (increased by 44.5%). At the bud stage, four GhCER1 genes (GH_A05G3496, GH_A06G1829, GH_A07G1203, and GH_D07G1182) were highly expressed in leaves. Their expression levels were slightly decreased by WD in XL22, but the expression levels of GH_A06G1829 and GH_D07G1182 were significantly induced by WD in XL17. From the seedling stage to the bud stage, the GhCER1 expression levels generally showed an increasing trend in both XL22 and XL17 (Fig. 7). Drought induced up-regulation of GhCER1 in both XL22 and XL17 is consistent with the enhanced accumulation of leaf wax content under the WD conditions.
In the wax and cutin biosynthesis pathways, the ω-hydroxylation reaction is typically catalyzed by cytochrome P450 monooxygenases (CYP), and peroxygenase (PXG) catalyzes the hydroperoxide-dependent epoxidation of unsaturated fatty acids to produce epoxy-fatty acids. In the seedling stage, the expression level of GhCYP86 genes (GH_A11G0933, GH_D11G0960, GH_A08G1579, and GH_D08G1598), which encode a very long chain fatty acid hydroxylase specifically involved in cutin biosynthesis, was significantly down-regulated by WD in both XL22 and XL17 (Fig. 7). However, in the bud stage, the expression levels of those GhCYP86 genes were slightly induced by WD in XL22 and significantly repressed by WD in XL17. Of the GhPXG genes, a pair of homoeologs (GH_A12G1707 and GH_D01G1678) was highly expressed in leaves of both XL22 and XL17 at the seedling stage, but displayed an opposite response to WD, with GH_A12G1707 being significantly induced by WD. The transcript abundance of GH_A12G1707 was also significantly induced by WD in both XL22 (> 4-folds) and XL17 (~ 2-folds) at the bud stage. GhALDH encodes an aldehyde dehydrogenase that catalyzes ω-oxo fatty acids to produce α, ω-dicarboxylic fatty acids. Four pairs of homoeologous GhALDH genes (GH_A12G2939 and GH_D12G2962, GH_A11G0436 and GH_D11G0455, GH_A06G1679 and GH_D06G1699, GH_A07G0707 and GH_D07G0693) exhibited high transcript abundance in leaves and their expression levels, especially GH_A12G2939 and GH_D12G2962, were apparently induced by WD (Fig. 7). This was consistent with the significant induction of dicarboxylic acid by WD. Fatty acyl reductase (FAR) is a important enzyme involved in the synthesis of primary alcohols . In Gossypium hirsutum, eight genes encoding FAR have been identified , and three of them, GH_A07G1194 (GhFAR3A.1), GH_D07G1174 (GhFAR3A.2), and GH_D09G1060 (GhFAR3D.2) had high transcript abundance in leaves. Notably, the expression of those GhFARs was significant repressed by WD at both developmental stages, except GH_A07G1194 at the bud stage (Fig. 7).
Gene network analysis with WGCNA
To identify the specific genes that were highly correlated with WD in cotton, the co-expression networks were generated by WGCNA using the all non-redundant DEGs. A total of 13 gene modules associated with the specific expression profiles of different samples were identified (Fig. 8A). Of the 13 modules, 5 (salmon, red, turquoise, black and green) were significantly associated with the two growth stages of XL22 and XL17 under the WW or WD conditions (Fig. 8B). The red module with 627 genes was highly associated with the seedling stage of XL22 under the WD conditions. The salmon with 54 genes was significantly associated with the seedling stage of XL17 under the WW conditions. The black (5440 genes) and green (690 genes) modules were significantly associated with the bud stage of XL22 and XL17, respectively, under the WD conditions (Fig. 8B). The 368 DEGs of the turquoise module were enriched with 23 GO terms, including microtubule binding, tubulin binding, cytoskeletal protein binding, movement of cell or subcellular component, and microtubule-based movement (Fig. S7A). The 5,440 DEGs of the black module were enriched with 43 GO terms, including oxidoreductase activity, dioxygenase activity, methionine adenosyltransferase activity, endopeptidase inhibitor activity, and peptidase inhibitor activity (Fig. S7B). KEGG analysis of the DEGs of the turquoise and black modules found that the turquoise module was enriched with genes involved in DNA replication (Padj = 1.29 × 10–6, 31 genes) for the seedling stage of XL22 under the WW conditions (Fig. S7C), while the black module was enriched with genes involved in the linoleic acid metabolic pathway (Padj = 5.78 × 10–5, 6 genes) for the bud stage of XL22 under the WD conditions (Fig. S7D).
Co-expression gene networks and hub genes
According to the different plant phenotypes and wax compositions and cutin monomers between XL22 and XL17 under WW and WD conditions, we further analyzed the co-expression networks of DEGs of the black and turquoise modules to identify the hub genes. Based on the criteria of eigengene-based connectivity (KME) value ≥ 0.84 and high weight value ≥ 0.19, 18 genes were found to be co-expressed in the black module. One of the GhCER1 genes, GH_D04G0747, was identified to be the hub genes of the black module (Fig. 9A). Based on the criteria of eigengene-based connectivity (KME) value ≥ 0.95 and high weight value ≥ 0.14, 17 genes were found to be co-expressed in the turquoise module. Four GhCYP86 genes (GH_A08G1579, GH_D08G1598, GH_A11G0933, and GH_D11G0960) encoding cytochrome P450 were identified to be the hub genes of the turquoise module (Fig. 9B). We found that these genes are related to the biosynthesis of unsaturated fatty acids, wax and cutin, and their expression might be sensitive to WD.
Expression profiles of the genes involved in biosynthesis of unsaturated fatty acids
Plants could adapt to environmental stresses by changing the proportion of unsaturated fatty acids in their membrane lipids. Therefore, the effect of water deficiency on the expression levels of the key genes of the unsaturated fatty acid synthesis pathway was further analyzed. C18:0-ACP is desaturated by Δ9-stearoyl-ACP desaturase (SAD) to form monounsaturated C18:1n-9-ACP. Six GhSAD genes (GH_A10G1563, GH_D10G1329, GH_A05G3299, GH_D05G3475, GH_A02G0927, and GH_D02G1157) identified in the cotton genome had high expression levels in leaves. Of those GhSAD genes, the expression levels of a pair of homoeologs (GH_A10G1563 and GH_D10G1329) were much higher than that of their homologs, which means that this pair of GhSAD could take a leading role in C18:0-ACP desaturation in leaves. WD treatment reduced the abundance of GH_A10G1563 and GH_D10G1329 transcripts in leaves at the seedling stage of XL22 and XL17, whereas WD had no significant effect on the expression of the two genes in leaves at the bud stage. The major polyunsaturated fatty acid (PUFA) in cottonseed is synthesized by fatty acid desaturase 2 (FAD2, C18:1n-9 to C18:2n-6 desaturation) and polyunsaturated linolenic acid (C18:3n-3) is desaturated from C18:2n-6 by fatty acid desaturase 3 (FAD3). Four pairs of homoeologous GhFAD2 with different expression patterns have been identified in the cotton genome . The homoeologous pair of GH_A11G3588 and GH_D11G3608 were highly expressed in cotton leaves, indicating that they might be the major contributors in catalyzing C18:2n-6 biosynthesis in cotton leaves. Under the WD conditions, the expression levels of GH_ A11G3588 and GH_D11G3608 decreased significantly in leaves of XL17 and XL22 at the seedling stage, resulting in a lower accumulation of C18:2n-6 in leaves. In the WW plants, the expression levels of GH_ A11G3588 and GH_D11G3608 in both XL17 and XL22 increased slightly from the seedling stage to the bud stage, consistent with a higher accumulation of C18:2n-6 in leaves of both varieties at the bud stage. Compared with the WW plants, the WD-treated plants had a significantly decreased expression level of GH_A11G3588 and a similar level of GH_D11G3608 at the bud stage. Consequently, the accumulation of C18:2n-6 in leaves under the WD conditions was lower than that under the WW conditions at the bud stage. In cotton, FAD3 is mainly encoded by a pair of homoeologous GhFAD3 genes (GH_A10G2629 and GH_D10G2734). Their expression level in leaves was significantly repressed by WD at the bud stage, consistent with decreased accumulation of C18:3n-6 in leaves of XL22 and XL17 under the WD conditions. In Gossypium, the low temperature could induce an increase in the mRNA levels of GhFAD2 to produce more C18:2n-6 . Different from low temperature stress, our results showed that under WD condition, the significant decrease of PUFA was consistent with those reported under salt stress [46,47,48].
Drought is increasingly becoming one of the serious abiotic stresses inhibiting crop productivity. Drought stress significantly affects plant growth and development by reducing plant height and leaf area . Understanding the physiological and molecular responses of plants induced by WD is the key to find solutions for mitigating the effect of drought stress on crop productivity. In this study, using two Upland cotton varieties tolerant (XL22) or sensitive (XL17) to drought stress, we compared their stoma-dependent (stomal density and aperture) and stoma-independent (leaf wax and cutin) responses to WD, and did comparative transcriptomic analysis to identify drought responsive genes and networks.
Under the WW conditions, there was no significant difference in stomatal density per unit of leaf area, leaf cuticular wax content, and RWC between the drought tolerant XL22 and the drought sensitive XL17, although XL17 seemed to consume slightly more water daily than XL22 during the period from 30 to 50 days after seedling emergence. WD stress increased stomata density per unit of leaf area and reduced stomatal aperture in both XL22 and XL17. The drought-induced increase in stomatal density per unit of leaf area had also been observed in Triticum aestivum . In contrast, reduced stomatal density was observed in Populus balsamifera under WD condition . Xu and Zhou reported an increase in the stomatal density on leaves of Leymus chinensis under moderate drought stress, while stomatal density tended to decrease under more severe drought stress conditions . Under drought conditions, there was no significant change in stomatal density of Arachis hypogaea . Taking into account the decrease of leaf area under WD conditions, compared with that of the controls, the total stomatal density of XL22 and XL17 leaves decreased respectively by 32% and 40% in our study, although the stomatal density per unit leaf area increased. The negative relationship between stomatal density and aperture has been previously documented in other plants under WD conditions [50, 54, 55]. The balance between stomatal density and size is coordinated that is related to the limitation of the leaf area allocated to stomata . Bi et al.  found that the response of stomatal density to drought varies with different varieties and the difference in stomatal density under drought could not explain the difference of water loss rate among varieties, and thought that the cuticle composition is important for water loss, rather than a simple correlation with stomatal density under WD conditions. In our study, WD stress increased stomata density per unit leaf area and decreased stomatal aperture in both XL22 and XL17, also suggesting the importance of stomata in regulation of drought stress.
Previous studies have shown that WD-induced increase of cuticle wax produces a cuticle layer with poor water permeability to limit transpiration and delay the occurrence of cell dehydration stress [10, 57, 58]. It has also been reported that WD increases the amounts of nearly all leaf cutin monomers . We also found that WD stress had a significant effect on the accumulation of cuticular wax and cutin in cotton leaves. Accumulation of two of the major leaf wax components, the ultra-long chain C29 and C31 alkanes, and one of the major cutin monomers, palmitic acid (C16:0), were significantly induced by WD stress in both XL22 and XL17, whereas accumulation of one of the major cutin monomers, linolenic acid (C18:3n-3) was significantly repressed by WD stress in both XL22 and XL17. Two cutin monomers, hexadecane-1,16-dioic acid (C16:0 DCA) and octadecanedioic acid (C18:0 DCA), were induced by WD stress in both XL22 and XL17, with a more significant induction observed in the drought sensitive XL17. These results suggest that drought stress induced both stoma-dependent and stoma-independent responses in drought tolerant or sensitive cotton varieties and the role of C16:0 DCA and C18:0 DCA in stoma-independent drought tolerance is yet to be further investigated.
Leaf epidermal wax is the protective layer covering the outermost layer of plants. It can effectively prevent the non-stomatal loss of water in plants and reduce transpiration, and plays an important role in maintaining leaf water content [1, 33]. After water deficiency, the increases of leaf wax amount in XL22 and XL17 indicate that their wax biosynthesis pathways were enhanced by drought stress. Similar results have been previously reported in other plants [10, 11, 13, 58,59,60,61]. Our results also showed that alkanes are the main waxy components in cotton leaves and the most notable change in the wax constituent profile in WD cotton leaves was the increase of alkanes, consistent with the results reported in Arabidopsis . CER1 plays an important role in alkane synthesis [15, 62]. Overexpression of CER1 in plants showed decreased cuticle permeability and reduced sensitivity to soil WD . In this study, GhCER1 was identified to be the hub genes of the black module. Transcriptome analysis showed that water deficiency induced the expression of GhCER1, which might contribute to the increase of the total content of alkanes, the major components of leaf wax, in leaves of the WD cotton plants.
Fatty alcohols are important components of aliphatics and have a crucial role in cuticle permeability and biosynthesis of wax, suberin and cutin [37, 63]. Studies presented here demonstrated that the relative content of fatty alcohols decreased significantly under the WD conditions. RNA-seq analysis showed that genes encoding the long chain fatty acid hydroxylase, such as GhCYP86 and GhFAR, were less abundant in cotton leaves experienced WD. AtFAR3/CER4 encodes FAR involved in the production of long chain primary fatty alcohols of Arabidopsis leaf cuticular waxes . The expression of several genes including AtFAR3 could be activated by MYB94 transcription factor and the amounts of hydroxy fatty acids increased by 39% in MYB94-overexpressing line relative to the WT . Under WD conditions, compared to that in WT leaves, cuticular transpiration occurred more slowly in MYB94-overexpressing line leaves due to the increase in cuticular wax amount . Silencing GhFAR in cotton leads to increased susceptibility to drought stress  and complete elimination of fatty alcohols is lethal to plants . During biosynthesis of the cutin monomers, the ω-hydroxylation reaction can be catalyzed by cytochrome P450 monooxygenases, e.g. CYP86 [22, 37]. It had also been reported that the mutants in P450-genes encoding fatty acid ω-hydroxylases presumably involved in the pathway to α, ω-diacids, show increased permeability to water vapor  or become permeable for the lipid boundary of the outer wall of epidermal cells . Consistent with the gene expression data, the total contents of hydroxyacids notably decreased in response to WW conditions. Interestingly, under the WD conditions, a significant expression difference of the GhCYP86 genes (GH_A11G0933, GH_D11G0960, GH_A08G1579, and GH_D08G1598) was observed between XL22 and XL17 at both the seedling and the bud stages. Decrease of the GhCYP86 (GH_A11G0933, GH_D11G0960, GH_A08G1579, and GH_D08G1598) transcription levels caused by water shortage was much greater in XL17 than in XL22. GH_A11G0933, GH_D11G0960, GH_A08G1579, and GH_D08G1598 were also identified to be the hub genes of the turquoise module. Consistent with the transcriptional abundance of GhCYP86, a more significant decrease of the fatty alcohol content was observed in drought sensitive XL17 than in drought tolerant XL22. Whether the different expression changes of these genes in XL22 and XL17 might contribute to the drought sensitivity of the two varieties is an open question.
In conclusion, WD stress changes stomatal density and aperture as well as accumulation of wax and cutin in leaves in both drought tolerant and sensitive cotton varieties. WD-induced changes of certain wax components and cutin monomers specific to XL22 or XL17 were observed, implying their potential role in drought tolerance or sensitivity, which is worth of further investigation. Transcriptomic analysis identified DEGs between WW and WD as well as between XL22 and XL17, and hub genes and their associated gene networks related to response to drought stress. GhCYP86 (GH_A11G0933, GH_D11G0960, GH_A08G1579, and GH_D08G1598) encoding long chain fatty acid hydroxylase are particularly of interest, as they showed consistent differential expression between XL22 and XL17 under the WW and WD conditions, implying their potential contribution to the difference of drought tolerance of the two cotton varieties.
Availability of data and materials
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/bioproject/;PRJNA776142.
Kosma DK, Bourdenx B, Bernard A, Parsons EP, Lü S, Joubès J, Jenks MA. The impact of water deficiency on leaf cuticle lipids of Arabidopsis. Plant Physiol. 2009;151(4):1918–29. https://doi.org/10.1104/pp.109.141911.
Zhou A, Liu E, Liu J, Feng S, Gong S, Wang J. Characterization of increased cuticular wax mutant and analysis of genes involved in wax biosynthesis in Dianthus spiculifolius. Hortic Res. 2018;15:40. https://doi.org/10.1038/s41438-018-0044-z.
Wu H, Shi S, Lu X, Li T, Wang J, Liu T, Zhang Q, Sun W, Li C, Wang Z, Chen Y, Quan L. Expression analysis and functional characterization of CER1 family genes involved in very-long-chain alkanes biosynthesis in Brachypodium distachyon. Front Plant Sci. 2019;10:1389. https://doi.org/10.3389/fpls.2019.01389.
Samuels L, Kunst L, Jetter R. Sealing plant surfaces: cuticular wax formation by epidermal cells. Annu Rev Plant Biol. 2008;59:683–707. https://doi.org/10.1146/annurev.arplant.59.103006.093219.
Lee SB, Suh MC. Recent advances in cuticular wax biosynthesis and its regulation in Arabidopsis. Mol Plant. 2013;6(2):246–9. https://doi.org/10.1093/mp/sss159.
Pollard M, Beisson F, Li Y, Ohlrogge JB. Building lipid barriers: biosynthesis of cutin and suberin. Trends Plant Sci. 2008;13(5):236–46. https://doi.org/10.1016/j.tplants.2008.03.003.
Li Y, Beisson F. The biosynthesis of cutin and suberin as an alternative source of enzymes for the production of bio-based chemicals and materials. Biochimie. 2009;91(6):685–91. https://doi.org/10.1016/j.biochi.2009.03.016.
Kurdyukov S, Faust A, Trenkamp S, Bär S, Franke R, Efremova N, Tietjen K, Schreiber L, Saedler H, Yephremov A. Genetic and biochemical evidence for involvement of HOTHEAD in the biosynthesis of long-chain alpha-, omega-dicarboxylic fatty acids and formation of extracellular matrix. Planta. 2006;224(2):315–29. https://doi.org/10.1007/s00425-005-0215-7.
Bessire M, Chassot C, Jacquat AC, Humphry M, Borel S, Petétot JM, Métraux JP, Nawrath C. A permeable cuticle in Arabidopsis leads to a strong resistance to Botrytis cinerea. EMBO J. 2007;26(8):2158–68. https://doi.org/10.1038/sj.emboj.7601658.
Kosma DK, Jenks MA. Eco-physiological and molecular-genetic determinants of plant cuticle function in drought and salt stress tolerance. Advances in Molecular Breeding toward Drought and Salt Tolerant Crops. 2007;91–120. https://doi.org/10.1007/978-1-4020-5578-2-5.
Shepherd T, Wynne GD. The effects of stress on plant cuticular waxes. New Phytol. 2006;171(3):469–99. https://doi.org/10.1111/j.1469-8137.2006.01826.x.
Cameron KD, Teece MA, Smart LB. Increased accumulation of cuticular wax and expression of lipid transfer protein in response to periodic drying events in leaves of tree tobacco. Plant Physiol. 2006;140(1):176–83. https://doi.org/10.1104/pp.105.069724.
Kim KS, Park SH, Kim DK, Jenks MA. Influence of water deficit on leaf cuticular waxes of soybean (Glycine max [L.] Merr.). Int J Plant Sci. 2007;168(3):307–16. https://doi.org/10.1086/510496.
Kim KS, Park SH, Jenks MA. Changes in leaf cuticular waxes of sesame (Sesamum indicum L.) plants exposed to water deficit. J Plant Physiol. 2007;164(9):1134–43. https://doi.org/10.1016/j.jplph.2006.07.004.
Aarts MG, Keijzer CJ, Stiekema WJ, Pereira A. Molecular characterization of the CER1 gene of Arabidopsis involved in epicuticular wax biosynthesis and pollen fertility. Plant Cell. 1995;7(12):2115–27. https://doi.org/10.1105/tpc.7.12.2115.
Bourdenx B, Bernard A, Domergue F, Pascal S, Léger A, Roby D, Pervent M, Vile D, Haslam RP, Napier JA, Lessire R, Joubès J. Overexpression of Arabidopsis ECERIFERUM1 promotes wax very-long-chain alkane biosynthesis and influences plant response to biotic and abiotic stresses. Plant Physiol. 2011;156(1):29–45. https://doi.org/10.1104/pp.111.172320.
Chen X, Goodwin SM, Boroff VL, Liu X, Jenks MA. Cloning and characterization of the WAX2 gene of Arabidopsis involved in cuticle membrane and wax production. Plant Cell. 2003;15(5):1170–85. https://doi.org/10.1105/tpc.010926.
Rowland O, Lee R, Franke R, Schreiber L, Kunst L. The CER3 wax biosynthetic gene from Arabidopsis thaliana is allelic to WAX2/YRE/FLP1. FEBS Lett. 2007;581(18):3538–44. https://doi.org/10.1016/j.febslet.2007.06.065.
Pruitt RE, Vielle-Calzada JP, Ploense SE, Grossniklaus U, Lolle SJ. FIDDLEHEAD, a gene required to suppress epidermal cell interactions in Arabidopsis, encodes a putative lipid biosynthetic enzyme. Proc Natl Acad Sci USA. 2000;97(3):1311–26. https://doi.org/10.1073/pnas.97.3.1311.
Schnurr J, Shockey J, Browse J. The acyl-CoA synthetase encoded by LACS2 is essential for normal cuticle development in Arabidopsis. Plant Cell. 2004;16(3):629–42. https://doi.org/10.1105/tpc.017608.
Wellesen K, Durst F, Pinot F, Benveniste I, Nettesheim K, Wisman E, Steiner-Lange S, Saedler H, Yephremov A. Functional analysis of the LACERATA gene of Arabidopsis provides evidence for different roles of fatty acid omega-hydroxylation in development. Proc Natl Acad Sci U S A. 2001;98(17):9694–9. https://doi.org/10.1073/pnas.171285998.
Xiao F, Goodwin SM, Xiao Y, Sun Z, Baker D, Tang X, Jenks MA, Zhou JM. Arabidopsis CYP86A2 represses Pseudomonas syringae type III genes and is required for cuticle development. EMBO J. 2004;23(14):2903–13. https://doi.org/10.1038/sj.emboj.7600290.
Liu F, Zhao YP, Zhu HG, Zhu QH, Sun J. Simultaneous silencing of GhFAD2-1 and GhFATB enhances the quality of cottonseed oil with high oleic acid. J Plant Physiol. 2017;21(5):132–9. https://doi.org/10.1016/j.jplph.2017.06.001.
Ma L, Cheng X, Wang C, Zhang X, Xue F, Li Y, Zhu Q, Sun J, Liu F. Explore the gene network regulating the composition of fatty acids in cottonseed. BMC Plant Biol. 2021;13, 21(1):177. https://doi.org/10.1186/s12870-021-02952-4.
Ullah A, Sun H, Yang X, Zhang X. Drought coping strategies in cotton: increased crop per drop. Plant Biotechnol J. 2017;15(3):271–84. https://doi.org/10.1111/pbi.12688.
Hejnák V, Tatar Ö, Atasoy GD, Martinková J, Çelen AE, Hnilička F, Skalický M. Growth and photosynthesis of Upland and Pima cotton: response to drought and heat stress. Plant Soil Environ. 2015;61:507–14. https://doi.org/10.17221/512/2015-PSE.
Zhang L, Xi D, Li S, Gao Z, Zhao S, et al. A cotton group C MAP kinase gene, GhMPK2, positively regulates salt and drought tolerance in tobacco. Plant Mol Biol. 2011;77(1–2):17–31. https://doi.org/10.1007/s11103-011-9788-7.
Shi J, Zhang L, An H, Wu C, Guo X. GhMPK16, a novel stress-responsive group D MAPK gene from cotton, is involved in disease resistance and drought sensitivity. BMC Mol Biol. 2011;16:12–22. https://doi.org/10.1186/1471-2199-12-22.
Liang C, Meng Z, Meng Z, Malik W, Yan R, Lwin KM, Lin F. GhABF2, a bZIP transcription factor, confers drought and salinity tolerance in cotton (Gossypium hirsutum L.). Sci Rep. 2016;7(6):35040. https://doi.org/10.1038/srep35040.
Chen T, Li W, Hu X, Guo J, Liu A, Zhang B. A cotton MYB transcription factor, GbMYB5, is positively involved in plant adaptive response to drought stress. Plant Cell Physiol. 2015;56(5):917–29. https://doi.org/10.1093/pcp/pcv019.
Hanson B. Field Estimation of Soil Water Content: A Practical Guide to Methods, Instrumentation and Sensor Technology: International Atomic Energy Agency. Vienna, Austria. Vadose Zone J. 2009;8(3):628. https://doi.org/10.2136/vzj2008.0171.
Drake JE, Aspinwall MJ, Pfautsch S, Rymer PD, Reich PB, Smith RA, Crous KY, Tissue DT, Ghannoum O, Tjoelker MG. The capacity to cope with climate warming declines from temperate to tropical latitudes in two widely distributed Eucalyptus species. Glob Chang Biol. 2015;21(1):459–72. https://doi.org/10.1111/gcb.12729.
Lu Y, Cheng X, Jia M, Zhang X, Xue F, Li Y, Sun J, Liu F. Silencing GhFAR3.1 reduces wax accumulation in cotton leaves and leads to increased susceptibility to drought stress. Plant Direct. 2021;5(4):e00313. https://doi.org/10.1002/pld3.313.
Dunn J, Hunt L, Afsharinafar M, Meselmani MA, Mitchell A, Howells R, Wallington E, Fleming AJ, Gray JE. Reduced stomatal density in bread wheat leads to increased water-use efficiency. J Exp Bot. 2019;70(18):4737–48. https://doi.org/10.1093/jxb/erz248.
Djanaguiraman M, Perumal R, Jagadish SVK, Ciampitti IA, Welti R, Prasad PVV. Sensitivity of sorghum pollen and pistil to high-temperature stress. Plant Cell Environ. 2018;41:1065–82. https://doi.org/10.1111/pce.13089.
Li-Beisson Y, Shorrosh B, Beisson F, Andersson MX, Arondel V, Bates PD, Baud S, Bird D, Debono A, Durrett TP, Franke RB, Graham IA, Katayama K, Kelly AA, Larson T, Markham JE, Miquel M, Molina I, Nishida I, Rowland O, Samuels L, Schmid KM, Wada H, Welti R, Xu C, Zallot R, Ohlrogge J. Acyl-lipid metabolism. Arabidopsis Book. 2013;11:e0161. https://doi.org/10.1199/tab.0161.
Liu F, Ma L, Wang Y, Li Y, Zhang X, Xue F, Nie X, Zhu Q, Sun J. GhFAD2-3 is required for anther development in Gossypium hirsutum. BMC Plant Biol. 2019;19(1):393. https://doi.org/10.1186/s12870-019-2010-9.
Franke R, Briesen I, Wojciechowski T, Faust A, Yephremov A, Nawrath C, Schreiber L. Apoplastic polyesters in Arabidopsis surface tissues - a typical suberin and a particular cutin. Phytochemistry. 2005;6(22):2643–58. https://doi.org/10.1016/j.phytochem.2005.09.027.
Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009;10(1):57–63. https://doi.org/10.1038/nrg2484.
Hu Y, Chen J, Fang L, Zhang Z, Ma W, Niu Y, Ju L, Deng J, Zhao T, Lian J, Baruch K, Fang D, Liu X, Ruan YL, Rahman MU, Han J, Wang K, Wang Q, Wu H, Mei G, Zang Y, Han Z, Xu C, Shen W, Yang D, Si Z, Dai F, Zou L, Huang F, Bai Y, Zhang Y, Brodt A, Ben-Hamo H, Zhu X, Zhou B, Guan X, Zhu S, Chen X, Zhang T. Gossypium barbadense and Gossypium hirsutum genomes provide insights into the origin and evolution of allotetraploid cotton. Nat Genet. 2019;51(4):739–48. https://doi.org/10.1038/s41588-019-0371-5.
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R STAT SOC B. 1995;57:289–300. https://doi.org/10.1111/j.2517-6161.1995.tb02031.x.
Cheng XQ, Zhang XY, Xue F, Zhu SH, Li YJ, Zhu QH, Liu F, Sun J. Characterization and transcriptome analysis of a dominant genic male sterile cotton mutant. BMC Plant Biol. 2020;20(1):312. https://doi.org/10.1186/s12870-020-02522-0.
Kanehisa M, Goto S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28:27–30. https://doi.org/10.1093/nar/28.1.27.
Wang M, Wu H, Xu J, Li C, Wang Y, Wang Z. Five fatty acyl-Coenzyme A reductases are involved in the biosynthesis of primary alcohols in Aegilops tauschii leaves. Front Plant Sci. 2017;8:1012. https://doi.org/10.3389/fpls.2017.01012.
Kargiotidou A, Deli D, Galanopoulou D, Tsaftaris A, Farmaki T. Low temperature and light regulate delta 12 fatty acid desaturases (FAD2) at a transcriptional level in cotton (Gossypium hirsutum). J Exp Bot. 2008;59(8):2043–56. https://doi.org/10.1093/jxb/ern065.
Kerkeb L, Donaire JP, Venema K, Rodríguez-Rosales MP. Tolerance to NaCl induces changes in plasma membrane lipid composition, fluidity and H+-ATPase activity of tomato calli. Physiol Plant. 2001;113:217–24. https://doi.org/10.1034/j.1399-3054.2001.1130209x.
Mansour MMF, van Hasselt PR, Kuiper PJC. Plasma membrane lipid alterations induced by NaCl in winter wheat roots. Physiol Plant. 1994;92:473–8. https://doi.org/10.1111/j.1399-3054.1994.tb08838.x.
Wu J, Seliskar DM, Gallagher JL. Stress tolerance in the marsh plant Spartina patens: Impact of NaCl on growth and root plasma membrane lipid composition. Physiol Plantarum. 1998;102(2):307–17. https://doi.org/10.1034/j.1399-3054.1998.1020219.x.
Zahid Z, Khan MKR, Hameed A, Akhtar M, Ditta A, Hassan HM, Farid G. Dissection of drought tolerance in upland cotton through morpho-physiological and biochemical traits at seedling stage. Front Plant Sci. 2021;12: 627107. https://doi.org/10.3389/fpls.2021.627107.
Bi H, Kovalchuk N, Langridge P, Tricker PJ, Lopato S, Borisjuk N. The impact of drought on wheat leaf cuticle properties. BMC Plant Biol. 2017;17(1):85. https://doi.org/10.1186/s12870-017-1033-3.
Hamanishi ET, Thomas BR, Campbell MM. Drought induces alterations in the stomatal development program in Populus. J Exp Bot. 2012;63(13):4959–71. https://doi.org/10.1093/jxb/ers177.
Xu Z, Zhou G. Responses of leaf stomatal density to water status and its relationship with photosynthesis in a grass. J Exp Bot. 2008;59:3317–25. https://doi.org/10.1093/jxb/ern185.
Clifford SC, Black CR, Roberts JA, Stronach IM, Singleton-Jones PR, Mohamed AD, Azam-Ali SN. The effect of elevated atmospheric CO2 and drought on stomatal frequency in groundnut (Arachis hypogaea (L.)). J Exp Bot. 1995;46(7):847–52. https://doi.org/10.1093/jxb/46.7.847.
Franks PJ, Drake PL, Beerling DJ. Plasticity in maximum stomatal conductance constrained by negative correlation between stomatal size and density: an analysis using Eucalyptus globules. Plant Cell Environ. 2009;32:1737–48. https://doi.org/10.1111/j.1365-3040.2009.02031.x.
Lawson SS, Pijut PM, Charles H, Michler CH. Comparison of Arabidopsis stomatal density mutants indicates variation in water stress responses and potential Epistatic Effects. J Plant Biol. 2014;57:162–73. https://doi.org/10.1007/s12374-014-0017-1.
Mansouri S, Radhouane L. Is Stomatal density an ideal tool to explain drought effect on photosynthesis in tunisian barley varieties? J Chem Bio Phy Sci Sec B. 2015;5(3):2853–64.
Williams MH, Rosenqvist E, Buchhave M. Response of potted miniature roses (Rosa × hybrida) to reduced water availability during production. J Hortic Sci Biotechnol. 1999;74:301–8.
Jenks MA, Andersen L, Teusink RS, Williams MH. Leaf cuticular waxes of potted rose cultivars as affected by plant development, drought and paclobutrazol treatments. Physiol Plant. 2001;112(1):62–70. https://doi.org/10.1034/j.1399-3054.2001.1120109.x.
Jefferson P, Johnson D, Rumbaugh M, Asay K. Water stress and genotypic effects on epicuticular wax production of alfalfa and crested wheatgrass in relation to yield and excised leaf water loss rate. Can J Plant Sci. 1989;69:481–90.
Bondada BR, Oosterhuis DM, Murphy JB, Kim KS. Effect of water stress on the epicuticular wax composition and ultrastructure of cotton (Gossypium hirsutum L.) leaf, bract, and boll. Environ Exp Bot. 1996;36:61–7. https://doi.org/10.1016/0098-8472(96)00128-1.
Samdur MY, Manivel P, Jain VB, Chikani BM, Gor HK, Desai S, Misra JB. Genotypic differences and water-deficit induced enhancement in epicuticular wax load in peanut. Crop Sci. 2003;43:1294–9.
Jenks MA, Tuttle HA, Eigenbrode SD, Feldmann KA. Leaf epicuticular waxes of the eceriferum mutants in Arabidopsis. Plant Physiol. 1995;108:369–77.
Vishwanath SJ, Kosma DK, Pulsifer IP, Scandola S, Pascal S, Joubès J, Dittrich-Domergue F, Wang R, Gao M, Ji S, Wang S, Meng Y, Zhou Z. Carbon allocation, osmotic adjustment, antioxidant capacity and growth in cotton under long-term soil drought during flowering and boll-forming period. Plant Physiol Biochem. 2016;107:137–46. https://doi.org/10.1016/j.plaphy.2016.05.035.
Rowland O, Zheng HQ, Hepworth SR, Lam P, Jetter R, Kunst L. CER4 encodes an alcohol-forming fatty acyl-Coenzyme A reductase involved in cuticular wax production in Arabidopsis. Plant Physiol. 2006;142(3):866–77. https://doi.org/10.1104/pp.106.086785.
Lee SB, Suh MC. Cuticular Wax biosynthesis is up-regulated by the MYB94 transcription factor in Arabidopsis. Plant Cell Physiol. 2015;56(1):48. https://doi.org/10.1093/pcp/pcu142.
We thank all the members of The Key Lab of Oasis Eco-agriculture, College of Agriculture, Shihezi University, Xinjiang, China, for their support throughout the study.
This work was financially supported by the State Key Program of National Natural Science Foundation of China (Grant No. U1903204) and the National Natural Science Foundation of China .
Ethics approval and consent to participate
All methods were carried out in accordance with relevant institutional, national, and international guidelines and legislation. G. hirsutum L. varieties XL17 and XL22 was legally obtained from Cotton Research Institute, Shihezi University and licensed for scientific research.
Consent for publication
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Yield in field trials of XL22 and XL17 under the WW and WD conditions.
Primers used in qRT-PCR analysis.
Transcriptome sequencing data quality and genome mapping.
Analysis of GO enrichment for 15,303 common DEGs.
Analysis of GO enrichment from the seedling stage and bud stage of XL22.
Analysis of KEGG enrichment of the DEGs.
Fig. S1. Consumed water analysis of XL22 and XL17 under the WW and WD conditions. A. Plant phenotyping of XL22 and XL17 at bud stage under the WW and WD conditions. B. Consumed water content of XL22 and XL17 under the WW and WD conditions. WW, well-watered; WD, water deficit. The scales bars are indicated with white lines. Error bars are standard errors. Values represent the means ± SE, n = 3. Different letters above the bars indicate statistically different from each other as determined by the Student’s t test: p < 0.05. Fig. S2. SEM observation of the leaves surface wax of XL22 and XL17 under the WW and WD conditions. a and b. XL22 leaf at the seedling stage under WW and WD condition, respectively. c and d. XL22 leaf at the bud stage under WW and WD condition, respectively. e and f. XL17 leaf at the seedling stage under WW and WD condition, respectively. g and h. XL17 leaf at the bud stage under WW and WD condition, respectively. i, j, k and l, the adaxial surface of XL22 leaf corresponds to a, b, c and d. m, n, o and p, the adaxial surface of XL17 leaf corresponds to e, f, g and h. The waxy crystals as indicated by the arrows. WW, well-watered; WD, water deficit. SS-22, XL22 leaves at seedling stage; SS-17, XL17 leaves at seedling stage; BS-22, XL22 leaves at bud stage; BS-17, XL17 leaves at bud stage. Fig. S3. Total wax content of cotton leaves of XL22 and XL17 under WW and WD conditions. WW, well-watered; WD, water deficit. SS, the seedling stage; BS, the bud stage. Error bars are standard errors. Values represent the means ± SE, n = 3. Different letters above the bars indicate statistically different from each other as determined by the Student’s t test: p < 0.05. Fig. S4. Heatmap showing the relative expression levels of the 9 selected genes in the two cotton varieties determined by RNA-seq analysis (A) and qRT-PCR (B). Of the selected genes, 4 genes are involved in fatty acid synthesis and 5 genes are related to wax and cutin biosynthesis. The enzymes encoded by these genes are: Very-long-chain enoyl-CoA reductase (GH_D10G2517), Very-long-chain 3-oxoacyl-CoA reductase 1 (GH_D03G1424), 3-oxoacyl-acyl-carrier-protein synthase II (GH_A13G2186), Stearoyl-acyl-carrier-protein 9-desaturase (GH_D10G1329), Delta(12)-fatty-acid desaturase (GH_A11G3588), Probable peroxygenase 4 (GH_A12G1707), Cytochrome P450 86A22 (GH_A08G2095), Very-long-chain -3-hydroxyacyl-CoA dehydratase (GH_D11G2259), Very-long-chain-3-hydroxyacyl-CoA dehydratase 2 (GH_A03G0493). S/W: well-watered plants at the seedling stage, S/D: water deficit plants at the seedling stage, B/W: well-watered plants at the bud stage, B/D: water deficit plants at the bud stage. GhUBI was used as a reference gene. All qRT–PCR reactions were performed in triplicate. Fig. S5. KEGG analysis of DEGs at the seedling stage of XL22 and XL17. A. KEGG categories of DEGs at the seedling stage of XL22. B. KEGG categories of DEGs at the seedling stage of XL17. Fig. S6. Heatmap comparison of DEGs associated with two cotton varieties. A. Venn diagram showing the different KEGGs at the two stages of XL22 and XL17. B and C show the DEGs at the seedling stage of XL22. D and E show the DEGs at the seedling stage of XL17. Red represents high expression, and blue represents low expression. Each row represents a DEG. S/W-22: well-watered XL22 at the seedling stage; S/D-22: water deficit XL22 at the seedling stage; S/W-17: well-watered XL17 at the seedling stage; S/D-17: water deficit XL17 at the seedling stage; B/W-22: well-watered XL22 at the bud stage; B/D-22: water deficit XL22 at the bud stage; B/W-17: well-watered XL17 at the bud stage; B/D-17: water deficit XL17 at the bud stage. Fig. S7. GO analysis and KEGG pathways of the black module and turquoise module genes. A. Analysis of GO enrichment of the black module genes, B. Analysis of GO enrichment of the turquoise module genes, C. KEGG categories of DEGs of the black module, D. KEGG categories of DEGs of the turquoise module. The horizontal axis represents rich factor, the vertical axis represents statistics of the enriched pathways. Circle size represents the number of genes. Red represents high KEGG enrichment, and blue represents low KEGG enrichment.
About this article
Cite this article
Yang, F., Han, Y., Zhu, QH. et al. Impact of water deficiency on leaf cuticle lipids and gene expression networks in cotton (Gossypium hirsutum L.). BMC Plant Biol 22, 404 (2022). https://doi.org/10.1186/s12870-022-03788-2
- Gossypium hirsutum
- Wax components
- Cutin monomers
- Transcriptomic analysis
- Differentially expressed genes
- Co-expression gene network