Understanding population structure in an evolutionary context: population-specific FST and pairwise FST

Abstract Populations are shaped by their history. It is crucial to interpret population structure in an evolutionary context. Pairwise FST measures population structure, whereas population-specific FST measures deviation from the ancestral population. To understand the current population structure and a population’s history of range expansion, we propose a representation method that overlays population-specific FST estimates on a sampling location map, and on an unrooted neighbor-joining tree and a multi-dimensional scaling plot inferred from a pairwise FST distance matrix. We examined the usefulness of our procedure using simulations that mimicked population colonization from an ancestral population and by analyzing published human, Atlantic cod, and wild poplar data. Our results demonstrated that population-specific FST values identify the source population and trace the evolutionary history of its derived populations. Conversely, pairwise FST values represent the current population structure. By integrating the results of both estimators, we obtained a new picture of the population structure that incorporates evolutionary history. The generalized least squares estimate of genome-wide population-specific FST indicated that the wild poplar population expanded its distribution to the north, where daylight hours are long in summer, to coastal areas with abundant rainfall, and to the south where summers are dry. Genomic data highlight the power of the bias-corrected moment estimators of FST, whether global, pairwise, or population-specific, that provide unbiased estimates of FST. All FST moment estimators described in this paper have reasonable processing times and are useful in population genomics studies.


