Soil Respiration Is Inﬂuenced by Seasonality, Forest Succession and Contrasting Biophysical Controls in a Tropical Dry Forest in Northwestern Mexico

: Soil respiration (R S ) is an important component of the C cycle because it contributes signiﬁcant CO 2 emissions to the atmosphere that result from metabolism and respiration of its autotrophic and heterotrophic components. However, the relative importance of different biophysical controls that drive the variability of this ﬂux and their inﬂuence along forest succession pathways is still unknown. We incorporate multiyear R S , ecosystem ﬂux and meteorological measurements in old-growth (OG), mid-secondary (MS) and early-secondary (ES) tropical dry forests (TDFs) with the goal of assessing the temporal variation of R S and identifying the biophysical controls at each site by applying structural equation models (SEM). Along forest succession, R S followed the pattern of precipitation events; we identiﬁed by the end of the wet season that R S was sustained by a longer period at OG, while in MS and ES, R S decreased according to the soil moisture availability. According to SEM, soil moisture and soil temperature exert an effect on the variability of R S in all sites. However, we found that R S was also controlled by the vapor pressure deﬁcit at MS and gross primary production at OG and ES. Our results suggest that seasonality has a different impact on R S along forest succession in TDFs found in northwestern Mexico and highlights the relevance of considering additional biophysical controls of R S for a better understanding this critical process of the C cycle


Introduction
Soil respiration comprises the CO 2 efflux resulting from metabolism and respiration of plants (i.e., roots and mycorrhizae), soil fauna and microorganisms (i.e., decomposers) that determine biogeochemical processes within the soil [1,2] and it is accompanied by nonbiological CO 2 sources, namely carbonate mineral weathering [3], CO 2 dissolution and gas diffusion [4]. This flux has received wide attention for its tight relationship with ecosystem productivity, soil fertility and the carbon (C) balance in terrestrial ecosystems [5,6]. R S is the second-largest flux in magnitude after photosynthesis as a global C cycle component [7,8] and represents the major component of ecosystem respiration (R eco ) [9] ranging from 30% to 80% of the total R eco in forests and seasonally dry ecosystems at annual scales [10,11].
The understanding of the intra-annual and interannual variability of R S and the mechanisms underlying this important process are still poorly explored in water-limited ecosystems [12,13], because seasonality plays a critical role in determining the temporal patterns and the interactions between biotic and abiotic drivers [14,15]. Furthermore, spatial heterogeneity of soils and vegetation in drylands is a defining feature that also influences ecosystem C fluxes such as R S [16,17].
Generally, at different temporal scales, R S rates are affected by environmental factors such as soil temperature and soil moisture, which are usually considered as the two main controlling factors [18]. However, soil physicochemical properties, substrate supply, soil enzymatic activity and soil microbial community also participate in this process [19][20][21]. Furthermore, the crucial role of plant phenology, plant functional types and gross primary production (GPP) in the temporal variability of R S has been recognized [22][23][24][25][26]. For instance, partitioning R S into its autotrophic and heterotrophic components to inquire about their contribution is still limited due to methodological challenges, leading to significant uncertainties in estimating the global C budgets [27,28], although some studies have been able to estimate the autotrophic contribution of R S , which can represent a significant fraction of the total annual R S in ecosystems [29]. Despite that in seasonally dry ecosystems soil moisture (SWC) exerts an overriding influence on ecosystem processes [30,31], there is still debate about the direct and indirect effects of SWC and the relative importance of different environmental (atmospheric moisture, air temperature, solar radiation, evapotranspiration) and biological drivers (GPP) on R S [32][33][34].
Globally, tropical dry forests (TDF) are considered an important hot spot of biodiversity, endemism and a fundamental source of goods and ecosystem services, implying high productivity and C storage rates [35,36]. The most notable feature of this forest is the seasonality caused by the variability of precipitation that leads to well-defined dry periods of low biological activity [37] and active periods with very dynamic eco-physiological activity from the vegetation [38] and soil organisms [39]. Therefore, TDFs support a variety of water-limited tolerant deciduous trees and shrubs with specific adaptive traits such as leaf drop and regrowth, along with thorny and succulent species and a few evergreen species holding their leaves and maintaining ecophysiologically active during dry periods [40][41][42].
Similarly, soils from these systems are influenced by seasonal drying-rewetting transitional events that strongly regulate the dynamic of soil biogeochemical processes [43]. For example, high rates of litterfall occur during the dry season, resulting in a forest floor with slow decomposition rates and a low CO 2 release [44,45]. At the onset of the wet season, litterfall promotes the rapid decomposition and mineralization of forest floor from the previous growing seasons [46], and in conjunction with microbial activity and nutrient transformations triggers a "Birch effect", where a large pulse of CO 2 is released upon soil rewetting [47,48].
TDFs are also one of the most threatened ecosystems due to land-use and land-cover changes due to anthropogenic activities, such as agriculture and grazing. These landscape transformations result in an important loss of vegetation cover [49,50], creating a mosaic of different stages of forest succession spanning from abandoned lands to secondary and old-growth forests [51]. Forest succession from TDFs in northwestern Mexico is triggered by land abandonments, followed by forest recovery, where seedling recruitment, resprouting and dominance of shallow-rooted fast-growing pioneer species occur, and then transition to the recruitment of primary forest species until the representative vegetation from TDF is settled and matures, including well-defined understories [52]. This forest succession path in conjunction with the precipitation variability and the effect of extreme climatic events has led to contrasting effects on ecosystem structure, species composition and soil properties [53][54][55]. Nonetheless, the effects of successional changes in carbon, water and energy fluxes have been less attended [56].
Despite the key contribution of R S to ecosystem C balance, there is still a lack of mechanistic information on the seasonal and interannual variability of R S from TDF soils, and on the relative importance of different biophysical controls and the variation of their influence across forest succession. In this study, we incorporate the analysis of monthly R S measurements and accompanying ecosystem flux and meteorological observations along a forest succession in a tropical dry forest in northwestern Mexico between 2015 and 2019, with the main goals of determining the temporal variations of R S and identifying the key biophysical controls of this flux in a TDF in northwestern Mexico. Based on this framework, we postulate the following three main hypotheses: (1) the magnitude of R S along successional sites will have significant differences at annual scales, and the variation would be strongly influenced by seasonality; (2) the relative importance of biophysical controls on R S such as solar radiation, air temperature, vapor pressure deficit, soil moisture and soil temperature will be different across TDF succession sites; and (3) since GPP represents the autotrophic activity of the system, it should be an important proxy for root respiration. Therefore, a different effect on R S will be expected across succession as the variations of GPP may be influenced by the metabolic activity of fast-growing pioneer species with shallow roots or mature vegetation with well-defined understories.

