ORIGINAL RESEARCH article

Front. Plant Sci., 15 April 2022

Sec. Plant Systematics and Evolution

Volume 13 - 2022 | https://doi.org/10.3389/fpls.2022.873788

Nucleotide Evolution, Domestication Selection, and Genetic Relationships of Chloroplast Genomes in the Economically Important Crop Genus Gossypium

  • 1. Shaanxi Key Laboratory for Animal Conservation, Key Laboratory of Resource Biology and Biotechnology in Western China (Ministry of Education), College of Life Sciences, Northwest University, Xi’an, China

  • 2. State Key Laboratory of Cotton Biology, Institute of Cotton Research, Chinese Academy of Agricultural Sciences, Anyang, China

  • 3. Gansu Provincial Key Laboratory of Aridland Crop Science, College of Life Science and Technology, Gansu Agricultural University, Lanzhou, China

Abstract

Gossypium hirsutum (upland cotton) is one of the most economically important crops worldwide, which has experienced the long terms of evolution and domestication process from wild species to cultivated accessions. However, nucleotide evolution, domestication selection, and the genetic relationship of cotton species remain largely to be studied. In this study, we used chloroplast genome sequences to determine the evolutionary rate, domestication selection, and genetic relationships of 72 cotton genotypes (36 cultivated cotton accessions, seven semi-wild races of G. hirsutum, and 29 wild species). Evolutionary analysis showed that the cultivated tetraploid cotton genotypes clustered into a single clade, which also formed a larger lineage with the semi-wild races. Substitution rate analysis demonstrated that the rates of nucleotide substitution and indel variation were higher for the wild species than the semi-wild and cultivated tetraploid lineages. Selection pressure analysis showed that the wild species might have experienced greater selection pressure, whereas the cultivated cotton genotypes underwent artificial and domestication selection. Population clustering analysis indicated that the cultivated cotton accessions and semi-wild races have existed the obviously genetic differentiation. The nucleotide diversity was higher in the semi-wild races compared with the cultivated genotypes. In addition, genetic introgression and gene flow occurred between the cultivated tetraploid cotton and semi-wild genotypes, but mainly via historical rather than contemporary gene flow. These results provide novel molecular mechanisms insights into the evolution and domestication of economically important crop cotton species.

Introduction

Since Darwin’s time, biologists have recognized that investigating the human domestication of wild plants can help to improve our understanding of the evolutionary process (Yoo et al., 2014). Generally, domesticated forms of cultivated species differ from their wild counterparts in numerous traits (; ). Insights into the evolution of chloroplast genome’s domestication and selection are made possible by comparative studies of wild and domesticated representatives of individual cultivated species. In the previous study, scholars used chloroplast genome data to analyze the genetic variation and evolution of olive. As a control, the cultivated species were employed to analyze genome variation and genetic association among olive chloroplasts (). Meanwhile, some other study have also examined the evolutionary mechanism of the chloroplast genome of cultivated Camellia sinensis and its relatives (). In recent studies, comparisons of wild and domesticated plants have provided important insights into the developmental mechanisms that underlie traits affected strongly due to targeted selection by humans (Yoo et al., 2014). In general, domesticated plants are characterized by reduced genetic variation and relaxed selection pressure compared with their wild counterparts. Several studies also found high levels of continuous gene flow from wild to cultivated genotypes (; ; ; ). Thus, the domestication process may provide a basis for studying the overall evolutionary relationships associated with wild crop transformation and identifying the genes under selection (; ).

Cotton (Gossypium) is one of the most important crops worldwide (Wendel, 1989; ) and a major source of natural fiber for the textile industry. Allopolyploid cotton originated in the New World and diverged into at least six species throughout the tropical and subtropical Americas: G. hirsutum (AD1), G. barbadense (AD2), G. tomentosum Nuttalex Seemann (AD3), G. mustelinum Miersex Watt (AD4), G. darwinii Watt (AD5), and G. ekmanianum (AD6) (Wendel and Cronn, 2003; Wendel and Grover, 2015). The diploid species comprise eight monophyletic genome groups: A, B, C, D, E, F, G, and K (Wendel and Cronn, 2003; ; Wendel et al., 2010). These groups can be separated into three main lineages in three continental regions: 13 D-genome species from the American continents, 15 species from the Asian and African continents (A-, B-, E-, and F-genomes), and 18 species (C-, G- and K-genomes) from Australia (Wendel and Cronn, 2003). Hence, cotton species provide a fascinating model system for studying evolution, domestication selection, genetic introgression, and gene flow among different continents (, ; Wendel, 1989; Wendel and Grover, 2015; , ,). Four species in the genus Gossypium are cultivated for the production of spinnable fiber, i.e., two allotetraploid species comprising G. hirsutum L. and G. barbadense L. (2n = 4x = 52), and two diploid species comprising Gossypium herbaceum L. (A1) and Gossypium arboretum L. (A2) (2n = 2x = 26) (Wendel and Cronn, 2003; Wendel and Grover, 2015). Allopolyploid cottons were considered to be about 1.5 million years old and were domesticated by humans 4,000 to 5,000 years ago (Wendel, 1989; Wang et al., 2017), which were originally domesticated from tree cotton in the Mesoamerican and Caribbean regions, and then further domesticated and improved in the southern United States (). And two diploid cotton species, G. arboretum and G. herbaceum, have been cultivated for several millennia (Simon et al., 2016), which were initially domesticated on Madagascar or in the Indus Valley (Mohenjo Daro), and was subsequently dispersed to Africa and other areas of Asia (Wendel and Grover, 2015; ; ). Due to the high-yield characteristics of allopolyploid cottons, the American upland cottons have been introduced and replaced by two diploid cotton species (G. arboretum and G. herbaceum) (; ). Up to now, the Upland cotton (G. hirsutum) accounts for more than 95% of the worldwide production of cotton (Yoo et al., 2014; ; ; Wang et al., 2019; Zhang et al., 2020).

