Biodiversity conservation in the sacred groves of north-west Ethiopia: diversity and community structure of woody species

of of of in fenced church forests. Seedling abundance depended on the local presence of large, mature conspeci ﬁ c trees and on the geographic distance to potential source populations of seeds, strongly suggest that metapopulation dynamics likely are important. We conclude that church forest conservation and minimizing further degradation of the landscape matrix are needed to help sustaining the ecological and socio- economic potential of this unique network of remnant forests. 0 ) and Pielou evenness (J) were calculated using the R ‘ vegan ’ package, once based on tree basal area (most strongly re ﬂ ecting a tree species ’ current contribution to forest structure and ecosystem services) and once using tree abundance (indicative of likely prospective importance). To elucidate which patch-and landscape level variables most strongly shape biodiversity patterns, we used linear models with a Gaussian error structure, specifying the different biodiversity indices as dependent variable. We applied a model selection procedure based on Akaike ’ s Information Criterium AIC (Burnham and Anderson 2004) and calculated AICc values for all possible models, using the R ‘ MuMIn ’ package. Models were ranked based on their AICc values, and the relative importance of variables was assessed by summing the AICc weights of all models in which the variable under consideration was included. Variables characterized by AICc weights > 0.5 and with model-averaged estimates exceeding their standard errors, were considered important predictors (Anderson 2008). Model residuals were normally distributed (i.e. all Shapiro-Wilk W > 0.90). Note that all (continuous) explanatory variables were ﬁ rst scaled by subtracting the mean and dividing by the standard deviation.


