Fine-scale structure of a whole regional population through genetics and genealogies

Domain: chemistry
Published: 2026-03-02
Editorial outcome: major_revision

Abstract
Population stratification can confound genetic association studies and often persists despite adjustment using principal components of common variants. Demographic history and rare variants can also contribute to this confounding. But to what extent does demographic history impact fine structure and can be detected in human populations? To address this question, we analysed the Saguenay–Lac-Saint-Jean region of Quebec, a recent founder population long assumed to be homogeneous. Integrating genotype data with genealogical records, we show a strong concordance between realised (genetic) and expected (genealogical) kinship. Using a time-efficient algorithm capable of computing billions of pairwise kinship coefficients for all individuals married in the region between 1931 and 1960, we reveal fine structure at the municipal level, including an east–west genetic gradient shaped by differential founders’ contribution, migration patterns and socioeconomic factors. These findings challenge the assumption of regional homogeneity and suggest that similar recent structure likely exists in other populations worldwide. These signals may be obscured by coarse, ancient structure under standard stratification corrections used in genome-wide association studies and polygenic risk score analyses which can lead to biases and false associations when realised allele frequencies correlate with phenotypic variation.

Introduction
Population stratification occurs when genetic structure (non-random allele distribution) and environmental differences between subpopulations confound genetic association studies. This can make environmental effects appear genetic, leading to misleading results in genome-wide association studies (GWAS) and polygenic scores1,2. Family-based GWAS, which controls for a shared familial environment, reveals that standard population GWAS often overestimates genetic effects3. This provides evidence that population stratification may persist despite conventional population structure controls, challenging assumptions in large-scale GWAS. When these GWAS are used to construct polygenic scores, even small biases can accumulate into large errors. Critically, it has been shown that the impact of population stratification is also shaped by demographic history which is often overlooked even after adjusting for 100 principal components of common variants4. This study suggests that using principal components derived from rare variants, which capture most of the genetic diversity, or from identity-by-descent (IBD) segments may correct the stratification for some types of effects.

IBD analysis has proven effective for uncovering fine-scale population structure in large cohorts. For example, a study using the UK Biobank5 revealed subtle genetic substructure among British individuals, while research on New York City residents identified 16 ancestry groups, including seven founder populations6, through the All of Us Research Program and the Mount Sinai BioMe biobank. Fine-scale structure has also been documented in the Quebec founder population7–9. These findings raise an important question: to what extent can fine-scale structure be detected in populations? In this study, we address this question by examining a small, recent founder population, the Saguenay–Lac-Saint-Jean (SLSJ) region, leveraging extensive genotype and genealogical data to characterise its genetic architecture.

The SLSJ region of Quebec, Canada (Fig. 1), represents a remarkable model for population genetics research due to its well-documented founder effect, recent settlement and unique demographic history. From the seventeenth century, European settlers of mostly French origin have migrated along the Saint Lawrence River, then founded the Charlevoix region, before settling in Saguenay and Lac-Saint-Jean from 183810. This serial founder effect was accentuated by a 25-fold population surge in SLSJ within a century, mainly due to high birth rates11,12. The region’s demographic characteristics are documented through the BALSAC population register13, a comprehensive database documenting civil records, primarily Catholic marriages, of French-Canadian individuals from the 1660s (New France) to the 1960s (Quebec), encompassing over three centuries of demographic history. This demographic history is known to have led to several variants being more frequent in the region14 and to the increase of some rare diseases15–18. Some researchers have characterised the genetic background of SLSJ as relatively homogeneous16,19, suggesting that, unlike the broader provincial context, the region might lack meaningful structure. Yet, demographic studies indicate that population structure does exist within SLSJ20,21, challenging the assumption of homogeneity.

Previous research in the French-Canadian population has shown that both genotype and genealogical data correlate substantially7–9,22,23. This shows that genealogy may be used as a proxy to study the genetic structure of populations. Yet, despite the wealth of genealogical and demographic data available for the SLSJ region, significant challenges remain in adequately visualising the region’s population structure at a large scale. This limitation primarily stems from algorithms that repeatedly compute kinship for redundant pairs of individuals within a genealogy24, making analysis of complete populations computationally intractable. Although efforts have optimised these computations25,26, implementing these methods at the scale needed for whole-population analysis, particularly for visualization purposes, has remained challenging, restricting samples to hundreds or a few thousand individuals at a time. The absence of appropriate analytical methods has thus limited comprehensive investigation of fine-scale structure in populations.

This study aims to comprehensively investigate the fine-scale population structure within the SLSJ region by integrating genotype and genealogical datasets. We also retrace the population’s origins through the expected genetic contribution of the region’s founders. By doing so, this article reveals a fine-scale population structure shaped by geographical boundaries, differential founders’ genetic contribution and socioeconomic factors in the recent SLSJ founder population.

Results
Comparison between genetic and genealogical structure

The CARTaGENE cohort comprises individuals with both genotype and genealogical data, of which 7,970 probands have sufficient genealogical completeness and are not siblings (see methods). Uniform Manifold Approximation and Projection (UMAP) with two dimensions based on realised and expected kinship (see methods) display strikingly similar population structure patterns (Fig. 2); a linear regression between realised and expected kinship further supports this consistency with a Pearson correlation coefficient of 0.78 (p < 0.001). Hierarchical Density-Based Spatial Clustering of Applications with Noise (HDBSCAN) on a UMAP with ten dimensions identified a subset of 938 individuals most likely originating from the SLSJ region (Supplementary Fig. 1). Within this subset, UMAP projections derived from both realised and expected kinship (Fig. 3) reveal partially overlapping gradients of individuals from the diverse subdivisions of SLSJ (Supplementary Fig. 2) with a significant correlation between realised and expected kinship of 0.83 (p < 0.001).

