\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[T2A]{fontenc}
\usepackage[ngerman,greek,polish,english]{babel}



\usepackage{float}








\begin{document}

\title{Common-garden study of introgression at loci associated with traits
adaptive to coastal environment from Quercus dentata into Q. mongolica
var. crispula}



\author[1]{Teruyoshi Nagamitsu}%
\author[1]{Kentaro Uchiyama}%
\author[1]{Ayako Izuno}%
\author[2]{Hajime Shimizu}%
\affil[1]{Forest Research and Management Organisation}%
\affil[2]{Hokkaido Research Organisation Forest Research Department}%


\vspace{-1em}



  
  \date{June 09, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
Adaptive introgression has been proposed in oaks (genus Quercus). In
northern Japan, Q. mongolica var. crispula (Qc) is common in inland
habitats, and Q. dentata (Qd) occurs in coastal habitats. At the
northern distributional limit of Qd, Q. \selectlanguage{ngerman}× angustilepidota (Qa), a hybrid
taxon between Qc and Qd, occurs in coastal habitats. The three taxa, Qc,
Qa, and Qd, were transplanted to inland and coastal common gardens
(sites). Genotypes at 27495 loci, phenotypes of eight traits of leaves
and shoots, and 30-year-old tree size were measured for 224 individuals
in both sites. Genotypic variation revealed a hybrid zone between Qc and
Qd, including both northern-edge Qd admixed with Qc and coastal Qa
backcrossed to Qc. Phenotypes of Qa trees were intermediate between
those of Qc and Qd trees. Size of Qa and Qd trees was smaller than that
of Qc trees in the inland site but was larger in the coastal site,
suggesting adaptation of Qa and Qd to coastal environment. Local
ancestry was estimated from phased genotypes of admixed trees using
reference genotypes of 47 Qc and 25 Qd trees, indicating heterogeneous
ancestry along chromosomes. Association mapping of genotypes and
admixture mapping of ancestry suggested that some loci potentially
associated with four traits were related to stress response and were
located at introgressed genomic regions. Further studies are necessary
to show the genetic basis of adaptive introgression resulting in Qd-like
phenotypes of Qa in coastal habitats.%
\end{abstract}\selectlanguage{ngerman}%



\sloppy


\textbf{Title}

Common-garden study of introgression at loci associated with traits
adaptive to coastal environment from \emph{Quercus dentata} into\emph{Q.
mongolica} var. \emph{crispula}

\textbf{Authors}

Teruyoshi Nagamitsu,\textsuperscript{1} Kentaro
Uchiyama,\textsuperscript{1} Ayako Izumo,\textsuperscript{1} and Hajime
Shimizu\textsuperscript{2}

\textbf{Affiliations}

\textsuperscript{1} Department of Forest Molecular Genetics and
Biotechnology, Forestry and Forest Products Research Institute, Forest
Research and Management Organization, Tsukuba 305-8687, Japan

\textsuperscript{2} Donan Station, Forestry Research Department,
Hokkaido Research Organization, Hakodate 041-0801, Japan

\textbf{Author for Correspondence}

Teruyoshi Nagamitsu orcid.org/0000-0003-3366-7197

nagamit@affrc.go.jp

+81-29-829-8264

\textbf{Contacts of co-authors}

Kentaro Uchiyama orcid.org/0000-0003-4083-3034

kruchiyama@affrc.go.jp

Ayako Izumo orcid.org/0000-0002-9339-2845

izuno@ffpri.affrc.go.jp

Hajime Shimizu

shimizu-hajime@hro.or.jp

\textbf{Abstract}

Adaptive introgression has been proposed in oaks (genus \emph{Quercus}
). In northern Japan, \emph{Q. mongolica} var. \emph{crispula}
(\emph{Qc} ) is common in inland habitats, and \emph{Q. dentata}
(\emph{Qd} ) occurs in coastal habitats. At the northern distributional
limit of \emph{Qd} ,\emph{Q.} \selectlanguage{ngerman}× \emph{angustilepidota} (\emph{Qa} ), a
hybrid taxon between\emph{Qc} and \emph{Qd} , occurs in coastal
habitats. The three taxa,\emph{Qc} , \emph{Qa} , and \emph{Qd} , were
transplanted to inland and coastal common gardens (sites). Genotypes at
27495 loci, phenotypes of eight traits of leaves and shoots, and
30-year-old tree size were measured for 224 individuals in both sites.
Genotypic variation revealed a hybrid zone between \emph{Qc} and
\emph{Qd} , including both northern-edge \emph{Qd} admixed with
\emph{Qc} and coastal \emph{Qa}backcrossed to \emph{Qc} . Phenotypes of
\emph{Qa} trees were intermediate between those of \emph{Qc} and
\emph{Qd} trees. Size of\emph{Qa} and \emph{Qd} trees was smaller than
that of \emph{Qc} trees in the inland site but was larger in the coastal
site, suggesting adaptation of \emph{Qa} and \emph{Qd} to coastal
environment. Local ancestry was estimated from phased genotypes of
admixed trees using reference genotypes of 47 \emph{Qc} and 25 \emph{Qd}
trees, indicating heterogeneous ancestry along chromosomes. Association
mapping of genotypes and admixture mapping of ancestry suggested that
some loci potentially associated with four traits were related to stress
response and were located at introgressed genomic regions. Further
studies are necessary to show the genetic basis of adaptive
introgression resulting in \emph{Qd} -like phenotypes of \emph{Qa} in
coastal habitats.

\textbf{Keywords}

adaptive introgression, genome-wide association study, hybridization,
local ancestry, phenotypic plasticity

\textbf{1 \textbar{} Introduction}

Introgression is the transfer of alleles from one taxon (donor) into
genomes of another taxon (recipient) through hybridization,
recombination, and backcrossing (Anderson 1953; Barton and Hewitt 1985).
Generally, alleles introgressed from the donor taxon can be deleterious,
because of conflicting interactions in the genomes and unfavorable
phenotypes in the environment of the recipient taxon (Burke and Arnold
2001; Goulet et al. 2017). Thus, introgressed alleles will be purged
from the recipient genomes, and the genetic integrity of each taxon can
be maintained (Rieseberg and Carney 1998). However, in specific
environment, where the donor-like phenotypes are more adaptive than the
recipient phenotypes, introgressed alleles can be selected and fixed in
the recipient genomes (Anderson and Stebbins 1954; Suarez-Gonzalez et
al. 2018c). Such environment-dependent adaptive introgression has been
found in some plants, resulting in the expansion to novel habitats and
the creation of different ecotypes in the recipient taxa (Martin et al.
2006; Rieseberg et al. 2007; Whitney et al. 2010; Arnold et al. 2016).
Adaptive introgression tends to occur at distributional range margins of
either taxon, where hybridization with the other taxon is frequent due
to reduced opportunity of mating within marginal populations, and
introgressed populations of the recipient taxa can colonize habitats
suitable for the donor taxa (Chhatre et al. 2018; Ma et al. 2019; Menon
et al. 2021; Rendón-Anaya et al. 2021). This process is thought to
contribute colonization and range shift of forest trees during past
climate change and offer an option of forest management under future
climate change (Petit et al. 2004; Hamilton and Miller 2016).

Oaks (genus \emph{Quercus} , Fagaceae) are temperate forest trees in the
northern hemisphere and are often interfertile among species (Denk et
al. 2017). Introgression between oak taxa has long been proposed
(Stebbins et al. 1947; Sork et al. 2016), and their genomes are regarded
as phylogenetic mosaics of different origins (Kim et al. 2018; Hipp et
al. 2020). Oaks inhabiting heterogeneous environments usually have
different phenotypic and genetic variations (Ortego et al. 2014; Riordan
et al. 2016; Cavender-Bares 2019) and often exhibit local adaptation to
indigenous environments (Sork 2018; Leroy et al. 2020). In northern
Japan, there are two species of white oaks (section \emph{Quercus}
),\emph{Q. mongolica} Fischer ex Ledebour var. \emph{crispula} (Blume)
H. Ohashi (\emph{Qc} ) is common in inland habitats, and \emph{Q.
dentata}Thunberg ex Murray (\emph{Qd} ) occurs in coastal habitats
(Matsumoto et al. 2009). In northern Hokkaido, \emph{Qd} trees are rare
because this area is the northern distributional limit of \emph{Qd} ,
whereas\emph{Qc} trees are abundant because \emph{Qc} is distributed to
more northern area, Sakhalin. In the northernmost area of Hokkaido, a
coastal\emph{Qc} ecotype with unique traits, which are similar to
\emph{Qd}phenotypes and tolerant to coastal stress, occurs in coastal
forests (Nagamitsu et al. 2019). Some taxonomists regarded this
coastal\emph{Qc} ecotype as a hybrid taxon \emph{Q.} ×
\emph{angustilepidota}Nakai (\emph{Qa} ) between \emph{Qc} and \emph{Qd}
(Ohba 2006; Aizawa et al. 2021). Multi-locus nuclear microsatellite
genotypes supported this hybrid origin (Nagamitsu et al. 2019), and
genome-wide single nucleotide polymorphism (SNP) genotypes demonstrated
environment-dependent introgression from \emph{Qd} to \emph{Qc} ,
resulting in \emph{Qa} that included hybrids after the first generation
of backcross to \emph{Qc}(Nagamitsu et al. 2020). Thus, introgression of
\emph{Qd} alleles to\emph{Qc} genomes at loci associated with traits
adaptive to coastal environment is expected but has not been confirmed
yet.

Oak trees in coastal habitats suffer from various stress, such as wind,
salinity, drought, heat, substrate instability, and nutrient scarcity
(Hesp 1991). For example, strong wind and salt spray in winter cause the
mortality of buds in the upper parts of shoots, resulting in frequent
dieback of shoots, slow elongation of stems, and low canopy height of
coastal oak forests (Asai et al. 1986). These stresses are likely to
cause natural selection in leaf and shoot traits, leading to local
adaptation (Ramírez-Valiente et al. 2010; Ciccarelli and Bona 2022).
Common-garden experiment using inland and coastal test sites is
effective to show local adaptation in \emph{Qc} and \emph{Qd} .
Relatively higher performance of \emph{Qa} trees with \emph{Qd} -like
phenotypes of leaf and shoot traits in coastal sites than in inland
sites will indicate the adaptation of these traits to coastal habitats.
Difference in phenotypes among admixed trees between inland and coastal
habitats depends not only on genetic and environmental variations but
also on their interaction, namely phenotypic plasticity
(Ramírez-Valiente et al. 2010). Common-garden experiment is useful to
elucidate phenotypic plasticity and to discriminate genetic variation
from environmental variation in phenotypes.

Phenotypes of admixed recombinants in common gardens enable us to detect
loci associated with focal traits (Rieseberg and Buerkle 2002). This
approach is referred to as admixture mapping and is useful in long-lived
and large-sized organisms, for which pedigrees from artificial crosses
are difficult to obtain (Buerkle and Lexer 2008). In poplar, one of the
forest trees with rich genomic information, admixture mapping was
applied to various traits (Suarez-Gonzalez et al. 2018a; Bresadola et
al. 2019). Admixture mapping requires the inference of locus-specific
(local) ancestry, the number of alleles that originate from either
ancestral population involved in the admixture at individual loci
(Lindtke et al. 2013). Local ancestry usually deviates from genome-wide
ancestry of individuals, indicating the excess or deficit of admixture
at individual loci (Buerkle and Lexer 2008). A genome-wide pattern of
introgression is often illustrated using various approaches (Martin and
Van Belleghem 2017), for example, the Paterson's \emph{D} statistics
from the ABBA-BABA test that discriminates introgression from incomplete
lineage sorting (Martin et al. 2015). These approaches require a
sufficient number of genome-wide loci, which can be obtained from
sequences of reduced genomic libraries mapped to a whole-genome
reference sequence that has been available in white oaks (Plomion et al.
2018; Sork et al. 2022). Recently, reference sequences of \emph{Q.
mongolica} var. \emph{mongolica} and \emph{Qd} also has been published
(Ai et al. 2022; Wang et al. 2023). Using genotypes, phenotypes, and
performance of trees admixed between \emph{Qc} and \emph{Qd} in common
gardens, we can verify adaptive introgression of \emph{Qd} alleles at
loci associated with adaptive traits into \emph{Qc} genome.

In this study, we tried to obtain evidence for the adaptive
introgression using inland and coastal common gardens (sites), where
seedlings of the focal white oak taxa, \emph{Qc} , \emph{Qa} ,
and\emph{Qd} , were planted from various provenances in Hokkaido and had
now reached about 30 years old. First, we examined their genetic
variation using genome-wide SNP genotypes to estimate genomic
compositions of hybrids and genome-wide patterns of introgression. Next,
we measured phenotypes of leaf and shoot traits and performance of trees
in each site to demonstrate adaptation of \emph{Qd} -like phenotypes of
these traits to coastal environment. Finally, we searched for loci, at
which SNP genotypes and local ancestry were associated with those
traits. We expected that these trait-associated loci were located at
introgressed genomic regions and were close to genes involved in
adaptation to coastal environment.

\textbf{2 \textbar{} Materials and methods}

\textbf{2.1 \textbar{} Common gardens}

Two common gardens were established in an inland site (44.7\selectlanguage{english}@N, 142.1@\selectlanguage{english}E,
20 m elevation) and a coastal site (45.3@\selectlanguage{english}N, 141.6@\selectlanguage{english}E, 5 m elevation) in
northern Hokkaido (Figure S1 in Supporting Information). In each site,
the three taxa, \emph{Quercus mongolica} var. \emph{crispula}(\emph{Qc}
), \emph{Q. dentata} (\emph{Qd} ), and their hybrids, \emph{Q.}\selectlanguage{ngerman}×
\emph{angustilepidota} (\emph{Qa} ), were transplanted from natural
forests in 32 provenances in Hokkaido (Figure 1, Table S1). Acorns
(seeds) of several maternal trees were collected from a forest stand in
each provenance, and their seedlings were reared in a nursery in the
Dohoku Station of the Forestry Research Institute of the Hokkaido
Research Organization. In the inland site located in the Dohoku Station,
seedlings were planted at 1.4 m intervals along plantation rows spaced
2.8 m apart in three blocks in 1992--1993, and half of trees in two
blocks and a few trees in one block were thinned in 2006, when tree
crowns became crowded (Figure S2). In the coastal site on sand dunes at
500 m distance from the seashore, seedlings were planted at 0.8 m
intervals along plantation rows spaced 2.8 m apart in two blocks in
1986, and trees were not thinned because tree crowns had been sparse
(Figure S3).

We observed living trees in both sites and identified the taxa of trees
in 2017--2019 (Figures S2, S3). Because taxonomic identification
of\emph{Qc} and \emph{Qd} is ambiguous due to intermediate phenotypes
(Ohba 2006; Aizawa et al. 2021), we classified the taxa using the
following criteria based on hairs on shoots according to our previous
studies (Nagamitsu et al. 2019; Nagamitsu et al. 2020). Trees with
pubescent one-year-old shoots with trichomes (stellate hairs) were
identified as \emph{Qd} . Among trees with glabrous one-year-old shoots,
trees from provenances in coastal habitats, sand dunes and coastal
cliffs, were identified as \emph{Qa} , which was recognized as the
coastal ecotype of \emph{Qc} in the previous studies (Nagamitsu et al.
2019; Nagamitsu et al. 2020), and trees from provenances in inland
habitats, hills and mountains, were identified as \emph{Qc} (Table S1).
In addition to living trees, we observed dead trees that had died
recently. We were not able to record individuals that had died when they
were saplings. We confirmed the provenances of living and dead trees
based on plantation records of the common gardens.

The following measurement and collection were conducted in both sites in
2017--2019. We measured the girth (cm) of every stem of each living tree
to evaluate tree size. The girth at the breast height in the inland site
and that at the ground level in the coastal site were measured, because
trees in the coastal site were crooked and branching near the ground due
to strong wind. The locations of living and dead trees and the girths of
living trees of \emph{Qc} , \emph{Qa} , and \emph{Qd} were recorded in
both sites (Figures S2, S3). We collected some branches of trees sampled
from living trees in July to obtain DNA and to measure phenotypes. Some
flesh leaves on current-year shoots were stored at --20\selectlanguage{english}@C, and DNA was
extracted from leaf tissues using DNeasy Plant Mini Kit (Qiagen, Hilden,
Germany). Three leaves on current-year shoots and three one-year-old
shoots selected from collected branches of each sampled tree were dried
and preserved for morphological measurement.

\textbf{2.2 \textbar{} Genotyping}

A double-digest restriction-site-associated DNA (ddRAD) library was
prepared using a modified protocol of the standard method (Peterson et
al. 2012) described in our previous study (Nagamitsu et al. 2020).
Extracted DNA was digested with \emph{Pst} I and \emph{Sau} 3AI
restriction enzymes, ligated with Y-shaped adaptors, and amplified by
polymerase chain reaction (PCR) with the KAPA HiFi polymerase (Kapa
Biosystems, Woburn, MA, USA). After PCR amplification with
adapter-specific primer pairs (Access Array Barcode Library for
Illumina, Fluidigm, South San Francisco, CA, USA), an equal amount of
DNA from each tree was mixed and size-selected with the BluePippin
agarose gel (Sage Science, Beverly, MA, USA). Approximately 450-bp
library fragments were retrieved. The library quality was checked using
a 2100 Bioanalyzer with a high sensitivity DNA chip (Agilent
Technologies, Waldbronn, Germany). The ddRAD library was sequenced using
Illumina MiSeq and HiSeq to generate paired-end reads with a 150-bp
length.

Obtained reads were mapped to the reference sequences consisting of 12
pseudomolecules (chromosomes) of the \emph{Q. robur} genome assembly
PM1N (Plomion et al. 2018). Although some reference sequences of white
oak species including \emph{Q. mongolica} var. \emph{mongolica}
and\emph{Qd} are available (Sork et al. 2022; Ai et al. 2022; Wang et
al. 2023), we selected PM1N because \emph{Q. robur} is a species
closely-related to the common ancestor of \emph{Qc} and \emph{Qd} . The
read mapping and variant calling were conducted using dDocent 2.7.8
(Puritz et al. 2014). The subsequent filtering procedures in detail are
explained in our previous study (Nagamitsu et al. 2020). From the
variant loci with various types of polymorphism, we selected loci that
were biallelic without indels and polymorphic with high sequencing
quality, \textless{} 5\% missing genotypes, and \textgreater{} 5\% minor
allele frequency, using VCFtools 0.1.14 (Danecek et al. 2011). We
removed loci mapped to the regions of transposable elements determined
in \emph{Q. robur} (Plomion et al. 2018). In addition, we removed loci
that extremely deviated from the Hardy--Weinberg equilibrium
(\emph{P}\textless{} 0.01) probably due to null alleles within a
population of each taxon, \emph{Qc} , \emph{Qa} , or \emph{Qd} , using
VCFtools.

\textbf{2.3 \textbar{} Genomic composition}

To summarize genetic structure of sampled trees in the inland and
coastal sites and to obtain allele frequencies in ancestral populations
that represent \emph{Qc} and \emph{Qd} , we conducted Bayesian
clustering for SNP genotypes at ddRAD loci using Admixture 1.3
(Alexander and Lange 2011). The log-likelihood was obtained, when the
number (\emph{K} ) of ancestral populations (clusters) was 1, 2, or 3.
Allele frequencies at ddRAD loci in each cluster (\emph{P} ) and
ancestry proportions of each cluster in individual trees (\emph{Q} )
were estimated for each \emph{K} . Genetic differentiation between
clusters (pairwise\emph{F} \textsubscript{ST}) was also estimated. We
regarded two clusters (\emph{K} = 2) as ancestral populations that
represented \emph{Qc} and\emph{Qd} and estimated allele frequencies
\emph{P} for the two clusters.

To describe genomic compositions of sampled trees in both sites, we
estimated the proportion of alleles derived from \emph{Qd} ancestral
population (\emph{S} ) and the proportion of heterozygous loci with both
alleles derived from \emph{Qc} and \emph{Qd} ancestral populations
(\emph{H} ) using the package HIest 2.0 (Fitzpatrick 2012) in R 3.3.2
(R\_Core\_Team 2019). We regarded allele frequencies \emph{P} for the
two clusters obtained from the Bayesian clustering as\emph{P}
\textsubscript{1} and \emph{P} \textsubscript{2} allele frequencies of
ancestral populations of \emph{Qc} and \emph{Qd} , respectively, and
estimated \emph{S} (\emph{Qd} ancestry) and \emph{H}(inter-ancestry
heterozygosity) of individual trees using the function HIest with the
options, method = SANN, iterations = 10000.

\textbf{2.4 \textbar{} Phenotyping}

We measured eight morphological traits of leaves and shoots, which were
thought to be distinctive between \emph{Qc} and \emph{Qd} and associated
with the tolerance to coastal stress (Nagamitsu et al. 2019).
Measurements for three leaves and three shoots of branches collected
from each tree were averaged and regarded as a phenotypic value of each
tree.

To describe leaf shape, three traits were recorded: relative leaf width,
lateral vein interval, and tooth angle (Figure S4). Among the three
traits, the two latter traits are important to identify \emph{Qc}
and\emph{Qd} taxa (Ohba 2006; Aizawa et al. 2021). The length between
the apex and base of each leaf and the width at the widest part of the
leaf (mm) were measured, and the number of lateral veins were counted.
Relative leaf width was calculated from (leaf width)/(leaf length), and
lateral vein interval (mm) was calculated from (leaf length)/(number of
lateral veins + 1). Tooth angle (@\selectlanguage{english}) was measured on a serration at the
central part of each leaf.

The density and size of stellate hairs on leaves vary between
\emph{Qc}and \emph{Qd} (Ishida et al. 2003) and may have ecological and
physiological functions (Bickford 2016). Stellate hairs on the lower
surface of each leaf were observed in a circular leaf area with a 3.1 mm
diameter between lateral veins using a stereo microscope (Figure S4).
The number of stellate hairs in each area were counted, and the density
(mm\textsuperscript{--2}) was obtained. As the size of stellate hairs,
the length of a radial filament (mm) of a typical stellate hair was
measured (Figure S4).

Leaf mass per area (LMA) and shoot diameter may be related with
adaptation to intense light, severe drought, and strong wind in coastal
habitats (Sancho-Knapik et al. 2021). A leaf disc with a 9 mm diameter
was collected from the central part of each leaf (Figure S4), and the
dry weight (mg) of each disc was measured after dehydration at 60@\selectlanguage{english}C for
48 h to calculate LMA (mg mm\textsuperscript{--2}). Shoot diameter (mm)
was measured at the upper part of shoots. Salt spray and harsh wind in
winter cause the mortality of buds, particularly those in the upper part
of shoots (Asai et al. 1986). Bud production in the lower part of shoots
is likely to convey the tolerance to coastal stress, because this bud
production can compensate for high bud mortality in the upper part of
shoots. Thus, the number of axillary buds at bud-scale scars located in
the lowest part of shoots was recorded (Figure S4).

\textbf{2.5 \textbar{} Phenotypic variation}

To summarize phenotypic variation of sampled trees in the inland and
coastal sites, we performed principal component analyses (PCA) using the
function prcomp in R 3.3.2. Among the eight traits, the frequency
distributions of the density and size of stellate hairs and the scale
bud number, were skewed. Thus, values of the size, the density + 1, and
the number + 1 were log-transformed. The phenotypic values of the eight
traits were scaled (mean = 0 and variance = 1) and applied to PCA.
Contributions of each principal component (PC) to the total phenotypic
variation and loadings of each trait to the first and second PCs were
obtained. Statistical differences in PC values among taxa were verified
using the function kruskal.test in R 3.3.2.

To evaluate phenotypic plasticity in response to inland and coastal
environments, we plotted the \emph{Qd} ancestry and the first phenotypic
PC values of individual trees and fitted smoothed lines in each of the
inland and coastal sites using the function lowess with the
parameter\emph{f} = 1 in R 3.3.2. We applied a linear model to the data
using the function lm in R 3.3.2,

\emph{y} \textasciitilde{} \selectlanguage{greek}α \selectlanguage{english}+ \selectlanguage{greek}β\selectlanguage{english}\textsubscript{1}\emph{x}
\textsubscript{1} + \selectlanguage{greek}β\selectlanguage{english}\textsubscript{2}\emph{x} \textsubscript{2} +
\selectlanguage{greek}γ\selectlanguage{english}\emph{x} \textsubscript{1}\emph{x} \textsubscript{2},

where \emph{y} is the first phenotypic PC value,\emph{x}
\textsubscript{1} is the \emph{Qd} ancestry, and\emph{x}
\textsubscript{2} is 0 in the inland site and 1 in the coastal site; \selectlanguage{greek}α
\selectlanguage{english}is an intercept, \selectlanguage{greek}β\selectlanguage{english}\textsubscript{1} and \selectlanguage{greek}β\selectlanguage{english}\textsubscript{2} are
coefficients of \emph{x} \textsubscript{1} and \emph{x}
\textsubscript{2}effects, respectively, and \selectlanguage{greek}γ \selectlanguage{english}is a coefficient of their
interactions. Significant \selectlanguage{greek}γ \selectlanguage{english}indicates different slopes between the
sites, suggesting different reaction norms of phenotypic plasticity in
response to coastal environment depending on the \emph{Qd} ancestry.

\textbf{2.6 \textbar{} Tree performance}

To show local adaptation of trees with different phenotypes to inland
and coastal environments, we plotted the mean values of the first
phenotypic PC and the mean values of the stem basal area of all living
and dead trees in individual plantation rows and fitted smoothed lines
in each of the inland and coastal sites using the function lowess
with\emph{f} = 1 in R 3.3.2. We interpreted different slopes of fitted
lines between the sites as local adaptation. The stem basal area
(cm\textsuperscript{2}) of living trees was calculated from the measured
girths of their stems. The stem basal area of dead trees was treated as
0 cm\textsuperscript{2}. Thinned trees in the inland site were excluded
from the calculation. Statistical differences in the mean stem basal
area among taxa in each site were verified using the function
kruskal.test in R 3.3.2.

\textbf{2.7 \textbar{} Local ancestry and genome-wide patterns of
introgression}

Haplotypes were estimated from SNP genotypes at ddRAD loci in each
chromosome through imputation and phasing using Beagle 5.2 (Browning et
al. 2018; Browning et al. 2021). Sampled trees in both sites were
divided into \emph{Qc} reference individuals (the \emph{Qd}
ancestry\emph{S} = 0), admixed individuals (0 \textless{} \emph{S}
\textless{} 1), and \emph{Qd} reference individuals (\emph{S} = 1). We
converted a haplotype matrix in each of three vcf files of the \emph{Qc}
reference individuals, the \emph{Qd} reference individuals, and the
admixed individuals into a Numpy array using the function vcf2npy. Local
ancestry at ddRAD loci of the admixed individuals was described as
homozygote of alleles derived from either \emph{Qc} or
\emph{Qd}ancestral population and heterozygote of alleles derived from
both\emph{Qc} and \emph{Qd} ancestral populations. The three ancestry
types were inferred from the haplotype arrays using the Python package
Loter with the functions lc.loter\_smooth and
loter.locanc.local\_ancestry (Dias-Alves et al. 2018). The package
assumed that haplotypes of admixed individuals originated from
hybridization and recombination of ancestral populations of reference
individuals (Dias-Alves et al. 2018). To visualize genome-wide patterns
of introgression, we calculated the mean number of alleles derived from
\emph{Qd} at each locus in the admixed individuals.

To depict genome-wide introgression patterns in another method, we
estimated the Patterson's \emph{D} statistic at ddRAD loci (Martin et
al. 2015). The \emph{D} statistic requires allele frequencies of four
populations with phylogenetic relationship (((P1, P2), P3), O), and
positive \emph{D} values indicate introgression from P3 to P2 as
discriminating incomplete lineage sorting. We assigned the
\emph{Qc}reference individuals to P1, the admixed individuals to P2,
the\emph{Qd} reference individuals to P3, and the \emph{Q. robur}
reference sequence to O. We calculated the number of derived alleles
that are different from the \emph{Q. robur} reference (ancestral)
alleles for individuals of P1, P3, and P3 from the vcf files and
obtained the frequency of derived alleles, \emph{p}
\textsubscript{1},\emph{p} \textsubscript{2}, and \emph{p}
\textsubscript{3}, respectively, using the package gaston 1.5.9 in R
3.3.2. Using the ABBA-BABA statistics, \emph{C} \textsubscript{ABBA} =
(1 --\emph{p} \textsubscript{1}) \emph{p} \textsubscript{2}\emph{p}
\textsubscript{3} and \emph{C} \textsubscript{BABA} =\emph{p}
\textsubscript{1} (1 -- \emph{p} \textsubscript{2})\emph{p}
\textsubscript{3}, at each locus, we obtained \emph{D} = (\selectlanguage{greek}Σ\selectlanguage{english}\emph{C}
\textsubscript{ABBA} -- \selectlanguage{greek}Σ\selectlanguage{english}\emph{C} \textsubscript{ABBA}) / (\selectlanguage{greek}Σ\selectlanguage{english}\emph{C}
\textsubscript{ABBA} + \selectlanguage{greek}Σ\selectlanguage{english}\emph{C} \textsubscript{ABBA}) at sliding
windows of 101 neighboring loci, which sufficiently reduced errors in
the calculation of \emph{D} .

