A Tracing Model for the Evolutionary Equilibrium of Octoploids.

Testing Hardy-Weinberg equilibrium (HWE) is a fundamental approach for inferring population diversity and evolution, but its application to octoploids containing eight chromosome sets has not well been justified. We derive a mathematical model to trace how genotype frequencies transmit from parental to offspring generations in the natural populations of autooctoploids. We find that octoploids, including autooctolpoids undergoing double reduction, attach asymptotic HWE (aHWE) after 15 generations of random mating, in a contrast to diploids where one generation can assure exact equilibrium and, also, different from tetraploids that use 5 generations to reach aHWE. We develop a statistical procedure for testing aHWE in octoploids and apply it to analyze a real data set from octoploid switchgrass distributed in two ecologically different regions, demonstrating the usefulness of the test procedure. Our model provides a tool for studying the population genetic diversity of octoploids, inferring their evolutionary history, and identifying the ecological relationship of octoploid-genome structure with environmental adaptation.

Introduction

As an evolutionary force of the organism to buffer against environmental perturbations, the evolutionary mechanisms of polyploidy have been a long-standing subject of population and evolutionary genetic research (Bever and Felber 1992; Ramsey and Schemske 1998; Otto and Whitton 2000; Soltis et al., 2004; Fawcett et al., 2009). While a number of population genetic studies are focused on tetraploids (Haldane 1930; Moody et al., 1993; Butruille and Boiteux 2000; Arnold et al., 2012; Dufresne et al., 2014; Meirmans et al., 2018; Sun et al., 2020), there is little knowledge about population variation in polyploids at a higher ploidy level. Octoploids, whose cells have a chromosome number as many as eight times the basic haploid chromosome number, are an even less-explored group of polyploids (Edger et al., 2019). As compared to tetraploids, octoploids have wide genetic diversity that can better adapt to changing environment (Johnson and Vance-Borland, 2016; Grabowski et al., 2017). Also, octoploids may have a larger body or organ size than tetraploids, implying their more desirable opportunity to be used in breeding programs. For example, basin wildrye have larger leaves, longer culms, and greater crown circumference for octoploids than tetraploids, although the numerical ranges of plant traits and their source climates overlap between ploidy types (Johnson and Vance-Borland, 2016). By creating a ploidy series from tetraploids to octoploid for althea (Hibiscus syriacus), Lattier et al. (2019) found an increase in guard cell length and rDNA signals with the level of ploidy.

Despite their evolutionary and economic value, it is unclear how octoploids vary and evolve across time and space scales, largely owing to the lack of suitable population genetic analysis tools. Hardy-Weinberg equilibrium (HWE) test has been widely used as an approach for inferring the evolutionary processes of natural populations (Waples 2015), but this approach was established on diploids, leaving its use to polyploids a mystery. Although the diploid-driven HWE test procedure has been modified to accommodate tetraploids (Moody et al., 1993; Meirmans et al., 2018; Sun et al., 2020), the use of such a modification to evaluate population variation in inherently more complicated octoploids is not justified. For example, double reduction has little impact on the asymptotic process of HWE for autotetraploids, but we do not know if this is true for autooctoploids that produce three types of gametes, characterized by two, one, or no double reduction, respectively. Sun et al. (2020) proposed a gamete-based approach for testing autotetraploid aHWE by estimating three equilibrium diploid gamete frequencies. This approach has a power to test the equilibrium of dosage-unknown markers because three gamete frequencies just can be estimated by three observable types of genotypes. However, this approach has no sufficient degrees of freedom to test aHWE for dosage-unknown markers in autooctoploids that produce five tetraploid gametes but still have three observable types of genotypes.

In this article, we propose an approach for HWE testing in octoploids using any type of molecular markers. We derive a system of recursive equations that transmit individual octoploid genotype frequencies from the parental to offspring generation under random mating. We find that, as opposed to diploids in which exact equilibrium can be reached after one generation of random mating, both allooctoploids and autooctoploids can only gradually approach aHWE after 15 generations, different from tetraploids that use 5 generations to reach aHWE. We propose specific statistical procedures for testing aHWE using dosage-known and dosage-unknown markers in octoploids. We investigate the power of octoploid HWE detection through computer simulation. By analyzing population genetic data of allooctoploid switchgrass (Grabowski et al., 2017), we validate the usefulness of our testing procedure.

Mathematical Model