Following human-mediated selection and agronomic improvement, the ability of cotton species to adapt to various environments was enhanced and the production of fiber from cotton improved significantly (). The domestication process also resulted in other morphological changes in other crops such as sorghum, rice and soybean (), including early flowering, larger and/or more fruits, annualized habit, plant height reduction, and loss of seed dormancy (Yoo et al., 2014). When plants undergo artificial domestication, the relaxation of certain features is inevitable (), that is, when plants undergo relatively large changes, such as from the transition from nature to domestication, certain characteristics important for survival in nature lose much of their adaptive significance under artificial directional selection. Hence, one would expect natural selection for such characteristics to lose its intensity (; ). Many studies have shown that the genetic diversity of upland cotton varieties is low, mainly due to several bottlenecks in the domestication process (; ; ; Wendel and Cronn, 2003). In addition, previous studies based on whole-genome resequencing of upland cotton have indicated that the genomic diversity of upland cotton decreased under the stress of artificial selection (; ). Thus, in the current era of genomic big data, high-throughput “omics” sequencing techniques allow detailed analyses of the genetic changes associated with artificial domestication, as well as providing new, accurate, and targeted genome-based crop breeding strategies (Wang et al., 2017; ; Yang et al., 2020). For example, in maize and rice, the use of high-quality backbone parents can obtain notable improvements in breeding efficiency (). The whole-genome sequences of allotetraploid cotton and its ancestors have been completed, and the high-quality allotetraploid upland cotton genome is an effective tool for systematically exploring the genomic mysteries of polyploidy (, ; Zhang et al., 2015). Compared with whole-genome sequencing, the chloroplast genome is single-copy, maternally inherited, and there is no chain exchange or free combination phenomenon. It has a relatively independent evolutionary route. In addition, the highly conserved characteristics of the chloroplast genome make them useful for the rapid analysis of species evolution (; ; Wang et al., 2013; ). However, the whole-genome resequencing (WGR) is parental inheritance, and there may be genetic recombination (; Wang et al., 2020).

In the current study, to better understand the evolution, domestication selection, and genetic relationships of cotton, we analyzed the chloroplast genomic variation in 72 cotton genotypes comprising G. hirsutum and its 29 cultivated upland cotton accessions, G. barbadense and its three cultivated accessions (Gossypium barbadense cultivar zhonghai 7, Gossypium barbadense cultivar Kaiyuan, and Gossypium barbadense cultivar yuanmou), G. africanum, G. arboretum, seven semi-wild races of G. hirsutum, and 29 wild cotton species. We also estimated molecular dating, genetic introgression, nucleotide substitutions, and indel variation.

Materials and Methods

DNA Extraction and Plant Materials

The fresh leaves of seven semi-wild races of upland cotton, i.e., punctatum, latifolium, richmondi, morrilli, marie-galante, palmeri, and yucatanense, were collected from the National Wild Cotton Nursery in Sanya, China. In addition, 29 cultivated upland cotton accessions were also obtained from different ecological geographic regions, with three accessions from the United States, eight from the Yellow River region, 12 from the Yangtze River area, four from northwest China, and two from north China (Table 1). Leaf tissues were dried with silica gel and genomic DNA was extracted using the modified CTAB method (). Approximately 5 μg of purified DNA was used to construct paired-end libraries with an insert size of 350 bp and sequencing was performed with the Illumina HiSeq 2500 platform by Novogene (Beijing, China). Additionally, we have also downloaded the 36 chloroplast genomes of cotton species from NCBI (National Center for Biotechnology Information) for further combination analysis.

TABLE 1

