Phylogenomic analyses highlight innovation and introgression in the continental radiations of Fagaceae across the Northern Hemisphere

Abstract

Northern Hemisphere forests changed drastically in the early Eocene with the diversification of the oak family (Fagaceae). Cooling climates over the next 20 million years fostered the spread of temperate biomes that became increasingly dominated by oaks and their chestnut relatives. Here we use phylogenomic analyses of nuclear and plastid genomes to investigate the timing and pattern of major macroevolutionary events and ancient genome-wide signatures of hybridization across Fagaceae. Innovation related to seed dispersal is implicated in triggering waves of continental radiations beginning with the rapid diversification of major lineages and resulting in unparalleled transformation of forest dynamics within 15 million years following the K-Pg extinction. We detect introgression at multiple time scales, including ancient events predating the origination of genus-level diversity. As oak lineages moved into newly available temperate habitats in the early Miocene, secondary contact between previously isolated species occurred. This resulted in adaptive introgression, which may have further amplified the diversification of white oaks across Eurasia.

Fig. 1. Phylogenetic relationships and divergence time estimation of Fagaceae inferred from analyses of 2124 nuclear genes.

A The global climate curve during the last 82 million years (modified from ref. ). Major climate events were indicated. B Rate-through-time plot showing the net diversification rate (species/million years) of Fagaceae. Red line is the median and the blue shadow represents the 95% confidence interval. C Chronogram derived from ASTRAL-III tree based on concatenated nuclear data. Nodes showing consistent relationships between ASTRAL-III, SVDquartets, maximum likelihood, and MrBayes are marked with red (phylogenetic support ≥ 95% in all four analyses) and blue (support < 95% in any one of the four analyses). Nodes showing conflicting relationships among analyses are marked with black dots. Light blue bars on nodes represent 95% confidence intervals of divergence time estimates and dashed vertical red line represents the age of the Cretaceous-Paleogene boundary (66 million years ago). Geological timescale is shown at bottom. Fossil calibration nodes are indicated with C1–C8 (stem calibration node; Supplementary Table S1). S1–S4 indicate four nodes where shifts in diversification rate were identified. Taxonomic labels of genera, subgenera and sections follow refs. ^{, ,} . Illustrations: lax catkins indicate the placement of the change from insect-pollination to wind-pollination that diagnoses the genus Quercus; hypogeous seed and seedling marks the origin of the HS clade. Images: representative cupule types are shown on the right. A consistent color scheme was used for taxonomic labels and image borders. Ma, million years ago; Pli, Pliocene; Ple, Pleistocene.

Fig. 2. Conflicts between nuclear (left) and plastome (right) species trees.

Pie charts on nodes indicate the geographic distribution of the clade (black = Old World, white = New World). The hypogeous seed (HS) clade consists of six genera divided into two major plastome clades: New World (light gray) and Old World (dark gray). Lineage colors are consistent with the color scheme in Fig. 1. Protobal, section Protobalanus; Cyclobalanop, section Cyclobalanopsis; Notholith, Notholithocarpus; Chryso, Chrysolepis; Trigono, Trigonobalanus.

Fig. 3. Gene flow between Fagaceae species revealed by D-statistic test.

A Number of species-pairs with significant D-value between sections of Quercus and other genera. Numbers on diagonal line indicate gene flow within each section or genus. Cells are colored based on the ratio of species-pairs with gene flow, with warmer colors indicating a higher proportion of species-pairs showing gene flow. For example, significant gene flow was detected for 12 species-pairs between sections Quercus (white oak) and Ponticae, representing 41% of tested species-pairs between these two sections. B Distribution of D-values for white oaks vs. Q. pontica (left), and two species of section Virentes, Q. virginiana (middle) and Q. fusformis (right). Each line summarizes a set of D-statistic tests performed on trios in the format ((H1,H2),H3) with different H1 species and fixed H2 and H3 species (one of the three species above; sample size n = 4, 4, and 7 species for H1 as NA, EU, and AS, respectively). Both H1 and H2 were white oaks, but represent different lineages. For example, if H2 was a North American white oak, then H1 was sampled from European or Asian white oaks. In each panel, points represent mean D-values and error bars represent minimum and maximum D-values across multiple tests. EU = European white oak; AS = Asian white oak; NA = North American white oak. A negative D-value indicates gene flow between H1 and H3 while a positive D-value indicates gene flow between the H2 taxon and H3. Q. pontica shows a clear pattern of gene flow with EU and AS white oaks but not with NA white oaks while the opposite pattern is recovered for Q. virginiana and Q. fusiformis. The significance of the D-value was tested by a two-sided standard block-jackknife procedure implemented in Dsuite v0.3 with default parameters, and adjusted by Bonferroni correction for multiple comparisons.

Fig. 4. Shared IBD blocks between Quercus species.

A Heatmap indicating the total length of identity-by-descent (IBD) blocks for each pair of comparisons. B, C box plots show shared total length of IBDs between sections Ponticae and Quercus, and between sections Virentes and Quercus. In these box plots, the horizontal lines indicate the median value, the bottom and top of each box represents the first and third quartiles, and the whiskers extend to 1.5 times the interquartile range (the sample size n = 4, 5, and 14 individuals for NA, EU and AS, respectively). NA = North American white oak; EU = European white oak; AS = Asian white oak. D Kernel distribution of the length of shared IBD blocks between sections. Vertical black line (at 11,724 bp) indicates the shortest IBD block that is significantly longer than the expectation for selectively neutral introgressed fragments maintained in a population under a constant recombination rate of 10^–8 per site per year, assuming an average divergence time of 3 million years (P = 0.0476; two-sided probability estimated under a Gamma distribution function, and adjusted by Bonferroni correction; see details in Methods).