Our analysis is based on the segregation of biallelic single nucleotide polymorphisms (SNPs). At such a SNP, there are a total of nine genotypes, AAAAAAAA (8A), AAAAAAAa (7A1a), AAAAAAaa (6A2a), AAAAAaaa (5A3a), AAAAaaaa (4A4a), AAAaaaaa (3A5a), AAaaaaaa (2A6a), Aaaaaaaa (1A7a), and aaaaaaaa (8a) in an octoploid population. These genotypes will produce five types of tetraploid gametes during meiosis, AAAA, AAAa AAaa, Aaaa, and aaaa, with different frequencies determined by Mendel’s first law and the rate of double reduction (α) defined as the probability of two sister chromatids occurring in the same gamete (Darlington 1929; Mather 1935; Haynes and Douches 1993). For 8A or 8a, the same gamete type is identified although its formation results from either double reduction or non-double reduction. Table 1 tabulates the frequencies of tetraploid gametes produced by each octoploid genotype. We use AAAAAAAa (written as A 1 A 2 A 3 A 4 A 5 A 6 A 7 a for the identification of allele A) as an example to demonstrate how we derive these gamete frequencies. This genotype produces 8 diploids A 1 A 1, A 2 A 2, A 3 A 3, A 4 A 4, A 5 A 5, A 6 A 6, A 7 A 7, aa, with the total frequency of α, through double reduction, and 28 diploids A 1 A 2, A 1 A 3, A 1 A 4, A 1 A 5, A 1 A 6, A 1 A 7, A 1 a, A 2 A 3, A 2 A 4, A 2 A 5, A 2 A 6, A 2 A 7, A 2 a, A 3 A 4, A 3 A 5, A 3 A 6, A 3 A 7, A 3 a, A 4 A 5, A 4 A 6, A 4 A 7, A 4 a, A 5 A 6, A 5 A 7, A 5 a, A 6 A 7, A 6 a, A 7 a, with the total frequency of 1–α, without double reduction. Together, we have the frequencies of AA as (1/8)α×7 + (1/28)(1–α)×21 = 3/4 + 1/8α, the frequencies of Aa as (1/8)α×0 + (1/28)(1–α)×7 = 1/4 – 1/4α, and the frequencies of aa as (1/8)α×1 + (1/28)(1–α)×0 = 1/8α. Tetraploid gametes produced by genotype AAAAAAAa include AAAA, with frequency calculated as (3/4 + 1/8α)2 = 9/16 + 1/64α2+3/16α, AAAa, with frequency as 2(3/4 + 1/8α) (1/4 – 1/4α) = 3/8−1/16α2 −5/16α, AAaa, with frequency as 2(3/4 + 1/8α)(1/8α) + (1/4 – 1/4α)2 = 1/16 + 3/32α2 +1/16α, and Aaaa, with frequencies as 2(1/4 – 1/4α)(1/8α) = (−1/16)α2+1/16α, and aaaa, with frequency as (1/8α)2 = 1/64α 2. All other gamete frequencies can be derived in a similar way.

TABLE 1

Gamete frequencies generated by different autooctoploid genotypes.

Gamete
Genotype	AAAA	AAAa	AAaa	Aaaa	aaaa
AAAAAAAA	1	0	0	0	0
AAAAAAAa	9/16 + 1/64α2+3/16α	3/8−1/16α2	1/16 + 3/32α2 +1/16α	(−1/16)α2	1/64α2
		−5/16α		+1/16α
AAAAAAaa	225/784 + 9/196α2 +45/196α	45/98−9/49α2	87/392 + 27/98α2	3/98−9/49α2	1/784 + 9/196α2+3/196α
		−27/98α	−6/49α	+15/98α
AAAAAaaa	25/196 + 225/3136α2 +75/392α	75/196	285/784	45/392	9/784
		−225/784α2	+675/1568α2	−225/784α2 +135/784α	+225/3136α2+45/784α
		−75/784α	−255/784α
AAAAaaaa	9/196 + 4/49α2	12/49−16/49α2 +4/49α	41/98 + 24/49α2	12/49−16/49α2	9/196 + 4/49α2+6/49α
	+6/49α		−20/49α	+4/49α
AAAaaaaa	9/784 + 225/3136α2 +45/784α	45/392−225/784α2 + 135/784α	285/784+	75/196	25/196
			675/1568α2	−225/784α2−75/784α	+225/3136α2+75/392α
			−255/784α