Site Description
The study site is located in the Area de Protección de Flora y Fauna Sierra de Álamos-Río Cuchujaqui (SARC) in northwestern Mexico. The SARC is registered to the Man and Biosphere Program Biosphere Reserve from UNESCO and to Comisión Nacional de Áreas Protegidas [57]. This study focuses on the Monte Mojino Natural Reserve, a private protected area owned and managed by Nature and Culture International which lies within SARC. The study sites belong to the MexFlux network, which is a system of ecohydrological monitoring sites equipped with Eddy Covariance technique at terrestrial, marine and coastal ecosystems in Mexico [58].
The climate in the area is classified as warm and semiarid with summer rains [59], a mean annual air temperature of 24 • C, and mean annual precipitation of 729 mm yr −1 ; more than 80% of the annual precipitation is recorded during the wet season due to the influence of the North American Monsoon System, with just a few rainfall events during winter [60,61]. For this study, we use a seasonal scale framework of two periods: (i) a dry season from November to May and (ii) a wet season from June to October, which includes the months of maximum growth rate and biological activity.

Forest Succession Sites
Our study sites are located in the northernmost limit of the TDF in America [46]. Three sites were strategically selected according to the conceptual TDF chronosequence from Álvarez-Yépiz et al. [57]. Following the criteria from [57] and advice from local inhabitants and forest ranges, the successional stages of the three sites were selected by the relative importance of the pioneering species vs. species representative of mature stages and the time since abandonment after clearcut. Currently, all sites are free from large-scale management practices and free of ranching activities.
The first site is an early-secondary forest (ES; 26.99 • N, −108.78 • W) that has a notable dominance of Acacia cochliacantha, considered a fast-growing legume and pioneer species. This site was previously used for local agriculture after being abandoned from management practice activities and now has a recovery period of~10 years. The soil organic matter content is 2.8 ± 0.001% and the litterfall production is 3.66 Mg ha −1 year −1 [62,63]. The canopy openness (%) of this site is 0.27 ± 0.11 with a mean leaf area index (LAI) of 1.51 ± 1.69. On the other hand, we selected a mid-secondary forest (MS; 27.00 • N, −108.77 • W) that has been in recovery for over~40 years after the last clearing. The main feature of this forest is the recruitment of primary forest species and a lesser dominance of Acacia cochliacantha [52]. For this site, litterfall production is estimated at 2.70 Mg ha −1 year −1 [62] and a soil organic matter content of 4.2 ± 0.1% [63]. This site presents a canopy openness of 0.32 ± 0.11% and a mean LAI of 1.44 ± 1.53.
Finally, the old-growth forest (OG; 26.99 • N, −108.78 • W) is characterized by a high biomass content and shows a greater floristic composition with a multispecies assemblage of representative vegetation from TDF such as Lysiloma divaricatum, Lysiloma watsonii, Croton lindquistii, Ceiba acuminata, Tabebuia impetiginosa, and Bursera spp. and a well-developed understory [52,57]. This forest has never been cleared, according to local residents and landowners. Litterfall production in this site is estimated at 4.32 Mg ha −1 year −1 [64] with a soil organic matter content of 5.5 ± 0.2% [63], a canopy openness of 0.25 ± 0.11% and a mean LAI of 1.55 ± 1.71. According to the total relative importance values, the most representative species at each site were: L. divaricatum and Croton lindquistii account for 30% at OG, while L. divaricatum with C. flavescens were 27% at MS, and A. cochliacantha represented 20% at ES [65]. . Within each transect, sixteen PVC soil collars (10 cm diameter and 8 cm in height) were installed with a 20 m distance among them. Due to site access conditions to sampling points over different years, between eight and twelve soil collars were measured at each site during each sampling campaign. All measurements were performed within a time span of three hours between 9 a.m. and 12 p.m. to avoid drastic temperature changes and potentially being representative of daily means [66,67]. R S measurements were carried out with a soil portable infrared gas analyzer (LI-8100, LI-COR Biosciences, Lincoln, Nebraska, USA) coupled to a 10 cm survey chamber (model: 8100-102, LI-COR Biosciences, Lincoln, NE, USA). Measurement cycles started by measuring the depth of each soil collar from the soil surface to the top, in order to adjust the total volume of the system that is required for flux calculations. Changes in soil CO 2 concentration were measured and recorded for 180 s; in each soil collar a 20 s prepurge prior to each measurement with the soil chamber open and a 20 s postpurge after each measurement was completed in order to remove possible remaining gas inside the soil chamber and the system [68,69]. Finally, R S rates were derived by fitting a linear equation in SoilFluxPro ® (version 4.0.1, LI-COR Biosciences, Lincoln, NE, USA); the first 30 s recorded of the dataset were discarded to allow stabilization conditions inside the chamber in the field.

