<figure><img src="image_1.png" /></figure>

HAL
open science

# Asymmetric character displacement in mixed oak stands

Fang K Du, Min Qi, Yuan-ye Zhang, Rémy J. Petit

► To cite this version:

Fang K Du, Min Qi, Yuan-ye Zhang, Rémy J. Petit. Asymmetric character displacement in mixed oak stands. New Phytologist, 2022, pp.1-13. 10.1111/nph.18311. hal-03745435

HAL Id: hal-03745435

https://hal.inrae.fr/hal-03745435

Submitted on 4 Aug 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Asymmetric character displacement in mixed oak stands

Fang K. Du¹ ID, Min Qi¹ ID, Yuan-Ye Zhang² ID and Rémy J. Petit³ ID

¹School of Ecology and Nature Conservation, Beijing Forestry University, Beijing 100083, China; ²Key Laboratory of the Ministry of Education for Coastal and Wetland Ecosystems, College of the Environment and Ecology, Xiamen University, Xiamen, Fujian 361102, China; ³Université de Bordeaux, INRAE, BIOGECO, F-33610 Cestas, France

Author for correspondence:
Fang K. Du
Email: dufang325@bjfu.edu.cn

Received: 22 October 2021
Accepted: 28 May 2022

New Phytologist (2022)
doi: 10.1111/nph.18311

**Key words:** allopatry, asymmetric gene flow, ecological character displacement (ECD), leaf shape, *Quercus aliena*, *Quercus dentata*, sympatry.

## Summary

* Ecological character displacement (ECD) refers to a pattern of increased divergence at sites where species ranges overlap caused by competition for resources. Although ECD is believed to be common, there are few in-depth studies that clearly establish its existence, especially in plants.

* Thus, we have compared leaf traits in allopatric and sympatric populations of two East Asian deciduous oaks: *Quercus dentata* and *Quercus aliena*. In contrast to previous studies, we define sympatry and allopatry at a local scale, thereby comparing populations that can or cannot directly interact.

* Using genetic markers, we found greater genetic divergence between the two oak species growing in mixed stands and inferred that long-term gene flow has predominantly occurred asymmetrically from the cold-tolerant species (*Q. dentata*) to the warm-demanding later colonizing species (*Q. aliena*). Analysis of leaf traits revealed greater divergence in mixed than in pure oak stands. This was mostly due to the later colonizing species being characterized by more resource-conservative traits in the presence of the other species. Controlling for relevant environmental differences did not alter these conclusions.

* These results suggest that asymmetric trait divergence can take place where species coexist, possibly due to the imbalance in demographic history of species resulting in asymmetric inter-specific selection pressures.

## Introduction

Ecological character displacement (ECD) is a process that enhances phenotypic differences between species through resource competition in sympatry (Brown & Wilson, 1956; Grant, 1966). This important process, one of the most thoroughly studied phenomena in evolutionary biology, plays key roles in speciation, divergence, diversification, and ultimately community assembly (Beans, 2014; Germain et al., 2018). Whereas ECD is often framed around competition between reproductively isolated species, it is not limited to such situations and can be investigated in species still exchanging genes (Goldberg & Lande, 2006; Reifová et al., 2011). However, rigorously demonstrating ECD's presence in empirical field studies is difficult, for several reasons (Schluter & McPhail, 1992; Stuart & Losos, 2013). First, assessing whether species can actually compete is not always trivial. Second, environmental conditions can vary across sites, potentially confounding effects of sympatry on studied traits. Third, for species still undergoing gene exchanges, increased divergence in sympatry can result either from greater competition for resources (ECD) or from reinforcement of reproductive barriers (Noor, 1999; Hopkins, 2013). Therefore, a rigorous test of ECD in sympatry, controlling for all these potentially confounding effects, is essential. As demonstrated by Stuart & Losos (2013), such tests are very rare, leading them to

advocate for 'more complete, rather than simply more, cases studies'. This applies especially to plants, which have received less attention in these respects than animals (Beans, 2014) even though they should be good models for such studies due to their sessile character.

Traditional field methods used to study divergence between species generally involve comparison of traits or mating barriers in sympatric and in allopatric regions, broadly defined at the range-wide scale. A limit of this approach is that, without precise and rigorous characterization of the degree of sympatry, other unstudied factors might complicate the interpretations. A review by Butlin et al. (2008) concluded that in speciation studies it is more productive to study the current balance between local adaptation and gene flow than to focus on broad and overly simplistic classification of spatial context. In line with this recommendation, we argue here that a simple solution for field studies is to concentrate on local scale patterns, by comparing populations in 'microsympatry' and 'microallopatry' within broadly sympatric parts of their ranges (Fig. 1), so the potential for direct interaction at each site is clear and the comparison is more relevant.

Even with such improved sampling scheme, environmental gradients can obscure patterns of character displacement (Goldberg & Lande, 2006). Therefore, in such studies, accounting for phenotypic variation along environmental gradients is essential (Adams & Collyer, 2007).
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Schematic model showing the classical and local (microgeographical) approaches to test character displacement hypotheses. (a) Classical approach: divergence of traits is tested by contrasting allopatric vs sympatric populations at a broad regional scale. (b) Local approach: divergence of traits is tested by contrasting allopatric vs sympatric populations defined at a local scale, i.e. in situations where the two species can potentially directly interact. Note that for the selection of microallopatric populations, the allopatric parts of the range, broadly defined, are avoided, to make results more comparable with microsympatric populations.</figcaption></figure>

Another challenge when studying ECD in plants is to deter-
mine the character traits to investigate, because character dis-
placement can occur in both subtle and clearly visible traits
(Beans, 2014). Geometric morphometric methods (GMMs) for
landmark-based analysis of shape variation have important
advantages (Adams & Rohlf, 2000). Leaf shape, which is often
used for taxonomic purposes, is now recognized as a trait
with great functional significance (Nicotra et al., 2011). With a
GMMs, all aspects of leaf shape variation among a set of land-
marks can be captured, so specific changes can be detected with-
out having to be specified and explicitly measured a priori
(Dryden & Mardia, 1998; Klingenberg, 2011). Other leaf func-
tional traits related to resource acquisition and hence to competi-
tion, such as leaf area and specific leaf area (SLA), are also of
great potential for ECD studies in plants (Cavender-Bares &
Pahlich, 2009; Cavender-Bares et al., 2015; Arenas-Navarro
et al., 2020).

Competition between closely related plants can lead not only
to ECD but also to reinforcement of prezygotic barriers (Jiggins
& Mallet, 2000). Focusing on leaf traits that are not directly
involved in reproduction and jointly monitoring populations
with molecular genetic markers should help clarify the underlying
mechanisms. For instance, it could help clarify if asymmetric
selection due to demographic imbalance could lead not only to
asymmetric gene flow but also to asymmetric ECD (Fig. 2).

The oak genus (Quercus spp.) is one of the most diverse and
ecologically important tree genera in the Northern Hemisphere,
with high species diversity in southeast Asia and both Central
and North America (Denk et al., 2018). Oaks are good models
for studies of community structure and niche construction, with
surveys showing that co-occurring species are often more dis-
tantly related than expected by chance (Cavender-Bares
et al., 2004). Oak leaves are hugely variable in sizes and shapes.
The possible adaptive response of this variation to environmental
conditions has been extensively investigated, perhaps more
so than in any other plant genus (e.g. Zwieniecki et al.,
2004; Niinemets, 2015; Hipp et al., 2018; Cavender-Bares,
2019; Royer et al., 2019; Li et al., 2021; Sancho-Knapik
et al., 2021; Skelton et al., 2021). The oak genus is also a classic
example for hybridization and introgression studies (Dar-
win, 1872; Rieseberg et al., 2006). Oak species are typically well
delimitated despite occasional hybridization thanks to high
fecundity combined with strong disruptive selection (e.g. Gailing
& Zhang, 2018; Cannon & Petit, 2019). In particular, there is
documentary evidence of often strong but plastic and context-
dependent reproductive barriers between species of oaks in
Europe (Abadie et al., 2012; Lepais et al., 2013); North America
(Cavender-Bares & Pahlich, 2009); and East Asia (e.g. Liu
et al., 2018; Lyu et al., 2018). Interestingly, recent studies have
suggested that oak species colonizing stands already occupied by

