We have collected samples from five reference hospital centers belonging to three health regions expanding the catalan Pyrenees. Recruitment period was extended from July 2014 to October 2016. In total, we recruited 435 samples (51.3% females) of subjects born and with all their grandparents born in the Pyrenean region (see insert Figure 1).

The data set represented an oldest extract of the population, with an average age of 55.6 years and equal proportion of sexes (Table 1). Therefore, this sample should be considered moderately affected by the demographic changes that occurred during the 19th and 20th centuries. All subjects signed an informed consent, and the study had the approval of the Ethics Committee of the University of Barcelona.

For each participant, we picked up demographic data such as sex, age and place of birth of the participant and place of birth their father and mother and for all their four grandparents. Birth places were assigned to geographic coordinates using Vissir3, the official viewer of the Cartographic and Geological Institute of Catalonia (Vidal, 2012). Assignment of centroid location of participants was done by calculating the mean longitude and latitude of grandparents (see Supplementary Table S1). Regional cartographic base maps and the terrain elevation model were from the Cartographic and Geological Institute of Catalonia. Maps were made with ArcGis 8.1 software (ESRI, Redlands, CA: Environmental Systems Research Institute), using facilities at the service of Cartography & SIG, of the University of Lleida (UdL).

We collected 1 ml of peripheral blood from each participant in EDTA-coated tubes. All participants were genotyped using Applied Biosystems Axiom™ Precision Medicine Diversity Array (PMDA) on GeneTitan™ Multi-Channel Instrument (Thermo Fisher, CA, United States) at Centro Nacional de Genotipado facilities (CEGEN, nodo Santiago de Compostela, Spain). Genotype clustering was analyzed using the Axiom Analysis Suite 6.0 software adjusting for BestWorkflow settings. Quality control (QC) was tuned to sample calling rate >0.95, marker calling rate >0.90, and Hardy-Weinberg equilibrium (HWE) >10^–5. After removing Indel and copy number variation (CNV), and QC, 780,310 SNPs (726,718 at autosomes, 49,319 at chromosome X, 406 at chromosome Y, 3,457 at pseudoautosomal region and 410 at mitochondrial genome) on 435 subjects were available for subsequent analyses (thereafter GENPIR data set).

We estimated relatedness among samples by using SNPRelate R package (Zheng et al., 2012). Identical by descent (IBD) analysis reveals a degree of relatedness among several individuals (Supplementary Figure S3A). In order to consider unrelated individuals for further analyses, we selected one individual from each pair of individuals having a kinship ≥0.09 (Supplementary Figure S3B). This yields a final set of 397 unrelated individuals.

Data was phased using the Michigan-imputation server-1.5.7 pipeline (Das et al., 2016) with the eagle-2.4 phasing (Loh et al., 2016), using TopMed reference panel (Version R2 on GRC38) as reference.

IBD tracts between pairs of chromosomes from different individuals and homozygote by descent (HBD) fragments between chromosomes from the same individual were inferred with Beagle 4.0 (IBD_Beagle) (Browning and Browning, 2013), using the human recombination map provided by the same software (Browning and Browning, 2013).

To make comparisons with other populations, GENPIR genotype data was merged with genotype data from samples of the 1,000 Genomes Phase 3 project (1KG3) (2,504 individuals) (Auton et al., 2015) and Human Genome Diversity Panel (HGDP) (929 individuals) (Bergström et al., 2020). After merging with the GENPIR dataset, a total of 3,868 individuals and 738,061 single nucleotide variants (SNVs) were selected with a genotyping rate of 0.989 (690830 at autosomes, 46,947 at chromosome X, 280 at chromosome Y). We also applied unrelatedness criteria to the merged data set resulting in a final worldwide set of 3,805 unrelated individuals that were used for further analyses.

We first analyzed the genetic relationships of the sampled individuals at worldwide level using a Principal Component Analysis (PCA) (Figure 2, Supplementary Figure S4).

