Phenomics-based GWAS analysis reveals the genetic architecture for drought resistance in cotton.

Drought resistance (DR) is a complex trait that is regulated by a variety of genes. Without comprehensive profiling of DR-related traits, the knowledge of the genetic architecture for DR in cotton remains limited. Thus, there is a need to bridge the gap between genomics and phenomics. In this study, an automatic phenotyping platform (APP) was systematically applied to examine 119 image-based digital traits (i-traits) during drought stress at the seedling stage, across a natural population of 200 representative upland cotton accessions. Some novel i-traits, as well as some traditional i-traits, were used to evaluate the DR in cotton. The phenomics data allowed us to identify 390 genetic loci by genome-wide association study (GWAS) using 56 morphological and 63 texture i-traits. DR-related genes, including GhRD2, GhNAC4, GhHAT22 and GhDREB2, were identified as candidate genes by some digital traits. Further analysis of candidate genes showed that Gh_A04G0377 and Gh_A04G0378 functioned as negative regulators for cotton drought response. Based on the combined digital phenotyping, GWAS analysis and transcriptome data, we conclude that the phenomics dataset provides an excellent resource to characterize key genetic loci with an unprecedented resolution which can inform future genome-based breeding for improved DR in cotton.

Introduction

Cotton (Gossypium hirsutum Linn.) is a relatively drought‐resistant economic crop, and the fibre is a valuable source material for textiles (Chen et al., 2007). Nevertheless, owing to the effects of climate change and extreme lack of fresh water (Li et al., 2009, 2009a,b), drought has been one of the most critical abiotic stresses influencing the growth of crops (Zhu, 2002; Zhu, 2016). This leads to reduced yield in the field, even in cotton (Deeba et al., 2012). Thus, there is an urgent need to improve drought resistance of cotton, not only through the use of irrigation but also through genetic improvement programmes.

Drought resistance (DR), a complex quantitative trait, is controlled by multiple interacting genes to induce morphological and physiological responses (Blum, 2011; Fukao and Xiong, 2013; Guo et al., 2018; Zhou et al., 2007). Plants respond differently at different growth stages under unpredictable hostile environments (Kooyers, 2015). In order to describe the complex mechanism of drought, an index defined as ‘stress elasticity’ has been used in barley, which shows a good correlation with drought tolerance (Chen et al., 2014). In a previous study, natural variation in the gene ZmVPP1 has been found to contribute to drought tolerance in maize seedlings, using survival rate as an indicator of drought tolerance (Wang et al., 2016). Another study in rice led to the identification of OsPP15 using a high‐throughput automated phenotyping platform (Guo et al., 2018). These studies support the idea that plants might have their own phenotypic characters when subjected to water deficit, which would be useful to dissect drought resistance in cotton.

Phenomics represents a new interdisciplinary field that has attracted more attention in recent years, but still lags far behind genomics research (Houle et al., 2010). Conventional phenotyping is labour‐intensive, time‐consuming, costly and low‐throughput, and acquiring accurate phenotypes is one of the major bottlenecks in modern crop breeding (Furbank and Tester, 2011). However, high‐throughput phenomics may avoid these drawbacks through the use of new technologies such as automatic phenotyping systems (Arvidsson et al., 2011; Granier et al., 2006; Nagel et al., 2012; Walter et al., 2007), noninvasive imaging (Al‐Tamimi et al., 2016; Araus and Cairns, 2014; Golzarian et al., 2011; Yang et al., 2014), and even artificial intelligence (Madec et al., 2017; Potgieter et al., 2017; Ubbens et al., 2018). Multidimensional image data based on high‐performance inputting have improved precision and saved time, and can provide scientists with more comprehensive information on plants. Massive plant phenomics datasets can be reanalysed by other users for integrating their omic data to verify and develop novel mechanisms and models (Cooper et al., 2018; Fiorani and Schurr, 2013; Neveu et al., 2019; Tardieu et al., 2017). Integrated analysis of phenomics and genomics has contributed to the discovery of candidate genes in crops in recent years (Dingkuhn et al., 2017; Prado et al., 2018; Spindel et al., 2018). Although phenomics has been applied to many crop and fruit species (Aquino et al., 2018; Sandmann et al., 2018; Singh et al., 2018; Yang et al., 2015), its use in relation to drought stress has been less reported.

In previous work, an automatic phenotyping platform (APP) was utilized to achieve high‐throughput screening of rice and maize for genetic studies and drought response (Guo et al., 2018; Yang et al., 2014; Zhang et al., 2017). In the present study, the APP was expanded for high‐throughput phenotyping of a natural cotton population consisting of 200 representative upland cotton accessions during the seedling stage of development once a week across six time points, either under standard watering or drought conditions. One hundred and nineteen different traits, with some traditional agronomic traits and more novel digital phenotypic traits, were obtained at each time point. Novel digital traits, as well as traditional DR‐related traits, were identified as indicators for cotton DR. Combined with next‐generation sequencing for genotyping, a total of 390 associated quantitative tarit loci (QTLs), including some hotspots, were identified based on phenome and genome data. These QTLs included some previously reported genes: GhRD2, GhNAC4, GhHAT22 and GhDREB2 (Hou et al., 2018; Li et al., 2017), as well as novel genes.

