Study Sites

The three study sites are located in the central Mexico, in natural T. mucronatum Ten. populations, located in Estado de México (MEX), Querétaro (QRO), and Morelos (MOR); between 18° 43’ 8.4” and 20° 33’ 00” N, and 98° 44’ 06” and 100° 05’ 38.4” W (Figure 1). The moisture regime is ustic in MEX and QRO and udic in MOR. The soil temperature regime of MEX is thermal and the regime of QRO and MOR is isothermal (Table 1 and Figure 1).

Figure 1 Study sites and climographs. Estado de México (MEX), Querétaro (QRO) and Morelos (MOR). The climographs were elaborated based on the information included in the ERIC III database Version 3.2 (Rapid Extractor of Climatic Information) (^{IMTA, 2013}), taking into consideration 3-5 nearby meteorological stations. 

Sampling and dendrochronological dating

Nine trees per site were sampled, using circular cross-section trunks as possible. Three increment cores were taken from several points around the trunk, at a height of 1.3 m above the ground level, with a 12-mm internal diameter Pressler increment core (HAGLÖF, model BS009). The wood samples were dried at room temperature and sanded with different grades of sandpaper, in order to improve the contrast of the growth rings. The total ring width was measured with a 0.001-mm precision VELMEX system (^{Robinson and Evans, 1980}). The dating was done with standard dendrochronological techniques (^{Stokes and Smiley, 1968}). The quality of the dating was verified with the COFECHA program (^{Holmes, 1983}).

Estimation of the growth in basal area

The basal area increment (IAB) was estimated from the measurements of the ring width, assuming regular concentric rings. The following formula was used:

where R is the radius and n, the year in which the ring was formed. IAB is a better representation of the tree growth, because (unlike ring width) it remains stable for a longer period, during the development of the tree and it enables a direct interpretation of the growth rate, expressed as area in cm² year^-1 (^{Poage and Tappeiner, 2002}). The IAB trend through time was estimated using polynomials of order five and the exponential model.

Identification of extreme humidity events (dry and wet periods)

Geometric and biological trends not related to the climate (competition, suppression, and release) were removed with the ARSTAN program (^{Cook and Holmes, 1984}), which generated a series of standardized annual growth in each population. To highlight low frequency extreme events (dry and wet periods), a spline was inserted in each series of ring width indexes (RWI) (^{Fritts, 1976}; ^{Cook and Peters, 1981}). Afterwards, dry and wet periods detected in growth series of T. mucronatum Ten. were compared with other dendroclimatic series developed in Mexico, in order to corroborate common responses to extreme climatic events at the regional level (^{Fritts, 1991}; ^{Stahle et al., 2009}; ^{Cook et al., 2010}).

Analysis of isotopes in wood