AAaaaaaa	1/784 + 9/196α2 +3/196α	3/98−9/49α2 + 15/98α	87/392 + 27/98α2	45/98−9/49α2	225/784
			−6/49α	−27/98α	+9/196α2
					+45/196α
Aaaaaaaa	1/64α2	(−1/16)α2	1/16 + 3/32α2	3/8−1/16α2	9/16
		+1/16α	+1/16α	−5/16α	+1/64α2 + 3/16α
aaaaaaaa	0	0	0	0	1

The newly-produced gametes combine randomly to generate offspring genotypes. We use P (t) (j = 1 for 8A, 2 for 7A1a, 3 for 6A2a, 4 for 5A3a, 5 for 4A4a, 6 for 3A5a, 7 for 2A6a, 8 for 1A7a, and 9 for 8a) to denote the frequencies of nine octoploid genotypes in the tth (parental) generation. Random mating of these parental genotypes produces 45 possible combination types, each of which forms different offspring genotypes, with frequencies depending on Table 1’s gamete frequencies and the frequencies of 45 parental combinations (Supplementary Table S1). We derive the mathematical expression of each genotype frequency in the (t+1)th (offspring) generation derived from those of the tth (parental) generation after random mating. Let denote the frequency of offspring genotype derived from parental mating type ( ). Then, the total frequency of offspring genotype j is the sum of these offspring genotype frequencies, weighted by mating frequencies, expressed as This explicit expression of Eq. 1 is a group of recursive equations, indicating that the frequency of genotype j in the (t+1)th generation is jointly determined by the mating frequencies of relevant genotypes in the tth generation, Mendelian segregation, and double reduction rate. As can be seen from Table S1, P (t+1) has a complicated but explicit relationship with P (t).

To numerically explore how P (t) transmits to P (t+1), we randomly sample nine parental genotype frequencies from their space, substitute these sampled values under different double reduction rates into the expression of offspring genotype frequencies given in Eq. 1, use these estimated offspring genotype frequencies to calculate the subsequent offspring genotype frequencies, and repeat this process until 20 generations have passed. By plotting genotype frequencies against generation (Figure 1A), we find that genotype frequencies will not reach absolutely stable values even after many generations of random mating. However, genotype frequencies will quickly approach the stationarity after 15 generations. We repeat the above procedure by sampling 1,000 sets of initial genotype frequencies, from which a similar conclusion is reached, i.e., octoploid genotype frequencies do not attach exact HWE but aHWE over generation, a similar phenomenon detected in tetraploids. Generation 15 can be regarded as one after which the octoploid population is at aHWE.

FIGURE 1

Generation-dependent change of genotype frequencies under different levels of double reduction (α) in a panmictic octoploid population. (A) Change trends of individual genotype frequencies, initiated with (0.10, 0.10, 0.15, 0.10, 0.20, 0.10, 0.05, 0.10, 0.10). (B) Change trends of homozygote and heterozygote genotype frequencies derived from (A). (C) Change trends of homozygote and heterozygote genotype frequencies, initiate with extremely high homozygote frequencies (0.99) and extremely low heterozygote frequencies (0.01). (D) Equilibrium genotype frequencies of homozygote and heterozygote change as a function of α. Legends of α are indicated by color metrics.

We plot the overall frequencies of genotypes 8A and 8a and those of genotypes 7A1a, 6A2a, 5A3a, 4A4a, 3A5a, 2A6a, and 1A7a over generation (Figure 1B), from which we find that homozygotes and heterozygotes change their frequencies differently during the approximation of aHWE. While homozygote frequency tends to decrease with generation, heterozygote frequency displays a generation-dependent increase. In an extreme example, we start the homozygote frequency of near one and the heterozygote frequency of near zero as initial values, in which the same trends are observed (Figure 1C). This phenomenon suggests that by approaching aHWE, octoploids are equipped to increase genetic diversity in natural populations.

We also investigate how double reduction influences the approximation of aHWE in octoploids. Under a range of double reduction rates, we find that octoploids are always close to stabilize their genotype frequencies in 15 generations, suggesting that double reduction is neutral for aHWE (Figure 1A). The case of no double reduction in the autopolyploid model reduces to allopolyploids if no preferential pairing is assumed. Thus, it is postulated that both autopolyploids and allopolyploids follow the same rule of aHWE. In autooctoploids, double reduction does influence the distribution of equilibrium genotype frequencies to some extent; the frequency of double reduction is positively associated with equilibrium frequencies of homozygotes but negatively associated with the equilibrium frequencies of heterozygotes (Figure 1D). Thus, double reduction may be regarded as a determinant of genetic diversity in autopolyploids.