<figure><img src="image_5.png" /><figcaption>**Fig. 2** Trait differences between two species in sympatry and allopatry: (a) trait divergence does not rely on ecological variables and is symmetrical: classical ecological character displacement (ECD); (b) trait divergence varies along an environmental gradient; (c) trait varies along an environmental gradient in late-successional (thick dashed lines) species but not in resident species (thick solid lines) (asymmetric ECD). (a, b) come from Goldberg & Lande (2006), with permission. Note that in all three examples an extreme case is displayed where species traits are identical in allopatry under similar environmental conditions (thin gray line).</figcaption></figure>
related oak species initially face conditions conducive to massive hybridization while they are still at low density (Lepais *et al.*, 2009; Lagache *et al.*, 2014). This can result in strong and asymmetric selection against inter-specific mating events, as described in *Drosophila* (Yukilevich, 2012).

In China, oaks are the most abundant trees in naturally regen-
erated forests, accounting for c. 13% of total numbers of trees in
them, with 35 species described, 15 of which are endemic to the
country (Huang *et al.*, 1999). The two species selected for this
study (*Quercus aliena* Blume and *Quercus dentata* Thunberg)
belong to a small monophyletic group of oak species nested
within Section *Quercus* (white oaks) (Hubert *et al.*, 2014; Hipp
*et al.*, 2020). They have widely overlapping distributions and co-
occur side by side in some forests. *Quercus aliena*, a large tree
species growing in mixed forests that can tolerate shade and com-
petition, is considered a characteristic species of temperate and
warm-temperate forests, while *Q. dentata*, an often-stunted tree
found in open habitats, grows further north and is considered a
characteristic species of temperate forests (Huang *et al.*, 1999;
Fujiwara & Harada, 2015). Both species can be discriminated by
leaf shape. Leaves of *Q. aliena* have 10–13 mm petioles with ellip-
tic leaves characterized by many shallow lobes. Leaves of *Q. den-
tata* are obovate with much shorter petioles (2–5 mm) and
fewer but deeper lobes than *Q. aliena* (Huang *et al.*, 1999; Liu
*et al.*, 2018).

For ECD studies, leaf morphology could provide suitable
functional traits (Zwieniecki et al., 2004; Royer et al., 2008;
Nicotra et al., 2011; Kusi & Karsai, 2019; Torres-Ruiz et al.,
2019; Sancho-Knapik et al., 2021). To help interpret differences
in leaf morphology, more integrated traits, such as leaf area and
SLA, can also be used, as their relationships with the environment
have been extensively investigated (Wright et al., 2004). In partic-
ular, associations between SLA and vulnerability to embolism are
widely reported in oaks (Cavender-Bares et al., 2004; Nardini
et al., 2012; Lobo et al., 2018). More drought tolerant species
tend to have lower SLA (e.g. Alexandre et al., 2020). Similarly, a
reduction in leaf size helps oaks withstand water deficits (e.g.
Peguero-Pina et al., 2015).

In this study, we used leaf traits and molecular markers to compare the divergence between the two species in mixed (locally sympatric) and pure (locally allopatric) populations, rather than in broadly defined, large-scale sympatry and allopatry, as in a classical survey of inter-specific divergence. We specifically addressed the following questions. First, is there more divergence in molecular markers and leaf traits between the two closely related oak species in mixed stands than in pure stands? Second, if so, is this divergence symmetric? Third, what are the underlying mechanisms?

**Materials and Methods**

Materials

We sampled leaves from 778 trees in 47 natural populations
throughout the entire geographic ranges of the two selected
deciduous oak species in mainland China (Fig. 3). We prospected
each site during about an hour, covering about 3–5 km and using
morphological traits evaluated in the field as a first attempt to
identify species. Sites in which we detected individuals of both
species were considered sympatric whereas those in which we
detected only one species were considered allopatric. Using this
procedure, there were eight sympatric and 39 allopatric popula-
tions. In sympatry, we collected on average 47% Q. aliena and
53% Q. dentata (with proportions varying from 30–70% to 66–
34%). Subsequently, in the laboratory, we reassessed the status of
each stand using genetic and morphological (GMM) approaches
(Supporting Information Table S1). During sampling, we walked
randomly and ensured the individuals were located at least 10 m
apart from each other. In sympatry, the two species often grow
together, sometimes at close distance. We collected four to five
mature leaves along the four cardinal directions in the middle
layer of the canopy for leaf trait analyses, and one to two young
leaves or new branches for DNA isolation. We dried all leaf sam-
ples in silica gel immediately and we recorded the latitude, longi-
tude and altitude of each sampling site using a 621sc global
positioning system (GPS) device (Garmin, Beijing, China).

<figure><img src="image_6.png" /><figcaption>Fig. 3 Geographical distribution and sampling location of *Quercus aliena* (a) and *Quercus dentata* (b) in China. Green colors indicate the species range. For the details of sampling locations see Supporting Information Table S1.</figcaption></figure>
Microsatellite genotyping

Total genomic DNA was isolated from 15 to 20 mg samples of leaf tissue or cambium using a Plant Genomic DNA Extraction Kit (Tiangen, Beijing, China) according to the manufacturer's instructions. We first checked the quality of the DNA samples electrophoretically, using a 1% agarose gel. Then we measured their DNA concentrations with an ultra-microspectrophotometer (Thermo Fisher, Waltham, MA, USA). For genotyping, we screened 25 nuclear simple sequence repeat (nrSSR) primers used in studies of other oak species (Ueno et al., 2008; Ueno & Tsumura, 2008; Durand et al., 2010; Lyu et al., 2018) (see Table S2 for primer details). We excluded loci harboring null alleles as identified using MICRO-CHECKER v.2.2 (Oosterhout et al., 2004). Finally, we used 12 polymorphic loci for genotyping all 778 individuals. For details of the PCR mixture and amplification conditions, see Lyu et al. (2018). We analyzed the PCR products using an ABI Prism 3730 Genetic Analyzer (Applied Biosystems, Foster City, CA, USA). We scored alleles using GENEMARKER v.2.2 (Softgenetics LLC, State College, PA, USA) and checked twice the genotypes.

Landmark configuration and specific leaf area estimation

For leaf GMM analyses, we flattened, dried and scanned mature
leaves at 600 dpi resolution using a CanoScan 5600 F scanner
(Canon, Tokyo, Japan). In total, we used 3391 leaves from 683
individuals (four or five leaves per individual). We used the
scanned images to extract landmarks and leaf contours with
IMAGEJ v.1.5 (Abràmoff et al., 2005). For each leaf, we selected
13 landmarks, focusing on the primary veins, symmetric struc-
tures, and homologous parts. Three landmarks were distributed
along the middle axis of leaves (LM 1–LM 3) and nine others
were symmetrically distributed in other parts of the leaves (LM
4–LM 13) (Viscosi et al., 2009; Savriama & Klingenberg, 2011;
Viscosi & Fortini, 2011; Liu et al., 2018) (Fig. 4a). We organized
the raw data for all leaf landmark configurations into 13 pairs of
Cartesian coordinates (x, y) using IMAGEJ v.1.5. Then we
imported all the xy coordinates as input data into MORPHOJ

program. We measured leaf area and SLA in a subset of the popu-
lations (Table S3). We estimated the area of each fully extended
leaf using IMAGEJ v.1.5 software. We calculated SLA (leaf area per
unit of dry leaf weight; cm² g⁻¹) as the ratio of leaf area to dry
leaf mass.

Tree phenology

We selected one sympatric site on Shangfang Mountain located in the southwest of Beijing Municipality, China (115°48'E, 39°39'N) to observe phenological characters of the two oak species. In this site, *Q. dentata* is growing in open conditions on the top of the hills, whereas *Q. aliena* is growing in the bottom and on the hillside, at higher density, resulting in straighter trees than *Q. dentata*. We recorded bud burst, leaf-expansion, flowering, and fruit stage periods of 15 *Q. aliena* and 10 *Q. dentata* trees, as described by Crawley & Akhteruzzaman (1988). During the 2018 flowering period, we observed nearly daily almost all flowers on each tree and recorded the duration of pollen release. We used SIGMAPLOT v.12.5 (SYSTAT Software GmbH, Erkrath, Germany) to illustrate the flowering phenology of the two oak species.

Data analysis