Fine-scale population structure of a whole generation

The great correlation between the population structure observed through genetic and genealogical data7–9,22,23 allows us to use the genealogy as a proxy for examining the fine-scale structure of the whole SLSJ population of the last generation available in BALSAC. The two-dimensional UMAP projection of non-siblings probands married between 1931 and 1960 in the SLSJ region reveals clear concentrations of individuals whose parents were married in the same subdivision (Fig. 4). The observed structure follows an east–west gradient, extending from the southeast (East of Ha! Ha! subdivision) to the northwest (Ashuapmushuan–Mistassini subdivision) (see Fig. 1). This gradient exhibits a circular pattern centred around individuals originating from Charlevoix and concludes with a distinct grouping of individuals whose ancestry lies outside both Charlevoix and SLSJ. More remote, rural municipalities exhibit higher levels of concentration in genealogical relationships, while urban centres display greater dispersion (Supplementary Fig. 3).

Expected genetic contribution of founders

To refine our understanding of the migratory patterns underlying the observed structure, we calculated the expected genetic contribution of the SLSJ founders (ancestors married in SLSJ whose parents married elsewhere, see methods) to the SLSJ probands of the last generation. Mean genetic contribution from the SLSJ founders originating from Charlevoix is the highest in the eastern part of SLSJ and progressively declines along an east–west axis (Fig. 5a)21. Individuals from the northwesternmost municipalities of Lac-Saint-Jean exhibit more diverse ancestral origins, as well as individuals originating from urban areas, especially south of Saguenay River that exhibit high diversity. A detailed analysis (Fig. 5b and 5c) reveals that this east–west genetic gradient is primarily driven by founders from La Malbaie and to some extent Les Éboulements, which both contributed most significantly to the Saguenay area (see also Supplementary Fig. 4). In contrast, Baie-Saint-Paul contributed more to the Lac-Saint-Jean region (Fig. 5b and 5c, see also Supplementary Fig. 4). Notably, Les Éboulements contributed disproportionately to the east of Ha! Ha! which also presents the highest contribution from Charlevoix founders overall (Supplementary Fig. 5). Individuals from the area between Ashuapmushuan and Mistassini largely trace their ancestry to regions outside Charlevoix (Fig. 5b and 5c), with this diversification becoming prominent around 1885 (Supplementary Fig. 5).

Methods
This study was approved by the ethics board of the Université du Québec à Chicoutimi (UQAC) and complies with all relevant ethical regulations.

Genotype data and cleaning