Statistical Testing of aHWE

Recursive Test

This test was proposed by Sun et al. (2020) to test aHWE in tetraploids. It is straightforward to extend it to test octoploid aHWE. As shown from recursive Eq. 1, an octoploid population reaches asymptotic equilibrium at approximately generation 15 (Figure 1A). Thus, the genotype frequencies at generation 15 are regarded as a proxy of equilibrium frequencies, denoted as for a SNP. Let N (totaling to N) and P denote the observations and frequencies of nine genotypes in the current octoploid natural population. We formulate the likelihood of genotypic observations as where c is the constant and the maximum likelihood estimates (MLEs) of P can be solved as . We test whether equilibrium genotype frequencies are not different from observed genotype frequencies by using Pearson’s chi-square testing approach, expressed as which is thought to follow a chi-square distribution with eight degrees of freedom. If test statistics is larger than the critical value , then we claim that the population deviates from aHWE at the SNP considered. Otherwise, the population is at equilibrium for this locus. Multiple comparisons will be needed if a set of SNPs are tested in a population genetic study.

Gamete-Based Model

Under HWE, genotype frequencies are expressed as the frequency products of the gametes that form the genotypes. Let P , P , P , P , and P denote the frequencies of five tetraploid gametes produced by octoploid genotypes in a natural population. Equilibrium genotype frequencies (Q ) are expressed as

We formulate a likelihood under HWE, i.e., which includes multiple mixture terms with component proportions determined by gamete frequencies. We implement the derivative-free EM algorithm to obtain the MLEs of gamete frequencies ( ) under HWE. Plugging the values into (4) obtains the likelihood value under the null hypothesis that the population is at HWE. The likelihood under the alternative hypothesis, i.e., the population deviates from HWE, is calculated by Eq. 2. The log-likelihood ratio (LR) is then calculated as which follows a chi-square distribution with four degrees of freedom. By comparing the LR test statistics with critical value , we determine whether the SNP considered is at HWE in the octoploid population.

Allele-Based Model

Because of technical and economic reasons, many studies may genotype octoploids at a limited resolution, in which case seven heterozygotes 7A1a, 6A2a, 5A3a, 4A4a, 5A3a, 6A2a, and 7A1a cannot be distinguished from each other. Thus, for such dosage-unknown markers, there are only three distinguishable genotypes, i.e., two homozygotes each for an alternative allele and one mixed heterozygote. The frequencies and sizes of two homozygotes and the heterozygote at a dosage-unknown SNP are denoted as P 8A, P 8a, and P_ and N 8A, N 8a, and N_, respectively. It is impossible that five gamete frequencies are estimated from three distinguishable genotypes. To make HWE testable for dosage-unknown markers, we make an assumption; i.e., alleles randomly unite to generate gametes during meiosis. Let p and q denote allele frequencies of A and a, respectively. Under HWE, we express genotype frequencies of homozygotes and heterozygote in terms of allele frequencies, i.e.,

Under the null hypothesis of HWE, the likelihood is formulated as in a comparison with the likelihood under the alternative hypothesis of no HWE,

Likelihood Eq. 8 contains a mixture term, whose solution can be made by implementing the derivative-free EM algorithm. We calculate the LR from likelihoods Eqs 8, 9 as a test statistic to test whether the dosage-unknown SNP considered is at HWE. This test statistic is thought of as being following a chi-square distribution with one degree of freedom.

Results

Example

In a genome-wide association study of allopolyploid switchgrass, Grabowski et al. (2017) collected samples from the southern-adapted upland ecotype and northern-adapted upland ecotype of this species in the United States. The sampled plants include allotetraploids from both ecotypes and allooctoploids from the upland ecotype. Here, we test whether upland octoploids, distributed in the east and west regions, deviate from aHWE. East and west populations include 66 and 101 samples, for which 24,859 and 23,795 quality SNPs are, respectively, available for our analysis. All markers are dosage-unknown, at each of which there are three distinguishable genotypes. Thus, the allele-based model is used to test whether they are at aHWE in the populations.