Soil Respiration Measurements
In addition, soil temperature was measured with a 15 cm-long thermocouple thermometer inserted in a 45-degree angle (Barnant Co., Barrington, IL, USA.) and 12 cm-long soil moisture probe (Theta Probe ML2x, Delta Services, Cambridge, UK) was used for soil volumetric water content; these measurements were conducted simultaneously in an adjacent area to the soil collars. Carbon (%C) and nitrogen (%N) concentrations in soil were determined from eight composite samples per site. Samples were weighted (3.5 mg) and loaded in tin pressed capsules (5 × 9 mm) and the C and N contents were obtained by flash combustion in an elemental analyzer (Flash 2000; Thermo Fisher Scientific, Walthman, MA, USA). Soil texture and pH were estimated with the Bouyoucos method and a pH meter (1:2 H 2 O), respectively. Soil physicochemical variables for each site are presented in Table 1.

Ecosystem Flux and Meteorological Measurements
Net ecosystem CO 2 exchange (NEE) and water (ET) flux measurements in the study sites were estimated with the eddy covariance technique (EC). In brief, the EC system in the ES site was coupled with a three-dimensional sonic anemometer (WindmasterPro, Gill Instruments, Lymington, UK) and an infrared gas analyzer (LI-7500A, LI-COR Biosciences, Lincoln, NE, USA). Data were collected at 10 Hz through an analyzer interface unit (LI-7550, LI-COR Biosciences, Lincoln, NE, USA). On the other hand, meteorological variables measured were relative humidity (RH), air temperature (T air ) (HMP45AC, Vaisala Inc., Helsinki, Finland) and net radiation (Q*; CNR1-L, Kipp and Zonen, Delft, The Netherlands), and precipitation was measured with a tipping bucket rain gauge (TE525-L, Texas Electronic, Dallas, TX, USA).
The flux measurements in the MS site were performed using a three-dimensional anemometer (CSAT3, Campbell Scientific, Logan, UT, USA) and an open-path gas analyzer (EC150, Campbell Scientific, Logan, UT, USA). Additionally, meteorological variables included a temperature-humidity probe sensor for relative humidity (RH) and air temperature (T air ) (CS215-L, Campbell Scientific, Logan, UT, USA), and net radiation was estimated with a four-component radiometer (Q*; CNR4-L, Campbell Scientific Inc., Logan, UT, USA). Precipitation was measured with a tipping bucket rain gauge (TE522MM, Texas Electronic, Dallas, TX, USA).
Finally, the EC system for the OG site consisted of a three-dimensional sonic anemometer (Windmaster, Gill Instruments, Lymington, UK) and an infrared gas analyzer (LI-7500A, LI-COR Biosciences, Lincoln, NE, USA). These measurements were made at 10 Hz, storing the data through an analyzer interface unit (LI-7550, LI-COR Biosciences, Lincoln, NE, USA). Ancillary, meteorological variables were collected in parallel. These included a temperature-humidity probe sensor for relative humidity (RH) and air temperature (T air ) (HMP45C-L, Vaisala Inc., Helsinki, Finland) and a four-component radiometer to estimate net radiation (Q*; CNR4-L, Kipp and Zonen, Delft, The Netherlands), and precipitation was measured with a tipping bucket rain gauge (TE525-L, Texas Electronic, Dallas, TX, USA).
At each site, soil volumetric water content was measured with water content reflectometers (CS616, Campbell Scientific Inc., Logan, UT, USA), and soil temperature was estimated with a thermocouple probe (TCAV-L, Campbell Scientific Inc., Logan, UT, USA). All meteorological data in the sites were collected and stored in a data logger and averaged to 30 min intervals, except for precipitation, which was accumulated. A further and detailed description of the main layout for EC instrumentation and measurements can be found in Rojas-Robles et al. [70].
Corrections and quality controls of EC data were carried out through a filtering of the friction speed coefficient (u*) for the periods of low turbulence and instrument failure; the data that did not meet these criteria were discarded. Missing data were filled, and then NEE was partitioned into its components of GPP and R eco using the REddyProc based on the procedures described by Reichstein et al. [71] and Wutzler et al. [72]. Additional information about the data quality control, flux partitioning and gap-filling is reported in Rojas-Robles et al. [70].