Introduction
Quantifying genetic relationships among populations is of substantial interest in population biology, ecology, and human genetics (Weir and Hill 2002). Appropriate estimates of population structure are the basis of our understanding of biology and biological applications, which vary from evolutionary and conservation studies to association mapping and forensic identification (Weir and Hill 2002). For such objectives, Wright's F ST (Wright 1951) is commonly used to quantify the genetic divergence of populations, and there have been many informative reviews of F ST estimators (e.g., Balloux and Lugon-Moulin 2002;Weir and Hill 2002;Rousset 2004Rousset , 2007Beaumont 2005;Excoffier 2007;Holsinger and Weir 2009;Gaggiotti and Foll 2010;Bhatia et al. 2013). The traditional F ST estimators have been defined as the ratio of the between-population variance to the total variance in allele frequencies (Wright 1965;Cockerham 1969Cockerham , 1973Weir and Cockerham 1984;Balloux and Lugon-Moulin 2002;Excoffier 2007;Holsinger and Weir 2009). An alternative approach for estimating population differentiation is to use population-specific F ST estimators (Balding and Nichols 1995;Nicholson et al. 2002;Weir and Hill 2002;Weir et al. 2005;Gaggiotti and Foll 2010;Weir and Goudet 2017). Model-based Bayesian approaches, based on beta and/or Dirichlet distributions, for estimating population-specific F ST have been proposed (Balding and Nichols 1995;Nicholson et al. 2002;Falush et al. 2003;Beaumont and Balding 2004). In addition to model-based methods, moment estimators of population-specific F ST have been derived (Weir and Hill 2002;Weir and Goudet 2017). A large number of approaches exist for estimating F ST that have different underlying assumptions (global, pairwise, or population-specific F ST ) and the framework used, such as frequentist and/or Bayesian. There have been many comprehensive reviews of traditional and populationspecific F ST estimators, as indicated above, most of which were written from the viewpoint of theoretical issues. Conversely, there has been no formal comparative study to describe their differences in terms of the evolutionary scenarios that best describe the data. Although these issues are well understood among statistical geneticists and theoretical population geneticists, empirical researchers, particularly those working in non-model organisms, could benefit from studies that address the problem.
In practice, the F ST value estimated from a set of population samples is called the global F ST , which measures population differentiation overall populations (e.g., Pé rez-Lezaun et al. 1997). Additionally, F ST values between pairs of population samples (pairwise F ST , Reynolds et al. 1983) are routinely used to estimate population structure in human genetics (Pé rez-Lezaun et al. 1997;Ramachandran et al. 2005), conservation biology (Palsbøll et al. 2007;Schwartz et al. 2007), and evolutionary biology and ecology (e.g., Hemmer-Hansen et al. 2013a;Therkildsen et al. 2013a;Geraldes et al. 2014;McKown et al. 2014a;Jorde et al. 2015;Rougemont et al. 2020). Divergent selection across an environmental gradient can influence population structure (Nosil et al. 2009;Orsini et al. 2013), and researchers have examined geographic distance and habitat differences between populations as explanatory variables that affect population structure estimated based on genome-wide (average over loci) pairwise F ST values (e.g., Rousset 1997;Bradbury and Bentzen 2007;Petrou et al. 2014;Jorde et al. 2015;Kitada et al. 2017). To identify the adaptive divergence of a gene among populations, locus-population-specific F ST was developed using empirical Bayes (Beaumont and Balding 2004) and full Bayesian methods (BayeScan) (Foll and Gaggiotti 2008). The Bayesian methods have been applied in many cases across various species to identify outlier single-nucleotide polymorphisms (SNPs) (e.g., Limborg et al. 2012;Therkildsen et al. 2013a;Geraldes et al. 2014). However, genome-wide populationspecific F ST is new among biologists. Despite the expected usefulness of genome-wide population-specific F ST in evolutionary biology (Weir and Goudet 2017), applications have been sparse to date (e.g., Nicholson et al. 2002;Weir et al. 2005;Foll and Gaggiotti 2006;Buckleton et al. 2016;Rougemont et al. 2020).
Traditional F ST estimators were originally developed to estimate F ST over a metapopulation (global F ST ) based on a set of population samples (Cockerham 1969;Nei and Chesser 1983;Weir and Cockerham 1984;Excoffier 2007;Rousset 2007). In this study, we use Nei and Chesser's (1983) bias-corrected G ST moment estimator (hereafter, NC83) as a pairwise F ST estimator (Supplementary Note). The pairwise F ST between populations (i, j) is defined as: where H ij T is total heterozygosity over all populations and H ij S is within-population heterozygosity.
We apply Weir and Goudet's (2017) bias-corrected moment estimator of population-specific F ST (hereafter, WG) (Supplementary Note). When only allele frequencies are used, the WG population-specific F ST at a locus is defined as: where M i W is the within-population matching of two distinct alleles of population i and M B is the between-population-pair matching average over pairs of populations i; i 0 . M B is homozygosity over pairs of populations. We represent heterozygosity over all pairs of populations as 1 À M B ¼ H B , and 1 À M i W ¼ H Si . Therefore: This formulation is reasonable because WG populationspecific F ST uses "allele matching, equivalent to homozygosity and complementary to heterozygosity as used by Nei (1973), rather than components of variance" (Weir and Goudet 2017). H B is heterozygosity for all pairs of populations, whereas H ij T in Equation (1) is heterozygosity for the pair of populations. Equation (2) shows that WG population-specific F ST measures population-specific genetic diversity (H Si ) under the framework of the relatedness of individuals and identifies the population with the greatest genetic diversity as the ancestral or oldest population. Because populations close to the ancestral population have had more opportunities for mutations than recently founded populations (Liu et al. 2006), they are likely to have the highest heterozygosity and low values of population-specific F ST . Thus, WG population-specific F ST works to infer evolutionary history through genetic diversity in terms of heterozygosity under the assumption that the ancestral population had the highest genetic diversity. By combining population-specific and pairwise F ST estimates, we can infer the present population structure, which reflects evolutionary history.
In this study, our objective is to demonstrate to empirical population geneticists and biologists how the two types of genomewide F ST estimators can be combined to help elucidate the population structure (pairwise F ST ) in the evolutionary context (population-specific F ST ). In our approach, the current population structure is represented by an unrooted neighbor-joining (NJ) tree (Saitou and Nei 1987) and a multi-dimensional scaling (MDS) plot based on pairwise F ST values, with population history being inferred by overlaying population-specific F ST values on the population structure. The colors of the populations (names and/or sampling points) based on the WG genome-wide population-specific F ST values enable the inference of the historical order of population colonization. We also present a representation on a geographical map, which is useful for visually understanding population history in a distribution range.
First, we examine the usefulness of our procedure using stepping-stone simulations that mimic population colonization from a single ancestral population for five scenarios of population range expansion. We then apply our approach to three datasets of human, Atlantic cod (Gadus morhua), and wild poplar (Populus trichocarpa). Human evolutionary history, migration, and population structure have been particularly well studied (e.g., Diamond 1997;Rosenberg et al. 2002;Ramachandran et al. 2005;Liu et al. 2006;Pickrell and Pritchard 2012;Lipson et al. 2013;Hellenthal et al. 2014;Rutherford 2016;Nielsen et al. 2017). These patterns are well known by statistical/theoretical population geneticists and biologists; therefore, testing our integrative F ST analysis on this dataset could provide a good example of the usefulness of this practical approach. Although dense human SNP datasets are currently available, we used microsatellite data for illustrative purposes, because the quality of human microsatellites has been examined thoroughly, and they are fundamentally neutral (Kanitz et al. 2018). Also, microsatellites are highly polymorphic and have more information at a locus than SNPs (Schlö tterer 2004).
The Atlantic cod SNP were genotyped from mature fish samples collected from the North Atlantic from the northern range margin of the species in Greenland, Norway, and the Baltic Sea. Two ecotypes (migratory and stationary) that were able to interbreed were genetically differentiated Berg et al. 2016). The inclusion of both types of data may improve the understanding of the demographic history of highly migratory marine fish. The wild poplar SNP data were collected from the American Pacific Northwest. Male poplar trees produce pollen and female trees produce small seeds with fine hairs, which enable dispersal of this species by wind (Geraldes et al. 2014). The samples covered various regions over a range of 2500 km near the Canadian-US border in British Columbia (BC) at alt between 0 and 800 m, where the variations in photoperiod and temperature have a north-south cline, while the variations in temperature, rainfall, and drought have an east-west (coastal to inland) cline (Geraldes et al. 2014). Each sampling location was associated with 11 environmental and geographical parameters. The analysis of environmental variables and population-specific F ST values may provide a good example for understanding the history of the range expansion of a wind-dispersed species. By analyzing different types of data with species-specific ecology and migration history, the usefulness of our approach may be identified to enable us to understand the current population structure in an evolutionary context.