\textbf{2.8 \textbar{} Association of loci with traits}

To detect loci associated with each of the eight traits, we examined
genome-wide association (GWA) of SNP genotypes (association mapping) and
GWA of local ancestry (admixture mapping) at ddRAD loci with the
phenotypic values using mixed linear models. Despite the similarities
between association mapping and admixture mapping, some differences
between them in terms of aims and methods result in both merits and
demerits of each (Buerkle and Lexer 2008). Association mapping can deal
with complicated background of ancestry and variation within ancestral
lineages underlying phenotypes using more loci densely distributed
across a genome, while admixture mapping can work with fewer loci
sparsely distributed owing to larger blocks of linkage disequilibrium in
recombinants of distinct ancestries (Buerkle and Lexer 2008).

Phenotypic values of individual trees with the same \emph{Qd} ancestry
were expected to differ between the inland and coastal sites due to
phenotypic plasticity. To exclude such confounding effects with
environment and plasticity and to conduct GWA in both sites together, we
scaled phenotypic values (mean = 0 and variance = 1) in each site
separately and applied a mixed linear model to pooled data of the scaled
values in both sites. Environmental heterogeneity within the sites were
incorporated as random errors in the models.

To conduct GWA, we used single-locus GWA (Yu et al. 2006) implemented in
the package rrBLUP 4.5 (Endelman 2011) in R 3.3.2. To discriminate
confounding factors with genetic structure of sampled trees in both
sites, we used the Q + K model with a population structure matrix Q
calculated from probabilistic principal component analysis using the
function ppca in the package pcaMethods (Stacklies et al. 2007) and a
kinship matrix K calculated from realized additive relationship using
the function A.mat in rrBLUP. The first four principal components were
used for the Q matrix. All parameters in the GWA were set at the default
values. Significant associations were regarded as \textless{} 0.05 false
discovery rate (FDR) of --log\textsubscript{10}\emph{p} calculated using
the function p.adjust in R 3.3.2. Q-Q plots of expected and observed
--log\textsubscript{10} \emph{p} values were examined for associations
of SNP genotypes and local ancestry.

