1. Introduction

Cotton is the most important renewable fiber and edible oil crop worldwide. The genus Gossypium L. contains more than 50 recognized species belonging to eight genome groups, of which four species, Gossypium herbaceum L. (A₁), Gossypium arboretum L. (A₂), Gossypium hirsutum L. (AD₁), Gossypium barbadense L. (AD₂), have been domesticated and cultivated in different regions worldwide (Wang et al., 2012). Cultivation of G. hirsutum, or upland cotton, contributes to over 95% of cotton production worldwide (Iqbal et al., 1997). Cotton traits such as high-yield and agroecological adaptability lead to significant economic benefits. Accordingly, most cotton breeding efforts have focused on improving these traits (e.g., fiber quality, lint yield, boll weight, seed weight) in upland cotton (Brubaker and Wendel, 1994; Ulloa and Meredith, 2000; Zhu et al., 2008).

Of the four cultivated Gossypium species, G. hirsutum shows the highest levels of genetic variation and gene flow (Wendel et al., 1992; Abdurakhmonov et al., 2009). However, some studies have indicated that the germplasm of G. hirsutum currently cultivated lacks genetic diversity and shows little sign of genetic introgression (Esbroeck et al., 1999; Campbell et al., 2010). This is likely because in addition to the initial bottleneck encountered during the domestication process, cotton breeding has frequently involved crossing and re-selection within small sets of breeding materials (May et al., 1995; Bradbury et al., 2007). The narrow genetic base of upland cotton has become a serious concern as limited genetic variation corresponds to limited allelic availability for continued genetic gain (Brown, 1983; Xu et al., 2001; Tyagi et al., 2014).

Cotton yield and fiber quality are typical quantitative traits, the phenotypes of which are the result of quantitative trait loci (QTLs), their interaction, the environment, and the interaction between QTLs and the environment. Understanding how QTLs and environmental conditions interact to influence phenotypes has generated numerous genetic models (Song et al., 2005; Dong et al., 2010). Traditionally, plant QTLs have been identified by linkage mapping, i.e., constructing a high-density molecular linkage map through limited parent hybridization, and localizing specific linkage segments associated with target traits (Fang et al., 2013). This approach was first used over twenty years ago to mine QTLs associated with agronomic and fiber traits of upland cotton (Shappley et al., 1998). In recent years, numerous QTLs in cotton varieties have been excavated through segregation analyses (Shen et al., 2005; Wang et al., 2006; Li et al., 2008). However, because of the limited number of markers used in upland cotton populations, the accuracy of QTLs is relatively low, which impedes further studies using mapbased cloning and marker-assisted selection.

One alternative approach to studying complex quantitative traits and mapping target genes is association analysis (Thornsberry et al., 2001; Flint-Garcia et al., 2005; Shen et al., 2007; Li et al., 2016a, b). Previous studies have used association analysis in upland cotton to map QTLs to agronomic traits, including disease resistance (Mei et al., 2014; Zhao et al., 2014) and seed quality (Liu et al., 2015). Linkage disequilibrium-based association mapping and model-based Bayesian analysis of upland cotton has been used to investigate the population genetics and significant quality traits. For example, association analysis of G. hirsutum accessions has revealed the existence of potential genetic variation for primary fiber quality traits and shown that the development of fiber quality traits is influenced by environmental factors (Abdurakhmonov et al., 2009). Additional studies in cotton have used association analysis to detect associations between elite alleles and fiber quality traits (Cai et al., 2014; Islam et al., 2016; Nie et al., 2016; Ademe et al., 2017). Investigating both population structure and kinship has revealed that in upland cotton interspecific gene combinations appear to improve fiber length traits (Handi et al., 2017; Wang et al., 2017a, b). Despite this progress, association mapping studies have been limited to single traits, e.g., lint yield or fiber quality and/or focused on either genetic structure or association analysis of the cotton quality traits (Islam et al., 2016; Ademe et al., 2017), while ignoring the effects of human transference on genetic exchange of superior cotton accessions.