Genotype data was drawn from the CARTaGENE cohort44, which comprises 43,032 participants aged 40-69 (in 2005) recruited in the main urban areas of the Quebec Province, notably in the city of Saguenay. Of those participants, 29,337 individuals were genotyped using a variety of platforms, including Omni 2.5, GSAv1 with a multi-disease panel, GSAv1, GSAv2 with a multi-disease panel, GSAv3 with a multi-disease panel, GSAv2 with a multi-disease panel and add-on (see https://cartagene.qc.ca/files/documents/other/Info_GeneticData3juillet2023.pdf for details). Each dataset was processed independently using PLINK v1.945, retaining individuals with at least 95% genotyping across all SNPs. At the SNP level, we retained autosomal variants with a call rate of at least 95% across individuals and that conform to Hardy–Weinberg equilibrium (p > 10-6, calculated within each dataset). All chips were merged, and the final dataset comprised 148,200 SNPs across 28,358 individuals.

Genealogical data

Genealogical data was obtained from the BALSAC population register13. Two different genealogical datasets have been extracted for this study: First, 9,405 Quebec residents from the CARTaGENE cohort have been matched to BALSAC records. Second, to analyse a whole generation, we also selected 80,348 probands (i.e. non-parent individuals) married in SLSJ from 1931 to 1960. For both probands’ sets, a multigenerational genealogy extending up to 19 generations was reconstructed. When available, information on the municipality, region, and year of marriage was included; to ensure confidentiality, marriage years are rounded up to the nearest 5. Individual anonymity is rigorously maintained through the use of unique, non-identifying codes.

Traditional geographic subdivisions split SLSJ into Lower Saguenay, Upper Saguenay, and Lac-Saint-Jean20. However, it has been demonstrated that geographic features and natural barriers, such as watersheds, can influence patterns of population structure7. We carefully carved eight subdivisions along primary watercourses to consolidate smaller communities and municipalities were grouped based on these watercourse boundaries (Supplementary Table 1 and Fig. 1). On average, we observed reduced migration in the last generation within our defined geographical subdivisions than between them (Supplementary Fig. 7).

IBD sharing (realised kinship)

Pairwise IBD segments were inferred on phased genotypes using Refined IBD (version 17Jan20)46 and Beagle (version 18May20)47. A matrix of realised kinship estimation was computed by summing shared segment lengths of at least 2 centiMorgans (cM) using a Python 3 ported version of relatedness_v1.py from Sharon R. Browning (see https://faculty.washington.edu/sguy/ibd_relatedness.html). The diagonal was set to 1.

Genealogical kinship coefficients (expected kinship)

Expected kinship (\phi)48 was used to infer and visualise population structure through genealogical data. It corresponds to the probability that, at one locus, one randomly picked allele from individual i and one randomly picked allele from individual j are IBD or come from the same ancestor. Expected kinship coefficients were computed between all pairs of probands from whole SLSJ and CARTaGENE genealogical datasets. Of the 80,348 SLSJ probands, 26,445 are non-siblings (\phi > 0.2) and have a genealogy deemed sufficiently complete, i.e. they are related (\phi > 0) with at least half of the probands. Of the 9,405 individuals in CARTaGENE data, 7,970 are neither parents nor siblings, and have a genealogy deemed sufficiently complete using the same criterion.

Hybrid algorithm for time-efficient computation of billions of kinship coefficients

To efficiently compute kinship coefficients across all pairs of individuals, we implemented a hybrid of two existing algorithms. Karigl’s algorithm24 estimates kinship coefficients on a per-proband-pair basis, enabling parallel processing of the kinship matrix, as in the GENLIB package for R49. However, it reprocesses each ancestor every time they appear, either multiple times across genealogies (kinship) or repeatedly within the same genealogy (inbreeding), which makes it computationally costly. Kirkpatrick’s algorithm26 addresses this inefficiency by partitioning the genealogy into successive sub-genealogies (e.g., one generation at a time) and computing kinship coefficients in a top-down manner, with the probands of one sub-genealogy treated as the founders of the next, thereby minimizing redundant ancestor processing. Yet, Kirkpatrick’s approach is not designed for parallel execution. Our method combines the strengths of both: we divide the genealogy into generational sub-genealogies, as in Kirkpatrick’s approach, but compute each sub-genealogy’s kinship matrix in parallel, achieving substantial time gains. For instance, computing kinship for all 6,455,801,104 pairs of 1931–1960 BALSAC SLSJ probands required only two minutes with our algorithm. It was not possible to compare directly with Kirkpatrick’s implementation due to the fact that this software is no longer available online, nor with the GENLIB implementation because its execution would require more resources than can be allocated. Following the pre-publication of this work, it was brought to our attention that an indirect method by Colleau50 can also be applied to compute kinship coefficients at a large scale51. However, this method incorporates all generations in the genealogy, which was not the objective of our analysis.

Our implemented hybrid algorithm is provided as part of GeneaKit 0.1.0, which is freely available on GitHub52 (https://github.com/Genopop/geneakit) under the MIT license, using the phi() function for kinship coefficients. A pseudocode of the algorithm is available in the supplementary material (Supplementary Note 1).

Expected genetic contributions

We also computed the expected genetic contribution of the region’s founders as a way to measure the region’s migratory origins and to better understand the fine structure. Region’s founders are defined as individuals who married in the SLSJ whereas their parents married elsewhere. This computation was done using GeneaKit 0.1.0’s gc() function, from all regions’ founders to all 80,348 genealogical probands who married between 1931 and 1960 in the SLSJ region. The expected genetic contribution of an ancestor corresponds to the probability that an allele is passed down to a given descendant. A founder’s region of origin is defined as the region of marriage of their parents.

Clustering and visualisation

CARTaGENE individuals originate from all regions of Quebec, it is thus necessary to identify those who are most likely to originate from SLSJ. UMAP (version 0.5.9)42 was performed on both the realised and expected kinship matrices, using precomputed distances of (1 – \phi). UMAP was run with ten dimensions, 15 neighbors, a minimum distance of 0.0, and a fixed random seed (random state = 0), to maximise clustering and reproducibility. UMAP was initialised from a classical multidimensional scaling (MDS), also known as principal coordinate analysis (PCoA), of the precomputed distances, implemented as pcoa from scikit-bio 0.6.353, using eigendecomposition and a seed of 0. HDBSCAN from scikit-learn 1.8.054 was then fitted to the UMAP outputs to identify individuals most likely originating from the SLSJ region, using a min_cluster_size of 25 and a cluster_selection_epsilon of 0.3. Those parameters were chosen as they allowed separation of regional populations in a previous study using the same cohort44. From the clustering results for both the realised and expected kinship matrices, the cluster containing the highest number of individuals of known SLSJ origin (i.e. whose four grandparents married in the region) was identified. 585/1141 individuals in the realised kinship cluster and 583/1036 individuals in the expected kinship cluster are of known SLSJ origin. The final set of 938 individuals was then determined as the intersection of these two best clusters (one from the realised kinship matrix clustering and one from the expected kinship matrix clustering). Of those 938 individuals, 566 are of known SLSJ origin.

For visualisation, UMAP (initialised with classical MDS/PCoA) was applied to the kinship matrices using their precomputed distances (1 – \phi), two dimensions, and a minimum distance of 0.5.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Morin et al, analyze the fine scale structure of the Saguenay–Lac-Saint-Jean (SLSJ) region of Quebec, known for its high frequency of some rare genetic disease. Although the founder events and demographic processes that lead to the increase of frequency of these disorders have been extensively studied, the fine scale genetic structure is less known. This is important if subregion have a different genetic makeup for genetic counseling. Further it is now demonstrated that substructure can lead to biases estimation in Gwas studies.\nTheir study use the extensive demographic genealogical database Balsac that covers the whole region. More precisely, they study the genealogies of 80 348 probands married between 1921 and 1960. And they use the genetic database the CARTaGENE cohort that includes whole genome genotyping of 29,337 individuals, among them 9405 matched in the Balsac genealogical database.\nTheir study show a clear genetic substructure. This is clearly important since this population has always been described as quite homogenous.\nThey first demonstrate a high correlation for these 9405 individuals between genetic and genealogical kinship. This result enables them to use the whole demographic database to infer genetic substructure in the region. In order to achieve these results, they develop a new method to compute kinship on their big genealogical dataset.\nTheir results show a clear fine scale structure.\nThe paper is well written, the methodology is impressive.\nMy remarks are minors:\n1. some maps would be usefull.\n- a map of the East of Quebec province in the main text would be usefull.\n- maps that present the substructure differences in term of founder genetic contribution\n- map that present the kinship level.\n2. It could be useful to address a little bit more the importance of substructure for GWAS studies", "reviewer_2": "Reviewer comment for\nFine-scale structure of a whole regional population through genetics and genealogies\nSummary: The authors investigate the fine-scale population structure of a French Canadian population by computing two kinship metrics, one expected (from a pedigree) and the other realized. They convert the kinship matrix into a distance matrix, then apply non-linear dimension reduction and clustering algorithms to obtain smaller subpopulations from the data. In particular, they developed a novel kinship coefficient computing algorithm that offers substantial improvement over existing methods. Although I think the novel algorithm is helpful, I’m not convinced by the clustering/dimension reduction results. Below are my comments.\nMajor comments\n1. Can the authors try out PC initialization? Unlike other branches of genomics where non-linear dimension reduction (hereafter, DR) is widely used, literature in applying non-linear DR to DNA sequence data has been relatively sparse. A few examples I can remember are the All-of-Us paper and another Quebec paper that used a dataset that overlaps with the paper currently being reviewed. Therefore, I believe there is little guidance on how to properly apply non-linear DR or if it’s even appropriate to apply these types of DRs. At this point, it would be informative to validate the wisdom from other fields in genomics, such as single-cell RNA sequencing. It is more common to apply PCA first and then feed the PC coordinates to non-linear DR as described in this paper (https://www.nature.com/articles/s41587-020-00809-z). Another paper working with DNA sequence data takes similar steps (https://www.cell.com/ajhg/pdfExtended/S0002-9297(22)00112-4), PCA first and then UMAP. It seems that initializing the non-linear DR with PCA makes a substantial difference. I would perhaps do the following two analyses:\na) Initialize with PCA, then feed PCA to the non-linear DR\nb) Initialize with PCA, but feed the original distance matrix to non-linear DR. i.e., do the exact same analysis as in the draft but add PC initialization.\nFor pedigree PC, there are several options including the paper in which the authors cited (https://gsejournal.biomedcentral.com/articles/10.1186/s12711-025-00994-y).\n2. How robust is the HDBSCAN/UMAP based classification? Visually, the points are very crowded and putting lines between them feels very artificial. One could try to adjust the parameters to check robustness, which the authors have already done to some extent, but we have a better alternative in this case. Since the migration patterns of Quebecois are already documented to a certain extent, it will be possible to run a population genetic simulation using msprime or SLiM (https://tskit.dev/) (or perhaps another suitable simulation scheme) . Then, authors could apply their pipeline to the simulated data to see if it recovers the demographic pattern codified into the simulation. The resulting tree sequence is amenable for kinship extraction (e.g. https://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1011537) which should be amenable to the paper’s pipeline.\n3. What is the time scaling of the new algorithm? Since it has to materialize an N-by-N kinship matrix, it should be at least O(N^2). I wonder if this is indeed the case.", "reviewer_3": "In this work Morine et al investigate the fine population structure of the Saguenay-Lac-Saint-Jean region of Quebec (Canada) integrating extensive genotype data and genealogical records. To do this they optimize an algorithm to efficiently compute kinship coefficients from genealogies. The study is thorough, methodologically sound, and well written; however, I have some reservations regarding its suitability for publication in Nature Communications. Indeed, while the study provides a detailed and quite convincing analysis of the Saguenay–Lac-Saint-Jean population (see below for specific comments on the presented analyses), the topic may be too specialized for the broad readership of Nature Communications. Moreover, the extent to which the proposed approach can be generalized remains unclear. The authors benefit from exceptionally rich datasets, both in terms of genotype information and genealogical records, which may not be available in other cases. Consequently, the broader applicability of their findings and methods might be limited. It is also worth noting that the authors do not present new data, but rather reanalyze previously generated genotype and genealogical datasets. While the methodological optimization is valuable, the lack of novel data somewhat limits the overall originality and impact of the study.\nDespite my view that this manuscript may not be appropriate for publication in Nature Communications, I include below some remarks and suggestions concerning the analyses carried out, specifically on the assessment of genetic structure and dimensionality reduction methods.\nThe authors investigate the genetic structure of the population through IBD sharing (realized kinship), aiming to efficiently compare genetic and genealogical relationships. It would, however, be informative to also explore the genetic structure directly from genotype data (e.g. SNP-based analyses), for instance through principal component analysis (PCA) or estimates of FST between clusters. Such analyses could provide complementary insights into the fine-scale population structure shaped by watercourses and the east-west migration patterns described in the study. In addition, the authors employ UMAP projections to summarize both genetic and genealogical data. While UMAP is a powerful dimensionality reduction method, it can sometimes introduce distortions that exaggerate the apparent separation between clusters. I therefore recommend that the authors also perform a PCA and subsequent clustering on the PCA coordinates (e.g. using k-means clustering) to verify whether the genetic structure observed in the UMAP plots is robust.\nOverall, this is a carefully executed study that makes good use of both genetic and genealogical data to characterize the fine-scale population structure of the SLSJ region. The methodological work is solid, and the manuscript is clearly written. However, the study mainly reanalyzes existing datasets and focuses on a highly specific regional case. As a result, the broader generalizability and potential impact of the findings appear limited. While the methodological contribution and the integration of rich genealogical resources is appreciable, I am not convinced that the current version of the manuscript meets the breadth and impact expected for publication in Nature Communications. I would, however, encourage the authors to consider submission to a more specialized genetics or population genomics journal.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "We appreciate that the editor sent our work for peer review and are pleased that the reviewers recognised its technical soundness and the value of integrating genealogical and genetic data with an optimised kinship algorithm. We respectfully request an appeal of the decision, as we believe we can address the reviewers’ concerns regarding the broader applicability of our approach.\nGeneral comments:\nIn particular, Reviewer 3 emphasised that the method has only been demonstrated in a uniquely rich, region-specific case with exceptional depth of genealogical and genotype data, making it difficult to judge whether the approach generalises beyond Saguenay-Lac-Saint-Jean.\nWe have substantially revised the manuscript to present the work with broader generalisation and relevance for populations with less comprehensive data resources. We believe that access to such rich data offers a unique opportunity to understand the complexity of fine structure in human populations. We show in our work that we can detect fine structure within what was considered a ‘homogeneous’ and very recent founder population like SLSJ. It is reasonable to expect that similar recent structure exists in all populations, signals that are often obscured by coarse, ancient structure under the standard stratification correction methods used in most GWAS and polygenic risk score analyses. Relying solely on PCA of common SNPs has been shown to leave residual stratification, which can lead to serious biases and false associations. Our study now emphasises why addressing this issue is critical. In addition to the revised scope of the manuscript asked by reviewer #3, we made the following changes:\n● Add maps of eastern Quebec and genetic contribution as requested by reviewer #1.\n● Replace UMAPs of expected and realised kinship matrices by UMAPs initialised on MDS as requested by reviewer #2.\n● Replace supplementary Fig. 3 to reflect the 10-dimensional UMAP used to identify the SLSJ cluster as raised by reviewer #2.\nWe have addressed these points and all others in our following detailed responses to the reviewers. We hope that the revised manuscript will now be considered of broad interest to the Nature Communications readership.\nAdditionally, Reviewers 2 and 3 highlighted the need for more robust validation of the clustering and dimensionality reduction strategy (e.g., PCA initialisation, simulation-based validation, and orthogonal assessments of structure from genotype-derived SNP analyses).\nWe have addressed these points in detail in our responses to the reviewers and implemented the corresponding changes in the manuscript. Please refer to the specific point-by-point answers for further clarification.\nReviewers' comments:\nReviewer #1 (Remarks to the Author):\nTheir study show a clear genetic substructure. This is clearly important since this population has always been described as quite homogenous.\nTheir results show a clear fine scale structure.\nThe paper is well written, the methodology is impressive.\nMy remarks are minors:\n1. some maps would be usefull.\n- a map of the East of Quebec province in the main text would be usefull.\nWe added this as Fig. 1 in the main text. During this process, we identified that St-Prime had been misclassified and corrected this across all figures and results.\n- maps that present the substructure differences in term of founder genetic contribution\nWe also added this as Fig. 5C in the main text.\n- map that present the kinship level.\nBecause kinship is inherently pairwise, we explored several visualisation approaches to represent it effectively. However, none of the configurations provided a clear or interpretable view without appearing cluttered and confusing. Therefore, we decided not to include this map in the revised manuscript.\n2. It could be useful to address a little bit more the importance of substructure for GWAS studies\nWe have incorporated this into the main text, and the scope of the work is now significantly broader, as also requested by reviewer #3. We sincerely thank the reviewer for this valuable comment, which has greatly improved our manuscript.\nReviewer #2 (Remarks to the Author):\nMajor comments\n1. Can the authors try out PC initialization? Unlike other branches of genomics where non-linear dimension reduction (hereafter, DR) is widely used, literature in applying non-linear DR to DNA sequence data has been relatively sparse. A few examples I can remember are the All-of-Us paper and another Quebec paper that used a dataset that overlaps with the paper currently being reviewed. Therefore, I believe there is little guidance on how to properly apply non-linear DR or if it’s even appropriate to apply these types of DRs. At this point, it would be informative to validate the wisdom from other fields in genomics, such as single-cell RNA sequencing. It is more common to apply PCA first and then feed the PC coordinates to non-linear DR as described in this paper (https://www.nature.com/articles/s41587-020-00809-z). Another paper working with DNA sequence data takes similar steps (https://www.cell.com/ajhg/pdfExtended/S0002-9297(22)00112-4), PCA first and then UMAP. It seems that initializing the non-linear DR with PCA makes a substantial difference.\nWe appreciate the reviewer’s suggestion and apologise for the lack of clarity in our methods. The reviewer is correct that initialisation plays an important role in UMAP. While our original analysis used spectral embedding (the default in the umap-learn package) for initialisation, we have now applied MDS initialisation as suggested. Because IBD represents pairwise similarity rather than raw genotype features, MDS provides a more appropriate initialisation than PCA, as it is applied to pairwise relationships directly. This approach is consistent with previous work that applied UMAP to MDS of IBD matrices1. UMAP figures were changed accordingly in the main text.\nI would perhaps do the following two analyses:\na) Initialize with PCA (MDS), then feed PCA (MDS) to the non-linear DR:\nThis analysis was performed using varying numbers of MDS components (see Fig. 1 below). To fully capture the population structure represented by the expected kinship matrix, approximately 1,000 MDS components were required. The first components primarily reflect the Charlevoix gradient and the East of Ha! Ha! cluster, which explain most of the variance. Consequently, when too few components are used, fine-scale structure remains masked behind these dominant signals.\nFig. 1: UMAPs of varying numbers of dimensions of an MDS calculated on genealogical kinship of individuals with parents married in SLSJ in 1931–1960, initialised on MDS 1 and 2.\nb) Initialize with PCA (MDS), but feed the original distance matrix to non-linear DR. i.e., do the exact same analysis as in the draft but add PC initialization.\nThe resulting UMAP is very similar to our previous version initialised with spectral embedding (see Fig. 2 below). Therefore, we included the new figure initialised with MDS in the manuscript, as suggested by Reviewer #2. Considering this and our focus on capturing recent fine-scale structure, which appears easier when using the kinship matrix with UMAP initialised by the first two MDS dimensions, we decided to retain this approach in our manuscript.\nFig. 2: UMAP of pairwise expected kinship matrix of individuals with parents married in SLSJ in 1931–1960, initialised on MDS 1 and 2.\nFor pedigree PC, there are several options including the paper in which the authors cited (https://gsejournal.biomedcentral.com/articles/10.1186/s12711-025-00994-y).\nWe did not use randPedPCA, as it appears to incorporate all ancestors in the genealogy and to add a genealogical depth in the response PCA, which was beyond the scope of this study. Our focus was on comparing contemporary probands, i.e. the last generation in the genealogy, based on their genealogical and genotype data.\n2. How robust is the HDBSCAN/UMAP based classification? Visually, the points are very crowded and putting lines between them feels very artificial. One could try to adjust the parameters to check robustness, which the authors have already done to some extent, but we have a better alternative in this case. Since the migration patterns of Quebecois are already documented to a certain extent, it will be possible to run a population genetic simulation using msprime or SLiM (https://tskit.dev/) (or perhaps another suitable simulation scheme) . Then, authors could apply their pipeline to the simulated data to see if it recovers the demographic pattern codified into the simulation. The resulting tree sequence is amenable for kinship extraction (e.g. https://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1011537) which should be amenable to the paper’s pipeline.\nThe reviewer is right that the HDBSCAN clustering feels artificial and can’t be visually validated on our actual Supplementary Fig. 3. We apologise for the confusion caused by this figure. It gives the impression that the regions are very close to each other, but this is due to the use of two-dimensional UMAPs for visualisation purposes. In reality, the clustering was performed using UMAPs with ten dimensions combined with HDBSCAN, this was stated in the text (line 161). When visualising these ten dimensions, it becomes clear that the SLSJ region is indeed distinct from the other regions in both our IBD data (realised kinship) and genealogical data (expected kinship). We have replaced Supplementary Fig. 3 with a visualisation showing all ten UMAP dimensions used for the HDBSCAN clustering. As stated in the text, the SLSJ cluster was chosen as the cluster containing the highest number of individuals of known SLSJ origin i.e. with grandparents married in SLSJ.\nPlease also refer to our response to reviewer #3 for a related question, and see Fig. 3 below, which illustrates the overlap between the HDBSCAN cluster based on the 10-dimensional UMAP and the SLSJ region identified through PCA of common SNPs.\n3. What is the time scaling of the new algorithm? Since it has to materialize an N-by-N kinship matrix, it should be at least O(N^2). I wonder if this is indeed the case.\nWe have added a paragraph under the pseudocode in the supplementary materials to describe the time complexity. It is the same as for Kirkpatrick’s recursive-cut exact algorithm (https://doi.org/10.1093/bioinformatics/bty725), i.e. O(N2G) where N is the number of individuals in the generation and G is the number of generations. Our improvement lies in the parallelisation of this process.\nWhile reassessing the time complexity of our hybrid algorithm, we noticed that by sorting the individuals during the preprocessing, the production of vertex cuts had a time complexity of O(N log N) instead of O(N) as we expected, where N is the number of individuals in the genealogy. Since the ordering of the individuals does not matter within the intermediate generations, we have removed this sorting from our pseudocode and from the GeneaKit implementation.\nReviewer #3 (Remarks to the Author):\nThe study is thorough, methodologically sound, and well written; however, I have some reservations regarding its suitability for publication in Nature Communications. Indeed, while the study provides a detailed and quite convincing analysis of the Saguenay–Lac-Saint-Jean population (see below for specific comments on the presented analyses), the topic may be too specialized for the broad readership of Nature Communications. Moreover, the extent to which the proposed approach can be generalized remains unclear. The authors benefit from exceptionally rich datasets, both in terms of genotype information and genealogical records, which may not be available in other cases. Consequently, the broader applicability of their findings and methods might be limited.\nIt is also worth noting that the authors do not present new data, but rather reanalyze previously generated genotype and genealogical datasets. While the methodological optimization is valuable, the lack of novel data somewhat limits the overall originality and impact of the study.\nWe respectfully disagree with the reviewer that our work is a reanalysis of existing datasets. To our knowledge, this is the first comprehensive characterization of fine-scale structure across the entire SLSJ population, made possible by our exceptionally rich genealogical dataset. Previous studies were limited to subsets of the region due to the absence of tools capable of handling such large-scale genealogical data. Although we use the same genotype Quebec dataset as the McClelland et al. study2 (recently accepted in Nature Communications), our objectives and contributions are different. The McClelland study highlights the broad diversity of Quebec and the opportunities offered by the biobank for health research, whereas our work uncovers the fine-scale genetic structure of the contemporary SLSJ population, a population shaped by a strong founder effect and long assumed to be homogeneous. Our work is therefore complementary, not redundant.\nMoreover, we believe that the richness of Quebec’s data is not a limitation but a strength: it serves as a model for understanding genetic fine structure linked to recent demographic history, which occurs not only in Quebec, but possibly in many populations worldwide3,4. Therefore, the standard stratification correction strategies used in most GWAS and polygenic risk score analyses might not correct for this fine structure which can lead to biases and false associations.\nOur goal is not to claim universal applicability of our method, but to demonstrate that fine-scale structure exists even at very small geographic scales and that accounting for it is relevant and necessary for genetic studies. Nevertheless, our approach can be applied to other populations wherever genotype data are available. Since our results are consistent on genealogical kinship and genomic IBD, this allows for the generalisation in any biobank or single population at low costs since only genotyping is required.\nWhat is new in our study:\n● Kinship coefficient calculation across the entire SLSJ population, described at a level of fine-scale structure that surpasses any previous study within Quebec’s regional divisions.\n● This was recently made possible to calculate pairwise kinship on billions of pairs of individuals with the new tools available, including the one that we present here, GeneaKit.\n● Detecting such fine-scale structure within a relatively small and ‘homogeneous’ population suggests that fine structure exists in all populations. It underscores the need to control for population stratification in all GWAS and polygenic analyses, requiring adjustment beyond the first 20 principal components generally used and the continental population structure5.\nThe authors investigate the genetic structure of the population through IBD sharing (realized kinship), aiming to efficiently compare genetic and genealogical relationships. It would, however, be informative to also explore the genetic structure directly from genotype data (e.g. SNP-based analyses), for instance through principal component analysis (PCA) or estimates of FST between clusters. Such analyses could provide complementary insights into the fine-scale population structure shaped by watercourses and the east-west migration patterns described in the study.\nWe do agree with the reviewer #3 that SNP-based analyses such as PCA and F have been ST widely used to infer population structure. However, these methods excel at capturing ancient divergence and broad continental patterns, but often lack the resolution needed to detect subtle, recent demographic events such as fine-scale structure. This limitation arises because differences in SNP allele frequency accumulate slowly over time, making them less sensitive to recent demographic events6. In contrast, IBD has been shown to directly reflect recent shared ancestry, with segment length correlating to the number of generations since a common ancestor. This property enables IBD-based analyses to uncover finer geographic and genealogical structure within the last few hundred years, as demonstrated in other large-scale studies on the UK Biobank4 and on New York City residents from the All of Us cohort3. Indeed, it was shown that when population structure is recent, it cannot be corrected using principal components of common variants because they are uninformative about recent history5. This can result in biases in GWAS, and even stronger biases in polygenic risks assessments.\nConsidering all these factors, we decided to continue using IBD so that we could apply the same methodology consistently across both genealogical and genotype data types.\nIn addition, the authors employ UMAP projections to summarize both genetic and genealogical data. While UMAP is a powerful dimensionality reduction method, it can sometimes introduce distortions that exaggerate the apparent separation between clusters. I therefore recommend that the authors also perform a PCA and subsequent clustering on the PCA coordinates (e.g. using k-means clustering) to verify whether the genetic structure observed in the UMAP plots is robust.\nThe question on the reliability of UMAPs has been addressed above, please refer to reviewer #2 answers to questions 1 and 2 for details.\nRegarding clustering: We apologise for the confusion caused by our original Supplementary Fig. 3, which did not accurately represent the analysis. The figure suggested that regions were very close, but this was an artifact of using two-dimensional UMAP for visualisation. In reality, clustering was performed with UMAP using 10 dimensions combined with HDBSCAN. We now provide an updated Supplementary Fig. 3 showing all 10 UMAP dimensions.\nTo illustrate the differences between HDBSCAN clustering on UMAP and the proposed k-means clustering on PCA of SNPs, Fig. 3 below compares: PCA of common SNPs colored by parents’ marriage region (top), k-means clustering on this PCA (k = 3 since this PCA clearly distinguishes a northeastern cluster (SLSJ) from an eastern cluster (Bas-St-Laurent and Côte-du-Sud)) (middle), and HDBSCAN clustering on the 10-dimensional UMAP derived from realised and expected kinship matrices (bottom).The SLSJ cluster (shown in red) is clearly visible along PC1, consistent with Quebec’s previously described genetic structure based on common SNPs and genealogical data7–9. The agreement between methods is strong: 75% of the 938 individuals identified by HDBSCAN were also captured by k-means clustering on the PCA of common SNPs. This PCA in Fig. 3 below likely does not provide an optimal shape for k-means, as it fails to capture the full continuum of SLSJ ancestry. In contrast, the HDBSCAN clustering consists almost entirely of individuals located along PC1 in the region with the highest concentration of parental marriages within the SLSJ.\nFig. 3: PCA calculated on common SNPs of more than 7,000 Quebec individuals with genotype and genealogical data colored according to parents’ region of marriage (top); K-means clustering (middle); and HDBSCAN clustering on a 10-dimensional UMAP (bottom).\nOverall, this is a carefully executed study that makes good use of both genetic and genealogical data to characterize the fine-scale population structure of the SLSJ region. The methodological work is solid, and the manuscript is clearly written. However, the study mainly reanalyzes existing datasets and focuses on a highly specific regional case. As a result, the broader generalizability and potential impact of the findings appear limited. While the methodological contribution and the integration of rich genealogical resources is appreciable, I am not convinced that the current version of the manuscript meets the breadth and impact expected for publication in Nature Communications. I would, however, encourage the authors to consider submission to a more specialized genetics or population genomics journal.\nWe thank the reviewer for this valuable comment. We have substantially revised the manuscript to present the work with broader generalisation and relevance for populations with less comprehensive data resources. We believe that access to such rich data offers a unique opportunity to understand the complexity of fine-scale structure in human populations globally. Indeed, if we can detect fine-scale structure within what was considered a ‘homogeneous’ founder population like SLSJ, it is reasonable to expect that similar recent signals exist in all populations, signals that are often obscured by coarse, ancient structure under the standard stratification correction methods used in most GWAS and polygenic risk score analyses5. Relying solely on PCA of common SNPs has been shown to leave residual stratification, which can lead to serious biases and false associations. Our study emphasises why addressing this issue is critical.\nReviewer #4 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nWe thank the reviewer #4 for their valuable comments."}, {"round": 2, "reviewer_comments": {"reviewer_2": "The authors have made substantial efforts to revise their manuscript and I only have one major concern left. It's an overreaching claim to say that any magnitude of population structure can induce spurious associations in association studies. It is in fact the realized allele frequency stratification that creates excess false positives. For instance, Edge et al. is a good reference to this point (https://pubmed.ncbi.nlm.nih.gov/24481204/). Also, see Rosenberg et al. (https://academic.oup.com/genetics/article/173/3/1665/6060953).\nIf the authors want to keep up with this claim, they should run simulations to confirm their findings: since they have a pedigree, simulate genomes conditioned on the pedigree and run GWAS on them. I'll recommend using msprime (https://tskit.dev/msprime/docs/stable/pedigrees.html), which can easily generate genomes conditioned on the pedigree. Then, using the genotypes from the data, one can further simulate traits and conduct GWAS to check the claim.", "reviewer_3": "The authors have provided convincing arguments addressing the most delicate point raised in my previous report, namely the specificity of the results and their limited generalizability. Their response is detailed and, this time, persuasive. In light of this, considering the opinions of the other referees and the revisions introduced in the new version of the manuscript, I now believe that the work is appropriate for publication in Nature Communications.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "We appreciate that the editor resubmitted our work for peer review and are pleased that the reviewers recognised the broader applicability of our approach through the updated manuscript.\nReviewers’ comments:\nReviewer #2 (Remarks to the Author):\nThe authors have made substantial efforts to revise their manuscript and I only have one major concern left. It's an overreaching claim to say that any magnitude of population structure can induce spurious associations in association studies. It is in fact the realized allele frequency stratification that creates excess false positives. For instance, Edge et al. is a good reference to this point (https://pubmed.ncbi.nlm.nih.gov/24481204/). Also, see Rosenberg et al. (https://academic.oup.com/genetics/article/173/3/1665/6060953).\nWe thank the reviewer for their comment. As rightly pointed out, we have adjusted the manuscript to clarify that population structure may lead to spurious associations only through correlation between allele frequencies and phenotypic variation.\nIf the authors want to keep up with this claim, they should run simulations to confirm their findings: since they have a pedigree, simulate genomes conditioned on the pedigree and run GWAS on them. I'll recommend using msprime (https://tskit.dev/msprime/docs/stable/pedigrees.html), which can easily generate genomes conditioned on the pedigree. Then, using the genotypes from the data, one can further simulate traits and conduct GWAS to check the claim.\nWe recognise the value of performing simulations to confirm the presence of population stratification when running GWAS on the SLSJ data. We have added in the manuscript Supplementary Fig. 7 (copied below) and its corresponding Supplementary Methods and References. We first simulated genomes using msprime conditioned on the genealogy for all unrelated SLSJ individuals (above first cousins). Knowing that the correlation between our PCoA 1 of expected kinship coefficients and expected genetic contributions (GC) from Charlevoix founders is high (Pearson’s r = 0.97, p < 0.001), we simulated a continuous trait based on the Charlevoix GC using this formula:\n2 2\n𝑦 ~ 𝐺𝐶 × ℎ + 𝑁(0,1) × 1 − ℎ ,\nwhere y is the continuous trait, GC is the normalized Charlevoix GC, h2 is the simulated heritability (0.1), and N(0, 1) is a random noise of normal distribution.\nWe then ran a GWAS and observed inflation due to stratification on the QQ plot when no correction was applied. These results demonstrate that standard PCA correction may be insufficient for GWAS in populations with fine-scale structure like SLSJ, if a trait happens to be correlated with this fine structure.\nMoreover, since SLSJ’s fine-scale structure is highly predicted by its Charlevoix ancestry, and this Charlevoix ancestry has been shown to highly correlate with the presence of rare diseases in the region, we think it is reasonable to conclude that some phenotypic variations may correlate with SLSJ’s population structure and cause spurious associations if left uncorrected, as shown in our simulation.\nSupplementary Fig. 7: PCA does not adequately control for inflation of GWAS on a simulated phenotype correlated with SLSJ fine-scale structure. QQ plots for simulated data (N = 10,000) across 100 Mb of neutral genomic sequence using msprime conditioned on the genealogy. The phenotype was simulated to be partially correlated with the expected genetic contribution (GC, computed on the genealogy) from Charlevoix founders (h2 = 0.1) since this GC was shown to be a major factor influencing the SLSJ fine structure. Panels display association results for: (A) a baseline null phenotype (pure noise); (B) the correlated phenotype with no correction; (C–D) the correlated phenotype corrected using the top 10 and 100 genetic principal components (PCA) of common SNPs (MAF < 0.05); and (E) the correlated phenotype corrected using the top 2 coordinates from a PCoA derived from the genealogical kinship matrix. The genomic inflation factor (λ ) is reported for each model to quantify test GC statistic inflation. The red dashed line represents the null expectation (y = x), and the gray shaded region indicates the 95% confidence interval.\nReviewer #3 (Remarks to the Author):\nThe authors have provided convincing arguments addressing the most delicate point raised in my previous report, namely the specificity of the results and their limited generalizability. Their response is detailed and, this time, persuasive. In light of this, considering the opinions of the other referees and the revisions introduced in the new version of the manuscript, I now believe that the work is appropriate for publication in Nature Communications.\nWe thank the reviewer #3 for their valuable contribution that improved our manuscript.\nReviewer #4 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nWe thank the reviewer #4 for their valuable contribution."}, {"round": 3, "reviewer_comments": {"reviewer_2": "The new benchmark looks excellent. I do not have any concerns left. It's a great paper, and best to my knowledge, one of the few paper that clearly demonstrates the significance of regional population structure in human populations."}, "author_response": ""}]

Decision letter