Population colonization simulations
To test the performance of our visual representation, we conducted simulations that mimicked the colonization of populations from a single ancestral population (population 1). We modeled five types of stepping-stone colonization: one, two, and three-directional population expansion; three-directional grid colonization from an edge; and eight-directional grid colonization from the center, with 24 demes (populations 2-25) (Figure 1, A-E). Our objective is to describe current population structure using an unrooted NJ tree and an MDS plot, and to infer population history by overlaying population-specific F ST values on the population structure. When gene flow is limited between adjacent populations, as shown in our simulations, the estimated population structure corresponds to population history.
We set the effective population size of the ancestral population to N e ¼ 10 4 (twice the number of individuals in diploid organisms in a random mating population). At the beginning of colonization, 1% of N e migrated into the adjacent vacant habitat once every 10 generations. For convenience, we considered one simulation cycle to be one generation. The effective population size of the newly derived population increased to the same size as the ancestral population (N e ¼ 10 4 ) after one generation, and the population exchanged 1% of N e genes with adjacent population(s) in every generation. Like the ancestral population, 1% of N e individuals migrated into the adjacent vacant habitat once every 10 generations. We simulated the allele frequencies of SNPs in the ancestral and 24 derived populations. We also examined the cases in which the effective population size of the ancestral population was 10 times greater (N e ¼ 10 5 ) than that of the newly derived population (N e ¼ 10 4 ).
We generated the initial allele frequencies in the ancestral population, q, at 100,000 neutral SNP loci from the predictive equilibrium distribution, f q ð Þ / q À1 1 À q ð Þ À1 (Wright 1931). Additionally, we introduced 10 newly derived SNPs to each existing population in each generation. When a new SNP emerged in a population, we set the initial allele frequency of the newly derived SNP to 0.01 in the population and 0 in the other populations. This mimicked new mutations that survived in the initial phase after their birth. We considered these 100,000 ancestral SNPs and newly derived SNPs to be "unobserved." We changed the allele frequencies of these SNPs using random drift under a binomial distribution in every generation. The frequencies of the derived alleles decreased for many of the SNPs over simulation generations and lost their polymorphism. After 260 simulation generations, we randomly selected SNPs that retained their polymorphism as "observed" SNPs. For grid colonization models (Figure 1, D-E), we randomly selected polymorphic SNPs after 100 simulation generations.
In our simulations, we selected 10,000 ancestral SNPs and 500 newly derived SNPs. Then, we generated 50 individuals for each population. We randomly generated the genotypes of these 10,500 SNPs for each individual following the allele frequencies in the population to which each individual belonged. Thus, we obtained "observed" genotypes of 1,250 individuals (¼ 50 individuals Â 25 populations) at 10,500 SNP loci. To examine the effect of generations on genetic diversity in newly derived SNPs, we also selected 9000 and 7000 ancestral SNPs, and 1000 and 3000 newly derived SNPs, respectively. We converted the simulated genotypes into Genepop format (Raymond and Rousset, 1995;Rousset, 2008). We then computed genome-wide population-specific and pairwise F ST values between the 25 populations.

Visual representation of population structure and demographic history
We integrated genome-wide population-specific and pairwise F ST estimates on a map of sampling locations on an NJ tree and MDS plot. We visualized the magnitude of the genome-wide population-specific F ST values using a color gradient based on rgb This conversion represents the standardized magnitude of a population-specific F ST value at the sampling point, with colors between red (for the smallest F ST ) and blue (for the largest F ST ). We drew the F ST map using the sf package in R, where we plotted sampling locations based on the longitudes (lon) and lat. We connected sampling points with genome-wide pairwise F ST values smaller than a given value using yellow lines to visualize the connectivity between populations. Under the assumption of Wright's island model at equilibrium between drift, mutation, and migration (Wright 1931 where N e is the effective population size and m is the average rate of migration between populations (Slatkin 1987). For example, F ST ¼ 0:02 refers to 4N e m % 49 migrants per generation (see Whitlock and Mccauley 1999;Waples and Gaggiotti 2006). The value was arbitrarily used in our case studies. We plotted the genome-wide population-specific F ST values on a dot chart with standard errors estimated using Equation (S5) (Supplementary Note). We drew the NJ tree based on the distance matrix of the genome-wide pairwise F ST values using the nj function in the R package ape. We performed MDS analysis on the pairwise F ST distance matrix using the cmdscale function in R. We used the cumulative contribution ratio up to the kth axis j ¼ 1; . . . ; k; . . . ; K À Á as the explained variation measure, which we computed using the R function as We colored the sampling locations on the F ST maps, dot charts, NJ trees, and MDS plots using a color gradient of the magnitude of genome-wide population-specific F ST values. We also examined a diverging color palette instead of blue to red to test the resolution using RColorBrewer on the F ST maps.

Computing F ST values
We converted the genotype data into Genepop format (Raymond and Rousset 1995;Rousset 2008) (Foll and Gaggiotti 2006) to compute the genome-wide population-specific F ST values. We examined the shrinkage effect of the Bayesian population-specific F ST estimator on inferring the ancestral population using a set of subsamples (37 populations) chosen from 51 populations.

