Algebraic equivalence and sampling properties of F-statistics in the genomic era

Abstract

This note evaluates variance- and heterozygosity-based genomic estimates of F-statistics through algebraic derivation and simulation. Genetic variance was partitioned among populations, among individuals within populations, and within individuals. The method-of-moments estimator of F_IS was evaluated in 100 Monte Carlo replicates of 1,000 independent SNPs incorporating population structure, sample-size imbalance, and missing genotypes. In the principal scenario (F_IS = 0.25; F_ST = 0.38), ANOVA, Weir-Cockerham, allele-sharing, and corrected Nei estimates of F_IS were 0.2498, 0.2498, 0.2497, and 0.2502, respectively, whereas the pooled statistic was 0.4844. Across the simulations, ANOVA and corrected Nei F_IS remained near 0.25 as FST increased from 0 to 0.50, while the pooled measure increased from 0.2481 to 0.5667. The mean within-population FIS was not much affected by missing data (up to 30%) and sample-size imbalances. The formulations therefore agree when they target the same hierarchical parameter and reference set.

.

Keywords:
F-statistics; ANOVA; variance components; population structure; reference population; Wahlund effect; genomic data

INTRODUCTION

The quantification of genetic structures, inbreeding, and population differentiation remains a foundational objective in quantitative genetics and plant breeding (Wright 1951). Two primary analytical paradigms have coexisted while estimating parameters using codominant molecular data since Wright's fixation indices ($F_{ST}$, $F_{IT}$, $F_{IS}$) were conceptualized (Wright 1951, Cockerham 1969): the allele/genotype-counting framework, exemplified by Nei's gene diversity statistics (Nei 1973, Nei 1978), and the variance-component framework, formalized through the analysis of variance (ANOVA) by Cockerham (1969, 1973) and expanded by Weir and Cockerham (1984).

At the true parametric level, variance-component and heterozygosity formulations can be reconciled when they refer to the same hierarchical level and reference population (Cockerham 1973, Long 1986). Empirical analyses, however, use finite numbers of populations and individuals, unequal sampling, and locus-specific missing data (Weir and Cockerham 1984). Finite-sample corrections for gene diversity and fixation indices have long been available in the heterozygosity framework (Nei and Chesser 1983, Nei 1986, Nei 1987), so heterozygosity-based estimators should not be treated as intrinsically uncorrected counting procedures. Cockerham and Weir (1986) further showed that apparent disagreements can arise because alternative formulations condition on population structure differently, including whether sampled populations are treated as fixed extant groups or as realizations from an evolutionary population process. More recent work has made the reference-set issue explicit: measures of inbreeding, relatedness, and population structure are relative to a specified set of allele pairs, and estimates can change when the reference population or weighting scheme changes (Weir and Goudet 2017, Goudet and Weir 2023, Weir and Goudet 2026). In plant breeding, this distinction is relevant whenever germplasm collections, breeding populations, landraces, or provenances contain differentiated groups that may be pooled during genomic analysis.

The algebraic transition from a discrete, binary allelic state to Mean Squares in an ANOVA is frequently treated as a black box in pedagogical and applied contexts. This technical note provides a step-by-step derivation of variance-based F-statistics from a binary allelic indicator and relates the resulting components to heterozygosity-based formulations. Variance-component estimators for F-statistics are well established (Cockerham 1973, Weir and Cockerham 1984), while revised heterozygosity and allele-sharing approaches provide alternative moment estimators when the reference set is defined appropriately (Nei 1986, Weir and Goudet 2017, Goudet and Weir 2023). The specific contribution tested here is therefore not a claim that ANOVA is universally superior to heterozygosity-based estimators. Instead, we ask whether formulations that target the same within-population FIS recover the same parameter, and how strongly a statistic based on pooled allele frequencies can diverge when differentiated subpopulations are treated as a single reference population. This distinction also provides a direct diagnostic for separating estimator behavior from a change in the biological quantity being estimated.

Theoretical framework and the linear model

Let $Y_{ijk}$ be a discrete binomial indicator variable denoting the presence or absence of a target allele at a specific locus, for the $k$-th allelic copy $\left(\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{}k=1,2\right)$ for diploid organisms) of the $j$-th individual, sampled within the $i$-th subpopulation (Cockerham 1969, Weir 1996):

$Y_{ijk}=\left\{1,ifallele1ispresent;0,otherwise\right\}$.