**Genetic diversity, population cluster analysis and gene flow** We estimated the average expected heterozygosity (H<sub>E</sub), observed heterozygosity (H<sub>O</sub>) and fixation index (F) for each population across all 12 nrSSR loci with GENALEX v.6.5 (Peakall & Smouse, 2012). We then evaluated the significance of the differences in diversity between the allopatric and sympatric populations with two-group Mann-Whitney *U* tests implemented in SPSS v.22 (SPSS Inc., Chicago, IL, USA).

We estimated genetic differentiation among populations and between species of sympatric and allopatric populations using hierarchical analysis of molecular variance (AMOVA), implemented in ARLEQUIN v.3.5 (Excoffier & Lischer, 2010). We then evaluated the significance of genetic differentiation using 10 000 permutations in ARLEQUIN v.3.5.

<figure><img src="image_7.png" /></figure>

**Fig. 4** Landmark configuration of leaves (a): 1, junction of the petiole and branch; 2, junction of the blade and petiole; 3, apex of the leaf; 4, 9, the first sinus of the right-hand and left-hand leaf apex, respectively; 5, 10, the first lobe immediately above the apex of the right-hand and left-hand leaf apex, respectively; 6, 11, tip of the lobe at the largest width of the right-hand and left-hand leaf, respectively; 7, 12, the sinus immediately above the lobes of 6 and 11, respectively; 8, 13, the first basal lobe of the right-hand and left-hand leaf starting from the petiole, respectively. Plot of leaf shape (b) of *Quercus aliena* and *Quercus dentata* in sympatry vs allopatry. The box in each boxplot shows the lower, median and upper quartile values, whiskers indicate the ranges of the variation in each species. *P* < 0.05 according to a linear mixed-effect model.
We assigned individuals to homogenous clusters, without con-
sideration of sampling information, by model-based Bayesian
clustering analysis using STRUCTURE v.2.3.4 (Pritchard
et al., 2000). We used 20 independent replicates for each value of
K (1–10), where K is the (unknown) number of genetic clusters,
and 100 000 Markov Chain Monte Carlo (MCMC) repetitions
followed by 100 000 burn-in iterations for each run. We subse-
quently estimated the most likely number of clusters with Pr(X|K)
and ΔK implemented in the program STRUCTURE HARVESTER
(Earl & VonHoldt, 2012), as described by Evanno et al. (2005).
We used admixture coefficient (Q) values to distinguish whether
the sampled individuals were purebreds or hybrids with a thresh-
old value of 0.9, as suggested in other oak population studies
(Lepais et al., 2009; Peñaloza-Ramírez et al., 2010). To visualize
the genetic clusters, we used DISTRUCT software (Rosen-
berg, 2004) and plotted the admixture proportion of each indi-
vidual on a map based on the best K value. We also performed a
single principal coordinate analysis (PCoA) with all genotypes
together based on the genetic distance matrix with GENALEX
v.6.5 (Peakall & Smouse, 2012) and then plotted separately the
sympatric and allopatric individuals to visualize individuals'
genetic proximities. In addition, we used a new threshold-free
approach based on the results of the PCoA between species in
sympatry and allopatry to extract information on inter-specific
divergence using the *MIXSMSN* package in R (Prates et al., 2013).
This allows us to measure the degree of overlap of the two species
in sympatry vs allopatry.

We also applied a Bayesian approach, based on the continuous
Brownian motion model, to estimate directions of historical gene
flow between species in both sympatry and allopatry. Before anal-
ysis, we distributed all individuals in four groups, corresponding
to *Q. aliena* and *Q. dentata* in sympatry or in allopatry. Species
assignments were based on tree identification in the field. We cal-
culated the parameters *θ* (four times effective population size
multiplied by mutation rate per site per generation) and *M* (mi-
gration rate divided by the mutation rate) with MIGRATE-N
v.4.4.3 (Beerli & Felsenstein, 2001; Beerli, 2006). Each run had
three long chains, with 10 000 recorded genealogies, a sampling
increment of 20, and the first 10 000 genealogies discarded as
burn-in, under a constant mutation model. This was followed by
a static heating with default temperatures. We assessed the mode
and 95% highest posterior density after checking for data conver-
gence. We selected the best model by comparing the marginal
likelihoods from three models with different gene flow directions.
To estimate asymmetrical gene flow, we used the posterior distri-
bution estimates to compute the mean and standard deviation of
mutation-scaled migration rate (*M*) for *Q. aliena* → *Q. dentata*
and *Q. dentata* → *Q. aliena* in sympatry and in allopatry
(Beerli, 2006). Differences among these parameters were tested
by analyses of variance (ANOVAs) followed by *post hoc* analyses
using Tukey's HSD (honestly significant difference) test for mul-
tiple comparisons implemented in SPSS v.22 (SPSS Inc.).

**Morphological analysis** We first performed a generalized pro-
crustes analysis (GPA) (Rohlf & Slice, 1990) to minimize differ-
ences between landmark configurations by translation, scaling

and rotation to maximize the coincidence of leaf coordinate data
using the MORPHOJ program (Klingenberg, 2011). Next, we took
the mean of the landmarks across four to five leaves of the same
tree as suggested by Klingenberg (2011). Finally, we exported the
normalized matrix of 26 procrustes coordinates of mean land-
marks per tree for leaf shape analyses (Viscosi et al., 2009). Using
mean landmark values per tree rather than individual leaf land-
mark data does not result in significant loss of information in
subsequent analyses, as shown by Viscosi & Cardini (2011). Dis-
criminant analysis is considered the most powerful statistical
method for investigating taxonomic differences between groups
(Klingenberg, 2011). In this study, we subjected procrustes coor-
dinates of the 683 individuals to linear discriminant analysis
(LDA), using the MASS package in R v.4.0.2 (Venables & Rip-
ley, 2002), to determine the proportion of purebreds and hybrids
in each population, based on the resulting linear discriminant
(LD) scores (Albarrán-Lara et al., 2018), and graphically visual-
ized the results on a map. In addition, we performed a principal
component analysis (PCA) and displayed the distribution fre-
quency of principal component (PC) scores for all individuals to
visualize the differences between species in leaf shape, using the
VEGAN package in R (Oksanen et al., 2022). We also used the
threshold-free approach on the results of PC1 for leaf shape
between species in sympatry and allopatry using the MIXSMSN
package in R, as described earlier for genetic markers (Prates
et al., 2013). We further compared the differences in leaf shape,
leaf area and SLA between species in sympatry and allopatry
using the GGPLOT2 R package (Wickham, 2016).

**Patterns of inter-specific phenotypic divergence** To character-
ize leaf traits' climatic responses, we first extracted 19 bioclimatic
variables at 30 s (c. 1 km²) resolution from WORLDCLIM v.2 raster
layers (Fick & Hijmans, 2017). We transformed these variables
using z-scores set to increase with temperature and precipitation
(Table S4), using a PCA in the VEGAN R package (Oksanen
et al., 2022). We then studied SLA variation for all, sympatric
and allopatric populations as a function of climatic factors (pre-
cipitation seasonality, annual precipitation, temperature seasonal-
ity and annual mean temperature). We also estimated the
divergence of leaf shape by GMMs between species in sympatry
and in allopatry using data acquired from all sampled individuals
or only from the genetically and morphologically purebred indi-
viduals. This was performed both across the whole climatic range
(the climatic range of all studied populations) and across the cli-
matic range common to both sympatric and allopatric popula-
tions, i.e. excluding allopatric populations experiencing divergent
climates not encountered by sympatric populations. Restricting
comparisons between sympatric and allopatric populations expe-
riencing the same climatic range should make the results more
comparable. We estimated mean leaf shapes of the two species in
the sympatric and allopatric populations (Dsym and Dallo, respec-
tively), as well as the difference between them (Dsym-dallo) as
described by Adams (2004). A greater divergence in sympatry
than in allopatry yields a positive Dsym-dallo value. Second, we used
a linear mixed model to detect effects of species type, community
type (sympatric or allopatric), and their interaction on each trait,
as suggested by Adams (2004) and Bates *et al.* (2015), using the LME4 package in R v.4.0.2. We assessed the significance of the effects using an ANOVA test implemented in the CAR R package (Fox *et al.*, 2013). A significant interaction between species type and community implies species-specific trait divergence (Adams, 2004). In addition, we used linear mixed modeling to test the hypotheses (1) that the trait's environmental response differs in sympatry and in allopatry and (2) that the divergence between the two species in sympatry and allopatry depends on environmental factors, focusing on purebred individuals (either morphologically purebred and genetically purebred). We also assessed the significance of the differences by ANOVA.

