Critical transitions and the fragmenting of global forests

Forests provide critical habitat for many species, essential ecosystem services, and are coupled to atmospheric dynamics through exchanges of energy, water and gases. One of the most important changes produced in the biosphere is the replacement of forest areas with human dominated landscapes. This usually leads to fragmentation, altering the sizes of patches, the structure and function of the forest. Here we studied the distribution and dynamics of forest patch sizes at a global level, examining signals of a critical transition from an unfragmented to a fragmented state. We used MODIS vegetation continuous field to estimate the forest patches at a global level and defined wide regions of connected forest across continents and big islands. We search for critical phase transitions, where the system state of the forest changes suddenly at a critical point in time; this implies an abrupt change in connectivity that causes an increased fragmentation level. We studied the distribution of forest patch sizes and the dynamics of the largest patch over the last fourteen years, as the conditions that indicate that a region is near a critical fragmentation threshold are related to patch size distribution and temporal fluctuations of the largest patch. We found some evidence that allows us to analyze fragmentation as a critical transition: most regions followed a power-law distribution over the 14 years. We also found that the Philippines region probably went through a critical transition from a fragmented to an unfragmented state. Only the tropical forest of Africa and South America met the criteria to be near a critical fragmentation threshold. This implies that human pressures and climate forcings might trigger undesired effects of fragmentation, such as species loss and degradation of ecosystems services, in these regions. The simple criteria proposed here could be used as an early warning to estimate the distance to a fragmentation threshold in forest around the globe and a predictor of a planetary tipping point.


