1. Introduction

Interspecific hybridization, by combining divergent genetic background, often generates genome-wide novel epistatic interactions and enhances phenotypic diversity and ecological specialization at a rate surpassing that of mutation alone (Lewontin and Birch, 1966; Rieseberg et al., 2003). Consequently, hybridization is recognized as a catalyst for speciation and evolutionary innovation (Anderson and Stebbins Jr, 1954; Arnold, 1997; Rieseberg, 1997). Hybrid speciation can proceed either via allopolyploidization or through homoploid hybrid speciation, where the hybrid retain the same ploidy level as its parents. While allopolyploidization could create rapid reproductive isolation between the hybrids and parents, homoploid hybrid speciation is a multi-generation process involving additional mechanisms, such as ecological divergence or chromosomal rearrangements, to establish reproductive isolation, making both the documentation and mechanistic understanding of this mode of hybrid speciation more challenging (Abbott et al., 2010; Rieseberg, 1997; Schumer et al., 2014). As such, the evolutionary significance of homoploid hybrid speciation remains debated, particularly regarding the extent to which hybridization directly contribute to reproductive isolation (Nieto Feliner et al., 2017; Schumer et al., 2014, 2018).

Theoretical expectations posit that ecological divergence, an important form of reproductive isolation, is a key factor in homoploid hybrid speciation. Without it, hybrid lineages risk being assimilated by parental species at early stages before advanced-generation hybrids can form (Buerkle et al., 2000; Gross and Rieseberg, 2005). Early-generation hybrids thus must exhibit phenotypic variation that facilitates niche transitions under new selective regime, otherwise homoploid hybrid speciation is unlikely to proceed successfully. Indeed, most of the documented homoploid hybrid species occupy either novel or intermediate niches relative to progenitors (Arnold et al., 2012; Gross and Rieseberg, 2005; Mao and Wang, 2011; Nevado et al., 2020; Wang et al., 2021). In such settings, barriers to genetic exchange evolve as a result of ecologically based divergent natural selection (Rundle and Nosil, 2005). However, knowledge about selection agents and corresponding phenotypic responses that drive postzygotic isolation and hybrid-specific adaptation remain limited.

The response to ecological selection can be expressed at morphological, physiological and gene expression levels. While morphological traits have traditionally been the focus of studies on hybridization, physiology remains underrepresented—likely due to the practical difficulties measuring them in natural settings. Similarly, the role of gene expression in hybrid and ecological speciation is still relatively unexplored, despite its potential to act as a crucial interface between genotype and phenotype and its ability to uncover ecologically relevant traits not apparent from morphology alone (Pavey et al., 2010). Comparative transcriptomics offers a powerful approach for identifying gene actions involved in speciation and provides unique insights into divergence associated with the colonization of new environments (Pavey et al., 2010; Runemark et al., 2024). However, despite its promise, few studies have applied this method to homoploid hybrids and their parental species—with notable exceptions in Helianthus (Lai et al., 2006), Senecio (Hegarty et al., 2009) and Argyranthemum (White et al., 2023), highlighting the need for broader application of this approach in speciation research.

Pinus densata is a keystone forest species on the southeastern Tibetan Plateau, where it forms extensive forests at elevations ranging from 2700 to 4200 m above sea level (Mao and Wang, 2011). The species represents the highest altitudes of all species of the subgenus Pinus in Asia (Mirov, 1967). Genetic evidence suggests it evolved from hybridization between P. tabuliformis and P. yunnanensis in the late Miocene (Song et al., 2002; Wang et al., 2011; Wang and Szmidt, 1994; Zhao et al., 2020). The three species are now largely allopatric, with P. tabuliformis distributed in northern to central China, P. yunnanensis in southwest China, and P. densata in between them (Mao and Wang, 2011). Ecological niche modeling has revealed distinct niche divergence in P. densata as reflected on key environmental factors such as frost frequency, diurnal and seasonal range of temperatures, growth season length, vapor pressure and UV radiation. These variables are characteristics of high alpine habitat and are strongly associated with habitat-specific adaptation in the three pine species (Mao and Wang, 2011). Transplantation experiments further demonstrated distinct habitat preferences as each species performed better in its native environment (Zhao et al., 2014). At the P. densata site, although all populations suffered high mortality, the survival rate of P. densata local population was more than 2-fold that of P. tabuliformis and P. yunnanensis by the third year, illustrating its significant adaptive advantage under the high plateau conditions (Zhao et al., 2014). While existing research has shed light on the patterns of genetic and ecological divergence in P. densata, the physiological and molecular mechanisms underpinning its superiority and direct evidence linking hybridization to environmental adaptation have yet to be fully explored.

The most direct approach for validating the role of hybridization in ecological transitions involves the creation of artificial hybrids and common garden testing. If artificial early-generation hybrids demonstrate traits similar to natural hybrids, this would support the hypothesis that hybridization facilitates ecological adaptation. Such approaches have successfully reconstructed the diverse phenotypes and their selective agents in several plant homoploid hybrid species, including Helianthus (Gross and Rieseberg, 2005; Rieseberg et al., 2003; Rosenthal et al., 2005), Senecio (Hegarty et al., 2009; Nevado et al., 2020), Iris (Arnold et al., 2012), Ostryopsis (Wang et al., 2021), and Argyranthemum (White et al., 2023). In this study, we investigate the role of hybridization in P. densata's ecological transition through comparative analyses involving the two putative parental species and their artificial F1 hybrids in a common garden within P. densata's natural habitat on the Tibetan Plateau. After the seedlings reached 5-years age, we conducted growth and physiological measurements across four seasons, along with RNA sequencing and gene expression analysis. By examining the patterns of phenotypic expression across a large set of physiological traits, we aim to address the following questions: 1) Does hybridization drive the adaptation of P. densata to high-altitude environment? 2) What are the specific physiological traits that enable P. densata's superior growth on the Tibetan Plateau? and 3) How does gene expression in P. densata differ from its parental species, and what role does hybridization play in this divergence? We hypothesize that adaptive physiological traits and gene expression patterns arising from hybridization will be retained and potentially enhanced during hybrid speciation, thereby contributing to hybrid fitness and promoting ecologically driven reproductive isolation.