## Results

### Genetic diversity and divergence

**Genetic diversity, differentiation and gene flow** All 12 nrSSR loci were polymorphic, with 5–28 alleles per locus. The genetic diversity was lower in sympatric populations ($H_O = 0.58$, $H_E = 0.63$) than in allopatric populations ($H_O = 0.63$, $H_E = 0.64$) ($P<0.05$) and heterozygote deficit much higher in the former than in the latter (Table S5). We also found that most of the genetic variation resided within populations (86% and 88% in the sympatric and allopatric populations, respectively) (Table S6). The fixation index ($F_{ST}$) values between species were higher for sympatric populations than for allopatric populations ($P<0.001$). MIGRATE-N analysis based on Bayesian factor values indicated that asymmetric historical gene flow had occurred between the two closely related oaks, mainly from *Q. dentata* to *Q. aliena* (63.7 and 48.7, respectively) (Table 1). The mutation-scaled migration rates ($M$) was higher from *Q. dentata* to *Q. aliena* (Fig. S1b) than in the reverse direction (Fig. S1a; Table S7), significantly so in sympatry ($P=0.05$, Table S7). Furthermore, we inferred that the rate of historical gene flow between species has been significantly lower in sympatric than in allopatric populations ($P<0.05$, Table S7; Fig. S1).

<table><caption>**Table 1** Historical gene flow as estimated by MIGRATE-N between Quercus aliena and Quercus dentata in sympatry vs in allopatry.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">0</th><th colspan="2">N<sub>e</sub> M</th></tr><tr><th>Q. aliena →</th><th>Q. dentata →</th></tr></thead><tbody><tr><td><strong>Sympatry</strong></td><td></td><td></td><td></td></tr><tr><td>Q. aliena</td><td><strong>2.6</strong> [0.7–4.7]</td><td></td><td><strong>49.9</strong> [40.0–51.2]</td></tr><tr><td>Q. dentata</td><td><strong>2.9</strong> [1.1–3.6]</td><td><strong>34.3</strong> [15.9–43.7]</td><td></td></tr><tr><td><strong>Allopatry</strong></td><td></td><td></td><td></td></tr><tr><td>Q. aliena</td><td><strong>3.1</strong> [1.6–4.4]</td><td></td><td><strong>67.8</strong> [55.4–81.8]</td></tr><tr><td>Q. dentata</td><td><strong>3.4</strong> [1.9–4.7]</td><td><strong>58.9</strong> [46.2–66.1]</td><td></td></tr><tr><td><strong>Total</strong></td><td></td><td></td><td></td></tr><tr><td>Q. aliena</td><td><strong>3.2</strong> [1.6–4.5]</td><td></td><td><strong>63.7</strong> [45.5–79.4]</td></tr><tr><td>Q. dentata</td><td><strong>3.3</strong> [1.5–4.8]</td><td><strong>48.7</strong> [45.6–50.4]</td><td></td></tr></tbody></table>

The mode of the posterior distribution is shown in bold and the values in square brackets give the 95% credibility interval; 0, 4× effective population size × mutation rate per site per generation; $N_eM$, effective number of migrants per generation; →, source populations.

**Genetic structure** The Pr($X|K$) and $\Delta K$ statistics obtained from Bayesian clustering strongly support the presence of two major clusters in the dataset, one corresponding to *Q. aliena* and the other to *Q. dentata* (Fig. S2). In addition, analysis with the threshold value $Q$ set to 0.9 indicated that there is a slightly lower percentage of admixed individuals in sympatric than in allopatric populations (29 and 33%, respectively, Table S1). The geographical distribution of genetic clusters also revealed a higher level of inter-specific genetic divergence in the sympatric than in the allopatric populations (Fig. 5; Table S1). PCoA results based on the genetic distance matrix at the individual level were largely concordant with the STRUCTURE analysis. They point to significant genetic differentiation between *Q. aliena* and *Q. dentata* in sympatry, with clear separation of the species along PC1 and small proportions of admixed individuals and a bimodal frequency distribution of individuals (Fig. 5c). In contrast, in allopatric populations *Q. aliena* and *Q. dentata* purebreds and admixed individuals were intermingled with each other to some extent, and there was a unimodal frequency distribution of individuals (Fig. 5d).

### Patterns of phenotypic divergence

**Leaf morphological variation** As recommended by Klingenberg (2011), we manually removed 28 leaf outliers, i.e. specimens deviating strongly from the average (Table S1 column H-I). The final dataset included 3363 leaves corresponding to 99.2% of the sampled leaves. For *Q. aliena*, we found a clear difference between mean landmark values in sympatry and in allopatry. For *Q. dentata*, the difference was not obvious (Fig. S3). The morphological data provided stronger discrimination between species than the genetic markers: 87% of the individuals were ‘pure’ according to leaf shape GMM analysis (with LD scores < −1 or > 1), compared to only 66% when genetic markers were used for assigning species ($Q$ scores < 0.1 or > 0.9) (Figs 5, 6; Table S1). Overall, morphological and genetic taxonomic criteria matched for 401 (65%) of the individuals (Table S8). Leaf shape PC1 accounted for 47% of the total leaf shape variance, so we used it as a proxy for leaf shape variation.

Conclusions regarding the allopatric or sympatric status of each population based on either morphological (GMM) or genetic analyses were the same as those arrived at in the field. However, hybrids could only be identified in the laboratory using dedicated approaches (GMMs or genotyping) followed by statistical analyses. We provide details on field identification, genetic and phenotypic assignment for each individual and population in Table S1.

We detected greater differences in all leaf traits (leaf shape by GMMs; leaf area and SLA) between species in sympatry than in allopatry, mostly because *Q. aliena* was more divergent from *Q. dentata* in sympatric than in allopatric populations (Figs 4b, S4). LDA using data acquired for all individuals detected significant differences in leaf shape between *Q. aliena* ($LD < -1$, $N=242$) and *Q. dentata* ($LD > 1$, $N=348$) and smaller proportions of admixed individuals in the sympatric than in the allopatric populations (8 and 16%, respectively; Fig. 6a,b;
<figure><img src="image_8.png" /><figcaption>**Fig. 5** Genetic differentiation between *Quercus aliena* and *Quercus dentata* in sympatry and allopatry. The geographic display of the STRUCTURE results with *K* = 2 at population level for sympatric (a) and allopatric populations (b). For the details of sampling locations see Supporting Information Table S1. Principal component analysis (PCA) plot at individual level in sympatry (c) and allopatry (d), with the distribution frequency and box plots of the first principal component (PC1) and second principal component (PC2) values plotted on the top and right sides of the scatter plot. The box in each boxplot shows the lower, median and upper quartile values, whiskers indicate the ranges of the variation. Percentage of total variance explained by each axis are noted in brackets.</figcaption></figure>

Table S1). Moreover, in some allopatric populations (KA, SC, ZS and ZY) more than 50% of the individuals were admixed (Fig. 6b). Clearly, oak trees' leaf morphology is more bimodal in sympatric than in allopatric populations (Fig. 6c,d). In addition, the degree of overlap of the two species in sympatry is significantly lower than in allopatry for both types of markers, and phenotypic divergence significantly greater than molecular divergence in sympatry in the threshold-free analysis (Fig. S5; Table S9).

The PC1 of the 19 extracted bioclimatic variables accounted for 53% of the total climatic variance in the species' total range. There were used as a proxy for climate (Table S10). Leaf traits of *Q. dentata* appear to be quite stable along this climatic gradient. In contrast, we detected significant correlation between *Q. aliena* leaf shape or SLA and climate (Figs 7, S6). In allopatry, leaf shape variation was correlated with climatic factors in both species, with especially large divergence for *Q. aliena* in more humid and warmer environments (Figs 7c, S6c). However, in sympatry, trait divergence did not depend on the environment (Figs 7b, S6b). In

addition, we found that SLA of *Q. aliena* but not of *Q. dentata* decreased when precipitation seasonality increased (Fig. S7).