NumberSpeciesAccession numberSourceLogogram
1Gossypium punctatumMK792868Sanya, ChinaJBM
2Gossypium richmondiiMK792869Sanya, Chinalqmd
3Gossypium morrilliMK792866Sanya, ChinaMLE
4Gossypium marie-galanteMK792865Sanya, ChinaMLJ
5Gossypium palmeriiMK792867Sanya, ChinaPME
6Gossypium yucatanenseMK792870Sanya, ChinaYKT1
7Gossypium hirsutum cultivar 06G415MK792871Yellow riverS32
8Gossypium hirsutum cultivar antongSP21MK792837United StatesS24
9Gossypium hirsutum cultivar chuanmian45MK792838Yangtze riverS47
10Gossypium hirsutum cultivar CJL-233MK792839Yangtze riverS252
11Gossypium hirsutum cultivar difenmian168MK792840Yangtze riverS64
12Gossypium hirsutum cultivar ekangmian7MK792841Yangtze riverS273
13Gossypium hirsutum cultivar emian12(4947)MK792842Yangtze riverS263
14Gossypium hirsutum cultivar gaochanbukangchong RRMMK792843Yangtze riverS246
15Gossypium hirsutum cultivar guangyedaizimianMK792844United StatesS59
16Gossypium hirsutum cultivar guokang12 (GK12)MK792845Yellow riverS156
17Gossypium hirsutum cultivar hanmian802MK792846Yellow riverS162
18Gossypium hirsutum cultivar humian204MK792847Yangtze riverS257
19Gossypium hirsutum cultivar Jan-86MK792848Yellow riverS211
20Gossypium hirsutum cultivar liaomian10MK792849North ChinaS234
21Gossypium hirsutum cultivar lumianyan21(lu1138)MK792850Yellow riverS163
22Gossypium hirsutum cultivar shan401MK792851Yellow riverS10
23Gossypium hirsutum cultivar simian4MK,792852Yangtze riverS272
24Gossypium hirsutum cultivar sizimian4MK792853United StatesS38
25Gossypium hirsutum cultivar sumian5MK792854Yangtze riverS45
26Gossypium hirsutum cultivar xinluzhong7MK792855Northwest ChinaS275
27Gossypium hirsutum cultivar xinluzhong9 (1318136-160)MK792856Northwest ChinaS277
28Gossypium hirsutum cultivar xinluzhong10MK792857Northwest ChinaS278
29Gossypium hirsutum cultivar xinluzhong19MK792858Northwest ChinaS281
30Gossypium hirsutum cultivar xuzhou209MK792859Yangtze riverS13
31Gossypium hirsutum cultivar yanmian48MK792860Yangtze riverS265
32Gossypium hirsutum cultivar youLU272⊕MK792861Yellow riverS175
33Gossypium hirsutum cultivar yumian1MK792862Yangtze riverS271
34Gossypium hirsutum cultivar zhong053MK792863Yangtze riverS8
35Gossypium hirsutum cultivar zhongzhimian GD89MK792864Yellow riverS185
36Gossypium barbadense cultivar zhonghai7HQ901199NCBIAD2_ 99
37Gossypium barbadense cultivar kaiyuanHQ901200NCBIAD2_200
38Gossypium barbadense cultivar yuanmouHQ901198NCBIAD2_98
39Gossypium darwiniiNC_016670NCBIAD5_70
40Gossypium tomentosumNC_016690NCBIAD3_90
41Gossypium mustelinumNC_016711NCBIAD4
42Gossypium hirsutumNC_007944NCBIAD1_44
43Gossypium barbadenseNC_008641NCBIAD2_41
44Gossypium africanumNC_016692NCBIA1_a
45Gossypium arboreumNC_016712NCBIA2
46Gossypium longicalyxJF317354NCBIF1
47Gossypium anomalumJF317356NCBIB1
48Gossypium capitis-viridisNC_018111NCBIB3
49Gossypium sturtianumJF317353NCBIC1
50Gossypium nandewarenseMG779276Sanya, Hainan, ChinaC1–n
51Gossypium robinsoniiNC_018113NCBIC2
52Gossypium bickiiJF317352NCBIG1
53Gossypium australeNC_033401NCBIG2
54Gossypium populifoliumNC_033398NCBIK2
55Gossypium thurberiJF317353NCBID1
56Gossypium armourianumMG891801Sanya, Hainan, ChinaD2–1
57Gossypium harknessiiNC_033333NCBID2–2
58Gossypium klotzschianumNC_033394NCBID3–k
59Gossypium davidsoniiNC_033395NCBID3–d
60Gossypium aridumNC_033396NCBID4
61Gossypium raimondiiNC_016668NCBID5
62Gossypium gossypioidesNC_017894NCBID6
63Gossypium lobatumMG891802Sanya, Hainan, ChinaD7
64Gossypium trilobumMG800783Sanya, Hainan, ChinaD8
65Gossypium laxumKF806549NCBID9
66Gossypium turneriNC_026835NCBID10
67Gossypium schwendimaniiMG891803Sanya, Hainan, ChinaD11
68Gossypium stooksiiJF317354NCBIE1
69Gossypium somalenseNC_018110NCBIE2
70Gossypium areyiabumNC_018112NCBIE3
71Gossypium incanumNC_018109NCBIE4
72Gossypium latifoliumMG800784Sanya, Hainan, Chinakym

List of taxa sampled in this study and species accession numbers (GenBank).

Chloroplast Genome Assembly, and Annotation

The raw sequencing reads obtained by the company (Novogene, Beijing, China) were filtered through the “AmbiguityFiltering.pl” script in the NGSQCToolkit software (), and removed the fragments with fuzzy bases greater than 2% and those with bases less than 50 bp. The clean reads were assembled by the MIRA 4.0.2 program () where the complete chloroplast genome of G. hirsutum (AD1) (NC_007944) was used as the reference sequence in this process. In order to further assemble the whole chloroplast genomes, some ambiguous regions were extended using the MITObim v1.7 program with a baiting and iteration method (). The contigs obtained were used to generate consensus sequences with Geneious v8.0.2 (). The chloroplast genomes were then annotated using the Dual Organellar Genome Annotator (DOGMA, Wyman et al., 2004) program and manual corrections were made for some specific genes. All tRNA genes were further confirmed using the online tool tRNAscan-SE (Schattner et al., 2005). All of the newly generated genome sequences were submitted to GenBank (accession numbers MK792837MK792871 and MG800784).

Genetic Clustering Analysis

To evaluate the genetic relationships among cotton genotypes, molecular phylogenetic analysis was conducted using 72 complete chloroplast genome sequences (Table 1) and two outgroups comprising Bombax ceiba (NC_037494) and Theobroma cacao (NC_014676). First, all of the sequences were aligned using the MAFFT program () and the best-fit model was then selected with Modeltest v3.7 () based on Akaike’s information criterion. Finally, a maximum likelihood tree was constructed using RAxML v7.2.8 (Stamatakis, 2006) where the best model was GTR + G based on 1000 bootstrap replicate tests.

Estimation of Divergence Times

Previously estimated dates of speciation events (fossil records) were used to calibrate the phylogenetic tree (). In BEAST v1.8.0 (), we used the Yule process speciation prior and the uncorrelated lognormal model of rate change with a relaxed clock to estimate the divergence times among cotton lineages. The divergence time was calculated based on 74 chloroplast protein-coding sequences shared by the cotton genotypes, and we used three fossil records: AD1 (G. hirsutum) and A2 (G. arboreum) diverged 1–2 Mya (Wendel, 1989), A2 (G. arboreum) and D5 (G. raimondii) diverged ∼ 5–10 Mya (Senchina et al., 2003), and Theobroma-Gossypium diverged 60 Mya (). A normal prior probability distribution was used to account for the uncertainty of prior knowledge. The analyses were run for 50,000,000 generations and the parameters were sampled every 5,000 generations. Tracer v 1.6 () was used to determine the effective sample size (>200) and the first 20% of the samples were discarded as burn-in. Tree Annotator v.1.8.0 () was used to summarize the set of post-burn-in trees and their parameters were used to produce a maximum clade credibility chronogram, which illustrated the mean divergence time estimates in the 95% highest posterior density (HPD) intervals. Finally, FigTree V1.3.1 () was used to visualize the molecular dating estimates.

