Does Cathaya argyrophylla, an ancient and threatened Pinaceae species endemic to China, show eco-physiological outliers to its Pinaceae relatives?

Abstract Cathaya argyrophylla is an ancient and threatened Pinaceae species endemic to China, but its eco-physiological traits are rarely reported. We hypothesized that Cathaya showed eco-physiological outliers to its Pinaceae relatives, which lead to its current endangered status. Here we collected the photosynthetic capacity (Pn, maximum photosynthesis rate) and branchlet hydraulic safety (P50, the water potential at which a 50% loss in conductivity occurs) of Pinaceae species globally, including our measurements on Cathaya. We applied the phylogenetic comparative methods to investigate: (i) the phylogenetic signal of the two key traits across Pinaceae species, and (ii) the trait–climate relationships and the photosynthesis–cavitation resistance relationship across Pinaceae species. We applied the polygenetic quantile regression (PQR) method to assess whether Cathaya showed eco-physiological outliers to its Pinaceae relatives in terms of cavitation resistance and photosynthetic capacity. It was found that P50, and to a less extent, Pn, had a strong phylogenetic signal consistent with niche conservation among Pinaceae species. Hydraulic safety largely determined non-threatened Pinaceae species’ distribution across moisture gradients at the global scale. There was also an adaptive trade-off relationship between Pn and P50. Cathaya is a high cavitation resistant, low photosynthetic capacity species. It showed eco-physiological outliers to its Pinaceae relatives because it had lower P50 and Pn below the 10% quantile boundaries along moisture and/or temperature gradients; also, it was above the 90% quantile boundary of the Pn and P50 relationship across non-endangered Pinaceae species. The PQR output demonstrated that in the subtropical area of China characterized by abundant rainfall, Cathaya has extra high hydraulic safety, suggesting inefficiency of carbon economy associated with either competition or other life history strategies, which lead to its current endangered status.