In this study, we have evaluated genetic structure and gene flow pattern of 285 domestication cotton superior germplasm accessions and performed association analysis of nine yield components and six fiber quality traits using 19 polymorphic microsatellite markers. Our results will provide useful information for understanding the genetic base of upland cotton superior varieties and will also promote future high yield and excellent fiber quality breeding.

2. Materials and methods 2.1. Collection of upland cotton superior accessions

A total of 285 G. hirsutum accessions were selected for this study (Table S1). All selections were derived from five sources: Yangtze River region, China (136 varieties), Yellow River region, China (123 varieties), Northern China (8 varieties), Northwestern China (9 varieties), and the United States (9 varieties) (Table S1).

2.2. Field experiments and trait phenotyping

All 285 accessions of upland cotton were planted in each of four locations that have distinct environmental conditions: Jingzhou, Hubei Province, China in 2015 and 2016 (designated environments E1 and E2, respectively); Jiujiang, Jiangxi Province in China in 2015 and 2016 (designated E3 and E4, respectively); Xinjiang in China in 2015 and 2016 (designated E5 and E6, respectively); and Anyang, Henan Province in China in 2015, 2016 and 2017 (designated E7, E8, and E9, respectively). Each accession was grown in a plot having 40-45 plants in two rows, with 0.10 m between plants in each row and 0.45 m between rows. Field planting followed a randomized complete block design with three replications in each environment. Field management followed conventional standard field practices.

The ten plants in the middle of each row were tagged for scoring and harvesting cotton seed. The yield traits evaluated included growing period, plant height (cm), number of fruit branches per plant, height of first fruit (cm), number of first fruit node, number of bolls per plant, seed butter (g), lint percentage (%), boll weight (g), and unit area yield (g·m^-2). The fiber quality traits evaluated were as follows: mean length of upper half fiber (mm), uniformity index (%), specific breaking strength (cN·tex^-1), micronaire value, and elongation (%).

2.3. Marker screening and genotyping

Total genomic DNA was extracted from collected samples following the methods of Doyle and Doyle (1987) and stored at -20 ℃. Amplification efficiency and polymorphism were randomly tested with 30 randomly selected microsatellite primer pairs (Nie et al., 2016). We found that 19 primer pairs generated polymorphic markers across all 285 G. hirsutum accessions (Table S2). Polymerase chain reaction (PCR) was performed in a volume of 10 μL containing 1 μL of DNA template, 0.3 μL of each primer (1000 ×), 5 μL of 2 × PCR Master Mix and 3.4 μL of ddH₂O under the following conditions: 3 min at 94 ℃, followed by 32 cycles of 30 s at 94 ℃, 40 s at 53 ℃, and 30 s at 72 ℃, and then a final extension of 5 min at 72 ℃. PCR products were firstly separated by 10% polyacrylamide gels and visualized by silver staining. Fragment sizes of each locus were estimated using Quantity One (Bio-Rad Laboratories, Berkeley city, CA, USA) and a 50 bp DNA ladders size standard (pBR322 DNA/MspI marker; Tiangen, Beijing, China). We then verified the polyacrylamide results by regenotyping all samples using an ABI3730 DNA analyzer (Applied Biosystems, Foster City, CA, USA) and fluorescently labeled primers. We performed the PCR amplification using a Veriti 96-Well Thermal Cycler. The upper primers were labeled with 6-FAM, HEX, TAMRA, or ROX (Sangon, Shanghai, China) (Table S2). The PCR reaction included 1 μL DNA template, 0.3 μL each primer (1000 ×), 10 μL2 × PCR Master Mix and 8.4 μL ddH₂O. ABI3730 DNA analyzer (Applied Biosystems) was used to score genotypes and GeneMarkers v. 2.0 was used for binning (Holland and Parson, 2011).

2.4. Data analysis