The first two PCA components (PCs), explaining 11.3% of variance, show that the considered catalan Pyrenees samples are genetically related to other samples with European ancestry (Figure 2A). Within Europe, the catalan Pyrenees samples fall within the first two PCs (explaining 0.7%) close to samples from Southern Europe, close to individuals from Iberian and French Basque origin (Figure 2B). Scatterplots of the top 4 PCs with proportional variance explained included are also presented in Supplementary Figure S4.

Given these results, we run an unsupervised ADMIXTURE analysis focusing on the European region for K ranging from one to four in order to quantify the European ancestry proportions in GENPIR samples, revealing a minimum cross validation error value for K= 2–3 (Supplementary Figure S5).

As shown in Figure 3A, for K = 2 a North (red) to South (blue) clinal component was observed that was extended to an additional Central (green) component at K = 3. Finally, for K = 4, Southern component was split into a Southeastern (orange) and Southwestern (blue) component. Considering these four major components maximizing in present-day Southwestern European populations as “Northern (N)”, “Southeastern (SE)”, “Southwestern (SW)” and “Central (C)” (Figure 3A), we estimated their proportion in our samples. Figure 3B, C shows a decreasing gradient for the contribution of N, C and SE components being highest at the East, intermediate at the Mid and lowest at the West Pyrenean regions. In contrast, the SW European component shows an increasing gradient being lowest at the East, intermediate at the Mid and highest at the West Pyrenean regions. These results agree with the longitudinal gradient observed on the same region by other studies using classical markers (Calafell and Bertranpetit, 1994), as well as with the main axes of genetic differentiation within the Iberian Peninsula (Bycroft et al., 2019). Therefore, all these results support that the catalan Pyrenees fall within the genetic continuum within the Iberian Peninsula. These overall patterns have been previously interpreted in terms of ancient migrations (Calafell and Bertranpetit, 1994), as well as on the more recent migrations of Christian kingdoms and linguistic differences during the middle ages (Bycroft et al., 2019). In particular, the West to East pattern observed in the Iberian Peninsula has been explained by the admixture between Paleolithic and Neolithic populations and further geographic isolation, independently of which was the cause of that isolation (i.e., Catholic Kingdoms or linguistic differences) (Ferreiro et al., 2021).

At a local level, the presence of a higher proportion of SW European ancestry in the West region of the catalan Pyrenees supports that this region has been more accessible for migrations than other populations from the Iberian Peninsula. In addition, the higher contribution of N, C and SE European components at the East region of the catalan Pyrenees reflects a major permeability of this Pyrenean region, where geographical barriers are much weaker and with a much richer communications network, that has allowed contact with the nearby regions of northern and eastern Europe. We noted that the SE component, being significant on the Mid and East Pyrenean regions, reflects potential migratory events from Southeast Mediterranean coasts.

Given the observed differences in ancestry among the GENPIR samples, we conducted a PCA only using GENPIR samples (Figure 4A; 0.5% of variance).

We observed that the genetic map defined by the two PCs supports the presence of a non-random distribution of GENPIR individuals. Individuals from West, East and Middle Pyrenean regions tend to be grouped into three prominent clusters (Figure 4A). Since shared haplotype tracts are more sensitive for identifying subtle population substructure (Lawson et al., 2012), we wondered if the pattern observed by the PCA could be enhanced by projecting the IBD matrix between pairs of individuals estimated from Beagle 4.0 (Browning and Browning, 2013). We applied an ordinal SMACOF approach, in which the distance between pairs of individuals is interpreted in rank terms rather than by its magnitude. This decision was taken due to the fact that the overall genetic differences between all individuals are expected to be limited. Therefore, the presence of few outliers showing a higher degree of relationship (even after excluding highly related individuals) could distort the projection towards these relationships rather than describing the patterns observed among all samples. The first two dimensions from the SMACOF analysis considering 15 dimensions (stress = 0.4) produce a similar picture as observed by PCA (Figure 4B). However, it enhances the differences between the three different clusters of populations, supporting micro substructure in the Pyrenees.