Statistical Analyses
For all data analyzed, a Kolmogorov-Smirnov test with Lilliefors correction was used to assess normality. Then, we explored the annual and seasonal dynamics of R S , TS, SWC and PPT. Monthly means of R S , TS and SWC were calculated from the replicated sampling points available in the transects at each site and a Kruskal-Wallis test following by a Dunn multiple comparison post hoc test (p < 0.05) was used to determine the differences of R S across study sites between dry and wet seasons. Additionally, annual C emission estimates were calculated for each soil collar using a time course fitting and area under the curve integration function (SigmaPlot v.12, Systat Software, Inc., Chicago, IL, USA) [73]. Finally, for ecosystem flux and meteorological variables, Spearman pairwise correlations were used to identify the relationships among all the variables and a Kruskal-Wallis test with Dunn multiple comparison post hoc test (p < 0.05) was used to determine the differences among sites.
Structural equation modeling (SEM) is a multivariate statistical tool that allows testing causal or correlative relationships among different variables. The general SEM model considered a full set of hypotheses based on literature, precedent exploratory analyses and our own experience [11,70]. Briefly, we hypothesize that the global (K) and net (Q*) radiation would control the evapotranspiration (ET, negatively), vapor pressure deficit (VPD, positively), and air temperature (T air, positively) regardless of the successional stage [74]. In turn, the VPD and ET would be correlated to gross primary production (GPP, negatively and positively, respectively) at all sites [75]. GPP would be controlled by either K or Q*, depending on the successional stage [76]. A higher ET reduces the residence time of soil moisture, therefore negatively affecting the soil water content (SWC) [77] and a cooling effect that constrains R S [78,79]. In consequence, SWC together with soil temperature (TS) would be the most important controlling of R S [18]. The variation of GPP across forest succession would be reflected on R S , since GPP represents the autotrophic activity of the ecosystem [27]. The general model was applied to each successional stage (i.e., old growth, mid-secondary or early-secondary).
The SEM model was tested for the combined years of R S and half-hour ecosystem flux and meteorological observations, which match with the dates and hours when R S measurements were carried out during the sampling campaigns (Table 2). Prior to the SEM analysis, the variables of interest from the dataset of each site were log transformed. Several models were run, and the best-fitted were selected according to the goodness of fit [80]; this procedure was performed to the dataset that integrates all three study sites (general model) and was also performed individually for each study site (i.e., OG, MS and ES). In all cases, to ensure a good model fitting, goodness of fit was assessed according to Grace [81] through the following parameters: (1) a Chi-square (χ 2 ) test-this parameter must be nonsignificant χ 2 (p > 0.05; df = 1); (2) a root-mean-square error of approximation statistic (RMSEA), whose values must be <0.10; and (3) a goodness-of-fit index (GFI) and the Bentler and Bonett's normed-fit index (NFI), where both estimates must be >0.9. The standardized path regression weights (SRW) were obtained with the maximum likelihood estimation method [82].
All SEM analyses were performed in the software IBM ® SPSS ® Amos™ (version 20.0, IBM Corporation, Armonk, NY, USA), while the remaining statistical tests were performed in Statistica (version 7, Stat Soft Inc., Tulsa, OK, USA).