Results

High‐throughput automatic phenotyping platform for cotton

Phenotyping experiments were carried out in Wuhan using the APP (Yang et al., 2014) with a natural population of 200 upland cotton accessions collected from 20 different regions (Fig. 1a; Table S1), grown either under standard watering or drought conditions (Fig. 1b). Each plant was inspected from the twenty side‐view images at six time points with drought conditions at seedling stage (Fig. 1c), producing a total of 1 440 000 images, yielding ca. 1.5TB of data. After image processing, 119 image‐based digital traits (i‐traits) were collected at each time point. These could be divided into two categories – 56 morphological i‐traits and 63 texture i‐traits (Fig. 1d‐e). These i‐traits could be further used in GWAS (Fig. 1f‐g).

Fig. 1

High‐throughput automatic phenotyping platform for cotton. (a) Geographic origin of 200 collected upland cotton. (b) Drought experiments were carried out in the automatic phenotyping platform, plants in control were irrigated regularly by sprinkler facility. (c) Images were collected at six time points either under mild drought, severe drought and rewater treatment during seedling stage. (d) Image processing. Binary image and G, H, I components of plant image were extracted from RGB images. (e) Two types of i‐traits were obtained: morphological traits and texture trait. (f‐g) Combined with genome sequencing, these i‐traits could be used in GWAS.

The 56 morphological i‐traits could be divided into four modules: plant geometry [plant height (PH), plant width (PW), height‐width ratio (HWR), side total projected area – bounding rectangle area ratio (SBR), perimeter‐area ratio (PAR), fractal dimension without and with image cropping (FD1, FD2)], plant compactness [plant density (PD) (i = 1, 2, 3, 4, 5, or 6) and PD derived features (PD_1, PD_2, PD_3, PD_R, PD_K and PD_S), relative frequencies (RF) (0 < j < 21) (j = n) and RF derived features (RF_1, RF_2, RF_3, RF_4, RF_5, RF_R1, RF_R2, RF_R3, RF_11, RF_21, RF_31, RF_R4, RF_K and RF_S)], plant colour [green area (GA), green projected area ratio (GPAR)] and plant biomass [side total projected area (SA)].

The 63 texture i‐traits were computed from green (G), hue (H), and intensity (I) components of the RGB images. For each component, 21 texture i‐traits were calculated, including 6 histogram traits [mean value (ME), standard deviation (SD), smoothness (ST), third moment (TM), uniformity (UF) and entropy (ET)] and 15 traits calculated from grey gradient co‐occurrence matrix (GGCM1‐GGCM15).

In order to evaluate the consistency of image‐based measurement and manual measurement, plant height and fresh weight (biomass) were measured manually immediately after imaging. Strong positive correlations (R 2 > 0.93, P‐value < 0.001 in plant height, and R 2 > 0.93, P‐value < 0.001 in fresh weight) were observed between image‐based measurement and manual measurement (Fig. S1).

Phenotypes of cotton were affected by genotype and environment, especially affected by the latter at the seedling stage, when plants grew and changed rapidly every day, except for compactness traits: RF and middle‐level PD (Table S2). However, the correlation between morphological i‐traits revealed a similar trend in years 2015 and 2017 (Fig. S2a; Table S3), only some of texture i‐traits were inconsistent in different years (Fig. S2b). A high positive correlation (r = 0.81, P‐value < 0.001, in year 2015; r = 0.77, P‐value < 0.001, in year 2017) was found between PH and PW, both i‐traits being representative of plant geometry. Similarly, GA and GPAR showed an obvious correlation. For plant compactness, an obvious correlation was obtained between plant density PD and relative frequencies RF. There were high positive correlations between low‐level PD (PD1 and PD2) and low‐level RF (RF1, RF2, RF3 and RF4), which were the important indices in showing the degree of the compactness at seedling stage for cotton. High correlations were also found for different types of i‐traits, such as GA and SA (Fig. S2).

Variation of i‐traits under drought stress in cotton

To monitor the drought response of cotton during seedling stage (T1–T6) using the APP, we collected images of the 200 cotton accessions, withdrawing water from T2 (mild drought) to T3 (severe drought) time points, and then re‐watering from T4 (early‐rewatering) to T6 (late‐rewatering). Plants of some accessions showed wilting, while plants of some accessions showed little reduced growth at T2 and T3. These i‐traits changed dramatically under drought conditions compared with under normal water condition, especially in geometric‐related and compactness‐related traits.