Furthermore, this micro-population substructure can be interpreted in geographic terms. In particular, a procrustes analysis between the geographic coordinates of the samples and the genetic map defined by the first two PCA dimensions showed a symmetric correlation of 0.61 (p-value after 10,000 permutations <9.999e-05).

Taking into account the orography of the region, the Cadi Massif constitutes a clear geographic barrier detaching the Western to the Middle and Eastern regions while no clear geographic barrier appears to detach the Middle to the Eastern regions (Supplementary Figure S1). When applying actual (comarques) and historical (bishoprics) administrative-religious boundaries, the West cluster comprises Pallars, Ribagorça, Cerdanya and Alt Urgell comarques as well as Lleida and Urgell bishoprics; the Middle cluster comprises Solsonès and Berguedà comarques and Solsona bishopric and the East cluster includes Ripollès and Garrotxa comarques as well as Vic and Girona bishoprics (Figures 4C, D).

Given these results supporting a geographic interpretation of the genetic variation, we wondered whether models incorporating geographic information could identify genetic barriers present between the samples. The ancestry map generated by tess3r at K = 2 to K = 4 (Figure 5) identifies a strong longitudinal component that agrees with ADMIXTURE analyses and procrustes, and splits individuals following a geographic partition that matches the genetic differentiation observed in previous analyses.

Interestingly, these divisions do not fully agree on orographic events, such as the presence of mountains, that could act as genetic barriers. Particularly, the division that is observed between East and West at K = 2 is set in a region lacking clear geographic barriers. Therefore, as previously considered, other factors correlating with geography must be shaping the observed genetic differentiation across the region. We also point out that at K = 3 the algorithm mainly reproduces the West, Middle and East clustering, which we previously observed with both PCA and SMACOF analysis.

We wondered if such genetic clusters identified in previous analyses were shaped by current administrative divisions (comarques), or if they were influenced by other long term ongoing cultural factors such as the limits defined by middle age-founded bishoprics. We applied a Mclust analysis using the first two dimensions from SMACOF in order to identify clusters showing a similar probability distribution function. The optimal number of clusters identified by Mclust was six (BIC = 406.33) that accurately reflect a hardest effect of geographic barriers that is shaped by administrative divisions, both comarques and bishopric (Figure 6). Cluster I includes subjects from the comarques of Pallars and bishopric of Lleida, cluster II subjects from the comarques of Pallars and Alt Urgell as well as bishopric of Urgell, cluster III mainly includes subjects from the comarca of Alt Urgell and bishopric of Urgell, cluster IV subjects from comarca of Solsones and Berguedà as well as bishopric of Solsona, cluster V subjects from comarca of Ripollès and bishopric of Vic and cluster VI mainly includes subjects from comarca of Garrotxa and bishopric of Girona (Figure 6).

Interestingly, such bishopric division from the SMACOF bi-plot is also observed within the geographic region (depicting ∼70 km far away) that contains the intersection of four bishoprics (Figure 7).

Since bishoprics and comarques depend on geography, the identified geographic clusters by tess3r could reflect cultural and/or political barriers. We conducted a partial Mantel test between the Euclidean distance estimated from the first two PCs, and the assignment to the same administrative division or bishoprics, controlled by the geographic distance. Pearson’s correlation between PC based distances and bishoprics after controlling by geography in a partial Mantel test was 0.434 (p-value <1e-04 using 9,999 permutations). This degree of association is higher than the observed between genetics and administrative divisions of comarques (Pearson’s correlation = 0.324, p-value <1e-04 using 9,999 permutations). When considering a classification based on bishoprics and administrative divisions at the same time, Pearson’s correlation decreased to 0.294 (p-value <1e-04 using 9,999 permutations). Similar results supporting higher correlation between bishopric divisions and genetic diversity were obtained when using the first two dimensions from the SMACOF analysis conducted with the IBD distance matrix. The degree of correlation between genetics and bishoprics after controlling by geographic distance was 0.563 (p-value <1e-04 using 9,999 permutations). This correlation decreases to 0.382 (p-value <1e-04 using 9,999 permutations) when considering administrative divisions and it is similar when considering both bishopric and administrative divisions (R = 0.361, p-value <1e-04 using 9,999 permutations). Overall, these results suggest that bishoprics play a role in shaping the genetic variation of GENPIR samples by modifying their mating patterns.

