\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[ngerman,greek,english]{babel}



\usepackage{float}








\begin{document}

\title{Environmental and landscape genomics of Picea sitchensis (Sitka spruce)
reveals clusters of loci involved with both resistance to wintering and
increases in height.}



\author[1]{Tomas Byrne}%
\author[1]{Niall Farrelly}%
\author[2]{Colin Kelleher}%
\author[3]{Trevor Hodkinson}%
\author[1]{Stephen Byrne}%
\author[1]{Susanne Barth}%
\affil[1]{Teagasc}%
\affil[2]{National Botanic Gardens}%
\affil[3]{Trinity College Dublin}%


\vspace{-1em}



  
  \date{May 23, 2024}


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





\selectlanguage{english}
\begin{abstract}
Picea sitchensis (Sitka spruce) occupies a large number of host
environments across the Pacific Northwest, stretching from California to
Alaska. Populations have become locally adapted to their environmental
niche rather than undergoing speciation and these local adaptions can be
discovered and utilised for breeding programs. These local adaptions may
contain traits which can confer increased resilience against climate
change. Here we employed Genome Wide Association Studies (GWAS),
Genotype Environment Analysis (GEA) and Canonical correlation analysis
to discover adaptive loci. A total of 593 adaptive loci were discovered
from the GEA. The GWAS found eight loci involved with height and one
with diameter at breast height. The canonical correlation analysis
allowed for the discovery of positive and negative correlated loci with
traits or environmental variables. We discovered 121 loci positively
correlated with environmental values and 10 loci positively correlated
with phenotypic traits. We discovered one cluster of adaptive loci
involved with snowy winter conditions and two clusters showing
adaptation to southern locations characterised with milder conditions.
The allelic distribution of positive and negatively correlated loci
involved with height showed a north south divide with loci positively
correlated with height being distributed below a latitude of 50oN. This
shows adaptations for height are generally found south of the 50th
latitude. We have discovered signs of local adaptation, with Minor
Allele Frequency (MAF) differing between northern and southern
populations.%
\end{abstract}%



\sloppy


\textbf{Introduction}

Low linkage disequilibrium (LD) is characteristic of conifer genomes and
the independence of genes contributes to the adaptation patterns seen in
conifers (Prunier et al., 2016). Angiosperms may undergo allopatric
speciation by adaptation leading to reproductive isolation, however,
conifers are less likely to do this due to low LD and large genome size
(Ahuja \& Neale, 2005; Pavy et al., 2012). For this reason, local
adaptation to host environments is quite common in conifers with many
alleles with small effects driving the adaptation (De La Torre et al.,
2019; Hornoy et al., 2015). In \emph{Picea sitchensis} (Bong.) Carr.
(Sitka spruce) this has been demonstrated using phenotype data where bud
set timing was adaptive to local climate types (Mimura \& Aitken, 2010).
In Loblolly pine, \emph{Pinus taeda} L., climate adaptation at a local
level is due to these subtle changes in the frequency of alleles of
moderate/small effect (De La Torre et al., 2019) and in Lodgepole
pine,\emph{Pinus contortus} Douglas, climactic adaptations have been
shown to be characterised by a slow rate of recombination (Lotterhos et
al., 2018). Conifer is generally slow and characterised by low changes
in frequency of many alleles, which is likely to affect the diversity of
local adaptations found in these types of conifer studies. The long
lifecycle of conifers results in low recombination cycles, contributing
to slow adaptation.

Evergreen conifers have multiple mechanisms to survive colder harsh
winter environments (Chang et al., 2021; Senser \& Beck, 1984). Low
temperatures can reduce the rate of enzymatic activity, disrupt
metabolic processes and inhibit protein repair systems (Chang et al.,
2021). To overcome this, some conifers have developed cold hardiness,
which can stabilise and change membrane lipid composition which function
to reduce the build-up of ice crystals in extracellular spaces (Chang et
al., 2021; Senser \& Beck, 1984). Cold hardiness is reported in some
white spruce populations, a close relative of Sitka spruce
(Sebastian-Azcona et al., 2018). Growth cessation during colder periods
is an adaptation linked to shorter diurnal ranges. This is seen in Sitka
spruce with delayed bud set due to wintering suggesting growth cessation
in buds (Mimura \& Aitken, 2010).

Low water stress, heat stress and solar radiation all combine to create
drought stress. Thus, drought stress is a complex environmental stressor
which requires a suite of differing adaptations (Brodribb et al., 2014).
Conifers exhibit two diverging pathways to adapt to low water stress,
one involves raised levels of abscisic acid (ABA) to close stomata and
another method involves leaf desiccation (Teskey et al., 2015).
Transpiration cooling is important for heat stress however, but stomatal
closure can increase the severity of heat stress (Teskey et al., 2015).
It is known that reductions in total seasonal solar radiation correlate
to a reduction in tree height, so generally it is not considered a
stressor on its own unless it exacerbates drought stress through
evapotranspiration (Strand et al., 2006).