Field studies on conifers in remote locations are challenging due to their long reproductive cycles, slow growth, and logistical constraints, which often resulting in high mortality and small sample sizes. Despite these challenges, our study yielded valuable insights into physiological differentiation among closely related pine species and the underlying genetic regulation driving high-altitude adaptation. These findings enhance our understanding of hybridization-derived novelty and the broader mechanisms shaping local adaption in conifer trees.

2. Materials and Methods

A more detailed Materials and Methods is provided in Appendix Supplementary data.

2.1. Plant materials and common garden experiment

This study was conducted in a common garden we established within the experimental station of Xizang Agricultural and Animal Husbandry University in Linzhi (29°40′N, 94°20′E, 2970 m above sea level), Tibet. The site was selected to represent the typical growth environment of Pinus densata on the Tibetan Plateau. The region has a semi-humid monsoon climate, and differs from the central regions of P. tabuliformis and P. yunnanensis on several climate variables, most distinctly the higher ground frost frequency, mean diurnal range and soil organic carbon, and lower heat energy (annual mean temperature, maximum temperature of warmest mouth, and growing degree days) and vapor pressure, which are characteristics of alpine environment (Mao and Wang, 2011).

The common garden hosted F1 hybrids from controlled crosses of Pinus tabuliformis x P. yunnanensis (hereafter F1s), and multiple provenances of P. densata, P. tabuliformis and P. yunnanensis. A series of controlled crossing experiments were initiated in 2008, involving P. tabuliformis of Ningcheng origin (42°16′N, 118°58′E) as maternal parents and P. yunnanensis of Kunming origin (24°58′N, 102°37′E) as paternal parents. The resulting F1 hybrid seeds were sown in the common garden in 2010, along with seeds collected from natural stands of the other three pine species. The planting followed a randomized complete block design. However, due to high rate of mortality, only a small portion of the seedlings survived to 5-years age when we started the 1st physiological measurements in Sep. 2015 (Table 1). After verifying the hybrid nature of all F1s using single nucleotide polymorphisms (SNPs) (see section 2.4 on RNA-seq), we included four trees of F1s (see Results 3.1), each from a different parental combination, four trees of P. densata, eight trees of P. tabuliformis and three trees of P. yunnanensis as the final sample set for comparative analyses. A suit of growth and physiological measurements were carried out over four seasons including Sep. 2015, Jan., May and Aug. 2016 to represent performance in autumn, winter, spring and summer, respectively (Table 1).

2.2. Growth and physiological performance

To characterize the physiological properties of the four groups of pines (i.e. the three pine species and F1s), we scored 44 traits (Table 2), including seven related to growth and soluble substances, 20 associated with photosynthetic activity, six related to isotopes, four of stress tolerance substances, and seven parameters associated with diurnal variations in photochemical efficiency. Measurements were taken across multiple time points of a day and/or a season for four seasons (Table S1). Detailed descriptions of each measurement are provided in Appendix Supplementary data. Here we briefly list the main measurements.

2.2.1. Primary metabolites and photosynthetic traits

To assess the physiological status of the plants, we measured soluble sugars, total protein content, and chlorophyll a and b and total chlorophyll contents from fresh needles, with five replicates for each individual (see Appendix Supplementary data). Soluble sugars were quantified as indicators of carbon assimilation and osmotic balance. Total protein content reflected metabolic activity and protein biosynthesis. Chlorophyll pigments were assessed to determine light-harvesting capacity. Moreover, to assess the photochemical efficiency and potential stress-induced damage to photosystem Ⅱ (PSII) under high-altitude conditions, chlorophyll fluorescence was measured by using a portable modulated fluorometer (Mini-PAM, Heinz Walz, Germany). The maximum quantum yield of photochemistry of PSII was determined in needles that had been subjected to a minimum of 2-h of darkness. Initial fluorescence (F_o) was measured under dim illumination, and the peak fluorescence yield (F_m) was measured subsequent to exposure to a 0.8 s pulse of saturating light at an intensity of 8000 μmol m⁻² s⁻¹. The maximum potential quantum efficiency of PSII was defined as F_v/F_m═(F_m-F_o)/F_m (Genty et al., 1989).

The rapid light curves (RLCs) were established to characterize the photosynthetic response and acclimation to ambient light of the samples following the procedure of Eilers and Peeters (1988). RLCs were constructed by exposing the needles to nine progressively increasing levels of photosynthetically active radiation (PAR), without prior dark acclimation. The PAR levels applied were 0, 36, 107, 215, 343, 509, 697, 1, 027, and 1387 μmol photons m⁻² s⁻¹, and measurements including electron transport rate (ETR) and maximum fluorescence under light-adapted conditions (F_m^') were taken using the Mini-PAM fluorometer. The non-photochemical quenching (NPQ) at each PAR was calculated as: NPQ = (F_m - F'_m)/F_m^' (Maxwell and Johnson, 2000). The duration of the irradiance steps was 30 s. All measurements were conducted on sunny days between 9:00 to 12:00 am, and with three replicates of each sample. Based on the RLCs, several photosynthetic parameters were derived including the maximum electron transport rate (ETR_m), representing the plateau of the RLC; the initial slope (alpha) of the light response curve, indicative of the efficiency of light utilization at low irradiance levels; and the light saturation point (E_k), derived as the ratio of ETR_m to alpha.

2.2.2. Gas exchange

To characterize diurnal variation in photosynthesis, we measured gas exchange of needles with an LI-6400XT photosynthesis system equipped with an LI-6400-01 CO₂ injector (LI-COR, Inc., Lincoln, Nebraska, USA). The measurements were performed for 1–3 days in each season during 9:00–18:00, at 1-hr intervals, to establish the light and CO₂ response curves. Measurements were made on 3–5 replicates for each sample with a 2 × 3 cm² cuvette head and a LI-6400-02B red-blue LED light source (10% blue light) under 13 different light irradiances and CO₂ concentrations. Alongside, steady-state rates of leaf net photosynthesis (A_n, μmol m⁻² s⁻¹), stomatal conductance to water vapor (g_s, μmol m⁻² s⁻¹), transpiration rate (TR, μmol m⁻² s⁻¹), leaf-air temperature difference (dT) and ratio of intercellular CO₂ to ambient CO₂ concentration (C_i/C_a), were determined. These indicators help characterize diurnal and seasonal variations in photosynthetic capacity and stomatal regulation.

