Abstract
Determining the genetic rearrangement and domestication footprints in Gossypium hirsutum cultivars and primitive race genotypes are essential for effective gene conservation efforts and the development of advanced breeding molecular markers for marker-assisted breeding. In this study, 94 accessions representing the 7 primitive races of G hirsutum, along with 9 G hirsutum and 12 Gossypium barbadense cultivated accessions were evaluated. The genotyping-by-sequencing (GBS) approach was employed and 146 558 single nucleotide polymorphisms (SNP) were generated. Distinct SNP signatures were identified through the combination of selection scans and association analyses. Phylogenetic analyses were also conducted, and we concluded that the Latifolium, Richmondi, and Marie-Galante race accessions were more genetically related to the G hirsutum cultivars and tend to cluster together. Fifty-four outlier SNP loci were identified by selection-scan analysis, and 3 SNPs were located in genes related to the processes of plant responding to stress conditions and confirmed through further genome-wide signals of marker-phenotype association analysis, which indicate a clear selection signature for such trait. These results identified useful candidate gene locus for cotton breeding programs.
Introduction
Cotton is the most important fiber crop in the world and comprises 52 Gossypium species including 7 allotetraploid species (amphidiploids [AD] genome [2n = 52]) and 45 diploid species (2n = 26).1,2 The allopolyploid species, Gossypium hirsutum L (AD1 genome) and Gossypium barbadense (AD2 genome), have been domesticated and G hirsutum cultivars alone account for more than 90% of global cotton fiber production.3-5 However, the high genetic similarity among G hirsutum accessions has hindered opportunities for breeding new cotton cultivars with improved agricultural traits such as higher yield, ease of harvest, and stronger resistance to pest, diseases, and environmental stresses.6-10
Seven genetically related accessions of G hirsutum, including Latifolium, Palmeri, Marie-Galante, Richmondi, Yucatanense, Morrilli, and Punctatum, have been identified based on their locations of origin. 11 These races have distinct characteristics that are common to wild cotton but not to cultivated G hirsutum, such as sensitivity to a short-term light cycle, greater disease resistance and drought tolerance, hard seed coats, and variable seed size, and are genetically compatible with domesticated cottons. All of these advantageous traits can be potentially applied to improve cotton yield and quality as well as tolerance to environmental stresses.8,12-14
In the past, the classification of G hirsutum has been primarily based on morphology, geographical distribution, and cytological markers.15,16 The classification of G hirsutum has now advanced significantly as a multitude of molecular markers have been identified,17,18 but simulation and empirical studies have shown that simple sequence repeat (SSR) markers are likely to result in a significant downward bias for FST estimation due to the mutational characteristics of highly polymorphic microsatellites.19-21 GBS (genotyping-by-sequencing) is a practical and low-cost single nucleotide polymorphism (SNP) marker identification platform which can be utilized for genetic variation screening, genome-wide association analyses, and genetic recombination studies. Application of a large number of genome-wide markers for genotyping across multiple populations enables the establishment of the adaptations that have taken place during evolution and the detection of novel trend during natural selection.22-25
The primary achievement of this study is the establishment of a fine scale genome-wide map of the distributions of SNPs and the determination of the phylogenetic relationships of 115 cotton genotypes including 94 G hirsutum primitive race accessions and 21 domesticated cotton cultivars. Selection-scan analyses and genome-wide association study (GWAS) signals were conducted, based on the phenotype association analysis, to correlate the linkage between the molecular markers for early seedling development with the evolutional and domestication of G hirsutum. The results obtained from this study will facilitate future investigations on the genetic structure of G hirsutum races and expand the marker resources available for breeding programs.
Materials and Methods
Plant materials and phenotypic evaluations
The study evaluated 115 cotton genotypes, including 94 accessions representing 7 G hirsutum primitive obtained from Wild Cotton Nursery located in Sanya City, Hainan Island, and supervised by the Institute of Cotton Research (28 Latifolium, 16 Marie-Galante, 14 Morrilli, 19 Punctatum, 8 Richmondi, 7 Palmeri, and 2 Yucatanense accessions; Supplemental Table S12), 9 G hirsutum cultivars were included to represent modern domesticated upland cotton genotypes, and 12 of the G barbadense cultivars as an outgroup.
The 115 genotypes were planted in the field at the Institute of Cotton Research of Chinese Academy of Agricultural Scientists in Sanya, in 5 m plots with 3 replications. All samples were planted on April 23, 2016, and the field was managed according to traditional production practices. During the harvest season, the number of fruiting branches and the number of bolls per plant were measured on 10 plants and then averaged. The single boll weight (average weight of 30 bolls), lint (lint weight obtained from the 30 bolls/weight of 30 bolls (g) × 100), and seed index (weight of 1000 seeds) were measured. Fiber quality (fiber length, uniformity index, strength, micronaire value, and elongation) were measured by HVI instrument. For the germination test, 30 de-linted seeds for each accession that had been stored for 3 months were placed in the germination box with 3 replicates, then cultured in a germination chamber at 28°C, and 10-hour daylight. The germination rate is the percentage of seeds that germinate at 8 days, the germination potential is the number of germinated seeds/30 seeds; the seedling weight, embryonic axis length, and root length were measured on 5 seedlings (Supplemental Table S1).
To determine phenotypic differences between G hirsutum races and G hirsutum cultivars and their association with genetic structure, phenotypic data comprising 15 morphological traits were subjected to principal components analysis (PCA) and agglomerative hierarchical cluster (AHC) analysis. In addition, differences between cultivars and G hirsutum races were determined by subjecting all morphological traits to analyses of random variance followed by the Tukey honest significant difference post hoc test at a significance level of P < .05. All calculations were performed with XLSTAT version 2013.
DNA extraction, GBS library preparation, sequencing, and data analysis
Young leaf tissues from a single plant for each genotype before flowering were collected, and DNA was extracted using a Qiagen DNeasy Plant Mini Kit following the manufacturer’s instructions. The concentration of DNA was determined by fluorimetry (Life Invitrogen Qubit 3.0, Qubit 3.0 Fluorometer; Thermo Fisher Scientific, Waltham, MA, USA) and confirmed by gel electrophoresis on a 1% (w/v) agarose gel. Genomic DNA at a concentration of at least 100 ng/µL was used to prepare the libraries for each genotype.
The library construction for GBS was conducted according to a previous report. 26 In brief, genomic DNA was digested with the restriction enzyme ApeKI,27,28 followed by ligation with a barcode adaptor and a standard Illumina sequencing adaptor. DNA fragments were pooled for polymerase chain reaction (PCR) amplification. Finally, 100 bp fragments were single end-sequenced on an Illumina HiSeq 2000 platform.
The high-quality FASTQ read sequences generated for each accession were aligned to the reference TM-1 cotton genome. 4 We applied Samtools 29 to produce BAM files for removing unmapped reads based on the mapping outputs. Vcflib packages (https://github.com/vcflib/vcflib.git) were then used to filter SNPs with a mapping quality score <30.
Population structure analysis
Population structures were determined in 2 steps. First, we applied principal coordinate analysis (PcoA) to investigate genetic relationships using a dissimilarity matrix obtained by DARwin 6.0 (http://darwin.cirad.fr). The PcoA results were plotted using the ggplot2 package in R studio. We also applied the discriminant analysis of the principal components (DAPC) using adegenet package, 30 which can determine relationships for redefined groups without requiring an a priori population genetics model. 31 In brief, the data were first transformed using PCA, and then the number of genetic clusters was assessed using the find clusters function. The Bayesian information criterion (BIC) was calculated for K = 1 to 10. For K-means clustering, all of the principal components were retained, and the K value with the lowest BIC was selected as the optimal number of clusters. DAPC was implemented using the optimized number of principal components as determined by the optim.a.score function. We further used the fast STRUCTURE tool to determine the most probable number of genetic clusters, which was run at K = 1 and K = 10 with default parameters. 32 Finally, we conducted the TreeMix (http://treemix.googlecode.com) to estimate population differentiation among all G hirsutum races and G hirsutum cultivar group by constructing the Maximum likelihood tree with m value (0-6) and block 1000, setting the G barbadense as the outgroup
Evidence of selection footprints in G hirsutum races
Population structure analysis prompted us to perform population genomic FST scans between G hirsutum cultivars and G hirsutum race groups to identify SNP-specific high FST outliers using both BAYESCAN version 2.1 33 and Arlequin v3.5. 34 For BAYESCAN, the “snp” option was applied using SNP genotype matrix as input data. The analyses were run using default settings, including 20 pilot runs of 5000 steps each, followed by 50 000 burn-in and 5000 sampling steps with a thinning interval of 10. The prior odds parameters were set to the default of 10. The false discovery rate (FDR) was set to 0.1 with the PLOT_BAYESCAN R function for outlier detection. For Arlequin, 50 000 simulations were run on the same data set with default parameters, using both the “neutral mean FST” and “force mean FST” options. Loci outside the 95% confidence interval and those with FST = 1 were considered outliers. For Arlequin, 20 000 simulations were run with 10 simulated groups and 100 demes per group to identify candidate loci under selection. High FST outlier SNPs were considered candidates for evidence of positive selection under population divergence. We identified all genes containing outlier SNPs (outlier genes) based on the TM-1 reference genome annotation and analyzed their functions based on known functions of Arabidopsis thaliana orthologous genes.
The phenotype analysis (Table 1) indicated significant differences with respect to seed germination and lint traits between G hirsutum cultivars and most of the G hirsutum race accessions tested. Therefore, we further tested for genome-wide signals of marker-phenotype associations using these differing phenotypes (seed germination and lint) to determine the selective footprint during the evolution of G hirsutum races to hirsutum cultivars. A mixed linear model (MLM) was used to analyze marker-trait associations with TASSEL 5.0, 35 which first discards heterozygous sites and then generates different sets of marker data. SNP data were further filtered using filter data with 0.05 site minor allele frequency (MAF) thresholds before association mapping. To correct the population structures, kinship analysis and PCA were carried out to obtain K and Q matrices, respectively, which were then used as a covariance matrix and integrated into the MLM. Correction for multiple tests was performed based on an FDR of 0.05 to identify significantly associated markers.36,37 The sequences of significant markers within genes were then used as queries for BLAST searches in the National Center of Biotechnology Information gene database based on the TM-1 genome sequence. Known genes linked to the significant loci were assigned as putative candidates based on the functions of A thaliana orthologous genes.
Phenotypic variations between Gossypium hirsutum races and cultivar populations.
Results
Phenotypic characterizations of G hirsutum races and G hirsutum cultivars
The biodiversity analysis demonstrated significant differences in several traits including seedling weight, embryonic axis length, bolls per plant, and single boll weight between the 7 races groups and cultivars (Table 1; Supplemental Table S2). Cultivars had a higher germination rate and lint percentage compared with races groups. However, no significant differences were detected for fiber quality traits (fiber length, uniformity index, and micronaire value) and fruiting branch traits. Only the Palmeri race genotypes were significantly different for the seed index from cultivars.
In the PCA, the first 2 components accounted for approximately 48.54% of the variation observed between G hirsutum races and G hirsutum cultivars (Figure 1; Supplemental Tables S3), and all accessions and cultivars could be simply clustered into 5 groups, which were the G barbadense cultivars group, G hirsutum cultivars group, Latifolium, Richmondi, and Marie-Galante group, Richmondi, Morrilliand, Punctatum, and Palmeri group along with the first component (32.26%). The second component (16.38%) separated Richmondi, Morrilliand, Punctatum, and Palmeri accessions group from the cultivar group. The PC analysis scatter plot suggests that most of these G hirsutum races occupy the most distant part of the ellipse including G hirsutum cultivars (Figure 1). The dendrogram of the AHC analysis further supported the PC analysis by classifying all samples into 3 major groups when the dissimilarity is at 200 (Figure 2), in which Punctatum, Latifolium, Morrilli, and Marie-Galante accounted for the more extended part of morphological variation along with G hirsutum. Furthermore, Latifolium, Marie-Galante, Morrilli, and Punctatum accessions were closer to G hirsutum cultivars group; meanwhile, these groups were clearly separated from the AD1 and AD2 samples that formed singular groups (Supplemental Table S4).
Scatter plot distribution of the first and second principal components (PCs) based on 15 morphological characteristics of 94 Gossypium hirsutum races and 21 cotton cultivars. The different color circles represent different clusters. The orange circle mainly contained Gossypium barbadense cultivars, the green circle mainly contained G hirsutum cultivars, the purple circle mainly contained Latifolium and Marie-Galante, the blue circle mainly contained Richmondi, Morrilli, and punctatum, and the light brown circle mainly contained Palmeri.
Agglomerative hierarchical clustering of 15 morphological characteristics of 94 Gossypium hirsutum races and 21 cultivated cotton accessions using the Ward’s agglomeration method and Euclidean distance dissimilarities; the horizontal dotted line represents the dissimilarity value is at 200; all cotton species can be divided into 3 groups.
Genome-wide detection of SNPs using GBS
A total of 877.41 million (98.4% of total reads) high-quality sequence reads (containing enzyme sites) were obtained. Sequence reads varied between 2.36 and 17.46 with a mean of 7.63 million reads per accession (Supplemental Figure S1, Supplemental Table S5). Approximately, 79.87% (75.7%-89.9%) of the sequenced nucleotides were evenly distributed across the 94 G hirsutum races and 21 cotton cultivars. The uniquely mapped sequence reads from the accessions or cultivars showed coverage from 8.27% (minimum) to 63.92% (maximum) of the G hirsutum acc.TM-1 reference genome (2189.14M) and the unique mapping ratios which ranged from 75.76% to 89.08% are presented in Supplemental Table S5.
A total of 146 558 SNPs was identified using Stacks tool (Supplemental Table S6). The MAF values varied from 0.005 to 0.5 with an average of 0.15 (Supplemental Table S6). According to the reference genome sequence, the detected SNPs were physically mapped to 26 chromosomes with an average density of 64 SNPs per Mb (Supplemental Table S6, Supplemental Figure S2). The leveraged SNP density is highest on chromosome 2 (117.25 Mb) with 8419 SNPs and lowest on chromosome 7 (23.27 Mb) with 1754 SNPs (Supplemental Figure S2). Tajima’s D, Theta, and Pi were calculated for the filtered SNPs with a mean of −0.22 (ranging from −0.97 to 0.28), 0.22 (ranging from −0.97 to 0.28), and 0.20 (ranging from 0.16 to 0.23), respectively (Supplemental Table S7). The transition/transversion ratio was calculated as 1.89. Most of the identified SNPs (62.9%) were transitions (A/G or T/C), whereas transversion events (A/C, A/T, C/G, or G/T) accounted for 37.1% of all SNPs. We also determined the physical locations of 146 558 SNPs based on the reference genome annotations, 35 499 SNPs were localized to 11 935 genes (Supplemental Table S8), and 111 316 (44.3%) SNPs were localized in the intergenic regions. Overall, 29 028 (11.6%) SNPs mapped to exons (coding sequences), 23 623 (9.4%) SNPs mapped to introns, and 42 261 (16.8%) mapped in the downstream regulatory regions (3′ untranslated region [UTR]). The SNPs detected in the upstream regulatory regions (promoter and 5′UTR) accounted for 14.4% (36 253) of all the SNPs (Supplemental Table S8).
Phylogenetic relationships of the cultivated and the wild relatives of G hirsutum
Initial assessment of the phylogenetic relationships was conducted using individual-based PC analysis with the identified 146 558 SNPs. PC1 analysis divided the selected cotton species into 2 groups associated with the AD1 and AD2 genomes (Figure 3). According to PC2 analysis, there were also 2 main groups: the first group mainly comprised Punctatum races and another group contains the cultivars and 6 G hirsutum races. However, the AD1 cultivars slightly overlapped with those 6 G hirsutum races except Punctatum (Figure 3; Supplemental Figure S3, Supplemental Table S9). Nine distinct clusters were identified as a result of the discriminant analysis of principal component (DAPC) analysis (Figure 4A). Visualization of DAPC results clearly clustered the first 50 principal components. AD2 samples were still in a separate cluster, whereas G hirsutum races were more diverse and could not be clearly separated, with some accessions (Latifolium, Richmondi and Marie-Galante) appearing to be more closely related to the AD1 cultivars (Figure 4B; Supplemental Table S10). The population splits and migration events from TreeMix analysis are shown in Figure 4C and D. In the model, the sampled populations in the selected cotton species were seemed to be related to their common ancestor through a graph of ancestral populations (Figure 4D). Using genome-wide allele frequency data and a Gaussian approximation to genetic drift, the historical distance among the population was as follows: Latifolium >Richmondi > Marie-Galante > Morrilli > Palmeri > Yucatanense > Punctatum (Figure 4D). The data also suggested that G hirsutum cultivars have the closest relationship with Latifolium, Richmondi, and Marie-Galante races and have the most distant relationship with Punctatum and Yucatanense races (Figure 4D; Supplemental Figure S3). A similar classification pattern was observed from the fast STRUCTURE analyses, 9 genetic clusters were visible, and AD1 cultivars showed a common genomic background with that of G hirsutum races (Figure 5; Supplemental Table S11).
Principal components analysis of 115 cotton samples, including 94 Gossypium hirsutum race and 21 cultivar accessions based on 146 588 SNPs. The different color circles represent different clusters. The purple circle only contained Gossypium barbadense cultivars, and the green circle mainly contained Punctatum.
Discriminant analyses of principal components (DAPC). (A) The optimal number of clusters (K) as determined by “K-means” clustering. The graph shows an apparent decrease of the Bayesian information criterion (BIC) until K = 9, red dot, which is the most likely value of K. After this value, BIC increases. (B) Scatter plot based on the DAPC output for 4 assigned genetic clusters indicated by different colors. Dots represent different individuals. pop1 represents group 1 containing almost all races accessions; pop2 represents group 2 containing mainly Latifolium accessions; pop3 represents group 3 containing mainly Punctatum accessions; pop4 represents group 4 containing mainly Punctatum and Marie-Galante accessions; pop5 represents group 5 containing mainly Latifolium, Richmondi, and all G hirsutum cultivars; pop 6 represents group 6 containing 2 Marie-Galante, 2 Palmeri accessions, Morrilli, and Marie-Galante accessions; pop 7 represents group 7 containing 2 Yucatanense accessions; pop 8 represents group 8 containing Richmondi, Morrilli, and Marie-Galante accessions; pop9 represents group 9 which only contained G barbadense cumainly cultivars (Table S10). (C, D) Plotted is the structure of the graph inferred by TreeMix for cotton populations, allowing 10 migration events. The scale bar shows 10 times the average standard error of the entries in the sample covariance matrix.
Phylogenetic structures of 115 cotton species estimated using the fast structure admixture model at K = 9.
Footprints of positive selection during G hirsutum domestication
In FST outlier scans, Arlequin yielded a significant number of high outlier SNPs than did BAYESCAN (Table 2). Cultivars of G hirsutum and Latifolium race pair produced fewer outliers than other pairs, whereas G hirsutum cultivars and Punctatum race pair had the highest number of outliers. In general, 54 outlier SNPs were located in the coding sequence, consistent with the proportion (~88%) among all outlier SNPs tested (Table 2), with 2 of the same SNPs being located in each of the genes LOC107896563, LOC107963906, LOC107913656, and LOC107927053. Thus, a total of 50 genes with outlier SNPs (designated as outlier genes) were considered as possible footprints of positive selection during G hirsutum domestication. This conclusion was reached without the need to analyze the outliers from each group of G hirsutum races separately. Gene Ontology terms (biological process) and the established functions of A thaliana orthologous genes are presented in Supplemental Table S12 for the outlier SNPs. Several outlier genes are predicted to be associated with pollen germination and tube growth and 3 genes (LOC107896563, LOC107927053, and LOC107913656) with the regulation of flower development. Other processes shared by more than 1 pair of comparisons include hormone pathways, biotic and abiotic stress responses. Intriguingly, several candidate genes were closely clustered in a specific region on the chromosome. This suggests that some outlier SNPs reside in areas that may have undergone sweeps of selection, which would make it more difficult to identify the specific genes targeted by selection.
Summary of high FST SNP outliers from BAYESCAN and Arlequin analyses using 16 288 SNPs.
Abbreviations: SNP, single nucleotide polymorphism; G hirsutum, Gossypium hirsutum.
Footprints of positive selection were investigated based on the genome-wide-association phenotype of early seedling development (stress responses) using 92 G hirsutum races, which contained 98 436 SNPs. Although no significant signal could be detected throughout the whole genome or on a specific chromosome, this trait did correlate with an SNPs-per-chromosome with the largest −log10(3e−06) value (14 SNP-containing genes) being found on G hirsutum cultivars (Figure 6; Supplemental Table S13). Three of these genes (LOC107914109, LOC107922201, and LOC107921406 are predicted to be involved in the biological processes including cellular response to water deprivation, defense responses to bacterial attack, reductive pentose-phosphate cycle, response to cold, and response to nematode infection (Table 3).
Genome-wide association study of early seedling development rate-related traits based on the genotyping-by-sequencing single nucleotide polymorphisms (SNPs), an SNPs-per-chromosome with the more than −log10(3e−06) value was selected to be involved in stress resistance–related traits. The red line means the cut off standard with −log10(3e−06); the number in X-axis represents chromosome number of Gossypium hirsutum.
Selective footprint of stress resistance–related traits between 94 Gossypium hirsutum races and 21 cultivated cotton populations.
Discussion
The GBS assay is a robust, simple, and affordable tool for SNP discovery and genome mapping. By using appropriate restriction enzymes and PCR amplifications, it can sufficiently reduce genome complexity, avoid repetitive regions of genomes, and target lower copy regions. In addition, it has the advantages of minimizing alignments problems in genetically highly diverse species and dealing with a large number of genome samples at a lower cost than other methods. Theoretically, the method can be used to map any plant species, as it does not require a reference genome. It has been frequently used to analyze large segregating progenies and marker trait association studies based on linkage disequilibrium and even the evolutionally genomic selection. 38 When the method was originally developed, it was used to analyze a high-resolution maize mapping population and doubled haploid barley lines for GBS accuracy and efficiency. The results demonstrated that 25 185 biallelic tags could be mapped to the maize genome and 24 186 sequence tags to the barley genome and the GBS reads were in 99% agreement with the reference markers. As a consequence, it has become a practical and convenient approach for genotyping in crop plants,39-43 including the analysis of cotton for genetic diversity, genetic mapping and GWAS analysis.44-47
In the current study, a total of 146 558 high-quality SNPs were identified in G hirsutum races and cultivated cotton accessions, which is consistent with the genome-wide SNP discovery capability previously reported in other crop plants.48-52 Moreover, 35 499 SNPs were located in 11 935 cotton genes correlating with some important agricultural traits such as fertility, germination, and fiber quality; thus, they can serve as useful marker-based resource or molecular design tools for genomic selection, genome-wide association analyses, and genetic diversity of G hirsutum. The obtained SNPs also have the potential for breeding programs aimed at improving cotton quality and yield.
Our results confirm those of others who have proposed that G hirsutum and other races were derived from a Yucatanense progenitor in the Old World. 11 Similarly, the Punctatum, Palmeri, Marie-Galante, and, especially, Latifolium races are considered to be progenitors of the modern G hirsutum cultivars in the New World.53-56 Our genotyping results are also consistent with the classification based on phenotypic similarities, that Latifolium, Punctatum, and Marie-Galante should be included in the G hirsutum cultivars group. 57 Overall, the results are consistent with the established history of cotton domestication and support the suggestion that G hirsutum cultivars originated from 2 distinct places in the Old World and New World from our TreeMix analysis.4,54,58
Genome scans and FST outlier SNP discovery approaches are effective strategies for identifying genes under evolutionary or domestication selection, and this approach is significantly different from phenotype-based GWAS analysis. 59 Extreme differentiation in allele frequencies between genetic groups or populations in different geographic zones, as measured by the FST statistic, provides a signature of recent positive selection, 60 and different populations of ancestral accessions/species can further serve as a useful biological model to study adaptive or directional selection in nature. 61 Our study indicates that, during cotton evolution and domestication, a large portion of the outlier genes are involved in reproductive processes such as the regulation of pollen germination and tube growth, the regulation of flower development, hormone signaling processes, biotic and abiotic stress responses. Outlier genes were closely clustered in a specific region on the chromosome suggesting that the whole chromosome fragment, and not just specific genes, has undergone sweeps of selection. As a result, it is difficult to isolate or characterize the genes contributing to specific traits from this region. This finding revealed direct genetic evidence for a positive selection from G hirsutum races to G hirsutum cultivars. Indeed, these genes appear to be involved in the regulation of range of important agricultural traits; therefore, they may be able to serve as candidate molecular markers for G hirsutum cultivars breeding programs using G .hirsutum races.
The most strongly differentiated trait between G hirsutum races and G hirsutum cultivars was their resilience during the domestication. G hirsutum races underwent some phenotypic adaptations including exhibiting more sensitivity to photoperiod changes; loss of perennial growth habits; reduced seed dormancy; and greater resistance to various stress conditions.62,63 Three SNPs were found in genes corresponding to potential Arabidopsis orthologous involving in responses to stress. These results also suggest that genes with the functions contributing to abiotic and biotic stress responses are conserved in the G hirsutum cultivars as a consequence of artificial selection for improved environmental adaptations.
In our study, the outlier SNPs, as indicators of positive selection, did not associate closely with genome-wide marker-phenotype association signals. This may be due to the possibility that the GBS mapping approach may miss some sequences due to the requirement for genome cleavage by specific restriction enzymes. Alternatively, the outlier genes were detected using 146 588 SNPs, whereas the genome-wide marker association analysis was conducted using only 98 436 SNPs, and therefore, it might not reflect genomic coordination between the selection loci and SNP markers of the genome-wide association. Another explanation is that SNPs may not be able to adequately represent the major signature of selection on coding variants as comprehensively as gene expression analyses. Further studies using transcriptomics or higher SNP density or larger number of SNPs could help to improve the resolution of differentially selected loci and increase the concordance between the SNPs identified by the selection footprint analysis, transcription analysis, and genome-wide phenotype association studies.
Conclusions
The current study adopted a GBS approach to analyze 94 G hirsutum races accessions and 21 cotton domesticated cotton cultivars. We concluded that the Latifolium, Richmondi, and Marie-Galante accessions were more genetically related to the selectively domesticated G hirsutum cultivars, and 54 outlier SNPs were identified and 3 SNPs located in genes related to plant responding to stress were isolated based on their orthologous genes function. These findings provide a preliminary indication of adaptation and selection footprints during cotton domestication and offer candidate DNA markers that could be used for cotton breeding programs.
Supplemental Material
Supplemental_File_Tables0908_xyz27463d4de7434 – Supplemental material for Genotyping-by-Sequencing of Gossypium hirsutum Races and Cultivars Uncovers Novel Patterns of Genetic Relationships and Domestication Footprints
Supplemental material, Supplemental_File_Tables0908_xyz27463d4de7434 for Genotyping-by-Sequencing of Gossypium hirsutum Races and Cultivars Uncovers Novel Patterns of Genetic Relationships and Domestication Footprints by Shulin Zhang, Yaling Cai, Jinggong Guo, Kun Li, Renhai Peng, Fang Liu, Jeremy A Roberts, Yuchen Miao and Xuebin Zhang in Evolutionary Bioinformatics
Footnotes
Funding:
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the National Natural Science Foundation of China (31671745), Innovation Scientists and Technicians Troop Construction Projects of Henan Province (184200510009), and The National Key Research and Development Program of China (2016YFD0100203-1-4, 2018YFD0100300).
Declaration of conflicting interests:
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Author Contributions
Y.M., X.Z., and R.P. contributed to the research design; S.Z. and J.G. conducted experiments and investigated the study; and S.Z. and X.Z. wrote the manuscript.
Data Availability
All the genotyping-by-sequencing (GBS) data during the current study are available in the NCBI Sequence Read Archive (SRA) under project accession number PRJNA498359. The authors state that all data necessary for confirming the conclusions stated are represented fully within the article and in Supplemental Materials.
Supplemental material
Supplemental material for this article is available online.
References
Supplementary Material
Please find the following supplemental material available below.
For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.
For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.
