Genetic diversity and population structure in *Nothofagus pumilio*, a foundation species of Patagonian forests: defining priority conservation areas and management

Patagonian forests are highly endemic, species-poor temperate forests isolated from other continental forests for millions of years and provide critical ecological, economic, and social services. Nothofagus pumilio (lenga) is a dominant foundation tree spanning ~20° latitude (36–55°S) across heterogeneous environmental gradients, exhibiting marked phenotypic variation and adaptability. These forests face increasing threats from climatic change (warming, drought, decreasing precipitation) and anthropogenic impacts (logging, increased fire frequency). Effective management requires understanding the distribution and structure of genetic diversity across the species’ range to delineate operational genetic management units (OGMUs) such as Genetic Zones (GZs) that allow safe movement of reproductive material and targeted conservation. Prior efforts in Argentina defined OGMUs for other native species, and national law emphasizes conserving native forests, making N. pumilio a priority. Given evidence for historical lineage breaks, hybridization with N. antarctica, and Quaternary glaciation shaping genetic patterns, this study aims to delineate genetically homogeneous regions (GZs) across the Argentine range of N. pumilio, identify priority conservation areas based on allelic richness and exclusive alleles, and provide guidance for conservation and management.

Previous research has shown that Quaternary climatic variability and glaciation shaped current distributions and genetic structures of Patagonian forests. In Nothofagus pumilio, phylogeographic studies based on cpDNA reported a north–south lineage break around ~42°S and multiple glacial refugia, with pollen dispersal exceeding seed dispersal. Natural hybridization with N. antarctica occurs in sympatry, with chloroplast capture documented. Operational genetic management units have been delineated for other Argentine native trees (Austrocedrus chilensis, N. nervosa, N. obliqua), and the concept of Genetic Zones (genetically homogeneous regions enabling safe seed transfer) contrasts with Provenance Regions (based on adaptive trait variation). Microsatellite diversity is often correlated with genomic diversity at the population level and has been used to define management units, although neutral markers may not reflect adaptive variation, highlighting the need for complementary adaptive data when defining Provenance Regions and Conservation Units. Argentina’s legal framework (Law 26331) emphasizes conserving native forests, with most N. pumilio forests in high or moderate conservation value zones, underscoring the need to identify priority conservation areas based on genetic criteria such as allelic richness and the presence of locally common or private alleles.