To evaluate the photosynthetic capacity and light-use efficiency of the plants, important light response parameters, such as the light compensation point (LCP); maximum net photosynthetic rate (A_max), photosynthetic capacity (P_m), and light saturation point (LSP) were calculated from the light response curve using a mathematical model as in Ye (2007) (see Appendix Supplementary data for details).

To assess the biochemical limitations and regulatory mechanisms of photosynthesis, CO₂ response parameters, such as the maximum rates of carboxylation (V_cmax), maximum ribulose 1, 5-bisphosphate regeneration rate (J_max), ratio of ribulose 1, 5-bisphosphate regeneration rate to carboxylation rate (J_max/V_cmax), CO₂ compensation point (Г), mesophyll conductance (g_m), mitochondrial respiration in light (R_d) and photorespiration (R_p), were calculated from the CO₂ response curve using a non-rectangular hyperbola model of photosynthesis (Ethier and Livingston, 2004; Ethier et al., 2006) (see Appendix Supplementary data for details). These indicators provide critical insights into carbon assimilation potential, enzymatic activity, and internal CO₂ diffusion processes.

2.2.3. Carbon concentrations and water utilization

Long-term water utilization analysis relied on data from carbon isotope ¹³C in needles using an isotope ratio mass spectrometer (DELTA V Advantage; Thermo Scientific, Inc., Waltham, USA). Following combustion in the elemental analyzer, needle carbon content (C) was determined simultaneously with the elemental analyzer (Flash, 2000 EA-HT, Thermo Scientific, Inc., Waltham, USA). For needle ¹³C composition (δ¹³C_leaf), carbon isotopic values were expressed in δ notation relative to the Vienna Pee Dee Belemnite (VPDB) standard. Needle photosynthetic carbon isotope discrimination (Δ¹³C) was calculated as: Δ¹³C = (δ¹³C_air-δ¹³C_leaf)/(1+δ¹³C_leaf/1000) (Farquhar and Richards, 1984; Farquhar et al., 1989). The assimilation weighted ratio of intercellular (C_i) to atmospheric (C_a) CO₂ mole fractions of photosynthesizing needles (C_i/C_a) was estimated by a linear relationship of Δ¹³C to C_i/C_a (Farquhar et al., 1982). This relationship was used to calculate the intrinsic water use efficiency (WUE_long) of needles (Farquhar and Richards, 1984), which is indicative of long-term trends in internal regulation of carbon uptake and water loss of plants.

2.2.4. Nitrogen concentrations and δ¹⁵N analysis

Nitrogen concentrations and δ¹⁵N in needles were determined with an elemental analyzer (Flash, 2000 EA-HT) coupled to an isotope ratio mass spectrometer. δ values of isotope ¹⁵N were expressed in units of per mil (‰), and calculated as: δ (‰) = (R_sample/R_standard −1) x 1,000, where R_sample and R_standard are the nitrogen isotopic ratios ¹⁵N/¹⁴N of the sample and standard, respectively. A positive or negative δ¹⁵N value indicates enrichment or depletion, respectively, of the heavy isotope relative to the light isotope, with respect to atmospheric nitrogen. Photosynthetic nitrogen-use efficiency (PNUE, μmol CO₂ g⁻¹ N s⁻¹) was calculated as the ratio of A_max to needle nitrogen content on a projected area basis (g [N] m⁻²). This parameter serves as an indicator of the carbon gain per unit nitrogen investment, thereby reflecting the efficiency of nitrogen utilization in the photosynthetic apparatus.

2.2.5. Antioxidant enzyme activities and flavonoid content