Precipitation
The mean historical annual precipitation for study sites is 724 mm year −1 [61]. Overall, rainfall recorded for study sites denote that 2015 was above average with 805 mm year  (Figure 1a).

Soil Respiration
The temporal variations of R S in all three sites are displayed in Figure 1c, which depicts a strong seasonality during the study period across all sites. R S showed a relationship with precipitation patterns (Figure 1d), although the seasonal responses and magnitudes at sites were different (Table 3). R S magnitudes during the dry season were similar across sites with minor responses to precipitation events during the winter and early wet season. During the wet season, R S varied from 0.06 to 6.25 µmol CO 2 m −2 s −1 in the OG site, while for MS and ES sites it ranged from 0.05 to 4.25 µmol CO 2 m −2 s −1 and from 0.10 to 4.25 µmol CO 2 m −2 s −1 , respectively. As expected for seasonally dry forests, the peak of R S occurred in the wet season (June-October). During this period, R S rates in the OG site were consistently higher than in the MS and ES sites, reaching a mean annual rate of 5.20 ± 0.25 µmol CO 2 m −2 s −1 , while for the MS and ES sites the mean annual rates were 4.52 ± 0.26 µmol CO 2 m −2 s −1 and 4.56 ± 0.32 µmol CO 2 m −2 s −1 , respectively (Table 3). In 2015, the R S peak for the OG sites was in October, whereas for the MS and ES sites it was in August; in 2016 and 2017 all sites presented their peaks in August and July, respectively. By the end of the wet season (October) of 2015 and 2016, R S rates in the OG site were sustained for a longer period, while those in the MS and ES sites decreased soon after the moisture availability in the system halted ( Figure 1c).
Annual C emission estimates throughout the study period in all the sites are displayed in Table 4. The mean annual budget for the entire observation years were 868.80 ± 70.95 g C m −2 in the OG, followed by 731.96 ± 30.40 g C m −2 for the ES site and 718.66 ± 48.54 g C m −2 in the MS site. Meanwhile, the main contribution of the seasonal C emissions to annual budgets was in the wet season, which ranged from 54% to 67% in the OG site, while for the MS site it was from 57% to 72%, and from 54% to 70% in the ES site.

Ecosystem Flux and Meteorological Variations
Ecosystem flux and meteorological means during the study period showed significant differences among the three sites (p < 0.05, Table 5). For example, K, VPD and ET showed significant difference among sites with the greatest values in the ES site (p < 0.05). VPD showed a decreasing pattern from the ES to MS site, followed by OG, while ET and K presented a gradient from the ES to the OG site, followed by MS. Conversely, GPP and T air at the MS and ES sites were lower than OG (p < 0.05), whereas Q* showed higher values at MS followed by OG and ES. Moreover, TS and SWC showed significant differences along the forest succession sites (p < 0.05); both variables showed a decreasing pattern from ES > MS > OG. The SWC recorded at ES was significantly higher than at MS and OG. TS was significantly lower at OG, and no significant differences were shown between ES and MS (p < 0.05). Table 5. Annual means of global radiation (K), net radiation (Q*), air temperature (T air ), vapor deficit pressure (VPD), evapotranspiration (ET), gross primary production (GPP), soil volumetric water content (SWC) and soil temperature (TS) along a forest succession from tropical dry forest in northwestern Mexico. Data are presented as mean ± standard error. Capital letters denote significant differences among sites (p < 0.05).