Sitka spruce occupies a large number of geographic and climatic niches
along the Pacific Northwest of North America. This large gradient in
environments requires specific local adaptations across the population
(Gapare et al., 2005). The IUFRO (International Union of Forest Research
Organisations) population was established in 1968 and 1970 by Matthews
and Fletcher, representing the natural range of Sitka, collected at
intervals of 50km and representing \textgreater{}80 provenances. (Sype,
1990). The collection includes samples over 3000 km in extent from
Alaska to California. Elevation ranges from sea level to 671m, and an
maximum distance from coast of 164m. Annual precipitation ranges from
4015 mm to 621 mm (Fick \& Hijmans, 2017). Mean annual temperature
ranges from -7.3\selectlanguage{ngerman}°C to 17.4°C (Fick \& Hijmans, 2017). Mean diurnal range
is 10.8hrs to 4.43hrs per day (Fick \& Hijmans, 2017). This highlights
the range of environments present in this data set making it key for the
study of local adaptation in conifers. This diverse and large range of
environments likely leads to adaptations at the local level driven by
shifts in allele frequencies, a common pattern in conifers (De La Torre
et al., 2019). To discover the genetic bases of these adaptations, a
number of methods have been developed. Genome wide association studies
take trait data and genotype data and associated molecular markers with
a trait (Casola, 2019; Chen et al., 2021; Hiraoka et al., 2018). This is
useful for the discovery of the molecular function of traits or for the
identification of molecular markers for the deployment in selection
programs as has been done in \emph{Picea abies} (L.) H. Karst. (Chen et
al., 2021). Similar techniques can be used to associate genes with an
environment using genotype-environment association. In conifers, trait
association has been completed in Lodgepole pine and Loblolly pine by
associating climactic variables with loci, elucidating mechanisms of
adaptation (De La Torre et al., 2019; Lotterhos et al., 2018). The
ability to associate genotype, phenotype and environment is important in
these studies, as it indicates a pathway towards adaptation. While not
previously undertaken in conifers, Canonical correlation analysis
(CANCOR) has allowed for the discovery of loci that correlate with
phenotypic growth due to environmental conditions in the forage
grass\emph{Lolium perenne} (Blanco-Pastor et al., 2021) (also
outbreeding and wind pollinated). This analysis merges data to discover
loci responsible for improved phenotypes for the environment, which is
key for breeders.

In this chapter we utilised data from a diverse collection of Sitka
spruce, available in an arboretum collection in Ireland, representing
the entire natural range of the species. This collection consists of 80
provenances with roughly 15 trees per provenance. In total 1177
genotypes were available for the study. Here we have utilised GEA, GWAS
and CANCOR to discover loci responsible for specific traits and
adaptations. We implemented CANCOR analysis to see the phenotypic
effects of environments on genotypes. These loci can then be traced
through populations, elucidating the driving forces of adaptation in
Sitka spruce.

\textbf{Materials and Methods}

\textbf{Genotype Data}

We had previously generated GBS data to assess genetic diversity and
structure of a diverse population of \emph{Picea sitchensis} (Byrne et
al., 2022) so were able to utilize that dataset for this study. GBS data
from 1177 genotypes of the IUFRO population, representing 80 sampling
points, were therefore obtained from NCBI BioProject PRJNA852515 (Byrne
et al., 2022). This data was generated with 150bp paired-end GBS using
Pst1-Apek1 following (Elshire et al., 2011). Sequences were aligned to
the Q903-v1-1000 plus Sitka spruce genome (GCA\_010110895.2) (Gagalova
et al., 2022) using BWA-mem with default parameters (Li \& Durbin,
2009). Variant calling was completed using SAMtools v1.9 mpileup (Li et
al., 2009). VCFtools was used for filtering (Danecek et al., 2021).
Indels were removed, multi-allelic SNPs were removed, minimum genotype
quality was filtered to 20, read depth was filtered to 5, max missing
data was filtered to 0.7 and MAF was set to 0.05. This resulted in 66974
loci which were transformed into an allelic matrix for further study.

\textbf{Phenotype Data}

Phenotypic data were collected in November 2022 from each of the 1177
genotypes in the IUFRO population at John F. Kennedy Arboretum in
Wexford, Ireland (52.315N, -6.941W). Diameter at breast (DBH) height was
taken at 1.3m from the base of the tree using a steel DBH tape. Height
was recorded for each tree using a vertex transponder (Hagloff V) to
represent the distance from ground level to tree tip.

\textbf{Environmental data}

All data (table 1) was downloaded in raster format and data for each
location of the IUFRO population was extracted using the Point Sampling
Tool in qGIS (Jurgiel, 2020; QGIS.org, 2022). Climatic data was
downloaded from worldclim.org (Fick \& Hijmans, 2017). This data
consisted of readings from 1970-2000 with a raster resolution of
1km\textsuperscript{2}. Annual snow cover data was obtained from the
Moderate Resolution Imaging Spectroradiometer (MODIS) on NASA's Terra
satellite for the years 2000-2020 (Hall et al., 2006). Soil data were
obtained from the Unified North America Soil Map which is a combination
of soil mapping from Canada and global soil maps (Liu et al., 2014).
Topographical maps were obtained from United States Geological Survey
(USGS) eros (EROS, 2010). Forest fire risk maps were obtained from
Genonode CSIRO (GFDRR, 2022). Data for each location of the IUFRO
population is made available as a supplemental table (Supplementary
table 1)

\textbf{GWAS}

GAPIT3 was used to complete the GWAS to identify adaptive loci (Wang \&
Zhang, 2021). GAPIT3 model selection was used to select the appropriate
model for performing the GWAS. Fixed and random model Circulating
Probability Unification (FarmCPU) was chosen as it was the most
effective for the size of the matrix used and is effective in the
control of false positives {[}33{]}. The GAPIT3 FarmCPU model was run
with an optimal 3 principal components to explain the variation. Regions
500bp from each loci were highlighted for future functional analysis.

\textbf{CANCOR}