Introduction
The 'giant panda' of the plant kingdom, Cathaya argyrophylla Chun et Kuang, was formally described in 1958 by Chun and Kuang (Chun, 1958;Ying et al., 1983). It is a monotypic genus of Pinaceae endemic to China and has been categorized as a paleoendemic species, with a fossil history dating at least to the Cretaceous (Ferguson et al., 1997;Liu et al., 2000). Cathaya is one of the most vulnerable conifers in the world with total number of mature individuals less than 5000 (Xie et al., 1999b). Having its first occurrence in North America and being a formerly widespread Asiatic species, Cathaya has now become severely restricted to only three sites in subtropical mountains of China as a consequence of Late Tertiary climatic deterioration and Quaternary glaciation (Sup. Fig. 1) (Liu and Basinger, 2000). Extant Cathaya is found in a few mesothermal, physiological arid (high water deficit), infertile and exposed slopes in this area (Taggart et al., 2009;Xie et al., 1999a). It was also found that the abilities of competition and recolonization of this species are very low and that the habitats suitable for this species are undergoing rapid deterioration and fragmentation (Fu, 1992;Xie and Chen, 1999); as such, identification of abiotic and biotic factors contributing to its endangered status is vital for its conservation purposes.
How has Cathaya been vulnerable to extinction? Probably, such limited adaptability to regional climates and environments are due to its unique eco-physiological traits (Hunter Jr et al., 2006;Lusk, 2008). A species' uniqueness can be shown in three ways: (i) outlier to certain traits of its relatives (Hunter Jr and Gibbs, 2006); (ii) outlier to trait-environment relationships among its relatives (Newsome et al., 2020); and (iii) outlier to trait-trait relationships among its relatives (Brodribb et al., 2008). In the first case, blind, unpigmented, cave-dwelling invertebrates, fishes and amphibians are good examples of special characteristics in response to a rare type of habitat (Culver et al., 2000). In the second case, a good example is vertebrate species near the lower boundary of the relationship between habitat ranges and body sizes were more often in higher extinction risk categories (Newsome et al., 2020). In the third case, a good example is the tropic unique flat-leaved pine, Pinus krempfii, showing departure from the trade-off pattern between leaf hydraulic conductivity and mesophyll path length across pine species. This ecophysiological outlier probably is a key factor attributing to its limited adaptability and endangered status (Brodribb et al., 2008).
Although there is a large body of knowledge on the palaeoclimate, palaeogeology, anatomy, ecology, population genetics and evolutionary history of this species (Fu, 1992;Hu et al., 1984;Liu and Basinger, 2000;Qian et al., 2016;Ran et al., 2018;Xie et al., 1999b), the eco-physiological traits of Cathaya are rarely reported Zhang et al., 2005). In the present study, we chose two key functional traits to test whether Cathaya shows eco-physiological outliers to its Pinaceae relatives. One is the drought tolerance assessed by the water potential at which a 50% loss in hydraulic conductivity occurs [P 50 , MPa (Choat et al., 2012)]; another is the competition capacity using the maximum photosynthetic capacity [P n , μ mol CO 2 m −2 s −1 , (Augusto et al., 2014)] as a surrogate. Theoretically, low P 50 (high drought tolerance) means high carbon investment on xylem (e.g. high wood density, small tracheid size). The synthesis of xylem tissue is costly [6.5 and 11.8 mmol glucose per g of cellulose and lignin, respectively (Lambers et al., 1992)]; as a result, less carbon investment to leaves and inefficiency of carbon economy (lower P n ) are expected (Augusto et al., 2014;Pittermann et al., 2012). We thereby hypothesize that there exists a trade-off between P 50 and P n among non-threatened Pinaceae species, which has already been observed for the Cupressaeae species (Pittermann et al., 2012).
It is well known that moisture gradient largely determines the distribution of species with different drought tolerances (Choat et al., 2018;Choat et al., 2006;Willson et al., 2006). In the subtropical area of China characterized by abundant rainfall, common species should evolve trait suits with low drought tolerance. Therefore, it is possible that the restricted range of Cathaya is due to its extra-high drought tolerance, which would deviate from the general relationship between species drought tolerance and regional rainfall/water deficit gradient. Further, the high drought tolerance of Cathaya possibly comes at a significant cost on competition capacity, if there exists a trade-off between P 50 and P n among nonthreatened Pinaceae species. As such, low competition capacity would make Cathaya hard to overcome local species' competitive advantages (ecological barrier) for its survival and dispersion. However, other possibilities cannot be excluded; for example, Cathaya could deviate from the P 50 -P n tradeoff. Such departure means the participation of some other life history trade-offs and/or the regulation of some unusual structure characteristics (Brodribb et al., 2008).
The supposed unusual suites of traits of Cathaya, either due to evolutionary conservatism (e.g. more ancient Pinaceae species have higher drought tolerance) or due to developmental constraints (e.g. the rarity of major mutations) (Brodribb et al., 2008;Hunter Jr and Gibbs, 2006;Lusk, 2008;Wiens et al., 2005), could be evaluated preliminarily by phylogenetic niche conservatism (PNC). PNC is the tendency of lineages to retain their niche-related traits through speciation events (Harvey et al., 1991;Ricklefs, 2010). Although it is a pattern, not a process, the initial general test can lead to specific tests addressing the hypothesized drivers (Crisp et al., 2012). Additionally, as the cross-species correlations may be biased by a lack of statistical independence among closely related species with shared evolutionary history, phylogenetic signals are also needed to be known a priori, in order to minimize the Type I error rates of non-phylogenetically corrected crossspecies correlations (Ackerly, 2000;Blomberg et al., 2003).
In the present study, we addressed three questions: to detect (i) the phylogenetic signal of P 50 and P n across Pinaceae species, based on a global database and our measurements (ii) which environmental variables drive the distribution of P 50 and P n across common Pinaceae species and if there exists a trade-off between P 50 and P n ; and (iii) if Cathaya shows outliers to the general relationships among common Pinaceae species between P 50 /P n and environmental variables and the P 50 vs. P n trade-off relationship.

Materials and methods
Field surveys were carried out at Jiaopengliao, Zixing City, Hunan Province, South China, with a geographic location of 25 • 56 54 N, 113 • 29 24 E and an altitude of 1250 m, on the southwest part of Bamian Mountain Reserve. The forest was a natural mixed stand of conifer and broadleaved trees. Besides Cathaya, there were also Castanopsis and Altingia tree species in the main canopy storey. The climate is a subtropical monsoon climate area with four distinct seasons and mild climate; the meteorological data could be seen in Supplementary Table 1. The main soil type is mountain yellow-brown earth, with a total nitrogen content of only 0.0061% . A total of 26 sunlit branches from six mature trees were selected for cavitation resistance measurement, and leaves from one branch from the six mature trees were selected for photosynthesis measurement.