The three most representative trees of each site (according to their chronological significance( were selected out of the nine trees sampled in the field. A sawdust sample was obtained from each tree and from each identified extreme period (dry and wet). Samples were taken during six-year periods, in order to collect enough sample material and perform the isotope analysis. Four dry periods were analyzed: I: 1895-1900, III: 1947-1952, VI: 1996-2001, and VII: 2007-2012; likewise, three wet periods were analyzed: II: 1910-1915, IV: 19651970, and V: 1988-1993. The sawdust samples of each period were weighed on a microbalance (Sartorius, 0.01 mg) and placed in tin and silver capsules (Elemental Microanalysis Limited, Okehampton, Devon, UK) to analyze δ¹³C and δ¹⁸O, respectively. Determination of Δ¹³C was performed using an Europa ANCA GSL elemental analyzer connected to a PDZ Europa 20-20 isotope ratio mass spectrometer (Sercon, Cheshire, UK). Determination of δ¹⁸O was performed in a HEKAtech HT Oxygen Analyser (Wegberg, Germany) with interface to the abovementioned spectrometer. All the analyzes were performed at the UC Davis Stable Isotope Facility and the values were expressed in terms of their respective standards -δ¹³C: VPDB (Vienna Pee Dee Belemnite); δ¹⁸O: VSMOW (Vienna Standard Mean Ocean Water).

Estimation of physiological variables (δ¹³C, C_i, and iWUE)

The evolution of the carbon isotope composition (¹²C and ¹³C) in the atmosphere was taken into account in order to estimate the physiological variables. The abundance of ¹³C in the plant is expressed in relative source-product terms (Δ¹³C, discrimination of carbon), taking into account the atmospheric variability (2) (^{Farquhar et al., 1982}).

where Δ is the discrimination against ¹³C, δ¹³C_air is the carbon isotope composition in the air (source), and δ¹³C_plant is the carbon isotope composition of the product (plant biomass). Discrimination of carbon is used in environmental physiology for long periods with changes in the composition of ¹³C in the atmosphere (pre-industrial vs. present times) (^{Mc Carroll and Loader, 2004}).

Intercellular CO₂ concentration (C_i) was estimated with Equation 3.

where a is the discrimination against ¹³CO₂ while the diffusion through the stoma takes place (≈4.4 ‰), b is the net discrimination resulting from carboxylation (≈27 ‰), C_i and C_a are the CO₂ internal and environmental concentrations, respectively. Taking into consideration Fick’s First Law (A = gCO₂ (Ci /Ca)), this equation can be expressed as follows (4).

where A is the assimilation rate of CO₂, and gCO₂ is the stomatal conductance to CO₂. Given that gH₂O (or leaf conductance to water vapor) is 1.6 gCO₂, the values of Δ¹³C can be connected with the A/gH₂O quotient; which represents the iWUE of the plant (^{Farquhar and Richards, 1984}; ^{Mc Carroll and Loader, 2004}).

Statistical analysis

Data were analyzed using a variable-spacing model of repeated samples in time, using the MIXED procedure of SAS 9.3 (^{SAS, 2011}). Intercepts for each tree and the site were considered to be random; the effect of time, moisture regime, and interaction were tested to (p≤0.05) (^{Fitzmaurice et al., 2004}). To test the differences between sites, a multivariate analysis of variance (MANOVA) was performed, complemented by contrast tests between sites.

Interpretation of trends in physiological variables

Changes in physiological variables were interpreted considering the scenarios proposed by ^{Saurer et al. (2004)}, and the relations between δ¹³C and δ¹⁸O regarding physiological variables (^{Scheidegger et al., 2000}). The three scenarios suggested by ^{Saurer et al. (2004)} consider that: 1) C_i is constant (increased iWUE); 2) C_i / C_a is constant (increased iWUE); and, 3) C_i - C_a is constant (unchanged iWUE). Those scenarios are interpreted as follows: C_i is independent of the changes in C_a ; C_i proportionally adjusts to the changes in C_a ; and the changes in C_i and C_a occur at the same rate. ^{Linares and Camarero (2012)} have provided a detailed explanation of these scenarios. The method proposed by Scheidegger et al. (2000) uses the relations between δ¹³C and v¹⁸O to infer changes in stomatal conductance (gl), photosynthetic capacity (Amax), and iWUE.

Chronological series and importance of low frequency periods

The chronological length varied from 109 to 233 years. Thus, the trees studied were not considered senile given that T. mucronatum Ten.’s longevity can reach several centuries or go beyond one-thousand years (^{Villanueva-Díaz et al., 2010}). The intercorrelation of the series was significant for all sites, with values ranging from 0.36 to 0.51; meanwhile mean sensitivity varied from 0.34 to 0.43 (Table 2). Both indicators show the influence and variation of a common climatic signal recorded in the ring width; therefore, the growth series are suitable to detect high and low frequency events (^{Fritts, 1976}; ^{Constant et al., 2010}). For example, a 0.32 degree of intercorrelation is considered suitable for general trend purposes (^{Holmes, 1983}). The degree of correlation between series tends to increase in dry places (^{Stahle et al., 2009}), but the higher inter-correlation was found in the intermediate humidity site (QRO). The mean sensitivity is a relative index of how much ring width changes from one year to the next: in this case, the 0.34-0.43 values are expected, because mean sensitivity varies from 0.15 to 0.65, depending on the type of weather. The autocorrelation parameter in chronologies indicates the relationship between the measurements of subsequent years. In this study, the 0.64-0.78 values are considered high and are the result of the annual humidity and drought cycles in study sites (Table 1). The greatest ring-width variation in the MEX site may be related with its lowest precipitation and the highest reduction in its phreatic zone (^{Correa-Díaz et al., 2014}).

