putative functional gene classes. Our study suggests that gene flow between eastern North American white oak species may not be as rampant as previously assumed and that hybridization is most strongly predicted by proximity to a species' range margin.

Interspecific gene flow is common and evolutionarily significant across the tree of life, with up to 10% of animal species and 25% of plant species estimated to engage in hybridization (Mallet 2005). Botanists in particular have long recognised the evolutionary importance of hybridization and subsequent introgressive gene flow due to backcrossing (Anderson 1953;Barton 2001;Grant 1981;Stebbins et al. 1947). Hybridization has the potential to generate novel genotypic combinations more rapidly than mutation (Anderson 1949, 62, 102;De Carvalho et al. 2010;Gompert et al. 2014), particularly within syngameons, multispecies networks of partially interfertile species (Buck and Flores-Rentería 2022;Cannon and Lerdau 2022;Grant 1981;Lotsy 1917). Syngameons have been documented in reef-building corals (Mao 2020), African lake cichlids (Olave and Meyer 2020;Schliewen and Klee 2004), Heliconius butterflies (Mallet et al. 2007) and numerous plant clades (Buck and Flores-Rentería 2022;Ellstrand et al. 1996), including many distantly related and ecologically dominant boreal, temperate and tropical trees (e.g., Buck et al. 2023;Chhatre et al. 2018;Gardner et al. 2023;Larson et al. 2021;Linan et al. 2022;Percy et al. 2014;Sun et al. 2018;Tsuda et al. 2017). Understanding landscape, climate, habitat, and demographic contexts of interspecific gene flow is key to understanding the evolutionary importance of these complex, multispecies systems.