Introduction
Forests are one of the most important biomes on earth, providing habitat for a large proportion of species and contributing extensively to global biodiversity (Crowther et al., 2015). In the previous century human activities have influenced global bio-geochemical cycles (Bonan, 2008;Canfield et al., 2010), with one of the most dramatic changes being the replacement of 40% of Earth's formerly biodiverse land areas with landscapes that contain only a few species of crop plants, domestic animals and humans (Foley et al., 2011).
These local changes have accumulated over time and now constitute a global forcing (Barnosky et al., 2012).
Another global scale forcing that is tied to habitat destruction is fragmentation, which is defined as the division of a continuous habitat into separated portions that are smaller and more isolated. Fragmentation produces multiple interwoven effects: reductions of biodiversity between 13% and 75%, decreasing forest biomass, and changes in nutrient cycling (Haddad et al., 2015). The effects of fragmentation are not only important from an ecological point of view but also that of human activities, as ecosystem services are deeply influenced by the level of landscape fragmentation (Mitchell et al., 2015).
Ecosystems have complex interactions between species and present feedbacks at different levels of organization (Gilman et al., 2010), and external forcings can produce abrupt changes from one state to another, called critical transitions (Scheffer et al., 2009). These abrupt state shifts cannot be linearly forecasted from past changes, and are thus difficult to predict and manage (Scheffer et al., 2009). Such 'critical' transitions have been detected mostly at local scales (Drake & Griffen, 2010;Carpenter et al., 2011), but the accumulation of changes in local communities that overlap geographically can propagate and theoretically cause an abrupt change of the entire system at larger scales (Barnosky et al., 2012). Coupled with the existence of global scale forcings, this implies the possibility that a critical transition could occur at a global scale (Rockstrom et al., 2009;Folke et al., 2011).
Complex systems can experience two general classes of critical transitions (Solé, 2011). In so-called first order transitions, a catastrophic regime shift that is mostly irreversible occurs because of the existence of alternative stable states (Scheffer et al., 2001). This class of transitions is suspected to be present in a variety of ecosystems such as lakes, woodlands, coral reefs (Scheffer et al., 2001), semi-arid grasslands (Bestelmeyer et al., 2011), and fish populations (Vasilakopoulos & Marshall, 2015). They can be the result of positive feedback mechanisms (Martín et al., 2015); for example, fires in some forest ecosystems were more likely to occur in previously burned areas than in unburned places (Kitzberger et al., 2012).
The other class of critical transitions are continuous or second order transitions (Solé & Bascompte, 2006).
In these cases, there is a narrow region where the system suddenly changes from one domain to another, 3 with the change being continuous and in theory reversible. This kind of transitions were suggested to be present in tropical forest (Pueyo et al., 2010), semi-arid mountain ecosystems (McKenzie & Kennedy, 2012), tundra shrublands (Naito & Cairns, 2015). The transition happens at a critical point where we can observe a distinctive spatial pattern: scale invariant fractal structures characterized by power law patch distributions (Stauffer & Aharony, 1994).
There are several processes that can convert a catastrophic transition to a second order transitions (Martín et al., 2015). These include stochasticity, such as demographic fluctuations, spatial heterogeneities, and/or dispersal limitation. All these components are present in forest around the globe (Seidler & Plotkin, 2006;Filotas et al., 2014;Fung et al., 2016), and thus continuous transitions might be more probable than catastrophic transitions. Moreover there is some evidence of recovery in some systems that supposedly suffered an irreversible transition produced by overgrazing (Zhang et al., 2005) and desertification (Allington & Valone, 2010).
The spatial phenomena observed in continuous critical transitions deal with connectivity, a fundamental property of general systems and ecosystems from forests (Ochoa-Quintero et al., 2015) to marine ecosystems (Leibold & Norberg, 2004) and the whole biosphere (Lenton & Williams, 2013). When a system goes from a fragmented to a connected state we say that it percolates (Solé, 2011). Percolation implies that there is a path of connections that involves the whole system. Thus we can characterize two domains or phases: one dominated by short-range interactions where information cannot spread, and another in which long range interactions are possible and information can spread over the whole area. (The term "information" is used in a broad sense and can represent species dispersal or movement.) Thus, there is a critical "percolation threshold" between the two phases, and the system could be driven close to or beyond this point by an external force; climate change and deforestation are the main forces that could be the drivers of such a phase change in contemporary forests (Bonan, 2008;Haddad et al., 2015). There are several applications of this concept in ecology: species' dispersal strategies are influenced by percolation thresholds in three-dimensional forest structure (Solé et al., 2005), and it has been shown that species distributions also have percolation thresholds (He & Hubbell, 2003). This implies that pushing the system below the percolation threshold could produce a biodiversity collapse (Bascompte & Solé, 1996;Solé et al., 2004;Pardini et al., 2010); conversely, being in a connected state (above the threshold) could accelerate the invasion of forest into prairie (Loehle et al., 1996;Naito & Cairns, 2015).
One of the main challenges with systems that can experience critical transitions-of any kind-is that the value of the critical threshold is not known in advance. In addition, because near the critical point a small change can precipitate a state shift of the system, they are difficult to predict. Several methods have been developed to detect if a system is close to the critical point, e.g. a deceleration in recovery from perturbations, or an increase in variance in the spatial or temporal pattern (Hastings & Wysham, 2010;Carpenter et al., 2011;Boettiger & Hastings, 2012;Dai et al., 2012).
In this study, our objective is to look for evidence that forests around the globe are near continuous critical point that represent a fragmentation threshold. We use the framework of percolation to first evaluate if forest patch distribution at a continental scale is described by a power law distribution and then examine the fluctuations of the largest patch. The advantage of using data at a continental scale is that for very large systems the transitions are very sharp (Solé, 2011) and thus much easier to detect than at smaller scales, where noise can mask the signals of the transition.

Study areas definition
We analyzed mainland forests at a continental scale, covering the whole globe, by delimiting land areas with a near-contiguous forest cover, separated with each other by large non-forested areas. Using this criterion, we delimited the following forest regions. In America, three regions were defined: South America temperate forest (SAT), subtropical and tropical forest up to Mexico (SAST), and USA and Canada forest (NA). Europe and North Asia were treated as one region (EUAS). The rest of the delimited regions were South-east Asia (SEAS), Africa (AF), and Australia (OC). We also included in the analysis islands larger than 10 5 km 2 . The mainland region has the number 1 e.g. OC1, and the nearby islands have consecutive numeration (Appendix S4, figure S1-S6).