Coefficient of variation (CV) of drought‐resistant coefficient (DRC) for different i‐traits among 200 accessions was used to determine the discrete degree of observational data. The CV ranged from 0.01 to 1.78 across all morphological traits, showing large differences between different i‐traits. By comparing the three sets of data with control, drought and DRC at six early growth stages, the CV became larger in the group of DRC, especially at severe drought stage (T3) for high‐level compactness indicators (RF20 and PD6) (Fig. S3).

PH and PW served as traditional indicators of drought response, which presented a significant difference between normal and drought groups. Taking two genotypes ZY7 (grouped as a drought‐sensitive accession) and ZY116 (grouped as a drought‐resistant accession) as examples, PH and PW reduced significantly under severe drought condition for ZY7 compared with ZY116, and an obvious symptom of wilting could be seen at T3 under drought condition (Fig. 2a‐c). After mild drought stress at T2, the increment of plant biomass was restricted and this trend continued with time and was maintained two weeks after re‐watering (Fig. 2d). The GPAR showed little variation during seedling stage, but a large fluctuation was seen in compactness traits, such as RF20 and PD6. In total, ninety morphological and texture i‐traits showed a significant difference between treatment and control panels after drought stress (Paired t‐test, P‐value < 0.001; Fig. S4) and most i‐traits were in normal distribution.

Fig. 2

Variation of i‐traits under drought stress during seedling stage. (a) Growth status images of two representative cotton plants (ZY7 and ZY116) at seedling stage (T1–T6) under normal water and drought conditions, respectively. (b) Statistical data for two geometric traits (plant height and plant width) of ZY7 and ZY116 at seedling stage. Prefix with C_ and D_ stand for control and drought treatment, respectively. (c‐d) Statistical analysis for PW, PH and SA at group level. Error bars are calculated by standard deviation (s.d.) based on three biological replicates, statistical significance was determined by a two‐sided t‐test: * P‐value < 0.05, ** P‐value < 0.01.

Based on data from selected morphological i‐traits (PD1‐6, PAR, PW, PH, FD1, FD2, SBR, SA, HWR, GA and GPAR) under both standard watering and drought conditions, the population of 200 upland cotton accessions could be divided into three sub‐groups, whereby sub‐groups I to III represented the cotton accessions that were affected severely or slightly by drought stress (Fig. S5). This result indicated that the vegetative growth of cotton plants was clearly restricted by water deficit, and variable responses to drought stress were observed among the population.

Novel i‐traits could be used as indicators for DR

In addition to traditional agronomic traits, novel i‐traits might be used to reflect the drought response. Compactness, which could be described by PD and RF, was an index that reflects the compactness of branches and leaves of plant. The plant density PD1‐6 stood for the ratio of six shading degrees, where PD1 and PD6 value represented the ratio of lowest and highest shading degree, respectively, and six features derived from PD [PD_1, PD_2, PD_3, PD_R, PD_K and PD_S] were used to give a comprehensive interpretation for drought response. Comparing two representative individuals, the distribution of PD1‐6 was basically the same of ZY116 under both drought and control conditions, but the distribution regularities of ZY7 were completely different whereby the percentages of PD6 and PD5 were obviously higher under drought stress than control (Fig. 3a‐b). Under severe water deficit, the treatment group showed an obviously higher PD6 than control (P‐value < 0.05) at T3 (Fig. 3c), so that drought‐sensitive and drought‐resistant individuals could be identified easily using this index.

Fig. 3

The plant compactness trait of plant density (PD). (a) Original and pseudo‐colour images for cotton plants. a‐d The original images of ZY7 and ZY116 under control and drought stress. e‐h The right pictures illustrate the corresponding plant density images. (b) Plant density in different individuals. ZY7_D has a different percentage distribution than ZY7_C, while the distribution characteristics of ZY116_C and ZY116_D are similar. The suffix of _C and _D stands for control and drought treatment, respectively. (c) Scatter distribution map of indicator of PD6 for cotton population at different seedling stages (T1‐4) under control (blue filled dot) and drought stress (red filled dot).

RF1‐20 represented twenty degrees of distribution for plants through twenty different angles, in which RF1 and RF20 represented the proportion of the lowest and highest compactness degree, respectively, and 14 traits derived from RF [RF_1, RF_2, RF_3, RF_4, RF_5, RF_R1, RF_R2, RF_R3, RF_11, RF_21, RF_31, RF_R4, RF_K and RF_S] were also used to evaluate the degree of drought response. RF_R4 was found to be effective in discriminating drought‐resistant individuals and drought‐sensitive individuals. For drought‐sensitive individuals, they showed a higher RF_R4 value (RF_R4 > 0.2) than drought‐resistant individuals (RF_R4 < 0.1) (Fig. 4a). To study the distribution of RF1‐20 of two representative plants of ZY7 and ZY116, a similar distribution was found for both individuals under normal water condition, but a large difference was observed in the distribution of ZY7 and ZY116 under drought conditions (Fig. 4b‐c). These results suggested that RF might be another useful i‐trait to indicate drought response in cotton.