Popgene v. 1.32 (Yeh and Boyle, 1996) was used to calculate linkage disequilibrium (LD) of all SSR markers (Table S3). Genetic diversity parameters of G. hirsutum from each source region were evaluated per locus using the following descriptive summary statistics: number of alleles (N_A), observed (H_o) and expected (H_e) heterozygosity, and inbreeding coefficient (F_IS) using GenAlEx v. 6.5 (Peakall and Smouse, 2012). CERVUS v. 3.0 (Kalinowski et al., 2007) was used to calculate polymorphismin formation content (PIC). In addition, Arlequin v. 3.5 (Excoffier and Lischer, 2010) was used to test the Hardye-Weinberg equilibrium (HWE) for all SSR markers (Table S4). We also performed a hierarchical analysis of genetic differentiation using an analysis of molecular variance (AMOVA) with 1000 permutations as implemented in Arlequin v. 3.5 (Excoffier and Lischer, 2010). The significance of fixation indices was tested using 10, 000 permutations.

We used Migrate-n v. 3.6.11 (Beerli and Felsenstein, 1999) to investigate long-term effective population sizes (θ) and migration rates (Nm) of G. hirsutum accessions from five geographic regions. We also calculated these measures for two subpopulations (P1 and P2) that had higher likelihood at K = 2 according to our evaluation of population structure. Subpopulation P1 contained 80 lines, including 43 cultivars from Yellow River region, 36 cultivars from Yangtze River region, and one line from Northern China region; subpopulation P2 contained 205 accessions including 80 lines from Yellow River regions, 100 cultivars from Yangtze River region, nine lines introduced from US, andeight lines both from Northern China and Northwestern China regions. Migrate-n uses coalescence theory to model population sizes and migration rates, and mutation models to explain change of alleles at sites over time. Bayesian analysis was run using the infinite allele option, which is recommended when the mutation model is unknown. Uniform prior for θ was set to min: 0.0, max: 100.0, delta: 10.0. Uniform prior for migration was set to 0.0, max: 1000.0, delta: 100.0. STRUCTURE analysis was run using 100, 000 burn-in MCMC iterations, with a length of 1, 000, 000 iterations, and eight replicates per run for K = 1-8 clusters with admixture model (Pritchard et al., 2000). The website program STRUCTURE HARVESTER was used to calculate the optimal value of K (Earl and Vonholdt, 2012) using the delta K criterion (Evanno et al., 2005; Jakobsson and Rosenberg, 2007). The corresponding Q-matrix at K = 2 was obtained for further markeretrait association analysis.

Many cotton traits are complicated quantitative traits; therefore, the most stable markers are those that can be detected at the same time in multiple populations and multiple environments (Li et al., 2016a, b). General linear model (GLM) and mixed linear model (MLM) are two prevalent statistical models in association mapping (Cardon and Palmer, 2003). To generate more accurate correlations with less-inflated type Ⅰ errors (Yu et al., 2006), the MLM (+K + Q) method was employed in the present study (Cardon and Palmer, 2003). Considering the cultivation history of upland cotton and the relatively simple population structure in this panel, GLM (+Q) was also employed, and the results derived from the GLM and MLM were compared. TASSEL v. 2.1 was used to determine relative kinship among the individuals of experimental materials (Yu et al., 2006), and the relative coefficient matrix (K-matrix), which was obtained for subsequent association analysis. General Linear Model (GLM) and Mixed Linear Model (MLM) were used to construct markers-yield and fiber quality traits association tests within TASSEL v. 2.1 (Yu et al., 2006); we set the number of permutations in "define F tests" to 1000.

3. Results 3.1. Phenotypic variation of G. hirsutum accessions grown in different environments

G. hirsutum yield and fiber quality traits displayed broad variation when grown under different environmental conditions (Table 1). For example, the average phenotypic values of ten upland cotton yield traits for accessions grown innine environments are as follows: Growth period (GP) was 131.93 d (range 126.31-142.19 d), plant height (PH) was 93.04 cm (range 58.80-117.97 cm), the number of fruit branches per plant (NFB) was 11.75 (range 8.00-17.39), the height of first fruit (HF) was 18.74 cm (range 14.59-24.33 cm), the number of first fruit node (ND) was 6.74 (range 5.20-8.11), the number of bolls per plant (NB) was 19.22 (range 6.59-33.89), seed butter (SB) was 25.35 g (range 9.87-37.22 g), lint percentage (LP) was 37.61% (range 35.99-41.60%), boll weight (BW) was 5.22 g (range 4.38-6.00 g) and unit area yield (UAY) was 300.36 g m^-2 (range 149.05-541.71 g m^-2). The average phenotypic values of five upland cotton quality traits for accessions grown in nine environments are as follows: the mean length of upper half fiber (LF) was 28.97 mm (range 27.73-29.54 mm), uniformity index (UI) was 84.57% (range 82.20-85.62%), specific breaking strength (BS) was 28.31 cN·tex^-1 (range 26.24-30.39 cN·tex^-1), miconaire value (MV) was 4.98 (range 4.71-5.26) and fiber elongation (FE) was 6.71% (range 6.67-6.76%).