Sampling and genotyping: The study analyzed 35 natural populations of Nothofagus pumilio across Argentina (36–55°S; 71°56′–66°37′W), totaling 965 adult trees. At least 20 trees per population were sampled (except Cholila: 18; El Chaltén: 19), with minimum 50 m spacing to avoid close relatives, and only N. pumilio phenotypes were included. DNA was extracted from buds. Seven species-specific nuclear microsatellite loci (Npum1, Npum3, Npum9, Npum10, Npum13, Npum17a, Npum18) were amplified using standard PCR conditions (M13 labeling; GoTaq or Platinum Taq), genotyped on an ABI 3130xl, and scored with GeneMarker. Genetic diversity and inbreeding: For each population, the number of alleles (NA), effective number of alleles (Ne), expected heterozygosity accounting for null alleles (HE), and Nei’s unbiased heterozygosity were estimated (GenAlEx, Info-Gen). Differences among populations were assessed by Wilcoxon tests. Null alleles were detected with MICRO-CHECKER and their frequencies estimated by FreeNA. Inbreeding coefficients (FIS) accounting for null alleles were estimated using INEST via two approaches: Bayesian Individual Inbreeding Model and Maximum Likelihood, testing models including null alleles, genotyping errors, and inbreeding; the best model (lowest DIC) was used. Standard errors were jackknifed across loci, and significance assessed by Z-tests. Demographic history: Recent bottlenecks were tested with BOTTLENECK under IAM, SMM, and TPM models, using Wilcoxon sign-rank tests; a bottleneck was inferred if heterozygosity excess was significant under at least two models. Population differentiation and IBD: Genetic differentiation was quantified by FST (with and without ENA correction) and standardized G’ST, with 95% CIs via bootstrap. Isolation-by-distance was assessed with (a) Mantel tests correlating FST with log geographic distance (Info-Gen) and (b) spatial genetic autocorrelation (GenAlEx) over ten 200 km distance classes, with 1000 permutations and 10,000 bootstraps for confidence intervals. Clustering and zone delineation: Given the broad range and complex history, clustering considered multiple criteria. Prior cpDNA evidence identified a lineage break near ~42°S, and new haplotypes were obtained for Cholila. Initial broad BAPS analysis on all populations showed weak grouping, prompting subdivision into two datasets: ‘North’ (12 populations with northern cpDNA haplotypes; Epulauquen to Cholila) and ‘South’ (23 populations with southern haplotypes; Huemules to Bahía Lapataia). For each subset, spatially explicit clustering was performed with BAPS 6.0 to infer K, followed by STRUCTURE admixture analyses (500,000 MCMC steps; 100,000 burn-in; five replicates per K). UPGMA dendrograms (Nei’s distance) explored relationships among clusters. Grid-based spatial visualization used R (Raster) and ArcGIS, averaging admixture coefficients (>0.7) within a circular 20-minute neighborhood. Further stratification of the ‘South’ subset considered geomorphological history: 42°30′–44°S (distinct watershed orientations influencing glaciation) analyzed separately; 44–52°S (excluding Tierra del Fuego) analyzed with attention to admixture in Alto Río Senguer; Tierra del Fuego populations analyzed as a distinct unit due to island isolation and local recolonization dynamics. Genetic Zones (GZs) and validation: GZs were delineated by integrating BAPS clusters, STRUCTURE admixture patterns, phylogeographic context (glacial refugia, recolonization), topography, mountain height, and water bodies. AMOVA (Arlequin 3.5) evaluated hierarchical partitioning among and within proposed GZs to confirm the number of zones (maximizing among-group variance and minimizing within-group variance). For each GZ, genetic parameters were calculated: private alleles (P), locally common alleles (LCA; frequency >5% occurring in <25% of GZs), Nei’s unbiased heterozygosity (HE), and allelic richness (AR; rarefied to a common sample size using FSTAT). Maps of GZ distributions were generated in ArcGIS.

- Genetic diversity and inbreeding: Populations differed significantly in unbiased heterozygosity (p<0.05). Tromen (Tr) had the highest HE (0.757), significantly greater than PG (0.581), Pi (0.573), BL (0.537), and HP (0.398). BL and HP showed the lowest heterozygosity. Signs of inbreeding were detected in 12 populations, but only Cerro Nahuelpan (Np) had a significant FIS (0.1461; Z>1.96). A recent bottleneck was inferred in Sierra Colorada (SC) with significant heterozygosity excess under IAM and TPM (p<0.05). - Differentiation and IBD: Overall FST with ENA correction was 0.077 (95% CI: 0.050–0.113); without ENA, 0.078 (0.053–0.115). Standardized G’ST (from GST=0.075) was 0.201. No isolation-by-distance was detected by Mantel or SGS tests using SSRs; cpDNA studies show contrasting patterns due to differential pollen vs seed dispersal. - Major phylogeographic division: A north–south lineage split occurs between Cholila (42°30′S; northern haplotypes) and Huemules (42°50′S; southern haplotypes), used as the first criterion for structuring analyses. - Genetic Zones (GZs) delineated (18 total): • North (12 populations; 36°49′–42°30′S): Five GZs based on BAPS clusters and STRUCTURE admixture: - GZ1: Epulauquen (E) - GZ2: Caviahue (Cav), Batea Mahuida (BM), Tromen (Tr) - GZ3: Quilanlahue (Q), Paso Puyehue (Pu), Cerro Otto (Ot), Valle del Chalhuaco (V), Cerro Piltriquitrón (Pi) - GZ4: Ea. Herodina Parada (HP) - GZ5: Mina de Indios (M), Cholila (Ch) GZ2 and GZ3 exhibited high allelic richness (AR), locally common alleles (LCA), and private alleles. • 42°30′–44°S (seven/eight populations analyzed separately): Four GZs: - GZ6: Trevelin (Te), Cerro Nahuelpan (Np) - GZ7: Huemules (Hm) - GZ8: La Hoya (H), José de San Martín (JSM) - GZ9: Sierra Colorada (SC), Lago Guacho (G), Lago Engaño (Eg) GZ6 showed high genetic diversity; GZ9 had high LCA and private alleles; SC (in GZ9) showed a recent bottleneck. • 44°–52°S (excluding Tierra del Fuego; 10 populations): Six GZs: - GZ10: Lago Azul (LA) - GZ11: Lago Fontana (F) - GZ12: Arroyo Perdido (AP) - GZ13: Río Unión (U) - GZ14: Monte Zeballos (MZ), Ea. Tucu Tucu (Tu), El Chaltén (N), Punta Gruesa (Sa) - GZ15: Cancha Carrera (CC), Mina I (MI) High LCA was found in GZ14 (2.14); AR values across these GZs were similar and generally lower than northern GZs. • Tierra del Fuego Island (five populations): Three GZs: - GZ16: Tierra del Fuego Norte (TdFN), Tierra del Fuego Centro (TdFC), Tierra del Fuego Este (TdFE) - GZ17: Paso Garibaldi (PG) - GZ18: Bahía Lapataia (BL) GZ16 exhibited high LCA and no private alleles; it is categorized as of medium conservation value and under active management. - AMOVA validation: With 18 GZs, among-group variance was 5.71% (FCT=0.057, p=0.000), among populations within groups 2.20% (FSC=0.023, p=0.000), within populations 23.03% (FIS=0.250, p=0.000), and within individuals 69.05% (FIT=0.309, p=0.000), supporting the proposed grouping. - Priority conservation areas: Based on high allelic richness (AR) and/or exclusive alleles (private alleles and LCA), six GZs were prioritized: GZ2, GZ3, GZ6, GZ9, GZ14, and GZ16.