Fig. 4

The plant compactness trait of relative frequency (RF) and texture i‐traits (G). (a) RF_R4 for sensitive and resistant cotton individuals. a‐h ZY7, ZY283, ZY301, ZY303, ZY9, ZY116, ZY168 and ZY347 under normal water condition; i‐p ZY7, ZY283, ZY301, ZY303, ZY9, ZY116, ZY168 and ZY347 under severe drought condition. (b) The distribution of RF1 to RF20 for a resistant plant of ZY116. (c) The distribution of RF1 to RF20 for a sensitive plant of ZY7. (d) Statistical analysis for ET_G at group level. Error bars are calculated by standard deviation (s.d.) based on three biological replicates, statistical significance was determined by a two‐sided t‐test: ** P‐value < 0.01. (e) The value of DRC of three representative cotton accessions (ZY7, ZY116 and ZY168). Error bars are calculated by standard deviation (s.d.) based on three biological replicates.

Some texture i‐traits also could be used as indicators for drought response in cotton. Taking G as an example, G was the green component that to some extent reflects the chlorophyll content of plants. There were 21 traits in G texture. The DRC value of some texture i‐traits, such as ET_G, declined under drought conditions (Fig. 4d). Based on the above study, we take three representative accessions for further analysis. The results showed that the DRC value decreased under severe drought treatment to a lowest DRC value of ET_G at 0.07 in the drought‐sensitive accession ZY7, which suffered a more serious impact than two drought‐resistant accessions ZY116 (DRC = 0.2) and ZY168 (DRC = 0.47). Moreover, the DRC value of ET_G increased to a relatively high‐level after re‐watering in the highly drought‐resistant accession ZY168 (Fig. 4e).

These results show that dynamic growth is revealed by either traditional indicators or novel i‐traits under drought treatments, providing a useful basis by which to find DR‐related genes using GWAS.

Candidate DR‐related genes based on GWAS and RNA‐Seq

To analyse the molecular basis of DR‐related i‐traits, the cotton accessions were sequenced. A total of ca. 3.8 Tb of sequence data were integrated, with an average depth of ca. 7.4× (Table S1). These data were mapped against the TM‐1 genome to identify genomic variants, and a total of 12 724 791 SNPs were detected (Table S4). After filtering scaffold and MAF (minor allele frequency < 0.05) SNPs, 2 060 458 SNPs were used in factored spectrally transformed linear mixed models (FaST‐LMM), and the strict significance P‐value was set at 2.43 × 10−8 by Bonferroni correction. Population structure was calculated as the covariance to control the false positives (Fig. S6a). We explored the phylogenetic relationships among the 200 cotton accessions, using whole‐genome SNP analysis, and those cotton accessions were divided into four groups, as supported by population structure analysis (Fig. S6b).

GWAS was performed across those cotton accessions for the value of RCRW (recovery capability after re‐watering) of SA and DRC of 119 i‐traits at different time points during seedling stage. A total of 622 unique SNPs were detected and with the smallest number of unique significant SNPs detected on Chromosome (Chr) A05 (Table S5). Some loci identified by GWAS associated with some morphological and texture i‐traits, as shown in Fig. S7a. Among them, the signal of one SNP peak in a region on Chr A04 showed gradual strengthening by GWAS analysis using SA at T2 and T3, and finally exceeded the suggested threshold value at T4 (Fig. S7b). Moreover, this SNP peak was found to be co‐located by GWAS analysis using three i‐traits (PH, PW and SA; Fig. 5a).

Fig. 5

Genome‐wide association studies. (a) Manhattan plots in Chr A04 for three traits of PH, SA and PW. The SNP that arrow points to as the same SNP described in the context. (b‐c) Manhattan plots and QQ‐plots for four representative novel i‐traits. The QTLs including reported genes were marked by dotted boxes. Horizontal red full line indicates the genome‐wide significant threshold (−log10 P = 7.6).

To find candidate genes around peak SNPs on different chromosomes, we needed to narrow down QTL regions. In our study, LD decay distance was ca. 500 kb using all SNPs (Fig. S8). A total of 390 candidate loci were detected by GWAS with those i‐traits (Table S5), and 71 loci could be identified in both two years. Among them, 252 QTLs were associated with morphological i‐traits, with a similar distributed amount to A and D subgenome (132 versus 120; Fig. S9 and Table S5), although the size of the A subgenome (1477 Mb) is nearly twice that of the D subgenome (831 Mb) (Zhang et al., 2015).