Analysis of Nucleotide Substitutions

Transitions/transversions explain the substitution rates of nucleotides, so we determined the transition/transversion rates using single nucleotide polymorphism (SNP) loci in protein-coding sequences in the cotton chloroplast genome. These analyses were conducted based on two genetic groups obtained from the phylogenetic analyses. One group contained the diploid cotton species (including G. africanum and G. arboretum) and the other group comprised tetraploid semi-wild races and cultivated upland cotton genotypes (excluding G. barbadense). MEGA files generated from SNP data were analyzed with MEGA7 software () to obtain the transition/transversion rate. The following parameters were employed: statistical method, maximum likelihood; analysis, substitution pattern estimation (MCL); substitution type, nucleotides; scope, all selected taxa; model/method, Tamura–Nei (automatic selection); gaps/missing data treatment, partial deletion, and site coverage cut off (%), 95 (). Finally, we converted the transition/transversion rates for the two groups into two histograms. In addition, DnaSP v5.10 () was used to calculate the non-synonymous (dN) and synonymous (dS) mutations in coding regions for the two groups.

Estimation of Mutation Rates

The two cotton groups described above were also used to calculate the mutation rates. The rate of mutation per site per year (μ) was estimated using the formula: μ = m/(nT), where m is the number of observed mutations, n is the number of total sites, and T is the divergence time of a node (). The μ values for structural mutations were calculated using the method described by Saiton and Ueda (), where the total number of structural mutations was divided by the additive time based on the branch lengths and by the length of the nucleotide sequences. Finally, we calculated the evolutionary rates for nucleotide substitutions and indels. The indel rates were calculated for the two groups using DnaSP v5.10 ().

Selection Pressure Analysis

To identify domestication selected genes, we performed selection pressure analysis using the Codeml program (Yang et al., 2005) and two different groups of genotypes, where one group comprised the wild diploid cotton species with a total of 28 genotypes and the other group contained the upland cotton semi-wild races and cultivated varieties with a total of 37 cotton genotypes (excluding G. barbadense and its three cultivated accessions, i.e., G. tomentosum, G. mustelinum and G. darwinii, because these seven genotypes were not involved in the domestication selection process for upland cotton). In general, the non-synonymous (dN) and synonymous substitution (dS) rate ratio (ω = dN/dS) was sensitive to selection pressure during evolution at the protein level, and it was particularly useful for identifying positive selection. Geneious v8.0.2 () and MAFFT v7.0.0 () were used to extract and align 77 protein-coding chloroplast genes from the two groups. Maximum likelihood phylogenetic trees were constructed based on the complete chloroplast genome sequences using RAxML v7.2.8 (Stamatakis, 2006). This model allowed the ω ratio to vary among sites with a fixed ω ratio for the whole tree to test for site-specific evolution in the gene phylogeny (Yang and Nielsen, 2002). Log-likelihood values of every model were compared against a neutral model based on likelihood ratio tests in order to determine statistically significant differences. Only the candidate sites for positive selection with significant support based on the posterior probability (p of (ω > 1) ≥0.99; Bayes Empirical Bayes approach) identified by M2 and M8 were considered further.

Diversity and Genetic Structure Analysis

DnaSP v5.10 () was used to analyze the genetic diversity parameters based on the complete chloroplast genome sequences of seven semi-wild races and 29 cultivated upland cotton genotypes. We also calculated the haplotype diversity (Hd) (), nucleotide diversity (π) (), and the number of haplotypes (H) with DnaSP v5.10 software.

We also analyzed the genetic structure patterns using the Bayesian Markov chain Monte Carlo clustering analysis method implemented in STRUCTURE 2.3.3 (; ; ). The admixture model with correlated allele frequencies was implemented for each run without a prior placed on the population information (). We conducted eight independent runs for each value from K = 1–10 to estimate the “true” number of clusters in 200,000 Markov chain Monte Carlo cycles following a burn-in step of 500,000 iterations. The most likely number of clusters was defined using log probabilities [Pr(X| K)] () and the △K method () via the online website STRUCTURE HARVESTER (). Next, CLUMPP 1.1.2 and the Greedy algorithm were used to align multiple runs of STRUCTURE for the same K value (). Finally, we applied DISTRUCT 1.1 () to graphically visualize the individual probabilities of cluster membership.

Gene Flow

We calculated the historical gene flow in semi-wild races and cultivated upland genotypes using Migrate-n (). First, we generated five independent Markov chain Monte Carlo cycles, each with 5,000,000 generations. We then sampled every 100 steps under a constant variation model and discarded the first 1,000,000 records as a burn-in and the other settings were at their default values. After checking for data convergence, we estimated the mode and 95% HPD (). In addition, we applied BAYESASS v3.0 to detect contemporary gene flow in the two groups (Wilson and Rannala, 2003). In these calculations, the three parameters comprising the migration rates (ΔM), allele frequencies (ΔA), and inbreeding coefficients (ΔF) were used as references to ensure that the optimal acceptance rates for the three parameters fell within the range of 20–60%. After continuous calculations, the correlation values for the genetic components were finally determined as 0.03, 0.16, and 0.14, respectively. We then conducted the analyses based on 57 iterations after a burn-in of 56 iterations and set 1,000 as the sampling frequency. Ten separate runs were performed to minimize the convergence problem (). The method proposed by Meirmans was used to obtain the results with the lowest deviance ().

Results

Evolutionary Relationships