We supported such a hypothesis based on genetic data by taking advantage of the sampling strategy conducted in this study, which incorporated the geographic birthplace of the four grandparents (Figure 8, Supplementary Table S1).

Figure 8 shows the grandparent to propositus network in which each propositus is linked to their four grandparents generating a complex network centered on bishopric boundaries. We can also observe that grandparent and parent couples come from the same bishopric (over 80% of couples) or to a nearest bishopric (near 10% of couples) (Figure 9A). The overall relationship between bishopric and mating was summarized by a multiple correspondence analysis (MCA) using the bishopric of each grandparent, parent and the considered sample. The relationship defined by the bishopric origin of the grandparents and samples is summarized by a triangle in the first two dimensions (explaining 39.3% of total inertia; Figure 9B) of the MCA. The first vertex contained relationships from the bishoprics of Urgell and Lleida. The second vertex was determined by the relationship between members from the bishopric of Solsona and the third one by the bishopric of Girona. The grandparents and current samples from the bishopric of Vic fall within the axis defined by Solsona and Girona, indicating a higher number of bishopric-admixed grandparent couples with the bishoprics of Solsona and Girona (Figure 9B). The center of the triangle depicted relationships where at least one of the grandparents is out of the bishoprics. Overall, the observed mating relationships agree with the geographic proximity between bishoprics (Lleida and Urgell, Vic between Solsona and Girona), as expected given the dependence between geography and these geohistorical organizations and highlight the bishopric-based mating choice of the individuals from GENPIR.

The ecclesiastical division of the modern era responds to human spheres of multisecular relationships. For this reason, it is not surprising that the bundles of family relationships roughly coincide with the bishoprics: a central Pyrenean space ‒bishopric of Urgell‒with abundant interactions (although it is capable of being subdivided into comarques); an eastern area (Girona); and a southern territory or central Catalonia (bishopric of Solsona, perhaps undifferentiated from that of Vic).

Traditionally, the main force used for explaining biases in mating choice in Spanish rural areas compared to urban areas has been geography (Fuster and Colantonio, 2004). Poor transportation facilities and limited communications - which in turn depend on orography- enhance the practice of consanguineous marriages and inbreeding in Spain. In particular, in small rural and dispersed localities, a limited number of suitable local partners conditioned marriage, and marrying a distant relative was a likely option when few other partners were available. This situation was even more common when the population started increasing due to medical advances, but mobility remained restricted (Gamella and Nez-Negrillo, 2019). Furthermore, in a sociocultural context, aunt-nephew and uncle-niece marriages could have been practiced in rural areas in order to maintain or enlarge the family inheritance (Fuster and Colantonio, 2004). Similarly, the regional practice in Northern and Northwestern Spain from the 18th to the 20th century of a complex household system could have favored first cousin mating in these rural areas (Fuster and Colantonio, 2004). Overall, these socio-cultural processes could explain the overall patterns of isolation by distance that we detected in our different analyses. However, our study also supports that bishopric can act on shaping the genetic variation of rural areas. This observation can be interpreted in the socio-cultural context of the catalan Pyrenees. In this region, the old urban centers of cities and towns with episcopal seats are the first location of markets or fairs. Cities such as La Seu d’Urgell (bishopric of Urgell), Besalú (bishopric of Girona), Berga (bishopric of Solsona) or Vic (bishopric of Vic) have documented markets or fairs since XII and XIII centuries. In rural areas, markets and fairs act as a space for sociability promoting personal contacts between relatives and acquaintances. Social relationships were established or renewed, labor contracts were established, and marriages were arranged (Gallart, 2004).