Forest patch distribution
We studied forest patch distribution in each defined area from 2000 to 2014 using the MODerate-resolution Imaging Spectroradiometer (MODIS) Vegetation Continuous Fields (VCF) Tree Cover dataset version 051 (DiMiceli et al., 2015). This dataset is produced at a global level with a 231-m resolution, from 2000 onwards on an annual basis. There are several definition of forest based on percent tree cover (Sexton et al., 2015); we choose a 30% threshold to convert the percentage tree cover to a binary image of forest and non-forest pixels. This was the definition used by the United Nations' International Geosphere-Biosphere Programme (Belward, 1996), and studies of global fragmentation (Haddad et al., 2015). This definition avoids the errors produced by low discrimination of MODIS VCF between forest and dense herbaceous vegetation at low forest cover (Sexton et al., 2015). Patches of contiguous forest were determined in the binary image by grouping connected forest pixels using a neighborhood of 8 forest units (Moore neighborhood).

Percolation theory
A more in-depth introduction to percolation theory can be found elsewhere (Stauffer & Aharony, 1994) and a review from an ecological point of view is available (Oborny et al., 2007). Here, to explain the basic elements of percolation theory we formulate a simple model: we represent our area of interest by a square lattice and each site of the lattice can be occupied-e.g. by forest-with a probability p. The lattice will be more occupied when p is greater, but the sites are randomly distributed. We are interested in the connection between sites, so we define a neighborhood as the eight adjacent sites surrounding any particular site. The sites that are neighbors of other occupied sites define a patch. When there is a patch that connects the lattice from opposite sides, it is said that the system percolates. When p is increased from low values, a percolating patch suddenly appears at some value of p called the critical point p c .
Thus percolation is characterized by two well defined phases: the unconnected phase when p < p c (called subcritical in physics), in which species cannot travel far inside the forest, as it is fragmented; in a general sense, information cannot spread. The second is the connected phase when p > p c (supercritical), species can move inside a forest patch from side to side of the area (lattice), i.e. information can spread over the whole area. Near the critical point several scaling laws arise: the structure of the patch that spans the area is fractal, the size distribution of the patches is power-law, and other quantities also follow power-law scaling (Stauffer & Aharony, 1994).
The value of the critical point p c depends on the geometry of the lattice and on the definition of the neighborhood, but other power-law exponents only depend on the lattice dimension. Close to the critical point, the distribution of patch sizes is: where n s (p) is the number of patches of size s. The exponent α does not depend on the details of the model and it is called universal (Stauffer & Aharony, 1994). These scaling laws can be applied for landscape structures that are approximately random, or at least only correlated over short distances (Gastner et al., 2009). In physics this is called "isotropic percolation universality class", and corresponds to an exponent α = 2.05495. If we observe that the patch size distribution has another exponent it will not belong to this universality class and some other mechanism should be invoked to explain it. Percolation thresholds can also be generated by models that have some kind of memory (Hinrichsen, 2000;Ódor, 2004): for example, a 6 patch that has been exploited for many years will recover differently than a recently deforested forest patch.
In this case, the system could belong to a different universality class, or in some cases there is no universality, in which case the value of α will depend on the parameters and details of the model (Corrado et al., 2014).
To illustrate these concepts, we conducted simulations with a simple forest model with only two states: forest and non-forest. This type of model is called a "contact process" and was introduced for epidemics (Harris, 1974), but has many applications in ecology (Solé & Bascompte, 2006;Gastner et al., 2009). A site with forest can become extinct with probability e, and produce another forest site in a neighborhood with probability c. We use a neighborhood defined by an isotropic power law probability distribution. We defined a single control parameter as λ = c/e and ran simulations for the subcritical fragmentation state λ < λ c , with λ = 2, near the critical point for λ = 2.5, and for the supercritical state with λ = 5 (see supplementary data, gif animations).