The chloroplast genome sequences and concatenated protein-coding genes were used to reconstruct the maximum likelihood phylogenetic relationships for 72 Gossypium genotypes, and the cotton relationships generated from the data sets had the same topology, as shown in Figure 1. The six major genetic clades identified comprised the A + AD, F, E, D, B, and C + G + K genomic groups. Interestingly, all of the cultivated upland cotton genotypes clustered with the semi-wild race latifolium, which also formed a large evolutionary lineage with the other semi-wild races. The A-genome cotton species and G. barbadense genotypes also formed a single clade and they were closest to the upland cotton branch, whereas the 13 D-genome species formed a strong monophyletic lineage. The Australian species (C + G + K) clustered into a small branch, which clustered into a large branch with the B-genome species. Four species representing the E-genomic group also clustered into a large evolutionary branch. These results were in good agreement with the biogeographic distributions of cotton species from different continents.

FIGURE 1

Divergence Time Estimation

The molecular dating showed that the divergence time between the genus Gossypium and outgroups (B. ceiba and T. cacao) was about 58.15 Mya (95% HPD = 56.53–60.04 Mya), which are consistent with previous estimates (; Figure 2). The genus Gossypium originated about 11 Mya (95% HPD = 9.34–11.74 Mya) and most genomic groups in the genus diverged radially in a relatively narrow time range. Interestingly, the divergence time between the B-genome (African origin) and Australian clades (C + G + K) was estimated at 7.7 Mya (95% HPD = 6.3–9.8 Mya), which again supported the genetic relationship present in the B-genome, i.e., the B-genome branch and Australian branch were strongly grouped phylogenetically. The semi-wild races and cultivated upland cotton accessions diverged about 3.12 Mya and the ancestor of the D-genome originated at 5 Mya (95% HPD = 3.59–5.44 Mya). The divergence time of the allotetraploid AD clade was about 3.37 Mya (95% HPD = 2.44–4.93 Mya).

FIGURE 2

Nucleotide Substitutions

The ratios of transition/transversion were high among the semi-wild races and cultivated upland cotton genotypes (1.41), but low among the genotypes of the wild cotton species (1.16) (Table 2). There were significant differences in the proportions of two transition mutations and four transversion mutations between the two groups (Figure 3). Among the four transversion mutations, the proportion of A-C + T-G mutations was similar to that of C-G + G-C mutations in the groups. In addition, few A-T + T-A and C-A + G-T mutations were found in all combinations.

TABLE 2

Semi-wild and cultivated cotton accessions
Wild cotton species
From\ToATCGTs/TvATCGTs/Tv
A4.27274.265522.22411.41004.67926.465413.94431.1600
T3.388013.98848.56705.730715.67456.2325
C3.388014.01208.56705.730711.34416.2325
G8.78894.27274.265512.82164.67926.4654

Ratios of transitions and transversions for plastid protein-coding sequences in cotton accessions.

Rates of different transitional substitutions are shown in bold, whereas those of transversional substitutions are not shown in bold.

FIGURE 3

The 4,074 biallelic SNPs were subdivided into coding, intron, and intergenic spacer regions, and sorted into two groups comprising wild cotton species, and semi-wild and cultivated upland cotton genotypes (Table 3). In wild cotton group, there were 3,753 SNPs in total: 1,375 in coding regions, 264 in intron regions, and 2,693 in intergenic spacer regions. The percentages of SNP to the total lengths were 1.72, 1.22, 2.95, respectively, manifesting the intergenic spacer region sequences were more variable than the intron regions. In the coding regions, there were 1,027 non-synonymous mutations and 347 synonymous mutations, and the dN/dS was about 2.96. In the semi-wild and cultivated cotton genotypes, the sequences of the intergenic spacers and intron regions were more variable than the coding regions. The dN/dS ratio (3.5) was larger for this group than the wild cotton species (56 non-synonymous mutations and 16 synonymous mutations).

TABLE 3

Wild accessions
Semi-wild and cultivated accessions
Genome regionLength (bp)Value%Length (bp)Value%
Total substitutions163,4003,7532.3166,2373210.19
Coding region79,7041,3751.7279,968770.1
Non-synonymous/1,0271.29/560.07
Synonymous/3470.44/160.02
dN/dS/2.96//3.5/
Intron21,5812641.2221,292140.07
Intergenic spacer81,3812,6932.9577,5241300.17

Taxonomic and genomic distribution of biallelic single nucleotide polymorphic loci in wild, semi-wild, and cultivated cotton plastid genomes.

Estimation of Mutation Rate

The evolutionary rates were calculated based on the lengths of the genomes, number of substitutions, and times since divergence. In total, 1,375 substitutions were estimated in the wild species group and 77 in the semi-wild races and cultivated upland cotton group. The evolutionary rate of nucleotide substitutions was 1.2 × 10–9 per site per year in the wild species group compared with 0.18 × 10–9 per site per year in the semi-wild and cultivated group. In addition, 479 indels were identified in the wild cotton species and the evolutionary rate for indels was estimated at 0.4 × 10–9 per site per year. In the semi-wild and cultivated group, 24 indels were detected and the evolutionary rate was estimated at 0.05 × 10–11 per site per year.

Selection Pressures

We identified 16 genes with sites under positive selection in the wild species group (Supplementary Tables 1, 2). These genes comprised two ATP subunit genes (atpB and atpE), three ribosome small subunit genes (rps2, rps3, and rps12), three genes encoding cytochrome b/f complex subunit proteins (petB, petD, and petN), one NADH oxidoreductase gene (ndhG), one DNA-dependent RNA polymerase gene (rpoC2), one gene encoding ribosome large subunit protein (rpl16), and five other genes (ccsA, cemA, rbcL, ycf1, and ycf2). According to the M2 and M8 models, the rps12 gene harbored 28 sites under positive selection, as well as 34 sites in ycf2, six and four sites in ycf1, two and five sites in ndhG, and one site each in the ccsA, cemA, rpl16, rps3, and petB genes. The M8 model detected 15 sites under positive selection in the rps2 gene. However, sites under positive selection in the atpB (five), atpE (two), and rbcL (two) genes were only detected by the M2 model, and the other six genes had only one active site.