Inferring environmental effects on the observed population structure
To infer the geography and environment that were experienced by the population range expansion, we regressed the genomewide population-specific F ST values on the geographical and environmental variables (j ¼ 1; . . . ; s): Figure 1 Results from population colonization simulations. Schematic diagrams of the models: (A) one, (B) two, (C) three-directional colonization, (D) three-directional grid colonization, and (E) eight-directional grid colonization. Population 1 in red is ancestral, and the yellow arrows indicate the direction of colonization. Lines show opportunities for migration. The effective population size of the newly derived population increased to the same size as the ancestral population (N e ¼ 10 4 ) after one simulation generation, and each population exchanged 1% of N e genes with adjacent population(s) in every generation, as indicated by the arrows (see the text). Neighbor-joining (NJ) unrooted trees (F-J) and multi-dimensional scaling (MDS) plots (K-O) based on the pairwise F ST distance matrix overlaid with population-specific F ST values for each model. The color of each population shows the magnitude of population-specific F ST values between red (for the smallest F ST ) and blue (for the largest F ST ).
where X is the variance matrix of population-specific F ST . We correlated residuals because of the population structure; therefore, the effective sample size was lower than the actual sample size. In such circumstances, ordinary least squares overestimate the precision. To account for the correlation, we derived the components of the variance-covariance matrix of the populationspecific F ST estimator [Supplementary Equations (S5) and (S6)] for generalized least squares (GLS). We performed this analysis on the wild poplar dataset, for which 11 environmental/geographical parameters were available for each sampling location. We used the variance-covariance matrix for the components of the variance matrix X in Equation (3), and performed regression using the GLS function in FinePop2 v0.2.

