Methods and Results

The visual inspection of the residuals distribution shows that for MAP (Figure 1a) it is normaly distributed, with the mean near zero. For the MSI, residuals presented a tailed distribution and to deal with heterocedasticity of residuals we used the function exp to transform the data (Figure 1b).

Figure 1: Residuals inspection.

Figure 1: Residuals inspection.

MAP

Firstly, we used the gls function from the R package nlme (Pinheiro et al., 2019) to fit a linear mixed-effects model. This model considered MAP (and MSI) as the response variable, with the factors of interest being Biome (with two categories: Amazon, Caatinga and Cerrado) and Vegetation Type (with three categories: Savanna, Gallery Forest, and Dry Forest) (Table 1).

m1 <- nlme::gls(Rainfall ~ as.factor(Biome),
    data = data,
    weights = nlme::varIdent(form = ~ 1 | as.factor(Biome)),
    method = "ML", na.action = na.omit) 

summary(m1)$tTable %>% knitr::kable()
Value Std.Error t-value p-value
(Intercept) 2268.5796 1.053557 2153.2581 0
as.factor(Biome)Caatinga -1577.1555 1.445449 -1091.1179 0
as.factor(Biome)Cerrado -832.1394 1.516494 -548.7260 0
as.factor(Biome)Dry Forest -972.7419 19.925029 -48.8201 0
as.factor(Biome)Gallery Forest -867.0943 2.885341 -300.5171 0
as.factor(Biome)Savanna -859.5753 2.113434 -406.7198 0

Next, we performed an analysis of variance (ANOVA) using the anova function in R to assess the statistical significance of the main effects of factors.

an_m <- anova(m1)
an_m %>% data.frame %>% knitr::kable()
numDF F.value p.value
(Intercept) 1 6535985.7 0
as.factor(Biome) 5 238171.8 0

ANOVA revealed statistically significant differences (Table 2) and we proceeded with post hoc tests to conduct multiple comparisons between groups. For this purpose, we performed pairwise comparisons between factors with the emmeans function from the package emmeans (Lenth, 2022) to compute marginal means and confidence intervals for each combination of factors (Table 3).

means <- emmeans::emmeans(m1, specs = list("Biome"))
means %>% data.frame %>% knitr::kable()
Biome emmean SE df lower.CL upper.CL
Amazon 2268.5796 1.0535568 138088.10590 2265.8000 2271.3592
Caatinga 691.4241 0.9896169 28495.18816 688.8131 694.0352
Cerrado 1436.4402 1.0907662 59247.67705 1433.5624 1439.3180
Dry Forest 1295.8377 19.8971552 12.38576 1233.4788 1358.1965
Gallery Forest 1401.4852 2.6861143 11.10997 1392.8858 1410.0847
Savanna 1409.0042 1.8321083 12.08495 1403.2359 1414.7726

Subsequently, we conducted multiple comparison tests using the contrast function, also from emmeans package, with the argument method = “pairwise”. This test allowed us to identify significant differences between groups (or not) by comparing all factor levels, such as biomes and vegetation types. Adjusted p-values were considered to control the family-wise error rate by Tukey method (see contrast-methods from the package emmeans).

contrasts <- emmeans::contrast(means, method = "pairwise")

contrasts %>%
  data.frame %>%
  rename(Pair = "emmeans.of.Biome.contrast") %>% 
  rename(estimate = "emmeans.of.Biome.estimate") %>% 
  rename(se = "emmeans.of.Biome.SE") %>% 
  rename(df = "emmeans.of.Biome.df") %>% 
  rename(`t ratio` = "emmeans.of.Biome.t.ratio") %>% 
  rename(`p-value` = "emmeans.of.Biome.p.value")  %>% knitr::kable()