We only identified the ribosome large subunit protein (rpl2) gene with sites under positive selection in the semi-wild and cultivated group, where it harbored four sites under positive selection in the M2 model (Supplementary Tables 3, 4).

Diversity and Genetic structure

Seven chloroplast DNA haplotypes were identified in the semi-wild races and 22 in the cultivated upland cotton genotypes (Table 4). The haplotypes diversity (Hd) and π values were slightly higher for the semi-wild races than the cultivated genotypes. STRUCTURE analyses and the ΔK statistic indicated an “optimal” value for K (number of populations modeled) of 2 (Supplementary Figure 1), thereby supporting the existence of two major clusters in the data set (Figure 4). The semi-wild races were primarily assigned to cluster I and the cultivated genotypes to cluster II, whereas the races marie-galante and latifolium had notable fractions assigned to cluster II, thereby suggesting genetic introgression between the two groups.

TABLE 4

PopulationNumber of samplesNumber of haplotypes (H)Hd (SD)π (SD) × 100Number of segregation sitesTheta
Semi-wild races771.000 (0.076)0.00035 (0.00006)1570.196
Cultivated accessions29220.946 (0.035)0.00010 (0.00003)1700.132

Nucleotide diversity and haplotype frequencies for plastid genomes in semi-wild and cultivated accessions of upland cotton.

FIGURE 4

Gene Flow

Patterns of historical and contemporary gene flow were detected between the semi-wild and cultivated upland cotton genotypes. Migrate-n analysis showed that historical gene flow ranged from 149.77 (135.69–164.85) for the semi-wild group to 377.47 (344.25–413.03) for the cultivated group, thereby indicating asymmetric gene flow between the groups. Significant asymmetric contemporary gene flow was also found between the groups, where the values ranged from 0.1110 (0.0612–0.1608) for the semi-wild group to 0.0108 (0.0004–0.0212) for the cultivated group. These results suggest a higher level of historical gene flow during domestication compared with the low level of contemporary gene flow.

Discussion

Evolutionary Relationships

Some previous studies have explored the molecular phylogenetic relationships of cotton, mostly based on a small number of plastids and nuclear DNA markers, as well as the complete chloroplast genome sequence and mitochondrial genome data set of a limited number of cotton species (; Senchina et al., 2003; Wendel et al., 2009; Xu et al., 2012; Wendel and Grover, 2015; , ,; Wu et al., 2018). However, the relationship between cultivated accessions of upland cotton and other species of Gossypium is not clear now. Therefore, we built phylogenetic analyses on 72 cotton plastid genome sequences including wild species, semi-wild races and cultivated accessions of Gossypium, representing the largest number of known cotton species. In the phylogenetic tree, Gossypium species were primarily divided into three large genetic branches. The outer two branches mainly comprised diploid cotton species and the upland cotton clade formed the inner branch. One of the two outside branches included the Australian species with C, G, and K-genomes, American D-genome species, and African E- and B-genome species. Other studies have also shown that species with the G-genome have a common nested relationship with C-genome species, probably due to the frequent capture of chloroplasts in the G. bickii lineage (Seelanan et al., 1999; ). The other outside branch comprised the African F-genome species, Asian–African A-genome species, and American AD-genome wild species and cultivated G. barbadense genotypes. The large internal branch included all of the upland cotton cultivars and semi-wild races. The race latifolium clustered more closely with the upland cotton genotypes, which may suggest a classification error because the race yucatanense is considered the closest progenitor of cultivated upland cotton. Some studies have reported that the maternal donor of the chloroplast genome for the allotetraploid species was the A-genome progenitor (; , ; ), and this was supported by our phylogenetic analysis. The latest research showed that the two A-genome species (G. herbaceum and G. arboreum) have evolved independently with no ancestor-progeny relationship (). In addition, the phylogenetic tree showed that all 13 D-genome species clustered into a single lineage with high support and they were more distantly related to the upland cotton genotypes. Some D-genome species formed closely associated pairs, including G. klotzschianum (D3–k) with G. davidsonii (D3–d), G. harknessii (D2–2) with G. turneri (D10), G. thurberi (D1) with G. trilobum (D8), and G. raimondii (D5) with G. gossypioides (D6). These results are consistent with previous reports of phylogenetic relationships based on nuclear genetic markers and chloroplast genome sequences (; Ulloa et al., 2013; ; Wu et al., 2018; ). The difference in phylogenetic relationships may be caused by the different genetic characteristics of the DNA markers used.

Divergence Time Analysis

We estimated the divergence time of Gossypium species based on the plastid protein-coding sequences. The results showed that the diversification between Gossypium and T. cacao was found to have occurred about 58 Mya, which was consistent with previous inferred results (Wendel, 1989; Senchina et al., 2003; ; ). Interestingly, the divergence time was estimated at 7.7 Mya (95% HPD = 6.3–9.8 Mya) between the B-genome and Australian clade (C + G + K), which was similar to the rapid radiation time calculated for all other cotton branches after differentiation from Australian cotton species (). In addition, the evolutionary time of the cotton ancestors was 11 Mya and cotton species then rapidly differentiated radially, where the differentiation time of most branches was 5–6 Mya. These results were largely consistent with those obtained in other molecular studies (, ,). The differentiation time for the semi-wild races, cultivated upland cotton genotypes, and AD-genome was estimated at 6.25 Mya, and that estimated for the race latifoloum and Gossypium hirsutum cultivar difenmian168 was 0.45 Mya. We also found that the divergence time betweem semi-wild races and cultivated upland cotton accessions were about 3.12 Mya, thereby indicating that they may have differentiated recently. The evolutionary time for the allotetraploid upland cotton accessions was 6.25 Mya (6.4–9.7), which agrees with the results obtained in previous studies (Senchina et al., 2003; Wang et al., 2017; ; ), where it was domesticated at least 4,000 to 5,000 years ago and subsequently subjected to direct selection (Wang et al., 2017). To the best of our knowledge, the present study is the first to use the protein-coding sequences in the chloroplast genome to estimate the divergence dates of the whole Gossypium species including semi-wild races and cultivated upland cotton genotypes, although the results could be improved by larger phylogenetic analyses.