For an allelic copy sampled from the reference population, $Y_{ijk}$ remains a binary variable and its expectation equals the allele probability $p$: $E\left(Y_{ijk}\right)=p$. To model hierarchical genetic structures, the random variable is partitioned using an additive random-effects linear model (Cockerham 1973, Weir and Cockerham 1984): $Y_{ijk}=p+u_i+v_{ij}+e_{ijk}$, where $p$ represents the grand mean allele frequency across the metapopulation. The term $u_i$ is the random effect of the $i$-th subpopulation, $v_{ij}$ is the random effect of the $j$-th individual nested within the $i$-th subpopulation, and $e_{ijk}$ represents deviations between allelic copies within the same individual.

Under Cockerham's genetic design, the random effects have expectation zero and variances $Var\left(u_i\right)={\sigma }_b^2$, $Var\left(v_{ij}\right)={\sigma }_d^2$, and $Var\left(e_{ijk}\right)={\sigma }_w^2$, representing variation among populations, among individuals within populations, and within individuals, respectively. With $p$ defined for the metapopulation, the hierarchical components are ${\sigma }_b^2=p\left(1-p\right)F_{ST}$, ${\sigma }_d^2=p\left(1-p\right)\left(1-F_{ST}\right)F_{IS}$, and ${\sigma }_w^2=p\left(1-p\right)\left(1-F_{ST}\right)\left(1-F_{IS}\right)$. Therefore, ${\sigma }_b^2+{\sigma }_d^2+{\sigma }_w^2=p\left(1-p\right)$, and $1-F_{IT}=\left(1-F_{IS}\right)\left(1-F_{ST}\right)$. When $F_{ST}=0$, as in a single unstructured population, the hierarchy reduces to ${\sigma }_d^2=p\left(1-p\right)F_{IS}$ and ${\sigma }_w^2=p\left(1-p\right)\left(1-F_{IS}\right)$, so $V_Y=p\left(1-p\right)$.

Step-by-step algebraic derivation of the ${\mathit{F}}_{\mathit{I}\mathit{S}}$ estimate

Reduction of the binomial total variance

To establish the Bernoulli variance, note that $Y_{ijk}\in \left\{0,1\right\}$, so $Y_{ijk}^2=Y_{ijk}$ and $E\left(Y_{ijk}^2\right)=E\left(Y_{ijk}\right)=p$. Therefore, $Var\left(Y\right)=E\left(Y^2\right)-{\left[E\left(Y\right)\right]}^2=p-p^2=p\left(1-p\right)$. Replacing $p$ by the observed allele frequency $\hat{p}$ gives the plug-in Bernoulli variance $\hat{p}\left(1-\hat{p}\right)$; this plug-in expression should not be described as the unbiased finite-sample variance estimator.

Expected mean squares and parameter isolation

In empirical designs utilizing finite datasets, parameters are estimated from Mean Squares (MS) in the hierarchical ANOVA. For a balanced diploid design with $g$ subpopulations and $n$ individuals per subpopulation, let $MSP$ denote the mean square among populations, $MSI$ the mean square among individuals within populations, and $MSG$ the mean square among allelic copies within individuals. The expected mean squares are $E\left(MSP\right)={\sigma }_w^2+2{\sigma }_d^2+2n{\sigma }_b^2$, $E\left(MSI\right)={\sigma }_w^2+2{\sigma }_d^2$, and $E\left(MSG\right)={\sigma }_w^2$. Thus, the method-of-moments estimates are ${\hat{\sigma }}_b^2=\left(MSP-MSI\right)/\left(2n\right)$, ${\hat{\sigma }}_d^2=\left(MSI-MSG\right)/2$, and ${\hat{\sigma }}_w^2=MSG$.

Derivation of the Fixation Index ${\mathit{F}}_{\mathit{I}\mathit{S}}$

By definition, the within-population fixation index is $F_{IS}={\sigma }_d^2/\left({\sigma }_d^2+{\sigma }_w^2\right)$. Substituting the method-of-moments components gives ${\hat{F}}_{IS}=\left[\left(MSI-MSG\right)/2\right]/\left\{\left[\left(MSI-MSG\right)/2\right]+MSG\right\}$, which simplifies to ${\hat{F}}_{IS}=\left(MSI-MSG\right)/\left(MSI+MSG\right)$. This is the corresponding method-of-moments estimator for the balanced diploid hierarchy.