Relationship between R S to Flux and Biophysical Controls
Using Spearman correlations, our results identified relationships from biophysical controls with R S (Table 6); where VPD was negatively correlated (p < 0.0001) in all sites K showed a negative relationship only in the OG site. Meanwhile, ET and GPP were only correlated in the OG and MS sites, and T air showed a significant positive relationship at MS (p < 0.05) and ES (p < 0.001) sites but not in the OG site. Furthermore, a positive correlation between SWC site and R S was found in the OG and MS sites, while TS site did not show a relationship in any of the three sites. Table 6. Correlation coefficients of soil respiration with global radiation (K), net radiation (Q*), air temperature (T air ), vapor deficit pressure (VPD), evapotranspiration (ET), gross primary production (GPP), soil volumetric water content (SWC) and soil temperature (TS) in old-growth (OG), midsecondary (MS) and early-secondary (ES) tropical dry forests in northwestern Mexico. The "*", "**" and "***" represent the 0.0001, 0.001 and 0.05 significance probability of p-value, respectively. In turn, soil respiration (R S ) was consistently controlled by TS site , having a positive and direct relationship at all sites; however, this effect had a stronger significant correlation in the ES (0.22) and MS (0.13) sites in contrast to OG (0.10). Then, a positive and direct relationship of SWC site and ET was identified at the OG and MS sites. However, we found that the tightest effect of SWC site over R S tended to be in the OG site (0.23), with the least being in the MS forest (0.13); meanwhile, the relationship with ET was highly significant in both sites (OG = 0.25, MS = 0.27). On the contrary, in the ES site, SWC site and ET did not have a significant effect over R S (Figure 2c). Additionally, K , seems to have a minor direct effect over R S , only affecting it negatively in the OG site (Figure 2a). Notably, the SEM results showed that GPP exerted an overall direct effect over R S in the ES and OG sites, but not in the MS site. A strong control of GPP was identified in the ES site (0.70), in contrast to the OG site, where the least effect was detected (0.41), and a null effect was observed in the MS site (Figure 2b).  . Each site path diagram shows the Bentler and Bonett's normed-fit index (NFI), the goodness-of-fit index (GFI), the root-mean-square error of approximation statistic (RMSEA), a Chi-square (χ 2 ) value, the significance probability of p-value, the degrees of freedom (d.f.) and sample size (n). Arrow widths are proportional to p-value scheme.

Discussion
In this manuscript, we analyzed the temporal variations of soil respiration (R S ) in order to determine intra-annual and seasonal variations along three forest succession sites in a tropical dry forest (TDF) in northwestern Mexico by incorporating the analysis of R S measurements, continuous ecosystem flux and meteorological observations. Furthermore, we assessed the relative contribution of biophysical controls of R S relaying in structural equation models (SEM). We discuss the interactions and the influence of different environmental and biological drivers that regulate the interannual variability of R S , such as soil moisture, soil temperature, atmospheric moisture demand or the autotrophic contribution and their relative importance along a TDF succession in northwestern Mexico.

Interannual and Seasonal Variation of R S along Succession in TDF
In seasonally dry ecosystems such as TDF, carbon, water and energy fluxes are tightly linked to the intra-annual and interannual seasonal variations in precipitation and moisture availability [83] and successional changes in the forest affecting R S [56,84,85].
In the TDF of northwestern Mexico, R S followed a typical seasonal pattern for seasonally dry ecosystems where after a period of water stress and dormancy, a large CO 2 pulse is stimulated by the first precipitation events. Soil respiration then reaches a maximum rate coinciding with the highest overall biological activity of the TDF, and finally a drastic decrease as the soil dries towards the end of the growing season [16,56,86]. We hypothesized that R S would have significant variations across forest succession sites at annual and seasonal scales. The mean annual R S rate in the old-growth (OG) was higher than at the mid-secondary (MS) and early-secondary (ES) forests, but contrary to our expectations, during the study period, the mean annual R S rates did not reveal significant differences between the sites. However, at interannual scales, R S mean rates differed significantly between seasons and according to forest succession (Table 3). These temporal R S variations across sites suggest that these seasonal differences could be related to site-specific direct and indirect interactions between biotic and abiotic factors that include environmental conditions and vegetation attributes (i.e., root density, above-and belowground biomass) coupled with physical and chemical soil properties, substrate availability and microbial activity along the forest succession [87][88][89].
Noteworthy, we found that R S responded to precipitation inputs differently according to the state of forest succession, where the sites on early-and mid-secondary successional stages depict a faster increasing response in this flux than an old-growth forest. However, the ES and MS sites resulted in a much more ephemeral flux with a rapidly decreasing trend at the end of the wet season, while the OG site showed sustained high R S rates for a longer time during the summer-fall transition (Figure 1c). A likely explanation for the sustained R S rates in the OG site may be related to the amounts of soil organic matter (5.5 ± 0.2%) and litterfall production (4.32 Mg ha −1 year −1 ) from previous seasons in this site. Furthermore, accumulated organic matter still available for decomposition from tree mortality that occurred in 2011 due to an extreme freezing event in the region may play an additional role [61,90].