Genetic Mutation

Mutation is the ultimate source of genetic variation, the substrate of evolution (; Zhang et al., 2020). A previous study suggested that the mutation/substitution rates varied between and within genomes (), and that they were influenced by factors such as the nearest neighbor bases, chromosomal position, and the efficiency of the repair systems between the leading and lagging DNA strands. In general, the presence of similar bases or derivatives of similar bases facilitates the base replacement in the DNA repair process, and thus transitions occur more frequently than transversions (). Our results of nucleotide sequence evolution analysis showed that the transition rate was higher than the transversion rate for the cotton genotypes evaluated, which is consistent with previous reports (; ). SNP represents the most common form of polymorphism in biological genomes. Common polymorphisms are effective genetic markers related to biological evolution (Zhang et al., 2020). In the present study, we identified 4,074 SNPs in the Gossypium cp genomes. Among them, there were more SNPs in the intergenic region than the intron region, indicating that intergenic spacer sequences were more variable than intron regions in the plastid genome, which was consistent with the latest research results (Zhang et al., 2020). Furthermore, the dN/dS ratios were larger than 1, thereby indicating that non-synonymous mutations were fixed in the genomes, which may be due to component-driven mutation pressure (). The dN/dS ratios were higher for the semi-wild and cultivated upland cotton genotypes than those determined for the wild cotton species, which may suggest that upland cotton has been subject to very strong artificial selection during domestication. The results of evolutionary rates indicated that the rates of nucleotide substitutions and indels were higher in wild species than the upland genotypes, thereby suggesting that the semi-wild and cultivated upland genotypes might have evolved more slowly after speciation. Due to the influence of artificial domestication, the cultivated genotypes exhibited less variation with fewer mutations. Previous studies have shown that selection can act on the mutation rate (). Moreover, according to our results, the mutation rate was lower for indels than nucleotide substitutions, which is consistent with a previous report (Wu et al., 2018).

Domestication Selection

By the mid-18th century, the coastal colonies of the southeastern United States had developed upland and Sea Island cotton varieties, which showed a long history of cotton domestication and breeding (). Evidence suggested that the domestication and breeding of allotetraploid cotton were superior to A-genomic diploid cotton in yield and quality (). And the allopolyploid cultivated cotton was first domesticated about 5,000 years ago (Yoo et al., 2014). Generally, synonymous and non-synonymous nucleotide substitutions are important markers of gene evolution. In most genes, synonymous nucleotide substitutions have occurred more frequently than non-synonymous substitutions (). The rates of non-synonymous and synonymous substitutions are relatively slow in plant chloroplast genomes because of purifying and neutral selection (; ). In the present study, selection pressure analysis identified 16 genes with sites under positive selection in the wild species group, but only one of these genes (rpl2) was identified in the semi-wild and cultivated group. We conclude that the selection pressure on semi-wild and cultivated cotton species has fewer genes at positive selection sites, whereas the wild species retained adaptive genes and the selected sites increased. These results are generally consistent with those obtained in previous studies of the effects of artificial domestication on selection pressure (). When plants experience relatively large changes in the environment, such as artificial domestication or natural selection, the relaxation of selection for certain characteristics is inevitable (; ). Thereby, humans would expect that natural selection of these features would lose its strength (). The rpl2 domestication selection gene identified in semi-wild and cultivated cotton species may have played an important role in the adaptation of Gossypium to various environments (; ; Wu et al., 2018; ). Moreover, selection pressure analysis for wild and domesticated cotton species can provide novel insights into how human selection has affected duplicated genes in allopolyploids (Yoo et al., 2014; ). It is known that many important crops such as potato, wheat and soybean are obvious polyploids, so studying the genes of allopolyploid cotton may provide new insights into the role of polyploids in crop evolution (Yoo et al., 2014).

Genetic Diversity

Additionally, genetic diversity is the basis of crop improvement (). Therefore, understanding the genetic diversity, structure, and relationships between varieties of upland cotton is very important for breeding (). The semi-wild races exhibited higher nucleotide diversity (Hd = 1.000, π = 0.00035) than the cultivated genotypes (Hd = 0.946, π = 0.00010), thereby suggesting that artificial domestication reduced the chloroplast genetic diversity, which is consistent with a previous report (). The low level of genetic diversity determined in the cultivated upland cotton accessions was primarily due to several genetic bottlenecks during the domestication process (; Wang et al., 2017). Various studies have also suggested that the genetic basis of cultivated upland cotton genotypes is narrow (; ; ), although the diversity of derived cultivars obtained by various breeding methods is still evident. In addition, cotton breeding often involves hybridization and re-selection with a small number of breeding materials, thereby resulting in a loss of genetic diversity (Tyagi et al., 2014). The genetic structure is mainly affected by geographical isolation and genetic exchange isolation (; ). Genetic structure analysis showed that the semi-wild races and cultivated upland accessions were divided into two groups when K = 2. We observed that the seven semi-wild races and cultivated upland accessions exhibited significant admixture, that was, the two semi-wild races Marie-galante and latifolium had notable fractions assigned to cultivated accessions group, which indicated that the race latifolium had closest relationships with cultivated accessions, followed by the race marie-galante race, thereby indicating the introgression of a certain gene between the semi-wild races and cultivated accessions, or possibly germplasm sharing (Tyagi et al., 2014). These results were consistent with a previous study on increasing human-mediated effects leading to significantly genetic introgression (). A previous study also showed that the existence of this mixture may be related to the domestication history and the frequent appearance of superior genotypes in different breeding programs (). China is not a natural cotton-growing region, and thus many cotton genotypes, such as Foster, STV, DPL, Trice, King, and Uganda, have been introduced as extensive genetic sources for upland cotton varieties in China from several overseas sources for improving varieties (; ; ; ). It is important to study the diversity and genetic structure of upland cotton genotypes as well as their relationships to facilitate the conservation and improvement of cotton (). In addition, the genetic diversity and population structure of upland cotton germplasm resources can be effectively used for genetic breeding, and it is of great significance for the systematic utilization of long-term genetic variation of upland cotton (Tyagi et al., 2014).