We found positive *D*<sub>sym-allo</sub> values. This indicates greater morphological divergence in sympatry than in allopatry, regardless of whether or not genetically and morphologically admixed individuals are included in the analysis. After controlling for climate, by restricting the comparison to the climatic range common to both sympatric and allopatric populations (i.e. excluding allopatric populations growing under climatic conditions not experienced by sympatric populations), we still found higher divergence in sympatry. However, the support was weaker when we considered only genetically purebreds (Table 2). Similarly, we detected significant species × community interaction in all but one of the six cases investigated (see Table 2), indicating that leaf shape response in sympatry and in allopatry differs between the two species.

**Flowering phenology** Both *Q. aliena* and *Q. dentata* produced male flowers a week earlier than female flowers (Fig. S8), and
<figure><img src="image_9.png" /><figcaption>**Fig. 6** Leaf geometric morphometric analysis. Geographic display of leaf shape by linear discriminant analysis (LDA) in sympatric (a) and allopatric (b) populations. For the details of sampling locations see Supporting Information Table S1. Principal component analysis (PCA) plot of leaf shape in sympatry (c) and allopatry (d) at individual level, with the distribution frequency and box plots of the first principal component (PC1) and second principal component (PC2) values plotted on the top and right sides of the scatter plot. The box in each boxplot shows the lower, median, upper quartile and the outlier values, whiskers indicate the ranges of the variation. Percentage of total variance explained by each axis are noted in brackets.</figcaption></figure>

*Q. aliena* trees flowered about a week earlier than *Q. dentata* trees (Table S11). Hence, female *Q. aliena* flowers were potentially exposed to pollen from *Q. dentata* at the monitoring site in the Shangfang Mountain in the studied year (Fig. S8).

Discussion

In this study, we proposed a new ECD approach based on the
composition of local communities, by comparing plant traits
measured in mixed-species stands or in pure stands. We there-
fore tested for ECD by focusing on whether species can interact
locally or not. Considering the sessile nature of plant species,
this local approach makes things more comparable than in the
classical ECD approach focusing on broad geographic scales.
We also accounted for the effect of climate on the

corresponding traits and genetically characterized all the trees
studied to confirm species assignment and infer long-term inter-
specific gene flow.

Asymmetric ecological character displacement in sympatry

Morphological data allowed better species delimitation, with
fewer intermediates, than the genetic data. We used a direct test
for pattern associated with character displacement by calculating
the difference between species leaf trait means in sympatry and
allopatry (Dsym-allo), as suggested by Adams (2004). Despite
biomechanical constraints, leaf shape is a functionally significant
trait (Givnish & Kriebel, 2017), because it affects rates and other
parameters of photosynthesis, which play crucial roles in plant
growth, development and survival (Nicotra et al., 2011). The
<figure><img src="image_10.png" /><figcaption>**Fig. 7** Leaf shape variation and correlation coefficient (R) and significance (P) value in all (a), sympatric (b) and allopatric (c) populations as a function of climatic data. Morphologically pure individuals in sympatry (175 individuals, 93 *Quercus aliena* and 82 *Quercus dentata*) and allopatry (415 individuals, 149 *Q. aliena* and 266 *Q. dentata*) were used for the analysis. The P-values for each graph were calculated using a linear mixed-effect model.</figcaption></figure>

pattern of morphological evolution of *Q. aliena* in the presence of *Q. dentata* is towards more lobes per leaf, a phenotype previously associated with occurrence in drier sites (Cavender-Bares, 2019). We also found an evolution of *Q. aliena* in sympatry towards smaller leaves and lower SLA. These features have been previously shown to be associated with a more resource-conservative strategy (for either nutrients, water or light; Wright *et al.*, 2004; Cavender-Bares & Pahlich, 2009; Cavender-Bares *et al.*, 2015; Arenas-Navarro *et al.*, 2020), suggesting that *Q. aliena* experiences more competition in sympatry.

The greater morphological divergence in sympatry was mostly due to *Q. aliena* trees having more divergent trait values in sympatry than in allopatry (individuals of the other species, *Q. dentata*, had similar trait values in sympatry and allopatry). Models of ecological character divergence generally consider only symmetric divergence. However, there is no reason to expect divergence to be the same for each species, especially if they differ substantially in abundance, which might impose greater selective pressures on the rarer species. Interestingly, we also found indications of more introgression in *Q. aliena* than in *Q. dentata*, based on estimates of long-term inter-specific gene flow.

## Impact of ecological variables on trait divergence

Using a variety of approaches, we found less morphologically intermediate individuals in sympatric than in allopatric populations. Yet, if the parental species can hybridize, we would expect instead more intermediate individuals in sympatry. We therefore infer that stronger selection against morphologically intermediate individuals is taking place in sympatry, resulting in ECD. A first alternative hypothesis is that the two species have responded differently to environmental factors, thereby confounding interpretations of the effects of species' coexistence on traits. However,

we found that the signal of greater divergence in sympatry than in allopatry was still present after restricting the analysis to a common climatic range, both when considering all individuals assigned to each species and when restricting the analysis to purebreds (defined using leaf morphology). A second alternative hypothesis to ECD is that disruptive selection acts on other (unknown) traits, and that the studied traits (leaf morphology, leaf area or SLA) covary with these unknown traits. This interpretation was proposed recently by Alexandre *et al.* (2020) to explain the rapid divergence in leaf morphological traits between European oak species. However, this model does not appear particularly parsimonious, especially considering that credible adaptive scenarios exist explaining the direction of observed phenotypic changes, as previously stressed.

## Origin of asymmetric ecological character displacement

Asymmetric divergence in leaf traits and asymmetric long-term inter-specific gene flow between the two species might be due to the same cause: demographic imbalance (i.e. uneven sample sizes) during colonization of the less hardy species. As *Q. aliena* is adapted to a warmer climate, we infer that *Q. dentata* first colonized its current range following climate warming in the past, with *Q. aliena* arriving later and experiencing strong selective pressures at the time of invasion when penetrating the range of the other, already abundant, species. This would have resulted in both asymmetric introgression, stronger towards the colonizing species than towards the resident species (Currat *et al.*, 2008; Du *et al.*, 2011), and asymmetric ECD for leaf traits. A trade-off between competition for resources and colonization abilities is indeed possible, as shown in a comparative study of reproductive strategies of two European oak species (Lagache *et al.*, 2014). We also note that the different flowering phenology of the two oak
<table><caption>**Table 2** Direct tests of leaf shape divergence between sympatric and allopatric populations for whole and common climatic range.</caption><thead><tr><th rowspan="3">Climatic range</th><th rowspan="3">Individual type</th><th colspan="3">Leaf shape divergence</th><th colspan="4">Statistics</th></tr><tr><th rowspan="2">D<sub>sym</sub></th><th rowspan="2">D<sub>allo</sub></th><th rowspan="2">D<sub>sym-allo</sub><sup>c</sup></th><th colspan="2">Species<sup>d</sup></th><th colspan="2">Community<sup>e</sup></th><th colspan="2">Species × Community<sup>f</sup></th></tr><tr><th>F</th><th>P</th><th>F</th><th>P</th><th>F</th><th>P</th></tr></thead><tbody><tr><td rowspan="3">Whole</td><td>All</td><td>0.15</td><td>0.1</td><td>0.05</td><td>350</td><td>&lt;0.001</td><td>5.7</td><td><strong>0.02</strong></td><td>5.7</td><td><strong>0.03</strong></td></tr><tr><td>Genetically purebred<sup>a</sup></td><td>0.15</td><td>0.13</td><td>0.02</td><td>340</td><td>&lt;0.001</td><td>4.9</td><td><strong>&lt;0.001</strong></td><td>3.6</td><td><strong>0.05</strong></td></tr><tr><td>Morphologically purebred<sup>b</sup></td><td>0.16</td><td>0.12</td><td>0.04</td><td>427</td><td>&lt;0.001</td><td>8.1</td><td><strong>0.006</strong></td><td>3.9</td><td><strong>0.04</strong></td></tr><tr><td rowspan="3">Common</td><td>All</td><td>0.15</td><td>0.11</td><td>0.04</td><td>356</td><td>&lt;0.001</td><td>6.7</td><td><strong>0.01</strong></td><td>3.1</td><td><strong>0.05</strong></td></tr><tr><td>Genetically purebred<sup>a</sup></td><td>0.14</td><td>0.12</td><td>0.02</td><td>407</td><td>&lt;0.001</td><td>1.5</td><td><strong>0.12</strong></td><td>1.8</td><td><strong>0.09</strong></td></tr><tr><td>Morphologically purebred<sup>b</sup></td><td>0.14</td><td>0.1</td><td>0.04</td><td>385</td><td>&lt;0.001</td><td>7.6</td><td><strong>0.01</strong></td><td>2.7</td><td><strong>0.02</strong></td></tr></tbody></table>