Mathematical reconciliation with Wright's heterozygosities

To confirm that this variance-component partition is consistent with Wright's heterozygosity formulation, consider the within-population reference level. For a subpopulation with local allele frequency $p_i$, $H_E=2p_i\left(1-p_i\right)$, and $H_O$ is the observed heterozygote frequency. The within-population variance components map to ${\sigma }_w^2=H_O/2$ and ${\sigma }_d^2=\left(H_E-H_O\right)/2$. Substitution gives $F_{IS}={\sigma }_d^2/\left({\sigma }_d^2+{\sigma }_w^2\right)=\left(H_E-H_O\right)/H_E=1-H_O/H_E$. Thus, the ANOVA and heterozygosity formulations agree when $H_E$ and $H_O$ are defined at the same within-population reference level.

MATERIAL AND METHODS

Genomic extension and simulation design

High-throughput genotyping datasets, including RAD-seq and GBS, can contain locus-specific missingness, allelic dropout, and genotyping error (Attard et al. 2018). Unequal population sample sizes also require explicit treatment in a hierarchical analysis. For the among-population component, the corrected sample-size coefficient nc follows the standard variance-component framework (Cockerham 1973, Weir and Cockerham 1984). This coefficient pertains to estimation of the among-population component and should not be inserted into the within-population derivation of FIS above. Simulations were performed in R 4.5.1 (R Core Team 2025), and independent software checks used hierfstat 0.5-11 (Goudet 2005). In the simulations considered here, missing genotypes were generated completely at random; genotyping error and allelic dropout were not simulated.

Simulation parameters

A genomic architecture comprising 1,000 independent bi-allelic SNP loci scored across four subpopulations was simulated. Population-specific allele frequencies were drawn independently from a symmetric Beta distribution with $\alpha =\left(1-F_{ST}\right)/\left(2F_{ST}\right)$, giving $E\left(p_i\right)=0.5$ and target $F_{ST}=0.38$ in the principal scenario. Within each subpopulation, genotype probabilities were $P\left(AA\right)=p_i^2+F_{IS}{p}_i\left(1-p_i\right)$, $P\left(Aa\right)=2p_i\left(1-p_i\right)\left(1-F_{IS}\right)$, and $P\left(aa\right)={\left(1-p_i\right)}^2+F_{IS}{p}_i\left(1-p_i\right)$, with target within-population $F_{IS}=0.25$. The principal sampling design was $n=\left[100,80,70,50\right]$ ($N=300$), and 15% of genotypes were masked completely at random. For each locus, genotypes were expanded into two binary allelic indicators. The hierarchical ANOVA estimated $F_{IS}$ from the mean square among individuals within populations (MSI) and the mean square among allelic copies within individuals (MSG) as $F_{IS}=\left(MSI-MSG\right)/\left(MSI+MSG\right)$. A corrected heterozygosity estimator following Nei (1987) was also calculated within populations. A separate pooled diagnostic used the allele frequency and observed heterozygosity after combining all non-missing individuals across subpopulations; because this changes the reference level, it was not treated as a second estimator of within-population $F_{IS}$.

Sensitivity analyses and independent validation

For each sensitivity scenario, 100 independent Monte Carlo replicates were generated, with 1,000 SNPs per replicate. Population structure was varied across $F_{ST}=0$, $0.05$, $0.10$, $0.20$, $0.38$, and $0.50$ while keeping $F_{IS}=0.25$, the principal unbalanced sampling design, and 15% missing genotypes. Missingness was then varied from 0 to 0.30 under balanced sampling for both $F_{ST}=0$ and $F_{ST}=0.38$. Sample-size imbalance was evaluated with $n=\left[75,75,75,75\right]$, $\left[90,80,70,60\right]$, $\left[100,80,70,50\right]$, and $\left[150,75,50,25\right]$, with no missing genotypes and $F_{ST}=0$ or $0.38$. Calibration was assessed at $F_{IS}=0$, $0.10$, $0.25$, and $0.50$ under $F_{ST}=0$ and $0.38$. Multilocus ANOVA estimates were calculated as the ratio of summed locus-specific numerators and denominators. Across Monte Carlo replicates we recorded the mean, standard deviation, empirical 2.5% and 97.5% quantiles, bias, and RMSE. The principal simulated dataset was checked with hierfstat::wc, an independent software implementation of the Weir-Cockerham variance-component estimator (Weir and Cockerham 1984); hierfstat::basic.stats, which reports Nei's diversity and fixation statistics (Nei 1987); and hierfstat::fs.dosage, which implements reference-set allele-sharing $F$-statistics for dosage data (Weir and Goudet 2017; Goudet and Weir 2023). The comparison after pooling all individuals was retained only as a diagnostic of the effect of changing the reference population and was not treated as a second estimator of within-population $F_{IS}$.