Light-response curve measurement
The light-response curve was measured from 9:00 to 12:00 AM in July with an infrared gas analyzer (LI-6400 Li-Cor, Lincoln, NE, USA) performed using the autoprogram function. One-year-old leaves were chosen for measurements. We chose the sequence of desired light settings of 1500; 1200; 800, 500, 300, 100, 50, 30 and 0 μmol (photon) m −2 s −1 , a minimum wait time of 120 s and a maximum wait time of 200 s, and matching the infrared gas analyzers for the 50-μmol(CO 2 ) mol(air) −1 difference in the CO 2 concentration between the sample and the reference, which allowed them to be matched before every change in irradiance. Leaf temperature was kept between 24 and 26 • C, and CO 2 partial pressure of and relative humidity in the leaf cuvette was 60-80% during measurements. Leaf mass ratio (LMA) measurements were then taken on P n -measured specimens in the laboratory. Ten to 15 fully expanded healthy leaves, not including petioles, were selected (Pérez-Harguindeguy et al., 2013). The stomatal conductance (g s ) and transpiration rate (T r ) at saturating irradiance was also recorded.
The light-response curve can be fitted by a non-rectangular hyperbola (Lambers et al., 2008): where A is the net photosynthetic rate under specific irradiance; P n is the light-saturated net photosynthetic rate; I is the irradiance; φ is the quantum yield of assimilation, which is the initial slope for the light-response curve; is the convexity of the curve; and R d is the dark respiration rate. φ and R d were calculated using linear regression analysis (A to irradiance <150 μmol m −2 s −1 ), then P n and were derived empirically by fitting (1) to light response data. The light compensation point (LCP) was the irradiance when A = 0, below which there is insufficient light to compensate for respiratory carbon loss.

Cavitation vulnerability curve and tracheid diameter measurements
Three to five branchlets (40-70 cm long, 2-4 years old) near that used for the light-response measurement were cut from each sample tree before sunrise (before 8:00 AM) to measure the specific conductivity (k s ; kg m −1 MPa −1 s −1 ), leaf specific conductivity (k l ; kg m −1 MPa −1 s −1 ) and cavitation vulnerability curve (totally 26 samples). The branchlet was wrapped in moist paper towels immediately after being cut and transported to the laboratory. All the samples were recut to 22-24 cm underwater, and all the measurements were carried out in an air-conditioned laboratory (26 • C). The maximum flow rate was measured under 8 kPa hydrostatic pressure after air emboli were flushed out by perfusion with 110 kPa distilled water (flowing through 0.2 μm filter) for 30 min. Measurements were initiated after ∼ 2 min when the flow rate stabilized. The weight of the collected efflux was measured every 30 s with a precision balance (Sartorius, BP221S, Göttingen, Germany) to obtain the flow rate. k l was calculated by dividing the maximum flow rate by the total leaf area distal to the measured segment and by the pressure gradient. The leaf area was determined using a Win-FOLIA system (Regent Instruments, Quebec City, Canada). k s was calculated by dividing the maximum flow rate by the segment's cross-section sapwood area. The total acropetalend cross-section area of the branch segment was determined from its maximum and minimum diameters. The area of the pith was determined from its dimensions measured under a dissecting microscope equipped with a stage micrometer and subtracted from the above acropetal-end cross-section area to determine the cross-section sapwood area. The Huber value was calculated as the total cross-section sapwood area per unit leaf area (Tyree et al., 1988).
The vulnerability of xylem to cavitation was characterized using a vulnerability curve, which was measured using a cavitation pressure chamber (PMS Instrument, Corvallis, Oregon, USA) according to Sperry and Saliendra (1994). A branch segment was inserted into a collar and sealed with both ends protruding. Air was injected into the collar at a set pressure, which was maintained for 15 min and then slowly decreased to 0.1 MPa. The hydraulic conductivity was then re-measured at a higher pressure. This procedure was repeated until at least 85% loss of hydraulic conductivity was reached. pressurization was calculated as PLC = 100 × (K h -K hi )/K h , where K hi is the hydraulic conductivity measured at pressure i. The vulnerability curve for each sample was fitted with an exponential sigmoidal equation (Pammenter et al., 1998): where Ψ is the negative of the injection pressure and a and b are coefficients estimated using a non-linear regression model with gnls() function of R software (R Development Core Team, 2016). The coefficient b represents P 50 , a represents the steepness of vulnerability curve. Wood density (W) was measured on stem segments used in the measurement of vulnerability curves after the removal of pith and bark, and the fresh volume was measured by the Archimedes principle of water displacement. The dry mass was determined after drying at 104 • C for 24 h. W is expressed as dry mass per unit of fresh volume (g cm −3 ).
Three to five samples out of those used for VC curve construction were used for the anatomical measurements. Micro thin sections (12 μm thick) were taken, following the procedure descried in Domec et al. (2010). Those sections were stained with Safranin O and mounted in resin for image analysis using ImageJ (Media Cybernetics, Silver Spring, MD). The mean hydraulic diameter (Kolb et al., 1999)] was estimated. The equivalent circle diameter was taken as tracheid diameter based on the lumen area measured . For each sample on pit-cambium transects separately, three radial files of tracheids were randomly selected, and all lumen areas were measured.