For yield traits, the coefficientof variance ranged from 4.03% (growth period) to 34.19% (unit area yield). Phenotypic variation in fiber quality traits was lower, with the coefficient of variance ranging from 1.46% (fiber elongation) to 12.68% (miconaire value).

3.2. Genetic diversity

We used 19 SSR markers to examine the genetic diversity of G. hirsutum accessions from different geographic regions and their relationships. The expected heterozygosity (H_e) and observed heterozygosity (H_o) ranged from 0.417 and 0.715 in G. hirsutum accessions from Northwestern China to 0.493 and 0.782 in the accessions from the Yangtze River region, with an average of 0.469 and 0.759, respectively. Moreover, the effective number of alleles per locus (N_e) varied from 1.903 in the accessions from Northwestern China to 2.121 in the accessions from the Yangtze River region, with a mean of 2.032 (Table 2). The average polymorphism information content (PIC) was 0.427 (ranging from 0.071 to 0.712) (Table S5). The average fixation indices values (F) were lower than zero for all accessions. These negative values indicate an excess of heterozygotes. The relatively small values of PIC suggest that the genetic diversity of the upland cotton cultivars examined in this study is relatively low.

3.3. Population structure and gene flow pattern

Bayesian clustering analysis revealed that G. hirsutum accessions were not clustered together according to geographic area or pedigree origin. The most likely value showed a much higher likelihood at K = 2, suggesting that the total panel could be divided into two subpopulations, designated P1 and P2 (Fig. S1; Fig. S2, Fig. 1). The P1 subpopulation contained 80 lines, including 43 cultivars (43/123, 35.0%) from the Yellow River region, 36 cultivars (36/136, 26.5%) from the Yangtze River region, and one line (1/9, 11.1%) from the Northern China region. The P2 subpopulation consisted of 205 accessions, including 80 lines (80/123, 65%) from the Yellow River region, 100 (100/136, 73.5%) cultivars from the Yangtze River region, nine lines (9/9100%) introduced from US, and eight lines both (8/8, 100%; 8/9, 88.9%) from the Northern and Northwestern China regions (Fig. S1). The corresponding Q matrix (at k = 2) was further used for markeretrait association mapping (Fig. S3). In addition, based on the results of the relatedness analysis, a K-matrix was also constructed for the association mapping.

Population structure analysis showed that the among-groups component of genetic variance was 0.05%, -98.06% among individuals within groups, indicating an excess of heterozygotes in accessions of the two subpopulations. The within-individuals component was 198.01% (Table 3). These results indicate that the variation among different cotton individuals contributed most to the overall variation. In addition, the Migrate-n analysis of the five geographic groups produced θ and M values greater than zero (Table 4). θ values did not vary among range sectors. Moreover, scaled immigration rates (M) revealed the existence of extensive historical gene flow between all five sectors. Gene movements occurred predominantly from northernto northwestern China, northwestern China to northern China (40.3744 vs 40.1705) and the Yellow River region (33.6081 vs. 21.9127), followed by the Yangtze River region to the US (34.9493 vs. 13.0342). These results illustrate that frequent human-induced gene flow may have occurred between the northern and northwestern China regions, the northwestern China region and the Yangtze River region, as well as into the Yellow River region. Additionally, the population Migrate-n analysis showed that there was a subtler deviation in the direction of migration from subpopulation P2 to subpopulation P1 (11.2123) than from subpopulation P1 to subpopulation P2 (7.8593) (Table 5).

3.4. Association analysis of yield and fiber quality traits