Patch size distributions
We fitted the empirical distribution of forest patch areas, based on the MODIS-derived data described above, to four distributions using maximum likelihood estimation (Goldstein et al., 2004;Clauset et al., 2009).
The distributions were: power-law, power-law with exponential cut-off, log-normal, and exponential. We assumed that the patch size distribution is a continuous variable that was discretized by remote sensing data acquisition procedure.
We set a minimal patch size (X min ) at nine pixels to fit the patch size distributions to avoid artifacts at patch edges due to discretization (Weerman et al., 2012). Besides this hard X min limit we set due to discretization, the power-law distribution needs a lower bound for its scaling behavior. This lower bound is also estimated from the data by maximizing the Kolmogorov-Smirnov (KS) statistic, computed by comparing the empirical and fitted cumulative distribution functions (Clauset et al., 2009). We also calculated the uncertainty of the parameters using a non-parametric bootstrap method (Efron & Tibshirani, 1994), and computed corrected Akaike Information Criteria (AIC c ) and Akaike weights for each model (Burnham & Anderson, 2002). Akaike weights (w i ) are the weight of evidence in favor of model i being the actual best model given that one of the N models must be the best model for that set of N models.
Additionally, we computed the goodness of fit of the power-law model following the bootstrap approach described by Clauset et. al (2009), where simulated data sets following the fitted model are generated, and a p-value computed as the proportion of simulated data sets that has a KS statistic less extreme than empirical data. The criterion to reject the power law model suggested by Clauset et. al (2009) was p ≤ 0.1, but as we 7 have a very large n, meaning that negligible small deviations could produce a rejection (Klaus et al., 2011), we chose a p ≤ 0.05 to reject the power law model.
To test for differences between the fitted power law exponent for each study area we used a generalized least squares linear model (Zuur et al., 2009) with weights and a residual auto-correlation structure. Weights were the bootstraped 95% confidence intervals and we added an auto-regressive model of order 1 to the residuals to account for temporal autocorrelation.

Largest patch dynamics
The largest patch is the one that connects the highest number of sites in the area. This has been used extensively to indicate fragmentation (Gardner & Urban, 2007;Ochoa-Quintero et al., 2015). The relation of the size of the largest patch S max to critical transitions has been extensively studied in relation to percolation phenomena (Stauffer & Aharony, 1994;Bazant, 2000;Botet & Ploszajczak, 2004), but seldom used in ecological studies (for an exception see Gastner et al. (2009)). When the system is in a connected state (p > p c ) the landscape is almost insensitive to the loss of a small fraction of forest, but close to the critical point a minor loss can have important effects (Solé & Bascompte, 2006;Oborny et al., 2007), because at this point the largest patch will have a filamentary structure, i.e. extended forest areas will be connected by thin threads. Small losses can thus produce large fluctuations.
One way to evaluate the fragmentation of the forest is to calculate the proportion of the largest patch against the total area (Keitt et al., 1997). The total area of the regions we are considering (Appendix S4, figures S1-S6) may not be the same than the total area that the forest could potentially occupy, thus a more accurate way to evaluate the weight of S max is to use the total forest area, that can be easily calculated summing all the forest pixels. We calculate the proportion of the largest patch for each year, dividing S max by the total forest area of the same year: This has the effect of reducing the S max fluctuations produced due to environmental or climatic changes influences in total forest area. When the proportion RS max is large (more than 60%) the largest patch contains most of the forest so there are fewer small forest patches and the system is probably in a connected phase. Conversely, when it is low (less than 20%), the system is probably in a fragmented phase (Saravia & Momo, 2017).
The RS max is a useful qualitative index that does not tell us if the system is near or far from the critical transition; this can be evaluated using the temporal fluctuations. We calculate the fluctuations around the mean with the absolute values ∆S max = S max (t) − ⟨S max ⟩, using the proportions of RS max . To characterize fluctuations we fitted three empirical distributions: power-law, log-normal, and exponential, using the same methods described previously. We expect that large fluctuation near a critical point have heavy tails (lognormal or power-law) and that fluctuations far from a critical point have exponential tails, corresponding to Gaussian processes (Rooij et al., 2013). As the data set spans 15 years, we do not have enough power to reliably detect which distribution is better (Clauset et al., 2009). To improve this we performed the goodness of fit test described above for all the distributions. We generated animated maps showing the fluctuations of the two largest patches to aid in the interpretations of the results.
A robust way to detect if the system is near a critical transition is to analyze the increase in variance of the density (Benedetti-Cecchi et al., 2015). It has been demonstrated that the variance increase in density appears when the system is very close to the transition (Corrado et al., 2014), thus practically it does not constitute an early warning indicator. An alternative is to analyze the variance of the fluctuations of the largest patch ∆S max : the maximum is attained at the critical point but a significant increase occurs well before the system reaches the critical point (Corrado et al., 2014). In addition, before the critical fragmentation, the skewness of the distribution of ∆S max should be negative, implying that fluctuations below the average are more frequent. We characterized the increase in the variance using quantile regression: if variance is increasing the slopes of upper or/and lower quartiles should be positive or negative.
All statistical analyses were performed using the GNU R version 3.3.0 (R Core Team, 2015), using code provided by Cosma R. Shalizi for fitting the power law with exponential cutoff model and the poweRlaw package (Gillespie, 2015) for fitting the other distributions. For the generalized least squares linear model we used the R function gls from package nlme (Pinheiro et al., 2016); and we fitted quantile regressions using the R package quantreg (Koenker, 2016). Image processing was done in MATLAB r2015b (The Mathworks Inc.). The complete source code for image processing and statistical analysis, and the patch size data files are available at figshare http://dx.doi.org/10.6084/m9.figshare.4263905.