Responses of R S to Biophysical Controls
In this study, structural equation models (SEM) were applied to integrate multiyear observations of R S, ecosystem flux and meteorological variables, in order to identify the main biophysical controls of the variability of R S and their relative importance along a three TDF succession sites. Overall, we hypothesized that the relative importance of biophysical controls in the interannual variability of R S , such as global radiation (K), air temperature (T air ), vapor pressure deficit (VPD), evapotranspiration (ET) and gross primary production (GPP), would be different among OG, MS and ES TDF sites.
At the TDF of northwestern Mexico, R S is mainly driven by soil volumetric water content (SWC site ), and soil temperature (TS site ) represents a lesser role in this flux, similar to other studies in seasonally dry ecosystems [91][92][93][94].
According to our results, the influence of SWC site was higher in the OG than the MS site; even though SWC site and ET showed the highest mean annual values in the ES site during the study period, there was no significant effect on R S from these controls ( Figure 2). The strong influence of SWC site observed in the OG and MS sites is consistent with observations in other forests and seasonally dry ecosystems [11,95,96].
Contrary to our expectations, R S was consistently influenced by TS site , regardless of forest succession on the TDF sites. However, the relative importance over R S was different along the forest successions since the influence of TS site on the ES site was higher and decreased as the forest transitioned to more advanced successional stages at MS and OG and as observed in other types of forests [89].
The contrasting effects of GPP on R S in our TDF succession gradient suggests that legume pioneer tree species may be exerting a strong influence on R S in the ES site due to a high metabolic activity (high GPP) as a consequence of a fast canopy and root system development during these early stages of succession in TDF [97]. Pioneering trees species maintain an active metabolism under water-stressed and drought conditions to maximize their net photosynthetic rates and resource acquisition efficiency [98][99][100][101]. At the TDF in northwest Mexico, the pioneer species A. cochliacantha forms a dense canopy and displays functional traits for high and rapid water use, such as an extensive lateral shallow-rooted system with no tap roots [52,102]. This trait combination suggests an active autotrophic activity in the ES site that is reflected in a strong positive effect of GPP on R S in our SEM analysis (SRW = 0.70, Figure 2c). It is worth mentioning that the effect of GPP over R S in the MS site is absent (Figure 2b), having been significant in the OG site but not as strong as in the ES forest (SRW = 0.41, Figure 2a). The complexity of the rooting systems might explain these contrasting effects of GPP over R S along the TDF succession as vegetation structure and composition changes with forest succession [103,104]. In the mid-secondary forest (MS), the dominance of A. cochliacantha decreases as tree species with deeper root systems increase [52]. In contrast, the OG site presents mature vegetation and a dense herbaceous and shrub understory with a well-developed shallow-rooted system that might contribute to the effect of GPP over R S .
Other biophysical controls should be considered in order to improve our understanding in mechanisms involving the variability of R S . First, vapor pressure deficit (VPD) and gross primary production (GPP) were negatively correlated at all sites ( Figure 2). Gas exchange and carbon isotope studies have discussed the link between tree photosynthesis, respiration and VPD [105][106][107][108][109]. VPD is considered as an ecophysiological moisture stress proxy of vegetation for regulating stomata conductance that contributes to determining the carbon uptake and water loss dynamics in ecosystems [110]; this suggests that photosynthetic rates and stomatal conductance at high VPD values constrain GPP [111], and therefore should have an influence on root respiration [108,112]. This is demonstrated in our results, as soil respiration (R S ) in the MS site showed a positive correlation with T air and was negatively correlated with VPD in the MS and ES sites (Figure 2), but did not influence R S in the OG site. The strongest effect of VPD in the MS site suggests a higher evaporative demand in the soil [113], resulting in a shorter water time residence in the shallow soil, and thus might be constraining the metabolic activity from the heterotrophic component of R S at this site.
Secondly, global radiation (K) exerted a minor direct control in R S with a negative effect in the OG site only. Contrary to early and secondary forests, old-growth forests present a well-developed canopy and understory with a low canopy openness, a high foliage density and a greater leaf area index [114,115]. These features suggest that the influence of K over R S in the OG forest is due to more drastic seasonal changes in vegetation cover across the wet and dry seasons [116]. For example, a lower leaf area index as a consequence of litterfall during the dry season [117], affecting forest soil microclimate [118,119].