RNA‐Seq was performed using a highly drought‐sensitive cotton accession ZY7 and a relatively drought‐resistant cotton accession ZY168 to differentially analyse the candidate genes at those associated loci. A total of 16 827 DEGs were found to be differentially expressed under drought treatments (P‐adjust < 0.01, log2FoldChange > 1). Among them, 8074 DEGs were identified in both genotypes, with 2212 genes differentially expressed (P‐value < 0.05) between the two. In addition, 5090 DGEs were identified only in ZY7, while 3663 DEGs were identified only in ZY168 (Table S6).

Fifteen QTLs within 167 DEGs through four i‐traits of PD, RF, G and SA were selected for further analysis (Table S7). For these QTLs, QTL3 could be identified by GWAS with both RF and G. Four genes (Gh_D11G2874, Gh_D11G2876, Gh_D11G2884 and Gh_D11G2885) were overlapping in QTL1 and QTL14 by GWAS with different i‐traits. Some previously verified DR‐related genes were co‐localized within these regions. For example, GhNAC4 and GhRD2 were associated with the derived features of PD; and GhDREB2 and GhHAT22 were detected as the features of RF (Fig. 5b‐c). However, more DEGs with homologous annotation of some reported stress‐related genes in other plants were present at these loci, such as Gh_D11G2863 (ethylene responsive element binding factor 3), Gh_D13G2394 (tonoplast intrinsic protein 1;3), Gh_D13G2393 (MYB domain protein 48), Gh_A12G1594 (WRKY DNA‐binding protein 72), Gh_D11G2887 (transducin family protein/ WD‐40 repeat family protein), Gh_D11G2899 (receptor‐like protein 50) and Gh_D10G1405 (Calcium‐binding EF‐hand family protein; Table S7). Moreover, there were some DEGs without available annotation even for the homologous genes in Arabidopsis within the associated loci, such as Gh_D13G2320, Gh_D13G2344, Gh_D10G0952, Gh_A04G0377, Gh_A04G0378, Gh_A12G1597, Gh_A12G0141 and Gh_A12G0142.

Identification of two novel DR‐related negative regulative genes

As mentioned above, an associated locus on Chr A04 was identified by several i‐traits, such as SA (Fig. 6a). Two closely linked sub‐regions were presented in this region (Fig. 6b). Other than some annotated DEGs, two previously unreported tandem repeat genes, Gh_A04G0377 and Gh_A04G0378, named as drought negative regulators (GhDNRs), were present in one sub‐region (Table S7). These two genes each had two exons (Fig. 6c). The genotypes were either CC or TT alleles in the promoter regions of these two genes, and a significant difference in SA was found in two sets of plants – those that had the ‘TT’ allele show a higher biomass than the ‘CC’ variant under drought stress (Fig. 6d). The relative expression level of Gh_A04G0377 showed a higher up‐regulation in the drought‐sensitive accession ZY7 than in the drought‐resistant accession ZY168 under drought stress (Fig. 6e). These results suggested that this candidate gene could be transcriptionally induced by drought and repressed under normal water condition in the cotton.

Fig. 6

Identification of the DR genes on Chr A04. (a) Manhattan plot of the trait of SA. The red line indicates the threshold of −log10 P = 7.6, the red arrow points to the lead SNP of A04_12697203. (b) LD heat map of surrounding the peak on Chr A04. Two interlocking regions within 472.3Kb are divided in the green dotted triangle. (c) Gene model of Gh_A04G0377 and Gh_A04G0378. (d) Boxplots for SA, based on the genotypes of SNP A04_10697203, which the red arrow points to in a. (e) Relative expression of Gh_A04G0377 in the leaf of ZY7 and ZY168 under control and drought stress. Values with different letters (a–c) indicate significant differences (Duncan’s multiple comparisons, P‐value < 0.05) between each sample. (f‐g) qRT‐PCR analysis confirmed that Gh_A04G0377 and Gh_A04G0378 were silenced in the TRV:GhDNRs plants. (h) VIGS of GhDNRs (TRV:GhDNRs) in drought‐sensitive accession ZY7, and ZY7 with TRV:00 was used as control. (i‐j) Indicators of SA, manually measured PH, defoliation rate measured in control and silenced plants. (k) Soluble sugar content detected in leaves under drought stress. Error bars were calculated by S.D. based on three biological replicates, statistical significance was determined by a two‐sided t‐test: * P‐value < 0.05, ** P‐value < 0.01.

Further functional validation of GhDNRs in cotton was carried out by virus‐induced gene silencing (VIGS) technology. The drought stress was applied to the cotton plants at the seedling stage, with the two tandem repeat genes of Gh_A04G0377 and Gh_A04G0378 both silenced in the seedlings (Fig. 6f‐h). The image‐based biomass (SA) of TRV:GhDNRs plants was higher than the control TRV:00 plants under drought stress (ASWC, ca. 1%; Fig. 6i), and the manually measured PH was higher in TRV:GhDNRs plants than the control TRV:00 plants (Fig. 6j). The defoliation rate was decreased to 13.3% in the silenced plants compared with the 93.3% in control (Fig. 6k). Delayed leaf defoliation of TRV:GhDNRs plants was also observed in drought‐sensitive plants (Fig. S10). TRV:GhDNRs plants accumulated more soluble sugar than control plants, which might contribute to drought tolerance (Fig. 6l). On the other hand, we found that silenced GhDNRs in the drought‐resistant accession ZY168 showed only small differences in plant height and less difference in leaf defoliation in TRV:GhDNRs plants compared with control plants (Fig. S10).