Introduction
Habitat loss and fragmentation have been recognized as major conservation threats to terrestrial biodiversity (Rands et al., 2010;Haddad et al., 2015;Martín-Queller et al., 2017). Loss and fragmentation typically occur simultaneously, and while there is debate about the relative importance of each (Fahrig 2017;Fletcher et al., 2018;Fahrig et al., 2019), the division of natural areas into smaller and more isolated fragments generally causes a reduction in the amount of habitat available to species and increasing edge effects (Ewers and Didham 2006). Forests have been hit particularly hard by fragmentation, as about 70% of the world's remaining forests are currently within 1 km of the forest's edge (Haddad et al., 2015). Habitat loss and fragmentation affect both forest structure and species composition (Echeverría et al., 2007). Apart from the fact that reduced population sizes in smaller patches negatively influence species persistence (Shaffer 1981;Wilsey et al., 2005), decreases in forest patch size additionally lead to reductions in tree basal area (which reflects the stand density and the strength of the microclimate inside forest patches), and a lower recruitment of seedlings (further increasing the risk that species will go extinct in the future, i.e. extinction debts) (Benitez-Malvido 1998).
Biodiversity of remaining forest patches can be determined by both patch-and landscape-level characteristics. At the patch level, patch size, quality and shape strongly influence biodiversity dynamics (Fahrig 2003). Edge effects such as changes in microclimate are stronger in smaller and irregular-shaped patches, resulting in habitat degradation and increased loss of core habitat (Laurance et al., 1997;Fahrig 2017). Particularly in tropical forests, such changes are likely to negatively affect large and old trees in favor of pioneer trees (Haddad et al., 2015), ultimately impacting forest regeneration due to a lack of seeds, which are predominantly produced by large, mature trees (Jorge and Garcia 1997). Besides patch-level characteristics, the ecological integrity of forest fragments also depends on land use and land cover changes in the matrix (i.e. at the landscape level). Functional connectivity between patches depends on the geographical distance between patches and the structure of the landscape matrix in which patches are embedded ). An inhospitably matrix can hamper dispersal and lead to reduced connectivity among remaining fragments, negatively affecting metacommunity structure as diminished recolonization ultimately increases species extinction rates. Such landscape-level processes may be especially important for tropical forests, as more than 90% of tropical woody plant species depend on zoochory (Bregman et al., 2016) while population persistence of seed-dispersing animals codepends on the quality of the surrounding matrix as a (secondary) habitat resource (Antongiovanni and Metzger 2005;Kennedy et al., 2011).
In many regions of the world, habitats are so severely degraded that only mosaics of smaller forest patches remain, while many competing land use demands exist (Nagendra et al., 2004), hindering full restoration of their former continuous state (Gibson et al., 2011;Watson et al., 2018). It is therefore crucial to identify factors that promote biodiversity and provision of ecosystem services in networks of forest remnants, in order to halt further losses and inform recovery and conservation strategies. Such information is especially needed to avert irreversible degradation of the Ethiopian Afromontane forests, which harbor high biodiversity and are of considerable economic and socio-ecological importance (Wassie et al., 2010). Agricultural expansion, urbanization, population pressure, demand for fuelwood and construction material, and other development activities have led to drastic deforestation (Nyssen et al., 2004;Hailu et al., 2015;Mesfin et al., 2016) across northern Ethiopia, and the majority of native forests currently only manifest as small patches that surround churches and monasteries of the Ethiopian Orthodox Tewahido Church (EOTC, Fig. 1). In the Lake Tana Basin, for example, woody vegetation cover sharply declined during the mid-20th century (Frankl et al., 2019), and currently covers only~2.3% of the area (i.e. including 0.76% natural forests and 1.56% woodlands) (Song et al., 2018). Particularly in northwestern Ethiopia, churches and monasteries are easily identified even from a great distance by their preserved trees (Wassie et al. 2005(Wassie et al. , 2010. Despite their generally small size, typically only a couple of hectares, and their isolation within a mainly agricultural or agropastoral matrix, church forests deliver important ecosystem services and are crucial for biodiversity conservation. Many rare and indigenous trees and shrubs, which went extinct in most parts of Ethiopia, are still found in the compounds of these churches (Wassie et al., 2005;Dudley et al., 2009). For example, Wassie et al. (2010) reported a total of 160 indigenous woody species in 28 studied church forests across the South Gondar region, which corresponds to more than 60% of the 263 tree species that occur in tropical northeast Africa (Friis et al., 2010). Church forests also serve as refuges for wildlife (Wassie et al., 2005;Scull et al., 2017), for example, in the Lake Tana Basin alone 437 species of bird occur, which is more than half of the total bird fauna of the country (Shimelis and Abebe, 2017). Similarly, the area harbours a considerable diversity of amphibians, lizards, snakes, freshwater turtles, reptiles and mammals too (Shimelis and Mengistu 2017). This suggests that church forests may be as, or even more important, for biodiversity conservation as the few contiguous natural forests still left in Ethiopia (Wassie et al., 2005;Aerts et al., 2016).
Yet, the long-term conservation value of this network of church forests is threatened by deteriorating habitat quality in individual patches combined with changes in the surrounding landscape matrix (Wassie et al. 2009a(Wassie et al. , 2010Aerts et al., 2016). For example, Cardelu's et al. (2019) found that human disturbance was high in more than half of the 44 studied church forests in northwest Ethiopia, resulting in reduced tree species richness, biomass and density. Inside many church forests, degradation is ongoing, as logging and tree dieback of species such as Juniperus procera and Olea europaea spp. cuspidata lead to community shifts, from dry Afromontane forest towards shrubland . Additionally, plantations of commercially important exotic species inside the church forests, such as Cupressus lusitanica, Citrus spp., Eucalyptus spp., Grevillea robusta, Jacaranda mimosifolia, Melia azedarach, and Pinus spp. strongly affect church forest species composition (Aerts et al., 2016;Cardelu's et al., 2019). Landscape-level changes in agricultural practices result in (even) fewer trees and bushland in the landscape matrix, intensifying the impacts of edge effects and isolation. For example, declining woody biomass in the matrix could result in increasing wood gathering pressure inside the church forest area  while grazing reduces forest regeneration and negatively affects the seedling bank due to increased soil compaction in and around the church forests (Wassie et al., 2009b;Cardelús et al., 2013). Beyond land use practices, the presence of roads and distance to population centers can also contribute to forest degradation (Cardelús et al., 2019). Several local communities have started erecting stone walls around the perimeter of the forests in an attempt to protect the forest interior ). Yet, ongoing habitat degradation and isolation necessitates improved protecting of church forests in order to achieve effective biodiversity conservation (Aerts et al. 2006(Aerts et al. , 2016Wassie 2007;Wassie et al., 2010;Scull et al., 2017).
Formulating conservation strategies is however hindered by the fact that across much of northwest Ethiopia, little information is available on the structure and species composition of church forests, nor is it known which patch or landscape characteristics most strongly influence forest biodiversity. Future planning of effective conservation programs additionally requires information on forest community structure (e.g. seedling recruitment), yet such knowledge is virtually absent. Therefore, here, we aim to fill this knowledge gap by performing a spatially explicit assessment of church forests and their surroundings in the Lake Tana Basin in northwest Ethiopia. Because of their typically long lifespan and regeneration times, forest tree species often show a time-lagged response to habitat changes (Metzger et al., 2009;Muscolo et al., 2014). Species may thus survive for a long time in small and degraded forest patches, before eventually becoming extinct. Local (i.e. patchlevel) extinctions may be counterbalanced by (re)colonization when species are still present in other nearby forest patches (i.e. in the regional species pool). To account for such dynamics, aside from surveying tree species composition and vegetation structure of church forests, our assessment also considers the prevalence of regenerating species (i.e. species present as both mature trees and saplings or seedlings), potential extinction debts (i.e. tree species only present in mature stages) and recruitment credits (i.e. tree species only present as seedling or sapling and thus potentially complementing the future tree layer composition (Thijs et al., 2014). Note that introduced, non-native species were excluded from the analyses because their inclusion would interfere with the unique evolutionary history and phylogeny of the region's Afromontane forests (sensu Pauchard et al., 2018). Instead, factors promoting the establishment of non-native species were assessed separately. Thus, concluding, this study is therefore done with the following core questions in mind: 1) what is the current state of church forest biodiversity and which factors shape the observed tree communities, 2) which factors determine (a) tree community structure, (b) the presence or absence of tree seedlings, and (c) the invasion of forests by non-native species?

Description of the study area and site selection
The Lake Tana Basin has an area of 15,089 km 2 with an elevation ranging from 1.782 to 4.109 m above sea level in the northwest highlands of Ethiopia (Lemma et al., 2018). This study was conducted in the lower elevations of this Basin southeast of Lake Tana (1794e2204 m a. s. l.). Cultivated land is the predominant land-use, interspersed with natural forests, bushland, plantation forests, urban areas, villages, waterbodies, wetlands, woodlands and grasslands (Song et al., 2018). The study area experiences a strong seasonal rainfall regime (about 70e90% of the total rainfall occurs during JuneeSeptember), but a very even monthly temperature regime (Peel et al., 2007). The annual average rainfall varies between 1250 and 1500 mm, and the mean annual temperature varies between 15.3 C and 19.6 C (Lemma et al., 2018).
Using Google Earth satellite maps, a total of 24 church forests located between 37 27 0 E À 37 55 0 E and 11 39 0 N e 11 56 0 N were selected from three districts namely Bahir Dar Zuriya, Dera and Fogera (Table 1, Fig. 2). Of the 24 church forests, three, nine and twelve church forests belonged to Bahir Dar Zuriya, Dera, and Fogera districts, respectively. The surface area of the church forests ranged from 2 ha to 13 ha and the churches in the forests were established between the years 340e2010. Several church forests are neighbored by larger, natural, conserved forests, which may act as seed sources and may provide important habitat for seed dispersing species. However, generally, the matrix surrounding the church forests represents a variegated landscape of smallholder agricultural land, where 'Tef' (Eragrostis teff), wheat (Triticum aestivum), maize (Zea mays), sorghum (Sorghum bicolor), barley (Hordeum vulgare) and finger millet (Eleusine coracana) are the main crops (information derived from 21 key informant interviews. Livestock husbandry is the second means of livelihood for the local community.

Data collection: church forest vegetation composition and structure
In each church forest, four sampling plots were systematically selected in the four cardinal directions (north, east, west, south), but on different distances along the axis from church to edge. The starting direction was selected randomly and continued in a clockwise direction (Sup Fig S1). A plot size of 20 m Â 20 m was employed to identify mature (i.e. diameter at breast height (DBH) ! 5 cm) and sapling (DBH<5 cm) trees. In each of the four sampling plots, all tree species were identified, counted (i.e. the number of individuals present) and measured (i.e. height and DBH of each tree). Height and DBH were measured using a 'Nikon Forestry Pro' rangefinder and diameter tape, respectively. When the height of the tree was >1.6 m, DBH was measured at 1.3 m above ground level, but when the height was <1.6 m and >1 m, the diameter was measured at 10 cm above the ground. For buttress root trees, the DBH was measured above the buttress and for multiple stems, all stems were counted and measured (Wassie et al., 2010). Trees having <1 m height were considered as a seedling, and to measure the presence and abundance of tree seedlings, a 5 m Â 5 m plot was established in each corner of the sampling plots (thus totaling 4 such 5 m Â 5 m plots per church forest). Identification was done to species level using identification keys (Bekele 2007). Species that were difficult to identify were identified in the national herbarium at Addis Ababa University. Nomenclature was based on the published guidelines of the Flora of Ethiopia and Eritrea, Vol. 1e8 (1989e2009) (Edwards et al. 1995(Edwards et al. , 1997Edwards 1997;Hedberg et al., 2003). The recently updated nomenclature for Acacia spp. was also used (Kyalangalilwa et al., 2013). Lastly, the percent of canopy cover was measured using a convex mirror densiometer at five points within each of the main plotplots (Sup Fig S1), and the average recording was used for each plot. 2.3. Data collection: church forest structural characteristics A combination of field observations, interviews with church forest leaders and Google Earth Pro images were used to describe the church forest patch characteristics. For each church forest, data was collected for the presence or absence of stone wall (this permanent fence is expected to protect seedlings against cattle), year of establishment (possibly related to the age of the agricultural use of the matrix and the time since fragmentation of the forest), altitude (related to the species composition) and forest fragment area (related to the risk of extinction). To account for possible effect of soil conditions, based on the main report of the 'Abay river basin integrated development master plan', church forests were assigned one of the following soil types: haplic or chromic luvisol, eutric vertisol, leptosol or fluvisol (BCEOM 1999). Together with canopy cover (see above), these metrics are henceforward referred to as 'patch level' variables.

Data collection: landscape-level characteristics
To assess the status of woody vegetation in the matrix, we first created a 1 km buffer area around each church forest using google earth pro. Inside these buffer areas, crown size of each tree or each group of trees (i.e. when single crowns were not distinguishable) was mapped. Resulting tree cover maps were saved as.KMZ files and then converted to a raster.tiff using the QGIS 3.4 software and further processed using Landscape Ecology Statistics (LecoS), an open-source plugin integrated into the QGIS processing framework for landscape ecology analysis (Jung 2016). Considering each cluster of trees as a 'matrix tree group' then allows to derive a number of variables describing how remnant trees are distributed in the landscape. For each church forest, the following landscape-level metrics were obtained: total number of tree groups, smallest tree group area, median tree group area, edge density (i.e. total length of tree group edge divided by total landscape area, larger values indicate more fragmented tree cover) and the fractal dimension index (which measures the irregularity of matrix tree group configurations, higher values indicate potentially stronger edge effects). In addition to these tree cover related metrics, the following variables were considered as well: distance to Bahir Dar (which reflects both the influence of closeness to the region's largest urban center, as well as, inversely, the historical colonization and anthropogenic modification of the area starting from the mountainous areas in the north-west, Fig. 2), distance to the nearest river and to the nearest permanent road. Additionally, presence or absence of natural reserve forest within 1 km radius was noted as well (as these larger forest may serve as source areas for seed dispersal). Anthropogenic land-use in the 1 km buffer was classified as either 'agriculture' or 'agropastoral'. Lastly, human population density was approximated by asking church forest leaders about the number of Christian households in the vicinity of the church forests.

Data collection: biodiversity metrics
Apart from total (a) species richness, the Shannon diversity (H') and Pielou evenness (J) indices were computed for each church forest, using both the tree abundance and basal area data sets. Basal area (m 2 /ha) was obtained by first calculating the basal area per tree (using as P x (DBH) 2 /40.000), which was then summed across all trees measured and converted to m 2 /ha using a conversion factor of 10.000/1.600 (as four 20 Â 20m sampling plots were surveyed).

Data collection: community structure
Three different community structure metrics were obtained, based on presence/absence data for each of the three tree life stages (mature, sapling and seedling), following the definitions provided by Thijs et al. (2014). First, 'regenerating species' are those tree species that are present in both the mature tree layer and the seedling or sapling layer. Second, the 'potential extinction debt' was quantified as the number of tree species that only occurred as mature individuals in a fragment (i.e. no seedlings or saplings present). Third, a 'potential recruitment credit' was calculated as the number of tree species represented by seedlings or saplings only, without representatives in the mature tree layer.

Statistical analyses of floristic and structural composition of church forests
Biodiversity indices: Biodiversity indices (Shannon diversity (H 0 ) and Pielou evenness (J) were calculated using the R 'vegan' package, once based on tree basal area (most strongly reflecting a tree species' current contribution to forest structure and ecosystem services) and once using tree abundance (indicative of likely prospective importance). To elucidate which patchand landscape level variables most strongly shape biodiversity patterns, we used linear models with a Gaussian error structure, specifying the different biodiversity indices as dependent variable. We applied a model selection procedure based on Akaike's Information Criterium AIC (Burnham and Anderson 2004) and calculated AICc values for all possible models, using the R 'MuMIn' package. Models were ranked based on their AICc values, and the relative importance of variables was assessed by summing the AICc weights of all models in which the variable under consideration was included. Variables characterized by AICc weights >0.5 and with model-averaged estimates exceeding their standard errors, were considered important predictors (Anderson 2008). Model residuals were normally distributed (i.e. all Shapiro-Wilk W > 0.90). Note that all (continuous) explanatory variables were first scaled by subtracting the mean and dividing by the standard deviation.
Tree species distributions. First, a Bray-Curtis dissimilarity matrix summarizing the raw tree community data was created, once using tree basal area and once based on tree abundance data. Then, the geographical (Euclidean) distances between church forests was calculated. Next, distances matrices were calculated summarizing dissimilarities in patch-and landscape level church forest characteristics by relying on Gower's distance to account for the categorical nature of several variables. Lastly, we applied Mantel and partial Mantel tests to examine the correlations between tree community composition and the geographic and environmental distance (patch and landscape) predictor variables, using 9999 data permutations. Second, a Bray-Curtis dissimilarity matrix based non-metric multidimensional scaling (NMDS) was used to identify specific patch and landscape-level variables influencing church forest community compositions. Ordinations were performed using the 'met-aMDS' function of the 'vegan' package, with a k ¼ 2 dimensionality. Patch and landscape environmental variables were fitted to ordinations using 'envfit' and their significance ascertained (p < 0.05) using 9999 permutations.
Forest types and indicator species: Using either tree abundance or tree basal area, the 24 church forests were grouped using full hierarchical clustering based on a Bray-Curtis dissimilarity matrix and flexible beta linkage. The best beta linkage, and therefore the best clustering of the dissimilarity matrix, was found using a mantel test between the dissimilarity matrix and the cophenetic correlation coefficient (McCune and Grace 2002). The highest possible number of indicator species in conjunction with the lowest average p-values were used to decide the optimal number of clusters. For each clustering level, the significant indicator species were identified at p-value < 0.05. To test for multivariate differences in the community composition of the identified clusters, a multi-response permutation procedure test (MRPP) with natural group weighting factor n i /S n i was used (n i is the number of church forests in each group). An indicator species analysis (sensu Dufrêne and Legendre, 1997) was then performed on the resulting cluster (function indval of R package 'labdsv').
Community structure: Patch and landscape-level factors controlling the number of regenerating tree species, the extinction debt and the recruitment credit of church forests were identified using linear models and the AIC-based multi-model selection mentioned above. The number of tree species per community structure metrics was divided by church forest total species richness and analyzed using a Gaussian error structure.
Seedling status: To uncover processes fostering seedling establishment, seedling abundance was related to the set of patch and landscape-level characteristics described above, but with the inclusion of two extra variables. For each species present as seedling, the basal area of that tree species, and the geographical distance to the nearest church forest where mature trees of that species were present were obtained. Mixed linear models were run under the AIC-based multi-model selection framework, specifying tree species and church forest as random effects and using a Poisson error distribution.
Non-native species: Exotic species richness, basal area and abundance were similarly assessed as native tree species, using linear models and the same set of patch and landscape-level variables as above.

Results
Biodiversity indices: In total, 56 indigenous tree species representing 33 families were recorded across the 24 church forests, and average species richness equaled 14.33 ± 3.43 (range: 8e23). Considering mature and sapling trees together, the church forests have a Shannon diversity of 2.51 ± 0.23 (based on tree species abundance) and 1.72 ± 0.27 (based on tree basal area) and a Pielou diversity index of 0.95 ± 0.02 and 0.66 ± 0.10, respectively. Seedling species richness was on average 7.39 ± 2.61 (range: 3e11).
Tree species distributions. For mature tree species and saplings, abundance based analyses showed that differences in tree communities between church forests was correlated with the geographic distance between them (r ¼ 0.12, p ¼ 0.059), while partial mantel tests revealed that aside from the effect of geographical distance, church forests with similar patch (but not landscape) level characteristics had a similar vegetation composition (r ¼ 0.17, p ¼ 0.011; r ¼ 0.076, p ¼ 0.20 respectively). A NMDS ('stress' level 0.18; linear fit R 2 ¼ 0.82, non-metric R 2 : 0.966) found that similarity among church forest communities was mainly driven by their year of establishment (r ¼ 0.26, p ¼ 0.041) and altitude (r ¼ 0.53, p ¼ 0.0003), and by the presence of a stone-wall too, when allowing a 90% confidence interval (r ¼ 0.11, P ¼ 0.080). In contrast to the partial mantel test above, the NMDS did retrieve landscape-level variables indicating that vegetation similarity was also influence by distance to Bahir Dar (r ¼ 0.42, p ¼ 0.0038), distance to rivers (r ¼ 0.27, p ¼ 0.039), the presence of natural forests (r ¼ 0.16, p ¼ 0.032) and land use (r ¼ 0.29, p ¼ 0.0026). Basal area based mantel tests only found similarities between vegetation composition and patchlevel characteristics (r ¼ 0.16, p ¼ 0.013), but not with geographical distance or landscape-level variables (all p ! 0.22, Supplementary Material SM1). A NMDS ('stress' level 0.23; linear fit R 2 ¼ 0.691, non-metric R 2 : 0.945) found that altitude (r ¼ 0.37, P ¼ 0.0073), church forest surface area (r ¼ 0.25, P ¼ 0.048), and under a 90% confidence interval also the presence of a stone wall (r ¼ 0.10, p ¼ 0.088) correlated with church forest vegetation similarity. Additionally, despite the fact that the partial mantel test did not uncover relationships with landscape-level variables, similarity among church forests did correlate with the number of household (r ¼ 0.28, P ¼ 0.031).
Forest types and indicator species: For mature tree species and saplings, the optimal number of clusters was three when based on abundance, but ten for basal area based analyses. Given that ten forest types for only 24 church forests is an unrealistically high number (e.g. see Aerts et al., 2016), only results for the more realistic abundance based analyses are discussed. Significant MRPP tests (p-values ¼ 0.001) confirm that differences exist in the tree species composition between the types of church forests identified by the silhouette-based cluster analysis (Table 2). A test statistics of T ¼ 0.71 indicated a maximum separation between groups and a within group homogeneity in comparison with random expectation of A ¼ 0.10 (sensu McCune and Grace, 2002). Cluster membership and indicator species values can be found in Supplementary Material SM2 and SM3).
Community structure: All 56 indigenous tree species had mature stage individuals, while only 48 and 36 species were also represented by saplings and seedlings stages, respectively. The number of regenerating trees species per church forest varied from 3 to 15 (mean and standard deviation 8.25 ± 3.04, total across all church forests: 51 species). Potential extinction debt was on average 3.04 ± 2.63 (range 0e11). Most tree species had seedlings or saplings in a limited number of church forests only (Fig. 4, Supplementary Material SM8). Five tree species may even face extinction across all studied forests as they were only found as mature trees (i.e. Euphorbia abyssinica, Ficus ingens, Ficus sycomorus, Myrica salicifolia and Senna petersiana, all species with known medical use; Njung et al., 2002;Olayinka et al., 2017;Sobeh et al., 2017;Oghenesuvwe et al., 2018). Another 15 species only had seedlings or sapling in a single church forest (Fig. 4). Recruitment credits averaged 4.25 ± 2.86 species (range 1e11). Higher numbers of regenerating tree species (i.e. tree species that are present both in the mature tree Patch and landscape-level variables that most strongly correlate with church forest woody species diversity. Linear relationships are visualized using partial residual plots. In the combined violin and boxplot, each dots represents a church forest and plot outlines illustrate kernel probability density, i.e. the width of the area represents the proportion of the data located there. layer and the seedling or sapling layer) correlated positively with church forest surface area (estimate ± standard error, AIC c variable weight: 0.46 ± 0.33, 0.78), canopy cover (0.07 ± 0.06, 0.68), and the number of matrix forest patches (0.029 ± 0.020, 0.80), but negatively with distance to Bahir Dar (À0.11 ± 0.06, 0.90, Fig. 5). Recruitment credit (i.e. tree species only present as seedlings or saplings) was positively related to year of establishment (3.53 ± 1.91, 0.90) and was higher in church forests protected by a stone wall (0.33 ± 0.17, 0.89). Extinction debt (i.e. mature trees without regeneration) was not related to any of the variables tested here (Supplementary Material SM5a and SM5b). Seedling status: The Zahar Mikal (ZM, 11 in Fig. 2) church forest had the highest total number of seedlings (~3.000), while no seedlings at all were recorded in Delemo Tekelhayemanot (DT, 13 in Fig. 2) church forest. Overall, twenty tree species had no seedlings in any surveyed church forest while about 50% had less than 10 seedlings. Diospyros abyssinica had the highest number of seedlings (3.899), although 2.773 of those were found in a single church forest only (Zahar Mikal (ZM), 11 in Fig. 2). The second most abundantly regenerating seedling species overall was Mimusops kummel, with a total of 2773 seedlings. Seedling abundance was higher in larger church forests (estimate ± standard error, AIC c variable weight: 0.31 ± 0.25, 0.77) and correlated positively with the basal area of conspecific mature trees (0.73 ± 0.14, 1.00). The smaller the geographic distance to another church forest harboring mature conspecific trees, the more seedlings were present (À0.21 ± 0.20, 0.66), while the presence of natural forest in the landscape also contributed to seedling abundance (0.86 ± 0.55, 0.84; Supplementary Material SM6).
Exotic species status: A total of 14 exotic species were found inside the church forests. Of these, Eucalyptus camaldulensis with a total of 751 individuals was the most abundant species, followed by Cupressus lusitanica, having a total of 53 individuals  4. Histogram depicting the number of church forests harboring seedlings or saplings for each woody species found across the 24 church forests surveyed here. Five tree species already face a regional extinction debt as no juvenile stages are found in any church forest. Many more species are on the brink as saplings or seedlings were found only in a small number of church forests, as e.g. no less than 15 species were found as seedling or sapling in a single church forest only.
( Table 2). Of the 14 exotic species, eight were present in both the mature tree layer and the seedling or sapling layer, two occurred only as mature individuals and another four were found as seedlings or saplings only. The number of exotic tree species per church forest did not correlate with any of the patch or landscape-level variables tested here, nor did exotic tree basal area (Supplementary Material SM7). When expressed as tree abundance, more exotic tree species were found when church forests were closer to a permanent road (estimate ± standard error, AIC c variable weight: À0.33 ± 0.22, 0.76) but lower numbers were found in forests equipped with a stone-wall (À0.71 ± 0.40, 0.85).

Overall status of church forest biodiversity
Knowledge on the floristic composition of a vegetation community is a prerequisite to understanding its overall status, structure and function. The inventories conducted here highlight the importance of churches for forest conservation, as tree species richness was high. A total of 56 indigenous tree species was recorded, which rivals the number of woody species that are reported for the few contiguous natural forests left in the region (e.g. 43 woody species in Menagesha Suba forest and 66 woody species in Hugumburda forest, both national forest priority areas, Demissew, 1988;Aynekulu, 2011;Aerts et al., 2016). The 24 church forests studied here represent about 20% of the ca. 270 tree species that occur in tropical northeast Africa (Friis 1992). The value of networks of church forests is further exemplified by the generally high diversity indices found. Especially when based on tree abundance, the Shannon index, while high values for the Pielou index indicate comparatively little variation in abundance between species, presenting communities in which species are rather equally abundant. Basal area based estimates point in the same direction, yet the generally lower values suggest that tree communities are to some extent dominated by the presence of a few large, old trees that represent the bulk of the total tree basal area. While overall tree diversity is high, community unevenness and strong extinction debts signal a negative biodiversity trend. Higher human population densities in the landscape matrix are associated with reduced biodiversity. Up to five tree species are at risk of going extinct in the study area as no juvenile stages were detected across any of the forest inventoried, while permanent transport infrastructure brings the risk of invasive non-native species.

Patch and landscape level correlates of church forest biodiversity
The current state and vegetation composition of church forests in northern Ethiopia results from the interaction between natural and anthropogenic factors, both recent and historical. Church forests located closer to each other tended to have a more similar composition, possible reflecting the fact that neighbouring forest fragments may have experienced a more similar deforestation history, which started in the north-eastern part of the region and gradually expanded to the south (Pankhurst 1995). This is also reflected by the fact that church forests with a similar year of establishment tend to have more similar vegetation communities, whereby older forest likely are composed of unique assemblages of indigenous woody species that have long disappeared from the surrounding areas, whereas this may not be the case for more recently established church forests. Likewise, land-use patterns affected similarity among church forests. Scull et al. (2017) noted that the conversion of agropastoral landscapes to a more intensive agriculture often is accompanied with a reduction of bushlands, effectively removing the buffer church forests had, ultimately impacting seedling regeneration and growth through increased edge effects, including increased wood extraction. Furthermore, Frankl et al. (2019) noted that vegetation corridors that previously protected valley floors from erosion largely disappeared, and therewith their functioning as landscape connectivity elements. Indeed, forest similarity was associated with the number of households in the surrounding landscape, as species richness was lower and less even in more populated areas. Interestingly, church forests equipped with a stone wall tended to be more similar to each other, which could point to early successes in protecting the forests, although this result could also arise as a consequence of the selection of church forests where walls have recently been erected . Landscape-level effects on church forests also manifested as differences in altitude, distance to rivers and proximity of natural forests. While this study was conducted in the lower elevations southeast of Lake Tana and thus spans a limited altitudinal range only (~1800 to 2200 m), (remnants of) Afromontane forests are strongly structured along altitudinal gradients (Wassie et al., 2010). Such patterns are most likely undergirded by climatic differences in precipitation and temperature, together with varying (micro)climate conditions such as moisture availability and topographic water holding potential. The same is likely to be true for the influence of distance to rivers, indicating a separation between the once dominant dry evergreen Afromontane forests covering the region, characterized by species such as Millettia ferruginea and Teclea nobilis, and remnants of moist, riparian forest found along rivers, represented by species such as Cordia africana, Celtis africana and Grewia ferruginea (Aerts et al., 2006;Yirgu and Delvare, 2019, see below, Table 2). Lastly, church forests closer to protected natural forests shared similarities, suggesting that these larger, natural forest serve as source patches, from which species may emigrate out of and into the surrounding matrix, rescuing species in isolated patches.

Forest types and indicators species
In general, the church forests have important contribution to the regional conservation of original potential vegetation type, more specifically DAF. Church forests can be classified into three variants of the 'Dry evergreen Afromontane forest and grassland complex' (DAF), namely the subtypes 'Undifferentiated Afromontane forest' (DAF/U), 'Dry single-dominant Afromontane forest of the Ethiopian highlands' (DAF/SD) and 'Afromontane woodland, wooded grassland and grassland' (DAF/ WG) (Friis et al., 2010). The first tree community group (DAF/U) consists of two significant indicator species (Millettia ferruginea and Teclea nobilis) and four general indicator species (Pittosporum viridiflorum, Bersama abyssinica, Prunus africana, Ekebergia capensis). The second group (DAF/SD) is characterized by one significant indicator species only (Acokanthera schimperi) plus three general indicator species (Juniperus procera, Olea europaea and Vachellia abyssinica). This forest is also known as the 'Juniperus e Olea forest', dominated by late-successional and higher altitude species (Friis et al., 2010). The last group consists of lower altitude species Croton macrostachyus, Grewia ferruginea, Dichrostachys cinerea and Vachellia lahai, indicative of the DAF/WG forest (Friis et al., 2010), which also includes widely distributed species such as Cordia africana ( Table 2). The occurrence of these forest types across the study is at least partly governed by elevational gradients, with higher areas showing reduced soil moisture, strongly affecting region wide tree distributions (Engelbrecht et al., 2007). Local soil conditions can be important too, as, for example, Vachellia were most prominent in Kirekus (K, 12 in Fig. 2) and Qere Giweregis (QG, 16 in Fig. 2) church forests, most likely because of their vertisol soils, which are preferred by Vachellia species (Friis et al., 2010).

Patch and landscape level correlates of church forest biodiversity
All but five species were present as both mature trees and as saplings or seedlings, but this high (>90%) overall prevalence of regenerating species obscures the fact that many tree species only regenerate in a small number of church forests (Fig. 4), and that extinction debts at the level of individual church forests are substantial. Only four church forests did not have any potential extinction debt and, on average, church forests are set to lose about three species, though some forests face losing up to eleven trees species in the future. No direct correlations were found linked to extinction debts, but several factors influencing tree species regeneration and recruitment were identified, shedding light on the processes influencing church forest community structure (Fig. 5). Regenerating rates provide an opportunity to assess long term sustainability of the forests, and the fact that forests further away from Bahir Dar show fewer regeneration suggests these habitats have been more strongly degraded because of the longer time period of anthropogenic disturbances they experienced. Species regeneration was higher in church forests characterized by a denser canopy cover. Higher canopy cover likely contributes to a microclimate more favorable to seed germination and seedling growth, especially as many church forest trees are shade tolerant species. Closed canopies add more litter to the soil, increasing its fertility (Berhane et al., 2013) and also lower the erosive impacts of raindrops, reducing runoff and resulting in less soil erosion and consequently less soil carbon loss (Mekuria et al., 2011). Regeneration also correlated positively with church forest surface area, which can be linked to the edge effects, as smaller forest have relatively fewer core area and thus fewer opportunities for seedlings to grow and germinate under optimal microclimates (Cardelu's et al., 2019). Alternatively, smaller forests may be more likely to suffer from higher self-fecundation and inbreeding depression, negatively effecting seedling success (Barbeta et al., 2011). The importance of surface area is also reflected by its direct effect on the abundance of seedlings, as larger forests showed higher densities of seedlings. Conspecific basal area was an important predictor of seedling abundance too, pointing to the importance of large, mature trees for securing species persistence in isolated fragments (Lindenmayer and Laurance 2017). Larger trees not only produce seeds locally, but also provide an important habitat for dispersal agents such as birds, who for example forage, rest or seek refuge from predators within their canopies. It should be noted that many of the large trees found in church forests produce fleshy fruits (e.g. Ficus spp., Mimusops kummel), attracting frugivores, and anecdotal evidence suggests that seeds are dispersed through bird faeces, leading to enhanced recruitment of seedlings especially under such large trees (pers. obs. F.M.).
Indeed, apart from patch-level characteristics, seedling abundance was also affected by landscape-level factors such as the distance to church forests harboring conspecific mature trees and the presence of larger, natural protected forests in the surrounding landscape. The importance of geographical distance to potential source populations points to the importance of seed dispersal for dynamics of plant populations and communities. Church forests act as 'islands' within a largely agricultural matrix, and given their overall relatively small area and anthropogenic pressures they face, their degree of functional isolation may be especially important. If seed dispersal between fragments occurs sufficiently frequently to influence local demography while maintaining subpopulation independence, networks of church forest may be considered a metapopulation (Hanski 1998). Species may be able to persist in such habitat networks notwithstanding high rates of local extinction in individual fragments. The fact that all church forest have a recruitment credit (on average about four species), whereby species are only present as seedling or saplings and thus must have been dispersed into the forest strongly suggest that such metapopulation dynamics may indeed be at play here. Seed dispersal likely is dependent on the functional connectivity of the church forest network, which considers not only Euclidean distance between fragments but also how seed dispersing animal vectors respond to various elements of the landscape matrix (Tischendorf and Fahrig 2000). Little is known about the communities of seed dispersing organisms in church forest communities. A study on two common church forest trees (Olea europaea and Schefflera abyssinica) by frugivorous birds however found overall relatively low dispersal distances, suggesting that matrix patches of forest serving as 'stepping stones' may be needed to ensure effective seed dispersal between church forests (Abiyu et al., 2016). The relationship found here between recruitment and the number patches in the surrounding landscape can be interpreted as support for such functional landscape connectivity as an important factor for the persistence of church forest biodiversity.

Invasive species
Dispersal can, at least for some species, however also be human-assisted, and indeed, the abundance of introduced, nonnative species was higher in church forests that are close to permanent road infrastructure. Fourteen different non-native species were identified across 24 church forests, a diversity that is somewhat higher than -but comparable to-the twenty non-native species found in 46 church forest surveyed by Mekonnen (2019) in the same general region (Table 2). Mekonnen (2019) also reported that non-native species were mainly associated with church forests located closer to towns accessible roads create a good link for marketing. Indeed, previous studies also have reported that non-native Eucalyptus plantations are more frequently found in church forests closer to roads (Orlowska and Klepeis 2018). This is likely due to the economic importance of the planted non-native trees, and the presence of roads create easy access for buyers, who customarily come with trucks to collect the trees from the plantations. Only two invasive species faced a potential extinction in the region, causing cause for concern and fitting recent studies reporting that intensification of church forest utilization increasingly results in the dominance of disturbance-tolerant native but also invasive species (Berhane et al., 2015). Interestingly, nonnative species abundance was lower in church forests equipped with a protective stone wall, which fits Mekonnen (2019) observation that church forests without a stone wall along the forest periphery have a more non-native woody species.
Stone walls do not only prevent access to the forest by grazing cattle but also directs the people to use a trail instead of randomly entering and walking through the forest. Woods et al. (2017) report that such reduced anthropogenic pressure is associated with increased native tree regeneration and higher densities of native tree seedling, possible strengthen a 'biotic resistance' (Levine et al., 2004) hindering invasive species establishment success. Indeed, also here, native species recruitment credits were higher when stone walls were present, and it is worth to note that seedlings of Celtis africana, Diospyros abyssinica, Millettia ferruginea, Rothmannia urcelliformis and Teclea nobilis were only found in church forests with a stone wall, indicating the important role stone walls likely play for the persistence of these species.

Conservation implications
Church forests play an important role in conserving many economically and ecologically important tree species, yet due to various anthropogenic pressures, their integrity is threatened. Focusing on the state and drivers of seedling and sapling dynamics as indicators of future tree recruitment, we found that all church forests face potential extinction debts, and that five tree species already are at risk of disappearing from the study region altogether. Importantly, these five species are registered as traditional medicinal plants, highlighting that together with unique biodiversity, church forest degradation risks the loss of important ecosystem services too. To counteract such losses, not only in church forests specifically, but more broadly of regional ecosystems such as the dry Afromontane forests (DAF), conserving indigenous forest is crucial and conservation priority should be given to the church forests, which are critical for biodiversity conservation and provide numerous ecosystem services to the local people. This study allows to suggests several conservation measures that may be considered to help restore church forest integrity. Given the degree of potential extinction debts uncovered here, to enhance forest regeneration potential, conservation strategies such as active planting programs both inside the church and its matrix may be considered. Fencing the forests with a stone wall also likely contributes to forest regeneration and may help conserving a minimum forest size by preventing wood harvesting from the forest edge, as well as discourage the creation of new paths and clearings. Removal of exotic species and preventing any further human-assisted introductions of non-native species should be encouraged. At the landscape level, management should focus on strengthen the functioning of the metapopulation structure of the church forest network. This could be achieved by conserving the large, old trees still present as they function as important seed sources, by minimizing the further degradation of the landscape matrix and by increasing connectivity through rehabilitation gullies into green corridors. In many cases, conservation initiatives usually prioritize biodiversity and ecosystem function, while human use receives minimal attention (Reynolds et al., 2017). Therefore, for sustainable conservation, human use should be fully integrated into conservation planning -indeed, the interdependence between these forests and the church community highlights that human socio-economic development and nature conservation cannot be uncoupled.

Conflicts of interest/Competing interests
The authors are not aware of any conflicts of interest or competing interests.

Availability of data and material
The data used in this study have been uploaded as zip files for review and will be made available publicly on Figshare or Dryad upon acceptance of the ms.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.