The study addressed the need to define operational genetic management units for Nothofagus pumilio by integrating neutral genetic variation (nSSR), phylogeographic history, and geomorphological context to delineate 18 Genetic Zones across Argentina. The absence of isolation-by-distance with SSRs and the clear north–south lineage break guided spatially explicit clustering and justified stratified analyses. AMOVA confirmed that the proposed GZs capture meaningful genetic structure, enabling seed transfer guidelines that minimize disruption of local genetic composition. Identifying six priority conservation areas—those with high allelic richness and/or exclusive (private or locally common) alleles—directly supports conservation planning by highlighting genetic reservoirs and putative refugial signatures. These areas are essential sources for in situ conservation, ex situ collections, restoration seed sourcing, and low-intensity breeding. Zones with balanced admixture (GZ7, GZ13) require tailored management to preserve unique genetic mosaics, recommending local or balanced mixed seed sourcing if restoration is needed. While microsatellites are neutral and may not capture adaptive variation, their population-level correlation with genomic diversity provides a robust first step for defining OGMUs, to be complemented by adaptive trait data when establishing Provenance Regions and Conservation Units.

This work provides the first range-wide, genetics-based management framework for Nothofagus pumilio in Argentina. Eighteen Genetic Zones were delineated and mapped, validated by hierarchical variance partitioning, and six priority conservation areas (GZ2, GZ3, GZ6, GZ9, GZ14, GZ16) were identified based on allelic richness and exclusive alleles. These outputs supply actionable guidance for seed transfer, restoration, and in situ conservation, supporting compliance with national native forest protection law and ongoing restoration projects. The GZs constitute a foundational step toward defining Provenance Regions and Conservation Units that incorporate adaptive variation. Future work should integrate adaptive trait and genomic data, extend sampling to Chilean populations to enable coordinated trans-Andean management, and monitor demographically vulnerable zones (e.g., bottlenecked SC within GZ9).

- Geographic scope: Sampling was limited to Argentina; Chilean populations were not included, potentially omitting trans-Andean genetic structure relevant for coordinated management. - Marker set: Only seven nuclear SSR loci were used; while informative for neutral diversity and structure, they may not reflect adaptive genetic variation and provide limited genomic resolution compared to genome-wide markers. - Neutral vs adaptive variation: Management units based on neutral markers may not align with adaptive Provenance Regions; additional assessments of adaptive traits and genomic scans are needed to inform assisted migration and restoration under climate change. - Temporal inference: Bottleneck detection and phylogeographic inferences rely on model assumptions and limited marker sets; recent demographic changes may be under- or over-detected. - Mapping resolution: Admixture and zone boundaries are inferred from discrete sampling locations and spatial interpolation, which may not capture fine-scale structure in unsampled areas.