The CANCOR test was based on scripts by (Blanco-Pastor et al., 2021).
The CANCOR test was used to analyse both environmental and phenotypic
data to find loci responsible for adaptations. Environmental and
phenotypic data were scaled using scale() function in R and then run
against the genotypic data using the CCorA function in the R package
vegan {[}34{]}. Then the significance of outlier loci was tested using
an \selectlanguage{greek}χ\selectlanguage{english}\textsuperscript{2} test on Mahalanobis distances. Loci were
considered an outlier if the FDR=0.1. The best represented loci in the
first two canonical dimensions were selected. To simplify results we
removed soil parameters as their projection norms weren't above 0.9 for
any soil parameter, meaning few significant loci could be found. Loci
with an alternative allele frequency that was relativly high
(\textbar{}\emph{r} \textbar{}\textgreater{}0.25) were kept. This
resulted in lists of loci that were positively or negatively correlated
with a trait of environmental parameter.

\textbf{Allele Frequency's}

To investigate the Minor Allele Frequency (MAF) of the adaptive traits
VCFtools v0.1.16 was used to separate adaptive and non-adaptive loci and
report on MAF statistics. VCFtools was also used to separate loci from
traits of interest and report on their MAF. The provenances were split
into Northern (north of 50\textsuperscript{th} latitude) and Southern
(south of 50\textsuperscript{th} latitude). MAF for each trait of
interest for was investigated for differences between the Southern and
Northern populations.

\textbf{Loci Tracing in Populations}

In order to visualise the distribution of traits associated with height
loci were traced amongst provenances. Loci associated with tree height
from the CANCOR analysis were selected and were searched for in the 80
geographic provenances using VCFtools v0.1.16. Each locus was labelled
with a positive and negative effect on height based on the CANCOR
analysis. Mean number of loci per provenance associated with a negative
or positive change in height were counted using R dplyr. Number of
positive or negative loci per provenance were geographically mapped
using the R package naturalearth (South, 2022).

\textbf{Functional Analysis}

BLASTX was used to search for genes positioned within 500bp of loci
(Altschul et al., 1990). Blast results were migrated to OmicsBox where
gene ontology was analysed (G\selectlanguage{ngerman}ötz et al., 2008).

\textbf{Co-association Analysis}

The co-association analysis was used to investigate groups of loci
associated with traits of interest. The results from the CANCOR analysis
were used as the input matrix for the co-association analysis. The input
matrix consisted of correlations of loci and environmental traits. The
Euclidian pairwise distances between the correlations of loci and
environmental traits were calculated using the dist function in R and
clustered using hierarchical clustering in the stats package in R.
Co-associations were visualised using the ComplexHeatmap package in R
(Gu et al., 2016).

\textbf{Results}

\textbf{GWAS/GEA}

The GAPIT3 model selection was performed to compare the efficacy of five
models, namely General Linear Model (GLM), Mixed Linear Model (MLM),
Compressed MLM (CMLM), Fixed and random model Circulating Probability
Unification (FarmCPU), and Bayesian-information and
Linkage-disequilibrium Iteratively Nested Keyway (BLINK) (Wang \& Zhang,
2021). FarmCPU was selected based off model selection and for its false
positive control (Liu et al., 2016). It was then implemented to test
both environmental and phenotypic parameters.

From the GEA analysis, 573 adaptive loci passed the significance
threshold of a false discovery rate of 0.1. 139 of these involved soil
parameters, 393 involved climatic conditions and 41 of these were
landscape specific. All loci and positions have been reported
(Supplementary Information 2). The GWAS discovered eight adaptive loci
for height and one for DBH. Of the loci involved with height, three were
also adaptive loci for environmental parameters, namely minimum water
vapour, forest fire risk and maximum solar radiation. As forest fire
risk assessment is based on aridity this can be considered a variable.

\textbf{CANCOR Analysis}

The CANCOR analysis revealed positive and negative adaptive loci that
were correlated (\textbar{}\emph{r} \textbar{}=0.25) with phenotypic
traits (Fig.1). Ten adaptive loci with positive phenotypic effects on
height were found but no adaptive loci for positive increases in DBH
being discovered. Incidentally, two of these loci for height were
identified in the GAPIT analyse. The two loci were associated with
environmental traits for maximum solar radiation and forest fire risk
(aridity-based parameters). Twenty-six loci were associated with a
reduction in height and one locus was associated with a reduction in
DBH. One of these loci was shared between DBH and height and none of
these loci were shared between CANCOR and GAPIT . From the environmental
analyses, 121 loci were positively adaptive towards one or more
variables. Of these 62 were shared between the GAPIT and CANCOR
analyses. 131 loci were negatively correlated with to one or more
environmental parameters. Zero of these negatively adaptive loci were
shared with the GAPIT analyses.

Association between environmental and phenotypic traits was defined as
sharing a position in the first two canonical dimensions. All positive
associations were relating to height. Four associations were due to
maximum solar radiation (max\_srad) which is typical for spruce trees to
grow taller in the presence of solar radiation, however, these loci may
indicate more efficient use of solar radiation. Two were towards and two
involved precipitation seasonality. Interestingly one association linked
height increases with forest fire risk (aridity-based parameters) which
indicated aridity adaptations. There were twenty-six loci associated
with a negative reduction in height. Of these, 15 loci were associated
with maximum or mean snow cover indicating adaptive reduction in size to
cope with high snow cover. 11 loci were involved with minimum or mean
solar radiation (mean srad).

There were approximately two environmental gradients discovered in the
environmental analyses (Fig.2). The first, and main cluster, is found
between quadrants 2 and 3 with adaptive loci being found for solar
radiation (mean srad, min\_srad, max\_srad), aridity (forest fire,
isothemality, precipitation seasonality, precipitation of driest month)
and temperature extremes (max temp, min temp) (Fig.2a). The second
cluster, which was smaller, is found between quadrants one and four with
adaptive loci for snowfall (min snow cover, mean snow cover, max snow
cover) and temperature ranges (temperature seasonality, temperature
annual range) being discovered (Fig.2a). Height is mainly affected by
the first cluster and DBH is located in quadrant three but is not
directly affected by any environmental parameters (Fig.2b).

\textbf{Allele Frequencies and Loci distribution}

To characterise adaptation, minor allele frequency (MAF) of adaptive
loci and non-adaptive loci was measured. Mean MAF for adaptive loci was
0.439 and mean MAF for non-adaptive loci was 0.117. In figure 3 A and B
we can see a considerable shift in MAF between adaptive and non-adaptive
loci indicating that environmental is characterised by high minor allele
frequency. The same cannot be concluded from phenotypic however, because
only nine adaptive loci were discovered from the phenotypic GAPIT
analysis. However, the mean MAF for phenotypic adaptive loci was 0.403.

Loci with either a negative or positive correlation with height were
traced through the population and mean number of loci per sample per
provenance were represented geographically (Fig.4). Northern populations
contained loci with negative correlations with height mainly caused by
high snow cover and low solar radiation. More negative loci were
distributed in Alaska than further south. Loci with positive
correlations towards height were distributed in the southern ranges of
Sitka spruce.

\textbf{Functional Annotation}

Of the loci discovered from the GAPIT and CANCOR analyses, 169 of them
could be functionally annotated as being located within a gene. 135 of
these had GO terms, 12 were unknown proteins and 20 were hypothetical
proteins. The loci which correlated with height or DBH and an
environmental parameter were functionally annotated using BLAST. Of
these 14 were annotated and 13 had GO terms (Fig.5).

\textbf{Co-association Analysis}

The co-association analysis revealed clusters of loci that correlated
with traits of interest. The loci with a relatively high correlation
(\textbar{}\emph{r} \textbar{}=0.25) from the CANCOR analysis could be
grouped into three main clusters. Cluster one and three are mainly
associated with traits associated with good growing conditions in Sitka.
Cluster three is strongly associated with these traits. Cluster two is
strongly associated with traits associated with freezing and snow cover.
Cluster one contains 425 loci in total, cluster 2 contains 562 loci and
cluster 3 contains 224 loci.

\textbf{Discussion}

\textbf{Adaptation in Sitka Spruce}

Studies such as GWAS are ideally suited to conifers due to the
characteristics of their genome which is marked by a slow rate of
evolution and low diversity despite the large range of ecological niches
the populations occupy (Farjon, 2010). The low rates of LD and the large
range of physical distance slow the rate of evolution and lead to clines
of locally adapted genotypes. This, however, is ideal for GEA and GWAS
studies, as high LD is problematic for the discovery of loci via GWAS,
due to the linkage of adaptive loci causing false discovery
(Christoforou et al., 2012). High LD does however necessitate very high
marker density. However linkage of adaptive loci can lead to evolution
by adaption through divergence, so harsh LD filtering is also not
recommended. Here we have discovered 694 loci that have positive effects
on the growth of spruce. We have discovered a split between the northern
and southern populations shows loci involved with increases in height in
the southern population. This is similar to what has been presented
previously (Mimura \& Aitken, 2010), with reductions in growth rate in
northern populations and slowing of bud set and burst. The allelic
distribution pattern suggests local adaptation in the northern and
southern populations. The distribution of these adaptive loci suggests
these alleles occurred in the farthest extent north and are naturally
selected until it meets the southern range in central British Columbia.
In the southern range, natural selection favours loci with positive
correlations toward height. This highlights a range ideal for breeding
programs that focus on height.

Previous studies in conifers have shown subtle shifts in MAF responsible
for adaptation (De La Torre et al., 2019). Here we see adaptive loci
with very high MAF compared to overall MAF. This indicates recent local
adaption due to the non-conserved alleles (Günther \& Coop, 2013). This
would agree with what is seen previously with post glaciation spread of
Sitka spruce creating the need for local adaptation (Byrne et al., 2022;
Gapare et al., 2005). Here we found shifts in MAF associated with some
climactic adaptive alleles. Traits identified as associated with
southern conditions (Fig.6) had subtle to moderate increases in MAF in
southern ranges indicating recent adaption. The same is seen in northern
ranges where traits associated with cold and snowfall saw subtle
increases in MAF. This subtle to moderate change in MAF is
characteristic of adaptations in conifers, as shown in \emph{Picea
glauca} and \emph{Pinus taeda} (De La Torre et al., 2019; Hornoy et al.,
2015). This indicates local adaption across a cline.

\textbf{Adaptations to Northern Wintering Conditions}

The two main wintering stressors on Sitka spruce are increased pressure
involved with snowfall and freezing. Northern conifers already have
conserved adaptions to snowfall, with homeotic genes conferring a
triangular form factor and thin needles to reduce the build-up of snow
(Du et al., 2020; Uddenberg et al., 2015). Here we identified a
reduction in size as a major adaption to snowfall. The high minor allele
frequency of these adaptive alleles indicates recent local adaptions.
This is clear in the geographic distribution of these alleles with
northern populations having more alleles negatively affecting height.
From the co-association analysis, we see the traits associated with
freezing grouping together as cluster 3. These traits are temperature
seasonality and minimum, maximum and mean snow cover. From the CANCOR
analysis we see a negative correlation between height and these freezing
traits.

\textbf{Adaptations to Southern Conditions}

The CANCOR analysis revealed a southern range of Sitka that contains
loci that correlate strongly with increases in height. From previous
studies, Sitka spruce has been shown to recolonise, post glaciation,
from Kodiak Island, Haidi Gwaii and the southern range (Byrne et al.,
2022). Here we have identified southern range specific adaptions that
allow for optimal growth in this range. This southern range has
implications for the adaptive evolution of Sitka, showing that southern
recolonisation led to adaptive loci which spread through the populations
to the 50\textsuperscript{th} latitude. The traits identified from the
co-association study show that cluster 3 can be strongly linked to good
growing conditions, with cluster 1 being mildly linked. Minimum, maximum
and mean solar radiation are all associated with loci involved with
increased height. This is common as an increase in solar radiation
allows for increased tree growth (Strand et al., 2006). Increases in the
diurnal range allow greater amounts of solar radiation, increasing
growth. The temperate climate of the southern range is optimal for
growth, with freezing temperatures causing growth cessation, as seen in
the northern range being less likely. Mean vapour is grouped with these
traits as Sitka is well adapted to receiving moisture through sea fog.
Forest fire risk is seen as linked to the same cluster of traits as it
is based on the aridity index. This may indicate some adaptations in
this loci cluster to aridity in these conditions. It may also be linked
to forest fires themselves, as known adaptations in other conifers have
been found. This would take additional phenotyping to elucidate. For a
breeding context we have defined a cluster of loci that should be
included for better growing trees in ideal conditions. We have also
identified a range of trees that could be included in breeding programs
to increase height of trees. These loci could be included through
genomic selection models by adding weight to selected loci.

The methods employed here have allowed for the merging of phenotype,
genotype and environmental data to elucidate the relationships between
the three. Here we have grouped loci into three clusters which relate to
wintering traits, like snowfall, and southern growth conditions which
seem optimal for maximising tree size. We have traced the spread of
these loci and discovered a north south divide on loci involved with
height with the southern population having a higher proportion of loci
positively correlated to height. Many of these adaptions have occurred
recently and are not conserved throughout the population, as shown by
the high MAF of adaptive loci compared to non-adaptive loci. We have
also shown that subtle shifts in allele frequency between northern and
southern populations are responsible for adaption to certain traits in
Sitka spruce. With improvement of genome resources discovered loci can
be anchored to chromosomes.

\textbf{Acknowledgements}

The authors would like to acknowledge the financial support of the Irish
Department of Agriculture, Food and the Marine under grant award
17/C/297. The authors would also like to thank Derek Gibson for help
with phenotypic measurements.

\textbf{References}

Ahuja, M. R., \& Neale, D. B. (2005). Evolution of Genome Size in
Conifers. \emph{Silvae Genetica} , \emph{54} (1-6), 126-137.
https://doi.org/doi:10.1515/sg-2005-0020

Altschul, S. F., Gish, W., Miller, W., Myers, E. W., \& Lipman, D. J.
(1990). Basic local alignment search tool. \emph{J Molecular Biology}
,\emph{215} (3), 403-410. https://doi.org/10.1016/s0022-2836(05)80360-2

Blanco-Pastor, J. L., Barre, P., Keep, T., Ledauphin, T.,
Escobar-Gutierrez, A., Roschanski, A. M., Willner, E., Dehmer, K. J.,
Hegarty, M., Muylle, H., Veeckman, E., Vandepoele, K., Ruttink, T.,
Roldan-Ruiz, I., Manel, S., \& Sampoux, J. P. (2021). Canonical
correlations reveal adaptive loci and phenotypic responses to climate in
perennial ryegrass. \emph{Molecular Ecology Resourses} , \emph{21} (3),
849-870. https://doi.org/10.1111/1755-0998.13289

Brodribb, T. J., McAdam, S. A. M., Jordan, G. J., \& Martins, S. C. V.
(2014). Conifer species adapt to low-rainfall climates by following one
of two divergent pathways. \emph{Proceedings of the National Academy of
Sciences} , \emph{111} (40), 14489-14493.
https://doi.org/doi:10.1073/pnas.1407930111

Byrne, T., Farrelly, N., Kelleher, C., Hodkinson, T., Byrne, S., \&
Barth, S. (2022). Genetic Diversity and Structure of a Diverse
Population of Picea sitchensis Using
Genotyping-by-Sequencing.\emph{Forests} , \emph{13} (9), 1511.
https://www.mdpi.com/1999-4907/13/9/1511

Casola, C. (2019). Resequencing of massive pine genomes helps to unlock
the genetic underpinning of quantitative traits in conifer
trees.\emph{New Phytologist} , \emph{221} (4), 1669-1671.
https://doi.org/10.1111/nph.15655

Chang, C. Y., Brautigam, K., Huner, N. P. A., \& Ensminger, I. (2021).
Champions of winter survival: cold acclimation and molecular regulation
of cold hardiness in evergreen conifers. \emph{New Phytologist}
,\emph{229} (2), 675-691. https://doi.org/10.1111/nph.16904

Chen, Z.-Q., Zan, Y., Milesi, P., Zhou, L., Chen, J., Li, L., Cui, B.,
Niu, S., Westin, J., Karlsson, B., García-Gil, M. R., Lascoux, M., \&
Wu, H. X. (2021). Leveraging breeding programs and genomic data in
Norway spruce (Picea abies L. Karst) for GWAS analysis. \emph{Genome
Biology} , \emph{22} (1), 179.
https://doi.org/10.1186/s13059-021-02392-1

Christoforou, A., Dondrup, M., Mattingsdal, M., Mattheisen, M.,
Giddaluru, S., Nöthen, M. M., Rietschel, M., Cichon, S., Djurovic, S.,
Andreassen, O. A., Jonassen, I., Steen, V. M., Puntervoll, P., \& Le
Hellard, S. (2012). Linkage-disequilibrium-based binning affects the
interpretation of GWASs. \emph{Am J Hum Genet} , \emph{90} (4), 727-733.
https://doi.org/10.1016/j.ajhg.2012.02.025

Danecek, P., Bonfield, J. K., Liddle, J., Marshall, J., Ohan, V.,
Pollard, M. O., Whitwham, A., Keane, T., McCarthy, S. A., Davies, R. M.,
\& Li, H. (2021). Twelve years of SAMtools and
BCFtools.\emph{Gigascience} , \emph{10} (2).
https://doi.org/10.1093/gigascience/giab008

De La Torre, A. R., Wilhite, B., \& Neale, D. B. (2019). Environmental
Genome-Wide Association Reveals Climate Adaptation Is Shaped by Subtle
to Moderate Allele Frequency Shifts in Loblolly Pine. \emph{Genome
Biology Evolution} , \emph{11} (10), 2976-2989.
https://doi.org/10.1093/gbe/evz220

Du, H., Ran, J.-H., Feng, Y.-Y., \& Wang, X.-Q. (2020). The flattened
and needlelike leaves of the pine family (Pinaceae) share a conserved
genetic network for adaxial-abaxial polarity but have diverged for
photosynthetic adaptation. \emph{BMC Evolutionary Biology} ,\emph{20}
(1), 131. https://doi.org/10.1186/s12862-020-01694-5

Elshire, R. J., Glaubitz, J. C., Sun, Q., Poland, J. A., Kawamoto, K.,
Buckler, E. S., \& Mitchell, S. E. (2011). A robust, simple
genotyping-by-sequencing (GBS) approach for high diversity
species.\emph{PLoS One} , \emph{6} (5), e19379.
https://doi.org/10.1371/journal.pone.0019379

Farjon, A. (2010). \emph{A Handbook of the World's Conifers} . Brill.
https://doi.org/https://doi.org/10.1163/9789047430629

Fick, S. E., \& Hijmans, R. J. (2017). WorldClim 2: new 1-km spatial
resolution climate surfaces for global land areas. \emph{International
Journal of Climatology} , \emph{37} (12), 4302-4315.
https://doi.org/https://doi.org/10.1002/joc.5086

Gagalova, K. K., Warren, R. L., Coombe, L., Wong, J., Nip, K. M., Yuen,
M. M. S., Whitehill, J. G. A., Celedon, J. M., Ritland, C., Taylor, G.
A., Cheng, D., Plettner, P., Hammond, S. A., Mohamadi, H., Zhao, Y.,
Moore, R. A., Mungall, A. J., Boyle, B., Laroche, J., . . . Birol, I.
(2022). Spruce giga-genomes: structurally similar yet distinctive with
differentially expanding gene families and rapidly evolving
genes.\emph{Plant J} , \emph{111} (5), 1469-1485.
https://doi.org/10.1111/tpj.15889

Gapare, W. J., Aitken, S. N., \& Ritland, C. E. (2005). Genetic
diversity of core and peripheral Sitka spruce (Picea sitchensis (Bong.)
Carr) populations: implications for conservation of widespread
species.\emph{Biological Conservation} , \emph{123} (1), 113-123.
https://doi.org/10.1016/j.biocon.2004.11.002

GFDRR. (2022). \emph{WF-GLOBAL-CSIRO-30} .
https://www.geonode-gfdrrlab.org/layers/hazard:csiro\_wf\_max\_fwi\_rp30\#more

Gu, Z., Eils, R., \& Schlesner, M. (2016). Complex heatmaps reveal
patterns and correlations in multidimensional genomic
data.\emph{Bioinformatics} , \emph{32} (18), 2847-2849.
https://doi.org/10.1093/bioinformatics/btw313

Günther, T., \& Coop, G. (2013). Robust identification of local
adaptation from allele frequencies. \emph{Genetics} , \emph{195} (1),
205-220. https://doi.org/10.1534/genetics.113.152462

Hall, D. K., Riggs, G. A., \& Salomonson, V. V. (2006).\emph{MODIS/Terra
Snow Cover 5-Min L2 Swath 500m.} Version 5).
https://doi.org/http://dx.doi.org/10.5067/ACYTYZB9BEOS.

Hiraoka, Y., Fukatsu, E., Mishima, K., Hirao, T., Teshima, K. M.,
Tamura, M., Tsubomura, M., Iki, T., Kurita, M., Takahashi, M., \&
Watanabe, A. (2018). Potential of Genome-Wide Studies in Unrelated Plus
Trees of a Coniferous Species, Cryptomeria japonica (Japanese
Cedar).\emph{Frontiers Plant Science} , \emph{9} , 1322.
https://doi.org/10.3389/fpls.2018.01322

Hornoy, B., Pavy, N., Gerardi, S., Beaulieu, J., \& Bousquet, J. (2015).
Genetic Adaptation to Climate in White Spruce Involves Small to Moderate
Allele Frequency Shifts in Functionally Diverse Genes. \emph{Genome
Biology Evolution} , \emph{7} (12), 3269-3285.
https://doi.org/10.1093/gbe/evv218

Jurgiel, B. (2020). \emph{Point sampling tool} . Retrieved 07/09/2022
from https://github.com/borysiasty/pointsamplingtool

Li, H., \& Durbin, R. (2009). Fast and accurate short read alignment
with Burrows-Wheeler transform. \emph{Bioinformatics} , \emph{25} (14),
1754-1760. https://doi.org/10.1093/bioinformatics/btp324

Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N.,
Marth, G., Abecasis, G., Durbin, R., \& Genome Project Data Processing,
S. (2009). The Sequence Alignment/Map format and
SAMtools.\emph{Bioinformatics} , \emph{25} (16), 2078-2079.
https://doi.org/10.1093/bioinformatics/btp352

Liu, S., Y. Wei, W.M. Post, R.B. Cook, K. Schaefer, \& Thornton., M. M.
(2014). \emph{NACP MsTMIP: Unified North American Soil Map} .
https://doi.org/https://doi.org/10.3334/ORNLDAAC/1242

Liu, X., Huang, M., Fan, B., Buckler, E. S., \& Zhang, Z. (2016).
Iterative Usage of Fixed and Random Effect Models for Powerful and
Efficient Genome-Wide Association Studies. \emph{PLoS Genetics}
,\emph{12} (2), e1005767. https://doi.org/10.1371/journal.pgen.1005767

Lotterhos, K. E., Yeaman, S., Degner, J., Aitken, S., \& Hodgins, K. A.
(2018). Modularity of genes involved in local adaptation to climate
despite physical linkage. \emph{Genome Biology} , \emph{19} (1), 157.
https://doi.org/10.1186/s13059-018-1545-7

Mimura, M., \& Aitken, S. N. (2010). Local adaptation at the range
peripheries of Sitka spruce. \emph{J Evolutionary Biology} ,\emph{23}
(2), 249-258. https://doi.org/10.1111/j.1420-9101.2009.01910.x

Pavy, N., Namroud, M. C., Gagnon, F., Isabel, N., \& Bousquet, J.
(2012). The heterogeneous levels of linkage disequilibrium in white
spruce genes and comparative analysis with other conifers.\emph{Heredity
(Edinb)} , \emph{108} (3), 273-284. https://doi.org/10.1038/hdy.2011.72

Prunier, J., Verta, J. P., \& MacKay, J. J. (2016). Conifer genomics and
adaptation: at the crossroads of genetic diversity and genome
function.\emph{New Phytologist} , \emph{209} (1), 44-62.
https://doi.org/10.1111/nph.13565

QGIS.org. (2022). \emph{QGIS Geographic Information System} . In QGIS
Association. http://www.qgis.org

Sebastian-Azcona, J., Hacke, U. G., \& Hamann, A. (2018). Adaptations of
white spruce to climate: strong intraspecific differences in cold
hardiness linked to survival. \emph{Ecology Evolution} , \emph{8} (3),
1758-1768. https://doi.org/10.1002/ece3.3796

Senser, M., \& Beck, E. (1984). Correlation of chloroplast
ultrastructure and membrane lipid composition to the different degrees
of frost resistance achieved in leaves of spinach, ivy, and
spruce.\emph{J Plant Physiology} , \emph{117} (1), 41-55.
https://doi.org/10.1016/S0176-1617(84)80015-2

South, S. (2022). \emph{rnaturalearth: World Map Data from Natural
Earth.} . https://github.com/ropensci/rnaturalearth.

Strand, M., Löfvenius, M. O., Bergsten, U., Lundmark, T., \& Rosvall, O.
(2006). Height growth of planted conifer seedlings in relation to solar
radiation and position in Scots pine shelterwood. \emph{Forest Ecology
and Management} , \emph{224} (3), 258-265.
https://doi.org/https://doi.org/10.1016/j.foreco.2005.12.038

Sype, H. R.-A., B. (1990). \emph{Genetic variability of Sitka spruce of
the IUFRO collection} Quebec, Canada.

Teskey, R., Wertin, T., Bauweraerts, I., Ameye, M., McGuire, M. A., \&
Steppe, K. (2015). Responses of tree species to heat waves and extreme
heat events. \emph{Plant Cell Environ} , \emph{38} (9), 1699-1712.
https://doi.org/10.1111/pce.12417

Uddenberg, D., Akhter, S., Ramachandran, P., Sundström, J. F., \&
Carlsbecker, A. (2015). Sequenced genomes and rapidly emerging
technologies pave the way for conifer evolutionary developmental
biology. \emph{Frontiers Plant Science} , \emph{6} , 970.
https://doi.org/10.3389/fpls.2015.00970

Wang, J., \& Zhang, Z. (2021). GAPIT Version 3: Boosting Power and
Accuracy for Genomic Association and Prediction. \emph{Genomics,
Proteomics \& Bioinformatics} , \emph{19} (4), 629-640.
https://doi.org/https://doi.org/10.1016/j.gpb.2021.08.005

\textbf{Data Accessibility and Benefit-Sharing}

Data accessibility: Raw sequence data were submitted to NCBI (BioProject
PRJNA852515). The phenotypic data are accessible in the supporting
information submitted with this manuscript from supplementary Excel file
`Byrne et al\_Supplementary data'. This file contains detailed
geographical information on 79 provenances including meteorological data
and soil data used in the analysis of this work.

Benefit sharing: The research described in this work complies with Irish
national laws implementing the Convention on Biological Diversity and
Nagoya Protocol agreements.

\textbf{Author contributions}

Conceptualization, SLB, NF, TRH, CK and SB; methodology, SLB, NF, TRH,
SB, CK and TB; software, TB and SLB; validation, SLB, TRH, NF, CK and
SB; formal analysis, TB and SLB; investigation, TB; data curation, TB
and SLB; writing---original draft preparation, TB; writing---review and
editing, TB, SLB, TRH, NF, CK and SB; visualization, TB; supervision,
SLB, TRH, NF, CK and SB; project administration, NF; funding
acquisition, SLB, TRH, NF, CK and SB. All authors have read and agreed
to the published version of the manuscript.

\textbf{Tables and Figures}

\textbf{Table 1. List of Environmental traits used in GAPIT and CANCOR
analysis}\selectlanguage{english}
\begin{longtable}[]{@{}ll@{}}
\toprule
\emph{Landscape Features} &\tabularnewline
Latitude & Decimal\tabularnewline
Longitude & Decimal\tabularnewline
Altitude & Meter\tabularnewline
Aspect & Degrees\tabularnewline
Forest fire risk & 0-300 (300 is high risk)\tabularnewline
Distance to coast & Kilometres\tabularnewline
\emph{Climate} &\tabularnewline
Precipitation of coldest quarter & mm\tabularnewline
Precipitation of warmest quarter & mm\tabularnewline
Precipitation of driest quarter & mm\tabularnewline
Precipitation of wettest quarter & mm\tabularnewline
Precipitation seasonality (coefficient of variation) & mm\tabularnewline
Precipitation of driest month & mm\tabularnewline
Precipitation of wettest month & mm\tabularnewline
Annual precipitation & mm\tabularnewline
Mean temperature of coldest quarter & \selectlanguage{ngerman}°C\tabularnewline
Mean temperature of warmest quarter & \selectlanguage{ngerman}°C\tabularnewline
Mean temperature of driest quarter & \selectlanguage{ngerman}°C\tabularnewline
Mean temperature of wettest quarter & \selectlanguage{ngerman}°C\tabularnewline
Temperature annual range & \selectlanguage{ngerman}°C\tabularnewline
Minimum temperature & \selectlanguage{ngerman}°C\tabularnewline
Maximum temperature & \selectlanguage{ngerman}°C\tabularnewline
Temperature seasonality & \selectlanguage{ngerman}°C\tabularnewline
Isothermality & \selectlanguage{ngerman}°C\tabularnewline
Mean diurnal range & \selectlanguage{ngerman}°C\tabularnewline
Annual mean temperature & \selectlanguage{ngerman}°C\tabularnewline
Mean wind speed & m/s\tabularnewline
Min wind speed & m/s\tabularnewline
Max wind speed & m/s\tabularnewline
Mean solar radiation & kJ/m/day\tabularnewline
Minimum solar radiation & kJ/m/day\tabularnewline
Maximum solar radiation & kJ/m/day\tabularnewline
Mean vapour & kPa\tabularnewline
Minimum vapour & kPa\tabularnewline
Max vapour & kPa\tabularnewline
Mean snow cover & \%\tabularnewline
Minimum snow cover & \%\tabularnewline
Maximum snow cover & \%\tabularnewline
\emph{Soil} &\tabularnewline
Dominant component & \%\tabularnewline
Soil depth & cm\tabularnewline
Top soil sand & \% weight\tabularnewline
Top soil silt & \% weight\tabularnewline
Top soil clay & \% weight\tabularnewline
Top soil gravel & \% volume\tabularnewline
Top soil organic carbon & \% weight\tabularnewline
Top soil pH & scale\tabularnewline
Top soil bulk density & g cm-3\tabularnewline
Top soil cation & meq/100g\tabularnewline
Subsoil sand & \% weight\tabularnewline
Subsoil silt & \% weight\tabularnewline
Subsoil clay & \% weight\tabularnewline
Subsoil gravel & \% volume\tabularnewline
Subsoil organic carbon & \% weight\tabularnewline
Subsoil pH & scale\tabularnewline
Subsoil bulk density & g cm-3\tabularnewline
Subsoil cation & meq/100g\tabularnewline
\bottomrule
\end{longtable}

\textbf{Figures}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/image1}

\end{center}
\end{figure}

Figure 1. Outlier loci for specific environmental and phenotypic data.
The outlier status of loci is determined by the Canonical Correlation
(CANCOR) analysis where minor allele frequency is relatively highly
correlated (\textbar{}\emph{r} \textbar{}=0.25) with a trait. Loci
coloured grey represent all the loci tested. Loci coloured red are
considered negatively correlated and loci coloured green are considered
positively correlated.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image2/image2}

\end{center}
\end{figure}

Figure 2. (A) Representation of environmental variables on canonical
axis as correlated with alternate allele frequency. (B) Representation
of phenotypic variables on canonical axis as correlated with alternate
allele frequency. (C) Representation of loci on the first two canonical
axis. In (A-B) the inner circle marks 0.9 on projection norm and the
outer circle represents 1 on projection norm values.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image3/image3}

\end{center}
\end{figure}

Figure 3. (A-B) Minor allele frequency of adaptive alleles and
non-adaptive alleles based on GAPIT3 analysis. (C) Shifts in minor
allele frequency of adaptive alleles between Northern and Southern
ranges of the population based of GAPIT analysis.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image4/image4}

\end{center}
\end{figure}

Figure 4. Geographic distribution of adaptive loci positively and
negatively correlated with height across the Pacific Northwest. Red
circles refer to loci negatively correlated with height whereas green
circles refer to positive correlations. Circle size refers to number of
loci associated with height.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image5/image5}

\end{center}
\end{figure}

Figure 5. Gene ontology of regions associated with adaptive loci. Gene
ontology terms shown have a frequency of greater or equal to three.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image6/image6}

\end{center}
\end{figure}

Figure 6. Environmental co-association analysis from the results of
CANCOR analysis. Correlations between SNP and environmental factors were
used for clustering analyses and pairwise comparisons.

\textbf{Supplementary data}

The following supporting information can be downloaded in this
supplementary Excel file `Byrne et al\_Supplementary data'. This file
contains detailed geographical information on 79 provenances including
meteorological data and soil data used in the analysis of this work.

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