Linkage disequilibrium (LD) tests showed that the level of LD for G. hirsutum accessions from five ecological areas were 43, nine, zero, zero, and zero (at the significance level of p < 0.05). This finding indicates that most SSR markers did not exist linkage disequilibrium (Table S3).

For all 15 agronomic traits, including ten yield component traits and five fiber quality traits, we applied general linear model (GLM) and mixed linear model (MLM) to analyze nine environment datasets derived from the 285 accessions at four locations over two years in Jingzhou, Jiujiang, and Xinjiang as well as three years in Anyang. When we compared the results of the GLM to MLM, we found that a total of 82 makere-trait associations were detected between 17 SSR markers (p = 0.05) and 15 agronomic traits in nine environments (Table S6). The number of SSR markers associated with G. hirsutum yield traits and fiber quality traits are shown in Table 6.

Association analysis further showed that 17 SSR markers are associated with from 3.04% to 22.35% of the phenotypic variation in G. hirsutum accessions (Table S7). Interestingly, the marker NAU5077 was simultaneously associated with 11 traits, including GP, PH, NFB, ND, LP, UAY, LF, UI, BS, MV and FE. Among these traits, UI, BS and MV were closely related to fiber quality (Table 6). The marker NAU4860 was simultaneously associated with NB, SB, BW, LF, UI, BS, MV and FE, three of which are fiber quality traits. Further, NAU5233 was simultaneously associated with GP, UAY, BS and FE; NAU4951 was simultaneously associated with HF, NB, SB and MV. NAU5013 was simultaneously associated with PH, LP and SB; NAU5195 was simultaneously associated with HF, ND and LP. NAU5148 was simultaneously associated with HF, BW and UAY. NAU5260 was simultaneously associated with NB and SB. NAU4932 and NAU4956 were simultaneously associated with GP and FE. NAU5120 and NAU5088 were simultaneously associated with NFB and NB; NAU5017 was simultaneously associated with HF and UI, and NAU5227 was simultaneously associated with LP, LF and MV.

4. Discussion

The geographic distribution of genetic variation in species is significantly associated with their evolutionary potential and future fate (Wendel et al., 1992). Generally, the domesticated crops have less genetic variability than their wildrelatives (Wendel and Cronn, 2003; Cao et al., 2014). Most cotton varieties planted in China were derived from a limited number of founder parents, such as DPL (a cotton germplasm type), Stoneville, King, Uganda, Foster, and Trice (Chen and Du, 2006). Therefore, to create association maps, it is especially critical to select samples that encompass as much genetic diversity as possible. In this study, the population panel consisted of 285 cultivars, including lines from cotton germplasm resources, historical varieties from abroad, multiple lines derived from radiation breeding programs, and some progenies of intraand inter species. Our results showed that the level of diversity in upland cotton varieties was relatively low (PIC = 0.427), with an average number of alleles per locus of 3.2 (ranging from 2.2 and 6.6 alleles/ locus). These results are similar to those detected in the variations analysis of 241 G. hirsutum cotton cultivars (Zhao et al., 2014). The average number of alleles per locus, gene diversity, and PIC in our study were less than those detected in 35 cultivars and eight inbred lines of G. hirsutum from Africa, United States, and Brazil (Lacape et al., 2007). On possible explanation for this difference is that cultivars domesticated directly in a native cotton growing area usually preserve higher levels of polymorphism than those cultivated in a non-native cotton growing areas. In our study, the average genetic diversity and PIC values were higher than those detected in previous research (Qin et al., 2015). Our results indicate that the selected markers have sufficient polymorphic information to reveal the genetic relationships between these upland cotton inbred cultivars.

