Genome duplication (polyploidy) is among the most dramatic mutational processes in nature, causing myriad saltational changes at the cellular and organismal levels (Doyle and Coate 2019;Bomblies 2020;Fernandes Gyorfy et al. 2021), and is associated with consequential phenomena ranging from crop domestication (Renny-Byfield and Wendel 2014;Qi et al. 2021) to cancer progression (Matsumoto et al. 2021). Polyploidy is especially common in the angiosperms, with all extant species having experienced at least one or more polyploidy events during their evolutionary history (Jiao et al. 2011), and at least 30% of currently recognized species having a polyploidy event in the recent past (One Thousand Plant Transcriptomes Initiative 2019).
Novel evolutionary patterns created by polyploidy at the genic (e.g., neofunctionalization, subfunctionalization, and loss; Kuzmin et al. 2021) and genomic (e.g., homoeologous recombination; Mason and Wendel 2020) levels have been well documented across taxa, including the frequent asymmetry of these responses with respect to coresident genomes in a polyploid nucleus. Nonetheless, many questions remain regarding the effects of natural selection on polyploid relative to diploid genomes (Baduel et al. 2019;Monnahan et al. 2019) and the interplay between these novel evolutionary patterns and the long-term trajectories of genome evolution (Qi et al. 2021) following polyploidization (e.g., biased fractionation).
One of the earliest predictions about natural selection in polyploids relative to diploids is that putatively deleterious mutations may accumulate faster due to the masking effect of completely or partially recessive deleterious mutations in duplicated genes (Haldane 1932;Hill 1970;Bever and Felber 1992). Only recently, however, have these predictions begun to be evaluated in young allopolyploids such as Arabidopsis kamchatica (Paape et al. 2018) and Capsella bursa-pastoris (Douglas et al. 2015;Kryvokhyzha, Milesi, et al. 2019;Kryvokhyzha, Salcedo, et al. 2019), and autotetraploid Arabidopsis arenosa (Monnahan et al. 2019). Because the number of deleterious mutations is strongly influenced by shifts in demography and mating system (Brandvain and Wright 2016), which may coincide with polyploid formation (Grant 1981;Barringer 2007), a clear link between ploidy level and the accumulation of deleterious mutations is challenging to demonstrate in natural polyploid populations.
ß The Author(s) 2022. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons. org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.
The cotton genus (Gossypium) represents one of the best studied allopolyploid systems (Wendel and Grover 2015;Hu et al. 2021). The genus includes approximately 45 currently recognized diploid species classified into eight genome groups (A-G, and K), and seven allopolyploid species resulting from a single (Grover et al. 2012) allopolyploidization event approximately 1-2 million years ago (Mya) between members of the A and D genome groups (Wendel 1989). Although the most closely related extant species of these two progenitor diploids are found in southern Africa and Northern Peru, respectively, the polyploids are only found in the Americas (fig. 1). Most wild populations, including those of the two independently domesticated species G. hirsutum (AD 1 ) and G. barbadense (AD 2 ), occur in small, isolated populations on islands or in coastal regions. Subsequent to their initial domestication 4,000-8,000 years ago in the Yucatan Peninsula (AD 1 ) and NW S. America (AD 2 ), respectively, the ranges of the two domesticated species rapidly expanded to encompass much of the American tropics and subtropics and then spread globally with the rise of the international cotton fiber trade (Yuan et al. 2021).
Here, we describe the evolutionary trajectory of deleterious mutations in two wild diploid and six wild allopolyploid cotton species (all descended from a single allopolyploidization event), with a focus on how allopolyploidization and speciation have shaped the number and genomic distribution of deleterious mutations. We use two methods to estimate the strength of selection at the amino acid and nucleotide level and show support for a nearly century-old hypothesis that polyploids accumulate mutations faster than their diploid progenitors. We demonstrate that, in agreement with this hypothesis, polyploidy has the greatest influence on strongly, rather than moderately or mildly, deleterious mutations. We also find that deleterious mutations accumulate asymmetrically between the two coresident subgenomes in the allopolyploid nucleus, indicating that these masking effects may act unequally between the subgenomes. In total, our results support theoretical predictions that allopolyploidy can lead to a faster rate of deleterious mutation accumulation through masking of recessive deleterious variants, and that the relationship of the rate of deleterious mutation accumulation between subgenomes and their progenitor diploids is complex, even when comparing identical pairs of single-copy homoeologs among lineages.
To investigate patterns of deleterious mutations, we viewed our data at three phylogenetic depths: single nucleotide mutations that occurred anywhere on the global phylogenetic tree (fig. 2A-D), mutations that emerged since the divergence of each subgenome from its respective diploid progenitor Using the curated set of 8,884 pairs of homoeologous genes, we found no evidence for differences in the rate of synonymous mutation accumulation in diploids versus polyploids at any phylogenetic depth (fig. 2A andE), although differences can be found within the polyploid species (fig. 2I), with G. mustelinum (AD 4 , Orange) and G. darwinii (AD 5 , Yellow) having consistently lower rates than the rest of the clade, and in both subgenomes. When viewing mutations that have accumulated since the divergence of the earliest polyploid lineage (fig. 2A-I), there is an asymmetry between subgenomes in the rate of synonymous site changes, with the Dt ("t" denoting "tetraploid") subgenome containing a moderately higher number of synonymous mutations than the At subgenome for all species. This difference potentially indicates a higher mutation rate or relaxed background selection in genic regions of the Dt subgenome compared with homoeologous genic regions of the At sugenome, and is consistent with previous analysis finding that genes in the Dt subgenome are evolving faster than genes in the At subgenome in five allopolyploid cottons (Chen et al. 2020). E-H) shows mutations that are variable within each subgenome and its associated progenitor diploid species. For example, in panel (F), the yellow solid circles indicate the number of derived nonsynonymous mutations in the At subgenome of AD5 (G. darwinii) that have occurred since its divergence from the A1 diploid (left half of panel), and in the Dt subgenome of AD5 (G. darwinii) that have occurred since its divergence from the D5 diploid (right half of panel). The bottom row (panels I-L) shows mutations that originated postpolyploidy and are variable within the polyploids; for example, in panel (L), the purple solid circles indicated the number of derived deleterious mutations identified by BAD_Mutations that have occurred in the At subgenome of AD3 (G. tomentosum) since its divergence from the At subgenome of AD4 (G. mustelinum; left half of panel), and in the Dt subgenome of AD3 (G. tomentosum) since its divergence from the Dt subgenome of AD4 (G. mustelinum; right half of panel). Panels (A), (B), (E), (F), (I), (J): the y axis for both synonymous and nonsynonymous represents the sum of the derived allele frequencies, interpreted as the average number of derived mutations in that category in each species. Panels (C), (G), and (K): the y axis indicates the GERP Load of each species, calculated as the sum of (derived allele frequencyÂGERP Score) for all variant positions with GERP>0. Panels (D), (H), and (L): the y axis shows the total number of deleterious mutations in each species, calculated by BAD_Mutations with Bonferroni-corrected significance (see Materials and Methods)and indicates the average number of deleterious mutations in each species at a given phylogenetic depth. For panels (E-H), comparisons between subgenomes cannot be made because the D 5 diploid is more distantly related to the D subgenome than the A 1 diploid is related to the A subgenome. Therefore, we would expect a larger number of derived mutations in D than A simply due to evolutionary history rather than to polyploidization per se.
Deleterious Mutations Accumulate Faster in Allopolyploids . doi:10.1093/molbev/msac024
In contrast to the relative homogeneity in rates of synonymous substitution among diploids and polyploids, rates of nonsynonymous mutation accumulation differed significantly at all phylogenetic depths. Notably, at the deepest phylogenetic depth (fig. 2B), estimates for the number of derived nonsynonymous mutations in the diploids G. herbaceum (A 1 , Red) and G. raimondii (D 5 , Black) did not differ between subgenomes, indicating that any mapping biases or erroneous SNP calls in these samples were removed by our variant filtering criteria. Gossypium raimondii (D 5 ) contained more derived nonsynonymous mutations than did G. herbaceum (A 1 ), and this lineage-specific difference was reflected in the Dt and At subgenomes as well. When lineage-specific effects that arose from the long, shared ancestry between the subgenomes and their progenitor diploids were removed (fig. 2F), a clear distinction between the rates of nonsynonymous mutations between diploids and their respective subgenomes in the allopolyploids becomes apparent. In all polyploids, the At subgenome contained between 25% and 58% more nonsynonymous mutations than G. herbaceum (A 1 , Red), and the Dt subgenome contained 17-36% more than G. raimondii (D 5 , Black). These results demonstrate that even though the rates of synonymous mutation accumulation did not differ significantly between the diploids and polyploids, polyploidy significantly increases the rate of nonsynonymous substitution accumulation. Finally, when restricting our attention to only those mutations that have arisen following polyploid formation (fig. 2J), the lineagespecific patterns observed for nonsynonymous mutations were largely identical to the patterns of synonymous mutations. For example, G. mustelinum, (AD 4 , Orange) consistently had the lowest number of derived mutations in both subgenomes. Notably, however, the Hawaiian Island endemic G. tomentosum (AD 3 , Purple) has a higher number of derived nonsynonymous mutations than expected based on the patterns of synonymous mutations, potentially reflecting the population bottleneck associated with its origin following long-distance dispersal to the Hawaiian Islands (see Discussion). In summary, polyploidy significantly enhances the rate of nonsynonymous mutation accumulation in all Gossypium allopolyploids, and does so asymmetrically across coresident genomes.
Because the fitness effects of most nonsynonymous mutations can vary widely, from neutral to lethal, we asked if the elevated rate of nonsynonymous mutations observed in polyploid Gossypium reflects an increase in neutral or nearly neutral nonsynonymous mutations, or if instead this elevation is attributable to a greater accumulation of deleterious mutations. To address this, we used two approaches to estimate whether a mutation was deleterious: BAD_Mutations and GERPþþ. BAD_Mutations performs a likelihood ratio test from a gene-specific multispecies alignment to determine if a mutation at a particular nonsynonymous site is potentially deleterious, whereas GERPþþ uses a genome-wide multiple sequence alignment (i.e., agnostic to genic regions) to estimate the degree of conservation at a particular site in the genome (including noncoding and synonymous sites). Notably, because one of the hallmark long-term processes following polyploidy is pseudogenization (Wendel 2015), recently pseudogenized sequences may still display some degree of conservation across the multiple sequence alignment, but may not be inherently deleterious. Therefore, to avoid inflating estimates of deleterious mutations in polyploids compared with diploids, we used GERP only within the exonic regions of the 8,884 homoeologs. Additionally, although GERP can score the degree of deleteriousness of a mutation, BAD_Mutations can only classify variants into deleterious or not deleterious. Therefore, the values shown in figure 2D, H, and L represent the sum of the allele frequencies of derived deleterious mutations, similar to the values for figure 2A,E, and I and figure 2B, F, and J. For analysis with GERP, we used the GERP load, which incorporates the deleteriousness of each variant into the score, summing the frequency of each derived allele multiplied by its GERP score (see Rodgers-Melnick et al. [2015] and Wang et al. [2017]).
As shown in figure 2, both of the foregoing analyses demonstrate that deleterious mutations accumulate in polyploids in a manner similar to nonsynonymous mutations, suggesting that the difference in nonsynonymous sites cannot be wholly attributed to putatively neutral or nearly neutral alleles. For example, there is remarkable consistency in the patterns of deleterious mutations that have accumulated since the divergence of the diploid from its respective diploid progenitor in both the count of nonsynonymous substitutions (fig. 2F), the GERP load (fig. 2G), and the number of deleterious mutations (fig. 2H). In all three columns, the diploids show fewer accumulated alleles than the polyploids, G. tomentosum (AD 3 , Purple) shows the highest number of all the polyploids, and G. mustelinum (AD 4 , Orange) shows the fewest of all the polyploids.
An interesting pattern arises when comparing estimates of the GERP load (fig. 2G) and number of deleterious mutations (fig. 2H) between diploids and their closely related subgenomes: Although the total number of deleterious mutations in the At subgenomes was 52-99% higher in the polyploids than the diploids (fig. 2H), the GERP load in the polyploids was only 13-42% higher (fig. 2G). Similar patterns were found in the Dt subgenome, with 34-66% more deleterious mutations in the polyploids than the diploid, but only a 9-13% increase in GERP load. This discrepancy could reflect inherent differences in the types of sequences used and how deleteriousness is quantified between the two methods, suggesting that the use of multiple analytical tools for detection of genetic load may yield more nuanced insights than either method on its own (See Discussion).
Although deleterious mutations are accumulating faster in polyploids relative to diploids, it is not obvious whether this increased rate is different from the increased rate of accumulation of nonsynonymous mutations. To test this, we compared, among ploidy levels, the total proportion of Conover and Wendel . doi:10.1093/molbev/msac024 MBE nonsynonymous mutations that were considered deleterious by BAD_Mutations (fig. 3). For mutations that originated since the divergence of the A and D diploids (fig. 3A), the proportion of nonsynonymous sites that are deleterious is roughly 2% higher in polyploids than in diploids, despite the shared evolutionary history of more than 4 My between each subgenome and their respective diploid progenitors. Notably, as similarly shown in figure 2, the proportion of nonsynonymous mutations that are inferred to be deleterious in both diploids is equivalent when mapped to either subgenome, indicating that our filtering criteria did not differentially exclude deleterious or nondeleterious mutations with respect to which subgenome the diploid reads were mapped.
At shallower phylogenetic depths (fig. 3B), the difference between diploids and polyploids becomes even clearer, with polyploids exhibiting 3-4% higher proportions of deleterious mutations in the Dt subgenome and 5-12% higher in the At subgenome than their respective diploid progenitors. The most unbiased and straightforward comparison of the asymmetry in strength of the masking effect of duplicate genomes between the two subgenomes of allopolyploid cottons is provided by mutations that have occurred following polyploidization (fig. 3C). Here, the At subgenome of all species contain a 2-3% high proportion of deleterious mutations than the Dt subgenome, indicating that differences exist in the strength of the masking effect between the two homoeologous subgenomes that have resided in the same nucleus for over a million years. This pattern is also observed when deleterious mutations are mapped onto the phylogeny (supplementary fig. 3, Supplementary Material online). Additionally, there is more variation among species in the At subgenome than in the Dt subgenome, although the patterns in this respect are not simple. The amount of subgenomic asymmetry is smallest in G. darwinii (AD 5 , Yellow) from the Galapagos Islands, and largest in the Brazilian endemic and inland species G. mustelinum (AD 4 , Orange), indicating that asymmetries between subgenomes of the same species may vary within a single clade of allopolyploids.
Demography is a potential confounding factor in estimating the rate of deleterious mutation accumulation. Shifts in demography are known to complicate inferences of the strength of selection and genetic load (Brandvain and Wright 2016); for example, even in one of the best studied demographic shifts, the Out of Africa migration in humans, several papers (Lohmueller et al. 2008;Gazave et al. 2013;Simons et al. 2014;Henn et al. 2016;Simons and Sella 2016) have reached seemingly contradictory conclusions on whether genetic load has increased as a result of these shifts in demography (but see Lohmueller [2014]). The pattern of deleterious mutation accumulation has also been welldocumented in bottlenecks and population growth associated with domestication in crops such as maize (Wang et al. 2017), soybean and barley (Kono et al. 2016), sorghum (Lozano et al. 2021), cassava (Ramu et al. 2017), and rice (Liu et al. 2017).
Polyploidy is typically associated with a population bottleneck (Grant 1981;Barringer 2007), but because the genetic diversity of both the diploid and polyploid species in this study is low (table 1), demographic modeling of the depth or duration of population bottlenecks and range expansion following polyploid formation is not straightforward. Generalized patterns of the effects of demography on deleterious mutations, however, can serve as a null expectation to test if our data follow the same trends observed under varying demographic scenarios, as explained in the following. Demographic shifts, including population bottlenecks and expansions, have a large influence on the accumulation of deleterious mutations. According to the nearly neutral theory (Ohta 1992), the fate of deleterious mutations is determined by genetic drift instead of selection when the selection coefficient (s) of deleterious mutations is less than or equal to 1/ (2 N e ), where N e is the effective population size. The reduction of N e during a population bottleneck would therefore allow weakly deleterious mutations to escape purifying selection (i.e., to behave as if they were neutral), whereas strongly deleterious mutations with a selective coefficient greater than 1/ (2N e ) would still be removed by purifying selection. On the In both demographic scenarios, we expect that mildly or moderately deleterious mutations would be most differentially affected, whereas strongly deleterious mutations would consistently be removed by purifying selection. Based on this theory, if the differences in the number of deleterious mutations we see between diploids and polyploids are due to demography, then we would expect to see most of that difference reflected in mildly, rather than strongly, deleterious mutations. In contrast, if masking of deleterious alleles in polyploids is driving a higher rate of accumulation relative to diploids, this pattern will not be observed.
To test if our data were consistent with changes in demography, we first asked if there was a correlation between the degree of deleteriousness of a mutation (as measured by GERP) and its relative increase in the polyploids compared with the diploids. To answer this question, we plotted the relative change of deleterious mutations in each subgenome relative to its most closely related diploid progenitor. We plotted this relative change for three different degrees of deleteriousness-strongly deleterious mutations (4 < GERP 6), moderately deleteriousness (2 < GERP 4), and mildly deleteriousness (0 < GERP 2) deleterious (fig. 4). We found that in both subgenomes of all six polyploids, when comparing mutations that had originated after the divergence of the diploid from its respective subgenome in the allopolyploids, strongly deleterious mutations accumulated at a faster rate relative to diploids than did moderately or mildly deleterious mutations, which is inconsistent with expectations under a demographic change model alone. We also observed this change under both an additive and recessive model of dominance (supplementary fig. 5, Supplementary Material online), and also when evaluating deleteriousness using the "masked constraint" estimates generated by BAD_Mutations (supplementary fig. 6, Supplementary Material online, but see Kono et al. [2018] and Chun and Fay [2009] for cautions on interpreting this value as deleteriousness). In total, the rate of accumulation among mutations with different inferred degrees of deleteriousness do not suggest that the patterns we see can be explained solely by demographic changes, but that the masking effect of duplicated genes may play an important role in the determining the fate of deleterious mutations in allopolyploids.
One of the earliest hypotheses regarding mutation accumulation in allopolyploids dates back to Haldane (Haldane 1932) where he posits that in allopolyploids, "one gene may be altered without disadvantage, provided its functions can be performed by a gene in one of the other sets of chromosomes." Allopolyploids are therefore predicted be able to tolerate a higher mutational load than their diploid relatives, and putatively deleterious mutations may accumulate faster in polyploids than in their diploid relatives due to the masking effect of recessive or incompletely dominant deleterious alleles. Here, we demonstrate that these predictions are true in allopolyploid cottons. All polyploids in Gossypium harbor more mutations at phylogenetically conserved sites than do their closest diploid progenitors, as determined by two different methods of detecting deleterious mutations. We also find that the proportion of all nonsynonymous mutations that are inferred to be deleterious is higher in polyploids than in their diploid progenitors and that polyploidy has the greatest effect on strongly deleterious (and, inferentially, more recessive; Eyre-Walker and Keightley 2007; Huber et al. 2018) mutations. Thus, using the power of comparative phylogenetics and genomics combined with analytical methods for detection of deleterious mutations, we For mutations that originated since the divergence of each subgenome from its diploid progenitor, we plotted the relative increase in deleterious alleles across three GERP load categories: mildly deleterious (0<GERP 2; light gray), moderately deleterious (2<GERP 4; gray), and strongly deleterious (4<GERP 6; black). We used the diploid as the reference population, meaning that the relative increase of GERP load in the diploid is always equal to one for all categories. In both subgenomes of all polyploids, strongly deleterious mutations had the greatest relative increase compared with the diploids, followed by the moderately deleterious mutations, and finally, mildly deleterious mutations. This pattern does not fit the expected patterns under demographic models alone, where most of the changes between two populations should be seen in mildly or moderately deleterious mutations. However, under a model where recessive deleterious mutations are masked by their homoeologs, we would expect that strongly deleterious mutations would accumulate faster than moderately or mildly mutations (i.e., the pattern we see here) due to the correlation between the recessivity of a mutation (h) and its selection coefficient (s).
Estimating the strength of natural selection and genetic load is notoriously challenging (Lohmueller 2014) and is complicated by shifts in effective population size (including bottlenecks and expansions), mating systems, and effective recombination rates, among other life-history and demographic factors (Brandvain and Wright 2016). Here, we illuminate an additional relevant consideration, that is, wholegenome duplication. Yet many of the considerations for populations that are not in demographic equilibrium also apply to Gossypium. Diversification in the cotton tribe (Gossypieae) has been characterized by numerous long-distance dispersal events (Grover et al. 2017), including the one from Africa to the Americas 1-2 Mya that led to the evolution of allopolyploid Gossypium. We note that in the Hawaiian Islands endemic G. tomentosum, the total number of synonymous substitutions is not significantly different from the rest of the polyploids, but the number of nonsynonymous and deleterious mutations is significantly increased, suggesting that the genetic bottleneck associated with island dispersal has elevated the number of deleterious mutations compared with the rest of the polyploids.
Although demographic changes upon polyploid formation have been shown to change the number and frequency of deleterious mutations in other systems (Douglas et al. 2015;Paape et al. 2018;Baduel et al. 2019;Kryvokhyzha, Milesi, et al. 2019;Kryvokhyzha, Salcedo, et al. 2019), we show here that the patterns of mutation accumulation in Gossypium cannot be explained by demography alone, and that the data are more consistent with the nearly century-old hypothesis that recessive deleterious mutations can accumulate faster in allopolyploids due to the masking effect of duplicated genes and lack of recombination between subgenomes (Haldane 1932). Specifically, we show that strongly (and, hence, more recessive; Morton et al. 1956;Mukai et al. 1972;Eyre-Walker and Keightley 2007;Agrawal and Whitlock 2011;Huber et al. 2018) deleterious mutations accumulate faster in polyploids compared with diploids than moderately or mildly deleterious mutations, and that this pattern is inconsistent with demographic shifts or longterm change in population size (fig. 4 and supplementary fig. 5, Supplementary Material online).
One of the elegant attributes of a clade of allopolyploid genomes derived from a single polyploidization event is that they offer a remarkable natural experiment for comparing subgenomes that have resided within the same nucleus for, in the case of Gossypium, approximately 1.5 My. Once an allopolyploid is established, each subgenome is subjected to identical external or population-level factors, including demography, mating systems, and environmental and ecological conditions, as well as internal cellular processes, including identical DNA replication and recombination machinery. These features remove many of the confounding factors that may influence the genetic load and provide a simple comparative context for revealing evolutionary forces that might differentially affect coresident genomes or homoeologs.
An unexpected finding of our analyses is the striking asymmetry in the proportion of all nonsynonymous mutations that are inferred to be deleterious between the two subgenomes of all allopolyploid species in Gossypium. We found that the At subgenome of all species contains 2-3% more nonsynonymous mutations that are inferred to be deleterious (fig. 3) even when only considering mutations that have arisen following the earliest allopolyploid diversification events, and correcting for removing the biases of unequal phylogenetic distances to each subgenome's model progenitor diploid. Our work adds to a growing recognition that the two coresident subgenomes in cotton allopolyploids may be shaped asymmetrically by evolutionary processes, including interspecific introgression and selection under domestication (Fang, Guan, et al. 2017;Fang, Wang, et al. 2017;Chen et al. 2020;Yuan et al. 2021), and that this phenomenon also extends to other important allopolyploid crop plants, including wheat (Pont and Salse 2017;Jiao et al. 2018) and Brassica (Tong et al. 2020).
Teasing apart the genesis of differential subgenomic responses to selection is rendered challenging by several factors independent of phylogeny. We note, for example, the relevant example of the recently formed allopolyploid Capsella bursa-pastoris and its diploid progenitors, where consistent asymmetries in genetic load are reported between the subgenomes (Kryvokhyzha, Milesi, et al. 2019;Kryvokhyzha, Salcedo, et al. 2019) the differences likely reflect the dramatically different mating systems of the progenitors, in which the subgenome with the higher genetic load originated from an obligate outcrosser, C. grandiflora (N e ¼ 800,000), whereas the subgenome with the lower genetic load derives from the predominantly selfing C. orientalis (N e ¼ 5,000) (Douglas et al. 2015). In another recently formed (20-250 ka) allopolyploid, Arabidopsis kamchatica, no asymmetry in the distribution of fitness effects between subgenomes was found, although it was observed that each subgenome of the allopolyploid contained more neutral and fewer deleterious alleles than either of the diploid progenitors (Paape et al. 2018). It is unclear, however, whether this shift was due to allopolyploidy per se or if it reflects the transition from an obligate outcrossing to a mating system with some degree of inbreeding, with a concomitant purging of partially or completely recessive deleterious alleles, as shown in several other systems (Arunkumar et al. 2015;Roessler et al. 2019). In Gossypium, all species have similar mating systems and a canonical outcrossing floral morphology including highly exserted styles and stigmas. Population sizes often are small, however, likely leading to relatively high levels of generalized inbreeding. At present, however, no data exist that address these considerations.
One possible interpretation of our results is that Gossypium polyploids are less fit than their closely related diploid Deleterious Mutations Accumulate Faster in Allopolyploids . doi:10.1093/molbev/msac024 MBE progenitors because they harbor more deleterious mutations in their genomes, especially mutations that have already been driven to fixation. We note, however, that the fitness effects of a mutation may change as a result of the genetic (e.g., epistasis) or environmental (e.g., local adaptation, conditional neutrality) context in which it occurs (Huang et al. 2021). Comparative genomics techniques used to infer deleterious mutations at phylogenetically conserved sites, as employed here, cannot identify these shifts in fitness effects (Huber et al. 2018), and also frequently incorrectly identify beneficial mutations as deleterious (Chen et al. 2020). Therefore, an additional possibility is that mutations in polyploids that occur at phylogenetically conserved sites may not actually have a deleterious effect on fitness as they do in diploids. Inferring the genetic load of a population simply by counting the number of deleterious variants also assumes that all alleles contribute independently to the total genetic load of a population. However, because of the functional overlap of duplicated genes and, in most cases, absence of recombination between homoeologous chromosomes in an allopolyploid, a recessive deleterious mutation can never be present in all four copies of a gene and thus may be invisible to selection because of the masking effect of its homoeologous partner.
An important takeaway from this study is that recessive deleterious mutations in allopolyploids, at least at some loci, may actually accumulate in a manner more similar to neutral mutations, presumably because of the lack of recombination between subgenomes and, hence, the inability of purifying selection to "see" the negative effects of these mutations. Because these recessive deleterious mutations escape the effects of purifying selection, many traditional tests for detecting positive and negative selection (e.g., dN/dS, p N /p S ) may be biased when comparing a polyploid to diploid because the polyploid would be expected to accumulate putatively deleterious sites more quickly (and maintain a higher genetic diversity at nonsynonymous sites) than their diploid relatives. This increased dN/dS value in allopolyploids compared with diploid progenitors was recently shown in five allopolyploid systems in addition to Gossypium (Sharbrough et al. 2022), and duplicated genes associated with an ancient polyploid event in Brassica rapa were shown to contain higher amounts of genetic variation than nonduplicated genes (Qi et al. 2021), indicating this phenomenon is likely not restricted to Gossypium and should be taken into consideration for other allopolyploid/diploid comparisons. Additionally, although this bias will be most notable in genome-wide comparisons, it may also be evident at the individual gene level, given that both homoeologs are still present and retain some degree of functional overlap (although the robustness of attributing this bias to a ploidy effect acting on a single gene is expected to be low).
Another important implication of this finding is that allopolyploidy (or gene duplication in general) may play an important and underrecognized role in determining how selection acts on new mutations, notwithstanding the burgeoning literature on fates of duplicate gene evolution (Conant et al. 2014;Shi et al. 2020;Veitia and Birchler 2022). The evolutionary trajectory of new mutations will largely be dependent on the selection coefficient (s) acting on that locus, and the dominance coefficient (h), defined as the proportion of the fitness cost that a mutation harbors when in a heterozygous state. In allopolyploids, however, the evolutionary fate of new mutations may be determined not only by allelic dominance at that locus, but also by the interaction arising from the coexistence of its homoeologous locus, a term we call "homoeologous epistatic dominance." The relationships between this homoeologous epistatic dominance, allelic dominance, and the selection coefficient are likely complicated and potentially heavily influenced by other biological considerations such as biased expression of homoeologs, sub-or neofunctionalization, and homoeologous recombination, among others. Moreover, notwithstanding these polyploidy-specific effects, even the genome-wide relationships between two of these factors, allelic dominance and the selection coefficient, have only been modeled using genomic data in the past few years (Huber et al. 2018).
Nonetheless, understanding how this homoeologous epistatic dominance impacts the fitness effects of new mutations is an unexplored aspect of polyploid genome evolution, and it is not yet clear whether this will equally affect advantageous and deleterious variants. How homoeologous epistatic dominance operates with respect to functional properties arising from considerations such as gene balance (Veitia and Birchler 2022), dosage effects (Conant et al. 2014), structural and functional entanglement (Kuzmin et al. 2020(Kuzmin et al. , 2021)), and intersubgenomic cis-and trans-effects (Bottani et al. 2018;Hu and Wendel 2019) would seem to represent important avenues for understanding how natural selection operates differently in polyploids compared with diploids. From an applied perspective, these insights could be important in agriculture, particularly because so many of our most important crop plants have a recent history that includes polyploidy (Renny-Byfield and Wendel 2014), and segregating patterns of genome fractionation have the potential to serve as targets of selection in crop improvement (Hufford et al. 2021).
We used whole-genome sequencing data from 46 individuals in Gossypium, including between two and ten individuals from each of eight species. Included in our sampling was six polyploid species originating from a single polyploidization event 1-2 Mya (Wendel 1989;Hu et al. 2021), two diploid species representing models of the genome donors to the allopolyploids (A and D), and three species from Australia that served as outgroups for polarizing mutations into ancestral and derived states. These sequences were previously described (Yuan et al. 2021), and SRA codes for all 46 resequenced individuals are listed in supplementary table 1, Supplementary Material online. For G. hirsutum, we randomly chose ten accessions that were classified in the "Wild" population from Yuan et al. (Yuan et al. 2021), and for the other species, we chose all accessions available that did not show evidence of being mislabeled, as determined by a PCA plot. Conover and Wendel . doi:10.1093/molbev/msac024 MBE After the data were downloaded from NCBI, adapter sequence removal and quality score filtering of FASTQ reads were performed using Trimmomatic v0.36 (Bolger et al. 2014) using the parameters "LEADING:28 TRAILING:28 SLIDINGWINDOW:8:28 SLIDINGWINDOW:1:10 MINLEN:65 TOPHRED33." Trimmed reads from each polyploid sample were mapped to the 26 chromosomes of the G. hirsutum reference genome (Saski et al. 2017), and reads from each diploid sample were mapped to each subgenome separately to avoid competitive mapping of the diploid reads against a tetraploid reference genome. Reads from the three outgroup species were separately mapped to both subgenomes to ensure that reads were not filtered out for mapping to multiple parts of the genome. All mapping was done using bwa-mem v0.7.17 (Li and Durbin 2009) and only uniquely mapping paired reads (-F 260 flag) that were mapped in their proper orientation (-f 2 flag) were retained using Samtools v1.9 (Li et al. 2009) before the files were sorted and converted to bam files. Using the Sentieon (Kendig et al. 2019) SNP Calling program, gVCF files were generated, and joint genotyping was performed using the GVCFtyper algorithm (see Github repository for full scripts). Variant filtering was performed using GATK v4.0.4.0 using the filter expression "QD < 2.0 k FS > 60.0 k MQ < 40.0 k SOR > 4.0 k MQRankSum < À12.5 k ReadPosRankSum < À8.0." For each species (excluding the outgroup species, and treating G. stephensii, and G. ekmanianum as a single species), we nullified any variant call in which all individuals were heterozygous to remove any collapsed genomic region in the reference genome or paralogous regions that were not present in the reference genome. We treated G. stephensii and G. ekmanianum as a single species because we only sampled two individuals of G. stephensii, so removing any sites in which both individuals were heterozygous errantly removed real variants that were not due to paralogy mapping issues. All scripts for generating and filtering variant calls are located on our GitHub repository (https://github.com/conJUSTover/ Deleterious-Mutations-Polyploidy).
We used the pSONIC pipeline (Conover et al. 2021) to identify syntenically conserved homoeologs in the G. hirsutum reference genome, and kept only homoeologous pairs that were less than 5% different in their total annotated CDS length. To remove homoeologous pairs that may have experienced homoeologous exchange events (though there is scant evidence for this in Gossypium; Salmon et al. 2010;Flagel et al. 2012;Chen et al. 2020), we removed any pair in which the proportion of the reads from the two progenitor diploid genomes (termed At and Dt in the allopolyploid, the "t" indicating "tetraploid") did not meet the expected 2:2 ratio, similar to previous analyses of homoeologous exchange events in other allotetraploids (Bird et al. 2021). Average read depth of CDS regions was determined by bedtools2 v.2.27.1 (Quinlan and Hall 2010). Briefly, for a single homoeologous pair, we calculated the average read depth of the At homoeolog divided by the sum of the average read depth of both homoeologs for each individual. We removed any homoeologous pair in which this fraction was less than 37.5 or greater than 62.5 in any individual. We expect any HEs that result in a 0:4 At: Dt copy number to contain 0% At reads/ total reads; HEs that result in 1:3 At: Dt copy number should have a 25% At reads/total reads; HEs that result in a 3:1 At: Dt copy number should have a 75% At reads/total reads; HEs that result in a 4:0 At: Dt copy number should have a 100% At reads/total reads; and no HE (i.e., 2:2 At: Dt copy number) would result in a 50% At reads/total reads. We used the midpoints between the "No HE" and the 1:3 and 3:1 copy numbers as cutoff points. Because we compared the read depth for homoeologous pairs within each individual, we expect there to be no difference in read depth between homoeologs, and differing read depths among individuals should not introduce any bias in our gene filtering criteria. This filtering resulted in 8,884 homoeologous pairs (17,768 genes) being analyzed further.
Nonreciprocal homoeologous exchanges (i.e., homoeologous gene conversion) could also bias the estimates of the genetic load in a way that is not related to new mutation following polyploidization or speciation. To control for positions in these non-HE homoeologs that may be influenced by gene conversion, we linked variant positions between homoeologs in the following way. We first performed pairwise alignments of the CDS sequences using MACSEv2 (Ranwez et al. 2011(Ranwez et al. , 2018)), which aligns CDS sequences in accordance with their translated amino acid sequences, but allows for the possibility of frameshift mutations. We then used the aligned CDS sequences to identify where indels were present and found the corresponding genomic positions for every nucleotide in the alignment, inserting gaps where indels occurred. We then extracted the genomic positions for each SNP position as well as the genomic position for its aligned nucleotide. We retained only those homoeologous SNP positions in which both positions had a confidently called ancestral allele (described above) and in which the ancestral allele matched between the two homoeologs. Importantly, for homoeologs that were encoded in opposite orientations in the reference genome (i.e., one homoeolog was encoded on the forward strand of the reference genome, and the other homoeolog was encoded on the reverse complement), we ensured that the inferred ancestral states for the two SNP positions included both nucleotides of a purine/pyrimidine pair (e.g., the ancestral state for homoeologous SNP was "A," whereas the ancestral state of the other homoeologous SNP was "T"). We also removed any pair of homoeologous SNPs in which more than two alleles were present (while similarly treating homoeologous pairs encoded in opposite directions as described in the previous sentence).
In total, we only used those SNP sites that: 1) did not link to an indel in its homoeologous pair, 2) were biallelic and had consistently inferred ancestral states in the two subgenomes, 3) the derived allele was found in only one of the two subgenomes or their respective diploid progenitors, and 4) the derived allele was fixed in a diploid and segregating in its respective subgenome (or vice-versa).
Deleterious Mutations Accumulate Faster in Allopolyploids . doi:10.1093/molbev/msac024 MBE
Downloaded from https://academic.oup.com/mbe/article/39/2/msac024/6517786 by guest on 02 December 2022