The global datasets, tree phylogenies and phylogenetic comparative analyses
We collected three datasets in the present study: P 50 dataset, P n dataset and P 50 vs. d h dataset. The branchlet P 50 dataset used in the present study can be accessed from the TRY Plant Traits Database (https://www.try-db.org/TryWeb/Home. php) (Choat et al., 2012). We found 35 Pinaceae species (in addition to Cathaya) from the database with definite locations on the maximally resolved seed plant tree (Zanne et al., 2014). Based on the P 50 dataset, we selected the same Pinaceae species for the P n dataset from the TRY Plant Traits Database. P n with the same (or slightly different) latitude and longitude of the P 50 dataset, measured under growth temperature, ambient CO 2 and saturating irradiance, with VPD (vapour pressure deficit) ranging from 0.8 to 1.5 kPa and/or relative humidity from 55 to 80%, of sunlit leaves from mature trees, were chosen out of more than 14 000 records. P n values of four Pinaceae species (Abies concolor, Cedrus deodara, Pinus caribaea, Pinus echinata) could not be found based on the selection criteria, and we added another two Pinaceae species' data (Pinus krempfii, Pinus kesiya) from Brodribb et al. (2008), resulting in 34 Pinaceae species in the P n dataset (including Cathaya). In the P 50 dataset, there were 18 Pinaceae species with P n measured at the same location (or a slightly different location with an elevation difference of <300 m); this sub-dataset was analyzed for the relationship between P 50 and P n . At last, we collected eight literatures (Davis et al., 1999;Domec et al., 2010;Hacke and Jansen 2009;Linton et al., 1998;Maherali et al., 2006;Mayr et al., 2006;Oliveras et al., 2003;Sterck et al., 2012), in which both d h and P 50 of 14 Pinaceae species have been measured simultaneously, to construct the P 50 vs. d h dataset. Six environmental variables [latitude, altitude, mean annual temperature (MAT), mean annual precipitation (MAP), mean temperature in the coldest month (MTCM) and water deficit (WD)] were selected for trait-environment analysis. Climate data (means over 1950-2000) were extracted from the WorldClim2 database with a spatial resolution of ca. 1 km (Fick et al., 2017) based on the geographic coordinates of sampling sites.
We then constructed the phylogenies for our species with the program Phylomatic (Webb et al., 2008), using a maximally resolved seed plant tree (Zanne et al., 2014). Phylogenies were constructed for the P 50 , P n , P 50 vs. P n and P 50 vs. d h dataset (Supplementary Figs. 2 and 3 for P 50 and P n ). Afterwards, multiple approaches were applied to evaluate phylogenetic signals in P 50 and P n . The first approach is Pagel's λ calculation (Freckleton et al., 2002) (using the 'Phytools' package in R). Secondly, the phylogenetic signal representation curve (PSR) approach has been adopted (Felizola Diniz Filho et al., 2012;Hawkins et al., 2014). The PSR area, expressing deviations from Brownian motion across the curve, is the same as Blomberg's K statistic and can reveal whether traits are evolving at a slower (negative) or faster (positive) rate than expected under Brownian motion in different parts of the phylogeny (Felizola Diniz Filho et al., 2012). PSR curves for P 50 and P n were calculated by using the "PVR" package in R. The third approach (phylogenetic generalized least squares) followed Hawkins et al. (2014), Kozak et al. (2010) and Wiens et al. (2010). Three models of evolution for P n and P 50 were applied: a white noise (WN) model of random variation; a Brownian motion (BM) model of gradual and continuous drift in trait; and an Ornstein-Uhlenbeck (OU) model of constrained evolution. OU-constrained evolution (stabilizing selection) was regarded as the sign of phylogenetic signal to a more stringent sense (Crisp and Cook, 2012). The log-likelihood of phylogenetic generalized least squares fitting of the three models was calculated. The "GEIGER" package in R was used to calculate and compare the fit of each model using a likelihood test; the Akaike information criterion (AIC) derived from model likelihoods was used to determine whether each trait fits an OU model better than either a BM or WN model.  branch length transformations optimized for the maximum likelihood was applied, using the 'caper' package in R. In PGLS analysis, regression parameters are found by maximum likelihood (ML) and 'weighted' by the variance-covariance matrix that represents the phylogenetic relationships among the species. The relationship between P n and P 50 , d h and P 50 were also explored by PGLS. Assume that the pressure that causes the cavitation of water inside xylem tracheids is inversely proportional to the mean of the maximum (per tracheid) size of pit pores (d p ) (air-seeding hypothesis (Zimmermann, 1983)) P 50 ∝1/d p , and d h is linear to the size of its largest pore: d h ∝d p (Martínez-Vilalta et al., 2002), resulting in P 50 ∝1/d h . Therefore, for the P 50 vs. d h relationship investigation, − P 50 and d h were log-transformed beforehand.
To investigate whether Cathaya shows eco-physiological outliers to its Pinaceae relatives, the polygenetic quantile regression (PQR) approach was applied. Since P 50 /P n is likely determined by multiple unmeasured factors, and species with trait values outside the quantile regression lines are more likely to be threatened with extinction, it is useful to define an approximate upper/lower boundary for P 50 /P n at a given variable, rather than focusing on means (Hacke et al., 2017;Newsome et al., 2020). To conduct PQR across non-threatened species, we followed the method by Newsome et al. (2020). We firstly fitted the non-phylogenetic quantile regression model and obtained residuals; secondly, we constructed a phylogenetic "distance" matrix from the phylogenetic tree (one minus the correlation matrix); thirdly, we calculated the phylogenetic autocovariate using the inverse distance squared weighted mean of residuals; then, the phylogenetic autocovariate was included in the quantile regression model (we set the auto-covariate to its average value, thus the boundaries that we show are for species with average values of the residual auto-covariate). We fitted quantile regression models for the 90% and 10% quantiles to describe the upper and lower bounds, respectively. The PQR was carried out by using the 'ape' and 'quantreg' packages in R (Koenker et al., 2018).

Results
The P n and branchlet P 50 of Cathaya were 5.96 μmol m −2 s −1 (Table 1; Fig. 1A) and −5.64 MPa (Table 1; Fig. 1B), respectively. P n was in the lowest 25% quartile range among compiled Pinaceae species [Sup. Table 1, compiled P n data of Pinaceae species can also be seen in Brodribb and Feild (2008)]. Besides two Cedrus species (Cedrus libani, Cedrus brevifolia), Cathaya was the most cavitation resistant one among compiled Pinaceae species (Sup. Table 1) [compiled P 50 data of Pinaceae species can also be seen in Delucia et al. (2000) and Martínez-Vilalta et al. (2004)]. High resistance to cavitation indicated low hydraulic conductance (k s and k l ), as well as narrow tracheids (d and d h ). In fact, the tracheid diameter of Cathaya was the narrowest one among the investigated Pinaceae species (Sup. Table 2; Sup. Fig. 3) [see also Hacke et al. (2017) and Martínez-Vilalta et al. (2012)]. The inefficiency of the carbon economy could also be indicated by low g s , T r and R d and high LMA ( Table 1). The LCP (light compensation point) was 12.41 μmol m −2 s −1 ; however, we could not obtain a reliable LSP (light saturation point), as with the increase in irradiance, P n increased slightly even at a very strong irradiance (>1500 μmol m −2 s −1 ); furthermore, there is no consensus on the accuracy of LSP estimation by a non-rectangular hyperbola (Li et al., 2019).
We calculated the phylogenetic signal of P 50 and P n based on global datasets of Pinaceae species; Pagel's λ of P 50 and Pn were 0.99 and 0.53, respectively (Table 2). P 50 was different from the null hypothesis (no phylogenetic dependence), while P n was marginally different from the null hypothesis,  suggesting P 50 was more phylogenetically conservative than P n . This could further be verified by the PSR area (Fig. 2); the PSR area for P 50 was −0.16 while that for P n was −0.11 (more negative area means more conservative). The fittings by phylogenetic generalized least squares showed that the Ornstein-Uhlenbeck (OU) model was the best model for comparative fits of P 50 and P n ( Table 2), suggesting that both traits satisfied the most stringent definition of phylogenetic niche conservatism (Crisp and Cook, 2012;Hawkins et al., 2014;Losos, 2008). Table 3 shows the PGLS slopes of P 50 and P n with eight environmental variables (including the first two axes of principal component analysis on the six environmental variables, i.e. PC1 and PC2). For the non-threatened species, P 50 had significant associations with altitude, MAP, WD and PC2. There were no associations between P n and environmental variables across all the common species, and only PC1 showed a marginal association with P n (P = 0.091). PC1 was mostly represented by temperature variables (MAT and MTCM), while PC2 was mostly represented by moisture variables (water deficit index, negatively with MAP and positively with WD) (Table 3).
Among non-threatened species, P 50 was negatively related to altitude, WD and PC2 (water deficit index) and positively related to MAP (Fig. 3). If Cathaya and other five threatened species were positioned, Cathaya and two Cedrus (Cedrus libani, Cedrus brevifolia) showed significant deviations from the general trends: the three species had extraordinarily lower P 50 than PGLS predictions and were below the 10% quantile boundaries. Further, the four PGLS models predicted that P 50 of Cathaya should be −3.3∼−4.2 MPa, while the measured value was −5.64 MPa. P n was only related to PC1 (temperature index) (marginally, P = 0.091), indicating that higher photosynthetic capacity is associated with higher growth temperature (Fig. 4). Again Cathaya showed a sign of deviation from the general trend and was below the 10% quantile boundary. The PGLS model predicted that P n should be around 11 μ mol CO 2 m −2 s −1 , while the measured value was 5.96 μmol m −2 s −1 .
P 50 was positively related to P n (trade-off) (P = 0.02) among the 13 non-threatened Pinaceae species, when PGLS was applied (Fig. 5A). If Cathaya and other threatened species were positioned, Cathaya and the two Cedrus species deviated from the general trend and were above the 90% quantile boundary, which means they maintained considerable P n at such low P 50 levels. There existed no significant relationship between d h and P 50 among non-endangered species when PGLS was applied; despite that Cathaya had the narrowest  The first two axes of principal component analysis, PC1 and PC2, across all investigated species in P 50 and P n datasets were also demonstrated. MAT: mean annual temperature ( • C); MAP: mean annual precipitation (mm); MTCM: mean temperature in the coldest month ( • C); WD: water deficit (mm). Principal component analysis was conducted after the data have been Z-score standardized. Numbers in bold mean significant relationships between P 50 /P n and environmental variables of nonthreatened Pinaceae species revealed by PGLS.

Discussion
Question 1: the phylogenetic signal of P 50 and P n across Pinaceae species The hydraulic safety, which represents a realization of drought tolerance, had a strong phylogenetic signal consistent with niche conservation (Sup. Fig. 2A; Fig. 2A; Table 2). Moreover, to a less extent, the photosynthetic capacity, which represents a realization of carbon economy efficiency, also had a strong phylogenetic signal (Sup. Fig. 2B; Fig. 2B; Table 2). This could be validated by Pagel's λ, PSR and model fitting comparison results. As such, not only morphological traits (Nobis et al., 2012) but also the key physiological functions of hydraulic safety and photosynthetic capacity showed phylogenetic niche conservation (PNC) signals among Pinaceae species. On the other hand, Wang et al. (2018) found that P n showed more PNC than P 50 among nine Picea species in mid-China, which is opposite to our finding. This is probably due to the 'local scale effect' (Forrestel et al., 2017) or statistical 'funnel effect' caused by the insufficient species amount investigated (Blomberg et al., 2003;Wright et al., 2005). The strong PNC of P 50 and/or P n as revealed by the present study is consistent with the observation on other gymnosperms (Brodribb et al., 2010;Pittermann et al., 2010;Willson et al., 2008). Ancient Pinaceae species tend to have high drought tolerance (low P 50 ) and low photosynthetic capacity (low P n ), in accordance with their initial diversification in the Permian, which was characterized by physiological aridity, and consequently they have a number of traits that provide physiological drought tolerance, permitting them to survive on frozen soils, deep sand and steep slopes (Graham, 1999).

Question 2: the trait-environment relationship and the trade-off between P 50 and P n
There existed significant PGLS relationships (slope) between altitude (negative), WD (negative) and P 50 among nonthreatened species. The results demonstrated that moisture gradient (WD and/or MAP) largely drives the distribution of species with different cavitation resistance (Table 3; Fig. 3). This could be further demonstrated from the significant PGLS negative slope of the second component (PC2) vs. P 50 relationship, while PC2 was largely determined by WD and MAP. High altitude, high WD and low MAP are conducive to high cavitation resistance (low P 50 ) among nonthreatened species. This output supports other studies that hydraulic safety largely determines species distribution across moisture gradients (Choat et al., 2018;Choat et al., 2012;Niinemets and Valladares, 2006;Willson and Jackson, 2006), although most of which did not pay attention to phylogenetic dependence of inter-specific data. On the contrary, there existed very weak or no relationships between environmental variables and photosynthetic capacity among non-threatened species (Table 3), suggesting that hydraulic safety is more adaptive than photosynthetic capacity across Pinaceae species. Only temperature (PC1) executed a marginal impact on P n (Fig. 4), which indicated the temperature dependence of enzymes' kinetics responsible for carbon assimilation (Brooks et al., 1985;Yamori et al., 2005).
It has been demonstrated that the hydraulic efficiency of branch has a tight association with leaf photosynthetic capacity (Brodribb et al., 2002;Scoffoni et al., 2016;Xiong et al., 2019). If there exists an adaptive trade-off relation between hydraulic efficiency and cavitation resistance among Pinaceae species (Choat et al., 2011;Gleason et al., 2016;Maherali et al., 2004), it means photosynthesis could have a negative association with cavitation resistance among the non-threatened species globally. As expected by either the carbon economy theory (Augusto et al., 2014) or k l -P n coordination (Brodribb et al., 2002;Scoffoni et al., 2016;Xiong et al., 2019), the PGLS output of the 13 species did support the P 50 -P n trade-off (Fig. 5A, P = 0.02), which was consistent with some other studies on gymnosperm species (Maherali et al., 2006;Pittermann et al., 2012;Sterck et   al., 2012). On the contrary, the non-significant association between d h and P 50 among the 14 Pinaceae species (Fig. 5B) did not support the theoretical model (Hacke et al., 2004, Martínez-Vilalta et al., 2002Mayr et al., 2006), which means that the torus-margo structure of gymnosperms (Choat et al., 2009;Domec et al., 2009;Hacke et al., 2009) possibly masks the tracheid allometry deviated from the theoretical value based on relative simpler pit structure of angiosperms.
Question 3: does Cathaya deviate from the general relationship among non-threatened Pinaceae species between P 50 /P n and environmental variables?
Cathaya did show eco-physiological outliers of P 50 /P n to its Pinaceae relatives. Besides two Cedrus species (Cedrus libani, Cedrus brevifolia), Cathaya was the most cavitation-resistant Pinaceae species among the Pinaceae species reported. P n was in the lowest 25% quartile range among compiled Pinaceae species (Sup. Table 1). Cathaya also had the narrowest tracheid among the investigated Pinaceae species (Sup. Table 2; Fig. 5B). More importantly, Cathaya had P 50 under the low (10%) quantile boundaries of P 50 vs. environmental variables (altitude, MAP, WD and moisture index PC2) (Fig. 3) and P n under the 10% quantile boundary of P n along the temperature gradient (temperature index PC1, Fig. 4). The extra high hydraulic safety of Cathaya and low photosynthetic capacity are probably due to its ancientness (Sup. Fig. 2 A, B, Table 2, Fig. 2), not developmental constraint. Interestingly, unlike Cathaya, the extra high hydraulic safety of Cedrus libani and Cedrus brevifolia could not be explained by their recent evolutionary history and may largely be attributed by developmental constraint (Sup. Fig. 2 A, B). High cavitation resistance is associated with inefficiency of carbon economy, leading to Cathaya's lower growth competition ability with its coexisting angiosperm species (Augusto et al., 2014;Brodribb et al., 2007;Boyce et al., 2009;Körner, 1995;Lusk et al., 2003), and can only be confined to sites with arid (high water deficit), infertile and steep slopes in the subtropical mountains of South China.
However, it is worth noting that Cathaya was above the 90% quantile boundary of the P n and P 50 trade-off [although not as significant as the two Cedrus species (Fig. 5A)], indicating that it could maintain considerable P n at such low P 50 levels. The departure indicates the participation of some other life history trade-offs. For example, extra high cavitation resistance could come at a significant energy cost in terms of insufficient reproductive effort and does not necessarily compromise the competition capacity of mature individuals. Actually, it was found that the average seed weight of Cathaya was 0.0197 g [only 6% of the average seed weight of arbour species (Silvertown, 1982)], with over 60% of seeds without viability (Xie et al., 2000), indicating insufficient energy input, although other mechanism(s) could not be overlooked (Hunter Jr and Gibbs, 2006;Vaario et al., 2006;Wang et al., 2006). The other mechanism(s) include specific arbuscular mycorrhizal root interaction (Vaario et al., 2006) and low gene flow between isolated populations (Wang et al., 2006). Another possibility is that some unusual structure characteristics such as torus thickness and depth of the pit chamber (Hacke and Jansen, 2009;Pittermann et al., 2005) and tracheid wall reinforcement (Hacke et al., 2017), but not tracheid diameter (Fig. 5B), drive the deviation. It is worth noting that Cathaya does not have tubular sclereids in the leaf (Dayong Fan, unpublished data), unlike the tropic unique flat-leaved pine, Pinus krempfii (Brodribb et al., 2008). The torus-pit structure characteristics probably also led to the absence of the trade-off between P 50 and k s within Cathaya (Sup. Fig. 4A), but not the P 50 -wood density relationship (Sup. Fig. 4B). As such, the mechanism(s) underlying the deviation of Cathaya (and the two Cedrus species) from the P n and P 50 trade-off merits future investigation. Interestingly, in the relationships between hydraulic safety and environment, three threatened species (Pinus caribaea, Pinus albicaulis, Abies pinsapo) were between the 10% and 90% quantile boundaries (Fig. 3), which were quite different from Cathaya and the two Cedrus species. The result suggests different mechanisms of endangerment for the six threatened species. For example, Cathaya and the two Cedrus were restricted to narrow range and habitats (Vidaković, 1991;Xie et al., 1999b), since the site-specific environmental conditions might be quite different from the regional environmental conditions (the Climate data were extracted based on regional scale), and their deviations from trait-environment relationships could be expected, while for the other three species, they have broader distribution areas (Farjon, 1990;Thompson et al., 1999), so probably other mechanisms [e.g. inter-specific interaction, human disturbance (Hunter Jr and Gibbs, 2006)] but not the limited adaptability (hydraulic safety) underpinned their vulnerability to extinction. As such, the PQR method on functional traits and their relationships to the environment can offer a useful tool to reveal the mechanism of species endangerment, in support of the viewpoint by Newsome et al. (2020). It is also worth noting that although the phylogenetic comparative methods applied in the present study could minimize the Type I error rates of non-phylogenetically corrected cross-species correlations (Ackerly, 2000;Blomberg et al., 2003), other factors such as non-random sampling, intra-specific variation and measurement error could impact the output (Ives et al., 2007;Huey et al., 2019). As such, further investigation on the difference of eco-physiological traits between Cathaya argyrophylla and its Pinaceae relatives is warranted.

Conclusions
In the present study, we found that the hydraulic safety, and to a less extent, photosynthetic capacity, had a strong phylogenetic signal consistent with niche conservation among Pinaceae species. Hydraulic safety largely determined the distribution of non-threatened Pinaceae species across moisture gradients at the global scale. There was a phylogenetic trade-off between hydraulic safety and photosynthesis across non-threatened Pinaceae species. Cathaya is a shadeintolerant, high cavitation resistant, low photosynthetic capacity Pinaceae species. These traits were largely attributed by its ancientness. Cathaya also showed eco-physiological outliers to its Pinaceae relatives because it had lower P n and P 50 outside the 10% quantile boundaries along moisture and/or temperature gradients. The safety-photosynthesis trade-off could provide an explanation on how Cathaya is vulnerable to extinction, but other factors contributing to its endangered status could not be overlooked. The somewhat departure of Cathaya from the safety-photosynthesis tradeoff merits future investigation. We suggested that the PQR method could be a useful tool to reveal the mechanism of species endangerment of Pinaceae species under the global climate change.