RESULTS AND DISCUSSION

Cross-method validation of within-population ${\mathit{F}}_{\mathit{I}\mathit{S}}$

In the principal scenario, the ANOVA multilocus estimate was $F_{IS}=0.249776$ for a target of 0.25 (Table 1). The hierfstat::wc implementation returned the identical value, $F_{IS}=0.249776$. This exact agreement is an important software-level verification of our implementation, but it should not be interpreted as validation by an unrelated statistical theory because both calculations are based on the Weir-Cockerham variance-component framework. More methodologically distinct estimators converged on the same target: fs.dosage returned $F_{IS}=0.2497$ using the reference-set allele-sharing formulation, and Nei's corrected heterozygosity statistics returned $F_{IS}=0.2502$. This concordance is consistent with the theoretical links among Cockerham's intraclass-correlation formulation, revised Nei estimators, and allele-sharing statistics when the same hierarchical parameter and reference set are used (Cockerham 1973, Nei 1986, Weir and Goudet 2017, Weir and Goudet 2026).

Table 1
Cross-method and independent software validation of fixation indices in the principal simulation scenario (1,000 SNPs; target $F_{IS}=0.25$; target $F_{ST}=0.38$; $n=\left[100,80,70,50\right]$; 15% MCAR missing genotypes). The pooled statistic is shown only as a reference-level diagnostic and is not an alternative estimator of within-population $F_{IS}$.