Relatability of SEM to Assess Cause-Effects on R S
Structural equation models (SEM) provide a scientific framework to represent direct and indirect cause-effect relationships from highly related variables and test multivariate hypotheses about multiple ecological processes [81,120]. The SEM presented in this study shows the multiple causal relationships that control the variability of R S from a tropical dry forest (TDF) in northwestern Mexico and allowed us to discern the relative importance among biophysical controls at individual sites and across the succession gradient. These models underline how the paths change among the TDF successional stages, highlighting the greater complexity of interactions in an old-growth forest (OG) compared to early-(ES) and mid-secondary (MS) sites.
Furthermore, the path model performance for the OG, MS and ES sites showed a good fit since models explained a high proportion of the variance of R S (ES, R 2 = 0.76; OG, R 2 = 0.72; MS, R 2 = 0.68). These estimates are comparable with previous studies in temperate, tropical, subtropical, arid and semiarid ecosystems that tested significant cause-effect interactions among different biophysical controls and R S . For example, Guan et al. [121] observed R 2 = 0.79 in a steppe, and Flores-Rentería et al. [82] explained R 2 = 0.17 in a Mediterranean forest; meanwhile, for a humid subtropical forest, Tian et al. [122] reported R 2 = 0.63 and Li et al. [123] showed R 2 = 0.42 and R 2 = 0.39 for tropical and subtropical forests and grasslands, respectively. In addition, Campuzano et al. [11] reported R 2 = 0.50 in a semiarid shrubland. Therefore, SEM is a robust modeling tool to account for the effects of key controls of R S variability and create a functional framework depicting the interactions, mechanisms and processes in TDF soils from northwestern Mexico [124].

Soil Respiration from Tropical Dry Forests in a Global Context
Dryland ecosystems cover about 40% of the terrestrial surface [125]. These ecosystems have received a significant role in the global carbon cycle (C) due to dynamic biogeochemical processes that define their productivity and control the global CO 2 concentration in the atmosphere [126]. At the global scale, they represent 51% and 39% of the global land C sink and the net C flux interannual variability [127], respectively, and nearly three-quarters of C stocks are stored in soil [128].
Tropical dry forests (TDF) in northwestern Mexico can contribute an important amount of CO 2 emissions from soils. We found that it is possible to compare our seasonal R S rates across several drylands and seasonally dry ecosystems. For instance, Leon et al. [92] reported R S variations for Mediterranean shrubland from 0.4 to 0.8 µmol CO 2 m −2 s −1 during the dry season, and during the wet season these rates ranged from 0.9 to 2.1 µmol CO 2 m −2 s −1 ; Arellano-Martin et al. [86] for TDFs in southeast Mexico found a mean rate of 1.5 µmol CO 2 m −2 s −1 and 5.3 µmol CO 2 m −2 s −1 during the dry and wet season, respectively. Meanwhile, for a TDF in Thailand, Adachi et al. [129] reported mean R S rates from 1.98 to 3.40 µmol CO 2 m −2 s −1 in the dry season and from 3.20 to 4.49 µmol CO 2 m −2 s −1 in the wet season; Hanpattanakit et al. [130] observed seasonal ranges from 1.93 to 2.20 µmol CO 2 m −2 s −1 in dry periods, whereas in the wet periods it varied from 2.70 to 3.93 µmol CO 2 m −2 s −1 . Regarding R S in TDF sites with succession, during the dry season, Calvo-Rodriguez et al. [56] found mean rates of 0.75 ± 0.85 µmol CO 2 m −2 s −1 for an abandoned pasture (AP), 0.72 ± 0.28 µmol CO 2 m −2 s −1 in the early-secondary forest (ES) and 1.26 ± 0.79 µmol CO 2 m −2 s −1 for a mid-secondary forest (MS), then mean rates reported at the wet season were 4.15 ± 2.25 µmol CO 2 m −2 s −1 , 6.85 ± 2.82 µmol CO 2 m −2 s −1 and 6.53 ± 2.56 µmol CO 2 m −2 s −1 for AP, ES and MS, respectively.
Finally, our annual cumulative R S estimates included ranges from 712.87 to 1005.89 g C m −2 in the OG site, while those in the MS site ranged from 618.16 to 837.38 g C m −2 and from 681.19 to 786.33 g C m −2 in the ES site, which are consistent and within the ranges reported in other tropical, subtropical and seasonally dry forests around the world (Table 7).

Conclusions
Seasonality plays a critical role in the variability of R S along forest succession in TDFs from northwestern Mexico. The patterns of R S responded differently to precipitation inputs according to the forest succession, where the early-(ES) and mid-secondary (MS) forests showed a faster response at the early wet season, but with an ephemeral flux with rapid decreasing towards the end of the wet season compared to the old-growth forest (OG). Mean annual rates of R S did not show significant differences among sites, although the OG forest sustained the highest rates along the study period.
Our results revealed through structural equation models (SEM) that as the forest transitions to advanced successional stages, the path interactions between R S and the biophysical factors become more complex. We found a variable control of GPP on R S across succession, suggesting that the complexity of the rooting systems plays an important role on the soil respiration flux of the TDF.
Finally, soils of the TDF in northwestern Mexico showed annual net emissions between 618 to 1006 g C m −2 , varying according to environmental conditions and forest succession but within the range observed in other tropical dry forest.

Informed Consent Statement: Not applicable.
Data Availability Statement: Data will be available on request.