Figure 2 illustrates the significance test of deviation from aHWE for dosage-unknown SNPs distributed throughout the switchgrass genome in the east and west populations of the upland ecotype. The majority of SNPs are detected to deviate from aHWE, with the western population having a slightly larger proportion (93.3%) than the eastern population (90.7%). This result suggests that the segregating genes of allootoploid switchgrass plants, especially those from the western population, throughout the entire genome are on their way toward equilibrium. In a similar aHWE test for allotetraploids of the same species, only a small proportion of SNPs was detected to deviate from aHWE (Sun et al., 2020). This comparison suggests that while lower-ploidy switchgrass tends to be evolutionarily stable, higher-ploidy switchgrass is still experiencing a strong evolutionary change.

FIGURE 2

Manhattan plots of significance test for marker aHWE throughout the genome in an allooctoploid switchgrass upland ecotype collected from eastern (A) and western populations (B). SNPs with unknown chromosomes are given in the “unannotated” part. Horizontal line denotes the significance level after Bonferroni correction.

Power Analysis

We perform simulation studies to examine the power of HWE detection. Using given gamete frequencies for dosage-known markers and allele frequencies for dosage-unknown markers, we simulate genotype data by allowing simulated genotype frequencies to deviate from equilibrium genotype frequencies by 20%. Under this deviation, genotype data are simulated for different sample sizes (N = 50, 100, 200, 400). We use the gamete- and allele-based models to test HWE for dosage-known and dosage-unknown markers, respectively. This simulation and estimation process is repeated 1,000 times to empirically calculate HWE detection power, defined as the proportion of the simulations for which LR tests are significant.

We find that to detect aHWE for octoploid populations, a sample size of n = 200 is required, under which the power of nonequilibrium detection reaches 0.94 (Table 2). A sample size of n = 100 can only have about 0.60 power, whereas n = 50 fails to detect deviation from aHWE in most cases. We also analyze the false positive rate of the model by simulating the data with genotype frequencies deviating from equilibrium genotype frequencies by zero. In all cases, the model has a low type I error (Table 2).

TABLE 2

Power analysis and false positive rate of aHWE detection under different sample size.

Degree of deviation	50	100	200	400
0	0.03	0.02	0.01	0.01
20	0.37	0.58	0.94	1.00

Note: Degree of deviation describes how much genotypes frequencies deviate from equilibrium frequencies.

Discussion

It has been recognized that polyploids gradually reach HWE through random mating, but it is not very clear how many generations they mate to approach equilibrium (Geiringer 1949; Bever and Felber 1992). Through a mathematical derivation, we showed that tetraploids never attain absolute HWE but aHWE after 5 generations of random mating (Sun et al., 2020). This phenomenon, in contrast to diploids in which exact HWE can be achieved only after one generation of random mating, suggests that polyploids are a “perpetual machine” of evolution; i.e., they can evolve even without actions of evolutionary forces, such as natural selection, mutation, genetic drift, admixture and so on.

In this article, we show that octoploids require more generations (i.e., 15) to approach aHWE. As compared to low ploidy-level polyploids, high ploidy-level polyploids have more allelic combinations and, thus, larger genetic diversity (Johnson and Vance-Borland, 2016; Grabowski et al., 2017). More generations required to reach aHWE suggest that high ploidy-level polyploids can maintain a longer-standing time of evolution than low ploidy-level polyploids. This theoretical postulation is well in agreement with empirical observations for many plants whose high-ploidy polyploids have greater diversity than low-ploidy relatives in the same habitat (Johnson and Vance-Borland, 2016). Double reduction is an important phenomenon in autopolyploids. It has been thought to affect evolutionary processes due to genetic drift (Moody et al., 1993) and gametophytic selection (Butruille and Boiteux 2000). Double reduction can change population structure and diversity by increasing the frequencies of homozygotes, but it does not strikingly accelerate the evolution of autopolyploid populations because of its subtle impact on deviation from HWE.

Based on mathematical derivations of frequency transmission, we propose three approaches for testing aHWE in octoploids using dosage-known and dosage-unknown markers. By analyzing dosage-unknown marker data of octoploid switchgrass collected from its natural distribution (Grabowski et al., 2017), we validate the usefulness of the gamete-based equilibrium-detecting model. Computer simulation is performed to determine the sample size required to detect aHWE, as guidance for designing octoploid population and evolutionary genetic studies. An increasing number of studies have begun to investigate the ecological relationship between HWE deviation and the organisms’ adaptation to heterogeneous environments (Tyukmaeva et al., 2011; Tarazona et al., 2019; Carvalho et al., 2021). Our aHWE testing model provide a generic tool to infer the genetic variation and evolutionary status of octoploid populations.