Three empirical datasets
We retrieved the human microsatellite data used in Rosenberg et al. (2002) from https://web.stanford.edu/group/rosenberglab/in dex.html. We removed the Surui sample (Brazil) from the data because that population was reduced to 34 individuals in 1961 as a result of introduced diseases (Liu et al. 2006). We retained genotype data (n ¼ 1035) of 377 microsatellite loci from 51 populations categorized into six groups, as in the original study: 6 populations from Africa, 12 from the Middle East and Europe, 9 from Central/ South Asia, 18 from East Asia, 2 from Oceania, and 4 from America. We obtained the lon and lat of the sampling sites from Cann et al. (2002) (Supplementary Table S1).
We combined the Atlantic cod SNP genotype data of 924 SNPs common to 29 populations reported in Therkildsen et al. (2013aTherkildsen et al. ( , 2013b and 12 populations reported in Hemmer-Hansen et al. (2013a. We compared the genotypes associated with each marker in samples that were identical in the two studies, that is, CAN08 and Western_Atlantic_2008, ISO02 and Iceland_migratory_2002, and ISC02 and Iceland_stationary_2002, and standardized the gene codes. We removed cgpGmo.S1035, whose genotypes were inconsistent between the two studies. We also removed cgpGmo.S1408 and cgpGmo.S893, for which the genotypes were missing in several population samples in Therkildsen et al. (2013b). For simplicity, we removed temporal replicates from the Norway migratory, Norway stationary, North Sea, and Baltic Sea samples. The final dataset consisted of genotype data (n ¼ 1065) for 921 SNPs from 34 populations: 3 from Iceland, 25 from Greenland, 3 from Norway, and 1 each from Canada, the North Sea, and the Baltic Sea. All individuals in the samples were adults, and most were mature (Therkildsen et al.  Table S2).
We retrieved wild poplar SNP genotype data and environmental/geographical data from the original studies of McKown et al. (2014aMcKown et al. ( , 2014b. The genotype data contained 29,355 SNPs of 3,518 genes of wild poplar (n ¼ 441) collected from 25 drainage areas (McKown et al. 2014c). Details of the array development and selection of SNPs are provided in Geraldes et al. (2011Geraldes et al. ( , 2013. A breakdown of the 25 drainages (hereafter, populations) is as follows: 9 in northern British Columbia (NBC), 2 in inland British Columbia (IBC), 12 in southern British Columbia (SBC), and 2 in Oregon (ORE) (Geraldes et al. 2014). We combined the original names of the clusters and population numbers, and used them as our population labels (NBC1, NBC3,. . ., ORE30). We associated each sampling location with 11 environmental and geographical parameters: lat, lon, alt, longest day length (DAY), frost-free days (FFD), mean annual temperature (MAT), mean warmest month temperature (MWMT), mean annual precipitation (MAP), mean summer precipitation (MSP), annual heat-moisture index (AHM), and summer heat-moisture index (SHM) (Supplementary Table  S3). The AHM was calculated in the original study as (MAT þ 10)/ (MAP/1000); a large AHM indicates extremely dry conditions.

Simulations of population colonization
First, we examined the effect of the number of simulation generations on genetic diversity in newly derived SNPs using the eightdirectional grid simulation ( Figure 1E). WG population-specific F ST correctly identified the ancestral population and traced the population history, and population structure reflected the population history regardless of the numbers of ancestral SNPs (9000 and 7000) and newly derived SNPs (1000 and 3000) selected after 100 simulation generations (Supplementary Figure S1). The result was consistent with the case that used 10,500 SNP loci (10,000 ancestral SNPs þ 500 newly derived SNPs) (Figure 1, J and O and Supplementary Figure S2E). In the following analysis, we used the results based on 10,500 SNP loci to generate clearer results, even for limited numbers of simulation generations (260 and/or 100).
In the one-directional simulation ( Figure 1A), populationspecific F ST correctly identified the ancestral population with the highest genetic diversity (Supplementary Figure S2A), and populations were located in order from 1 to 25 on the NJ tree ( Figure 1F). The first axis of the MDS plot explained 93% of the variance of the pairwise F ST distance matrix and indicated population expansion from populations 1 to 25 ( Figure 1K). In the twodirectional simulation ( Figure 1B), our analysis correctly identified the ancestral population (Supplementary Figure S2B) and detected that populations were split at population 9 and expanded in two directions ( Figure 1G), which was consistent with the simulation scenario. The first axis of the MDS plot identified population expansion from populations 1 to 25 and explained 56% of the variance of the pairwise F ST distance matrix, whereas the second axis identified another manner of population expansion from populations 1 to 16 and explained 40% of the variation ( Figure 1L). In the three-directional simulation ( Figure 1C), the ancestral population was also correctly identified (Supplementary Figure S2C). It was closely located to the adjacent populations 2, 9, and 17, but correctly detected three directions ( Figure 1H). The first axis of the MDS plot identified population expansion from population 1 to populations 16 and 25, and explained 57% of the variance of pairwise F ST , whereas the second axis identified population expansion from population 1 to populations 8 and 16, and explained 38% of the variance ( Figure 1M).
In the three-directional grid colonization model from an edge ( Figure 1D), population-specific F ST correctly identified the ancestral population (Supplementary Figure S2D) and pairwise F ST detected that populations expanded in three directions ( Figure 1I), which agreed with the simulation scenario. The first axis of the MDS plot identified population expansion from population 1 to other edge populations (populations 11, 15, and 25), and explained 60% of the variance of the pairwise F ST distance matrix, whereas the second axis indicated genetic differentiation between populations 24 and 25, and 15, 19, and 22, and explained 17% of the variance ( Figure 1N). In the eight-directional grid colonization model ( Figure 1E), population-specific F ST identified the ancestral population (Supplementary Figure S2E) and pairwise F ST estimated that populations expanded in five directions from the center ( Figure 1J). The first axis of the MDS plot identified vertical population expansion from population 1 to populations 24 and 25 and explained 40% of the variance of the pairwise F ST distance matrix, and the second axis indicated horizontal population expansion from population 1 to populations 23 and 24, which explains 30% of the variance ( Figure 1O). We obtained similar results in the cases in which the effective population size of the ancestral population was 10 times greater (N e ¼ 10 5 ) than that of the newly derived population (N e ¼ 10 4 ) (Supplementary Figure S3). We obtained very similar results from different data for more than 20 simulations (figures not shown).

Humans
The F ST map (Figure 2A) shows integrated information from genome-wide population-specific and pairwise F ST , which visualizes population structure with the migration and demographic history of human populations in terms of genetic diversity. Interestingly, Bantu Kenyans had the smallest F ST value (shown in red). Figure 2B ordered population-specific F ST values from Africa to Central/South Asia, the Middle East, Europe, East Asia, Oceania, and America (Supplementary Table S4). As indicated by the sampling points connected by yellow lines with pairwise F ST values below 0.02 (Figure 2A), genetic connectivity from Africa was low. Conversely, migration was substantial within Eurasia but much smaller than that inferred from Eurasia to Oceania and America. H e was the highest in Africa, followed by the Middle East, Central/South Asia, Europe, and East Asia, but relatively small in Oceania and lowest in South America. The Karitiana in Brazil had the lowest H e . The NJ tree ( Figure 2C populations and explained 44% of the variance of the pairwise F ST distance matrix, whereas the second axis indicated genetic differentiation between Melanesian and Karitiana populations, and explained 19% of the variance ( Figure 2D).
The Bayesian population-specific F ST values estimated using the methods of Beaumont and Balding (2004) (empirical Bayes) and Foll and Gaggiotti (2006) Figure  S5A, Supplementary Table S4). However, in African populations, they were higher than the WG population-specific F ST values ( Figure  S5B). Our F ST map based on the empirical Bayesian populationspecific F ST values indicated that the Middle East, Europe, and Central/South Asia were centers of human origin ( Figure 3A), which was consistent with that from the full Bayesian population-specific F ST estimator (figure not shown). Our integrated NJ tree showed that the Hazara, Pakistan population was genetically closest to the human ancestors (Supplementary Figure S6A). The numbers of sampling locations of the 51 human populations were as follows: 6 from Africa, 12 from the Middle Figure 3 Population structure of humans based on Bayesian and moment estimators of population-specific F ST . Results from the Bayesian populationspecific F ST estimator using (A) 51 samples and (B) 37 subsamples, and from (C) the WG population-specific F ST moment estimator using 37 subsamples. The numbers of sampling locations of the subsamples were as follows: 3 from Africa, 6 from the Middle East/Europe, 9 from Central/South Asia, 18 from East Asia, 2 from Oceania, and 4 from America. Populations connected by yellow lines are those with pairwise F ST < 0.02. The color of each population indicates the magnitude of population-specific F ST values between red (for the smallest F ST ) and blue (for the largest F ST ).
East/Europe, 9 from Central/South Asia, 18 from East Asia, 2 from Oceania, and 4 from America. When we used a subsample of the 37 human populations (3 from Africa, 6 from the Middle East/ Europe, 4 from Central/South Asia, 18 from East Asia, 2 from Oceania, and 4 from America; Supplementary Table S5), the area with the highest population-specific F ST values shifted toward Central/South Asia and East Asia ( Figure 3B), whereas Bantu Kenyans had the smallest WG population-specific F ST value ( Figure 3C); this was consistent with the results from the full dataset (Figure 2A). The integrated NJ trees provided similar results (Supplementary Figure S6, B and C).

Atlantic cod
The F ST map ( Figure 4A) visualizes the population structure, migration, and genetic diversity of the Atlantic cod populations. The Canadian population had the smallest population-specific F ST value (shown in red) and the greatest H e . H e was also high in Greenland, low in other areas, and lowest in the Baltic Sea. Figure 4B shows the order of population-specific F ST values from Canada to the Baltic sea (Supplementary Table S6). Greenland west coast populations (green in Supplementary Figure S7) generally had small population-specific F ST values, whereas fjord populations (violet) had relatively higher population-specific F ST values. The population-specific F ST values were much higher for populations in Iceland, Norway, and the North Sea, and were highest in the Baltic Sea. Based on pairwise F ST values (< 0.02) between sampling points ( Figure 4A), substantial migration was suggested between Greenland, Iceland, and Norway. Conversely, migration could be low from Canada to Greenland and from Iceland and Norway to the North and Baltic Seas. Our integrated NJ tree with population-specific F ST values ( Figure 4C) inferred that Atlantic cod originated from Canada, migrated to the west coast of Greenland, and then expanded their distribution to Iceland, Norway, the North Sea, and the Baltic Sea. According to the ordinal NJ tree of the pairwise F ST distance matrix (Supplementary Figure S7), the populations were divided into four large clusters: (1) Canada; (2) Greenland west coast, (3) Greenland east coast, Iceland, and Norway; and (4) North and Baltic Seas. Fjord populations (in purple) formed a sub-cluster within the Greenland west coast, and migratory (orange) and stationary (magenta) ecotypes also formed a sub-cluster. The first axis of the MDS plot characterized the differentiation between Canadian and North Sea/Baltic Sea populations and explained 72% of the variance of the pairwise F ST distance matrix, whereas the second axis highlighted the differentiation between Norwegian migratory populations and Canadian and North Sea/ Baltic Sea populations, which explained 22% of the variance ( Figure 4D).

Wild poplar
The F ST map ( Figure 5A) indicated that population-specific F ST values were lowest in SBC27 and inner British Columbia (IBC15 and IBC16) (shown in red, Supplementary Figure S8). The sampling points connected by yellow lines (pairwise F ST < 0.02) indicated migration between all populations. H e was highest in SBC27, IBC15, and IBC16, and lowest in NBC5. Figure 5B shows that samples collected from areas close to the SBC coast had higher population-specific F ST values than other SBC samples (Supplementary Table S7). The NBC samples had populationspecific F ST values similar to those of SBC. Among the NBC samples, NBC8 had the smallest population-specific F ST , and NBC5 had the highest value, followed by NBC6 and NBC7. The pairwise F ST NJ tree integrated with population-specific F ST values ( Figure 5C) suggested that wild poplar originated from around SBC27 and interior BC, and expanded in three directions: to the southern coast of BC, NBC and south-western Alaska, and ORE. The ordinal NJ tree based on the pairwise F ST distant matrix divided populations into four large clusters: (1) IBC, (2) SBC, (3) NBC, and (4) ORE (Supplementary Figure S8). The population represented by sample ORE30 was isolated from ORE29. The first axis of the MDS plot characterized the differentiation between southern and northern populations and explained the 51% variance of the pairwise F ST distance matrix, whereas the second axis characterized the southernmost ORE30 population, and explained the 18% variation ( Figure 5D).
To avoid multicollinearity, we excluded 7 out of 11 environmental variables that were significantly correlated with each other: lat, lon, alt, FFD, MWMT, MSP, and AHM (Supplementary  Table S3). Our GLS of genome-wide population-specific F ST values on the four environmental variables (DAY, MAT, MAP, and SHM) indicated that DAY, MAP, and SHM were significant (Table 1). All estimates were positive, which indicated that higher populationspecific F ST values were expected for longer DAY (longer daylight time), higher MAP (abundant rain), and higher SHM (dry summers), and these values might reflect the directions of population expansion. The scatter plot of DAY and SHM (with each population colored according to the population-specific F ST value) ( Figure 6A) suggested three directions of population range expansion: the wild poplar that might have originated from around  SBC27 and IBC15 expanded its distribution to NBC, where daylight hours are long in summer, as well as expanding to coastal SBC with its lower SHM (humid summer and abundant rainfall), and to the south (ORE29 and ORE30) with its higher SHM (dry summer). This was consistent in the scatter plot of DAY and MAP ( Figure 6B).

Diverging color palette
RColorBrewer had 35 color palettes. Each palette had a minimum of eight colors. Two palettes had 12 colors (maximum) and nine had 11 colors. We chose a color palette with 10 colors (RdYlBu). The color gradient better identified the middle range of population-specific F ST values, but failed to detect the ancestral population (Supplementary Figure S9).

CPU times
With an Intel Core i7-8650U CPU, 89.8 s of CPU time were required to compute the WG population-specific F ST estimates and SEs of wild poplar (29,355 SNPs; 25 populations, n ¼ 441). To obtain the pairwise F ST (NC83) between all 300 (¼ 25 Â 24/2) population pairs, 120.7 s were required. CPU time is proportional to the number of SNPs analyzed.

Discussion
Genome-wide population-specific F ST traced population history as reflected by genetic diversity Our simulations demonstrated that the WG population-specific F ST estimator identified the source population and traced the evolutionary history of its derived populations based on genetic diversity (heterozygosity estimated from each population). The NC83 pairwise F ST estimator correctly estimated the current population structure. As explained in the Introduction, the population-specific F ST estimator is a rescaling of expected heterozygosity, and we expect a linear relationship between expected heterozygosity and population F ST . This shows that the population-specific F ST estimator implicitly assumes that populations closest to the ancestral population have the highest heterozygosity. In our three case studies, a linear relationship between H e of each population (¼ H Si ) and psF i ST was evident (Supplementary Figure S10). The coefficient of determination, R 2 , was 0.91 for 51 human populations (n ¼ 1035), 0.99 for 34 Atlantic cod populations (n ¼ 1065), and 0.82 for 25 wild poplar populations (n ¼ 441). The goodness of fit to the linear function should depend on the sample size (number of individuals). Our simulations evaluated the performance of the population-specific F ST estimator for such cases. However, in populations that experienced extensive admixture events, heterozygosity was enhanced, whereas a bottleneck in the ancestral population reduced heterozygosity. In such cases, the population-specific F ST estimator misidentifies the ancestral population.
In our analysis, the genome-wide WG population-specific F ST values successfully illustrated human evolutionary history, and indicated that humans originated in Kenya, expanded from the Middle East into Europe and from Central/South Asia into East Asia, and then possibly migrated to Oceania and America ( Figure 2). Kenya is located just below Ethiopia, where the earliest anatomically modern humans were found from fossils (Nielsen et al. 2017). Our results are also in good agreement with the highest levels of genetic diversity being detected in Africa (Rosenberg et al. 2002), the relationship uncovered between genetic and geographic distance (Ramachandran et al. 2005), the shortest colonization route from East Africa (Liu et al. 2006), and major migrations inferred from genomic data (Nielsen et al. 2017).
Our analysis indicated that Atlantic cod might originate in Canada (CAN08). Figure 4 suggested that the population expansion of Atlantic cod began by minimal gene flow from Canada. They might have first expanded to the west coast of Greenland before spreading to Iceland, the North Sea, Norway, and the Baltic Sea. This result was consistent with genomic evidence that Atlantic cod inhabit both sides of the Atlantic Ocean and evolved from a common evolutionary origin (Berg et al. 2017). The migratory ecotypes characterized by deeper and more offshore habitats  Figure 4C and Supplementary Figure S7. In our study, CAN08 had the highest H e , which was lower in Iceland than in Greenland; this result implies that Icelandic populations were the descendants of colonists from Greenland, which in turn originated in Canada. The BAS0607 sample from the Baltic Sea had the highest population-specific F ST and the lowest H e , suggesting that Baltic cod is the newest population. This result agrees with the findings of a previous study, which identified Baltic cod as an example of a species subject to ongoing selection for reproductive success in a low salinity environment (Berg et al. 2015). In the Atlantic cod case study, CAN08 had the highest H e and a very large negative population-specific F ST value of À0.21 6 0.019 compared with the maximum value of 0.22 6 0.014 in BAS0607 (Supplementary Figure S10, Supplementary Table S6). The WG population-specific F ST value can be negative (Weir and Goudet 2017). In the one and two-directional models of our simulations, the WG population-specific F ST value was significantly negative in the ancestral population, whereas H e was the largest (Supplementary Figure S2, A and B). Our consistent results between the simulations and Atlantic cod case study indicate that when gene flow from other populations into the source population is limited, a relatively large H e (Ĥ Si ) is maintained in the source population. In such cases withĤ Si >Ĥ B , Equation (2) produces negative values for population-specific F ST .
Although the wild poplar samples used in this study might not cover the entire distribution range of the species, which extends from southern California to northern Alaska, Montana, and Idaho (Geraldes et al. 2013), the genome-wide population-specific F ST values suggested three directions of population expansion of wild poplar: from SBC27 and IBC (IBC15, IBC16) to coastal British Colombia, southern ORE, and NBC ( Figure 5). The largest population-specific F ST value was found in the population with the smallest heterozygosity, SBC22, which may have resulted from a bottleneck (Geraldes et al. 2014).
Our continuous color gradient from red to blue successfully detected the ancestral population, whereas the R diverging color palettes, which had a limited number of colors (maximum of 12), better identified the middle range of population-specific F ST values (Supplementary Figure S9). Both color gradients may be useful.
Genome-wide pairwise F ST described current population structure Our stepping-stone simulations did not account for long-range dispersal (Hallatschek and Fisher 2014). Human samples were collected from present populations, but they might reflect history from the Age of Discovery, when humans were travelling far beyond their native continents (Diamond 1997). Wild animals primarily move locally, but occasionally disperse over long distances (Hallatschek and Fisher 2014). The migratory ecotype of the Atlantic cod is characterized by its long-distance migration . Sea currents also play a role for passive transportation of fertilized eggs and juveniles. The pollen and small seeds with fine hairs of poplar trees enable long-range dispersal by wind (Geraldes et al. 2014). Our stepping-stone simulations were conducted with the assumption that dispersal was step-by-step, but long-range dispersals were not taken into consideration. When gene flow is limited between adjacent populations, as in our simulation scenarios, estimated population structure reflects population history.
Genome-wide population-specific F ST detects key environments that relate to population expansion Our GLS of genome-wide population-specific F ST values revealed that long daylight hours, abundant rainfall, and dry summer conditions are the key environmental factors that relate to the demographic history of wild poplar (Table 1). Wild poplar could have expanded its distribution by its fluffy seeds being blown away by the wind. In the NJ unrooted tree, the root of the populations cannot be fixed without out-group populations. There is no direct evidence, but we can infer from the population-specific F ST values that wild poplar seems to have spread from SBC (SBC27) and IBC (IBC15, 16), and expanded its distribution to NBC, where daylight hours are long in summer, to coastal SBC with its rainy environment, and to southern ORE (ORE30) with its dry summer conditions (Figures 5 and 6). A previous study on wild poplar revealed that genes involved in drought response were identified as F ST outliers using BayeScan (Foll and Gaggiotti 2008;Geraldes et al. 2014). The F ST outlier test of Geraldes and colleagues also revealed that genes involved in the circadian rhythm and response to red/far-red light had high locus-specific global F ST values. The first principal component of SNP allele frequencies of the poplar tree was significantly correlated with day length, and a previous enrichment analysis for population structuring uncovered genes related to circadian rhythm and photoperiod. The circadian clock pathway might play a central role in adaptation, and clinal variation in phenological traits might be an adaptive response to the north-south climatic variation of P. trichocarpa (Geraldes et al. 2014;McKown et al. 2014a). Our GLS analysis does not detect a locus-specific effect of environmental adaptation like genome scan methods, but detects key environmental variables that affected the population history through genome-wide population-specific F ST . Our results confirmed the previous findings and the hypothesis of postglacial northward expansion of the poplar tree to refugia north of the southern margin of the ice sheet (Geraldes et al. 2014), showing the usefulness of applying the GLS estimate of genome-wide population-specific F ST to infer environmental effects on the population expansion of species.
The two types of F ST s can also be calculated for the genomic windows. By comparing the pairwise F ST among genomic windows, it is possible to identify the genomic regions that are largely differentiated due to local adaptation by using the population branch statistics (Yi et al. 2010). The local principal component analysis using lostruct software (Li and Ralph 2019) and MDS (Fuller et al. 2020) also identify genomic regions linked to local adaptation. Likewise, by comparing population-specific F ST among windows, it is possible to identify unique genomic regions for adaptation (Akey et al. 2002;Weir et al. 2005;Oleksyk et al. 2008).

Genome-wide F ST moment estimators converge to their true means
Previous studies have suggested or indicated that the "ratio of averages" works better than the "average of ratios" as the number of independent SNPs increases (Reynolds et al. 1983;Weir and Cockerham 1984;Bhatia et al. 2013). Because "the combined ratio estimate (ratio of averages) is much less subject to the risk of bias than the separate estimate (average of ratios)" (Cockran 1977), scholars recommend using the "ratio of averages" estimators (Bhatia et al. 2013). To explicitly show the underlying mechanism, we used the observed heterozygosity of population i (Ĥ Si ) as derived by Nei and Chesser (1983) (Supplementary Note). When the number of loci (L) increases, the average observed heterozygosity over all loci converges to its expected value according to the law of large numbers as: Similarly,Ĥ S andĤ T converge to their expected values. This example indicates that the numerators and denominators of bias-corrected F ST moment estimators, whether global, pairwise, or population-specific, converge to their true means and provide unbiased estimates of F ST in population genomics analyses with large numbers of SNPs.
Bayesian F ST estimators measure the deviation from the average of the sampled populations In the Bayesian framework, the population-specific F ST is the coefficient of the genetic drift that represents the amongpopulation variation of the allele frequencies at neutral loci from the allele frequencies of the ancestral population. The allele frequency of the ancestral population is assumed to be the amongpopulation mean allele frequency. The analysis of human populations ( Figure 3A) highlights the need to account for geographical heterogeneity in sampling fractions of populations. The unbiased estimate of the allele frequency of the ancestral population will be obtained as the weighted among-population average. The weights are inversely proportional to the sampling fractions.
The shrinkage effect on allele frequencies in Bayesian inference (Stein, 1956) may shift population-specific F ST values toward the average of the entire population. Because of the shrinkage toward mean allele frequencies, the maximum likelihood and Bayesian estimators of locus-specific global F ST improve the power to detect genes under environmental selection (Beaumont and Balding 2004;Foll and Gaggiotti 2008). An empirical Bayes genome-wide pairwise F ST estimator (Kitada et al. 2007) is useful in cases involving a small number of polymorphic marker loci, particularly in high gene flow scenarios, but it suffers from the shrinkage effect when larger numbers of loci are used. The shrinkage of allele frequencies should affect inference in genome-wide population-specific F ST , particularly in cases when samples (populations) were not representative of the populations.

Conclusions
The WG genome-wide population-specific F ST moment estimator can identify the source population and trace the evolutionary history of the derived populations based on genetic diversity under the assumption that populations closest to the ancestral population have the highest heterozygosity. Conversely, the NC83 genome-wide pairwise F ST moment estimator represents the current population structure. By integrating population-specific and pairwise estimates on F ST maps, NJ trees, and MDS plots, we obtain a picture of population structure by incorporating evolutionary history. Our GLS analysis of genome-wide population-specific F ST , which accounts for the correlation between populations, provides insights into how a species expanded its distribution in different environments. Given a large number of loci, bias-corrected F ST moment estimators-whether global, pairwise, or population-specific-provide unbiased estimates of F ST supported by the law of large numbers. Genomic data highlight the usefulness of the biascorrected moment estimators of F ST .

Data availability
The authors affirm that all data necessary for confirming the conclusions of the article are present within the article, figures, tables, and supplementary material. The R codes for our representation method exemplified by the human data and simulations of population colonization used in this study are available in the Supplementary material at figshare: https://doi.org/10. 25387/g3.14813490.