species, combined with their protandry, could facilitate asymmetric inter-specific gene flow.

## Reproductive isolation and ecological character displacement

In principle, the greater divergence observed in contact zones for both phenotypic traits and multilocus genetic markers might be due to higher postzygotic or prezygotic reproductive isolation (Jiggins & Mallet, 2000). However, we found numerous intermediate adult individuals, both in mixed and in pure stands, so the action of a postzygotic barrier (resulting in reduced hybrid viability) seems unlikely. Alternatively, reinforcement of prezygotic barriers could have caused the observed divergence in leaf traits. However, this seems also unlikely. Our flowering phenology data, obtained by observations of sympatric material, do point towards some prezygotic restriction to their reproductive interaction, as *Q. aliena* flowers earlier than *Q. dentata*. However, we found that the two species are more divergent at morphological traits than at presumably neutral molecular markers, unlike patterns often observed in other oak studies using similar approaches (e.g. Lagache *et al.*, 2014). This suggests that divergent selection acting on trees' phenotypic traits (i.e. ECD) rather than reinforcement of prezygotic barriers is the driving process. We also note that reinforcement, if present, is more likely to be a consequence than a cause of ECD. Indeed, the driving force of reinforcement is disruptive selection (Merrill *et al.*, 2012; Hopkins, 2013). Overall, we argue that, in these two oak species, the most plausible explanation for asymmetric divergence in leaf traits in sympatry is biased competition resulting in asymmetric selection. We think that, in the future, greater consideration for patterns of asymmetric character displacement triggered by inter-specific competition will benefit plant evolutionary studies.

## Acknowledgements

The authors thank Yang Xu, Wei Su, Ping liu from Beijing Forestry University (BFU) for leaf sampling, landmarking and flowering phenology record; Dr Sun Han from BFU for discussion on SLA analysis. Prof. Weining Bai from Beijing Normal University for gene flow analysis. The authors thank the Associate Editor, Prof. David Ackerly and two anonymous reviewers for their helpful comments on the manuscript. This research was supported by a grant from the National Science Foundation of China (no. 42071060) to FKD.

## Competing interests

None declared.

## Author contributions

FKD designed the research; MQ performed the experiments; MQ, Y-YZ and FKD performed the analyses, with input from RJP; FKD, MQ, and RJP wrote the manuscript and all authors contributed to its revision.

<footnote><sup>a</sup>Genetically purebred: excluding admixed trees with Q values, 0.1–0.9, <sup>b</sup>Morphologically purebred, excluding admixed trees with linear discriminant (LD) values: –1 to 1. <sup>c</sup>D<sub>sym-allo</sub>, positive values are indicative of patterns associated with greater divergence of trait in sympatry than in allopatry; <sup>d</sup>Species (Quercus aliena vs Quercus dentata), <sup>e</sup>Community (allopatry vs sympatry), and <sup>f</sup>Species × Community interaction effects via linear mixed-effect models similar to Adams & Collyer (2007). The P and F values were calculated by the linear mixed-effect models. Bold values indicate statistical significance at P < 0.05.</footnote>
ORCID

Fang K. Du ID https://orcid.org/0000-0002-7377-5259
Rémy J. Petit ID https://orcid.org/0000-0002-4809-1453
Min Qi ID https://orcid.org/0000-0002-6570-5795
Yuan-Ye Zhang ID https://orcid.org/0000-0002-8636-2242

Data availability

Photographs of sampling sites can be obtained in https://www.oakofchina.org/photo-of-sampling/, and leaf morphological data can be obtained in https://doi.org/10.6084/m9.fishshare.20000021. Genotyping data can be found in figshare: https://doi.org/10.6084/m9.fishshare.16821439.

References

Abadie P, Roussel G, Dencausse B, Bonnet C, Garnier-Géré P. 2012. Strength, diversity and plasticity of postmating reproductive barriers between two hybridizing oak species (*Quercus robur* L. and *Quercus petraea* (Matt) Liebl.). *Journal of Evolutionary Biology* 25: 157-173.

Abrámoff MD, Magalhães PJ, Paulo J, Ram SJ. 2005. Image processing with IMAGEJ Part II. *Biophotonics International* 11: 36-42.

Adams DC. 2004. Character displacement via aggressive interference in Appalachian salamanders. *Ecology* 85: 2664-2670.

Adams DC, Collyer ML. 2007. Analysis of character divergence along environmental gradients and other covariates. *Evolution* 61: 510-515.

Adams DC, Rohlf FJ. 2000. Ecological character displacement in *Plethodon*: biomechanical differences found from a geometric morphometric study. *Proceedings of the National Academy of Sciences, USA* 97: 4106-4111.

Albarrán-Lara AL, Petit RJ, Kremer A, Caron H, Peñaloza-Ramírez JM, Gugger PF, Dávila-Aranda PD, Oyama K. 2018. Low genetic differentiation between two morphologically and ecologically distinct giant-leaved Mexican oaks. *Plant Systematics and Evolution* 305: 89-101.

Alexandre H, Truffaut L, Klein E, Ducoussou A, Chancerel E, Lesur I, Dencausse B, Louvet JM, Nepveu G, Torres-Ruiz JM et al. 2020. How does contemporary selection shape oak phenotypes? *Evolutionary Applications* 13: 2772-2790.

Arenas-Navarro M, García-Oliva F, Terrazas T, Torres-Miranda A, Oyama K. 2020. Leaf habit and stem hydraulic traits determine functional segregation of multiple oak species along a water availability gradient. *Forests* 11: 894.

Bates D, Mächler M, Bolker B, Walker S. 2015. Fitting Linear Mixed-Effects Models using LME4. *Journal of Statistical Software* 67: 1-48.

Beans CM. 2014. The case for character displacement in plants. *Ecology and Evolution* 4: 862-875.

Beerli P. 2006. Comparison of Bayesian and maximum-likelihood inference of population genetic parameters. *Bioinformatics* 22: 341-345.

Beerli P, Felsenstein J. 2001. Maximum likelihood estimation of a migration matrix and effective population sizes in n subpopulations by using a coalescent approach. *Proceedings of the National Academy of Sciences, USA* 98: 4563-4568.

Brown WL, Wilson EO. 1956. Character displacement. *Systematic Zoology* 2: 49-64.

Butlin RK, Galindo J, Grahame JW. 2008. Sympatric, parapatric or allopatric: the most important way to classify speciation? *Philosophical Transactions of the Royal Society of London. Series B: Biological Sciences* 363: 2997-3007.

Cannon CH, Petit RJ. 2019. The oak syngameon: more than the sum of its parts. *New Phytologist* 226: 978-983.

Cavender-Bares J. 2019. Diversification, adaptation, and community assembly of the American oaks (*Quercus*), a model clade for integrating ecology and evolution. *New Phytologist* 221: 669-692.

Cavender-Bares J, González-Rodríguez A, Eaton DAR, Hipp AAL, Beulke A, Manos PS. 2015. Phylogeny and biogeography of the American live oaks

(Quercus subsection *Virentes*): a genomic and population genetics approach.
*Molecular Ecology* 24: 3668-3687.

Cavender-Bares J, Kitajima K, Bazzaz FA. 2004. Multiple trait associations in relation to habitat differentiation among 17 Floridian oak species. *Ecological Monographs* 74: 635-662.

Cavender-Bares J, Pahlich A. 2009. Molecular, morphological, and ecological niche differentiation of sympatric sister oak species, *Quercus virginiana* and *Q. geminata* (Fagaceae). *American Journal of Botany* 96: 1690-1702.

Crawley MJ, Akhteruzzaman M. 1988. Individual variation in the phenology of oak trees and its consequences for herbivorous insects. *Functional Ecology* 2: 409-415.

Currat M, Ruedi M, Petit RJ, Excoffier L. 2008. The hidden side of invasion: Massive introgression by local genes. *Evolution* 62: 1908-1920.

Darwin C. 1872. On the origin of species: by means of natural selection, or preservation of favoured races in the struggle for life. *Nature* 144: 971-972.

