Consistent with the genetic data, niche analysis of the first two axes of the environmental PCA and both overlap indices showed that the eastern and western clusters had a higher degree of niche overlap with each other than with HC (Fig. S7). HC was characterized by niches with higher mean annual temperature and lower precipitation in spring compared to those of the other clusters (Fig. S7). The equivalency test, background test, and ANOVA demonstrated that the three clusters occupied significantly different ecological niches (Fig. S7; Fig. S8).

 

The PCA and RDA of phenotypic traits were also congruent with the genetic analyses and showed that HC was morphologically and ecologically differentiated from the other clusters (Fig. 6a-b; Fig. S9; Fig. S10), which was classified as Q. longinux var. kuoi. Although the eastern and southern clusters had a high degree of overlap according to the first two PC axes, PERMANOVA indicated significant differentiation between the two clusters (Table S9). Partial RDA using all phenotypic traits as responses demonstrated that pure geography contributed more variation than pure environment, and a large intersection between their interactions was found (Fig. 6c-d; Table S10). The GLMs revealed that leaf traits were significantly associated with environment, but the directions of the relationships (i.e., positive or negative) differed depending on the trait and environmental variable (Fig. S11). For example, leaf thickness was positively correlated with annual temperature (r = 0.22, p < 0.01), whereas leaf length was negatively correlated with annual precipitation (r = −0.22, p < 0.01).

The high AUC value (average AUC = 0.81) suggested a good model fit for the predicted distribution of Q. longinux (Table S11; Fig. S12). As the predictions under different RCPs were highly correlated (Table S12), we inferred RONA based on the average values from both predictions. Substantial variations in the RONA estimates between populations and the top three representative environmental variables were observed (Fig. 7a-c; Table S10). The RONA estimates were larger in regions with greater differences between current and future climates (Fig. 7a-c). Whereas the eastern and western populations had relatively low RONA values (< 0.2) for the three variables, the northern populations were predicted to suffer from severe winter rainfall (precipitation in October, Fig. 7a) in the future and had much higher RONA values (> 0.6). 

 

The GF model constructed with five climatic variables suggested that precipitation in winter was the most influential variable with the highest weighted R² (Fig. S13). The focal species exhibited strong spatial patterns, indicating adaptation to local climate conditions (Fig. 7d). Consistent with the results of the RONA analysis, the GF model estimated highly correlated results between RCPs with similar genetic offset patterns (Fig. 7e-f). Both the RONA analysis and the GF model indicated that the northern populations were the most vulnerable to future climate change (Fig. 7e-f).