Although our selection of samples emphasized genetic diversity as much as possible, when we compared the genetic diversity of five ecological groups in our study, we found that genetic differentiation among these ecological areas was still very small. This finding might be explained by frequent gene exchange and germplasm domestication events in different ecotype collections, reflecting the probable extensive exchange of parental lines by breeders. Meanwhile, the low genetic differentiation of upland cotton in five ecological regions is also probably associated with the single origin of its domestication. Over 95% of cultivated cotton crop worldwide is allotetraploid upland cotton, which was possibly initially domesticated in northern Yucatan peninsula (Stephens, 1958; Wendel et al., 2009; Coppens d'Eeckenbrugge and Lacape, 2014; Fang et al., 2017). Therefore, the low genetic differentiation of allotetraploid cotton is likely the result of frequent gene exchange with a restricted domestication source. Similarly, numerous studies have shown that within Gossypium species genetic diversity is low (Abdalla et al., 2001; Iqbal et al., 2001; Rungis et al., 2005; Lacape et al., 2007).

Additionally, population structure is important for explaining the heterogeneity of genetic architecture and is mostly affected by geographic isolation and genetic exchange isolation (Guo et al., 1997; Gutiérrez et al., 2002). We found that instead of being separated in accordance with their geographic origins, all 285 accessions could be classified into two subpopulations. This classification indicates that when upland cotton germplasm is intersperesed or crossbred, genetic exchange may occur frequently, independent of geographic restriction. These results may provide important insight into evaluating the effects of cross-breeding in molecular breeding programs.

Our results provide strong support for the validity of traitassociation results. We identified nine markers that are associated with upland cotton yield traits and two markers that are associated with fiber quality within specific environments. These findings demonstrate that traits related to yield quality in G. hirsutum germplasm have potential genetic variation that may be useful for future breeding programs. We also found that yield and fiber quality traits are correlated with diverse environments. These findings suggest that agriculturalists attempting to select target traits with the same cotton cultivar but in different environments should consider using different practices (Zhang et al., 2005).

The major target traits of cotton during the breeding process are quantitative; thus, phenotypic variation of each trait is directly or indirectly affected by that of other traits (Huang et al., 2018; Keerio et al., 2018; Wen et al., 2018). Numerous QTL mapping studies aimed at improving cotton fiber quality and yield traits have been previously reported (Mei et al., 2013; Adhikari et al., 2017; Wang et al., 2017a, b; Dong et al., 2018a, b). In this study, a total of 16 SSR markers for yield and fiber quality traits were detected. Among these markers, 14 were simultaneously associated with more than two traits, which may have resulted from geneegene interactions or pleiotropism (Zeng et al., 2009; Lehner, 2011). For example, the SSR markers NAU5077 and NAU4860 were simultaneously associated with numerous fiber quality traits, including LF, UI, BS, MV and FE. In addition, SSR markers NAU4951, NAU5013, NAU5195, NAU5148, NAU5260, NAU5088, and NAU5017 were mainly associated with yield quality traits, including HF, NB, SB, PH, LP. These identified associations for different fiber and yield quality traits in domesticated upland cotton, along with those reported in previous studies, add toa rich cluster of yield and fiber quality QTLs (Cai et al., 2014; Adhikari et al., 2017; Ademe et al., 2017; Dong et al., 2018a, b).

In conclusion, we used 19 SSR markers to determine the genetic structure and gene flow patterns of 285 upland cotton accessions. We then identified which markers are associated with agronomic traits in upland cotton. Our results showed that the extensive gene flow occurred in different ecological and geographic regions for crossbreeding. Specific markers identified from association analysis are potential QTLs for the selected traits inselected cotton production regions and can provide more information for markerassisted breeding programs.

Author contributions

ZL designed the work. NZ, TZ, WL, and XM performed the experiments. WL, XZ, XP, YL, KH, WZ, KZ, DY, FZ, and ZR contributed materials/analysis tools. ZL and TZ wrote the manuscript. ZL, XM, DY, and TZ revised the manuscript.

Declaration of Competing Interest

The authors declare that they have no conflict of interests.

Acknowledgments

This research was co-supported by grants from National Key R and D Program for Crop Breeding (2016YFD0100306), National Natural Science Foundation of China (No. 31401431), the Shaanxi Science and Technology Innovation Team (2019TD-012), the Public health specialty in the Department of Traditional Chinese Medicine (Grants no. 2017-66 and 2018-43), and the Open Foundation of the Key Laboratory of Resource Biology and Biotechnology in Western China (Ministry of Education) (Grants no. ZSK2017007 and ZSK2019008).

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2020.03.001.