Framework	FIS	FST / FST'	Interpretation
Parametric target	0.2500	0.3800	Generating values
ANOVA, this study	0.249776	-	Within-population variance components
Weir & Cockerham 1984 (wc)	0.249776	0.383695	Independent software implementation of the same variance-component framework
Allele sharing (fs.dosage)	0.2497	0.3824	Reference-set allele-sharing estimator; distinct moment formulation
Nei 1987 (basic.stats)	0.2502	0.3167; 0.3819 (FST')	Heterozygosity formulation; FST and corrected FST' are distinct outputs
Pooled global diagnostic	0.484446	-	Total-sample heterozygote deficit; different reference level, not within-population FIS

Population differentiation was also recovered closely by Weir-Cockerham ($F_{ST}=0.3837$), allele sharing ($F_{ST}=0.3824$), and the corrected Nei $F_{ST}\text{'}$ (0.3819) for a target $F_{ST}=0.38$. The pooled global statistic, by contrast, was 0.484446. Because this quantity compares observed heterozygosity with expected heterozygosity calculated from an allele frequency estimated after pooling the four differentiated subpopulations, it is not an alternative estimate of the same within-population $F_{IS}$. Population structure drives the pooled-reference shift Across 100 Monte Carlo replicates, the within-population estimators remained centered near $F_{IS}=0.25$ as population structure increased (Figure 1).

Figure 1
Effect of population structure on within-population $F_{IS}$ and on the pooled global fixation measure. Points are means across 100 Monte Carlo replicates of 1,000 independent SNPs; error bars are $\pm 1$ SD. The target within-population $F_{IS}$ was 0.25 in all scenarios. Population sizes were $\left[100,80,70,50\right]$ and missingness was 15% MCAR.

The ANOVA means ranged from 0.2497 to 0.2503 across $F_{ST}=0-0.50$, and the corrected Nei means ranged from 0.2482 to 0.2501. In contrast, the pooled global measure increased monotonically from 0.2481 at $F_{ST}=0$ to 0.2762, 0.3050, 0.3653, 0.4824, and 0.5667 at $F_{ST}=0.05$, $0.10$, $0.20$, $0.38$, and $0.50$, respectively. At the principal $F_{ST}=0.38$, the Monte Carlo means were 0.2503 (SD 0.0026) for ANOVA, 0.2495 (SD 0.0028) for Nei, and 0.4824 (SD 0.0044) for the pooled statistic. This pattern is the expected consequence of changing the reference set rather than evidence that a properly defined heterozygosity estimator of $F_{IS}$ fails in structured data. Cockerham (1973) explicitly discussed the consequences of ignoring levels of population subdivision, and Weir and Goudet (2017) later formalized the general point that measures of allelic dependence are relative to a reference set. Goudet and Weir (2023) and Weir and Goudet (2026) extended this argument for genomic data, emphasizing that population-structure parameters should be interpreted relative to the populations used as the reference.

Our simulation provides a simple numerical illustration of that principle: the within-population target remains stable while the pooled-reference statistic tracks the amount of differentiation among populations. Missing data and sample-size imbalance Missing genotypes generated completely at random had little effect on the mean within-population $F_{IS}$ over the range tested (Figure 2). With $F_{ST}=0$, increasing missingness from 0 to 30% changed the ANOVA mean from 0.2497 to 0.2499 and the corrected Nei mean from 0.2497 to 0.2498. With $F_{ST}=0.38$, the corresponding values changed from 0.2504 to 0.2499 for ANOVA and from 0.2498 to 0.2490 for Nei; the pooled statistic remained near 0.485. Similarly, sample-size imbalance had little effect on within-population $F_{IS}$. Under $F_{ST}=0.38$, ANOVA means were 0.2500, 0.2500, 0.2501, and 0.2496 for the balanced, mild, current, and severe designs, while corrected Nei means were 0.2493, 0.2494, 0.2494, and 0.2486. These results are compatible with the purpose of finite-sample and unequal-sample corrections developed in the classical moment-estimation literature (Cockerham 1973, Nei and Chesser 1983, Weir and Cockerham 1984). They do not imply, however, that missing genomic data are generally harmless. Attard et al. (2018), studying genomic relatedness rather than $F_{IS}$, showed that missing data, allelic dropout, and genotyping error can affect estimator performance. Our conclusion is narrower: under the MCAR mechanism used here, missingness up to 30% produced little systematic displacement of the mean $F_{IS}$, whereas other error mechanisms were not tested. The pooled statistic was more sensitive to weighting under severe imbalance, declining from 0.4847 under equal sample sizes to 0.4633 for $n=\left[150,75,50,25\right]$, which is expected because changing sample weights changes the pooled allele frequency and therefore the reference itself. Calibration across $F_{IS}$ values Calibration across $F_{IS}=0$, $0.10$, $0.25$, and $0.50$ confirmed that both within-population estimators tracked the generating value over the tested range (Table 2). Across both $F_{ST}=0$ and $F_{ST}=0.38$, the absolute Monte Carlo bias of the ANOVA estimator never exceeded 0.0003, while the largest absolute bias of the corrected Nei estimator was 0.00124. This agreement is consistent with the finite-sample corrections developed by Nei and Chesser (1983) and with the equivalences among moment estimators discussed by Weir and Goudet (2026). When $F_{ST}=0$, the pooled statistic was also close to $F_{IS}$, with a small downward finite-sample deviation of approximately 0.0015-0.0022. Once population structure was introduced, however, the pooled means became 0.3097, 0.3788, 0.4821, and 0.6549 for generating $F_{IS}$ values of 0, 0.10, 0.25, and 0.50 at $F_{ST}=0.38$. The increasing separation is therefore not a calibration failure of either within-population estimator. It reflects the fact that the pooled statistic uses a total-sample allele frequency and consequently absorbs both within-population heterozygote deficit and among-population differentiation into a single global contrast.

Table 2
Monte Carlo calibration across target $F_{IS}$ values. Values are mean ± SD across 100 replicates of 1,000 independent SNPs. The pooled global measure is included only as a reference-level diagnostic.

Target FIS	Target FST	ANOVA FIS	Nei (1987) FIS	Pooled global measure
0.00	0.00	-0.0003 ± 0.0020	-0.0003 ± 0.0020	-0.0022 ± 0.0020
0.00	0.38	0.0003 ± 0.0021	0.0003 ± 0.0021	0.3097 ± 0.0054
0.10	0.00	0.1001 ± 0.0019	0.1001 ± 0.0021	0.0982 ± 0.0019
0.10	0.38	0.1000 ± 0.0025	0.0995 ± 0.0027	0.3788 ± 0.0050
0.25	0.00	0.2500 ± 0.0018	0.2501 ± 0.0019	0.2482 ± 0.0017
0.25	0.38	0.2501 ± 0.0026	0.2492 ± 0.0030	0.4821 ± 0.0042
0.50	0.00	0.4999 ± 0.0019	0.4999 ± 0.0020	0.4985 ± 0.0019
0.50	0.38	0.5000 ± 0.0024	0.4988 ± 0.0024	0.6549 ± 0.0031

Figure 2
Sensitivity to missing genotypes and sample-size imbalance. Points are means across 100 Monte Carlo replicates of 1,000 independent SNPs; error bars are $\pm 1$ SD. Top panels vary MCAR missingness under balanced sampling; bottom panels vary sample-size imbalance without missing genotypes. Left panels use $F_{ST}=0$ and right panels use $F_{ST}=0.38$. The target within-population $F_{IS}$ was 0.25.

Interpretation in relation to F-statistic theory

The revised analyses change the interpretation of the original comparison. Cockerham's variance-component framework and the estimators of Weir and Cockerham keep the hierarchy among populations, individuals within populations, and alleles within individuals explicit (Cockerham 1973, Weir and Cockerham 1984). Nei's framework provides heterozygosity-based fixation indices with finite-sample corrections (Nei and Chesser 1983, Nei 1986, Nei 1987). The present simulations show that these approaches are not in conflict for within-population $F_{IS}$ when they are aligned to the same reference: ANOVA, Weir-Cockerham, allele sharing, and corrected Nei estimates all recovered the generating value. This empirical result closely matches the modern synthesis of Weir and Goudet (2026), who showed that differences among several classical estimators are often more apparent than real once parameter definitions and reference sets are made explicit; under the assumptions of the Weir-Cockerham framework, their mean-square expressions coincide with allele-sharing formulations, and revised Nei estimators have corresponding allele-sharing forms. Thus, the correct contrast is not 'ANOVA versus counting', but rather estimators of the same hierarchical parameter versus a statistic whose reference population has been changed by pooling. The $F_{ST}$ results add a second caution. Weir-Cockerham $F_{ST}=0.3837$, allele-sharing $F_{ST}=0.3824$, and corrected Nei $F_{ST}\text{'}=0.3819$ were close to the target 0.38, whereas the uncorrected Nei $F_{ST}$ returned by basic.stats was 0.3167. We do not interpret this difference as evidence that the Nei framework is defective. Cockerham and Weir (1986) explicitly contrasted parameters defined for fixed extant populations with parameters that allow evolutionary variation among populations, and Nei (1986) likewise emphasized definition when comparing fixation indices.

The result reinforces a narrower requirement: the estimator, correction, and reference population should accompany the symbol $F_{ST}$ when values are compared across methods.

Scope and limitations

The simulations used four populations, a metapopulation mean allele frequency of 0.5, independent population-specific allele frequencies generated from a symmetric Beta distribution, independent SNPs, and missingness generated completely at random. These choices make the causal comparisons transparent but restrict the scope of inference. Goudet and Weir (2023) considered much more general population structures, including dependence among populations, whereas our simulations do not model migration-induced covariance, admixture, or shared population histories. Likewise, the present MCAR experiment does not cover genotype-dependent missingness, allelic dropout, or genotyping error, all of which can occur in reduced-representation genomic data (Attard et al. 2018). The SNPs were also independent. In real genomic data, linkage disequilibrium reduces the effective number of independent observations; Weir and Goudet (2026) therefore recommend resampling strategies such as block bootstrap when uncertainty is assessed across linked variants. Consequently, the numerical bias and precision reported here should be interpreted as properties of the simulated conditions, not as universal performance guarantees for genomic datasets.

Practical implications for plant breeding and genetic diversity research

For plant breeding and conservation genetics, the main practical requirement is to match the statistic to the biological level of interest. If $F_{IS}$ is the target, expected and observed heterozygosities should be defined within populations, or an equivalent hierarchical or reference-set estimator should be used. Pooling differentiated provenances, landraces, breeding groups, or conservation units can create a heterozygote deficit relative to the pooled allele frequency that should not be interpreted automatically as within-population inbreeding. This point is consistent with the long-standing hierarchical view of Cockerham (1973) and with the explicit reference-set perspective of Weir and Goudet (2017, 2026). Unequal sample sizes and missing genotypes still require transparent handling, but our simulations do not support universal thresholds for either factor. For genomic studies, the reference population, population weighting, missing-data mechanism, estimator, and multilocus aggregation rule should therefore be stated explicitly so that $F_{IS}$ and $F_{ST}$ values are biologically interpretable and comparable across studies.

CONCLUSIONS

The hierarchical variance-component formulation correctly separates variation among populations, among individuals within populations, and within individuals. In the revised simulations, the ANOVA estimator recovered within-population $F_{IS}$ with negligible bias. The exact agreement with hierfstat::wc verifies the implementation of the same variance-component framework, while the close agreement with the reference-set allele-sharing and corrected Nei estimators provides cross-method support that the formulations recover the same within-population parameter when the reference is aligned. The large value obtained after pooling differentiated populations was not a second estimate of the same $F_{IS}$; it arose because the pooled calculation used a broader reference population. Across the scenarios tested, increasing $F_{ST}$ strongly shifted the pooled statistic, whereas MCAR missingness up to 30% and the tested sample-size imbalances had little effect on the mean within-population $F_{IS}$. The revised evidence therefore does not support the original interpretation that heterozygosity-based estimation intrinsically produces a two-fold inflation under missing or unbalanced data. Instead, it supports the more limited and theoretically grounded conclusion that $F$-statistics must be defined, estimated, and interpreted at the correct hierarchical level and relative to an explicit reference set. Further work with linked variants, alternative missing-data mechanisms, dependent population histories, and empirical plant genomic datasets is required before broader quantitative generalizations are made.

ACKNOWLEDGEMENTS

Evandro V. Tambarussi is supported by the CNPq research fellowship grant no. 303789/2022-0.

Data Availability Statement

The R script generated and the simulation outputs used in this research are available at https://doi.org/10.48331/SCIELODATA.UF7OPH

REFERENCES

• Attard CRM, Beheregaray LB, Möller LM2018 Genotyping-by-sequencing for estimating relatedness in nonmodel organisms: Avoiding the trap of precise biasMolecular Ecology Resources 18:381-390
• Cockerham CC1969 Variance of gene frequenciesEvolution 23:72-84
• Cockerham CC1973 Analyses of gene frequenciesGenetics 74:679-700
• Cockerham CC, Weir BS1986 Estimation of inbreeding parameters in stratified populationsAnnals of Human Genetics 50:271-281
• Goudet J2005 Hierfstat, a package for R to compute and test hierarchical F-statisticsMolecular Ecology Notes 5:184-186
• Goudet J, Weir BS2023 An allele-sharing, moment-based estimator of global, population-specific and population-pair FST under a general model of population structurePLoS Genetics 19:e1010871
• Long JC1986 The allelic correlation structure of Gainj- and Kalam-speaking people. I. The estimation and interpretation of Wright's F-statisticsGenetics 112:629-647
• Nei M1973 Analysis of gene diversity in subdivided populationsProceedings of the National Academy of Sciences 70:3321-3323
• Nei M1978 Estimation of average heterozygosity and genetic distance from a small number of individualsGenetics 89:583-590
• Nei M1986 Definition and estimation of fixation indicesEvolution 40:643-645
• Nei M1987 Molecular evolutionary genetics. Columbia University Press, New York.
• Nei M, Chesser RK1983 Estimation of fixation indices and gene diversitiesAnnals of Human Genetics 47:253-259
• R Core Team2025 R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. Available at <Available at https://www.R-project.org/ >. Accessed on January 10, 2026.
  » https://www.R-project.org/
• Weir BS1996 Genetic data analysis II: Methods for discrete population genetic data. Sinauer Associates, Sunderland, 445p.
• Weir BS, Cockerham CC1984 Estimating F-statistics for the analysis of population structureEvolution 38:1358-1370
• Weir BS, Goudet J2017 A unified characterization of population structure and relatednessGenetics 206:2085-2103
• Weir BS, Goudet J2026 Allelic association analyses: estimation recommendations. Genetics 233: iyag088. doi:10.1093/genetics/iyag088.
• Wright S1951 The genetical structure of populationsAnnals of Eugenics 15:323-354

Correspondence

Corresponding author: [e-mail]

SCIENTIFIC EDITOR:

Luiz Antônio dos Santos Dias

Publication Dates

• Publication in this collection
  05 Oct 2026
• Date of issue
  2026

History

• Received
  28 May 2026
• Accepted
  06 July 2026
• Published
  17 July 2026