To search for genes involved in the adaptation to coastal environment
around the trait-associated loci, we selected genes (protein ID) nearest
to the loci detected by association mapping and within regions including
the loci detected by admixture mapping from PM1N v2.3 \emph{Q.
robur}annotation database
(https://urgi.versailles.inra.fr/OakMine\_PM1N) (Plomion et al. 2018).
We obtained gene ontology (GO) terms and protein descriptions of the
genes.

\textbf{3 \textbar{} Results}

\textbf{3.1 \textbar{} Genotypes}

In the inland site, 575 living and 197 dead trees were found, and 106
trees were sampled from the living trees (Figure S2, Data S1 in
Supporting Information). In the coastal site, 350 living and 46 dead
trees were found, and 118 trees were sampled from the living trees
(Figure S3, Data S1). Taxonomic assemblages of sampled trees were biased
to \emph{Qc} (41 \emph{Qc} , 40 \emph{Qa} , and 25 \emph{Qd} trees) in
the inland site and biased to \emph{Qd} (12 \emph{Qc} , 48 \emph{Qa} ,
and 58\emph{Qd} trees) in the coastal site probably due to different
mortality in inland and coastal habitats (Table S1).

SNP genotypes at 27495 ddRAD loci were obtained from the 224 sampled
trees in both sites. Bayesian clustering for SNP genotypes demonstrated
a larger increase in log-likelihood from one to two clusters (\emph{K} =
1--2) and an additional smaller increase from two to three clusters
(\emph{K} = 2--3; Figure S5). When two clusters representing
\emph{Qc}and \emph{Qd} were recognized, genetic differentiation based
on\emph{F} \textsubscript{ST} between these clusters was 0.183. Ancestry
proportion of the \emph{Qd} cluster (\emph{Q} ) was lower in
\emph{Qc}trees (0.00 \textless{} \emph{Q} \textless{} 0.17),
intermediate in\emph{Qa} trees (0.00 \textless{} \emph{Q} \textless{}
0.45), and intermediate or higher in \emph{Qd} trees (0.23 \textless{}
\emph{Q}\textless{} 1.00; Figure S5). These findings suggested that
genetic variation mainly existed between \emph{Qc} and \emph{Qd} and
that their admixture characterized \emph{Qa} and a part of \emph{Qd} .

For the sampled trees in both sites, the \emph{Qd} ancestry (\emph{S} )
and inter-ancestry heterozygosity (\emph{H} ) were estimated from allele
frequencies of the \emph{Qc} and \emph{Qd} clusters obtained from the
Bayesian clustering. Genomic compositions shown in the coordinates
of\emph{S} and \emph{H} indicated nearly random admixture between the
ancestral populations of \emph{Qc} and \emph{Qd} (Figure 2). In both
sites, \emph{Qd} trees from provenances near the distributional limit in
the northernmost area of Hokkaido showed lower \emph{Qd} ancestry (0.23
\textless{} \emph{S} \textless{} 0.91) than those from provenances in
the other southern areas (0.78 \textless{} \emph{S} \textless{} 1.00),
and \emph{Qa} trees from coastal provenances north to the distributional
limit of \emph{Qd} showed higher \emph{Qd} ancestry (0.00
\textless{}\emph{S} \textless{} 0.45) than \emph{Qc} trees from inland
provenances in Hokkaido (0.00 \textless{} \emph{S} \textless{} 0.06;
Figure 2). These findings suggested a hybrid zone between \emph{Qc} and
\emph{Qd}around the northern distributional limit of \emph{Qd} ,
including both northern-edge \emph{Qd} populations admixed with
\emph{Qc} and coastal\emph{Qa} populations backcrossed to \emph{Qc} .

\textbf{3.2 \textbar{} Phenotypes}

Phenotypic values of the eight morphological traits of leaves and shoots
were obtained from the 224 sampled trees in the inland and coastal sites
(Data S1). PCA for the eight traits summarized the first and second PCs
that contributed to 45.2\% and 16.7\%, respectively, of phenotypic
variation (Figure 3a). All the eight traits positively affected the
first PC, which represented \emph{Qd} -like phenotypes (Figure 3a).
Among the eight traits, four leaf traits (relative leaf width, lateral
vein interval, tooth angle, and stellate hair length), namely
taxon-specific traits, negatively affected to the second PC. Another two
leaf and two shoot traits (stellate hair density, LMA, shoot diameter,
and scale bud number), namely habitat-specific traits, positively
affected the second PC.

In both sites, \emph{Qc} , \emph{Qa} , and \emph{Qd} trees showed lower,
intermediate, and higher first PC values, respectively
(\emph{P}\textless{} 0.001), and \emph{Qa} and \emph{Qd} trees showed
higher and lower second PC values, respectively (\emph{P} \textless{}
0.001; Figure 3b, c). Thus, coastal \emph{Qa} trees mainly differed from
inland\emph{Qc} trees in the habitat-specific traits, while \emph{Qa}
trees mainly differed from \emph{Qd} trees in the taxon-specific traits.
Overall, both taxon-specific and habitat-specific traits contributed to
phenotypic variation between \emph{Qc} and \emph{Qd} , which was
represented by the first PC.

Relationship between the \emph{Qd} ancestry and the first phenotypic PC
value differed between the sites, indicating phenotypic plasticity
(Figure 4). The first phenotypic PC value was increasing with
the\emph{Qd} ancestry (\selectlanguage{greek}β\selectlanguage{english}\textsubscript{1} = 4.35, \emph{P} \textless{}
0.001) and higher in the coastal sites than in the inland site
(\selectlanguage{greek}β\selectlanguage{english}\textsubscript{2} = 2.20, \emph{P} \textless{} 0.001). A slope of the
first phenotypic PC value to the \emph{Qd} ancestry was more gentler in
the coastal site than in the inland site (\selectlanguage{greek}γ \selectlanguage{english}= --1.97,
\emph{P}\textless{} 0.001), indicating larger phenotypic response of
\emph{Qa}trees than in \emph{Qd} trees to coastal environment.

\textbf{3.3 \textbar{} Tree performance}

Relationship between the mean first phenotypic PC values and the stem
basal area of both living and dead trees in individual plantation rows
showed the opposite patterns between the inland and coastal sites
(Figure 5), suggesting local adaptation. While the tree size was
generally larger in the inland site than in the coastal site, the size
of \emph{Qc} trees was extremely reduced in the coastal site, but the
size of \emph{Qd} trees was similar in both sites (Figure 5). The size
of \emph{Qa} trees was smaller than that of \emph{Qc} trees (\emph{P} =
0.004) but was not significantly different from that of \emph{Qd} trees
(\emph{P} = 0.096) in the inland site, whereas the size of
\emph{Qa}trees was larger than that of \emph{Qc} trees (\emph{P} =
0.004) but was not significantly different from that of \emph{Qd} trees
(\emph{P}= 0.064) in the coastal site (Data S1). These results
suggested that\emph{Qd} -like phenotypes of \emph{Qa} trees in both
taxon-specific and habitat-specific traits were adaptive to coastal
environment.

\textbf{3.4 \textbar{} Local ancestry and genome-wide patterns of
introgression}

Haplotypes at 27495 ddRAD loci in each chromosome were obtained from
imputation and phasing for SNP genotypes of the 224 sampled trees in
both sites. Among the sampled trees, 47 \emph{Qc} trees (the
\emph{Qd}ancestry \emph{S} = 0) and 25 \emph{Qd} trees (\emph{S} = 1)
were selected as \emph{Qc} and \emph{Qd} reference individuals,
respectively (Table S1). Local ancestry (\emph{Qc} homzygote, \emph{Qc}
and \emph{Qd}heterozygote, and \emph{Qd} homozygote) at ddRAD loci in
each chromosome was inferred for the remaining 152 admixed individuals
based on the haplotypes of the \emph{Qc} and \emph{Qd} reference
individuals (Figure 6a). Genome-wide patterns of local ancestry in
admixed individuals indicated heterogeneous ancestry along chromosomes
(Figure 6b).

The Patterson's \emph{D} statistics obtained at sliding windows of 101
neighboring ddRAD loci were also heterogeneous along chromosomes (Figure
S7). Genomic regions showing excessive introgression (higher
\emph{D}values; Figures S7) were not always consistent with those with
higher local ancestry from \emph{Qd} (Figure 6b).

\textbf{3.5 \textbar{} Trait-associated loci}

Difference in phenotypic values between the inland and coastal sites and
dependence of phenotypic values on the \emph{Qd} ancestry varied among
the eight traits (Figures S8, S9). Phenotypic values of each trait were
scaled in each site separately. Scaled phenotypic values of the 224
sampled trees in both sites were pooled and applied to mixed linear
models with effects of genetic structure and kinship relationship. Based
on PCA of SNP genotypes at ddRAD loci, genetic variation
between\emph{Qc} and \emph{Qd} mainly determined genetic structure in
the first PC (Figure S10), and some northern-edge \emph{Qd} populations
had unique genetic properties in the second, third, and fourth PCs
(Figure S10b--d). Q-Q plots of association mapping of SNP genotypes and
admixture mapping of local ancestry showed consistent relationships
between expected and observed log-likelihood values at loci with a lower
range of --log\textsubscript{10} \emph{p} values, which represented the
majority of unassociated loci, indicating successful fitting (Figures
S11, S12).

We did not detect any loci significantly associated with the eight
traits (FDR \textgreater{} 0.05) but found some loci potentially
associated with one taxon-specific trait and three habitat-specific
traits (0.051 [?] FDR [?] 0.148; Figure 7, Table S2). Association mapping
detected three loci in chromosomes 2, 7, and 9 associated with lateral
vein interval (Figure 7c), and admixture mapping also supported one
locus in chromosome 2 of the three loci (Figure 7d). For loci associated
with stellate hair density, association mapping detected one locus in
chromosome 1 (Figure 7i), and admixture mapping detected loci in two
genomic regions in chromosomes 10 and 11 (Figure 7j). Admixture mapping
detected loci associated with LMA in chromosome 11 (Figure 7l), and
association mapping detected one locus associated with shoot diameter in
chromosome 1 (Figure 7m). Among them, a locus in chromosome 9
potentially associated with lateral vein interval showed a high value of
the Patterson's \emph{D} statistics (\emph{D} = 0.51; Figure S7), and a
region of loci potentially associated with stellate hair density in
chromosome 11 showed high local ancestry from \emph{Qd} (1.33
\emph{Qd}alleles; Figure 6b).

From 31 genes around the loci potentially associated with the four
traits, we obtained 39 GO terms, in which most frequently observed terms
were DNA binding and protein binding (Table S2). Among the 31 genes, we
obtained 10 protein names (Table S2), some of which were related to
drought stress resistance (Lorenzo et al. 2002; Kushiro et al. 2004),
disease resistance (Van Damme et al. 2009; Mokryakova et al. 2014), and
senescence (Engqvist et al. 2011).

\textbf{4 \textbar{} Discussion}

Adaptive introgression is known in some plant taxa, often resulting in
divergent ecotypes colonizing novel habitats (Suarez-Gonzalez et al.
2018c). In European and American oaks, divergence and introgression of
adaptive traits between interfertile species have been demonstrated. In
drought-averse \emph{Q. rubra} and drought-tolerant \emph{Q.
ellipsoidalis} (section \emph{Lobatae} ), parapatric populations of both
species are fixed for alternate alleles at a \emph{CONSTANS} -like gene,
which affected seedling survival in common gardens, and the \emph{Q.
ellipsoidalis} allele is introgressed into the \emph{Q. rubra} genome
depending on soil moisture conditions in their sympatric populations
(Lind-Riehl et al. 2014; Lind-Riehl and Gailing 2016; Khodwekar and
Gailing 2017). In \emph{Q. robur} and \emph{Q. petraea}
(section\emph{Quercus} ), introgression from \emph{Q. robur} contributes
to divergence among locally adapted \emph{Q. petraea} populations along
a gradient of climatic conditions (Leroy et al. 2020). In Asian oaks,
introgressed genomic regions were investigated in admixture
between\emph{Q. acutissima} and \emph{Q. variabillis} (section
\emph{Cerris} ), and admixed populations occupying similar habitats
tended to share the same introgressed regions (Fu et al. 2022). We also
found a candidate system for adaptive introgression, which consists of
\emph{Q. mongolica}var. \emph{crispula} (\emph{Qc} ) in inland habitats
and \emph{Q. dentata} (\emph{Qd} ) in coastal habitats in northern Japan
(Nagamitsu et al. 2019; Nagamitsu et al. 2020). In this system, we tried
to verify intensive introgression of \emph{Qd} alleles at loci
associated with traits adaptive to coastal environment.

Our previous study demonstrated environment-dependent introgression
from\emph{Qd} to \emph{Qc} in coastal habitats, resulting in \emph{Q.}
\selectlanguage{ngerman}×\emph{angustilepidota} (\emph{Qa} ), which are hybrids after the first
generation of backcross to \emph{Qc} (Nagamitsu et al. 2020). Our
present study using a comprehensive tree collection in common gardens
revealed a hybrid zone around the northern distributional limit
of\emph{Qd} , including both northern-edge \emph{Qd} populations admixed
with \emph{Qc} and coastal \emph{Qa} populations backcrossed to\emph{Qc}
. The observed genomic compositions indicated random admixture between
\emph{Qc} and \emph{Qd} , suggesting frequent gene flow and weak genetic
drift. The genotypic PCA indicated that some northern-edge\emph{Qd}
populations had unique genetic properties in some PCs with minor
contributions to genetic variations, suggesting the presence of weak
genetic structure in coastal habitats, probably due to fragmentation of
sand dunes and coastal cliffs in the distributional range margins. This
hybrid zone is likely to be maintained by the tension zone process
(Barton and Hewitt 1985). Near the northern distributional limit of
\emph{Qd} , gene flow between \emph{Qc} and\emph{Qd} seems relatively
frequent, because their flowering synchrony increases at higher
latitudes probably due to shorter periods of growing season (Shimizu et
al. 1995). Lower fitness of \emph{Qd} at higher latitudes in comparison
to fitness of \emph{Qc} coupled with both lower fitness of \emph{Qc} in
coastal habitats and lower fitness of \emph{Qd}in inland habitats may be
responsible for natural selection that maintains the tension zone.

Common garden experiments have revealed local adaptation in various
traits of oaks (Cavender-Bares and Ramírez-Valiente 2017). Local
adaptation of \emph{Qc} to climatic conditions was evident at regional
scales (Nagamitsu and Shuri 2021). The size of 30-year-old trees in our
common gardens indicated that the size was larger in \emph{Qc} trees
than in \emph{Qd} trees in the inland site but larger in \emph{Qd} trees
than in \emph{Qc} trees in the coastal site. This result demonstrates
both adaptation of \emph{Qc} to inland habitats and adaptation
of\emph{Qd} to coastal habitats. Although the plantation density and
forest management differed between the sites, smaller tree size in the
coastal site than in the inland site suggests higher stress in coastal
habitats than in inland habitats. In the coastal site, which was located
in north to the northern distributional limit of \emph{Qd} , the size
of\emph{Qd} trees seemed to be suppressed and to become similar to the
size of \emph{Qa} trees. This result implies that environmental
conditions at higher latitudes are maladaptive for \emph{Qd} , which is
distributed at lower latitudes than \emph{Qc} . These findings of local
adaptation are consistent with natural selection expected from the
tension zone process.

Morphological traits of leaves and shoots observed in each site
indicated more \emph{Qd} -like phenotypes as the \emph{Qd} ancestry
increased, suggesting that these traits have the genetic basis derived
from \emph{Qd} . Artificially crossed F\textsubscript{1} hybrids
between\emph{Qc} and \emph{Qd} exhibited pubescent shoots and
intermediate phenotypes of morphological traits of leaves and fruits
(Ubukata et al. 1996). These findings are consistent with our observed
phenotypes of admixed individuals. We categorized into observed traits
into taxon-specific and habitat-specific traits based on the results of
phenotypic PCA. The habitat-specific traits, of which phenotypes
of\emph{Qa} and \emph{Qd} trees in coastal habitats mainly differed from
those of \emph{Qc} trees in inland habitats, included the density of
stellate hairs (trichomes) on the lower leaf surface, leaf mass per area
(LMA), shoot diameter, and the number of axillary buds at bud-scale
scars. These traits seem relevant to the stress tolerance in coastal
habitats. Higher water absorption capacity of oak leaves with dense
trichomes is advantageous in soils with lower water storage capacity
(Fernández et al. 2014). Increased LMA in oaks are commonly observed in
response to higher intensity of drought and light (Sancho-Knapik et al.
2021). More buds at the lower part of shoots, such as bud-scale scars,
can compensate for bud mortality at the upper part of shoots caused by
salt spray and harsh wind in winter (Asai et al. 1986). These lines of
evidence suggest that \emph{Qd} -like phenotypes of the habitat-specific
traits are more adaptive in coastal environment. Thus, we expected
that\emph{Qd} alleles at loci associated with these traits are more
intensively introgressed than those in the genomic background.

The observed shift in phenotypic variation to more \emph{Qd} -like
phenotypes from inland to coastal habitats suggests phenotypic
plasticity in response to coastal stress. This phenotypic shift was
larger in \emph{Qa} trees than in \emph{Qd} trees, indicating lower
plasticity in more stress-tolerant \emph{Qd} trees. This result is
consistent with our previous findings that phenotypic plasticity in
natural populations in response to coastal stress, which was measured by
bud mortality, decreased with higher \emph{Qd} ancestry estimated from
nuclear microsatellite genotypes (Nagamitsu et al. 2019). Similar
features in phenotypic plasticity were evident in Mediterranean
oaks,\emph{Q. faginea} , which showed lower morphological and
physiological plasticity and more conservative resource-use strategy in
populations from drier and colder provenances than in those from mesic
and milder provenances (Solé-Medina et al. 2022). The phenotypic shift
to more\emph{Qd} -like phenotypes of \emph{Qa} trees may result from
induced expression of alleles introgressed from stress-tolerant
\emph{Qd} in response to coastal stress.

Introgression from \emph{Qd} to \emph{Qc} was evident from a positive
genome-wide estimate of the Patterson's \emph{D} statistic in our
previous study (Nagamitsu et al. 2020), and genome-wide patterns of
introgression in our present study using more loci indicated that the
extent of introgression were heterogeneous along chromosomes. However,
extensively introgressed regions were not always consistent with local
ancestry, probably due to incomplete linage sorting. To determine
introgressed genomic regions precisely, denser SNP loci mapped to
reference sequences of both parental taxa are necessary (Ai et al. 2022;
Wang et al. 2023). In spite of the low resolution in our present study,
the heterogeneity along chromosomes and various admixture patterns among
individuals may make our GWA analysis feasible owing to a diverse array
of recombinants.

We detected few loci consistently associated with traits in both SNP
genotypes and local ancestry, probably due to sparse distribution of the
available ddRAD loci, polygenic regulation of the examined traits, and
confounding factors with the measured phenotypes. In a model
plant,\emph{Oryza sativa} , for example, major genes of some agronomic
traits with alleles introgressed between rice varieties were clearly
detected from both association mapping using SNP genotypes and admixture
mapping using local ancestry (Zhao et al. 2010). In poplar, various
traits were measured in a common garden of \emph{Populus trichocarpa} ,
and associations of SNP genotypes with these traits were examined using
29k genome-wide loci (Mckown et al. 2014a; Mckown et al. 2014b).
Introgression of \emph{P. balsamifera} alleles to \emph{P.
trichocarpa}genome at some candidate genes around these trait-associated
loci were detected using denser loci around these candidate genes
(Suarez-Gonzalez et al. 2016). Furthermore, using 1169k loci, additional
introgressed regions were found from local ancestry inference and
admixture mapping (Suarez-Gonzalez et al. 2018a; Suarez-Gonzalez et al.
2018b). These studies imply that more loci are necessary to detect genes
responsible for adaptive traits in our oak system. In spite of the
limited power of our present study using 27k loci, we found several loci
potentially associated with one taxon-specific and three
habitat-specific traits. Two of these loci were located at intensively
introgressed regions in chromosomes 9 and 11.

GO terms of genes around the loci potentially associated with the four
traits suggest that these genes are related to signaling processes and
stress responses with functions of DNA binding and protein binding. Some
proteins of those genes were related to catabolic pathways of abscisic
acid (ABA) in response to drought stress (Kushiro et al. 2004) and ABA
transduction cascades in dormant seeds (Lorenzo et al. 2002),
immunophilin responses to bacterial invasion (Mokryakova et al. 2014)
and homoserine accumulation triggering pathogen resistance (Van Damme et
al. 2009), and catabolism of lysine during leaf senescence (Engqvist et
al. 2011). These findings imply that these genes can be related to the
response to coastal stress. In introgression between \emph{Q.
acutissima} and \emph{Q. variabilis} in China, not only genes involved
in signaling processes such as responses to abiotic stimulus, stress,
hormone, and pathogens but also \emph{cis} -regulatory elements in
untranslated regions were enriched in introgressed genomic regions (Fu
et al. 2022). To detect candidate genes responsible for the\emph{Qd}
-like phenotypes in response to coastal stress, not only morphological
and physiological traits but also gene expression patterns should be
investigated to select essential transcription factors and regulatory
elements underlying adaptive phenotypes.

In summary, we revealed a hybrid zone with random admixture
between\emph{Qc} and \emph{Qd} in the northernmost area of Hokkaido.
This hybrid zone is located in coastal habitats around the northern
distributional limit of \emph{Qd} , where \emph{Qd} -like phenotypes are
adaptive. We found heterogeneity of introgression along chromosomes and
various admixture patterns among individuals, but failed to detect loci
associated with the \emph{Qd} -like phenotypes. To elucidate adaptive
introgression suggested in our present study, more SNP genotypes at
denser loci in admixed populations and deeper knowledge of genetic basis
of adaptive traits are necessary.

\textbf{Acknowledgements}

We thank Kazuhito Kita and Nobuhiro Akashi (Hokkaido Research
Organization) for their help to investigate the inland and coastal
sites; Atsushi Nakanishi and Keiko Kitamura (Forest Research and
Management Organization) for their cooperation of the investigation;
Akiko Takazawa and Yasuyuki Komatsu for their assistance in the
laboratory work; Mineaki Aizawa (Utsunomiya University), and Ichiro
Tamaki (Gifu Prefectural Research Institute for Forests) for their
suggestions for improving the manuscript.

\textbf{Author contributions}

TN conceived the study. HS provided the inland and coastal sites. TN,
KU, and AI obtained and analyzed the data. TN drafted the manuscript.
KU, AI, and HS revised the manuscript.

\textbf{Conflict of Interest Statement}

The authors declare no conflict of interest.

\textbf{Funding Information}

This study was supported by JSPS KAKENHI grant no. 17K07859 to TN.

\textbf{Data Accessibility Statement}

The data of phenotypes and performance of trees in the common gardens
are given in Supporting Information (Data S1).

\textbf{References}

Ai W, Liu Y, Mei M, et al (2022) A chromosome\selectlanguage{english}-scale genome assembly of
the Mongolian oak (Quercus mongolica). Molecular Ecology Resources
22:2396--2410. doi: 10.1111/1755-0998.13616

Aizawa M, Maekawa K, Mochizuki H, Iizuka K (2021) Taxonomic revision of
quercus serrata subsp. mongolicoides. Acta Phytotaxonomica et
Geobotanica 72:113--123. doi: 10.18942/apg.202017

Alexander DH, Lange K (2011) Enhancements to the ADMIXTURE algorithm for
individual ancestry estimation. BMC Bioinformatics 12:246. doi:
10.1186/1471-2105-12-246

Anderson E (1953) Introgressive hybridization. Biological reviews of the
Cambridge Philosophical Society 28:280--307.

Anderson E, Stebbins GL (1954) Hybridization as an evolutionary
stimulus. Evolution 8:378--388.

Arnold BJ, Lahner B, DaCosta JM, et al (2016) Borrowed alleles and
convergence in serpentine adaptation. Proceedings of the National
Academy of Sciences of the United States of America 113:8320--8325. doi:
10.1073/pnas.1600405113

Asai T, Shinmura Y, Usui G (1986) Mortality factor of the winter buds of
Quercus dentata and Quercus mongolica var. grosseserrata in natural
coastal forests of northern Hokkaido. Journal of Japanese Forest Society
68:368--374. doi: 10.11519/jjfs1953.68.9\_368

Barton NH, Hewitt GM (1985) Analysis of hybrid zones. Annual Review of
Ecology and Systematics 16:113--148. doi:
10.1146/annurev.es.16.110185.000553

Bickford CP (2016) Ecophysiology of leaf trichomes. Functional Plant
Biology 43:807--814. doi: 10.1071/FP16095

Bresadola L, Caseys C, Castiglione S, et al (2019) Admixture mapping in
interspecific Populus hybrids identifies classes of genomic
architectures for phytochemical, morphological and growth traits. New
Phytologist 223:2076--2089. doi: 10.1111/nph.15930

Browning BL, Tian X, Zhou Y, Browning SR (2021) Fast two-stage phasing
of large-scale sequence data. The American Journal of Human Genetics
108:1880--1890. doi: 10.1016/j.ajhg.2021.08.005

Browning BL, Zhou Y, Browning SR (2018) A one-penny imputed genome from
next-generation reference panels. American Journal of Human Genetics
103:338--348. doi: 10.1016/j.ajhg.2018.07.015

Buerkle CA, Lexer C (2008) Admixture as the basis for genetic mapping.
Trends in Ecology and Evolution 23:686--694. doi:
10.1016/j.tree.2008.07.008

Burke JM, Arnold ML (2001) Genetics and the fitness of hybrids. Annual
Review of Genetics 35:31--52.

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. doi:
10.1111/nph.15450

Cavender-Bares J, Ramirez-Valiente JA (2017) Physiological evidence from
common garden experiments for local adaptation and adaptive plasticity
to climate in American live oaks (Quercus section Virentes):
Implications for conservation under global change. In: Gil-Pelegrin E,
Peguero-Pina J, Sancho-Knapik D (eds) Oaks Physiological Ecology.
Exploring the Functional Diversity of Genus Quercus L. Tree Physiology,
vol 7. Springer, Cham, pp 107--135

Chhatre VE, Evans LM, DiFazio SP, Keller SR (2018) Adaptive
introgression and maintenance of a trispecies hybrid complex in
range\selectlanguage{english}-edge populations of Populus. Molecular Ecology 27:4820--4838. doi:
10.1111/mec.14820

Ciccarelli D, Bona C (2022) Exploring the functional strategies adopted
by coastal plants along an ecological gradient using morpho-functional
traits. Estuaries and Coasts 45:114--129. doi:
10.1007/s12237-021-00945-y

Danecek P, Auton A, Abecasis G, et al (2011) The variant call format and
VCFtools. Bioinformatics 27:2156--2158. doi:
10.1093/bioinformatics/btr330

Denk T, Grimm GW, Manos PS, et al (2017) An updated infrageneric
classification of the oaks: review of previous taxonomic schemes and
synthesis of evolutionary patterns. In: Oaks Physiological Ecology.
Exploring the Functional Diversity of Genus Quercus.

Dias-Alves T, Mairal J, Blum MGB (2018) Loter: A software package to
infer local ancestry for a wide range of species. Molecular Biology and
Evolution 35:2318--2326. doi: 10.1093/molbev/msy126

Endelman JB (2011) Ridge regression and other kernels for genomic
selection with R package rrBLUP. The Plant Genome 4:250--255. doi:
10.3835/plantgenome2011.08.0024

Engqvist MKM, Kuhn A, Wienstroer J, et al (2011) Plant
D-2-hydroxyglutarate dehydrogenase participates in the catabolism of
lysine especially during senescence. Journal of Biological Chemistry
286:11382--11390. doi: 10.1074/jbc.M110.194175

Fernandez V, Sancho-Knapik D, Guzman P, et al (2014) Wettability,
polarity, and water absorption of holm oak leaves: Effect of leaf side
and age. Plant Physiology 166:168--180. doi: 10.1104/pp.114.242040

Fitzpatrick BM (2012) Estimating ancestry and heterozygosity of hybrids
using molecular markers. BMC Evolutionary Biology 12:131. doi:
10.1186/1471-2148-12-131

Fu R, Zhu Y, Liu Y, et al (2022) Genome-wide analyses of introgression
between two sympatric Asian oak species. Nature Ecology \& Evolution
6:924--935. doi: 10.1038/s41559-022-01754-7

Goulet BE, Roda F, Hopkins R (2017) Hybridization in plants: Old ideas,
new techniques. Plant Physiology 173:65--78. doi: 10.1104/pp.16.01340

Hamilton JA, Miller JM (2016) Adaptive introgression as a resource for
management and genetic conservation in a changing climate. Conservation
Biology 30:33--41. doi: 10.1111/cobi.12574

Hesp PA (1991) Ecological processes and plant adaptations on coastal
dunes. Journal of Arid Environments 21:165--191. doi:
10.1016/S0140-1963(18)30681-5

Hipp AL, Manos PS, Hahn M, et al (2020) Genomic landscape of the global
oak phylogeny. New Phytologist 226:1198--1212. doi: 10.1111/nph.16162

Ishida TA, Hattori K, Sato H, Kimura MT (2003) Differentiation and
hybridization between Quercus crispula and Q. dentata (Fagaceae):
Insights from morphological traits, amplified fragment length
polymorphism markers, and leafminer composition. American Journal of
Botany 90:769--776. doi: 10.3732/ajb.90.5.769

Khodwekar S, Gailing O (2017) Evidence for environment-dependent
introgression of adaptive genes between two red oak species with
different drought adaptations. American Journal of Botany
104:1088--1098. doi: 10.3732/ajb.1700060

Kim BY, Wei X, Fitz\selectlanguage{english}-Gibbon S, et al (2018) RADseq data reveal ancient,
but not pervasive, introgression between Californian tree and scrub oak
species (Quercus sect. Quercus, Fagaceae). Molecular Ecology
27:4556--4571. doi: 10.1111/mec.14869

Kushiro T, Okamoto M, Nakabayashi K, et al (2004) The Arabidopsis
cytochrome P450 CYP707A encodes ABA 8\selectlanguage{english}'-hydroxylases: Key enzymes in ABA
catabolism. EMBO Journal 23:1647--1656. doi: 10.1038/sj.emboj.7600121

Leroy T, Louvet JM, Lalanne C, et al (2020) Adaptive introgression as a
driver of local adaptation to climate in European white oaks. New
Phytologist 226:1171--1182. doi: 10.1111/nph.16095

Lind-Riehl J, Gailing O (2016) Adaptive variation and introgression of a
CONSTANS-like gene in North American red oaks. Forests 8:3. doi:
10.3390/f8010003

Lind-Riehl JF, Sullivan AR, Gailing O (2014) Evidence for selection on a
CONSTANS-like gene between two red oak species. Annals of Botany
113:967--975. doi: 10.1093/aob/mcu019

Lindtke D, Gonzalez-Martinez SC, Macaya-Sanz D, Lexer C (2013) Admixture
mapping of quantitative traits in Populus hybrid zones: Power and
limitations. Heredity 111:474--485. doi: 10.1038/hdy.2013.69

Lorenzo O, Nicolas C, Nicolas G, Rodriguez D (2002) Molecular cloning of
a functional protein phosphatase 2C (FsPP2C2) with unusual features and
synergistically up-regulated by ABA and calcium in dormant seeds of
Fagus sylvatica. Physiologia Plantarum 114:482--490. doi:
10.1034/j.1399-3054.2002.1140318.x

Ma Y, Wang J, Hu Q, et al (2019) Ancient introgression drives adaptation
to cooler and drier mountain habitats in a cypress species complex.
Communications Biology 2:213. doi: 10.1038/s42003-019-0445-z

Martin NH, Bouck AC, Arnold ML (2006) Detecting adaptive trait
introgression between Iris fulva and I. brevicaulis in highly selective
field conditions. Genetics 172:2481--2489. doi:
10.1534/genetics.105.053538

Martin SH, Davey JW, Jiggins CD (2015) Evaluating the use of ABBA-BABA
statistics to locate introgressed loci. Molecular Biology and Evolution
32:244--257. doi: 10.1093/molbev/msu269

Martin SH, Van Belleghem SM (2017) Exploring evolutionary relationships
across the genome using topology weighting. Genetics 206:429--438. doi:
10.1534/genetics.116.194720

Matsumoto A, Kawahara T, Kanazashi A, et al (2009) Differentiation of
three closely related Japanese oak species and detection of
interspecific hybrids using AFLP markers. Botany 87:145--153. doi:
10.1139/B08-121

Mckown AD, Guy RD, Klap\selectlanguage{polish}š\selectlanguage{english}t\selectlanguage{polish}ě \selectlanguage{english}J, et al (2014a) Geographical and
environmental gradients shape phenotypic trait variation and genetic
structure in Populus trichocarpa. New Phytologist 201:1263--1276. doi:
10.1111/nph.12601

Mckown AD, Kl\selectlanguage{ngerman}áp\selectlanguage{polish}š\selectlanguage{english}t\selectlanguage{polish}ě \selectlanguage{english}J, Guy RD, et al (2014b) Genome-wide association
implicates numerous genes underlying ecological trait variation in
natural populations of Populus trichocarpa. New Phytologist
203:535--553. doi: 10.1111/nph.12815

Menon M, Bagley JC, Page GFM, et al (2021) Adaptive evolution in a
conifer hybrid zone is driven by a mosaic of recently introgressed and
background genetic variants. Communications Biology 4:1--14. doi:
10.1038/s42003-020-01632-7

Mokryakova M V., Pogorelko G V., Bruskin SA, et al (2014) The role of
peptidyl-prolyl cis/trans isomerase genes of Arabidopsis thaliana in
plant defense during the course of Xanthomonas campestris infection.
Russian Journal of Genetics 50:140--148. doi: 10.1134/S1022795414020100

Nagamitsu T, Shimizu H, Aizawa M, Nakanishi A (2019) An admixture of
Quercus dentata in the coastal ecotype of Q. mongolica var. crispula in
northern Hokkaido and genetic and environmental effects on their traits.
Journal of Plant Research 132:211--222. doi: 10.1007/s10265-018-01079-2

Nagamitsu T, Shuri K (2021) Seed transfer across geographic regions in
different climates leads to reduced tree growth and genetic admixture in
Quercus mongolica var. crispula. Forest Ecology and Management
482:118787. doi: 10.1016/j.foreco.2020.118787

Nagamitsu T, Uchiyama K, Izuno A, et al (2020) Environment-dependent
introgression from Quercus dentata to a coastal ecotype of Quercus
mongolica var. crispula in northern Japan. New Phytologist
226:1018--1028. doi: 10.1111/nph.16131

Ohba H (2006) Fagaceae. In: Iwatsuki K, Boufford DE, Ohba H (eds) Flora
of Japan. Volume IIa. Kodansha, Tokyo, pp 42--60

Ortego J, Gugger PF, Riordan EC, Sork VL (2014) Influence of climatic
niche suitability and geographical overlap on hybridization patterns
among southern Californian oaks. Journal of Biogeography 41:1895--1908.
doi: 10.1111/jbi.12334

Peterson BK, Weber JN, Kay EH, et al (2012) Double digest RADseq: An
inexpensive method for de novo SNP discovery and genotyping in model and
non-model species. PLoS ONE 7:e37135. doi: 10.1371/journal.pone.0037135

Petit RJ, Bod\selectlanguage{ngerman}énès C, Ducousso A, et al (2004) Hybridization as a
mechanism of invasion in oaks. New Phytologist 161:151--164. doi:
10.1046/j.1469-8137.2003.00944.x

Plomion C, Aury JM, Amselem J, et al (2018) Oak genome reveals facets of
long lifespan. Nature Plants 4:440--452. doi: 10.1038/s41477-018-0172-3

Puritz JB, Hollenbeck CM, Gold JR (2014) dDocent: a RADseq,
variant-calling pipeline designed for population genomics of non-model
organisms. PeerJ 2:e431. doi: 10.7717/peerj.431

R\_Core\_Team (2019) A language and environment for statistical
computing. R Foundation for Statistical Computing

Ramírez-Valiente JA, Sánchez-Gómez D, Aranda I, Valladares F (2010)
Phenotypic plasticity and local adaptation in leaf ecophysiological
traits of 13 contrasting cork oak populations under different water
availabilities. Tree Physiology 30:618--627. doi:
10.1093/treephys/tpq013

Rendón-Anaya M, Wilson J, Sveinsson S, et al (2021) Adaptive
introgression facilitates adaptation to high latitudes in European aspen
(Populus tremula L.). Molecular Biology and Evolution 38:5034--5050.
doi: 10.1093/molbev/msab229

Rieseberg LH, Buerkle CA (2002) Genetic mapping in hybrid zones.
American Naturalist 159:S37--S50. doi: 10.2307/3078920

Rieseberg LH, Carney SE (1998) Plant hybridization. New Phytologist
140:599--624. doi: 10.1046/j.1469-8137.1998.00315.x

Rieseberg LH, Kim S-C, Randell RA, et al (2007) Hybridization and the
colonization of novel habitats by annual sunflowers. Genetica
129:149--165. doi: 10.1007/s10709-006-9011-y

Riordan EC, Gugger PF, Ortego J, et al (2016) Association of genetic and
phenotypic variability with geography and climate in three southern
California oaks. American Journal of Botany 103:73--85. doi:
10.3732/ajb.1500135

Sancho-Knapik D, Escudero A, Mediavilla S, et al (2021) Deciduous and
evergreen oaks show contrasting adaptive responses in leaf mass per area
across environments. New Phytologist 230:521--534. doi:
10.1111/nph.17151

Shimizu H, Kikuchi K, Yamada K (1995) Local variation of bud number on
axillary buds of bud-scales of Quercus dentata in coastal forest along
Japan sea of Hokkaido. Trans Hokkaido Branch Jpn For Soc 43:140--142.

Solé-Medina A, Robledo-Arnuncio JJ, Ramírez-Valiente JA (2022)
Multi-trait genetic variation in resource-use strategies and phenotypic
plasticity correlates with local climate across the range of a
Mediterranean oak (Quercus faginea). New Phytologist 234:462--478. doi:
10.1111/nph.17968

Sork VL (2018) Genomic studies of local adaptation in natural plant
populations. Journal of Heredity 109:3--15. doi: 10.1093/jhered/esx091

Sork VL, Cokus SJ, Fitz-Gibbon ST, et al (2022) High-quality genome and
methylomes illustrate features underlying evolutionary success of oaks.
Nature Communications 13:1--15. doi: 10.1038/s41467-022-29584-y

Sork VL, Riordan E, Gugger PF, et al (2016) Phylogeny and introgression
of California scrub white oaks (Quercus section Quercus). International
Oaks 27:1--14.

Stacklies W, Redestig H, Scholz M, et al (2007) pcaMethods: A
bioconductor package providing PCA methods for incomplete data.
Bioinformatics 23:1164--1167. doi: 10.1093/bioinformatics/btm069

Stebbins GL, Matzke EB, Epling C (1947) Hybridization in a population of
Quercus marilandica and Quercus ilicifolia. Evolution 1:79. doi:
10.2307/2405406

Suarez-Gonzalez A, Hefer CA, Christe C, et al (2016) Genomic and
functional approaches reveal a case of adaptive introgression from
Populus balsamifera (balsam poplar) in P. trichocarpa (black
cottonwood). Molecular Ecology 25:2427--2442. doi: 10.1111/mec.13539

Suarez-Gonzalez A, Hefer CA, Lexer C, et al (2018a) Introgression from
Populus balsamifera underlies adaptively significant variation and range
boundaries in P. trichocarpa. New Phytologist 217:416--427. doi:
10.1111/nph.14779

Suarez-Gonzalez A, Hefer CA, Lexer C, et al (2018b) Scale and direction
of adaptive introgression between black cottonwood (Populus trichocarpa)
and balsam poplar (P. balsamifera). Molecular Ecology 27:1667--1680.
doi: 10.1111/mec.14561

Suarez-Gonzalez A, Lexer C, Cronk QCB (2018c) Adaptive introgression: A
plant perspective. Biology Letters 14:20170688. doi:
10.1098/rsbl.2017.0688

Ubukata M, Koono K, Iizuka K (1996) Morphological characteristics of
Quercus crispula x dentata hybrids. Trans Hokkaido Branch Jpn For Soc
44:113--116.

Van Damme M, Zeilmaker T, Elberse J, et al (2009) Downy mildew
resistance in arabidopsis by mutation of Homoserine Kinase. Plant Cell
21:2179--2189. doi: 10.1105/tpc.109.066811

Wang W, He X, Yan X, et al (2023) Chromosome\selectlanguage{english}-scale genome assembly and
insights into the metabolome and gene regulation of leaf color
transition in an important oak species, Quercus dentata. New Phytologist
238:2016--2032. doi: 10.1111/nph.18814

Whitney KD, Randell RA, Rieseberg LH (2010) Adaptive introgression of
abiotic tolerance traits in the sunflower Helianthus annuus. New
Phytologist 187:230--239. doi: 10.1111/j.1469-8137.2010.03234.x

Yu J, Pressoir G, Briggs WH, et al (2006) A unified mixed-model method
for association mapping that accounts for multiple levels of
relatedness. Nature Genetics 38:203--208. doi: 10.1038/ng1702

Zhao K, Wright M, Kimball J, et al (2010) Genomic diversity and
introgression in O. sativa reveal the impact of domestication and
breeding on the rice genome. PLoS ONE 5:e10780. doi:
10.1371/journal.pone.0010780

Table and Figure Legend

\textbf{Figure 1.} Locations of provenances of \emph{Quercus
mongolica}var. \emph{crispula} (\emph{Qc} , red), \emph{Q.}
x\emph{angustilepidota} (green), and \emph{Q. dentata} (\emph{Qd} ,
blue), which were planted in the inland site (a) and the coastal site
(b). Provenances are categorized into southwestern (inverse triangles),
northwestern (triangles), northern (diamonds), eastern (squares)
districts of the northernmost area of Hokkaido and other areas of
Hokkaido (circles). Enlarged symbols indicate reference individuals
of\emph{Qc} and \emph{Qd} for inferring local ancestry. Locations of the
inland site (red arrow head) and the coastal site (blue arrow head) are
shown.

\textbf{Figure 2.} \emph{Quercus dentata} ancestry (\emph{S} ) and
inter-ancestry heterozygosity (\emph{H} ) of sampled trees in the inland
site (a) and the coastal site (b). Ancestral populations of \emph{Q.
mongolica} var. \emph{crispula} are \emph{S} = 0 and \emph{H} = 0, and
those of \emph{Q. dentata} are \emph{S} = 0 and \emph{H} = 0. Solid
lines indicate backcrossing from F\textsubscript{1} hybrids (\emph{S} =
0.5 and \emph{H} = 1) to ancestral populations, and a dotted curve
indicate random mating between the ancestral populations. Symbols are
shown in the same way as Figure 1.

\textbf{Figure 3.} Principal component analysis for phenotypes of eight
morphological traits of leaves and shoots. Proportion (\%) of variances
for the first five principal components (PCs) and vectors (arrows) of
the traits contributing to the first and second PCs are shown (a).
Coordination of sampled trees in the inland site (b) and the coastal
site (c) on the first and second PCs are shown. Symbols are shown in the
same way as Figure 1.

\textbf{Figure 4.} Relationship between \emph{Quercus dentata} ancestry
and the first principal component of phenotypic variation of sampled
trees in the inland site (a) and the coastal site (b). Solid and dotted
lines are fitted to observations in the focal site and the other site,
respectively, using the lowess method. Symbols are shown in the same way
as Figure 1.

\textbf{Figure 5.} Relationship between mean values of the first
principal component of phenotypic variation of sampled trees and mean
values of the stem basal area of recorded trees at plantation rows in
the inland site (a) and the coastal site (b). Solid and dotted lines are
fitted to observations in the focal site and the other site,
respectively, using the lowess method. Symbols are shown in the same way
as Figure 1.

\textbf{Figure 6.} Genome-wide patterns of local ancestry (a) and mean
number of \emph{Quercus dentata} alleles in admixed individuals (b).
Along vertical axis (a), \emph{Q. dentata} reference individuals
(white), admixed individuals sorted by \emph{Quercus dentata} ancestry,
and \emph{Q. mongolica} var. \emph{crispula} reference individuals
(black) are arranged. White, black, and grey bars (a) indicate
homozygote of \emph{Q. dentata} alleles, homozygote of \emph{Q.
mongolica} var. \emph{crispula} alleles, and heterozygote of both
alleles, respectively, in admixed individuals. Blue vertical lines
indicate borders of chromosomes, and red vertical lines indicate
positions of loci potentially associated with traits shown in Figure 7
and Table S2.

\textbf{Figure 7.} Manhattan plots of likelihood of genome-wide
association with eight morphological traits of leaves and shoots. Left
panels (a, c, e, g, i, k, m, o) are association of SNP genotypes
(association mapping), and right panels (b, d, f, h, j, l, n, p) are
association of local ancestry (admixture mapping). Taxon-specific
traits: relative lead width (a, b), lateral vein interval (c, d), tooth
angle (e, f), and stellate hair length (g, h); and habitat-specific
traits: stellate hair density (i, j), leaf mass per area (k, l), shoot
diameter (m, n), and scale bud number (o, p) are shown. Black and grey
circles indicate loci in different chromosomes, and green circles
indicate loci potentially associated with traits (FDR \textless{} 0.15).

Supporting Information

Table S1. Provenances of sampled trees in inland and coastal sites.

Figure S1. Photographs of inland site (a) and coastal site (b) in 2017.

Figure S2. Locations and sizes of trees planted in three blocks (a--c)
in inland site. Solid circles, circles with crosses, and open circles
indicate sampled reference, sampled non-reference, and non-sampled
living trees, respectively. Crosses, + and x, indicate dead and thinned
trees, respectively. Triangles indicate non-experimental trees. Colors
indicate taxa: \emph{Quercus mongolica} var. \emph{crispula}
(red),\emph{Q.} x \emph{angustilepidota} (green), and \emph{Q.
dentata}(blue).

Figure S3. Locations and sizes of trees planted in two blocks (a, b) in
coastal site. Solid circles, circles with crosses, and open circles
indicate sampled reference, sampled non-reference, and non-sampled
living trees, respectively. Crosses, + and x, indicate dead and thinned
trees, respectively. Triangles indicate non-experimental trees. Colors
indicate taxa: \emph{Quercus mongolica} var. \emph{crispula}
(red),\emph{Q.} x \emph{angustilepidota} (green), and \emph{Q.
dentata}(blue).

Figure S4. Measurements for morphological traits of leaves and shoots.

Figure S5. Bayesian clustering of SNP genotypes at ddRAD loci of sampled
trees in both inland and coastal sites. Log likelihood when the number
of clusters (\emph{K} ) are 1, 2, and 3 (a). Genetic differentiation
(\emph{F} \textsubscript{ST}) between clusters when \emph{K} = 2
and\emph{K} = 3 (b). Bar plots of ancestry proportions of trees sorted
by the ancestry proportion of clusters, when \emph{K} = 2 (c) and
\emph{K}= 3 (d). Colors correspond to clusters (b, c, d). Colors of
circles indicate taxa: \emph{Quercus mongolica} var. \emph{crispula}
(red),\emph{Q.} x \emph{angustilepidota} (green), and \emph{Q.
dentata}(blue), and enlarged symbols indicate reference samples (e).

Figure S6. Frequency distributions (scatter plots) and correlations
(histograms) for pairs of shoot hair (0: absent and 1: present) and
eight morphological traits: relative leaf width, lateral vein interval
(mm), tooth angle (\selectlanguage{english}@), loge length (mm) and log\textsubscript{e} density
(mm\textsuperscript{--2}) of stellate hair on lower leaf surface, lead
mass per area (mg mm\textsuperscript{--2}), shoot diameter (mm), and
log\textsubscript{e} number of axillary buds at bud-scale scars. Numbers
indicate Kendall's correlation coefficients, and red lines indicate
linear regression lines.

Figure S7. Patterson's \emph{D} statistics at neighboring 101 ddRAD loci
along chromosomes. Positive \emph{D} values indicate introgression
from\emph{Quercus dentata} to \emph{Q. mongolica} var. \emph{crispula} .
Red vertical lines indicate positions of trait-associated loci shown in
Figure 7 and Table S2.

Figure S8. Phenotypic values of four taxon-specific traits in inland
(left) and coastal (middle) sites and those values scaled in each site
and pooled together (right) in relation to \emph{Quercus
dentata}ancestry. Solid and dotted lines are fitted to observations in
the focal site and the other site, respectively, using the lowess
method. Symbols are shown in the same way as Figure 1.

Figure S9. Phenotypic values of four habitat-specific traits in inland
(left) and coastal (middle) sites and those values scaled in each site
and pooled together (right) in relation to \emph{Quercus
dentata}ancestry. Solid and dotted lines are fitted to observations in
the focal site and the other site, respectively, using the lowess
method. Symbols are shown in the same way as Figure 1.

Figure S10. Principal component analysis for SNP genotypes at ddRAD loci
of sampled trees in both inland and coastal sites. Proportion (\%) of
variances for the first 10 principal components (PCs) are shown (a).
Coordination of sampled trees on the first and second PCs (b), the first
and third PCs (c), and the first and fourth PCs (d) are shown. Symbols
are shown in the same way as Figure 1.

Figure S11. Q-Q plots of expected and observed
--log\textsubscript{10}\emph{p} values for four taxon-specific traits.
Red lines indicate identical expected and observed
--log\textsubscript{10} \emph{p} values.

Figure S12. Q-Q plots of expected and observed
--log\textsubscript{10}\emph{p} values for four habitat-specific traits.
Red lines indicate identical expected and observed
--log\textsubscript{10} \emph{p} values.

Table S2. ddRAD loci of SNP genotypes (G) and local ancestry (A)
associated with traits (FDR \textless{} 0.15) and proteins nearest to or
between trait-associated loci obtained from PM1N v2.3 \emph{Quercus
robur} annotation database.

Data S1. Data of phenotypes of morphological traits of leaves and shoots
and performance of trees in plantation rows in inland and coastal common
gardens.

\selectlanguage{english}
\FloatBarrier
\end{document}