Results
The figure 1 shows an example of the distribution of the biggest 200 patches for years 2000 and 2014, this distribution is highly variable, but the biggest patch usually maintains its place. Sometimes the biggest patch breaks and then big fluctuations in its size are observed, as we will analyze below. Smaller patches can merge or break more easily so they enter or leave the list of 200, and this is why there is a color change across years.
The power law distribution was selected as the best model in 92% of the cases (Appendix S4, Figure S7).
In a small number of cases (1%) the power law with exponential cutoff was selected, but the value of the parameter α was similar by ±0.02 to the pure power law. Additionally the patch size where the exponential tail begins is very large, thus we used the power law parameters for this cases (See Appendix S4, Figure S2, region EUAS3). In finite-size systems the favored model should be the power law with exponential cutoff, because the power-law tails are truncated to the size of the system (Stauffer & Aharony, 1994). Here the regions are so large that the cutoff is practically not observed.
There is only one region that did not follow a power law: Eurasia mainland, which showed a log-normal distribution, representing 7% of the cases. The log-normal and power law are both heavy tailed distributions and difficult to distinguish, but in this case Akaike weights have very high values for log-normal (near 1), meaning that this is the only possible model. In addition, the goodness of fit tests clearly rejected the power law model in all cases for this region (Appendix S4, table S1, region EUAS1). In general the goodness of fit test rejected the power law model in fewer than 10% of cases. In large forest areas like Africa mainland (AF1) or South America tropical-subtropical (SAST1), larger deviations are expected and the rejections rates are higher so the proportion is 30% or less (Appendix S4, Table S1).
Taking into account the bootstrapped confidence intervals of each power law exponent (α) and the temporal autocorrelation, there were no significant differences between α for the regions with the biggest (greater than 10 7 km 2 ) forest areas (Figure 2 and Appendix S4, Figure S8). There were also no differences between these regions and smaller ones (Appendix S4, Tables S2 & S3), and all the slopes of α were not different from 0 (Appendix S4, Table S3). This implies a global average α = 1.908, with a bootstrapped 95% confidence interval between 1.898 and 1.920.
The proportion of the largest patch relative to total forest area RS max for regions with more than 10 7 km 2 of forest is shown in figure 3. South America tropical and subtropical (SAST1) and North America (NA1) have a higher RS max of more than 60%, and other big regions 40% or less. For regions with less total forest area (Appendix S4, figure S9 & Table 1 and Java (OC7), all others were power laws (Appendix S4, Table S7).The goodness of fit test (GOF) did not reject power laws in any case, but neither did it reject the other models except in a few cases; this was  13 due to the small number of observations. We only considered fluctuations to follow a power law when this distribution was selected for both absolute and relative fluctuations.
The animations of the two largest patches (see supplementary data, largest patch gif animations) qualitatively shows the nature of fluctuations and if the state of the forest is connected or not. If the largest patch is always the same patch over time, the forest is probably not fragmented; this happens for regions with RS max of more than 40% such as AF2 (Madagascar), NA1 (North America), and OC3 (Malaysia). In the regions with RS max between 40% and 30% the identity of the largest patch could change or stay the same in time.
For OC7 (Java) the largest patch changes and for AF1 (Africa mainland) it stays the same. Only for EUAS1 (Eurasia mainland) we observed that the two largest patches are always the same, implying that this region is probably composed of two independent domains and should be divided in further studies. The regions with RS max less than 25%: SAST2 (Cuba) and EUAS3 (United Kingdom), the largest patch always changes reflecting their fragmented state. In the case of SEAS2 (Philippines) a transition is observed, with the identity of the largest patch first variable, and then constant after 2010.
The results of quantile regressions are identical for ∆RS max and ∆S max (Appendix S4, table S4). Among the biggest regions, Africa (AF1) has upper and lower quantiles with significantly negative slopes, but the lower quantile slope is lower, implying that negative fluctuations and variance are increasing (Figure 4). Eurasia mainland (EUAS1) has only the upper quantile significant with a positive slope, suggesting an increase in the variance. North America mainland (NA1) exhibits a significant lower quantile with positive slope, implying decreasing variance. Finally, for mainland Australia, all quantiles are significant, and the slope of the lower quantiles is greater than the upper ones, showing that variance is decreasing (Appendix S4, figure S10).
These results are summarized in Table 1.
The conditions that indicate that a region is near a critical fragmentation threshold are that patch size distributions follow a power law; temporal ∆RS max fluctuations follow a power law; variance of ∆RS max is increasing in time; and skewness is negative. All these conditions were true only for Africa mainland (AF1) and South America tropical & subtropical (SAST1).  Table 1: Regions and indicators of closeness to a critical fragmentation threshold. Where, RS max is the largest patch divided by the total forest area, ∆RS max are the fluctuations of RS max around the mean, skewness was calculated for RS max and the increase or decrease in the variance was estimated using quantile regressions, NS means the results were non-significant. The conditions that determine the closeness to a threshold are: power law distributions in patch sizes and ∆RS max ; increasing variance of ∆RS max and negative skewness.

Discussion
We found that the forest patch distribution of most regions of the world followed power laws spanning seven orders of magnitude. These include tropical rainforest, boreal and temperate forest. Power laws have previously been found for several kinds of vegetation, but never at global scales as in this study. Interestingly, Eurasia does not follow a power law, but it is a geographically extended region, consisting of different domains (as we observed in the largest patch animations, see supplementary data). It is known that the union of two independent power law distributions produces a lognormal distribution (Rooij et al., 2013). Future studies should split this region into two or more new regions, and test if the underlying distributions are power laws.
Several mechanisms have been proposed for the emergence of power laws in forest: the first is related self organized criticality (SOC), when the system is driven by its internal dynamics to a critical state; this has been suggested mainly for fire-driven forests (Zinck & Grimm, 2009;Hantson et al., 2015). Real ecosystems do not seem to meet the requirements of SOC dynamics (Sole et al., 2002;Pueyo et al., 2010;McKenzie & Kennedy, 2012). A second possible mechanism, suggested by Pueyo et al. (2010), is isotropic percolation, when a system is near the critical point power law structures arise. This is equivalent to the random forest model that we explained previously, and requires the tuning of an external environmental condition to carry the system to this point. We did not expect forest growth to be a random process at local scales, but it is possible that combinations of factors cancel out to produce seemingly random forest dynamics at large scales. If this is the case the power law exponent should be theoretically near α = 2.055; this is close but outside the confidence interval we observed (1.898 -1.920). The third mechanism suggested as the cause of pervasive power laws in patch size distribution is facilitation (Manor & Shnerb, 2008;Irvine et al., 2016): a patch surrounded by forest will have a smaller probability of been deforested or degraded than an isolated patch. We hypothesize that models that include facilitation could explain the patterns observed here. The model of Scanlon et al. (2007) showed an α = 1.34 which is also different from our results. Another model but with three states (tree/non-tree/degraded), including local facilitation and grazing, was also used to obtain power laws patch distributions without external tuning, and exhibited deviations from power laws at high grazing pressures (Kéfi et al., 2007). The values of the power law exponent α obtained for this model are dependent on the intensity of facilitation, when facilitation is more intense the exponent is higher, but the maximal values they obtained are still lower than the ones we observed. The interesting point is that the value of the exponent is dependent on the parameters, and thus the observed α might be obtained with some parameter combination.
It has been suggested that a combination of spatial and temporal indicators could more reliably detect critical transitions (Kéfi et al., 2014). In this study, we combined five criteria to detect the closeness to a fragmentation threshold. Two of them were spatial: the forest patch size distribution, and the proportion of the largest patch relative to total forest area (RS max ). The other three were the distribution of temporal fluctuations in the largest patch size, the trend in the variance, and the skewness of the fluctuations. Each one of these is not a strong individual predictor, but their combination gives us an increased degree of confidence about the system being close to a critical transition.
We found that only the tropical forest of Africa and South America met all five criteria, and thus seem to be near a critical fragmentation threshold. This confirms previous studies that point to these two tropical areas as the most affected by deforestation (Hansen et al., 2013). Africa seems to be more affected, because the proportion of the largest patch relative to total forest area (RS max ) is near 30%, which could indicate that the transition is already started. Moreover, this region was estimated to be potentially bistable, with the possibility to completely transform into a savanna (Staver et al., 2011). The main driver of deforestation in this area was smallholder farming.
The region of South America tropical forest has a RS max of more than 60%, suggesting that the fragmentation transition is approaching but not yet started. In this region, a striking decline in deforestation rates has been observed in Brasil, which hosts most of the biggest forest patch, but deforestation rates have continued to increase in surrounding countries (Argentina, Paraguay, Bolivia, Colombia, Peru) thus the region is still at a high risk.
The monitoring of biggest patches is also important because they contain most of the intact forest landscapes defined by Potapov et al. (2008), thus a relatively simple way to evaluate the risk in these areas is to use RS max index. The analysis of RS max reveals that the island of Philippines (SEAS2) seems to be an example of a critical transition from an unconnected to a connected state, the early warning signals can be qualitatively observed: a big fluctuation in a negative direction precedes the transition and then RS max stabilizes over 60% ( Figure S9). In addition, there was a total loss of forest cover of 1.9% from year 2000 to 2012 (Hansen et al., 2013) and deforestation rates were not substantially reduced in 1990-2014; this could be the results of an active intervention of the government promoting conservation and rehabilitation of protected areas, ban of logging old-growth forest, reforestation of barren areas, community-based forestry activities, and sustainable forest management in the country's production forest (Lasco et al., 2008). This confirms that the early warning indicators proposed here work in the correct direction.
The region of Southeast Asia was also one of the most deforested places in the world, but was not detected as a region near a fragmentation threshold. This is probably due to the forest conservation and restoration programs implemented by the Chinese government, which bans logging in natural forests and monitor illegal harvesting (Viña et al., 2016). The MODIS dataset does not detect if native forest is replaced by agroindustrial tree plantations like oil palms, that are among the main drivers of deforestation in this area (Malhi et al., 2014). To improve the estimation of forest patches, data sets as the MODIS cropland probability and others about land use, protected areas, forest type, should be incorporated Sexton et al., 2015).
Deforestation and fragmentation are closely related. At low levels of habitat reduction species population will decline proportionally; this can happen even when the habitat fragments retain connectivity. As habitat reduction continues, the critical threshold is approached and connectivity will have large fluctuations (Brook et al., 2013). This could trigger several negative synergistic effects: populations fluctuations and the possibility of extinctions will rise, increasing patch isolation and decreasing connectivity (Brook et al., 2013). This positive feedback mechanism will be enhanced when the fragmentation threshold is reached, resulting in the loss of most habitat specialist species at a landscape scale (Pardini et al., 2010). Some authors argue that since species have heterogeneous responses to habitat loss and fragmentation, and biotic dispersal is limited, the importance of thresholds is restricted to local scales or even that its existence is questionable (Brook et al., 2013). Fragmentation is by definition a local process that at some point produces emergent phenomena over the entire landscape, even if the area considered is infinite (Oborny et al., 2005). In addition, after a region's fragmentation threshold connectivity decreases, there is still a large and internally well connected patch that can maintain sensitive species (Martensen et al., 2012). What is the time needed for these large patches to become fragmented, and pose a real danger of extinction to a myriad of sensitive species? If a forest is already in a fragmented state, a second critical transition from forest to non-forest could happen, this was called the desertification transition (Corrado et al., 2014). Considering the actual trends of habitat loss, and studying the dynamics of non-forest patches-instead of the forest patches as we did here-the risk of this kind of transition could be estimated. The simple models proposed previously could also be used to estimate if these thresholds are likely to be continuous and reversible or discontinuous and often irreversible (Weissmann & Shnerb, 2016), and the degree of protection (e.g. using the set-asides strategy Banks-Leite et al. (2014)) than would be necessary to stop this trend.
The effectiveness of landscape management is related to the degree of fragmentation, and the criteria to direct reforestation efforts could be focused on regions near a transition (Oborny et al., 2007). Regions that are in an unconnected state require large efforts to recover a connected state, but regions that are near a transition could be easily pushed to a connected state; feedbacks due to facilitation mechanisms might help to maintain this state. Crossing the fragmentation critical point in forests could have negative effects on biodiversity and ecosystem services (Haddad et al., 2015), but it could also produce feedback loops at different levels of the biological hierarchy. This means that a critical transition produced at a continental scale could have effects at the level of communities, food webs, populations, phenotypes and genotypes (Barnosky et al., 2012). All these effects interact with climate change, thus there is a potential production of cascading effects that could lead to a global collapse. Therefore, even if critical thresholds are reached only in some forest regions at a continental scale, a cascading effect with global consequences could still be produced, and may contribute to reach a planetary tipping point (Reyer et al., 2015). The risk of such event will be higher if the dynamics of separate continental regions are coupled (Lenton & Williams, 2013). Using the time series obtained in this work the coupling of the continental could be further investigated. It has been proposed that to assess the probability of a global scale shift, different small scale ecosystems should be studied in parallel (Barnosky et al., 2012). As forest comprises a major proportion of such ecosystems, we think that the transition of forests could be used as a proxy for all the underling changes and as a successful predictor of a planetary tipping 19 point.

Supporting information
Appendix Table S1: Proportion of Power law models not rejected by the goodness of fit test at p ≤ 0.05 level. Table S2: Generalized least squares fit by maximizing the restricted log-likelihood.              Figure S9: Largest patch proportion relative to total forest area RS max , for regions with total forest area less than 10 7 km 2 . Figure S10: Fluctuations of largest patch for regions with total forest area less than 10 7 km 2 . The patch sizes are relativized to the total forest area for that year.

Data Accessibility
Csv text file with model fits for patch size distribution, and model selection for all the regions; Gif Animations of a forest model percolation; Gif animations of largest patches; patch size files for all years and regions used here; and all the R and Matlab scripts are available at figshare http://dx.doi.org/10.6084/m9.figshare. 4263905.