Genetic Introgression

Ancient gene flow between domesticated varieties and their wild relatives probably occurred historically through seed transmission, and it was possibly influenced by human activities and environmental events (Wegier et al., 2011). In the present study, asymmetric historical gene flow was determined between the semi-wild and cultivated upland genotypes, which is consistent with a previous study (). However, contemporary gene flow was greatly reduced, which may have been due to current isolation. Genetic studies of species in the early stages of domestication have identified multiple domestication origins or high levels of sustained gene flow between wild and cultivated genotypes (). A previous study also suggested that the genetic structure of upland cotton genotypes was weak or an admixture, which may have resulted in a strong historical gene flow (). In general, gene flow is an important factor that affects the population structure over time, where it may reduce local adaptation by homogenizing the populations found in different environments or by spreading harmful alleles between populations. Gene flow might also contribute to the introduction of potential adaptive alleles into populations and increased genetic variation (Sexton et al., 2011; ; Welt et al., 2015). Some studies have also indicated that gene flow from cultivated upland genotypes to wild cotton tetraploid species has increased the risk of extinction for these wild species (Wegier et al., 2011). I.e., some wild cotton species G. tomentosum (in Hawaii), G. mustelinum (in Brazil) and G. darwinii (in Galapagos) were in danger of extinction as a result of hybridization with domesticated tetraploid cotton (; Simard, 2010). In addition, numerous studies have shown that interspecific hybrids (G. hirsutum x G. barbadense) can serve as genetic links for gene transfer from domesticated cotton to other wild relatives (G. darwinii) (; Simard, 2010). This occurred during or after speciation lead to the retention of ancestral polymorphism due to incomplete lineage sorting (; Wilyard et al., 2009), or introgression or introgressive hybridization of previously geographically isolated species resulting from the genetic exchange after secondary contact (; ). Moreover, among the four cultivated Gossypium plants, upland cotton exhibits the highest level of gene flow (Wendel et al., 1992; ), which is related to the strong artificial domestication that it has undergone. The extensive gene flow and/or genetic introgression among cotton accessions might have provided the novel genetic resources of cotton breeding. Therefore, the suitable management and conservation of different cotton species accessions are important in the future.

Conclusion

In conclusion, our phylogenetic analysis confirms the evolutionary relationship within the whole Gossypium, especially the relationships between semi-wild races and cultivated accessions were well resolved. We also identified that the rpl2 gene was positively selected in semi-wild races and cultivated genotypes. Meanwhile, we found that the cultivated genotypes have experienced very strong selection pressure. In addition, we found that the genetic diversity of cultivated accessions was low compared to wild ones due to artificial domestication. Through the analyses of genetic structure and gene flow, we concluded that there was a certain gene introgression between semi-wild races and cultivated accessions. The present research provided novel genetic resources for cotton breeding, as well as novel molecular mechanisms insights for the evolution and domestication of cotton species.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Statements

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Author contributions

TZ, NW, and Z-HL: data curation and writing – original draft. YW: formal analysis. TZ, NW, and X-FM: investigation. YW, X-LZ, B-GL, WL, J-JS, C-XW, and AZ: methodology. X-FM: resources and validation. TZ and NW: software. Z-HL: supervision and writing – review & editing. All authors contributed to the article and approved the submitted version.

Funding

This research was funded by grants from the National Key R&D Program (2021YFF1000100), the Innovation Project of the Chinese Academy of Agricultural Sciences (CAAS-ASTIP-ICR-KP-2021-01), the Xinjiang Tianshan Talents Program (2021), the Central Public-interest Scientific Institution Basal Research Fund (Y2021XK12), the Project of Introduction High-level Talents in Xinjiang Uygur Autonomous Region Flexible Talents (2020), the Shaanxi Science and Technology Innovation Team (2019TD-012), and the Key Program of Research and Development of Shaanxi Province (2022ZDLSF06-02), and the Public Health Specialty in the Department of Traditional Chinese Medicine (2019-39).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2022.873788/full#supplementary-material

References

Summary

Keywords

cotton, domestication selection, gene flow, genetic relationship, nucleotide evolution

Citation

Zhou T, Wang N, Wang Y, Zhang X-L, Li B-G, Li W, Su J-J, Wang C-X, Zhang A, Ma X-F and Li Z-H (2022) Nucleotide Evolution, Domestication Selection, and Genetic Relationships of Chloroplast Genomes in the Economically Important Crop Genus Gossypium. Front. Plant Sci. 13:873788. doi: 10.3389/fpls.2022.873788

Received

11 February 2022

Accepted

24 March 2022

Published

15 April 2022

Volume

13 - 2022

Edited by

Lin-Feng Li, Fudan University, China

Reviewed by

Jie Qiu, Shanghai Normal University, China; Nian Wang, Huazhong Agricultural University, China; Xiongming Du, State Key Laboratory of Cotton Biology, Cotton Institute of the Chinese Academy of Agricultural Sciences, China

Updates

Copyright

*Correspondence: Zhong-Hu Li, Xiong-Feng Ma,

†These authors have contributed equally to this work

This article was submitted to Plant Systematics and Evolution, a section of the journal Frontiers in Plant Science

Disclaimer

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.

Outline

Figures

Cite article

Copy to clipboard


Export citation file


Share article

Article metrics