Pair estimate se df t ratio p-value
Amazon - Caatinga 1577.155475 1.445449 102390.78352 1091.117868 0.0000000
Amazon - Cerrado 832.139410 1.516494 160995.80547 548.725977 0.0000000
Amazon - Dry Forest 972.741893 19.925029 12.45527 48.820100 0.0000000
Amazon - Gallery Forest 867.094348 2.885341 14.79007 300.517127 0.0000000
Amazon - Savanna 859.575345 2.113434 21.39806 406.719782 0.0000000
Caatinga - Cerrado -745.016065 1.472791 81755.45753 -505.853311 0.0000000
Caatinga - Dry Forest -604.413582 19.921750 12.44713 -30.339382 0.0000000
Caatinga - Gallery Forest -710.061127 2.862613 14.33063 -248.046519 0.0000000
Caatinga - Savanna -717.580130 2.082297 20.16873 -344.609824 0.0000000
Cerrado - Dry Forest 140.602483 19.927031 12.46039 7.055867 0.0001242
Cerrado - Gallery Forest 34.954938 2.899135 15.07896 12.057025 0.0000000
Cerrado - Savanna 27.435935 2.132227 22.16928 12.867268 0.0000000
Dry Forest - Gallery Forest -105.647545 20.077649 12.71369 -5.261948 0.0017983
Dry Forest - Savanna -113.166548 19.981327 12.59476 -5.663615 0.0009654
Gallery Forest - Savanna -7.519003 3.251435 18.00246 -2.312518 0.2393606

We found differences in MAP between dry forests and both gallery forests and savannas (Table 4). Unexpectedly, the MAP of our samples are distinct from the MAP of both Biomes, including the Cerrado where they are located (Table 4).

Furthermore, we found that the MAP for more than half of all dry forests samples are outside the bistability zone (1300-2100 mm), i.e. the range of MAP in which forests and savannas may occur, defined by Ciemer et al. (2019) for the neotropics (Figure 2). However they are near the lower limit. Despite the unexpected differences in MAP for dry forests, most of our samples are inside the bistability zone for the neotropics (1300-2100mm; Ciemer et al., 2019) (Figure 2). By the other hand, Staver et al. (2011) defined a more wide range (1000-2500 mm), but using a global model. Thus, the limits defined by Staver et al. (2011) correctly includes dry forests into the bistability zone.

Figure 2: MAP distribution for Amazon, Caatinga and Cerrado Biomes, and for our sampled vegetation types. Red lines are the bistability zone defined by Staver et al. (2011).

Figure 2: MAP distribution for Amazon, Caatinga and Cerrado Biomes, and for our sampled vegetation types. Red lines are the bistability zone defined by Staver et al. (2011).

Markham Seasonality Index (MSI)

To compare between vegetation types and Biomes regardind MSI (i.e. seasonality lenght) (Markham, 1970), we performed the same method used to test for differences in MAP (Tables 5-8).

m2 <- nlme::gls(exp(Rainfall) ~ as.factor(Biome),
    data = data1,
    weights = nlme::varIdent(form = ~ 1 | as.factor(Biome)),
    method = "ML", na.action = na.omit) 

summary(m2)$tTable %>% knitr::kable()
Value Std.Error t-value p-value
(Intercept) 1.4846488 0.0004642 3197.9982 0
as.factor(Biome)Caatinga 0.3037148 0.0014085 215.6307 0
as.factor(Biome)Cerrado 0.2894743 0.0005671 510.4890 0
as.factor(Biome)Dry Forest 0.3248452 0.0031265 103.8992 0
as.factor(Biome)Gallery Forest 0.3262108 0.0012231 266.7144 0
as.factor(Biome)Savanna 0.3295360 0.0010822 304.5066 0

The ANOVA revealed significative differences among groups (Table 6).

an_m <- anova(m2)
an_m %>% data.frame %>% knitr::kable()
numDF F.value p.value
(Intercept) 1 47773529.45 0
as.factor(Biome) 5 59169.47 0

The range of values of MSI is very similar among groups. However, the range of MSI values for the Amazon biome clearly differs from the other groups and do not overlap with Cerrado and Caatinga nor with our samples (Table 7; Figure 3).

means <- emmeans::emmeans(m2, specs = list("Biome"))
means %>% data.frame %>% knitr::kable()
Biome emmean SE df lower.CL upper.CL
Amazon 1.484649 0.0004642 1.336860e+05 1.483424 1.485874
Caatinga 1.788364 0.0013298 2.444528e+04 1.784855 1.791872
Cerrado 1.774123 0.0003256 5.886448e+04 1.773264 1.774982
Dry Forest 1.809494 0.0030919 8.823026e+00 1.799037 1.819951
Gallery Forest 1.810860 0.0011315 6.742534e+00 1.806689 1.815030
Savanna 1.814185 0.0009776 7.525629e+00 1.810718 1.817652

The comparisons among groups showed differences between biomes and between biomes and all sampled vegetation types (Table 8). Dry forests, gallery forests and savannas at the PNCV do not differ in terms of length of dry season.