The chronology of this study was associated with 83 % of extreme dry events and with 44 % of extreme wet events detected in other chronologies of montezuma baldcypress in the states of Puebla, Guanajuato and Querétaro (^{Therrell et al., 2006}; ^{Cortés-Barrera et al., 2010}; ^{Stahle et al., 2011}), as well as for chronologies of Douglas-fir (Pseudotsuga menziesii Mirb. Franco) in Nuevo León (^{Arreola-Ortiz and Návar-Cháidez, 2010}). Since rings are narrower in dry years, a greater correlation exists with extreme dry events. Although wet events generate greater ring widths, the degree of influence is limited by the water storage capacity of the soil profile, and surpluses (resulting from overflow or evapotranspiration( do not influence growth, particularly in periods of physiological inactivity (^{Sarris et al., 2013}).

Changes in the composition of isotopes (δ¹³C and δ¹⁸O)

The mean values of δ¹³C were -24.05±0.12 ‰, -24.59±0.18 ‰, and -25.38±0.11 ‰, in MEX, QRO and MOR, respectively. The ¹³C signature in the growth rings is mainly controlled by the difference between the atmospheric CO₂ concentration and the tree; i.e., it depends on the photosynthetic rate (carboxylation) and stomatal conductance (^{Mc Carroll and Loader, 2004}; ^{Holzkämper et al., 2008}). The changes of δ¹³C through time or according to moisture regime were not significant (p≤0.05) in MEX and QRO. In MOR, the changes of δ¹³C through time depended on the moisture regime (dry or wet event) (Figure 2, Table 3). In this sense, the composition of ¹³C in the wood did not reflect the δ¹³C changes in the atmosphere, like in other forest species (^{García-G et al., 2004}; ^{Gómez-Guerrero et al., 2013}; ^{Silva et al., 2015}). Generally, the δ¹³C signature reflects the δ¹³C dilution caused by the use of fossil fuels (^{Mc Carroll and Loader, 2004}). Although the differences were not significant through time, the (24.05 ‰, (24.59 ‰, and (25.38 ‰ mean values do show significant differences (p≤0.0001), and reflect an expected humidity gradient from lowest to highest (MEX-QRO-MOR). The δ¹³C signature in wood does not show the δ¹³C changes in the atmosphere, possibly because the montezuma baldcypress trees grow near water reservoirs and the discrimination of ¹³C is stable throughout the year. However, the net change expressed as Δ¹³C -observed in all three sites- indicates that changes in C_a do have an influence in gas exchange (Figure 2).

Figure 2 Composition of carbon and oxygen isotopes (δ¹³C and δ¹⁸O) in growth rings, discrimination (Δ¹³C), intrinsic wateruse efficiency (iWUE), internal partial pressure of CO2 (Ci), and atmospheric CO2 concentration (Ca), in three natural populations of montezuma baldcypress trees (Taxodium mucronatum Ten.). The bars indicate standard errors. 

The analysis of contrasts indicated significant differences (p≤0.05) of δ¹³C between the MEX and MOR sites, during periods III, IV, and VII. There were no differences between the MEX and QRO periods. In QRO and MOR, only period IV was statistically different (Table 4). Sites with contrasting rainfall (MEX and MOR) are the only ones that show δ¹³C differences between them, and, if a marginal significance (p≤0.079) is taken into consideration, the differences are evident from periods III to VII (Table 4). The more negative values of δ¹³C in MOR indicate less discrimination and less hydric stress as a result of the site’s higher humidity (Figure 2).

The values of δ¹³C in our study are lower than those reported for Pinus tabulaeformis and P. taiwanensis in subtropical forests in China (^{Xu et al., 2014}); however, they are similar to those reported by ^{Silva et al. (2009)} for Araucaria angustifolia in Brazil, and by ^{Gómez-Guerrero et al. (2013)} and ^{Silva et al. (2015)} for Abies religiosa and Pinus hartwegii in Mexico.

The values of δ¹⁸O were different (p≤0.05) through time in all three sites, with a 4.0 ‰ mean increase over the last century. The time x regime interaction was not significant anywhere. The contrasts indicated differences between MEX and QRO in periods V and VII. For MOR and QRO, there was a statistical difference of δ¹⁸O in period V (Table 3 and 4).

The increase of δ¹⁸O through time in all three sites indicates possible changes in the water source or greater evapotranspiration (^{Kohn and Welker, 2005}). Since there is not data of δ¹⁸O in rainwater in study sites, it is difficult to separate the water source effect or demonstrate that evapotranspiration (and temperature) increased in the twentieth century. However, when nearby weather stations are taken into consideration, there is indeed evidence of an increase in average annual temperature over the last 50 years (1.2-1.4 °C in all sites). Therefore, the 4.0 ‰ increase in the three sites under study is likely a reflection of the temperature and evapotranspiration increase that enriched the ¹⁸O signature in the water used by trees.

The explanation of an evapotranspiration and temperature increase is congruent with the water supply crisis in the study sites. According to the National Commission for Water (CNA), there is an 111 755 911 m³ year^-1 reduction in MEX site aquifers and 12 929 028 m³ year^-1 in the QRO site (^{CNA, 2015}). The groundwater decrease forces the montezuma baldcypress trees to extend their root system to deeper layers and it limits their aerial growth (^{Enríquez-Peña and Suzán-Azpiri, 2011}).

Physiological variables (Δ¹³C), C_i and iWUE)

The Δ¹³C, C_i and iWUE differences over time were significant (p≤0.05) in the three sites. The moisture regime and the time x regime interaction were significant for C_i in MEX and MOR (Table 3). The decrease in Δ¹³C and increase in iWUE were more evident after 1950 in all three sites (Figure 2). The increase in C_i was noticeable after 1970 (Figure 2). There were statistical differences in Δ¹³C, C_i, and iWUE between MEX and MOR in periods III, IV and VII, but there were no significant differences in Δ¹³C, C_i, and iWUE between MEX and QRO. For MOR and QRO, only period IV was significant (Table 4).

During the twentieth century (p≤0.05), iWUE increased from 77 to 116 μmol mol^-1 (a 50 % increase). However, it is well-known that the annual rate increased from 0.35 to 0.65 μmol mol^-1 after the extreme event of 1950 (Figure 2). In recent decades, iWUE has increased in different forest ecosystems (^{Saurer et al., 2004}; ^{Peñuelas et al., 2011}; ^{Battipaglia et al., 2013}; ^{Xu et al., 2014}), but it has had a variable relationship with tree growth, including in our study sites. Theoretically, the increase of iWUE improves the drought resistance of vegetation, because it reduces the cost associated with the transport of water within the tree (^{Raven, 2002}), but tree growth does not increase in all sites.

Regardless of their locality, montezuma baldcypress trees respond to climate changes and to the increase of atmospheric CO₂ concentration in the twentieth century. The most sensitive variable was C_i , which varied according to the weather and moisture regime in the sites where weather is most contrasting (MEX and MOR). The periods with the largest significant differences between sites were IV and V, which match wet periods. This result indicates that -compared with rainfall- temperature exerted greater influence on the physiological changes in montezuma baldcypress trees. This result is similar to that of ^{Linares and Camarero (2012)}, who observed that the physiological responses of Abies alba had a deeper relationship with changes in temperature. Our study also shows that temperature is the variable that has a deeper relationship with changes in iWUE. The recent increase in average annual temperature matches an increase in iWUE in all sites (Figure 5).

The net reduction of Δ¹³C in all sites proves the physiological response of trees to the increase of C_a (^{Silva and Horwath, 2013}). Regarding extreme events, C_i has increased at a higher rate in dry events than in wet ones (Figure 3), but the change is greater in humid places with isothermal regime (such as MOR). The increase of C_i indicates that there is more carbon available within the leaves, although their assimilation does not always increase, since the growth rates in the basal area of the trees increased or decreased (Figure 6). The results indicate that, in recent decades, the temperature increase has had a differentiated impact in tree growth.

Figure 3 Moisture regime effect (dry and wet periods) in the composition of isotopes of carbon 13 (δ¹³C) and internal partial pressure of CO₂ (C_i). (A) MEX, (B) QRO y (C) MOR. The dotted line represents linear effect in variables. Bars indicate standard errors. 

Interpretation of physiological variables (Δ¹³C, C_i, and iWUE)

The trends of the physiological variables are connected with the scenarios posed by ^{Saurer et al. (2004)}, where C_i and C_i /C_a remain constant (Table 5). However, in the first case, an increase of Δ¹³C was expected, but this occurred only in wet periods (Figure 3). In the second case, although changes to Δ¹³C were not expected, they decreased (Table 5). The third scenario (C_i -C_a) shows less correspondence with physiological variables. There was no correlation between δ¹³C and δ¹⁸O (Figure 4) which (according to conceptual model of ^{Scheidegger et al. (2000)} (suggests that variation in iWUE are the result of the proportional adjustments of photosynthesis and stomatal conductance. For data representation purposes, we have included the δ¹³C and δ¹⁸O deviations in relation to the global mean, as suggested by ^{Barnard et al. (2012)}.

Table 5 Expected (black) and observed (gray) trend of the physiological variables. 

Escenario	Ci	Ci /Ca	Ci -Ca	δ13C	Δ13C	iWUE
Ci constante	—
Ci /Ca constante		—
Ci -Ca constante			—

Scenarios match those suggested by ^{Saurer et al. (2004)}. The trends found in this study are indicated in gray. Result that partially match or differ from the conceptual model are indicated with “Partial” and “NO”, respectively.

Figure 4 Deviation of δ¹³C and δ¹⁸O in growth rings in relation to the global mean. Data are for the three sites: MEX (triangle), QRO (square), and MOR (circle). 

Figure 5 Ratio of atmospheric CO₂ concentration (C_a), internal partial pressure of CO₂ (C_i), average annual temperature (average of extreme wet or dry years from 1947-1952), and intrinsic water-use efficiency (iWUE), in three natural areas of Taxodium mucronatum Ten. 

The response of montezuma baldcypress to recent C_a changes match theoretical results of two of the scenarios proposed by ^{Saurer et al. (2004)}, in which C_i and C_i /C_a are constant. However, the expected increase of δ¹³C in the first scenario was partial because it only happened in wet periods. The inconsistency with regard to the second scenario (C_i /C_a constant) is that Δ¹³C was expected to remain relatively constant; however, its decrease after 1950 was highly significant in all sites (Figure 2 and Table 3). The results indicate that in the case of montezuma baldcypress trees, C_i tends to remain constant in face off C_a increases, showing the active role of the stomatal opening. Specifically, this appears to occur during extreme wet events and in low and high humidity places (like MEX and MOR). The results indicate that, in extreme drought events, trees lean towards the second scenario, with C_i and C_a proportional adjustments. According to the conceptual model proposed by ^{Scheidegger et al. (2000)}, a positive correlation between δ¹³C and δ¹⁸O is not observed (Figure 4); therefore, the idea that the increase of iWUE was caused by the reduction in stomatal conductance (g_s) and stable photosynthesis (A) levels was dismissed. The horizontal trend of δ¹³C and δ¹⁸O deviations with regard to the general mean suggest proportional adjustments of A and support the second scenario proposed by ^{Saurer et al. (2004)}.

Relation between physiological variables and climate

The C_i and iWUE variables show a linear increase in relation to C_a, in the twentieth century (Figure 5). This could suggest that the increase of atmospheric CO₂ is reflected in higher gas concentration in the trees. Although C_i values are different in each site, the same slope or increase rate was found in both variables. The average annual temperature also plays an important role in iWUE and C_i in MEX and QRO; a positive relationship between temperature and C_i and iWUE increase was found; compared with temperature, rainfall did not explain C_i and iWUE changes. The results of this study are similar to those found by ^{Silva et al. (2009)} in Araucaria angustifolia and by ^{Tognetti et al. (2014)} in Fagus sylvatica and Nothofagus spp. This relation also depends on the steam pressure deficit in the area where trees grow; however, this may also reflect a regulation process (homeostasis) of forest species, in face of climate variation, because the correlation between iWUE and C_i tends to increase in recent decades (^{Tognetti et al., 2014}).

Growth in basal area

The growth in basal area for all three sites, using an equation with negative exponents for MEX and quintic polymorphic equations for QRO and MOR, show three trends: reduction in growth, stable growth, and increased growth, in MEX, QRO, and MOR (Figure 6). However, since the last decade, growth has decreased everywhere. In 1920 MEX’s IAB exceeded that of the other sites; afterwards, its growth rate showed a progressive decrease from 1920 to 2012, until it reached values four to five times lower than in the other two sites, even if the younger trees were found in this site. A probable explanation for this situation is that the woodland has been significantly affected by changes in water availability and the effects of land-use change (^{Correa- Díaz et al., 2014}).

Figure 6 Growth rate in basal area (IAB) in three natural areas of Taxodium mucronatum Ten. MEX trends were determined with a model with negative exponents, and MOR and QRO, with a quintic polynomial. The bars indicate standard deviation.