The oak genus (Quercus L.) comprises one of the best-known syngameons (Kremer and Hipp 2020;Lazic et al. 2021). This northern hemisphere clade of ca. 425 species, mostly trees, has served as a model of plant hybridization for at least 150 years (Engelmann 1878;MacDougal 1907;Palmer 1948;Wiegand 1935), syngameon dynamics for more than 75 (Dodd and Afzal-Rafii 2004;Muller 1952;Stebbins et al. 1947;Tucker 1961) and genomic patterns of adaptive introgression for the past several years (Fu et al. 2022;Leroy et al. 2020;O'Donnell et al. 2021;Zhou, Yuan, et al. 2022). The eastern North American white oaks (Quercus sect. Quercus) have been particularly important in shaping botanists' understanding of introgressive gene flow and species concepts (Burger 1975;Coyne and Orr 2004;Van Valen 1976). Nearly all eastern North American white oak species can hybridize, but they still remain ecologically, morphologically and genomically distinct (Hardin 1975;Hipp et al. 2019).

Bur oak (Quercus macrocarpa Michx.) has the longest latitudinal range of any eastern North American white oak, ranging from north of Winnipeg, Manitoba to south of Houston, Texas, and from Nova Scotia west to Montana (Figure 1). Bur oak is a dominant of dry, fire-prone upper Midwestern oak savannas; mesic, closed-canopy forests around the Great Lakes and Northeast; and bottomland forests along the floodplains of the southcentral United States (Johnson 1990). In most of these habitats, bur oak grows in close association with at least one other eastern North American white oak species. Bur oak is sympatric with more than 12 white oak species across its range, with hybridization described between it and at least 8 naturally co-occurring species (Palmer 1948). Prior population genetic studies of bur oak have demonstrated high stand-level diversity and weak isolation-by-distance, both at smaller scales among fragmented populations (Craft and Ashley 2007;Dow and Ashley 1998) and at larger scales among rangewide samples (Hipp et al. 2019;Schnabel and Hamrick 1990;Whittemore and Schaal 1991). However, previous studies have relied upon a limited number of loci to assess nuclear patterns of population structure in the species, and none have included sufficient sampling of other species to estimate introgression across the range of the species.

To evaluate patterns and determinants of interspecific gene flow in the eastern North American white oak syngameon, we sampled bur oak and the four most prevalent white oak species sympatric with it: eastern white oak (Q. alba L.), a dry-mesic forest species that ranges from near the Canada-U.S. border to northern Florida and Texas; swamp white oak (Q. bicolor Willd.), a wet forest species mostly restricted to north of 37°; post oak (Q. stellata Wangenh.), common in drier locales south of approximately 40°; and chinkapin oak (Q. muehlenbergii Engelm.), an upland limestone specialist mostly contained within the range of Q. alba, but extending a bit further into eastern Kansas, Oklahoma and Texas (Figures 1 and 2). We sequenced 465 target genes in 517 individuals across 56 sites in eastern North America. We addressed the following questions in our study: (1) How prevalent is interspecific gene flow in the syngameon, and what are the spatial and ecological drivers of introgression? (2) How do levels of introgression and/or genetic differentiation vary over candidate gene classes potentially implicated in environmental adaptation? (3) How is genetic variation geographically partitioned in Q. macrocarpa, and does genetic variation in the species follow the ecological gradient observed from predominantly upland sites in the north to increasingly bottomland sites in the south?

Sampling was designed to cover the majority of the Q. macrocarpa range and include co-occurring white oaks at as many sites as possible. A focused field sampling of wild populations was conducted in the summer of 2017; these wild populations were supplemented with botanical garden material from Starhill Forest Arboretum (Petersburg, IL, USA) grown from acorns of known wild provenance (thus representing wild populations our group did not visit in the field) (Figure 1). At each wild population, leaf samples were collected from up to three individuals of Q. macrocarpa and of each co-occurring white oak species found at the site (Figure 2). Conspecific samples were taken at least 30 m apart from each other in wild populations, when possible, to avoid genotyping close relatives. Botanical garden material was selected based on wild provenance, so spatial distance within gardens was disregarded. All trees collected from wild populations were georeferenced in the field using a hand-held GPS unit. Botanical garden samples from acorns of known wild provenance could not be georeferenced precisely and were consequently excluded from spatial analyses (Table S1).

Morphological intermediates with potential hybrid ancestry were neither intentionally targeted nor avoided in the course of sampling; clearly intermediate individuals were determined to the closest morphospecies. Species determinations made in the field were checked and, if necessary, updated by examining leaf vestiture, gross leaf morphology, and twig and bud characters according to descriptions provided in the Flora of North America (Nixon and Muller 1997). After harvesting, leaf specimens were placed on ice in the field, then frozen at -80°C within 1 day to 1 week of collecting for later genomic analysis. Voucher samples were collected for all individuals and deposited at the Morton Arboretum Herbarium (MOR). Duplicate vouchers were deposited at the Bell Museum at the University of Minnesota (MIN) (herbarium acronyms from Thiers, Updated continuously [accessed 2025-02-02] n.d.).

We selected target loci and designed sequence capture exon baits based on gene annotations of coding domain sequences in the Q. robur (Plomion et al. 2016(Plomion et al. , 2018) ) and Q. lobata (Sork, Fitz-Gibbon, et al. 2016) draft assemblies, as well as leaf transcriptomes of Q. robur, Q. alba and Q. rubra (Lesur et al. 2015) (WO454_v2, Hardwood Genomics Project, https:// doi. org/ 10. 25504/ FAIRs haring. srgkaf; RO454_v2, Hardwood Genomics Project, https:// doi. org/ 10. 25504/ FAIRs haring. srgkaf). We used MarkerMiner v1.0 (Chamala et al. 2015) to identify 403 longexon, single-copy nuclear genes shared across Quercus references that were expected to provide high-quality phylogenetic resolution. 30 candidate genes for bud-break phenology and waterlogging response, 30 for drought tolerance and 2 for freezing tolerance were included as targets (Lesur et al. 2015;Meireles et al. 2017;Oney-Birol et al. 2018;Ueno et al. 2013). Target loci are provided as a fasta file in Appendix S1 (Table S2).

DNA was extracted from frozen leaf tissue using a DNeasy Plant Mini Kit (Qiagen) with a modified DNeasy protocol as described in Hipp and Weber (2008). DNA was quantified using a Qubit Fluorometer and visualised on an agarose gel to confirm DNA quality. DNA was fragmented to approximately 550 bp using a Covaris M220 (Covaris, Woburn, MA, USA) or NEBNExt dsDNA fragmentase (New England Biolabs, Ipswich, MA, USA). Library preparation used the Adapterama iTru system and KAPA Hyper Prep Kit (Roche Diagnostics Corporation, Wilmington, MA, USA), with KAPA Pure Beads, KAPA HyperPure Beads, or Speedbeads (Gardner et al. 2021;Glenn et al. 2016;Hale et al. 2020;Rohland and Reich 2012). Pools of 12-16 libraries were enriched using our 465-gene custom Quercus MYbaits kit (Morales-Saldaña et al. 2024) and either the Mybaits v3 or v4 reagent kit (Arbor Biosciences, Ann Arbor, MI, USA). Hybridization was done at 65°C for 16 h with ½ volume of baits used per enrichment. QIAquick PCR purification kit (Qiagen) was used to purify PCR-amplified enriched libraries. Pooled libraries were quantified using a Qubit Fluorometer and run on a Bioanalyzer (Agilent, Santa Clara, CA, USA) to ascertain fragment size distributions before these were combined for sequencing. Samples were subsequently sequenced using the Illumina MiSeq (2 × 300 bp) or HiSeq 4000 (2 × 150 bp) platforms. Two-sample t-tests with sequencing platform as the predictor of latitude were used on all samples as well as Q. macrocarpa samples to test whether there was a latitudinal bias in the sequencing platform used.

Raw FASTQ files were cleaned with Trimmomatic 0.39 (Bolger et al. 2014) to remove ILLUMINA adapters, leading and trailing bases with Phred score < 15, portions of reads with an average Phred score < 20 across a 4 bp sliding window, and cleaned reads shorter than 40 bp. Subsequently, reads were mapped to the Q. lobata genome v3.2 (Sork et al. 2022) using BWA-MEM 0.7.17 (Li and Durbin 2009). Duplicate reads were marked using Samtools 1.17 (Danecek et al. 2021). Variants were called with FreeBayes 1.3.6 (Garrison and Marth 2012) using default settings. After decomposing multinucleotide polymorphisms and retaining only biallelic SNPs, we implemented an iterative filtering pipeline adapted from O'Leary et al. (2018, pipeline FS5). Site genotyping rates and missing data per individual filters were repeatedly applied within the same pipeline according to progressively more stringent thresholds, removing low-quality variants while maximising the number of retained samples and SNPs.

Filtering was performed through a custom pipeline of vcflib 1.0.3 (Garrison et al. 2022), VCFtools 0.1.16 (Danecek et al. 2011) and BCFtools (Danecek et al. 2021); see details in methods in Appendix S1. After filtering, SNPs were thinned to a minimum distance of 500 bp using PLINK 1.9 (Purcell et al. 2007), as prior studies in oaks have demonstrated rapid decay of linkage disequilibrium, both within and across individual genes (Kremer et al. 2012;Nocchi et al. 2022;Sork, Squire, et al. 2016). Scripts to conduct variant-calling and SNP-filtering are archived on Zenodo (v 0.9-3, https:// doi. org/ 10. 5281/ zenodo. 15216009). Sequence reads are archived on NCBI (SRA Bioproject #1223965; http:// www. ncbi. nlm. nih. gov/ biopr oject/ 1223965).

Ancestry estimation was performed using the maximumlikelihood software ADMIXTURE (Alexander et al. 2009), which jointly optimises the allele frequencies of a predefined number of ancestral clusters (K) and the proportional membership ("Q-value") of each individual in each cluster based on SNP genotype data. Fifty replicate runs were carried out for each value of K from K = 2 to K = 8, and the run with the highest log-likelihood per K was selected for visualisation and downstream analysis, as performed by Marcus et al. (2021). Evaluation of the optimal value of K was performed using the default ADMIXTURE cross-validation procedure enabled by the "--cv" flag. We visualised ADMIXTURE ancestry estimation and cross-validation results with ggplot2 (Wickham 2016) in R version 4.3.0 (R Core Team 2023). We also ran ADMIXTURE on downsampled datasets with 20 randomly sampled individuals per species and only putatively "pure" Q. macrocarpa individuals (Q > 0.95) to assess the robustness of our dataset to detect population structure within this species.

Ancestral clusters in the K = 5 and K = 6 runs corresponded almost exactly to the 5 morphospecies included in our study, with the exception of Q. macrocarpa, which was separated into two admixed subpopulations at K = 6. For subsequent analyses using ADMIXTURE Q-values, we defined the conspecific Q-value for non-Q. macrocarpa individuals as the maximum ancestry proportion in the K = 6 cluster, which almost invariably matched a sample's morphospecies determination. For samples identified as Q. macrocarpa, the conspecific Q-value was calculated as the sum of the Q-values for the two Q. macrocarpa clusters estimated in the K = 6 run. Heterospecific Q-values were calculated as 1-conspecific Q-value.

We tested for spatial autocorrelation of the northern Q. macrocarpa ancestral cluster (hereafter termed MAC N , in contrast to the southern genetic cluster subsequently referred to as MAC S ) using global Moran's I (Moran 1950). "Pure" (conspecific Qvalue > 0.95) Q. macrocarpa individuals were grouped by sampling site. An inverse distance matrix for spatial weighting was calculated using each site's centroid coordinate. Pairwise geographic distances between sites were calculated using the Haversine formula, assuming a spherical globe (Sinnott 1984). Moran's I was calculated with the MAC N Q-value for each Q. macrocarpa individual using adegenet 2.1.10 and spdep 1.3_3 in R (Bivand 2022;Jombart 2008;Jombart and Ahmed 2011). We evaluated the significance of Moran's I via permutation test with 99,999 resamples to obtain a one-tailed p-value, corresponding to an alternative hypothesis of positive spatial autocorrelation.

To further clarify whether genetic clustering within Q. macrocarpa reflected either continuous genetic variation or genetic discontinuity, we examined spatial genetic variation using the Python package Fast Estimation of Effective Migration Surfaces (FEEMS) (Marcus et al. 2021). FEEMS uses deme-aggregated allele frequencies at spatial nodes to fit edge-specific effective migration estimates, inferring landscape-wide patterns of heterogeneous isolation-by-distance via a penalised likelihood approach. We subsetted the thinned SNP dataset to only precisely georeferenced Q. macrocarpa samples with conspecific Q-value > 0.95 and removed any resulting invariant sites ("-maf" PLINK flag) from the PLINK binary files. Demes were constructed by aggregating individuals to their nearest node in a 100 km-resolution spatial grid. We performed leave-oneout cross-validation to determine the optimal value of the model tuning hyperparameter λ, evaluating 20 values of λ from 10 -6 (low smoothing) to 10 2 (high smoothing). All analyses were conducted using modified FEEMS tutorial scripts published by the package authors (https:// github. com/ Novem brelab/ feems ).

To more precisely characterise hybrids, we used NewHybrids, a Bayesian ancestry inference algorithm that uses patterns of observed genotype frequencies and inferred population allele frequencies to estimate the posterior probability of assignment to a predefined set of hybrid classes (Anderson and Thompson 2002), specifying expected genotype frequencies for multiple backcross generations (Table S3). Identifying hybrids based on genotype frequency is more tractable with loci that approach fixation between populations (Wringe et al. 2017). For this reason, and because NewHybrids can only consider pairwise species comparisons, we selected the 200 most highly differentiated loci per species pair using the getTopLoc function in the R package HybridDetective 4.3.1 (Wringe et al. 2017) based on putatively pure individuals identified in ADMIXTURE runs (conspecific Q-value > 0.95). The choice of 200 loci for analysis was based on computational limits of NewHybrids and supported by a previous finding that including more than 10 markers per linkage group (more than 120 loci in the case of oaks) provides minimally increased power to detect F1s and first-generation backcrosses (Chakraborty and Rannala 2023). All individuals with cumulative two-species Q-values > 0.95 across pairwise parental ancestral clusters were considered for each run of NewHybrids; individuals with substantial hybrid parentage from more than two species were therefore excluded, as NewHybrids does not model 3-way hybrids. For each species pair, four independent MCMC chains were run with an uninformative Jeffreys prior for a burn-in period of 300,000 iterations prior to 600,000 sampling sweeps. Traceplots and effective sample sizes, computed across all chains per pairwise analysis for the population-wide hybrid class membership proportion parameters via the R package coda 0.19-4 (Plummer et al. 2006), were inspected to evaluate MCMC convergence. As a nonparametric check on the hybrid class assignment inferred by NewHybrids, we plotted ancestry proportion against interclass heterozygosity for all individuals in the Q. alba-Q. macrocarpa and Q. bicolor-Q. macrocarpa comparisons, using triangulaR v.0.0.1 (Wiens and Colella 2024).

BIOCLIM rasters were downloaded at 30 arcsecond resolution from the WorldClim 2 database (Fick and Hijmans 2017). Moisture index I m was calculated as 100 × MAP -PET PET (Thornthwaite 1948), where MAP is mean annual precipitation (BIO12) and PET is potential evapotranspiration extracted from the Global Aridity Index and Potential Evapotranspiration Database v3 (Zomer et al. 2022). Little's (1971) species range maps, digitised and compiled by Prasad and Iverson (2003), were downloaded as shapefiles (https:// github. com/ andre w-hipp/ white -oak-synga meon/ tree/ master/ data/ littl e-maps; doi: 10.5281/zenodo.13150341). For all spatial analyses, coordinates were transformed to the North American Albers equal area conic projection to ensure accurate raster grid sizing. Projections and all subsequent spatial calculations, analyses and raster extractions were performed using the R packages sf v1.0-14 and terra v1.7-39 (Hijmans 2019(Hijmans -2025;;Pebesma 2018).

Given range-margin uncertainty and long-distance pollen dispersal (Ashley 2021), we extended all range limits by 20 km for analyses involving range calculations. Distance from range edge was calculated as the distance between a sample's coordinate and the nearest buffered range edge, with lake and ocean shorelines masked to avoid artificial range truncations. Number of sympatric species was calculated as the number of eastern North American white oak ranges (of the 5 species analysed in our study), aside from an individual's species determination, overlapping with the coordinates of a sample. Soil data were obtained from the USDA-NRCS Soil Survey Geographic Database (SSURGO; Soil Survey Staff, Natural Resources Conservation Service, United States Department of Agriculture, n.d.) accessed through the Web Soil Survey API using a custom R script.

Species occurrence records were obtained from the Global Biodiversity Information Facility (GBIF) database and previously curated USDA Forest Inventory and Analysis (FIA) data (Cavender- Bares et al. 2018) (https:// gitlab. com/ meire les/ compi le_ fia_ data/ ). GBIF records were filtered to remove erroneous coordinates and data points corresponding to county or state centroids, fossils and planted specimens, urban occurrences and coordinates over bodies of water using R packages rgbif 3.7.7 and CoordinateCleaner 2.0.20 (Chamberlain et al. 2024;Chamberlain and Boettiger 2017;Zizka et al. 2019). FIA data as analysed in Cavender- Bares et al. (2018) were used to summarise species co-occurrence, treating the presence of two species within any of the subplots of a single FIA plot as co-occurrence.

Preliminarily cleaned occurrence records from both sources were merged, downloaded and manually curated to remove spatial outliers in QGIS 3.28 (QGIS Development Team 2023) and thinned to a minimum distance of 20 km via spThin 0.2.0 (Aiello-Lammens et al. 2015). Study extent was defined by the 250 km-buffered convex hull of occurrence points for each respective species, with 20,000 unique, randomly sampled raster grid cells serving as a background dataset. We used all 19 BIOCLIM variables as environmental predictors in our SDMs, as Maxent's regularisation algorithm is robust to predictor multicollinearity (Merow et al. 2013). SDMs were fit using MAXENT 3.4.3 as implemented in the R package ENMeval 2.0.4 (Kass et al. 2021;Phillips et al. 2017). See methods in Appendix S1 for model-fitting details.

We employed Bayesian generalised linear mixed models (GLMMs) in brms v2.20.1 (Bürkner 2017;Carpenter et al. 2017) to examine associations between environmental and ecological factors and introgression and Q. macrocarpa population structure. Beta regression with a logit link and default flat priors on all parameters was used for all models, and sampling site was incorporated as a random effect to control for spatial autocorrelation and reduce the possible risk of pseudoreplication (Douma and Weedon 2019). See methods in Appendix S1 for details of model-fitting and MCMC diagnostics.

To test the association of genetic differentiation in Q. macrocarpa with the gradient from more southern bottomlands to more northern uplands, we included as fixed effects climate and soil predictors expected to vary along moisture and temperature gradients, with expectations as follows: temperature seasonality (BIO4) and moisture index (I m ), based on our previous work demonstrating that both vary strongly among North American oak species (Hipp et al. 2018); maximum temperature of the warmest month (BIO5), as a proxy for drought severity; mean annual precipitation (BIO12), as a predictor of overall plant water requirements; and categorical soil predictors for ponding frequency, flooding frequency and drainage class, as predictors of wetland status (which we expect would differentiate Q. alba and Q. stellata as upland species most strongly from Q. bicolor, Q. lyrata and bottomland ecotypes of Q. macrocarpa). For modelling population structure within Q. macrocarpa, we analysed only individuals with conspecific Q-value > 0.95, including site as a random effect. Whole-dataset admixture was modelled by regressing heterospecific Q-value against latitude, number of sympatric species, Maxent-inferred habitat suitability, and distance from range edge, with species as a fixed categorical predictor and a random effect for site. Because we were interested in understanding the drivers of gene flow between species, rather than merely first-generation hybridization, F1s identified by NewHybrids were excluded from the admixture GLMM. See methods in Appendix S1 for modelling details.

Single-nucleotide polymorphisms were annotated to indicate which candidate genes they fell within, and genes with fewer than three SNPs were excluded from jackknife analysis. To test whether genetic differentiation varied significantly between candidate gene classes, we obtained per-class F ST distributions by jackknife resampling (10,000 replicates) of 15 genes per class (bud break phenology, drought tolerance, and other) and three SNPs per gene. Pairwise F ST estimates (calculated per Weir and Goudet 2017) for each combination of two species were computed using the R package hierfstat v. 0.5-11 (Goudet 2005) for each jackknife replicate. We compared the means of the F ST distributions for the BBP and DT classes to that of the 'Other' genes to calculate two-tailed pvalues. See methods in Appendix S1 for details on establishing gene orthology, categorising genes, and analysis.

Our post-filtering SNP dataset consisted of 3412 biallelic SNPs (Figure 3, black dots) with an average SNP genotyping rate of 87% across all individuals and a mean read depth of 22.0 ± 12.1 SD across all SNPs, averaged over all samples. Variants were recovered from 445 of the 465 target genes (Figure 3, red bars); nine of the targets initially identified in Q. robur failed to BLAST to the Q. lobata genome and 11 lacked any associated SNPs. Target loci and SNPs were distributed across all (as mapped to Q. lobata genome v3.2). Of filtered and thinned SNPs, 22% fell outside of the target gene regions. Given that Hyb-Seq can recover high-copy, off-target nuclear regions (Weitemier et al. 2014), and because we applied stringent filtering parameters to remove erroneous and poorly resolved paralogous calls, we chose to retain off-target SNPs for subsequent analyses. ADMIXTURE results below are those from analyses conducted with the full dataset; results using only SNPs confined to target loci produced highly similar ancestry estimates (not shown

). Of the 517 individuals sequenced, 364 (224 of 319 Q. macrocarpa [70.2%], 64 of 78 Q. alba [82.1%], 30 of 60 Q. bicolor [50.0%], 25 of 31 Q. stellata [80.6%] and 21 of 29 Q. muehlenbergii [72.4%]) passed the iterative SNP filtering pipeline (Table S1). 318 individuals with associated GPS coordinates or manual georeferences were retained for downstream spatial analyses, including GLMMs (Table 1). Sequencing platform was unbiased with respect to sample latitude across species (t = -0.99667, df = 184.63, p = 0.3202) as well as within just Q. macrocarpa (t = -0.34215, df = 230.86, p = 0.7325).

Non-metric multidimensional scaling (NMDS) ordination of the SNP distance matrix revealed well-defined, discontinuous clusters corresponding to the 5 species sampled, with a much smaller number of admixed individuals found at varying distances between parental clusters (Figure 4a,b). ADMIXTURE's cross-validation procedure for determining the optimal number of ancestral clusters (K) returned a minimum prediction error at K = 6 (0.17007) and slightly higher errors for K = 5 (0.17232) and K = 7 (0.17214) (Figure 5). Cross-validation errors increased considerably for values of K smaller than 5 and larger than 7, supporting both specieslevel genomic delineations (K = 5; Figure 5a) and a population genetic divide within Q. macrocarpa (which appears only at K ≥ 6; Figure 5b, Figure S1). As these genomic clusters within Q. macrocarpa are strongly geographically structured between northern and southern populations (see results below, Spatial and environmental predictors of admixture and population structure), we refer to these clusters as MAC N and MAC S throughout the paper.

Of the 364 individuals included in the total-dataset ADMIXTURE runs, 19 (5.2%) were estimated as being F1s (17) or three-way hybrids (2) based on having their maximum Qvalue < 0.60; all but one has at least 0.30 assignment to Q. macrocarpa (the exception being an inferred F1 between Q. alba and Q. muehlenbergii). The two three-way hybrids inferred include Q. alba as the dominant genotype (Q = 0.52, 0.57) and Q. macrocarpa and Q. bicolor as the only other genomic components of Q > 0.05 (Q = 0.14-0.32). Some of these might also be F2, F3 or other generation hybrids, but NewHybrid analyses demonstrate that we have minimal power to distinguish among these classes. Of the remaining samples, 6 (1.6%) have a maximum Q-value < 0.9, 18 (4.9%) have maximum Q-value < 0.95 and 56 (15.4%) have maximum Q-value < 0.99 (Table S4).

The distribution of introgressed (not hybrid; i.e., mixed ancestry Site State Latitude Longitude Sample size Sakatah Lake MN 44.6243 -93.3961 3 Q. macrocarpa Twin Cities MN 44.9557 -93.1620 4 Q. macrocarpa Winona MN 43.8992 -91.6427 1 Q. alba; 1 Q. macrocarpa; 1 Q. muehlenbergii Eminence MO 37.1573 -91.3650 1 Q. alba; 3 Q. macrocarpa; 1 Q. muehlenbergii Eureka MO 38.5122 -90.5657 1 Q. alba; 1 Q. macrocarpa; 1 Q. muehlenbergii; 2 Q. stellata Ha Ha Tonka MO 37.9732 -92.7623 2 Q. alba; 1 Q. muehlenbergii; 1 Q. stellata Shaw MO 38.4639 -90.8100 3 Q. macrocarpa Sullivan MO 38.2226 -91.0855 3 Q. alba; 1 Q. bicolor; 2 Q. macrocarpa; 3 Q. stellata Bienville MS 32.3582 -89.5472 8 Q. stellata Daughmer OH 40.7299 -83.0936 3 Q. macrocarpa Goll Woods OH 41.5522 -84.3574 3 Q. alba; 2 Q. macrocarpa; 3 Q. muehlenbergii Pearl King OH 40.0440 -83.4794 2 Q. macrocarpa Hinton OK 35.4453 -98.3543 2 Q. macrocarpa; 2 Q. muehlenbergii Norman OK 35.1771 -97.4496 3 Q. macrocarpa Pawhuska OK 36.8455 -96.4232 1 Q. macrocarpa; 3 Q. stellata Tulsa OK 36.2177 -95.8987 1 Q. alba; 3 Q. macrocarpa Custer SD 43.7887 -103.3634 2 Q. macrocarpa Spearfish SD 44.4913 -103.8090 1 Q. macrocarpa Buescher TX 30.0764 -97.1817 2 Q. stellata Addison VT 43.9580 -73.1881 3 Q. macrocarpa Bennington VT 42.9102 -73.2193 3 Q. macrocarpa Burlington VT 44.4682 -73.2109 2 Q. alba; 2 Q. macrocarpa Ferrisburgh VT 44.2398 -73.2444 1 Q. alba; 3 Q. macrocarpa Grand Isle VT 44.6925 -73.3354 2 Q. macrocarpa Rutland VT 43.6341 -73.1529 1 Q. alba; 2 Q. macrocarpa Dodgeville WI 43.0187 -90.1162 3 Q. alba; 7 Q. macrocarpa Gotham WI 43.2451 -90.3034 5 Q. bicolor; 3 Q. macrocarpa Madison WI 43.0391 -89.4394 3 Q. alba; 3 Q. macrocarpa Ottawa WI 43.0168 -88.4384 3 Q. alba; 3 Q. macrocarpa Pleasant Valley WI 43.1068 -89.8063 1 Q. alba; 1 Q. macrocarpa Round Lake WI 46.0247 -91.1393 3 Q. macrocarpa Note: Sites represent the closest named locality to the population; state/province abbreviations are standard. Latitude and longitude are computed as the centroid of the coordinates for all individually georeferenced samples at a given site and rounded to 4 decimal places.

TABLE 1 | (Continued) 1365294x, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mec.17822 by University Of Chicago Library, Wiley Online Library on [29/09/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

other than inferred F1) individuals within each species at the 0.60 < Q < 0.95 threshold was 14 out of 217 samples (6.5%) in Q. macrocarpa, 1 out of 58 (1.7%) in Q. alba, 2 out of 29 (6.9%) in Q. bicolor, 1 out of 22 (4.5%) in Q. stellata, and 0 out of 19 (0%) in Q. muehlenbergii or 18 out of 345 (5.2%) for the whole dataset (Table S4). Within Q. macrocarpa, samples were much more likely to have assignment to both MAC N and MAC S than to exhibit admixture with different species: of the 203 samples with MAC N + MAC S > 0.95, 119 (58.6%) had Q < 0.95 membership in one of the two Q. macrocarpa genetic clusters. However, genetic differentiation between MAC N and MAC S is lower than interspecific variation: pairwise F ST from ADMIXTURE between inferred Q. macrocarpa genetic clusters was 0.103, compared to 0.292-0.468 for interspecific comparisons (Figure 5).

362 out of 364 individuals had a summed Q > 0.95 across at most two species clusters and were included in pairwise NewHybrids analyses. Using a 0.9 posterior probability threshold of hybrid class assignment, 314 individuals (86.7%) were confidently classified as "pure" (unadmixed) members of their determined species, 16 (4.4%) were classified as F1s, 7 (1.9%) were classified as backcrosses (BC1, BC3, or BC of undetermined generation), and 2 (0.55%) were classified as F2s (Table S5; Data Summary S1). Only three of the seven backcrosses (0.8%) could be confidently assigned to a specific backcross generation. Additionally, discrimination between putatively pure individuals and thirdgeneration backcrosses proved difficult for NewHybrids; 23 (6.3%) individuals could not be unambiguously placed into either class (Table S5). The percentage of admixed individuals in each pairwise comparison ranged from 0% (Q. stellata-Q. muehlenbergii, Q. bicolor-Q. stellata, Q. bicolor-Q. muehlenbergii, Q. alba-Q. stellata, Q. alba-Q. bicolor) to 11.6% (29 out of 248)

for the Q. macrocarpa-Q. bicolor pair (Table S6; Data Summary S1). NewHybrids analyses converged across all MCMC chains, with effective sample sizes of > 100,000 for all parameters across species pairs (Table S7). Triangle plots for Q. macrocarpa × bicolor and Q. macrocarpa × alba (Figure 6) demonstrate that F1s FIGURE 6 | Triangle plots. Interancestry heterozygosity plotted against hybrid index (i.e., ancestry proportion) for the (a) Q. macrocarpa-Q. alba and (b) Q. macrocarpa-Q. bicolor species pairs. Samples are coloured according to their NewHybrids consensus hybrid class assignments. Grey plus signs indicate the expected position of specified hybrid classes on the triangle plots. Expected heterozygosity as a function of ancestry fraction under Hardy-Weinberg conditions is given by the dashed parabola.

and F2s are strongly concordant with evidence from heterozygosity for both crosses. However, the inferred backcrosses beyond BC1 are not readily distinguishable from pure Q. macrocarpa in the cross with Q. bicolor, whereas the BC of undetermined generation toward Q. alba appears to be a plausible BC1.

Of the FIA plots subsetted to just those that had at least one oak, 39.5% had Q. alba, 14.3% had Q. stellata, 12.7% had Q. macrocarpa, 5.6% had Q. muehlenbergii and 2.2% had Q. bicolor. The most common associate of Q. macrocarpa was Q. alba, at 19.4%, followed by Q. muehlenbergii at 6.0%, then Q. stellata and Q. bicolor, both at 3.5%. The most common associate of Q. alba was Q. stellata (23.2%), followed by Q. muehlenbergii (7.5%), Q. macrocarpa (6.3%) and Q. bicolor (2.6%).

Proportional membership in ADMIXTURE-inferred Q. macrocarpa genetic clusters displayed significant spatial autocorrelation (Moran's I = 0.088391, p = 0.009), broadly corresponding to northern (MAC N ) and southern (MAC S ) genetic clusters (Figure 7). Cross-validation tuning of the FEEMS hyperparameter λ suggested that a modest degree of smoothing (λ = 2.069) yielded the best model fit (Figure S2). Congruent with ADMIXTURE results, analysis of spatial genetic connectivity via FEEMS revealed a 300-500 km latitudinal band of reduced effective migration stretching from southern Ohio in the east to northern Kansas and southern Iowa in the west (Figure 8). Effective migration rates to the north and south were approximately 1-2 orders of magnitude larger than rates within the band of reduced effective migration.

Model selection via AICc favoured species distribution models (SDMs) with small regularisation multipliers and many feature classes (Table S8). SDMs projected across eastern North America corresponded closely to published range maps and known occurrence data for each species in our dataset (Little 1971) (Figure S3). Range edges displayed steep habitat suitability FIGURE 8 | Bur oak effective migration landscape estimated by FEEMS. Black points represent sampled populations of Quercus macrocarpa. Line segments forming the lattice are coloured (and thickened) according to log-transformed migration rates (w). Migration estimates are reported as the log 10 -transformed edge weights (w ik ) standardised against a scaling parameter w 0 estimated under a pure isolation-bydistance model, i.e., a constant value of w ik across all edges.

gradients for most species, especially toward northern and western margins (Figure S3).

In the all-species admixture GLMM, distance from range edge exhibited a strong negative association with heterospecific Qvalue ([-0.48, -0.13], 95% CI) (Figure 9a). No other environmental predictors were significantly correlated with admixture rates, nor were any of the species fixed effects, though the small number of hybrids found may limit our power to detect any such effects. Conversely, the Q. macrocarpa population structure GLMM revealed little predictive power for any individual environmental covariate (Figure 9b). As assessed by traceplots, posterior predictive checks and visual examination of parameter posterior distributions, our GLMMs exhibited convergence across all MCMC chains, with R-hat = 1.00 and cumulative Bulk/Tail ESS > 82,000 across all model parameters for all models (Tables S9 and S10).

Rates of genetic differentiation across candidate gene classes for bud-break phenology and drought tolerance did not diverge significantly from the null distributions constructed from background (non-candidate) genes (Figure S4, Table S11).

Our study demonstrates that the five most prevalent species of the eastern North American white oak syngameon are genomic cluster species (Mallet 1995(Mallet , 2020)), largely predictable by their ecology and morphology. All, however, are connected by gene flow with the widespread bur oak or directly with one another. Although the species form clear, discontinuous genomic clusters, ca. 13% of the individuals we sampled are hybrids or introgressed. The eastern North American white oaks are a syngameon.

Our study also shows that genetic variation in Q. macrocarpa is due in part to introgression and in part to a substantial phylogeographic split separating the northern and southern populations. Our work thus supports the analysis of the eastern North American white oak syngameon made by Hardin (1975, 336) on the basis of morphology alone, but here supported by genomic data: while "a significant, although relatively minor, component of the variation is due to hybridization and localized introgression," most of the variation in the species is due to intraspecific variation that apparently has nothing to do with hybridization.

Understanding the relative importance of these sources of genetic diversity will be key to predicting how bur oak, one of eastern North America's foundational oak species, will adapt to evolving environmental conditions.

Both the NewHybrids and ADMIXTURE analyses demonstrated that the eastern North American white oaks in our study are genomically distinct species that nonetheless hybridize and, in some cases, introgress. Of 362 samples, the NewHybrids analyses identified 16 (4.4%) as F1 hybrids and 9 (2.5%) as backcrosses or F2, for a total of 6.9% confirmed hybrids. In addition, a further 23 (6%) are of undetermined ancestry (meaning they could be either unadmixed or backcross), making 315 samples (87%) confirmed "pure" (unadmixed).

The estimated frequency of hybrids in our study is much higher than estimates from historical, morphology-based studies, which generally estimated the frequency of hybrids in eastern North American oaks as less than 1% (Jensen 2002;Palmer 1948). This discrepancy may be due to both sampling bias-in at least one study, oak hybrids were found to be undersampled in herbaria relative to their prevalence in natural populations (Wu et al. 2023)-and the increased power of genomic data to detect admixture. Nonetheless, our finding recalls the claim made 50 years ago based on morphology alone that "the term 'syngameon' [in the eastern North American white oaks] does not imply unlimited and widespread gene exchange between the species. To the contrary, the gene 'flow' is but a 'trickle' in most cases" (Hardin 1975, 360).

Our admixture findings are in line with other DNA-based studies of North American oaks (e.g., Cavender- Bares et al. 2015;Craft et al. 2002;Kim et al. 2018;Ortego et al. 2017;Sullivan et al. 2016). Higher hybridization rates have been reported for oaks of Europe (e.g., Degen et al. 2023;Jensen et al. 2009;Reutimann et al. 2020;Valbuena-Carabana et al. 2005), East Asia (e.g., Liu et al. 2018; Pei et al. 2022; Qi et al. 2024; Shi et al. 2023) and Mexico (e.g., Albarrán-Lara et al. 2019; Castillo-Mendoza et al. 2019; González-Rodríguez et al. 2004; Morales-Saldaña et al. 2022; Peñaloza-Ramírez et al. 2010), but this may be influenced in part by human interplanting of species that would not grow together naturally in Europe (Stace et al. 2015) and substantially younger clade ages in Mexico (Hipp et al. 2018).

21 apparent backcrosses from Q. bicolor to Q. macrocarpa are scored as ambiguous in the NewHybrids analysis, but most make sense geographically: all Q. macrocarpa with ≥ 3.0% admixture from Q. bicolor are located within the range of Q. bicolor except for two individuals (QUE001882 from Oklahoma; QUE001635, cultivated from acorn collected in New Mexico); and four of the five Q. macrocarpa individuals with ≥ 2.0% admixture from Q. muehlenbergii are in the area of sympatry with that species (the exception being QUE002581, from South Dakota). However, ancestry estimation methods can fail to recover signals of ancient introgression, particularly when the magnitude of gene flow is low and occurred long enough in the past for drift, selection and recombination to fix variants and reshuffle species-specific combinations of allele frequencies (Kong and Kubatko 2021). The rates and spatial patterns of admixture observed in our dataset should be interpreted as a snapshot of one point in the history of the eastern North American white oak syngameon and possibly an under-representation of the importance of admixture.

Proximity to species' range edge most strongly predicted admixture in our dataset, concordant with prior rangewide findings on introgression in Q. alba based on morphological data (Hardin 1975). Similar geographic patterns have been observed in poplar, pinyon pine and other temperate forest trees (Buck et al. 2023;Chhatre et al. 2018;Mimura et al. 2024) as well as clades of butterflies, toads, grasshoppers and other taxa (Barton andHewitt 1985, 1989;Kawakami et al. 2009;Shepherd et al. 2022). These patterns align well with predictions and observations of environmental intermediacy favouring genetic admixture (Anderson 1948;Muller 1952). However, unlike these systems, which are largely characterised by parapatric hybridization, the eastern North American white oak syngameon contains many species occurring in widespread sympatry.

Range-marginal hybridization may be due to pollen swamping from larger heterospecific populations (Lepais et al. 2009), as populations reach their geographic limits, or reproductive failure due to increased stress in ecologically marginal populations, as demonstrated in Quercus gambelii and Q. grisea, where hybridization was driven in part by pollen abortion due to heat and drought stress (Williams et al. 2001). A disproportionate number of admixed individuals in our study were identified in the Hudson Valley-Lake Champlain corridor in New York and Vermont, where the eastern North American white oak species present consist of disjunct or range-marginal populations confined to lower elevation river valleys between climatically unsuitable mountain ranges. Under these circumstances, compressed distributions along topographically diverse river valleys might also facilitate higher rates of gene flow between species (Vallejo-Marín and Hiscock 2016). However, increased gene flow in the north also raises the question of whether adaptive introgression might contribute to the large geographic range and ecological variation in Q. macrocarpa. Adaptive introgression can facilitate species migration and persistence in novel environments (Burge et al. 2019;Cronk and Suarez-Gonzalez 2018;Dodd and Afzal-Rafii 2004;Nagamitsu et al. 2020;O'Donnell et al. 2021;Pfennig et al. 2016) and has been argued to play a key role in the northward extension of the widespread European species Q. petraea (Leroy et al. 2020). Climate suitability was not a significant partial predictor of admixture in our study (Figure 9a, "habitat suitability"), but edaphic factors at a finer scale than we were able to measure using the soil layers we had may be more important, as has been suggested, for example, in eastern North American Q. ellipsoidalis and Q. rubra (Khodwekar and Gailing 2017).

Even along range margins, hybrids were largely confined to a small subset of sites. Since early-generation hybridization appears rare, this heterogeneity could simply be explained by a low probability of a hybrid occurring at any given site. However, heterogeneity in local factors, such as habitat structure, disturbance (Anderson 1948;Bray 1960) or community composition, may exert strong effects on the generation and persistence of F1s and their backcrossed progeny. More detailed data on micro-scale habitat variation could help clarify whether habitat intermediacy promotes the survival of F1s and facilitates introgressive gene flow back into one or more parent populations.

The two most closely related species in our study, bur oak (Q. macrocarpa) and swamp white oak (Q. bicolor) (Hipp et al. 2020;Larson et al. 2024), have long been thought to exchange genes relatively freely (Bray 1960;Burger 1975;Stebbins 1950, 63-64;Van Valen 1976). While our linear mixed models (GLMMs) suggest that the partial effect of species identity on admixture is non-significant (Figure 9), this may be partially due to low within-population sampling. By contrast, ADMIXTURE and NewHybrids analyses both show exceptionally high interspecific gene flow between Q. bicolor and Q. macrocarpa: well-supported F2s, 4 of the 7 demonstrated backcrosses, and 21 of the 23 samples of undetermined ancestry (unadmixed or backcross) were derived from these two species. This is particularly remarkable in light of the fact that in U.S. Forest Inventory and Analysis data, 19.4% of Q. macrocarpa plots also include Q. alba, while only 3.5% of Q. macrocarpa plots include Q. bicolor. This suggests that contemporary gene flow is much more frequent between Q. macrocarpa and Q. bicolor than expected by co-occurrence alone.

Given the close evolutionary relationship between the two species and high gene tree discordance observed among white oaks in general (Crowl et al. 2020;Larson et al. 2024;McVay, Hauser, et al. 2017;McVay, Hipp, et al. 2017), it is conceivable that our inferred admixture estimates could represent either gene flow or unsorted ancestral polymorphism. Incomplete lineage sorting between recently diverged species can masquerade as admixture when phylogeny is not accounted for (Lexer et al. 2006;Muir andSchlötterer 2005, 2006;Thomson et al. 2015). However, in our dataset, admixed individuals of Q. macrocarpa and Q. bicolor were distributed unevenly across the landscape, clustering at specific sites and extending no further than the zone of sympatry between parent species. This pattern supports interspecific gene flow over incomplete lineage sorting: under the latter, inferred admixture is expected to be widespread across the ranges of parental populations (Muir and Schlötterer 2005).

The causes of high gene flow between Q. bicolor and Q. macrocarpa as well as the strong asymmetry in introgression, almost exclusively from Q. bicolor to Q. macrocarpa, are unknown. The two recently diverged species may have lower reproductive barriers due to insufficient time for intrinsic incompatibilities to evolve. Moreover, Q. macrocarpa and Q. bicolor co-occur closely in mixed stands throughout the Midwest and portions of the Northeast, though relatively infrequently (3.5% of Q. macrocarpa plots). Close co-occurrence could provide ample opportunity for heterospecific pollen flow and potentially relax postzygotic selection against hybrids with intermediate ecological niches. However, most of the eastern North American white oak species in our study commonly grow in codominant mixed stands, particularly Q. macrocarpa and Q. alba (Cavender-Bares et al. 2018; Cho and Boerner 1991; Hardin 1975); and as noted above, more than 5 times as many Q. macrocarpa plots in the U.S. FIA data include Q. alba as include Q. bicolor. Disentangling the mechanisms maintaining genomic coherence of Q. macrocarpa and Q. bicolor would require controlled crossing experiments to tease apart the relative contributions of pollen competition, gametic incompatibility, hybrid viability and ecological selection on hybrids, among other factors.

Genetic differentiation did not vary significantly across candidate functional gene classes for any of the species pairs in our study (Figure S4). We may have lacked adequate sampling of Quercus individuals and target loci-our post-filtering dataset includes only 445 of the ca. 31,000-42,000 gene models inferred in a variety of annotated oak genomes to date (e.g., Ai et al. 2022;Fu et al. 2022;Kapoor et al. 2023;Larson et al. 2024;Sork et al. 2022;Zhou, Liu, et al. 2022)-to detect differential rates of divergence or gene flow across classes: adaptive differentiation or introgression may occur at a small number of loci, yielding insignificant results when averaging over an entire functional gene class consisting mostly of loci not under positive selection. Additionally, the drought-tolerance candidate genes in our marker set were selected based on orthology and literature search for California oaks (Oney-Birol et al. 2018); the bud-break phenology candidate genes were selected based on a transcriptomic study in European white oaks (Lesur et al. 2015; Ueno et al. 2013). It is unknown whether these loci are all under selection, particularly in eastern North America. Alternatively, because we combined all samples of a given species for pairwise analyses, spatially heterogeneous patterns of gene flow or differentiation could obscure regionally significant divergence, gene flow or local adaptation. For example, adaptive introgression of drought tolerance alleles from Q. stellata into southwestern populations of Q. macrocarpa could be counterbalanced by introgression of flooding-tolerance alleles from Q. bicolor at northern sites. Further studies using larger marker sets and denser population sampling could provide the power to geographically partition jackknife analyses or conduct genotype-environment association tests to examine finer-scale, locus-or environment-specific patterns of (potentially) adaptive genetic differentiation across the eastern North American white oak syngameon. Finally, adaptive divergence and/or gene flow may be focused around non-coding regions, as has been demonstrated in east Asian Quercus (Fu et al. 2022), potentially producing null results when testing for adaptive introgression using target-capture methods for protein-coding genes.

Genomic coherence of Q. macrocarpa is strongly supported by NMDS and ADMIXTURE analyses. Divergence between the northern and southern populations of the species is supported by ADMIXTURE and FEEMS, the latter showing a reduction in gene flow near the middle of the species' range, thus suggesting that the phylogeographic structure within Q. macrocarpa is not simply an outcome of isolation-by-distance. Moreover, our inclusion of a phylogenetically diverse set of widespread eastern North American white oak species (Hipp et al. 2018) suggests that population structure within Q. macrocarpa, beyond evidence of admixture and hybridization with the included species, represents intraspecific variation rather than gene flow from an unsampled eastern North American white oak species. Are the bur oak north and south genomic clusters distinct enough to be good species? This seems unlikely, as genetic differentiation between inferred "pure" Q. macrocarpa north and south clusters (F ST = 0.103) is approximately one-third as great as that between the most closely related species in our dataset, Q. bicolor and Q. macrocarpa (F ST = 0.292-0.304). In combination, these analyses support our inference that, despite population structure, Q. macrocarpa is a single genomic cluster species, distinct from sympatric congeners.

The phylogeographic signal within Q. macrocarpa is not significantly predicted by climate or soil, leaving open the question of whether selection along an environmental gradient from northern upland savannas to southern bottomlands might explain genetic divergence within the species. Genetic clustering within Q. macrocarpa also does not correspond to known natural barriers to reproductive connectivity in eastern North America, nor to any widely documented phylogeographic patterns in the region (Soltis et al. 2006). It also does not align with a previous phylogenetic study of chloroplast haplotypes across the eastern North American white oak syngameon (Pham et al. 2017), which in fact showed weak geographic clustering, in contrast with the strong geographic structure in European oak chloroplast haplotypes (Petit et al. 2002). These findings are consistent with the previous inferences that North American oaks occupied more and smaller northern refugia and underwent less dramatic late-Pleistocene population bottlenecks than European oaks (Lumibao et al. 2017;McLachlan et al. 2005).

The genetic structure within bur oak does not correspond to any known morphological discontinuities. While latitudinal variation exists for bur oak acorn size and leaf lobing, these characters vary continuously (Desmond et al. 2021;Koenig et al. 2009). Additionally, none of the varieties previously diagnosed within Q. macrocarpa, which consist mostly of regional morphs throughout the Great Plains (Nixon and Muller 1997), align with the pattern of population structure we obtained. Rangewide common garden approaches hold the promise of evaluating what, if any, implications rangewide population structure holds for functional phenotypic diversity in Q. macrocarpa.

Our study is the first genomic investigation of hybridization, introgression and population genetic structure in bur oak and sympatric species of the eastern North American white oak syngameon, a model system for examining hybridization and species coherence in temperate trees. Our study confirms that the traditional species form largely discontinuous genomic clusters, despite the presence of interspecific hybrids. Several interesting trends are suggested by our results. Levels of hybridization vary across the bur oak range, with hybrids establishing at a higher level in range-marginal populations. Other trends suggested in this dataset require further testing, since the small percentage of hybrids prevented strong statistical confirmation. Patterns of hybridization vary among species pairs. For instance, Q. macrocarpa showed little backcrossing with Q. stellata and Q. alba, although F1 hybrids with these species were found. On the other hand, backcrosses and F2s are much commoner for Q. bicolor × Q. macrocarpa. Backcrossing in this combination is asymmetrical, however, with Q. bicolor genes becoming incorporated into Q. macrocarpa but little evidence for the reverse. The demonstration of strong phylogeographic structure within the species suggests the possibility of either secondary contact from multiple refugia or strong divergent selection. The high amount of introgression into Q. macrocarpa raises the question of whether the large ecological, geographic, and morphological range of the species may be due in part to adaptive gene flow.

Combined with the species distribution models presented, these data are a first step toward understanding how the interplay between genetic diversity within and among populations, genetic diversity among species and interspecific gene flow has shaped one of eastern North America's most widespread oak species.

and hybridization to managers of participating natural areas where we collected tissue (see Acknowledgments) and will inform conservation decisions in at least one state. Additional benefits accrue from the sharing of our data and results on public databases as described above.