contrasts <- emmeans::contrast(means, method = "pairwise")

contrasts %>%
  data.frame %>%
  rename(Pair = "emmeans.of.Biome.contrast") %>% 
  rename(estimate = "emmeans.of.Biome.estimate") %>% 
  rename(se = "emmeans.of.Biome.SE") %>% 
  rename(df = "emmeans.of.Biome.df") %>% 
  rename(`t ratio` = "emmeans.of.Biome.t.ratio") %>% 
  rename(`p-value` = "emmeans.of.Biome.p.value")  %>% knitr::kable()
Pair estimate se df t ratio p-value
Amazon - Caatinga -0.3037148 0.0014085 3.087989e+04 -215.630689 0.0000000
Amazon - Cerrado -0.2894743 0.0005671 1.921352e+05 -510.488963 0.0000000
Amazon - Dry Forest -0.3248452 0.0031265 9.225340e+00 -103.899241 0.0000000
Amazon - Gallery Forest -0.3262108 0.0012231 9.202904e+00 -266.714364 0.0000000
Amazon - Savanna -0.3295360 0.0010822 1.130368e+01 -304.506613 0.0000000
Caatinga - Cerrado 0.0142405 0.0013691 2.744074e+04 10.401580 0.0000000
Caatinga - Dry Forest -0.0211304 0.0033657 1.238590e+01 -6.278118 0.0003959
Caatinga - Gallery Forest -0.0224959 0.0017461 3.841572e+01 -12.883842 0.0000000
Caatinga - Savanna -0.0258212 0.0016504 6.115541e+01 -15.645018 0.0000000
Cerrado - Dry Forest -0.0353709 0.0031090 9.019871e+00 -11.377016 0.0000123
Cerrado - Gallery Forest -0.0367365 0.0011775 7.904815e+00 -31.199792 0.0000000
Cerrado - Savanna -0.0400618 0.0010304 9.287198e+00 -38.881078 0.0000000
Dry Forest - Gallery Forest -0.0013656 0.0032924 1.203992e+01 -0.414754 0.9980062
Dry Forest - Savanna -0.0046908 0.0032427 1.096499e+01 -1.446566 0.7012329
Gallery Forest - Savanna -0.0033253 0.0014953 1.135516e+01 -2.223782 0.2982834

MSI do not differ among our sample plots, but they differ biomes (Table 8). However, the lenght of the dry season (i.e MSI) of the sampled vegetation types is in the range of values found for the Cerrado Biome (Table 8; Fig. 3), as they are located within the biome.

Figure 3: MSI distribution for Amazon, Caatinga and Cerrado Biomes, and for our sampled vegetation types.

Figure 3: MSI distribution for Amazon, Caatinga and Cerrado Biomes, and for our sampled vegetation types.

Conclusion

Gallery forests and savannas do not differ in MAP, but dry forests differ from them. Despite differences showed in the analysis of MAP for dry forests, the values are very close compared with the range of values for the Cerrado biome. Furthermore, all vegetation types are inside the bistability zone for MAP defined by Staver et al. (2011). Vegetation types do not differ in dry season length.

REFERENCES

Ciemer, C., Boers, N., Hirota, M., Kurths, J., Müller-Hansen, F., Oliveira, R.S., & Winkelmann, R. (2019). Higher resilience to climatic disturbances in tropical vegetation exposed to more variable rainfall. Nat. Geosci. 12, 174–179. https://doi.org/10.1038/s41561-019-0312-z

Markham, C. G. (1970). Seasonality of precipitation in the United States. Annals of the Association of American Geographers, 60(3), 593–597. https://doi.org/10.1111/J.1467-8306.1970.TB00743.X

Pinheiro J, Bates D, DebRoy S, Sarkar D, R Core Team (2019). nlme: Linear and Nonlinear Mixed Effects Models. R package version 3.1-140, <URL: https://CRAN.R-project.org/package=nlme>

Lenth, R.V. (2022). emmeans: Estimated Marginal Means, aka Least-Squares Means. R package version 1.7.4-1. https://CRAN.R-project.org/package=emmeans

Staver, A. C., Archibald, S., & Levin, S. A. (2011). The global extent and determinants of savanna and forest as alternative biome states. Science, 230–232(334). https://doi.org/10.1126/science.1210465