To assess oxidative stress levels and the antioxidant defense capacity of the plants, we measured the activities of key antioxidant enzymes and quantified oxidative damage indicators. The activity of superoxide dismutase (SOD) and peroxidase (POD) were determined following the method of Nakano and Asada (1981) using fresh needles. The flavonoid content was measured at 510 nm with a spectrophotometer, and the content of malondialdehyde (MDA) was determined by the thiobarbituric acid (TBA) test following the protocol at Nanjing Jiancheng Bioengineering Institute (http://www.njjcbio.com/, Nanjing, China).

2.3. Feature selection of physiological trait

To gain an overview about the differentiation among the four groups of pine samples, we first conducted a principal component analysis (PCA) using R package ade4 (Chessel et al., 2004) on 100 measurements (Table S1), which revealed clear differences among the groups (see Results section). To identify the variables that contributed significantly to inter-groups differentiation, we proceeded with a feature selection using classification trees of the R-package rfPermute (Archer and Archer, 2016), which is an extension of the package randomForest (Liaw and Wiener, 2002). The random forest algorithm implements Breiman's forest algorithm for classification and regression. We divided the four groups of pine trees into three classification models (Model 1, 2 and 3) and one regression model (Model 4) (Table S2). Model 1 compares F1s vs. the two parental species. Model 2 compares P. densata with the two parental species. Model 3 considers F1s and P. densata as one group and compares them to the two parental species. Model 4 is a regression model for tree height that includes all four pine groups. Tree height is taken as the dependent variable, while other indices are considered as independent variables.

Following the identification of significant variables by the four models, we performed more detailed comparative analyses on them, including conducting PCA and multiple comparisons across groups based on the non-parametric Kruskal–Wallis test using agricolae package (de Mendiburu and de Mendiburu, 2019). The Kruskal–Wallis test was chosen over parametric alternatives (e.g., ANOVA) due to its lower sensitivity to assumptions of normality and equal variance, making it more robust for small sample sizes. Its rank-based approach is better suited for data with skewed distributions or outliers, which are common in our dataset.

2.4. RNA-seq and analyses

Fresh needles from 17 to 16 trees were collected in Jan. and May 2016 to characterize gene expression (Table 1). mRNA library preparation and paired-end sequencing (2 x 125 bp) on an Illumina Hi-Seq 2500 were performed by Beijing Yuanyi Gene Technology Co., Ltd. (Beijing, China). The average sequencing output per sample was 8.47 Gbp, with a range of 4.22–31.23 Gbp (Table S3). Raw reads were preprocessed with Fastp (v.0.21.0) (Chen et al., 2018) to remove adapters and low-base-quality sequences. The resulting clean reads were aligned to the reference genome of P. tabuliformis (Niu et al., 2022) using STAR (v.2.7.8a) (Dobin et al., 2013). More details are presented in Appendix Supplementary data.

We sorted the sequence alignment using SAMtools (Danecek et al., 2021). PCR duplicates were removed using SAMtools 'rmdup'. Variant calling was performed using BCFtools (Danecek et al., 2021). Variants with quality values below 25, SNPs within 3 bp of an indel were filtered out. SNPs with a minor allele frequency above 0.05 and genotype missingness below 0.1 were kept using VCFtools (Danecek et al., 2011).

The normalized transcripts per million reads (TPM) of each sample were estimated using RSEM (v.1.3.3) (Li and Dewey, 2011). Differences in expression between genes were evaluated using DESeq2 (Love et al., 2014), and compared at the seasonal level. To characterize gene expression patterns in Pinusdensata and F1 hybrids, we calculated the mid-parent expression values (MPV) based on the mean TPM from all pairwise comparisons between samples of parental species. We also performed pair-wise gene expression comparisons between F1s, P. densata, P. tabuliformis, and P. yunnanensis. These analyses were performed separately for two seasons, Jan. and May 2016 (Table 1 and Fig. S1). Differentially expressed genes (DEGs) were identified based on criteria of a false discovery rate (FDR) < 0.05. The DEGs identified between hybrids and MPV were further categorized into upregulated and downregulated genes. Genes had a fold change > 2 but an adjusted p-value (P_adj) > 0.05 were labeled as uncategorized.

2.5. Inheritance mode of gene expression and functional annotation

Based on patterns of differential gene expression between hybrids and parental species, genes were classified into five categories. Regardless of statistical significance, any pair of genotypes with expression differences of less than 2-fold were classified as exhibiting conserved inheritance, whereas genes whose expression in hybrids deviated by more than 2-fold (P_adj < 0.05) from either parent were categorized as having nonconserved inheritance. Nonconserved genes were further divided into additive, dominance, or overdominance inheritance in hybrids. Specifically, genes whose hybrid expression fell between the expression levels of the two parents (i.e., lower than one parent but higher than the other) were classified as additive. Genes with hybrid expression similar to that of one parent were categorized as dominance, while genes whose hybrid expression exceeded or was lower than both parents were classified as overdominance (or transgressive). For genes exhibiting nonadditive action, but the statistical power was insufficient to confidently assign them to any of the aforementioned categories, they were then labeled as uncategorized (McManus et al., 2010; Rapp et al., 2009; Swanson-Wagner et al., 2006).

Differentially expressed genes were further analyzed for gene ontology (GO) enrichment using KOBAS databases (Xie et al., 2011). The GO terms that exhibited a corrected P-value ≤ 0.001 were considered to be significantly enriched.

3. Results 3.1. Genetic verification of F1 hybrids

This study was conducted in a common garden that represents the typical habitat of Pinus densata on the Tibetan Plateau. Four groups of samples were included, P. densata, P. tabuliformis, P. yunnanensis and F1 hybrid of P. tabuliformis × P. yunnanensis (Table 1). To validate the hybrid status and exclude mislabeling of F1s from controlled crosses, we performed a PCA on the SNPs recovered from RNA-seq of 16 individuals sampled in May 2016. Following rigorous filtering, we obtained 29,491 high-quality SNPs, and genetic clustering by PCA showed clear separation of the four groups of samples, with the F1 hybrids genetically positioned between P. tabuliformis and P. yunnanensis (Fig. S2), supporting the hybrid nature of the F1 samples and their genetic uniqueness. This validates their use in the following comparative analyses.

3.2. Divergence in growth and physiological traits among Pinus densata, P. tabuliformis and P. yunnanensis

In this study, we measured 44 growth and physiological traits over four seasons in the common garden (Table 2). Altogether, 100 variables were scored for each sample (Table S1), and used in a preliminary analysis to understand the pattern of differentiation among the four groups of samples (Fig. S3). We then used four random forest models to identify significant traits contributing to the separation among groups, with each model selecting 11–20 parameters (Table S2). In total, 32 parameters were identified by the four models and used as the final dataset for multiple comparison analyses among groups. The overall pattern of divergence over the 32 parameters is shown in Fig. 1, in which P. densata and F1 hybrids are closer together in the PCA space and are separated from P. tabuliformis and P. yunnanensis along PC1 (growth, photosynthesis and oxidative stress tolerance) and PC2 (needle water content and photosynthetic activity), respectively. The first two PCs explained 75% of the total variance in this dataset. The clear advantage of P. densata in growth is reflected in height, ground diameter and needle number (Figs. 1 and 2A). This superiority could have resulted from differentiated physiological processes, which we examined, by main category, in details below.

3.2.1. Photosynthetic efficiency and photoprotection

We analyzed 12 photosynthesis-related traits to characterize differences in photoreaction in the pine species (Fig. 2B). These included traits that specify light-harvesting capacity (chl a, chl b, total chl), photosynthetic efficiency (ETR_m, alpha), tolerance to high light intensity (LSP, E_k), photosynthetic capacity (P_m, A_max), and Calvin cycle dynamics (J_max, J_max/V_cmax, Г). The three species showed significant differences on all these traits. Piuns densata is characterized by lower chlorophyll content, higher photosynthetic efficiency and tolerance to high-light conditions, and lower CO₂ compensation point than the other two species (Fig. 2B).

We further analyzed plants' photosynthetic efficiency and photoprotection from low to high photosynthetic active radiation (PAR) by examining the NPQ-PAR and ETR-PAR curves (Fig. 3), along with the maximal photosynthetic electron transport rate (ETR_m, Fig. 2B). In total, 207, 171, and 171 datapoints were collected for Jan, May and Aug. 2016, respectively, and used to construct the NPQ and ETR responses shown in Fig. 3. In summer season (Aug. 2016), Pinus densata had an optimal energy investment ratio, evidenced by a higher ETR_m than the other two species (Fig. 2B), a larger NPQ pool (Fig. 3), and utilization of more energy for photosynthesis rather than for photoprotection through heat generation (i.e., higher NPQ) as in the other two species (Fig. 3). In winter season (Jan. 2016), P. densata exhibited lower ETR and higher NPQ, indicating a greater energy investment in photoprotection and less energy used for photosynthesis. This pattern was particularly strong in P. yunnanensis. Generally, higher NPQ would result in a lower photosynthetic rate (A_n) during winter unless the plant possesses a strong CO₂ fixation capacity or its photoprotection mechanisms optimize energy utilization to sustain high A_n levels (see next section). Overall, P. densata displayed a balanced capacity between photosynthesis and photoprotection, reflecting its unique ability to optimize light energy utilization.

3.2.2. Diurnal variation in photosynthesis

Because of the pronounced diurnal variation in temperature and radiation on high plateau, we analyzed the daily variation in photosynthesis among the three pine species over three seasons. Measurements were taken with 1-hr intervals across 10 time points from 9:00–18:00, with a total of 190, 140 and 243 datapoints collected for Sept. 2015, Jan. and May 2016, respectively. In spring (May 2016), Pinus densata maintained higher net photosynthetic rate (A_n), stomatal conductance (g_s), and transpiration rate (TR) for most of the daylight period, a lower leaf-to-air temperature difference (dT), and a lower ratio of intercellular to ambient CO₂ concentration (C_i/C_a), indicating optimal physiological activity (Fig. 4). In autumn (Sep. 2015), P. densata and P. tabuliformis showed comparable photosynthetic induction rates (i.e. the initial slope of net photosynthetic rate), which were higher than P. yunnanensis. At milder light intensities (from 9:00 to 12:00, and after 16:00), the A_n of P. densata was highest, while the C_i/C_a ratio was correspondingly lower. Conversely, during high light intensity periods (from 12:00 to 15:00), a decline in A_n and an increase in the C_i/C_a ratio were observed, likely as a rapid response to light and temperature changes. Additionally, under high light intensity, P. densata exhibited a lower dT and a higher g_s relative to the other two species, thereby demonstrating superior A_n performance (Fig. 4). In winter (Jan. 2016), P. densata and P. tabuliformis had little photosynthetic activity during the cold morning hours (from 9:00 to 13:00), but P. densata was able to restart photosynthesis under suitable temperature conditions (after 13:00), whereas P. tabuliformis did not exhibit photosynthetic activity throughout the day. It is noteworthy that P. yunnanensis had slightly higher photosynthetic activity during most of the daylight hours, which might be a behavior less sensitive to cold conditions, making it prone to temperature damages, like frost, in alpine environments.

3.2.3. Stress tolerance related traits

The carbon (Δ¹³C) and nitrogen (δ¹⁵N) stable isotope compositions in plants hold significant implications for environmental adaptation. Variations in leaf Δ¹³C are used to study plant responses to drought and mechanisms underlying acclimation to changing conditions (Cernusak et al., 2013). Leaf δ¹⁵N values are useful for studying nitrogen cycling in ecosystems, plant nitrogen use efficiency, and the effects of nitrogen availability on plant (Spangenberg et al., 2021). Our analysis of Δ¹³C and δ¹⁵N in the three pine species showed lower Δ¹³C discrimination values in P. densata (Fig. 2C), indicating that the other two species experienced drought stress relative to P. densata in the common garden site. On the δ¹⁵N values, P. densata were intermediate between the other two species.

Malondialdehyde (MDA) is recognized as a biomarker for oxidative stress and for assessing the extent of lipid peroxidation in cellular membranes. The concentrations of MDA in Pinus densata were higher than those in the other two species (Fig. 2D), suggesting a heightened antioxidant capacity and potential for coping with environmental stress. The total protein content in plant is related to the plant's growth and metabolic activities. We detected a higher value in P. densata (Fig. 2A), potentially reflecting its enhanced growth vigor and metabolic efficiency.

3.3. Phenotypic similarity between Pinus densata and the F1 hybrids

Over all the 32 selected traits, the F1 hybrids exhibited a high similarity to Pinus densata on 26 traits, resulting in an overall affinity with P. densata in growth and physiological traits as shown in the PCA plot (Fig. 1). For example, with respect to the photoreactive properties of photosynthesis (Fig. 2B), both P. densata and the F1s demonstrated good low-light harvesting ability (chl a, chl b, total chl), higher photon energy use efficiency (ETR_m, alpha), high-light tolerance (LSP, E_k), and high photosynthetic capacity (P_m, A_max). Additionally, the F1s also showed similar performance to P. densata in photoprotective capacity (NPQ) and diurnal pattern of photosynthesis (Figs. 3 and 4). However, on water content, ETR_m, J_max, J_max/V_cmax and MDA of May 2016, the F1s differed significantly from P. densata (Fig. 2). The overall parallel performances of F1 hybrids and P. densata bolster the notion that hybridization can enhance adaptability to new environments.

3.4. Differential gene expression

Given the significant divergence in trait values among the four groups of samples, we delved into transcriptomic differences between the groups to investigate the potential regulatory basis for phenotypic divergence. Mapping and annotation recovered 33, 993 and 33, 431 expressed genes from the dataset of May and Jan. 2016, respectively. Analysis of gene expression over the two seasons revealed 13.41%–13.87% significant differentially expressed genes (DEGs) between P. densata and P. tabuliformis, 3.10%–3.36% between P. densata and P. yunnanensis, 6.64%–12.01% between P. tabuliformis and P. yunnanensis, and 19.14%–18.47% between P. densata and mid-parent expression values (MPV) of the other two species (Fig. 5A and B, Table S4). Among the DEGs between P. densata and MPV (6503 and 6106 genes for each season), 64.34%–70.41% were down regulated and 35.66%–29.59% up regulated, of which 149–108 (3.56%–2.51%) genes were lower and 9–25 (0.39%–1.38%) genes were higher than both P. tabuliformis and P. yunnanensis, respectively (Fig. 5C and E; Table S4).

For the F1s, 4.10%–6.07% and 2.15%–2.83% DEGs were detected when compared to Pinus tabuliformis and P. yunnanensis, respectively, and 8.25%–10.90% DEGs when compared to MPVs over the two seasons (Fig. 5A and B; Table S4). These values were much smaller than the corresponding values in P. densata. The DEGs between F1s and P. densata were 5.50%–2.43% (Fig. 5 A and B; Table S4). Similar to the pattern in P. densata, 61.75%–71.02% and 38.25%–28.98% of the DEGs in F1s when compared to MPV were down- and up-regulated, respectively, of which 21–23 (1.22%–0.91%) genes were lower and 3–4 (0.28%–0.38%) genes were higher than both parents, respectively (Fig. 5D and F; Table S4). The pronounced differences between F1 hybrids and the MPV suggest that hybridization can lead to rapid transcriptomic reprogramming.

In this study, all transcriptome data were mapped to the Pinus tabuliformis reference genome. While genomic differences due to interspecific divergence may introduce some systematic bias in transcriptome comparisons, mapping to the P. tabuliformis genome remains the most practical option, as the three pine species involved are the most closely related within the genus Pinus (Jin et al., 2021). We examined the mapping rates for all samples and found that most exceeded 90%, with only two samples had lower rates at 84%. Importantly, there was no apparent bias in mapping rate among species (Table S3). Therefore, we consider the mapping results to be appropriate for comparative analyses among our samples.

3.5. Mode of inheritance of gene expression in Pinus densata and F1 hybrids

Further analysis of the gene expression patterns in Pinus densata and F1s vs. P. tabuliformis (Pt) and P. yunnanensis (Py) classified 47.63%–52.02% as conserved, 3.23%–4.21% as additive, 4.27%–10.42% as dominance, and 0.26%–1.45% as overdominance (or transgressive, Table 3). These values were consistent in P. densata and F1s over the two seasons. One noticeable difference between P. densata and the F1s is that P. densata had about 2-fold more genes in the dominance class and more than 2-fold of genes in overdominance class (Table 3). We found substantial overlap between F1s and P. densata across different categories of gene actions, e.g. 384–417 additive, 893–1312 dominance, and 9–28 overdominance loci were shared between the two groups, supporting the presence of shared molecular mechanisms underlaying their phenotypic novelty. Further classification of the dominance genes into Py- and Pt-like groups revealed another distinct difference between P. densata and the F1s, in that 63.17%–66.05% of the dominance genes were Py-like and 36.83%–33.95% were Pt-like in F1s, by contrast 86.18%–83.13% of the dominance genes were Py-like and 13.82%–16.87% were Pt-like in P. densata (Table 3 and Fig. 6). These results indicate a stronger influence of P. yunnanensis genes in the two hybrid groups, and higher dominance and overdominance gene actions in the advanced generation natural hybrid P. densata.

3.6. Functional annotations of dominance and overdominance genes

Functional annotation of the Py- and Pt-like dominance genes in Pinus densata and F1s revealed significant enrichment in biological processes related to light response, oxidation reduction and tolerance to biotic and abiotic stress. P. densata showed a higher enrichment of Py-like genes (Fig. 6), whereas F1 hybrids had a relatively higher enrichment of Pt-like genes (Fig. S4), as expected from the retention of more P. tabuliformis genes in the F1s. Enrichment analysis of the overdominance genes in P. densata reveled their concentration in flavonoid biosynthetic process, UV response, drought response, and oxidative stress related processes, which are characteristic of adaptation to high plateau habitat (Fig. S5). This illustrates that non-additive gene actions in hybrids have a plausible role in ecological adaptation and could explain the physiological properties that enable P. densata to survive in the alpine environment that neither of its parents can tolerate. The generation of novel patterns of gene expression in hybrids, therefore, can result in evolutionary novelty in hybrid species.

4. Discussion

Previous transplantation experiments involving Pinus densata and its putative parents, P. tabuliformis and P. yunnanensis, have revealed its distinct fitness advantage on the Tibetan Plateau, as reflected in superior growth, survival and seed production potential (Mao et al., 2009; Mao and Wang, 2011; Zhao et al., 2014). However, direct evidence supporting hybridization-derived adaptation in P. densata has been lacking. In this study, we conducted a direct in situ comparison of physiological traits and gene expression in P. densata, F1 hybrids and the two parental species in a common garden on the Tibetan Plateau. Despite a limited sample size resulting from high mortality, our findings on surviving samples provide valuable insights into how hybridization has facilitated ecological transition in P. densata, enabling it to thrive in the alpine environment.

4.1. Ecophysiological adaptation to high plateau habitat in Pinus densata

Pinus densata exhibited superior growth compared to its parental species. We suspected that this growth advantage is driven by enhanced photosynthetic efficiency and stress tolerance. Indeed, physiological measurements revealed higher maximum electron transport rates (ETR_m) and light saturation points (LSP), reflecting improved capacity to utilize light energy under high irradiance (Côté and Platt, 1984; Gu, 2023). Additionally, the species maintained lower CO₂ compensation points (Г), suggesting more efficient carbon assimilation.

Survival in high-altitude environments requires balancing photochemical activity and photoprotection in response to intense radiation and temperature fluctuations (Demmig-Adams and Adams Ⅲ, 2006; Mao and Wang, 2011). Pinus densata showed a seasonal flexible adjustment of non-photochemical quenching (NPQ) and electron transport rates (ETR): in summer, it allocated more energy towards photochemical activity with higher ETR_m and lower NPQ; in winter, it shifted toward photoprotection with higher NPQ and lower ETR. This dynamic energy allocation enhances plants' resilience to thermal and irradiance-related stress (Adams et al., 2004; Ploschuk et al., 2014). Diurnal measurements further highlighted P. densata's physiological flexibility. Compared to its parents, it maintained higher net photosynthetic rates (A_n), stomatal conductance (g_s), and transpiration rates (TR), while exhibiting lower leaf-to-air temperature differences (dT) and intercellular-to-ambient CO₂ concentration ratio (C_i/C_a) throughout the day, particularly in spring and autumn. This combination supports more stable photosynthesis under fluctuating conditions (Kaiser et al., 2015). Notably, its capacity to resume photosynthesis during warmer afternoon hours in winter underscores its adaptive advantage in cold, variable climates. In contrast, P. tabuliformis showed limited photosynthetic activity in winter, while P. yunnanensis lacked precise regulation, potentially increasing frost vulnerability. These results suggest that an efficient CO₂ fixation capacity and optimization of energy utilization for photosynthesis and photoprotection contribute to the observed higher A_n levels in P. densata.

Furthermore, Pinus densata displayed traits associated with abiotic stress tolerance, including improved water-use efficiency, and heightened antioxidant response to oxidative stress from high light intensity and fluctuating temperatures. Such stress tolerance traits are essential for maintaining cellular function under abiotic stress (Ma et al., 2010; Tommasino et al., 2012). Taken together, these finding suggest that P. densata combines efficient low-light harvesting, enhanced photosynthetic performance, and robust photoprotection. The combined advantages of these traits likely support its vigorous growth and successful ecological transition from its parental species in the challenging environment of the Tibetan Plateau.

4.2. Hybridization contributes to the ecological transition

To determine whether the novel adaptation in Pinus densata originated from hybridization, we compared its performance with artificial interspecific F1 hybrids. Of the 32 growth and physiological traits identified as major contributors to intergroup divergence, 26 were highly similar between P. densata and F1s, reflecting a strong phenotypic affinity. The F1 hybrids displayed growth characteristics and seasonal/diurnal photosynthetic patterns closely mirroring P. densata, indicating that hybridization confers immediate physiological shifts relevant to high-altitude adaptation.

Our physiological and transcriptomic assessments were conducted at the 5-year age on individuals that have survived strong natural selection during early seedling establishment in the field—a phase marked by high mortality. While this limits sample size, it captures the effects of both hybridization and selection on trait variation. In contrast, if this study were done under mild greenhouse environment, greater trait variance would be expected in a larger F1 pool due to the absence of environmental filtering. The in situ design, although challenging, offers unique value in identifying trait combinations shaped by selection that are directly relevant to local adaptation.

Homoploid hybrid speciation is typically a multi-generation process, during which natural selection continuously adjusts the genetic and ecological properties of the hybrid population to its environment. Thus, although hybridization provides the initial genetic substance for ecological transitions, further refinement in the following generations, via e.g. divergent selection on segregating genotypes, would eventually buildup strong enough genetic differentiation and isolation to establish a new species. Such processes have been documented in other plant groups, including Helianthus (Gross and Rieseberg, 2005; Rieseberg et al., 2003; Rosenthal et al., 2005), Senecio (Nevado et al., 2020) and Ostryopsis (Wang et al., 2021). In outcrossing conifers, frequent introgressions is common between sympatric and parapatric sister species due to weak genetic incompatibility (Zhao et al., 2024). Evidence for this can be seen in crossing experiments between Pinus tabuliformis and P. yunnanensis, which yielded a seed rate as high as 24% (Zhao et al., 2014). Weak genetic barriers and frequent introgressions between species could potentially lead to genetic assimilation over time, unless there are geographical and ecological isolation barriers that facilitate pre- and postmating reproductive isolation. Our common garden transplantation experiments involving P. densata and the two parental species have revealed strong species-specific adaptation to local climate conditions (Mao et al., 2009; Mao and Wang, 2011; Meng et al., 2015; Xing et al., 2014; Zhao et al., 2014), suggesting that geographic isolation and ecological selection are both important in maintaining the species boundaries. Advancing from these findings, the present study provides direct evidence for the role of hybridization in generating phenotypic novelty and demonstrates how physiological adaptation to a novel high-altitude niche has been shaped by hybridization and subsequent selection in P. densata.

4.3. Hybridization promotes transcriptomic reprogramming

Ecological differentiation is typically polygenic and many of the traits of interest in speciation are shaped by complex multi-loci interactions, thus natural selection acting on admixed populations will likely take advantage of gene actions from both parents (Hegarty et al., 2009; Leal et al., 2020; Rieseberg et al., 2003; White et al., 2023). Indeed, gene expression patterns in Pinus densata and the F1s differentiated strongly from the parental species, with 19% of 34,000 examined genes in P. densata and 9%–11% in F1s deviated from MPV. The much higher proportion of DEGs in P. densata likely reflects the influence of natural selection in advanced-generation natural hybrids, driving enhanced divergence in gene expression during adaptive evolution. Furthermore, the deviation of F1 hybrids from the MPV indicates that hybridization does not merely produce intermediate phenotypes. Instead, it triggers complex gene interactions, leading to up- or down-regulation or dominance or transgressive expression. Depending on the function of a gene, increased stress tolerance or other adaptive performance can come from either increased or decreased expression regulations (Hegarty et al., 2009; Ng et al., 2019; Rieseberg et al., 2003; Rivera et al., 2021).

Analysis of inheritance mode of gene expression revealed that both additive (ca. 4%) and non-additive effects (5%–6% in F1s, 10%–12% in P. densata) are common genetic dynamics of the two hybrid groups, with dominance and overdominance loci more common in P. densata than in F1s. Relative to the whole 34,000 examined genes, dominance and overdominance (transgressive) loci in P. densata were ca. 10% and 1%, respectively, while in F1s the corresponding values were ca. 5% and 0.3%, respectively (Table 3). These proportions increase substantially if relative to the DEGs category. The dominance hypothesis posits that hybrid vigor arises from the complementation of favorable dominance alleles, masking the deleterious effects of recessive alleles (Thompson et al., 2021). This is consistent with findings in maize (Baldauf and Hochholdinger, 2023), rice (Shao et al., 2019), sorghum (Hashimoto et al., 2021), and oilseed rape (Ye et al., 2023), where dominance complementation has contributed to hybrid vigor and environmental adaptability.

We observed a parental-like expression bias in the dominance loci in both F1s and Pinus densata, with 65Py: 35 Pt in F1s and 85Py: 15 Pt in P. densata, suggesting that genes or gene regulation machine from P. yunnanensis contribute more to these hybrids. The predominance of P. yunnanensis is also detected in genetic composition of P. densata populations from its central to western range, reflecting high levels of introgression (Wang et al., 2011; Zhao et al., 2020). Parent bias in traits or gene expression has been commonly observed in plant homoploid hybrids, for example in Helianthus deserticola (Lai et al., 2006), Senecio squalidus (Hegarty et al., 2009), and Argyranthemum lemsii and A. sundingii (White et al., 2023), and is regarded as an operating mechanism improving fitness to a specific niche (Thompson et al., 2021).

Functional enrichment analyses revealed that these dominance genes are involved in biological processes such as light response, oxidation-reduction, and biotic and abiotic stress tolerance. Additionally, the overdominance (transgressive) genes identified in Pinus densata are enriched in pathways related to flavonoid biosynthesis, UV response, drought response, and oxidative stress. These biological properties are essential for coping with the intense UV radiation and diverse abiotic stressors typical of high-altitude habitats (Du et al., 2022; Rodríguez-Calzada et al., 2019; Sun et al., 2020). The enrichment of functions that reflect adaptation to high plateau environments strongly suggests that the dominance and transgressive expression observed in P. densata could be adaptive. While our findings provide a tentative link between gene actions and the physiological adaptation, direct causal relationships require further experimental validation. Nonetheless, these results enhance our understanding of how hybridization drives transcriptomic reprogramming and facilitates gene expression evolution in conifer genomes during ecological transition.

Gene expression profiles serve as powerful surrogate phenotypes, capable of revealing adaptive traits that are otherwise difficult to measure, and may play a key role in initiating ecological divergence during speciation (Pavey et al., 2010; Runemark et al., 2024). They also provide insight into the genetic architecture underlying ecological divergence, an essential component of homoploid hybrid speciation. However, this remains an underexplored area. A recent review (Runemark et al., 2024) identified only four studies comparing gene expression between stabilized plant homoploid hybrids and their parental species, all within Asteraceae (Helianthus, Senecio, Argyranthemum) (Brouillette and Donovan, 2010; Hegarty et al., 2009; Lai et al., 2006; White et al., 2023). Our study extends this limited body of work to conifers, a distant and previously unrepresented taxonomic group. A common pattern across these studies is the prevalence of non-additive gene expression—particularly complementary dominance and, to a lesser extent, transgressive expression. Notably, loci associated with ecological adaptation tend to be enriched among those showing transgressive expression, suggesting that novel gene expression patterns resulting from hybridization is an operating mechanism of ecological divergence in hybrid species.

5. Conclusion

Although adaptation to novel environments is often observed in homoploid hybrids, strict definition of homoploid hybrid speciation emphasizes whether such divergence arose at the species' origin or evolved after the hybrid species formation (Abbott et al., 2010; Schumer et al., 2014). While this distinction may be important from a theoretical standpoint, understanding the full evolutionary trajectory of novel adaptation requires a broader perspective. In this study, we conducted in situ physiological and transcriptomic comparisons of Pinus densata, its putative parental species, and synthetic F1 hybrids in a common garden on the Tibetan Plateau. We identified key physiological traits and gene expression patterns that underpin the superior growth and stress tolerance of P. densata, and provided clear evidence that hybridization played an important role in enabling its ecological transition. The strong phenotypic and molecular similarity between P. densata and F1s suggests that hybridization can rapidly generate novel trait combinations favorable for adaptation to new environments, which in turn allow for ecological diversification and, ultimately, speciation. However, speciation is a dynamic, multi-generational process. Continued selection and introgression further shape the genetic and phenotypic landscape of hybrid populations, leading to divergence between early and advanced generations. Finally, the widespread involvement of physiological and gene regulatory traits emphasizes the complexity of local adaptation in conifers, and highlights the power of integrative approaches in uncovering the mechanisms behind hybrid speciation and ecological innovation.

Acknowledgments

Genomic data processing and analyses were performed using resources provided by the Swedish National Infrastructure for Computing (SNIC), through the High Performance Computing Centre North (HPC2N) at Umeå University. This study was supported by the National Natural Science Foundation of China (32171816), and T4F program Sweden.

CRediT authorship contribution statement

Zhi-Chao Li: Writing – original draft, Visualization, Methodology, Investigation, Formal analysis, Data curation. Chao-Qun Xu: Writing – original draft, Visualization, Methodology, Investigation, Formal analysis, Data curation. Wei Zhao: Methodology, Formal analysis. Shuai Nie: Investigation, Formal analysis, Data curation. Yu-Tao Bao: Methodology, Investigation, Formal analysis. Hui Liu: Methodology, Formal analysis, Data curation. Zhen Xing: Resources, Investigation, Data curation. Jian-Feng Mao: Writing – review & editing, Supervision, Project administration, Methodology, Investigation, Funding acquisition, Conceptualization. Xiao-Ru Wang: Writing – review & editing, Project administration, Funding acquisition, Conceptualization.

Data availability statement

The raw transcriptome sequencing reads generated in this study have been deposited in the Genome Sequence Archive (GSA) at the National Genomics Data Center (NGDC), China National Center for Bioinformation, Beijing Institute of Genomics, Chinese Academy of Sciences, under accession number: CRA022689, and are publicly accessible at https://ngdc.cncb.ac.cn/gsa.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2025.05.009.