Denk T, Grimm GW, Manos PS, Deng M, Hipp AL. 2018. An updated infrageneric classification of the oaks: review of previous taxonomic schemes and synthesis of evolutionary patterns. In: Gil-Pelegrín E, Peguero-Pina J, Sancho-Knapik D, eds. *Oaks physiological ecology. Exploring the functional diversity of genus Quercus*. Cham, Switzerland: Springer, 13-38.

Dryden IL, Mardia KV. 1998. *Statistical shape analysis*. *Wiley series in probability and statistics Chichester*. New York, NY, USA: John Wiley & Sons.

Du F, Peng XL, Liu JQ, Lascoux M, Hu FS, Petit RJ. 2011. Direction and extent of organelle DNA introgression between two spruce species in the Qinghai-Tibetan Plateau. *New Phytologist* 192: 1024-1033.

Durand J, Bodénès C, Chancerel E, Frigeriol JM, Vendramin G, Sebastiani F, Buonamici A, Gailing O, Koelewijn HP, Villani F. 2010. A fast and cost-effective approach to develop and map EST-SSR markers: Oak as a case study. *BMC Genomics* 11: 570.

Earl DA, VonHoldt BM. 2012. STRUCTURE HARVESTER: A website and program for visualizing STRUCTURE output and implementing the Evanno method. *Conservation Genetics Resources* 4: 359-361.

Evanno G, Regnaut S, Goudet J. 2005. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. *Molecular Ecology* 14: 2611-2620.

Excoffier L, Lischer HEL. 2010. ARLEQUIN suite v.3.5: A new series of programs to perform population genetics analyses under Linux and Windows. *Molecular Ecology Resources* 10: 564-567.

Fick SE, Hijmans RJ. 2017. WORLDCLIM 2: New 1-km spatial resolution climate surfaces for global land areas. *International Journal of Climatology* 37: 4302-4315.

Fox J, Friendly M, Weisberg S. 2013. Hypothesis tests for multivariate linear models using the *CAR* package. *R Journal* 5: 39-52.

Fujiwara K, Harada A. 2015. *Warm-temperate deciduous forests around the northern hemisphere*. Cham, Switzerland: Springer.

Gailing O, Zhang R. 2018. Experimental evidence for selection against hybrids between two interfertile red oak species. *Silvae Genetica* 67: 106-110.

Germain RM, Williams JL, Schlüter D, Angert AL. 2018. Moving character displacement beyond characters using contemporary coexistence theory. *Trends in Ecology & Evolution* 33: 74-84.

Givnish TJ, Kriebel R. 2017. Causes of ecological gradients in leaf margin entirety: Evaluating the roles of biomechanics, hydraulics, vein geometry, and bud packing. *American Journal of Botany* 104: 354-366.

Goldberg EE, Lande R. 2006. Ecological and reproductive character displacement of an environmental gradient. *Evolution* 60: 1344-1357.

Grant V. 1966. Selective origin of incompatibility barriers in the plant genus *Gilia*. *American Naturalist* 100: 99-118.

Hipp AL, Manos PS, González-Rodríguez A, Hahn M, Kaproth M, McVay JD, Avalos SV, Cavender-Bares J. 2018. Sympatric parallel diversification of major oak clades in the Americas and the origins of Mexican species diversity. *New Phytologist* 217: 439-452.

Hipp AL, Manos PS, Hahn M, Avishai M, Bodénès C, Cavender-Bares J, Crowl AA, Deng M, Denk T, Fitz-Gibbon S et al. 2020. Genomic landscape of the global oak phylogeny. *New Phytologist* 226: 1198-1212.

Hopkins R. 2013. Reinforcement in plants. *New Phytologist* 197: 1095-1103.

Huang CJ, Zhang YT, Bartholomew B. 1999. Fagaceae. In: Wu ZY, Raven PH, eds. *Flora of China, vol. 4* [in English]. Beijing, China: Science Press, 370-380.
Hubert F, Grimm GW, Jousselin E, Berry V, Franc A, Kremer A. 2014.

Multiple nuclear genes stabilize the phylogenetic backbone of the genus *Quercus*. *Systematics and Biodiversity* **12**: 405–423.

Jiggins CD, Mallet J. 2000. Bimodal hybrid zones and speciation. *Trends in Ecology & Evolution* **15**: 250–255.

Klingenberg CP. 2011. MORPHO: An integrated software package for geometric morphometrics. *Molecular Ecology Resources* **11**: 353–357.

Kusi J, Karsai I. 2019. Plastic leaf morphology in three species of *Quercus*. The more exposed leaves are smaller, more lobated and denser. *Plant Species Biology* **35**: 24–37.

Lagache L, Klein EK, Ducoussou A, Petit RJ. 2014. Distinct male reproductive strategies in two closely related oak species. *Molecular Ecology* **23**: 4331–4343.

Lepais O, Petit RG, Guichoux E, Lavabre JE, Alberto F, Kremer A, Gerber S. 2009. Species relative abundance and direction of introgression in oaks. *Molecular Ecology* **18**: 2228–2242.

Lepais O, Roussel G, Hubert F, Kremer A, Gerber S. 2013. Strength and variability of postmating reproductive isolating barriers between four European white oak species. *Tree Genetics & Genomes* **9**: 841–853.

Li YJ, Zhang Y, Liao PC, Wang TR, Wang XY, Ueno S, Du F. 2021. Genetic, geographic, and climatic factors jointly shape leaf morphology of an alpine oak, *Quercus aquifolioides* Rehder & EH Wilson. *Annals of Forest Science* **78**: 1–18.

Liu Y, Li Y, Song J, Zhang RP, Yan Y, Wang Y, Du FK. 2018. Geometric morphometric analyses of leaf shapes in two sympatric Chinese oaks: *Quercus dentata* Thunberg and *Quercus aliena* Blume (Fagaceae). *Annals of Forest Science* **75**: 90.

Lobo A, Torres-Ruiz JM, Burlett R, Lemaire C, Parise C, Francioni C, Truffaut L, Tomášková I, Hansen JK, Kjær ED et al. 2018. Assessing inter-and intraspecific variability of xylem vulnerability to embolism in oaks. *Forest Ecology and Management* **424**: 53–61.

Lyu J, Song J, Liu Y, Wang YY, Li JQ, Du FK. 2018. Species boundaries between three sympatric oak species: *Quercus aliena*, *Q. dentata*, and *Q. variabilis* at the northern edge of their distribution in China. *Frontiers in Plant Science* **9**: 414.

Merrill RM, Wallbank RWR, Bull V, Salazar PCA, Mallet J, Stevens M, Jiggins CD. 2012. Disruptive ecological selection on a mating cue. *Proceedings of the Royal Society B: Biological Sciences* **279**: 4907–4913.

Nardini A, Pedà G, Rocca NL. 2012. Trade-offs between leaf hydraulic capacity and drought vulnerability: morpho-anatomical bases, carbon costs and ecological consequences. *New Phytologist* **196**: 788–798.

Nicota AB, Leigh A, Boyce CK, Jones CS, Niklas KJ, Royer DL, Tsukaya H. 2011. The evolution and functional significance of leaf shape in the angiosperms. *Functional Plant Biology* **38**: 535–552.

Niinemets Ü. 2015. Is there a species spectrum within the world-wide leaf economics spectrum? Major variations in leaf functional traits in the Mediterranean sclerophyll *Quercus ilex*. *New Phytologist* **205**: 79–96.

Noor M. 1999. Reinforcement and other consequences of sympatry. *Heredity* **83**: 503–508.

Oksanen J, Blanchet FG, Friendly M, Kindt R, Legendre P, McGlinn D, Minchin PR, O'Hara RB, Simpson GL, Solymos P et al. 2022. VEGAN: Community ecology package. R package v.2.5-7. [WWW document] URL https://CRAN.R-project.org/package=vegan [accessed 3 June 2022].

Oosterhout CV, Hutchinson WF, Wills DPM, Shipley P. 2004. MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. *Molecular Ecology Notes* **4**: 535–538.

Peakall R, Smouse PE. 2012. GENELEX 6: Genetic analysis in Excel. Population genetic software for teaching and research—an update. *Bioinformatics* **28**: 2537–2539.

Peguero-Pina JJ, Sisó S, Sancho-Knapik D, Díaz-Espejo A, Flexas J, Galmés J, Gil-Pelegrín E. 2015. Leaf morphological and physiological adaptations of a deciduous oak (*Quercus faginea* Lam.) to the Mediterranean climate: a comparison with a closely related temperate species (*Quercus robur* L.). *Tree Physiology* **36**: 287–299.

