Genetic Dissection of Growth and Eco-Physiological Traits Associated with Altitudinal Adaptation in Sakhalin Fir (Abies sachalinensis) Based on QTL Mapping

(1) Background: The genetic basis of local adaptation in conifers remains poorly understood because of limited research evidence and the lack of suitable genetic materials. Sakhalin fir (Abies sachalinensis) is an ideal organism for elucidating the genetic basis of local adaptation because its altitudinal adaptation has been demonstrated, and suitable materials for its linkage mapping are available. (2) Method: We constructed P336 and P236 linkage maps based on 486 and 516 single nucleotide polymorphisms, respectively, that were derived from double digest restriction site-associated DNA sequences. We measured the growth and eco-physiological traits associated with morphology, phenology, and photosynthesis, which are considered important drivers of altitudinal adaptation. (3) Results: The quantitative trait loci (QTLs) for growth traits, phenology, needle morphology, and photosynthetic traits were subsequently detected. Similar to previous studies on conifers, most traits were controlled by multiple QTLs with small or moderate effects. Notably, we detected that one QTL for the crown area might be a type-A response regulator, a nuclear protein responsible for the cytokinin-induced shoot elongation. (4) Conclusion: The QTLs detected in this study include potentially important genomic regions linked to altitudinal adaptation in Sakhalin fir.


Introduction
Local adaptation is an evolutionary process that is considered a home-site advantage, in which a genotype performs best in its native site [1]. As sessile organisms, plants are directly affected by the local climate, which is one of the most important drivers of adaptation [2]. However, a pairwise comparison between local and foreign plants revealed that local adaptation was observed in only 45% of the 1032 studied populations [3]. Hence, local adaptation is determined by the balance between gene flow and selection intensity [4]. For example, forest tree species with extensive gene flow generally exhibit local adaptation at a wide scale [5]. Selection drivers are usually complex because of a Table 1. The Sakhalin fir segregated population used for QTL analysis. Grandmother parent Clones (A33 and A39) are mother trees inhabiting high-altitude zones (1200 m). Clone 1-1 and Clone 1-2 are low-altitude genotypes derived from a pollen mixture collected from three genotypes living at low-altitudes (530 m). The saplings were grown for another two consecutive growing seasons at the Forestry Research Institute, Hokkaido Research Organization (HRO), Bibai, Japan (43 • 17 N, 141 • 51 E, 40 m asl). Considering the environmental heterogeneity among microhabitats (e.g., sun exposure), periodic rotation of the pots was performed. Additionally, the saplings were irrigated every morning to prevent drought stress. The soil was watered with Osmocote ® Exact Standard 3-4M (N16-P9-K12) (ICL Specialty Fertilizers, Geldermalsen, The Netherlands) at the start of each growing season and with HYPONeX liquid (N6-P3360-K5) (HYPONeX Japan Co., Ltd., Osaka, Japan) during the latter part of the growing season, following the manufacturers' protocols.

Pedigree
The saplings were subsequently transplanted to a common garden that was established at the UTHF (43 • 13 N, 142 • 24 E, 230 m asl) on 27 April 2017. The saplings were planted in a randomized order at 2.0 m intervals. No irrigation or fertilization was performed for the transplants in the common garden, only weed control.

