Adaptive Differentiation in Seedling Traits in a Hybrid Pine Species Complex, Pinus densata and Its Parental Species, on the Tibetan Plateau

Evidence from molecular genetics demonstrates that Pinus densata is a natural homoploid hybrid originating from the parent species Pinus tabuliformis and Pinus yunnanensis, and ecological selection may have played a role in the speciation of P. densata. However, data on differentiation in adaptive traits in the species complex are scarce. In this study, we performed a common garden test on 16 seedling traits to examine the differences between P. densata and its parental species in a high altitude environment. We found that among the 16 analyzed traits, 15 were significantly different among the species. Pinus tabuliformis had much earlier bud set and a relatively higher bud set ratio but poorer seedling growth, and P. yunnanensis had opposite responses for the same traits. P. densata had the greatest fitness with higher viability and growth rates than the parents. The relatively high genetic contribution of seedling traits among populations suggested that within each species the evolutionary background is complex. The correlations between the seedling traits of a population within a species and the environmental factors indicated different impacts of the environment on species evolution. The winter temperature is among the most important climate factors that affected the fitness of the three pine species. Our investigation provides empirical evidence on adaptive differentiation among this pine species complex at seedling stages.


Introduction
Natural hybridization plays an important role in plant evolution [1][2][3]. Hybridization facilitates adaptive evolution and generates new species either through allopolyploid or homoploid hybrid speciation (HHS) [4][5][6]. HHS is hybridization that occurs in the wild without a change experiment (April 2011-March 2013) were collected from the China Meteorological Administration (Table 1).
Twenty-four natural populations, 9 of P. tabuliformis, 8 of P. densata, and 7 of P. yunnanensis, were used to represent populations found in the primary geographical regions of each species (Table 2 and Fig. 1). Seeds were collected from each population when cones ripened naturally. In each sample population, 30 or more individuals (separated by at least 50 m) were selected for cone collection. More than 3000 seeds for each population were mixed in separate cloth bags and stored under low temperature and humidity conditions in darkness before  The meteorological data of each year were collected from China Meteorological Data Sharing Service System (http://www.cma.gov.cn/2011qxfw/ sowing. Only full and healthy seeds (determined by X-ray analyses) were used in the plantation experiment.

Experimental design
The experimental site was located in the nursery at Tibet University. Plants were cultivated in seedling beds with forest soil, and seedlings from the different populations were distributed according to a randomized complete block design. 50 seeds from each population were sown in a population plot with the spacing of 10cm×10 cm. Four replication blocks, each with the full set of populations, were exposed to natural selection factors without human intervention (except weeding in the autumn). We examined 16 quantitative traits involving growth, survival, phenology and needle color variation (Table 3). Surviving plants were measured from the approximate time of bud sprouts (April) to the time of bud set (October) over two years from early April 2011 to late March 2013. The stem height (SH, cm) and stem diameter 2 cm above the ground (D0, mm) were measured on 10 surviving individuals randomly from each plot (all plants in a plot were measured if the number of surviving plants was less than 10). Seedling phenology and viability traits that were measured for each experimental plot included the germination rate (GR), the survival rate (SR) in October, the bud set rate (BSR) in October, the survival rate after winter (SRAW) in April, and the second year-growth rate (SGR) in October. The ratios of seedlings with different needle colors (RLC) were determined after frost in October; the seedling colors observed were green, yellow, purple, and red.

Statistical analyses
To examine the effects of species and populations for each trait, we used a two-way nested analysis of variance (ANOVA) followed by Duncan's multiple range test with block and species as fixed factors and population as a random factor. Variables were logit-transformed to ensure data had normal distributions prior to ANOVA and Duncan's analyses. The general statistics for each trait in the population included the mean and the maximum and minimum values, which were used to describe the differences among species. The descriptions of each trait within each population were made using medians, quartiles, and amplitudes based on the average of each plot. All analyses were performed in the R (version 2.13.0) statistical and programming environment using the vegan package [26].
To evaluate the underlying dimensionality of the variables and to obtain an overview of the dominant patterns, we used PCA with a standardized matrix containing data on measured traits. Growth indicators were not considered in PCA for unavailable data from dead plants. Analysis using the ade4 package was performed using R statistics [27]. The PCA scores of the 1 st four components were used to calculate the species distance and a scatter diagram was generated using R-2.13.0 (ggplot2) to describe the dimensionality of the three pines in the 1 st two components.
Cluster analysis using arithmetic averages was performed for the 24 sample populations of the three pine species to estimate the similarity of populations for fitness related traits. Following correlation analyses among all indicators, 10 unrelated variables (SH of the 1st year, D0 of the 2nd year, GR and SR of the 1st year, SRAW of the 1st year, BSR of the 2nd year, RLC of the green, RLC of yellow, RLC of purple, and RLC of red) with correlation coefficients lower than 0.6 were chosen for cluster analysis. Cluster analysis was performed using the ward method in R statistics, and the Euclidean distances among populations were calculated from z-scores.
To examine habitat effects and adaptability variation, we calculated Pearson correlations between measured variables and primary habitat factors for the 24 sample populations of P. densata, P. tabuliformis, and P. yunnanensis. Data on three meteorological factors covering the last 50 years were collected (WorldClim, www.worldclimate.com), including three geographic factors (latitude, longitude, and elevation) and three meteorological factors (annual average temperature, annual precipitation, and average temperature in January). This multivariate correlation analysis was performed using SPSS 18.0 [28].

Results
Differentiation among the three pine species in the P. densata habitat on the Tibetan Plateau Because of the high elevation of the Tibetan Plateau, 15 traits differed significantly among P. densata, P. tabuliformis, and P. yunnanensis (Table 4). P. densata possessed higher values of SRAW, SR, SGR, RLC in purple and SH, but displayed intermediate values between parent species for BSR, D0, GR, RCL in yellow and red. P. yunnanensis had the lowest values for traits of GR, SR and SRAW but had the highest value of D0. P. tabuliformis possessed intermediate. values of GR, SR in the 1st year, RLC in green of the three pine apecies to conditions on the Tibetan Plateau. BSR over 2 years for P. tabuliformis were higher than in the two other pine species, while RLC in red, D0, SH and the SGR in the 2nd year were lower than in the other species.
Significant differences among populations were observed for P. densata (Table 4 and Fig. 2). This hybrid species had broader variations in SR in the 1 st year, SRAW in the 1 st year, BSR in the 1 st year, SGR in the 2 nd year, SH over 2 years, and RCL in four colors than parent species. Distributed in the northeast, population 15 had the highest GR, SR in the 1st year, and BSR in the 2nd year. Central populations 11 and 12 had the lowest survival rate but had high RCL in red seedlings. And the western populations 8, 9 and 10, had the best survival ratios and the fastest growth. Remarkably, population 15 had excellent seedling survivability, which was similar to the western populations.
For the two parental species, significant differences were observed among geographic populations within a species for 10 traits (Table 4 and Fig. 2). P. tabuliformis had less variation for each trait than P. densata. In population 24, BSR of the 1 st year was the highest. Population 18 had higher fitness and a higher SR in the 2 nd year. Populations 16 and 17 had the lowest BSR in the 1 st year and then grew faster with higher SH and D0 values over the next 2 years. P. yunnanensis had the broadest variation in GR, SR in the 2 nd year, BSR in the 2 nd year, and D0 over 2 years among the three pines. Populations 1 and 2 rarely did bud set in October, and all plants in population 1 died during the 1 st winter. population 4 displayed the highest value of SRAW in the 2 nd year. Population 7 possessed higher RLC in red, BSR over 2 years, and SGR for 2 years than the other populations.
Principal components analysis (PCA) on seedling tested traits revealed four principal components with eigenvalues greater than 1. The four principal components accounted for 74.7% of the variation. Individually, the traits of GR, SR over 2 years, and BSR over 2 years, had high loadings in the 1st component, while RLC in green had high loadings in the 2nd component ( Table 5). The first two components explained 52.1% of the variation (Fig. 3). Both P. tabuliformis and P. yunnanensis were clearly distinguished from each other. P. densata had a broad distribution that covered almost all of the area between the two parent species. Significant differences were observed between P. tabuliformis and P. yunnanensis. P. densata appeared to inherit adaptive characteristics of the parent species, suggesting that greater variation is present in the hybrid species Based on eight weakly related characteristics (Pearson correlation coefficient lower than 0.5), hierarchical cluster analysis using the ward method produced two clearly defined groups within the 24 populations (Fig. 4). Populations of P. tabuliformis and P. yunnanensis were completely separated. Populations of P. densata were divided into four clusters: P. densata east 1 (14,15), P. densata east 2 (13), P. densata central (11 12), and P. densata west (8,9,10). Populations of 13, 14, and 15 in P. densata were clustered in the P. tabuliformis group showing a  closer relationship to P. tabuliformis, while populations 8, 9, 10, 11, and 12 in P. densata were clustered in the P. yunnanensis group and showed a closer relationship to P. yunnanensis. Population 7, located in the northeastern-most region of the P. yunnanensis range, had little differentiation from P. densata. A relationship may exist between population 7 and P. densata.

Environmental association
For P. yunnanensis, annual average temperature and average temperature of the coldest month were associated with BSR over 2 years and RLC in red. Annual precipitation was significantly correlated with seedling survival (SR, SRAW), growth (D0 and SH), and needle traits (RLC) in the 1 st year. The average temperature of coldest month was positively correlated with the BSR over 2 years (Table 6). For P. densata, annual average temperature was significantly positively correlated with growth rate of the 2 nd years, annual precipitation was significantly correlated with seedling survival (SR, SRAW), growth (D0 and SH), and needle traits (RLC) of the 2 nd years. Similarly, the average temperature of the coldest month was significantly positively correlated with seedling survival (SR, SRAW), phenology (BSR) and growth rate (D0 and SH) ( Table 6).
For P. tabuliformis, the annual average temperature and the average temperature in January were significant correlation with phenology (BSR). Additionally, the average precipitation was significantly correlated to germination rate (GR) ( Table 6).

P. densata shows better fitness than parent species on the Tibetan Plateau
Speciation of P. densata was supposed to be the result of niche adaption [1,5,9]. Common garden experiment in hybrid habitat is one of the most direct approaches to investigation of juvenile fitness [29]. We found P. densata had the highest survival rate and the fastest growth rate compared to the two parental species which is indicative of an adaptive advantage over the parental species in the hybrid niche. Plant germination and seedling growth is environmentally demanding [30,31]. Traits for survival and development are the most direct indicators of seedling adaptation [32][33][34]. Although lower than P. tabuliformis in the GR, the survival rate (SR) of P. densata surpassed both parental species, Species with higher survival rate will have advantage in interspecific competitions [35,36]. Habitat factors on the Tibetan Plateau, such as precipitation, temperature, light, soil, and stress, significantly affected the fitness of seedling plants [9]. For traits associated with phenology, P. densata was intermediate between parent species. Bud set is a feature of woody plant for avoiding winter injury, which is induced by photoperiod and temperature [37][38][39][40][41]. The much earlier and higher ratio of bud set of P. tabuliformis is the  results of adaptation of higher latitude condition in north China [42], while the much later and lower bud set rate of P. yunnanensis is the results of adaptation of lower latitude condition in south China [43]. The phenology of bud set in P. densata showed adaptation of a lower latitude and higher altitude condition of the Tibetan Plateau. The higher BSR of P. densata than P. yunnannesis will provide better cold resistance, and the later bud set than P. tabuliformis will allow for more shoot growth. P. densata also showed the highest rate of the secondary growth, which suggests that most plants in this species have the latent growth in years with favorable climate conditions. This growth strategy is a benefit for the species to gain a competitive advantage in the hybrid niche.
The changes of seedling color reflect the status of plants responding to lower temperature in autumn [13]. Yellowing of seedling needles is suggestive of plant death, while needles changing color from green to red or purple may be a protective mechanism of plants against environmental stresses from ultraviolet radiation and cold [44][45][46][47][48][49][50]. For these three pine species, P. densata had higher value of RCL in purple, and P. yunnanensis had higher value of RCL in red, suggesting that the two species could make an effective response to strong ultraviolet rays. A higher yellowing rate was found in P. tabuliformis seedlings, reflecting their maladaptation to high altitude habitat. All lines of evidence from this study provided support to a better fitness of P. densata than parental species on the Tibetan Plateau through utilizing a longer growth season, better survival, moderate bud set time and the protection from physiological mechanism in in hybrid habitat. The fitness advantage of P.densata likely has contributed to the hybrid speciation.

Adaptation to different climate factors in the species complex
We included 16 representative populations of P. tabuliformis and P. yunnanensis in the field testing to provide a comprehensive assessment to the adaption of parental species in hybrid habitat. There were significant differences in precipitation, heat, light and other climate factors between the habitats of P. densata and the parental species. Low winter temperatures can cause chilling injuries for species which typically grows in the warmer and humid niches in the subtropical zone such as P. yunnanensis [43]. For P. yunnanensis, water depletion and hypothermy may have increased the mortality of its seedlings on the plateau. P. tabuliformis suffers selection pressure from a variety of ecological factors [40]. The growth time of P. tabuliformis should be determined by bud set phase and secondary growth rate [41]. We showed that bud set rate and secondary growth rate were significantly correlated with latitude (S1 Table) and the climate factors ( Table 4), suggesting that the climate factors could influence the adaptability of P. tabuliformis by affecting its bud set process. Blue and violet light on the Tibetan Plateau have strong inhibitory effects on plant growth [49,[51][52][53][54]. Purple rates of P. tabuliformis and P. yunnanensis seedlings were significantly correlated with altitude, indicating that the changes of seedling colors might affect the development of the parental species in the Tibetan Plateau.
A few populations of parent species (e.g. populations 16,17,18 of P. tabuliformis and populations 3,4,6,7 of P. yunnanensis) performed generally well in hybrid habitat, suggesting the fitness of each population may be related to its geographical distribution. The difference in fitness among population within parent species might come from (1) local adaption in parental habitat. Populations 16, 17, and 18, which are from the south distribution region, grew much faster than the other populations of P. tabuliformis. Population located in lower latitude of P. tabuliformis set bud at a later time than other populations, which maybe lengthen the growing period of the shoot on Tibetan Plateau. Similarly, P. yunnanensis is typically distributed in the subtropical zone, but the habitats of populations 3, 4, 6, and 7 were colder and drier than the other populations. Seedlings fitness of these populations was relatively high, which may be due to the cold tolerance of these populations. (2) Interspecific gene flow with P. densata. In our test, the GR, BSR, and SR in population 16 of P. tabuliformis were similar to the eastern population of P. densata, while the GR, BSR, and RCL in population 7 in P. yunnanensis were similar to central population of P. densata. Reproductive barrier between these species is weak [55]. Population 16 is located very close to the northern range of P. densata, bidirectional gene flow between the P. tabuliformis and P. densata may be possible through pollen exchange [10,56,57]. For P. yunnanensis, the genetic background of population 7 was similar to that of P. densata, which may have resulted from introgression and gene flow [10,23].
In general, adaptation to different climate factors may determine the adaptive differentiation at seedling stages of P. tabuliformis and P. yunnanensis on Tibetan Plateau. Both P. tabuliformis and P. yunnanensis occupy vast areas and different habitats from the hybrid. When transplanted out of the habitat, ecological factors could set strong selections on the fitness of each species. Due to the wide distribution of each species, differentiation in fitness traits within species is expected.

Differentiation in fitness traits among P. densata populations
We included 8 populations of P. densata in the field trial. We found adaptive differences among these geographic populations which reflected an altitude gradient of selection. The results from clustering analysis showed that populations from the west, the middle and the east of P. densata could be divided into three distinct groups, suggesting that the populations in each geographical region had different evolutionary process, corroborating with the molecular investigations on the evolutionary history of P. densata on the plateau [12]. Moreover, the east groups adapted to the western habitat generally well even though it is more than 1000 km away from its geographic origin. However, populations 11 and 12, from the central range of P. densata adapted to the western habitat poorly. In the field testing, continuous variation in SH and D0 was found associated with altitude. Similar altitude effects on the evolution of other species have been reported for the eastern Tibetan Plateau [58][59][60]. We also noticed similar performance between adjacent populations of different species, which might be explained by pollenmediated introgression [12; 16]. In our test the traits GR, BSR, and SR in eastern P. densata were similar to P. tabuliformis, while GR, BSR, and RCL in central P. densata were similar to P. yunnanensis. Molecular investigations show bidirectional gene flow between different species in the above regions [12,15,16]. Bidirectional gene flow between the hybrid species and either parent species may be occur via pollen flow when there is only weak crossability barrier among species [22,54,61]. The performance of the introgressed populations will be affected by location adaption as well as the gene flow.

Conclusions
As one of the most documented examples of a homoploid hybrid speciation, P. densata evolved a unique adaptation in the high plateau environment. The findings from the present study support the hypothesis of previous authors that habitat divergence promoted the speciation of P. densata on the Tibetan Plateau. Although results of our experiments can only partially identify the adaptive features of the three pine species placed in the habitat of P. densata, our data supports the general findings from molecular studies. Homoploid hybrid speciation can be driven by chromosomal and recombinational variation and ecological selection. Field test on fitness related traits provides direct evidence complementing the theory of homoploid hybrid speciation. Our present study together with more long-term field experiment will be helpful for understanding of speciation process of P. densata.
Supporting Information S1 Table. Correlation coefficient between measured traits and altitude for P. yunnanensis, P. densata and P. tabuliformis. (DOCX) S2 Table. Correlation coefficient between geographical variables and the climate factors. (DOCX)