Discussion

The accurate identification and description of crop phenotypes are a prerequisite for a deep understanding of the relationship between phenotype and genotype as an important aim in modern biology. There is also an urgent need to cultivate new breakthrough crop varieties and ensure national food production and food safety. Compared with the rapid development of genomic technology, the current methods for phenotyping are labour‐intensive, costly, time‐consuming and low‐throughput, so becoming a major bottleneck in crop improvement programmes (Chen et al., 2014; Furbank and Tester, 2011). Here, we utilized efficient crop phenotyping techniques and acquired multidimensional i‐traits for above‐ground parts based on images generated by a high‐throughput APP (Fig. 1), which has proved to be an efficient system in plant phenotyping (Guo et al., 2018; Yang et al., 2014; Zhang et al., 2017).

Drought resistance is a complex quantitative trait. Plants have evolved drought adaptation in part via changes in morphological phenotype (Hu and Xiong, 2014). With the increasing amount of phenotypic data, phenomics has become a mode of investigation and several important phenotype databases are integrated for public analysis (www.plantphenomics.com; http://www.plantphenomics.org.au/PODDProject). The problem remains of how to filter out these redundant i‐traits and pick out the data that meet our needs is another problem (Chen et al., 2014; Schadt et al., 2010). In our study, two kinds of i‐traits were identified to study the drought response of cotton. Using the index of DRC, we studied the correlation among all morphological and texture i‐traits to explain the change regularity of traits under drought conditions, which helped us to screen out some important candidate traits, such as PH, SA, PD and RF.

Plant height and some other traditional traits are vital indices for cotton DR (Fig. 2). In the current study, we paid more attention to novel image‐based indicators (PD, RF and G) to explain the DR of cotton in different dimensions. Drought‐sensitive individuals presented a higher percentage of PD6 under drought stress (Fig. 3). RF1 to RF20, and their derived feature from RF could also be used to evaluate DR and allow visualization of the drought‐resistant individuals and drought‐sensitive individuals (Fig. 4a‐c). One trait in texture i‐traits was also selected to display the DR response (Fig. 4d‐e).

We performed GWAS on a population of 200 upland cotton accessions with the phenotypes collected from the above APP for DR during seedling stage. We not only co‐located some previously reported DR‐related genes in cotton (Fig. 5), but also found some proposed meaningful stress‐related genes, such as Gh_D13G2393, and Gh_A12G1594 (Table S7). Moreover, due to the high positive correlation between PH, SA and PW, a QTL was co‐located on Chr A04. Within this region, two previously unreported tandem repeat genes, Gh_A04G0377 and Gh_A04G0378, were found to relate to drought stress.

To the best of our knowledge, this is the first study to conduct GWAS based on high‐throughput APP derived multiple i‐traits in a cotton drought study. During seedling growth, dynamic DR was revealed in cotton, which might be influenced by different genes. Among a large population, drought‐sensitive and drought‐resistant individuals presented different growth models. Yield and quality of fibre were not investigated in this study, which are the two important traits for the value of cotton. Accurate phenotypic identification data for field‐grown plants was hard to obtain because of exposure to a range of biotic and abiotic stresses. Nevertheless, we associated high‐density resequencing information with abundant phenotypic data and identified novel genes related to DR, as well as previously reported relevant genes (Table S7). Some candidate genes were verified as being differentially expressed under drought stress, and novel genes, GhDNRs were shown to negatively regulate drought tolerance at physiological and transcriptional levels through VIGS experiments, or gene expression tests that showed the drought‐sensitive cotton accession exhibited more significant up‐regulation of relative expression of the genes (Fig. 6). The GhDNRs genes might control the biomass and defoliation under drought stress, but the deeper mechanisms of drought tolerance remained to be clarified.

Materials and methods

Plant materials and experimental design

A natural population of 200 representative upland cotton accessions was planted in APP facilities (Yang et al., 2014) at the National Key Laboratory of Crop Genetic Improvement in Huazhong Agricultural University, Wuhan, Hubei Province, P.R.China (30°28′27″N, 114°21′1″E, 30 m elevation) with two year replicates and three replicates each year, either under standard watering or drought conditions during seedling stage. Among them, 186 cotton accessions were collected from the major Chinese cotton cultivation regions: the northwestern inland region (NIR), the northern specific early maturation region (NSEMR), the Yellow River region (YRR) and the Yangtze River region (YtRR), another 14 varieties were introduced from America and Soviet Union (AS), as described previously (Huang et al., 2017; Wang et al., 2017). Information on the cotton accessions is provided in Table S1. Seeds were sown in late April (Table S8). Plants were grown under controlled greenhouse conditions and were watered by an automatic irrigation system at germination stage, and then manually watered with about 900 mL/pot once a week three days before images capturing. Fertilization was carried out at sowing, and a compound boron fertilizer was added before T5 (60 kg of water + 100 g of boric fertilizer + 370 g of carbamide, fully dissolved and 150 mL applied per pot).

Cotton plants were subject to automatic phenotyping at the three‐leaf stage (ie 30 days after sowing in year 2015, and 32 days after sowing in year 2017). All the cotton individuals were screened at a one‐week interval to generate six time points (as T1 to T6) during the seedling stage. Drought stress was conducted by withdrawing water from the plants (as D) at T2 (mild drought) to T3 (severe drought) time points, while control plants were well watered. Absolute soil water content (ASWC) was used to evaluate the degree of drought in cotton [the ASWC of mild drought (T2) was ca. 5%, and 1% for severe drought (T3), while ca. 15% for the parallel control]. The experiments were conducted on three biological replicates. Each plant grown in a pot with 4.5 kg soil was used as a biological replicate. In total, 1200 plants were used and imaged in this study. The line interval size was 400 mm, and individual plant interval size was 250 mm with one line remaining empty between two planted lines. A small population of ten representative varieties was planted for manual measurement and digital modelling. The sowing dates, watering dates and detection dates for two years were provided in Table S8.

Phenotyping and image‐based digital traits (i‐traits) extraction

The details of the phenotyping facility were described previously (Yang et al., 2014). The automatic phenotyping of cotton was conducted as shown in Fig. 1. RGB images of the cotton plants were acquired using a modified high‐throughput rice phenotyping facility (HRPF). Twenty side‐view images (2452 × 2056 pixels) were taken from different angles at one time point for each plant. Image processing and features extraction were accomplished using LabVIEW as described previously, with minor modifications (Yang et al., 2014). Manual measurements were conducted at each time‐point for digital modelling. After image processing, 119 image‐based digital traits (i‐traits) were computed, including 56 morphological i‐traits and 63 texture i‐traits.

The textures calculated from green (G), hue (H), and intensity (I) component were denoted as G texture, H texture and I texture. For instance, mean value in G texture was denoted as ME_G. The detailed definition and equation of all the i‐traits were shown in the File S1.

Analysis of DR‐related i‐traits

CV using value of control (suffix with _C), drought (suffix with _D) among 200 accessions for different i‐traits was used to determine the discrete degree of observational data. We focused on mild drought condition (T2), severe drought condition (T3) and early‐rewatering treatment (T4). The changes observed in some i‐traits on these three time points were expected to reflect the drought response of cotton. These dynamic changes were investigated in two ways: (i) i‐trait changes following mild to severe drought conditions (such as traditional indicator, PH and some novel i‐traits as defined below); (ii) i‐trait changes representing phenotypic recovery following serious water deficiency (such as SA). DRC (ratio between i‐traitdrought and i‐traitcontrol) and RCRW [(i‐trait_T4 ‐ i‐trait_T3)/i‐trait_T3] were used to evaluate the phenotypic response. Analysis of variance (ANOVA) across two years was performed using QTL IciMapping software (Meng et al., 2015).

The DRC of all traits at T2, T3 and T4 was used to analyse the performance of the accessions under mild drought, severe drought and rewatering conditions. RCRW was defined to show the recovery capability by re‐watering after drought treatments. The values of DRC (suffix with _R) were also used to analyse the correlations among the morphological i‐traits and texture i‐traits, respectively. The distribution of all i‐traits was tested by Shapiro‐Wilk test.

Genome‐wide association studies for DR‐related i‐traits

Of the 200 cotton accessions, 170 were sequenced in a previous study (Wang et al., 2017) and another 30 accessions were sequenced in 2017 (Table S1). The allotetraploid cotton genome (Gossypium hirsutum L. acc. TM‐1) was used as reference genome (Zhang et al., 2015). Clean paired‐end reads (PE150) were mapped to the TM‐1 genome using BWA software (Li and Durbin, 2009) with the default parameters. Picard‐tools in Genome Analysis ToolKit (GATK) (McKenna et al., 2010) software were used to sort BAM format files, remove duplicated data and realign InDel reads. SNPs and InDels calling were calculated separately by two different softwares: GATK and SAMtools (Li et al., 2009, 2009a,b), and linkage disequilibrium (LD) for all pairs of SNPs was calculated using PopLDdecay with the default parameters (Zhang et al., 2019).

GWAS was performed with FaST‐LMM (version 2.02) (Lippert et al., 2011). The RCRW of SA and DRC of all i‐traits was used for GWAS. All filtered SNPs used in association analysis were compressed into binary formats using PLINK software (version 1.90b6.8) (Purcell et al., 2007). The variant annotations were annotated using the software SnpEff (version 4.3) (Cingolani et al., 2012). All SNPs were categorized as being in intergenic regions, upstream and downstream regions, in exons or introns by default parameters. The significant association threshold values were set at −log10 (0.05/n) (n, total filtered SNP number). The significant association regions were manually verified from the aligned resequencing reads against the TM‐1 genome with SAMtools.

Differential expression analysis and validation of candidate genes

Leaves of the accession of ZY7 (drought‐sensitive) and ZY168 (drought‐resistant) from treated and parallel control plants were sampled with two biological replicates for RNA‐Seq analysis at severe drought stage (T3) just as for phenotyping. Total RNA was extracted using the TG‐DP441 RNA Kit, according to the manufacturer's instructions. A total of 2 μg RNA were used for library construction, and RNA sequencing was performed on an Illumina Novaseq™ 6000 system. Clean reads (PE150) were mapped to the TM‐1 genome with Hisat2 (version 2.1.0) (Kim et al., 2015). The expression level of each gene was determined by StringTie (v1.3.3b) (Pertea et al., 2015), and the differentially expressed genes (DEGs) (P‐adjust < 0.01, log2FoldChange > 1) were calculated by DEseq2 R package. t‐Test was implemented to check the significance between the changes of transcript abundance between the two cotton accessions under drought conditions. qRT‐PCR was conducted to verify the expression pattern of candidate genes.

VIGS was used to validate the candidate genes as described previously (Gao et al., 2011). The coding region for Gh_A040377 was cloned into the TRV2 (tobacco rattle virus) vector. The constructed vectors were transformed into A.tumefaciens GV3101 by electroporation and agroinfiltrated into the cotyledons of drought‐sensitive accession ZY7 and drought‐resistant accession ZY168 through syringe inoculation when its two cotyledons had fully spread (13 days after sowing). qRT‐PCR was conducted with the ABI 7500 system (Applied Biosystems, Foster City, CA) to verify the expression of the candidate genes, with the internal control GhUBQ7. The 2−ΔΔCT method was used to present relative changes in gene expression levels (Schmittgen and Livak, 2008). Primers are listed in Table S9. Drought stress was applied to the plants by withdrawing water after seed germination in water‐saturated soil. After plants geminated with no watering for 30 days (ASWC, ca. 1%), the phenotypes were collected. Leaves were sampled to verify the silence effect when the second leaf had fully emerged. Nondestructive biomass (SA) was extracted just as for phenotyping. The plant height and defoliation of leaves were measured manually. Soluble sugar was determined by Anthrone‐Sulfuric acid colorimetry (Yemm and Willis, 1954). Significant test was performed by t‐test (P‐value < 0.05).

Conflict of interest

The authors have declared that no conflict of interest exists.

Author contributions

X.Y. and X.Z. conceived the project and designed the experiments; B.L. conducted the phenotyping and GWAS experiments; W.S. and D.W. participated in the phenotyping; W.Y. participated in designing the APP platform; G.C. contributed to the cotton cultivation; Z.L. and Y.Y. provided the seed materials; L.D. extracted imaging traits from raw images; L.C. conducted the VIGS experiments; M.W. offered writing advice; B.L. wrote the manuscript and managed the main data analysis, including statistical analysis, figure and table design; X.Y. revised the manuscript.

Figure S1. Correlation between manual measurement and image‐based measurement in trait of plant height and biomass.

Figure S2. Correlation analysis for morphological and texture i‐traits in two years, respectively.

Figure S3. Coefficient of variation (CV) among 200 accessions for three sets of at seedling stage.

Figure S4. T‐test between treatment and control panels at seedling stage.

Figure S5. Three sub‐groups i‐traits of 200 upland accessions.

Figure S6. Structure analysis of 200 upland cotton accessions.

Figure S7. GWAS results by morphological and texture i‐traits.

Figure S8. Linkage disequilibrium (LD) decay distance at whole genome level.

Figure S9. Distribution of QTLs associated with morphological traits in 26 chromosomes in cotton.

Figure S10. VIGS of GhDNRs (TRV:GhDNRs) in drought‐resistant accession ZY168 under drought stress.

Click here for additional data file.

Table S1. Cotton resources information and resequencing information used in this study.

Table S2. ANOVA for morphological traits across years in 2015 and 2017.

Table S3. Correlation analysis of morphological traits over two years.

Table S4. Variants rate details.

Table S5. Information of associated loci.

Table S6. DEGs in two cotton accessions under drought stress.

Table S7. QTLs with DEGs in four i‐traits.

Table S8. Date of two‐year seedling experiment of cotton.

Table S9. List primers used in this study.

Click here for additional data file.

SupFile S1. The image processing and i‐traits extraction.

Click here for additional data file.