Phenotyping
The phenotype data of 15 functional traits, including four growth and 11 eco-physiological traits, for the 252 saplings in the segregation population were recorded (Table 2). On 4 October 2016, the stem diameter (D16, mm) and height (H16, cm) of the saplings at the HRO were measured. After transplanting to the UTHF, the tree heights (H17, cm) were measured in the autumn of 2017, and images via photograph-based projection were captured from the top points of each tree. Then, the crown area (CR17, cm 2 ) was estimated using LIA32 software (http://www.agr.nagoya-u.ac.jp/~shinkan/LIA32/) ( Figure 1).

Characterization of Morphological Traits
One-year-old needles were randomly collected from the middle portion of each branch and scanned using a digital scanner (MSC20; King Jim, Tokyo, Japan). The mean length and width of the needles (Lw_ratio) were determined using ImageJ software (NIH, Bethesda, MD, USA; http://rsbweb.nih.gov/ij/). The leaf mass per area (LMA, g cm −2 ) was calculated by dividing the dry mass of 20 needles by the total needle area using ImageJ software. Following a previously used method [17], the stomatal density (SD, number mm −2 ) was estimated as the ratio of the number of stomata to the respective stomatal band area and the number of stomatal rows (SRN, number mm −1 ) was measured from the same images in a 0.88 mm × 0.66 mm area. Transverse sections of the current season shoots of the saplings were photographed using a digital camera (DP71; Olympus, Tokyo, Japan) mounted on a light microscope (BX50; Olympus). The images were used to measure the ratio of bark width to xylem width (Bark_xy) and the ratio of normal wood width to reaction wood width (Norm_reac) using ImageJ software.

Characterization of Morphological Traits
One-year-old needles were randomly collected from the middle portion of each branch and scanned using a digital scanner (MSC20; King Jim, Tokyo, Japan). The mean length and width of the needles (Lw_ratio) were determined using ImageJ software (NIH, Bethesda, MD, USA; http://rsbwe b.nih.gov/ij/). The leaf mass per area (LMA, g cm −2 ) was calculated by dividing the dry mass of 20 needles by the total needle area using ImageJ software. Following a previously used method [17], the stomatal density (SD, number mm −2 ) was estimated as the ratio of the number of stomata to the respective stomatal band area and the number of stomatal rows (SRN, number mm −1 ) was measured from the same images in a 0.88 mm × 0.66 mm area. Transverse sections of the current season shoots of the saplings were photographed using a digital camera (DP71; Olympus, Tokyo, Japan) mounted on a light microscope (BX50; Olympus). The images were used to measure the ratio of bark width to xylem width (Bark_xy) and the ratio of normal wood width to reaction wood width (Norm_reac) using ImageJ software.

Evaluation of Bud Phenology
As a trait associated with plant phenology, the timing of bud flush (Bud_fl) was also recorded. On 6-29 May 2016, the status of the terminal bud in each pot sapling was monitored every two days. Following a previous study [19], Bud_fl was defined as the time point at which the initial emergence of needles from the terminal bud was observed. Hence, the number of days after 1 May was recorded as the Bud_fl.

Evaluation of Freezing Tolerance
To measure freezing tolerance, a freezing test was conducted on November 4, following the method of a previous study on Sakhalin fir [25]. Before the test, three healthy current-year needles were collected in the morning, placed in a thin plastic bag (0.03 mm), and incubated for >3 h in a dark chamber maintained at 4 °C for dark acclimation. Each test was conducted by exposing the needles to the target temperature (−20 °C) for 5 h, and then decreasing temperature at a rate of 12 °C h −1 using a programmed freezer (SC-DF25K; Nippon Freezer Co., Tokyo, Japan). We measured the freezing damage to the needles after incubation for 2 days at 4 °C, followed by an increase in temperature at a rate of 8 °C h −1 using a pulse-amplitude modulated (PAM) chlorophyll fluorometer (Mini-PAM; Walz,

Evaluation of Bud Phenology
As a trait associated with plant phenology, the timing of bud flush (Bud_fl) was also recorded. On 6-29 May 2016, the status of the terminal bud in each pot sapling was monitored every two days. Following a previous study [19], Bud_fl was defined as the time point at which the initial emergence of needles from the terminal bud was observed. Hence, the number of days after 1 May was recorded as the Bud_fl.

Evaluation of Freezing Tolerance
To measure freezing tolerance, a freezing test was conducted on November 4, following the method of a previous study on Sakhalin fir [25]. Before the test, three healthy current-year needles were collected in the morning, placed in a thin plastic bag (0.03 mm), and incubated for >3 h in a dark chamber maintained at 4 • C for dark acclimation. Each test was conducted by exposing the needles to the target temperature (−20 • C) for 5 h, and then decreasing temperature at a rate of 12 • C h −1 using a programmed freezer (SC-DF25K; Nippon Freezer Co., Tokyo, Japan). We measured the freezing damage to the needles after incubation for 2 days at 4 • C, followed by an increase in temperature at a rate of 8 • C h −1 using a pulse-amplitude modulated (PAM) chlorophyll fluorometer (Mini-PAM; Walz, Effeltrich, Germany). The maximum photochemical quantum yield of photosystem II (PSII; Fv/Fm) is negatively correlated with the extent of freezing damage, and thus, represents a reliable parameter for the assessment of freezing tolerance during autumn [25,26]. The Fv/Fm values used as an index of freezing tolerance in this study were termed "Freez_tol".

Measurement of Chlorophyll Fluorescence
The saplings at the HRO were placed in a darkroom for 3 days. Chlorophyll fluorescence was then measured with a chlorophyll fluorometer (Junior-PAM; Walz, Effeltrich, Germany) using a detached leaf from the major axis of a current-season branch. The measurements were performed in a ventilated room, with 40 Pa CO 2 and 21 kPa O 2 at 5 • C. Saturation pulses from blue light-emitting diodes (LEDs; >6000 µmol photons m −2 s −1 , 800 ms duration) were applied to determine the maximum chlorophyll fluorescence, with closed PSII centers in the dark (Fm) and under actinic light (Fm ). The maximum photochemical quantum yield (Fv/Fm) and effective quantum yield (Φ II ) of PSII were calculated as (Fm-F0)/Fm and (Fm -Fs )/Fm , respectively [27], where F0 is the minimal chlorophyll fluorescence yield in the dark and Fs is the steady-state chlorophyll fluorescence level under actinic light from blue LEDs, with a wavelength peak at 445 nm. Non-photochemical quenching (NPQ) was calculated as (Fm-Fm )/Fm [28]. Additionally, the PSII quantum yields Φ NPQ and Φ NO [29] that represent the regulated and non-regulated energy dissipation at the PSII centers, respectively, and add up to the unity of the photochemical quantum yield (i.e., Φ II + Φ NPQ + Φ NO = 1), were also calculated. The values of Φ NPQ and Φ NO were calculated as Fs /Fm-Fs /Fm and Fs/Fm, respectively [30,31]. The Fv/Fm in dark-adapted leaves was measured overnight, and then the other parameters were measured under actinic light of 1500 µmol photons m −2 s −1 for 2 min.

Determination of Relationships between Phenotypic Traits
The relationships between the functional traits were analyzed using R version 3.6.1 [32]. The Pearson's correlation coefficients and p values for correlations between functional traits were calculated. To simplify and visualize the strong correlations, a network plot was generated using the R package 'corrr' version 0.4.3 [33]. The proximity of the points was determined using multidimensional clustering, and the minimum coefficient value (in absolute terms) was set to 0.1.

Construction of Double-Digest Restriction Site-Associated DNA Sequencing (ddRAD-seq) Library
A total of 252 individuals were included in ddRAD-seq experiment [34]. Compared with the original restriction site-associated DNA sequencing (RAD-seq) [35], ddRAD-seq effectively reduced the complexity of studying the large genome of conifers [36]. The total genomic DNA (250 ng) from each sample was digested with SphI and PstI, ligated with Y-shaped adaptors, and amplified by PCR using KAPA HiFi polymerase (KAPA Biosystems, Woburn, MA, USA). After PCR amplification with adapter-specific primer pairs (Access Array Barcode Library for Illumina Sequencers; Fluidigm, San Francisco, CA, USA), an equal amount of DNA from each sample was mixed and size-selected using BluePippin agarose gel cassettes (Sage Science, Beverly, MA, USA). The library fragments (~450 bp) were retrieved, and the quality of the library was checked using an Agilent 2100 Bioanalyzer with a high-sensitivity DNA chip (Agilent Technologies, Waldbronn, Germany). The library was sequenced using the Illumina ® HiSeq X platform (Illumina, San Diego, CA, USA) to generate 150 bp long paired-end reads (see details in Supplementary Materials S1). The raw ddRAD-seq data were deposited in the DNA Data Bank of Japan (DDBJ) under accession number DRA012397.

Genotyping
Quality control, adapter trimming, and quality filtering of the ddRAD-seq raw data were performed using fastp program version 0.20.0 [37]. A sliding window of 15 bases was used to filter reads with low-quality scores (Phred quality values < 10). The high-quality ddRAD-seq data were processed using the ustacks, cstacks, sstacks, tsv2bam, and gstacks programs of Stacks software version 2.5 [38], following the Stacks manual for a de novo genetic mapping cross (http://catchenlab.life.illinois.edu/stacks/manual/, accessed on 21 July 2021). This method builds catalogs from the parents of a mapping population only, which enables the exclusion of potential errors derived from the mapping progeny. The number of mismatches allowed between stacks within and between individuals was set to four (i.e., M and n options in ustacks and cstacks, respectively). The populations program in Stacks pipeline was used to retain the loci that were observed in >80% of individuals in the mapping progeny. The retained loci were exported in JoinMap format via a cross pollination (CP) mode using the same program to construct the linkage maps.

Linkage Map Construction
Two types of linkage maps (P336 and P236) were constructed based on the segregation patterns of the genetic markers, following the method of the previous study [19]. The P336 and P236 maps were comprised of the aa × ab and ab × aa type markers, respectively. Markers that contained a large amount of missing data (>20%) and showed significant deviation (chi-square test; p < 0.001) from the expected Mendelian segregation ratio were excluded. Linkage analysis was conducted for the populations of cross pollinators (CP) with the regression mapping method using JoinMap v5.0 software [39]. Default settings were used for determining the recombination frequency.

QTL Analysis
QTL analysis was conducted using linear models with a regularized horseshoe prior to distribution, which can detect and assess relevant variables from a large number of irrelevant variables [40]. First, to reduce collinearity and redundancy among variables, one marker from a given pair of markers that were in close proximity (<1 cM) and showed high correlation (r > 0.7) with other markers was removed prior to the analysis. Second, four chains were used to estimate the posterior distribution, and a total of 5000 iterations after the first 1000 iterations were discarded as burn-in. The convergence of the chains was confirmed using the Gelman-Rubin statistic (R < 1.1). Third, the posterior distribution was estimated using the R package 'rstanarm' [41]. The global shrinkage parameter of the regularized horseshoe prior was set using the equation proposed by a previous study [40]: where p0 is the expected number of relevant variables, D is the number of variables, and n is the number of observations. For all traits, the value of p0 was set to five because the QTLs of various tree species were known to exhibit small effects [42]. The sensitivity of the p0 value for QTL detection was confirmed by arranging p0 from 1 to 9 (Table S1). The default values of the regularized horseshoe prior in the R package 'rstanarm' were used for the remaining parameters.
Fourth, given that the variables in the aforementioned model were highly correlated, the estimated marginals of the posterior distribution of the marker effects may be biased because of collinearity. To address this problem, a model selection method, projectionpredictive variable selection [43], was applied to the model. The model selection method uses a model that has the best predictive power (a reference model) to find a simpler model that has similar predictions as the reference model (predictive projection). The linear model described above was used as the reference model. Model selection was performed when the variables in the reference model contained significantly more information than the null model, which was confirmed by leave-one-out cross-validation [44]. When the marginal posterior distribution of the marker effects in the selected model did not contain zero in the 95% and 80% credible intervals, the markers were defined as significant and suggestive of QTLs, respectively. The contribution of each QTL was estimated by the coefficient of determination (R 2 ) of the simple regression; and the contribution was defined as the percentage variance explained (PVE) for a QTL. The QTL analysis and model selection were both performed using R version 3.6.1 [32]. The leave-one-out cross-validation procedure was performed using the R package 'loo' [45], while the model selection method was conducted using the R package 'projpred' [43].

Candidate Gene Prediction
Genes in the target QTL regions were predicted in the Sakhalin fir transcriptome database TodoFirGene [20] using the basic local alignment search tool (BLAST) nucleotide algorithm. The transcripts of the candidate genes were obtained from TodoFirGene and subsequently annotated using the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG; https://www.kegg.jp/kegg/kegg2.html, accessed on 16 July 2021) databases.

Phenotypic Traits
Based on the results of previous studies, these traits are potentially associated with altitudinal adaptation ( Table 1). The mean values, standard deviations, and coefficients of variation for the 15 functional traits are presented in Table 1. The distribution of each phenotype was observed to have a modal shape ( Figure S1). The network analysis revealed the relationships among the functional traits investigated, in which two clusters were identified for the growth and photosynthetic traits ( Figure 2). Furthermore, growth-related traits were discovered to be positively correlated with Bud_fl, whereas Lw_ratio was negatively correlated with Bark_xy ratio. Pairwise correlations between functional traits using Pearson's correlation coefficients showed significant negative correlations between certain photosynthetic traits, such as Φ II vs. NPQ (r = −0.63, p < 0.001) and NPQ vs. Φ NO (r = −0.89, p < 0.001) ( Figure S2). However, a positive correlation was observed between Φ II and Φ NO (r = 0.24), whereas strong and highly positive correlations were observed between growth traits (r > 0.5).

Linkage Map Construction and QTL Detection
A total of 486 and 516 markers were mapped to 12 linkage groups (LG) in the P336 and P236 maps, respectively. For the P336 and P236 maps, the total lengths were 1986.2 cM and 1932.8 cM, respectively, whereas the average distances between adjacent markers were 4.1 cM and 3.7 cM, respectively. The maximum marker gaps for P336 map and P236 map were 29.8 cM and 27.3 cM, respectively (Table 3).

Linkage Map Construction and QTL Detection
A total of 486 and 516 markers were mapped to 12 linkage groups (LG) in the P336 and P236 maps, respectively. For the P336 and P236 maps, the total lengths were 1986.2 cM and 1932.8 cM, respectively, whereas the average distances between adjacent markers were 4.1 cM and 3.7 cM, respectively. The maximum marker gaps for P336 map and P236 map were 29.8 cM and 27.3 cM, respectively (Table 3).  Three significant and 11 suggestive QTLs were detected in the P336 map, whereas five significant and 10 suggestive QTLs were observed in the P236 map (Table 4; Figure 3). Significant QTLs for Φ NO and CR17 and for Φ NO , H16, H17, and CR17 were detected in the P336 and P236 maps, respectively. In addition, suggestive QTLs for NPQ, Bud_fl, Freez_tol, Lw_ratio, H16, H17, and CR17 and for NPQ, Bud_fl, Lw_ratio, D16, H17, and CR17 were detected in the P336 and P236 maps, respectively. A significant QTL (Locus #2055) was detected for NPQ and Φ NO in both linkage maps. In the P236 map, Locus #12865 was discovered as a significant QTL for CR17 and a suggestive QTL for D16, H16, and H17 ( Figure 3). In the P336 map, suggestive QTLs for Bud_fl and Freez_tol were co-located in LG5 within a 10 cM distance (Figure 3).

Candidate Gene Prediction
Out of 16 unique QTLs, eight were found in the transcripts of TodoFirGene (Table S2). Among the eight QTLs, locus #1970, which is associated with CR17, was annotated as a type-A Arabidopsis response regulator (type-A ARR) of phosphorelay signal transduction system and a type-A ARR of cytokinin signaling in the GO and KEGG databases, respectively. The sequence of locus #1970 was shown in Supplementary Materials S2. Genes 2021, 12, x FOR PEER REVIEW 10 of 16

Segregated Population and Linkage Maps
Long generation times and the outbred mating system of conifers were used to construct a suitable mapping population. Because several pedigrees are available for important traits, including height growth and wood properties, in a tree breeding program, secondgeneration populations were used for QTL mapping [46][47][48][49]. However, F 1 pedigrees have been used for QTL mapping of eco-physiological studies in the evolutionary biology of conifers [50,51]. In the present study, we used segregated populations based on control crosses with four different grandparent trees with distinct genetic backgrounds ( Table 1) [19]. Furthermore, altitudinal adaptation was clearly demonstrated by a reciprocal transplant experiment [18]. Our studies are limited by the number of pedigrees, locations and sample replications. Nevertheless, the segregated populations were constructed by a control cross between F 1 hybrids with two contrasting ecotypes (high-altitude and lowaltitude grandparents), which might be useful to clarify the genetic basis of altitudinal adaptation in conifers as described in a previous paper [19].
The linkage maps were saturated, with a 10 cM interval between nearest markers. The average distance between nearest markers detected in this study (approximately 4 cM; Table 3) was also comparable with that in previous studies [46,48,50,52]. The maximum marker gap distance in this study (approximately 20 cM) was also similar to that in a previous study [51]. Recently, high-density linkage maps with an average distance between nearest markers of less than 1 cM have been constructed with several thousand SNP markers [47,53]. In addition, gene-based markers have been applied to obtain functional information [50]. The reconstruction of the linkage map using gene-based markers might improve the accuracy of QTL mapping in future studies.

Growth Traits
The segregation populations reflected the broad growth variations along an altitudinal gradient. The coefficient of variation (CV) for the heights of 4-year-old Sakhalin fir seedlings was~0.20, which is similar to the CV value (0.19) reported for 4-year-old seedlings in the provenance trials [16]. The growth traits analyzed in the present study were mostly controlled by several QTLs with small or moderate effects (PVE: 2.6-9.3%) ( Figure 3; Table 4). This result corroborates that of a previous study on Sakhalin fir, in which several QTLs with small or moderate effects (PVE: 3.8-7.4%) were detected for H17 and D16 at the early stage. QTLs with small or moderate effects on growth traits have also been detected in other conifers [46][47][48][49].
The crown area and width seem to be a good indicator for growth. In the present study, CR17 was measured using the digital images for growth traits (Figure 1). The number of QTLs for CR17 was the highest among the traits examined in this study; five exhibited moderate effects (PVE > 5%; Table 2; Figure 3). The QTLs for crown width detected in this study have also been reported in other conifers [47]. Notably, Locus #1970 of the QTLs for CR17 was annotated as a type-A Arabidopsis Response Regulator (type-A ARR; Table S2), which functions as a key node in the integration of ethylene and cytokinin signaling to regulate plant responses and shoot elongation in Arabidopsis thaliana [54]. Type-A response regulators (Type-A RRs) play an important role in cytokinin-induced adventitious shoot formation in Pinaceae (Pinus pinea) [55]. Furthermore, it is known that crown size is determined by the expansion and division of cells in the apical and cambial meristems. Yang et al. [47] previously detected QTLs for growth traits, including crown width, in the interspecific hybrid clones of Taxodium species, which is widely grown in southern China for economic and ecological purposes, and reported that the three markers were annotated as a leucine-rich repeat receptor-like kinase (LRR-RLK) gene, which is specifically expressed in the vascular tissues and is vital for stem growth.

Eco-Physiological Traits
Among the morphological traits, three QTLs were detected for needle length/width ratio (Table 4; Figure 3). The provenance trial of Sakhalin fir demonstrated that the length/width ratio was negatively correlated with the altitude of seed source [17]. However, no studies related to the QTLs for needle length/width ratio were available for other conifers. For example, needle length was only investigated in Pinus elliottii var. elliottii × P. caribaea var. hondurensis hybrids [48]. In Populus and European beech, a major QTL for leaf length/width ratio was identified using linkage mapping [56,57]. No QTLs were detected for the other morphological traits (i.e., Norm_reac, LMA, SD, and SRN) examined in the present study. However, altitudinal clines have been detected in mature stands based on a provenance trial [17], suggesting that these morphological trends may have emerged at an advanced life-cycle stage.
Additionally, QTLs were detected for Bud_fl and Freez_tol (Table 4; Figure 3). A QTL for Bud_fl was also detected in a previous study of Sakhalin fir [17]. Bud phenology traits generally exhibit high heritability and are often controlled by several QTLs. A positive correlation between Bud_fl and growth traits was observed, indicating that the saplings with later Bud_fl tended to be larger in growth, which supports the data of a previous study on Sakhalin fir [19]. The Freez_tol was also found to be positively correlated with growth traits (Figure 2; Figure S2). The timing of acquiring Freez_tol was associated with the altitudinal adaptation of Sakhalin fir [25,58]. The locations of the QTLs for Bud_flush and Freez_tol were relatively close in the P336 map-LG5, with a distance of~10 cM (Figure 3). QTL clusters for height, density, and wood properties were also previously reported in Douglas fir [52]. These QTL clusters potentially represent pleiotropic effects or may be evidence for clusters of linked genes [59]. Because generating large amounts of meaningful data related to photosynthesis is easy, the data on chlorophyll fluorescence, including Φ II , NPQ, and Φ NO , can be directly applied to eco-physiological studies [60]. Φ II is an indicator of photosynthetic activity [29]. Because NPQ is associated with functional traits required to avoid photoinhibition [61] and Φ NO represents the non-regulated energy dissipation at the PSII centers, these parameters might be important for trees living in high-altitude zones. Although we did not detect any QTLs for Φ II , QTLs for NPQ and Φ NO were readily observed (Table 4; Figure 3). In general, the inhibition of photosynthesis by low temperature results in a considerable excess of energy, leading to photodamage. Specifically, selection to avoid photoinhibition should be stronger than that in low-altitude zones. Therefore, energy dissipation of heat and non-regulated energy might be essential for adaptation in high-altitudinal zones.

Conclusions
To elucidate the genetic basis of altitudinal adaptation in Sakhalin fir, a QTL analysis with linkage maps was conducted for 15 functional traits potentially associated with altitudinal adaptation. Similar to a previous study, QTLs with moderate PVE were detected for growth traits; ten were linked to CR17. Notably, one QTL for CR17 might function as a type-A RR, which plays an important role in cytokinin-induced shoot elongation. The QTLs for morphological traits (needle length/width ratio) and phenology (Bud_fl and Freez_tol) were also identified. Furthermore, our results suggest that two significant QTLs may be genetically controlling Φ NO in Sakhalin fir. Reproducibility of the QTLs detected in the present study should be examined in the future because we conducted the analysis with a single pedigree, in a single location, and there was no replication of segregated populations. Nevertheless, we believe that these findings might be an important initial step for detecting candidate genes for altitudinal adaptation in conifers because altitudinal adaptation has been demonstrated by a reciprocal transplant experiment [18], and a segregation population was produced by a control cross between F 1 hybrids with grandparents that were high-altitude and low-altitude genotypes [19].

Data Availability Statement:
The raw ddRAD-seq data generated in this study were deposited in the DNA Data Bank of Japan (DDBJ) under accession number DRA012397.