Ecophysiological Vulnerability to Climate Change in Mexico City’s Urban Forest

Urban forests play an important role in regulating urban climate while providing multiple environmental services. These forests, however, are threatened by changes in climate, as plants are exposed not only to global climate change but also to urban climate, having an impact on physiological functions. Here, we selected two physiological variables (stomatal conductance and leaf water potential) and four environmental variables (air temperature, photosynthetically active radiation, vapor pressure deficit, and water availability) to compare and evaluate the ecophysiological vulnerability to climate change of 15 dominant tree species from Mexico City’s urban forest. The stomatal conductance response was evaluated using the boundary-line analysis, which allowed us to compare the stomatal response to changes in the environment among species. Our results showed differential species responses to the environmental variables and identified Buddleja cordata and Populus deltoides as the least and most vulnerable species, respectively. Air temperatures above 33°C and vapor pressure deficit above 3.5 kPa limited the stomatal function of all species. Stomatal conductance was more sensitive to changes in leaf water potential, followed by vapor pressure deficit, indicating that water is a key factor for tree species performance in Mexico City’s urban forest. Our findings can help to optimize species selection considering future climate change by identifying vulnerable and resilient species.


INTRODUCTION
Urbanization drastically changes the natural environment with significant land-use changes (Berry, 2008). Urban ecosystems are also exposed to pollutants, high CO 2 , nitrogen deposition, and elevated temperatures, which affect urban vegetation (Gregg et al., 2003, Zhao et al., 2016. Urban tree growth, for example, can be stimulated by warmer temperatures, particularly in temperate climates. However, high temperatures can also accelerate senescence and aging (Moser et al., 2017;Smith et al., 2019) and additional challenging conditions in cities, such as limited soil volume, modified soils, and reduced water and nutrient availability affect tree growth and performance (Day and Bassuk, 1994;Hauer et al., 2020).
The imposition of large concrete and asphalt slab restricts natural evapotranspiration in cities and considerably increases heat storage (Gui et al., 2007;Göbel et al., 2013). This redistribution of energy results in the urban heat island (UHI) effect. The UHI effect consists of an increase in the air temperature of the urban area with respect to the rural surroundings (Oke, 1995). The UHI in Mexico City, for example, can be as high as 7 • C on average, but in some instances can reach intensities of up to 10 • C (Ballinas, 2011).
Global climate change exacerbates urban temperature differences (Brazel and Quattrochi, 2005). In coming decades, heat shocks or waves will increase in frequency and severity and be enhanced by the UHI effect (Luber and McGeehin, 2008;Manoli et al., 2019). These predicted increases in air temperature around the world will affect the thermal balance with consequences on natural and urban ecosystems (IPCC, 2021). Undoubtedly, all organisms living in cities, including human populations, will experience more heat stress, which might decrease productivity and severely affect human health and mortality (Laschewski and Jendritzky, 2002).
With proper selection and strategic planning, urban vegetation can mitigate the UHI because of the cooling potential plants have via transpiration (Barradas, 1991(Barradas, , 2000Susca et al., 2011;Ballinas and Barradas, 2016a;Kornaska et al., 2016). Some trees, for example, can transpire up to 5 L day −1 tree −1 , which is equivalent to a heat dissipation of up to 12.62 W day −1 tree −1 (Ballinas and Barradas, 2016b). For this reason, urban vegetation has gained recognition for management of the UHI effect, along with other mechanisms that contribute to the improvement of the urban environment, such as the capture of atmospheric particulate material and improvement of human health (García-Sánchez et al., 2019;Keeler et al., 2019).
Plants, however, are exposed not only to global but also urban climate change, affecting key physiological functions. As a result, there are three possible plant responses to changes in climate: migration, adaptation, and extinction (Aitken et al., 2008). As cities are fragmented landscapes, species' movement is limited (Dobbs et al., 2017); therefore, adaptation (i.e., evolutionary changes and/or physiological acclimatization) becomes crucial for species survival (Atkins and Travis, 2010).
Stomatal conductance (g S ) is a key plant trait in photosynthesis and plant productivity (Jones, 1992;Jarvis and Davies, 1998). This trait also responds rapidly to climate and environmental changes (Jones, 1992;Ackerly, 2004;Buckley, 2005). Some climatic changes, such as elevated temperatures, can negatively affect plant functions by reducing transpiration, nutrient assimilation and photosynthesis (Farquhar and Sharkey, 1982;Atkin and Tjoelker, 2003;Mathur et al., 2014). Therefore, understanding how climatic and environmental changes affect species functional responses (e.g., g S response) becomes necessary in face of global and urban climate change.
Global climate change represents a threat to the persistence of urban forests globally (Esperon-Rodriguez et al., 2021) and these changes increase species' vulnerability-i.e., the extent to which climate change (urban and global) can damage a system depending on its sensitivity and its ability to adapt itself to new climatic conditions (IPCC, 2001). Further, when these climate changes affect how species respond physiologically to the environment, their ecophysiological vulnerability increases, limiting their ability to adapt physiological functions to environmental changes (Esperon-Rodriguez and Barradas, 2015a). Thus, identifying vulnerable species in urban environments can help to improve species selection while ensuring the provision of desired ecosystem services toward sustainable development.
Here, we aimed to determine and evaluate the ecophysiological vulnerability of Mexico City's urban forest to climate change. For this, we used g S as a vulnerability indicator and four environmental variables were selected to assess the species ecophysiological vulnerability: air temperature, photosynthetically active radiation, vapor pressure deficit, and water availability. The objectives of this study were to (1) identify more and less vulnerable species to environmental changes and (2) develop a model summarizing the g S response to the environment (i.e., the four aforementioned variables). We used the boundaryline analysis to predict stomatal responses. This analysis has great precision and uses meteorological data to assess ecophysiological responses to the environment (Jarvis, 1976;Lambers et al., 1998;Barradas et al., 2011;Ballinas and Barradas, 2016b) and can be used to identify vulnerable species (Esperon-Rodriguez and Barradas, 2015a,b).

Study Area
Data were collected in five parks of Mexico City (19 • 19 N, 99 • 11 W, 2230 m above sea level): España (EP), Mexico (MP), Luis G. Urbina (LGUP), Bombilla (BP), and Ciudad Universitaria (CU) (Figure 1). The city's mean annual rainfall (average of 40 years) is 748 mm, where ∼94% occurs during the rainy season (June-November). Winds are light and predominantly from the Northeast. Extreme average temperatures occur in April  Values represent averages and the standard deviation is given in brackets (N = 10).
(26 • C) and January (5.3 • C) and currently the annual average temperature ranges from 14.6 to 17.3 • C (CONAGUA et al., 2021), with a distinct urban heat island effect (Jauregui, 1973) not only at night but also during daytime, with differences of up to 10 • C (T U −R ) compared to surrounding areas (Ballinas, 2011; Figure 1). The increase on average in air temperature since 1990 has been 0.052 • C per year and predictions of temperature suggest an increase to 24.5 • C by 2050 (Ballinas, 2011).

Selected Species
Fifteen dominant tree species were selected based on their presence and abundance in five parks from Mexico City and included:

Data Collection
All individual trees were identified within urban parks located toward the center and west of the city (Figure 1). Stomatal conductance (g S ) was measured with a steady-state diffusion porometer (LI-1600, LI-COR, Lincoln, NE, United States) on at least five sunlit and shaded expanded leaves per individual tree in the mid-level of the tree and for at least four individuals at each park. Photosynthetically active radiation (PAR), air temperature (T A ) and relative humidity (RH) were determined next to each measured leaf with a quantum sensor (LI-190SB, LI-COR Ltd., Lincoln, Nebraska, United States), a thermocouple and a humicap sensor (Vaisala, Helsinki, Finland), respectively. The quantum sensor was installed near and parallel to the leaf maintaining the orientation of the leaf when measuring. The vapor pressure deficit (VPD) was calculated from T A and RH measurements: VPD = e S (1-RH) and e S = 0.6108{exp[(17.27T A )/(T A + 237. 3)]} and 0 ≤ RH ≤ 1. Concomitantly, leaf water potential (ψ) was measured in two fulllight exposed leaves per individual tree, with a pressure chamber (PMS Corvallis, Oregon, United States) (Scholander et al., 1964). Measurements were made from 08:00 to 16:00 h (local time, h) during 16 days in April and May, the warmest and driest months of the year. We recorded the highest (i.e., less negative) ψ between 08:00 and 09:00 h. During measurements, the maximum Ψ µ FIGURE 2 | Scatter diagrams and probable boundary-line function of stomatal conductance (g S ) vs. air temperature (T A ), leaf water potential (ψ), vapor pressure deficit (VPD), and photosynthetically active radiation (PAR) for Acacia longifolia.
average temperature was 29 • C and precipitation was 10.2 mm was registered only during the last 6 days of March before measurements. No rain was recorded during the study period.

The Boundary-Line Analysis
The boundary-line analysis assesses the relationship between g S and climatic/environmental and physiological variables (i.e., T A , PAR, VPD, ψ). From this analysis, graphs are generated representing the optimal stomatal response to a given selected variable (Jarvis, 1976;Fanjul and Barradas, 1985;Jones, 1992;Barradas et al., 2004;Esperon-Rodriguez and Barradas, 2014). The analysis of the effect of each variable on stomatal conductance is determined with simple models that are referred to as boundary-line functions. These simple models consist of the selection of data from a probable upper limit represented by a cloud of points in each graph made by plotting g S as a function of any driving variable (i.e., T A , PAR, VPD, ψ). This analysis has three assumptions: (1) the boundary-line function represents the optimal response of g S to the selected variable (e.g., T A ); (2) the points that fall below the selected function are the result of a change in any of the other variables (e.g., PAR, VPD); and (3) it assumes no synergism at the boundary (Jarvis, 1976;Fanjul and Barradas, 1985;Jones, 1992; Ramos-Vázquez and Barradas, 1998;Barradas et al., 2004; Figure 2).
The relationship of g S as a function of air temperature (T A ) is given by the boundary-line values that conform to a quadratic equation: where a, b and c are the function parameters that allow calculating the optimum temperature (T O ) at which the maximum stomatal conductance (g SMAX ) is reached and the cardinal temperatures (minimum and maximum) when g S is zero. Also, from this relationship (g S vs. T A ), a thermal interval (T I ) can be obtained (considering this interval from the maximum stomatal opening up to 30% closure) and the maximum temperature (T M AX ) before g S decrease by 50%; this value is used because this decrease is considered to represent a potential vulnerability and stress that threatens stomatal performance (Esperon-Rodriguez and Barradas, 2015b). The boundary-line values of stomatal conductance (g S ) as a function of photosynthetically active radiation (PAR) are consistent with a hyperbolic function: 2 | Maximum stomatal conductance (g SMAX , mmol m −2 s −1 ) and averages of stomatal conductance (g S , mmol m −2 s −1 ), leaf water potential (ψ, MPa), air temperature (T A , • C), vapor pressure deficit (VPD, kPa) and photosynthetically active radiation (PAR, µmol m −2 s −1 ) for 15 tree species from Mexico City's urban forest. where A is the asymptotic value of g S or g SMAX , and B is g S sensitivity to changes in PAR.
While the function of gs vs. the vapor pressure deficit (VPD) generates a simple linear equation: where β is the sensitivity of g S to VPD and α is the zero shift. With these plots, the minimum VPD value (VPD MIN ) before g S decreases by 50% was calculated. Similarly, the g S response to leaf water potential (ψ) is also a simple linear equation: where b is g S sensitivity to the ψ, and a is the zero shift. Also, the maximum ψ value (ψ 50% ) before g S decrease by 50% was calculated. Boundary-line points are selected and adjusted to the functions and coefficients mentioned above.
Finally, to predict and analyze g S at any present and future time, a function composition model was used (i.e., multiple regression analysis). In this case, we used a multiplicative model weighted by the effect of each environmental variable (Jarvis, 1976;Ballinas and Barradas, 2016b), i.e., when all these variables act in concert on g S : whereg(T A ),g(PAR),g(VPD) andg(ψ) are the boundary-line functions and are weighted by the effect of each environmental variable taking a proportionality k from 0 to 1; these values change according to the studied species (kT A , kPAR, kVPD, kψ). The maximum value 1 coincides with g SMAX . Theoretically, g SMAX is obtained mainly from the function of g S vs. PAR, being the parameter A of equation (2), and all response variables have the same weight. We used the consistent estimated values of T A (25.0 • C), PAR (300 µmol m −2 s −1 ), VPD (2.0 kPa), and ψ (-1.2 MPa) for model development. These values were derived from data generated from weather stations of Mexico City (Barradas, unpublished dataset) to reduce uncertainty of future climate projections.
Here, we determined species vulnerability by identifying the sensitivity of g S responses to environmental conditions (i.e., T A , PAR, VPD, ψ); that is, a species that is more sensitive to environmental changes will close stomata as a response to small changes in a given environmental variable. For example, in the case of T A , the most sensitive species will be those that have a narrower interval (T I ) confining the optimal function of g S to a smaller interval. For the other predictor variables (i.e., VPD, PAR, ψ), the slope of the response variable indicates the sensitivity of each species. Highly vulnerable species have higher slopes (i.e., high sensitivity) compared to less vulnerable species (Esperon-Rodriguez and Barradas, 2015b). However, we highlight that the stomatal function is species specific, and the functional responses to environmental variables could reflect species plasticity. This response, therefore, can indicate species resilience to climate change (Frank et al., 2017;Ming et al., 2019). Differences among species reflect different evolutionary histories and adaptations to distinct niches, and although all the species studied here share a similar environment (i.e., Mexico City), they respond differently to environmental factors.

Statistical Analysis
To evaluate significant differences, the non-parametrical test Kruskal-Wallis was used to compare the variables of T A , PAR and VPD among sites. Statistical significance was considered at 95% for all cases. We used the software Statgraphics to perform this analysis. For the boundary line analysis, the upper points of the point cloud were chosen, which are assumed to be the maximum function with no restriction of the other environmental or physiological variables. For this analysis, we used the software Table Curve 2D.

RESULTS
The highest g S was recorded for L. lucidum, followed by E. camaldulensis, whereas the lowest g S corresponded to D. viscosa and B. cordata. The lowest ψ (i.e., more negative) belonged to B. cordata Q. rugosa and A. acuminata, and the highest ψ corresponded to P. deltoides and L. styraciflua. All species had the lowest ψ, values between 14:00 and 15:00 h ( Table 2). Similarly, all environmental variables (i.e., T A , PAR and VPD) varied among species ( Table 2). We found significant differences in air temperature when we compared species among sites, finding higher T A in the parks EP and MP (average of 24.8 • C) and lower T A in CU and VP (22.1 • C) (H = 14.85, df = 7, p = 0.05). This result probably reflects the effect of the urban heat island. In contrast, we found no significant differences when comparing VPD (H = 6.350, df = 7, p = 0.15) and PAR (H = 1.356, df = 7, p = 0.153) among parks.
The parameters for each curve and each variable among the 15 species were compared (one case study is shown in Figure 2 and the parameters of the boundary-line analysis for g S vs. T A are presented in Supplementary Table 1). All species presented different T O (ranging from 18.0 to 28.5 • C), optimal thermal interval, and different T I range. Populus deltoides had the lowest T O , whereas E. camaldulensis had the highest T O . Optimal temperature on average and across all species was 22.5 ± 3.2 • C (standard deviation). Regarding T I , P. alba (20.3 • C), D. viscosa (19.9 • C) and B. cordata (19.3 • C) had the widest thermal range, whereas C. occidentalis (11.0 • C) and E. camaldulensis (11.2 • C) had the narrowest intervals across all species (Table 3).
For VPD, P. deltoides and F. uhdei reached the highest rates of g S change, while the lowest rates were registered for B. cordata and D. viscosa (Supplementary Table 2). As for VPD 50% , for all species, the critical temperature and relative humidity values were above 33.0 • C and less than 30% of relative humidity, where D. viscosa and P. deltoides presented the lowest VPD 50% and B. cordata, E. camaldulensis and U. parvifolia had the highest VPD 50% (Table 4). Concerning PAR, the highest and lowest values were recorded in A. acuminata (119.7 ± 155 µmol m −2 s −1 ) and P. alba (410.5 ± 235 µmol m −2 s −1 ) ( Table 2).
Regarding vulnerability, for some species, we found wide thermal ranges (e.g., P. alba and D. viscosa); however, we highlight that although a species interval can be wide, this does not necessarily mean the species is not vulnerable. In the case of the increase in T A , the temperature at which stomatal conductance decreases by 50% is more useful to assess species vulnerability. Here, we identified P. deltoides as the most vulnerable species, since it reaches its T 50% at 27.1 • C ( Table 3). 6 | Predictions of stomatal conductance (g S ) using a multiplicative model (equation 5 from Methods, g SMOD ) and the contribution of each environmental (kT A , kVPD, kPAR) and physiological (kψ) variable that can restrict or enhance stomatal performance in 15 tree species from Mexico City's urban forest.

Speciesg(T A )g(VPD)g(ψ)g(PAR)
kT We considered less vulnerable, those species that present their T 50% at least at 30 • C. The species that presented the highest T 50% were E. camaldulensis and F. uhdei (Table 3).
Leaf water potential was used to analyze tree vulnerability to water availability. Based on the ψ when stomatal conductance decreased by 50% (ψ 50% ), the most vulnerable species were P. deltoides, L. styraciflua and C. occidentalis, as these species had the highest (i.e., less negative) values. In contrast, the least vulnerable species were B. cordata, Q. rugosa, and E. camaldulensis, which also coincide with their lowest ψ values ( Table 5). Concerning VPD and PAR, vulnerability was assessed using the slope of the hyperbola (or sensitivity) given by the constant β and B, respectively, in which the more vulnerable species were those with the highest β and B values. For VPD, the most vulnerable species was P. deltoides (β = 1367.42), whereas B. cordata (β = 136.01) was the least vulnerable. As for PAR, we identified B. cordata (B = 40.21) and L. lucidum (B = 2.08) as the most and least vulnerable species, respectively (Supplementary Table 2).
The performance of the g S response to all variables is presented in Table 6. For U. parvifolia, the proportionality kT A was higher than 1, probably caused by recording high g S at elevated T A . We found P. deltoides and F. uhdei presented negative proportionalities (kψ), indicating high vulnerability. Further, P. deltoides presented the lowest proportionality values in all environmental variables, except PAR, indicating that this is the most vulnerable species to climate change. According to the predictions of g S using the multiplicative model, survival of P. deltoides and F. uhdei might be not possible in urban environments. From our results, it became evident that the effect of ψ was restrictive in all species, except for Q. rugosa, which showed good performance ( Table 6).

DISCUSSION
We used the boundary-line analysis to identify vulnerable species to environmental changes based on their stomatal response. The approximation of the boundary line functions of each variable through the determination coefficient (r 2 ), indicated that these functions were adequate for determining the environmental effect on the g S response. This response showed a higher g S sensitivity to VPD and ψ compared to T A and PAR. Differences in stomatal sensitivity found here may be caused by different microclimates at each study site, which could cause acclimatization of stomata to the urban environment (Barradas et al., 2016). Stomatal responsiveness could also be an adaptation (i.e., protective of water loss) related to different management practices (e.g., irrigation) at each site.
The multiplicative model allowed us to identify the most (i.e., P. deltoides and F. uhdei) and least (i.e., E. camaldulensis, U. parvifolia and A. acuminata) vulnerable species to environmental changes. It is worth noting that P. alba had a better g S performance than P. deltoides, showing that other factors rather than phylogeny, may improve species performance, such as plasticity and climate of origin (Nicotra and Davidson, 2010;Du et al., 2019). However, studies in tropical areas showed that phenomenological models, such as the boundary line analysis, are site-specific (Barradas et al., 2004). Also, for the multiplicative model, we used input data derived from meteorological stations from Mexico City; a different approach might be using modeled meteorological data (e.g., PRISM). Our results, therefore, are site-and species-specific, and caution must be taken if the parameters are extrapolated to other sites and individuals under different conditions. As a response to climate change, global warming will continue toward the end of the century increasing air temperature in urban areas (Esperon-Rodriguez et al., 2021;IPCC, 2021). The average air temperature in Mexico City is expected to be 24.5 • C by 2050 (Ballinas, 2011) and daytime temperatures will likely exceed the species' temperature ranges found here. While the annual average air temperature in Mexico City keeps increasing, extreme temperatures will be aggravated by the UHI effect negatively affecting some species by limiting their stomatal function. In this study, we found that temperatures above 27 • C limited the stomatal opening, although some species could maintain stomata open (T 50% ) at higher temperatures. This was the case of E. camaldulensis and F. uhdei, which had g S associated with T 50% of up to 39 • C. In contrast, species such as P. deltoides, had the lowest T 50% at 27.1 • C. Additionally, changes in precipitation and cloudiness, which affect solar radiation, may also affect the stomatal function with effects on species growth and performance.
A direct effect of the increase in air temperature is the increase in VPD when water is limited and relative humidity is low. The response of g S to high temperatures can be uncertain, particularly when these temperatures are coupled with low relative humidity and high VPD (Medlyn et al., 2002;Esperon-Rodriguez and Barradas, 2014). The increase in VPD restricts g S , thus, species capable of maintaining stomata open when VPD is high are less vulnerable (i.e., Q. rugosa, E. camaldulensis and B. cordata). In contrast, species with low g S under high VPD conditions are more vulnerable (i.e., E. americana, A. negundo, L. styraciflua, and P. deltoides). From our analysis, we found all species were vulnerable to changes in VPD and ψ indicating sensitivity to water stress conditions. However, caution must be exercised when interpreting the proportionalities in the model because if two proportionalities are negative, for example, the value of g SMOD is positive. As for PAR, its effect was less significant compared to the other variables (i.e., VPD, ψ) because C 3 plants become saturated toward 300 µmol m −2 s −1 (Jones, 1992), hence, g S would be limited by PAR below the saturation point; that, however, did not occur in our study.
Finding a high vulnerability to VPD and ψ, two variables directly related to water availability, highlights the importance of water in urban environments to keep stomatal function. Water is crucial to maintain g S for leaf cooling via evapotranspiration, and therefore is associated with the cooling benefits of trees even under extreme high temperatures (Drake et al., 2018). Leaf temperature is an important parameter; a decrease in g S can be associated with increases in leaf temperatures (Lu et al., 1994) and leaf temperatures can exceed air temperature by several degrees (Miller et al., 2021) leading to leaf damage (Hueve et al., 2011). Thus, identifying resistant species to high temperatures in water-limited environments is beneficial in urban environments because canopy loss can have negative environmental and social impacts (e.g., decrease in cooling benefits) (Sanusi and Livesley, 2020).
Here, we used leaf water potential as an indicator of plant resistance to drought (Fanjul and Barradas, 1987;Jones, 1992). Furthermore, decreases in ψ can reduce g S Barradas, 1985, 1987;Esperon-Rodriguez and Barradas, 2015b). When water is limited, species that respond by decreasing their ψ are less vulnerable. This was the case of B. cordata, Q. rugosa, A. acuminata, and E. camaldulensis, which presented ψ lower than -2.0 MPa. In contrast, the most vulnerable species were P. deltoides and L. styraciflua. For L. styraciflua, previous research showed its vulnerability to water-limited natural environments (Esperon-Rodriguez and Barradas, 2015b), and here, we confirm its vulnerability also in urban environments. To confirm the species vulnerability to ψ, we recommend using the leaf turgor loss point as an additional trait. The less vulnerable species would be those that are able to maintain ψ above their turgor loss point under water stress conditions, allowing gas exchange and evaporative cooling to continue.
Our results can aid species selection in a changing climate to minimize economic losses related to urban tree failures. By working together, scientists, policymakers and urban planners can develop and implement strategies and management actions that not only guarantee the preservation of urban forests but also benefit local communities.

DATA AVAILABILITY STATEMENT
The datasets presented in this article are not readily available because when data is generated for an investigation, by law it is reserved for 10 years in Mexico. Requests to access the datasets should be directed to corresponding author.

AUTHOR CONTRIBUTIONS
VB and ME-R designed the research, collected, analyzed the data, and wrote the manuscript. Both authors contributed to the article and approved the submitted version.

FUNDING
This wok was supported by the PAPIIT Programa de Apoyo a Proyectos de Investigación e Innovación Tecnológica, DGAPA, UNAM, Grant/Award No. IT200620.