Penalosa-Ramírez JM, González-Rodríguez A, Mendoza-Cuenca L, Caron H, Kremer A, Oyama K. 2010. Interspecific gene flow in a multispecies oak hybrid zone in the Sierra Tarahumara of Mexico. *Annals of Botany* **105**: 389–399.

Prates MO, Lachos VH, Cabral CRB. 2013. MIXSMSN: Fitting finite mixture of scale mixture of skew-normal distributions. *Journal of Statistical Software* **54**: 1–20.

Pritchard JK, Stephens M, Donnelly P. 2000. Inference of population structure using multilocus genotype data. *Genetics* **155**: 945–959.

Reifová R, Reif J, Antczak M, Nachman MW. 2011. Ecological character displacement in the face of gene flow: evidence from two species of *nightingales*. *BMC Evolutionary Biology* **11**: 1–11.

Rieseberg LH, Wood TE, Baack EJ. 2006. The nature of plant species. *Nature* **440**: 524–527.

Rohlf FJ, Slice DE. 1990. Extensions of the Procrustes method for the optimal superimposition of landmarks. *Systematic Zoology* **39**: 40–59.

Rosenberg NA. 2004. DISTRUCT: A program for the graphical display of population structure. *Molecular Ecology Notes* **4**: 137–138.

Royer DL, McElwain JC, Adams JM, Wilf P. 2008. Sensitivity of leaf size and shape to climate within *Acer rubrum* and *Quercus kelloggii*. *New Phytologist* **179**: 808–817.

Royer DL, Moynihan KM, McKee ML, Londoño L, Franks PJ. 2019. Sensitivity of a leaf gas-exchange model for estimating paleoatmospheric CO₂ concentration. *Climate of the Past* **15**: 795–809.

Sancho-Knapik D, Escudero A, Mediavilla S, Scoffoni C, Zailaa J, Cavender-Bares J, Alvarez-Arenas TJ, Molins A, Alonso-Forn D, Ferrio JP et al. 2021. Deciduous and evergreen oaks show contrasting adaptive responses in leaf mass per area across environments. *New Phytologist* **230**: 521–534.

Saviряна Y, Klingenberg CP. 2011. Beyond bilateral symmetry: geometric morphometric methods for any type of symmetry. *BMC Evolutionary Biology* **11**: 280.

Schluter D, McPhail JD. 1992. Ecological character displacement and speciation in sticklebacks. *American Naturalist* **140**: 85–108.

Skelton RP, Anderegg LDL, Diaz J, Klinga MM, Pappera P, Lamarqued LJ, Delzone S, Dawsona TE, Ackerlya DD. 2021. Evolutionary relationships between drought-related traits and climate shape large hydraulic safety margins in western North American oaks. *Proceedings of the National Academy of Science* **118**: e2008987118.

Stuart YE, Losos JB. 2013. Ecological character displacement: glass mostly full or empty? *Trends in Ecology & Evolution* **28**: 1–7.

Torres-Ruiz JM, Kremer A, Carins-Murphy MR, Brodribb TJ, Lamarque LJ, Truffaut L, Bonne F, Ducoussou A, Delzon S. 2019. Genetic differentiation in functional traits among European sessile oak populations. *Tree Physiology* **39**: 1736–1749.

Ueno S, Taguchi Y, Tsumura Y. 2008. Microsatellite markers derived from *Quercus mongolica* var. *crispula* (*Fagaceae*) inner bark expressed sequence tags. *Genes & Genetic Systems* **83**: 179–187.

Ueno S, Tsumura Y. 2008. Development of ten microsatellite markers for *Quercus mongolica* var. *crispula* by database mining. *Conservation Genetics* **9**: 1083–1085.

Venables WN, Ripley BD. 2002. *Modern applied statistics with S*, 4^th^ edn. New York, NY, USA: Springer.

Viscosi V, Cardini A. 2011. Leaf morphology, taxonomy and geometric morphometrics: a simplified protocol for beginners. *PLoS ONE* **6**: e25630.

Viscosi V, Fortini P. 2011. Leaf shape variation and differentiation in three sympatric white oak species revealed by elliptic Fourier analysis. *Nordic Journal of Botany* **29**: 632–640.

Viscosi V, Fortini P, Slice DE, Loy A, Blasi C. 2009. Geometric morphometric analyses of leaf variation in four oak species of the subgenus *Quercus* (Fagaceae). *Giolnale Botanico Italiano* **143**: 575–587.

Wickham H. 2016. GPLOT: Elegant graphics for data analysis. New York, NY, USA: Springer-Verlag.

Wright IJ, Reich PB, Westoby M, Ackerly DD, Baruch Z, Bongers F, Cavender-Bares J, Chapin T, Cornelissen JHC, Diemer M et al. 2004. The worldwide leaf economics spectrum. *Nature* **428**: 821–827.

Yukilevich R. 2012. Asymmetrical patterns of speciation uniquely support reinforcement in Drosophila. *Evolution* **66**: 1430–1446.

Zwieniecki MA, Boyce CK, Holbrook NM. 2004. Hydraulic limitations imposed by crown placement determine final size and shape of *Quercus rubra* L. leaves. *Plant, Cell & Environment* **27**: 357–365.
**Supporting Information**

Additional Supporting Information may be found online in the Supporting Information section at the end of the article.

**Fig. S1** Posterior distributions for mutation-scaled migration rate (M) for *Quercus aliena* → *Quercus dentata* and *Q. dentata* → *Q. aliena* in sympatry and allopatry with corresponding 95% credible intervals.

**Fig. S2** Population clusters identified with STRUCTURE software.

**Fig. S3** The mean landmark output of *Quercus aliena* and *Quercus dentata* in sympatry and allopatry by geometric morphometric methods.

**Fig. S4** Plot of leaf area and specific leaf area of *Quercus aliena* and *Quercus dentata* in sympatry vs allopatry.

**Fig. S5** Inter-specific divergence along the first principal component (PC1) axis at genetic markers and morphological markers in sympatry and allopatry.

**Fig. S6** Leaf shape variation, correlation coefficient (R) and significance (P) value for genetically pure individuals in all, sympatric and allopatric populations as a function of whole climatic range.

**Fig. S7** Specific leaf area variation, correlation coefficient (R) and significance (P) value for all, sympatric and allopatric populations as a function of climatic factors (precipitation seasonality, annual precipitation, temperature seasonality and annual mean temperature).

**Fig. S8** Line chart of flowering phenology observations of *Quercus aliena* and *Quercus dentata* in one sympatric site.

**Table S1** Details of the 47 natural populations of *Quercus aliena* and *Quercus dentata*.

**Table S2** Detailed information for the 12 selected pairs of nuclear simple sequence repeat (nrSSR) primers.

**Table S3** The population information for leaf area and specific leaf area.

**Table S4** Loading values of climate principal component analysis.

**Table S5** Estimation of genetic diversity and heterozygote deficit in sympatric and allopatric oak populations based on 12 nuclear simple sequence repeat (nrSSR) loci.

**Table S6** Analysis of molecular variance (AMOVA) for sympatric and allopatric populations based on nuclear simple sequence repeats (nrSSRs).

**Table S7** Means and standard deviations of mutation-scaled migration rate (M) for *Quercus aliena* → *Quercus dentata* and *Q. dentata* → *Q. aliena* in sympatry and allopatry with corresponding 95% credible intervals computed using the posterior distribution estimates obtained from 12 nuclear simple sequence repeat (nrSSR) loci.

**Table S8** Number of ‘pure’ or ‘hybrid’ individuals assigned using morphological and genetic data. The bold indicate that the morphological and genetic criteria matched.

**Table S9** The parameters and standard deviation of the Gaussian distribution curves for the first principal component (PC1) of *Quercus aliena* and *Quercus dentata* in sympatry and allopatry for genetics vs morphological markers using a threshold-free approach.

**Table S10** Variance summary of climate using principal component analysis.

**Table S11** Phenology of *Quercus aliena* and *Quercus dentata* in one sympatric site on Shangfang Mountain.

Please note: Wiley Blackwell are not responsible for the content or functionality of any Supporting Information supplied by the authors. Any queries (other than missing material) should be directed to the *New Phytologist* Central Office.