Sacred sites as refuges for native plant diversity in an urbanized tropical landscape

Land-use change, especially the conversion of natural ecosystems into agricultural landscapes, has driven widespread biodiversity loss and declines in ecosystem functioning globally¹. However, many agricultural landscapes retain substantial biodiversity and ecosystem functions, often supported by remnant natural habitats and diverse land-management practices². Urbanization leads to further loss, fragmentation, isolation and degradation of these habitats, which generally results in simplified woody plant communities, with natural regeneration—the spontaneous recruitment of plants without active human intervention—becoming dominated by a subset of disturbance-adapted species^3,4. An increasing proportion of the regional species pool failing to disperse to, or establish in, remaining green spaces reduces long-term ecosystem resilience with consequences for ecosystem service provision^5,6. Understanding the complex, multiscale effects of urbanization on biodiversity has become a central theme in ecological research since the early 2000s, with substantial implications for conservation management and sustainable urban planning^7,8,9. Yet these effects remain poorly understood in developing tropical regions^10,11,12, where urban growth is predicted to accelerate most rapidly¹³, and which might be most affected by such developments^1,11.

Regeneration in plant communities is shaped by a hierarchical set of filters: vegetation cover and landscape configuration influence seed sources and dispersal; local environmental conditions affect establishment and survival and biotic interactions such as herbivory and pollination further modulate plant populations. Urbanization alters each of these filters. For example, the form and arrangement of buildings transform the landscape matrix in which the remaining vegetation patches are embedded, which may create or intensify dispersal and pollination barriers for local species pools^5,14,15. Urbanization also creates new social, economic and cultural contexts that change how green spaces are perceived, used and managed, with direct consequences for local plant and animal communities and environmental conditions. These changes usually favor colonization and persistence of a limited group of species with traits and life histories adapted to these novel urban conditions, resulting in homogenized local plant communities across the landscape^5,16.

The role of urban morphology in influencing such biodiversity patterns is increasingly discussed^5,14,17, including its effects on ecosystem functions and services^18,19. However, how specific building forms and configurations affect biodiversity variation, particularly within rather than across cities, remains understudied—especially in tropical countries^5,7. Tropical cities differ from most temperate ones not only in biogeography but also in the timing and nature of their—often postcolonial—development. Whereas urbanization in temperate regions typically occurred earlier, over longer periods and under more coordinated planning, many tropical cities have undergone rapid, often informal growth, resulting in heterogeneous land use, dense building structures and pronounced peri-urbanization¹¹. Whereas greening efforts have increased vegetation cover in many temperate cities since the early 2000s²⁰, tropical urban landscapes often remain in transitional phases marked by continued habitat loss and fragmentation, with little consideration for biodiversity aspects¹¹. In complex Afrotropical urban contexts, plant community assembly may further be shaped by wealth- and education-related species preferences^12,21,22 or customary and culturally driven land-management practices, for instance, in sacred forest groves²³ or church forests^24,25. The interplay of delayed and accelerated urbanization, contingent development histories and distinct socio-economic and cultural valuations of nature may lead to greater stochasticity and less predictable patterns of plant diversity in such landscapes, particularly for sensitive native species.

This study examines urbanization-related diversity patterns of naturally regenerating native and exotic woody plant species in green spaces across Antananarivo, the capital of Madagascar (Fig. 1). Antananarivo has followed—and continues to follow—decentralized urbanization trajectories rooted in 12 major ‘sacred’ hill settlements (12 collines sacrées de l’Imerina), distributed across the surrounding highlands and once forming the core of a pre-colonial kingdom²⁶. To this day, these hill sites—marking birthplaces, former residences or burial grounds of Malagasy royalty—often remain culturally and spiritually important places of worship (Doany)^{26,27,28,29,30}. This fragmented and peripheral early development history contrasts with remnants of planned, state-supervised city planning, densification and industrial expansion in Antananarivo’s core and surrounding areas during the French colonial era (1896–1960)³¹ and more recent uncontrolled urban sprawl driven by rapid population growth and sustained rural-to-urban migration, reaching wide into former rice plains, wetlands and agricultural lands²⁶. Combined, these processes have created a highly heterogeneous landscape with strong pressures on remnant vegetation and uncertain consequences for biodiversity and ecosystem functions (Methods provide more details on study area).

We studied both local and landscape-scale variation in natural regeneration to assess resilience and capacity of woody plant species and communities to maintain and recover biodiversity in heavily disturbed urban landscapes without human assistance. For our study, we define natural regeneration as unassisted establishment of woody plants through natural seed dispersal or vegetative resprouting and we exclude planted and larger established individuals from analyses (Methods). Here we combined a large-scale field inventory of natural regeneration of tree and shrub species, including over 25,000 individuals of 206 species (39% native), in 329 vegetation patches across the broader Antananarivo area with a comprehensive set of fine-scale urban spatial metrics, socio-economic, historical, cultural and environmental predictors—with the goal of identifying (1) variation in diversity and abundance of naturally regenerating native and exotic species along gradients of these predictors and (2) potential drivers of the observed patterns. Exploratory studies such as this are critical for developing strategies and policies for the conservation and restoration of native plant communities and their ecological roles in complex tropical urban landscapes. This represents one of the first comprehensive analyses of spontaneous urban vegetation in Madagascar and one of few in the context of a lower-income tropical urban context³² (Supplementary Section 0).

We recorded 25,261 naturally regenerating woody plant individuals in 329 vegetation patches (Fig. 1), of which 72% were exotic and 27% native (1% unidentified). These individuals represented 206 species across 55 families, including 81 native (39%), 107 exotic (52%) and 18 undetermined (9%). Across the landscape, native species were less abundant, averaging 84 individuals per species compared to 170 for exotics (Fig. 2), with many—43% of native species and 44% of exotic species—being recorded fewer than ten times across all patches and a substantial proportion (14% of native species, 19% of exotic) only once. The most abundant native species was Psiadia altissima (5% of all individuals in the study area, 20% of natives), whereas Fraxinus excelsior (ash) was most common among exotics (10% of all, 14% of exotics). The most common family among native species was Asteraceae (43% of all native individuals)—Madagascar holds the highest number of native island Asteraceae species worldwide³³—which is known for reproductive flexibility, wind dispersal and persistent seed banks^33,34. Exotic individuals were most frequently from Oleaceae (for example, ash) and Myrtaceae (for example, Eucalyptus robusta), families widely planted for ornamental purposes, fuel wood and other productive uses and among the world’s most common invasive woody taxa³⁵.

Across our 329 sampled patches (mean: 42.5 individuals per transect), species richness ranged from 0 to 14 species per patch for natives (mean: 1.9 ± 0.2 s.d.) and from 0 to 12 for exotics (mean: 3.7 ± 0.2). Native species were absent in nearly a quarter of patches (n = 76), whereas exotics were missing in only 3% (n = 9). Seven patches (2%) contained no recorded natural regeneration. Psiadia altissima (n = 149) and Psidium guajava (n = 200) were the species with the highest occurrence across all patches. Notably, despite the overall higher exotic species richness, the most species-rich patches were dominated by native species (Fig. 2 and Supplementary Section 2).

Species richness patterns of natural regeneration in relation to the studied variables differed between native and exotic species at both local (Fig. 3) and landscape scales (Fig. 4). We first explored whether variation in local (patch-level) diversity was associated with the 12 pre-selected predictor variables, which capture dimensions of urban form and configuration (UFC), socio-economic and cultural factors (SEC) and environmental heterogeneity (EL) (Supplementary Table 1a). To do so, we fitted a negative-binomial generalized linear model (NB-GLM) including all predictor variables, ranked the models using the corrected Akaike Information Criterion (AICc) and computed model-averaged coefficients across the top models (ΔAICc ≤ 2) (Methods). The fit of the averaged model was better for native diversity (McFadden’s R² = 0.15) than for exotics (R² = 0.03), indicating stronger associations between predictor variables and native diversity, with exotic distributions being more stochastic or unrelated to the variables included in our models. Native species richness was predicted to be highest in non-managed vegetation (mean of 1.9 ± 0.3 species per patch), which is about 40–90% higher than in other vegetation types. Native diversity was also higher in sacred (1.6 ± 0.6) than non-sacred (1.0 ± 0.3) sites. Contrastingly, native species richness was predicted to decrease with increasing population density by ~40% per s.d. (incidence rate ratio (IRR) = 0.6 ± 0.1) (Fig. 3). For exotic richness, mean values were generally higher across vegetation types, but associations with numeric predictors were weaker (Extended Data Fig. 1 and Supplementary Section 3).

To examine landscape-level diversity patterns, we selected four variables from each of the three thematic groups (UFC, SEC, EL). Continuous variables were categorized, and species richness was estimated for each category using coverage-based rarefaction and extrapolation at 95% sample coverage (iNEXT function in R; Methods). Estimated landscape-level native species richness was lowest across patches located in areas with greater building density and complexity (UFC variables), higher socio-economic index scores and population density and more recent development. Native species diversity differed most among the Socio-economic index category, with the estimated number of species per subcategory decreasing from 34.9 ± 2.2 to 15.1 ± 1.8 (low to high socio-economic development, respectively). For all variables, landscape-scale exotic species richness was higher than for natives in most categories, with no consistent directional trends (Fig. 4 and Supplementary Section 4).

In both local and landscape-level analyses, we found higher native species richness of natural regeneration in non-managed and sacred patches compared to more intensively managed and non-sacred ones. After accounting for spatial autocorrelation (Methods), native species richness was, on average, predicted to be about twice as high in sacred (n = 32) than in non-sacred (n = 297) patches (Fig. 5), with only a few (n = 8) having lower native diversity than the average per patch across the study area (1.9 ± 0.2). Moreover, six of the seven patches with highest native diversity (> 8 species) were sacred; the other being within 100 m of a sacred patch. Among the top ten patches, seven were sacred. Most of the high-diversity sacred patches were located in a few city-surrounding hill sites (Extended Data Fig. 2). For example, on the royal Ambohimanga hill—former residence of the last kings of Madagascar²⁷—native richness ranged from 7 to 12 species (average across landscape: ~1.9). The highest overall native (14 species) and mean total species richness (20.7 species) were recorded in Antsahadinta, a more isolated, non-urbanized sacred hill containing ancient royal tombs and houses of the noble class. Ambohijoky, another former royal hill site²⁹, held two of the top seven and three of the top ten patches with the most native species (Fig. 5b). Contrastingly, exotic species richness differed little between sacred and non-sacred sites (Fig. 5c), and none of the top-ranking patches in exotic diversity were sacred, with only one in the top ten (Fig. 5a,b). Patches with high native diversity outside sacred areas reflected other historical contingencies. They were typically adjacent to sacred sites, situated within (semi-)abandoned agroforestry systems, or in deep, unmanaged hadivory formations—trenches historically used to protect dwellings in this region³⁶ (Supplementary Table 5b).

Sacred patches were generally larger and contained unmanaged vegetation (69%, n = 22), typically located in pre-colonial neighborhoods (56%, n = 18) and embedded in hilly landscapes with low building density and high surrounding tree cover (Extended Data Fig. 3). The sacred versus non-sacred contrast (cultural–spiritual significance) was thus partly confounded with other variables. ‘Drop-one’ model comparisons (Methods) showed that removing vegetation type increased the sacred coefficient from 0.49 (± 0.14 SE) to 0.72 (± 0.15), indicating both variables jointly explain part of the variation in native diversity but also that the sacred versus non-sacred contrast independently explains part of this variation (likelihood ratio (LR) = 5.73, ΔAIC = 3.7, P = 0.02) (Supplementary Table 5c,d).

It is widely recognized that urbanization drives biodiversity loss, particularly of native species. Yet more recently, initiatives such as the 3–30–300 rule for urban tree-canopy design³⁷, green streets³⁸ and biophilic cities³⁹ reflect a growing narrative that cities may also mitigate such loss, while providing important ecosystem services^8,9,40. However, such nature-based solution narratives are largely shaped by more affluent temperate cities. Our results indicate that biodiversity patterns in heterogeneous, rapidly transforming tropical urban landscapes may be shaped by multiple, interacting socio-ecological factors that can generate variable barriers to the natural regeneration of diverse plant communities, highlighting the need for context-specific research approaches and conservation planning. To our knowledge, the only study examining similar questions in a lower-income tropical urban context is by Bigirimana et al. (2011) in Bujumbura, Burundi, who assessed spontaneous vegetation in 437 patches in relation to anthropogenic (urban land cover, trampling intensity) and environmental (soil type, shade) factors—finding exotic species dominant and native diversity declining toward more urbanized areas³². Our study expands on this by assessing natural regeneration of woody plants across multiple spatial scales and a wide range of urbanization-related variables, representing a first in Afrotropical urban ecological research. Three principal findings emerged: (1) species richness declines along urbanization gradients for natives but not for exotics, (2) native diversity is typically higher in little-managed ancient, sacred hill sites and (3) urbanization effects on natural regeneration are complex and interconnected.

Species richness declines along urbanization gradients for natives but not for exotics

Diversity of natural regeneration was associated with UFC at the landscape scale but not the local scale, with native species richness across sites declining from less to more complex and dense building arrangements. We did not observe elevated diversity in moderately urbanized zones, such as peri-urban areas, a pattern often reported for affluent temperate cities and attributed to the temporary coexistence of disturbance-adapted native and exotic species⁴¹. We found similar patterns when examining the effect of SEC factors, namely decreasing landscape-level native species diversity towards regions with higher socio-economic development, characterized by modern buildings and infrastructures expected in more urban areas. Similar to other studies from Africa^21,22, we observed no ‘luxury effect’, a common phenomenon in more developed cities where wealthier neighborhoods with better education support more diverse plant communities^42,43. Likewise, we did not find elevated native species richness in private gardens, as reported in other African cities^12,21,22,44. In those cases, elevated biodiversity often appears to be driven by actively planted ornamental species^21,42 and is therefore not indicative of plant community resilience through natural regeneration—which was the focus of our study. Exotic species richness remained comparatively stable along UFC and SEC gradients, suggesting greater tolerance to urbanization-associated environmental variability.

Native diversity is typically higher in little-managed ancient, sacred hill sites

Prior research has highlighted how land-use dynamics drive patterns of natural regeneration in tropical agricultural landscapes^2,45,46. The recent and rapid urbanization in such regions implies that legacy effects from past land use and landscape change may still influence present-day species assemblages¹¹. Biodiversity studies in such contexts should account for both current urban conditions and historical land-use trajectories. In Antananarivo, this includes both the continued, multi‑nuclear influence of pre‑colonial village settlements²⁶ and the spatial legacy of colonial‑era city planning³¹. Whereas a global meta-analysis found a positive relationship between city age and plant biodiversity⁷, within-city analyses have reported inconsistent associations between neighborhood age and biodiversity^42,43,47—highlighting context-dependence and idiosyncrasy of such relationships. In our study, native diversity tended to be higher in pre-colonial neighborhoods compared to newer districts, which could partly be due to more unmanaged and sacred sites in these areas (Extended Data Fig. 3). Such findings align with evidence from other African regions, including church forests in Ethiopia²⁴, urban churchyards in South Africa²⁵ and sacred groves in Morocco²³. Whereas higher native biodiversity in minimally managed vegetation fragments may be directionally expected, it is notable that the most species-rich patches across our study area were almost exclusively located within or adjacent to the sacred royal hill sites surrounding the city. This pronounced concentration of native species may partly reflect lower urban pressures in these peripheral sites but could also be driven or reinforced by other factors, such as customary conservation practices associated with Malagasy ancestral rules and socio-cultural taboos (fady) against cutting specific endemic plants in such sites, including forests around tombs^30,48. Because these alternative explanations operate over historical and socio-cultural dimensions not captured by our spatial and compositional data, incorporating additional information on land-use history, education and socio-cultural context may help disentangle how surrounding populations influence nature in such historical sacred sites.

Urbanization effects on natural regeneration are complex and interconnected

Our results suggest complex associations between variables of urban form and configuration and the diversity of natural regeneration at local and landscape scales. For instance, the four UFC variables used in this study may capture different spatial effects influencing natural regeneration (for example, anthropogenic dispersal barriers) or, perhaps more likely, reflect broader socio-economic and cultural factors that influence how people perceive, manage and use local urban green spaces. Also, many variables and underlying drivers may covary and interact, making isolating single potential drivers of diversity inherently challenging. Hence our explanatory variables are likely to collectively capture the overarching ‘the more urban, the fewer native species’ trend observed in other Afrotropical studies^32,44,49. Therefore, while urbanization might seem like an intuitive term, it encompasses spatially and conceptually multifaceted processes that are difficult to capture in complex tropical landscapes—especially when relying on observational data and correlation analysis. Our study points to the resilience of native plant communities in rapidly urbanizing landscapes and emphasizes that, to further advance our understanding of how different aspects of urbanization alter the drivers and barriers to natural regeneration, we need long-term, landscape-scale data—including on spatial-temporal (seasonal) variation in seed dispersal and plant recruitment, growth and survival—combined with additional information on (changing) landscape context (for example, degree and quality of connectivity), local environmental conditions (for example, soil properties) and socio-economic context (for example, land-use dynamics), replicated across diverse city contexts.

On the basis of our findings, we recommend exploring three avenues for conservation and restoration of diverse and resilient plant communities in the broader study area: (1) foster conditions favorable for natural regeneration, (2) focus on species of high conservation value and (3) sacred sites as socio-cultural conservation strategy.

Foster conditions favorable for natural regeneration

Our extensive survey of woody vegetation across the broader Antananarivo area identified a strongly depleted regional native species pool. The small size, fragmentation and degradation of remaining forested patches constitute strong barriers to the natural regeneration of most species and ongoing, unregulated urbanization further erodes the resilience of native tree communities, leading to biotic homogenization. Most studies on the diversity and functioning of urban trees and forests concentrate on mature, established individuals, whose composition largely mirrors legacy conditions and historical planting decisions and thus offers limited insight into the community’s resilience—that is, capacity to sustain diversity and ecosystem functions through natural regeneration⁴⁵. Research from tropical secondary forests^2,45 and urban forests^6,16,50 suggests that promoting natural regeneration can be a viable complement to active revegetation. This is especially relevant in lower-income tropical countries, where municipalities often lack resources for direct ecological interventions¹¹. In Antananarivo, green space restoration on the city outskirts that reflects historical landscape characteristics, such as open or mosaic forest ecosystems⁵¹, could help establish peripheral seed sources that support long-term biodiversity conservation through natural dispersal. Here Psiadia altissima—the native species we recorded most frequently in our study area (Supplementary Table 2c)—could be a promising, disturbance-adapted nurse species for supporting natural regeneration. Research from Madagascar found elevated mycorrhizal soil infectivity in the rhizosphere of Psiadia altissima, which may facilitate the establishment of later-arriving species⁵². Alternatively, planting trials in the country’s Central Highlands showed high survival of the endemic Dodonaea madagascariensis (recorded in two of our patches) in plots close to existing forest fragments, suggesting that it may be a promising candidate for reforestation plantings intended to facilitate natural succession in grassland openings⁵³.

Focus on species of high conservation value

The Central Highlands surrounding Antananarivo are considered a center of paleo-endemism⁵⁴, yet afforestation efforts in the region often focus on exotic species such as Eucalyptus spp., Pinus spp. and Acacia spp.⁵¹. Whereas these deliberately introduced species have important socio-economic value for local populations⁵⁵, they are also known to increase crop damage and management costs⁴⁶, disrupt soil and water cycles⁵⁶ and reduce biodiversity and regeneration capacity of native species^46,51. Razafimanantsoa et al. (2025) advocated shifting the focus from exotic mono-plantations to restoring forest patches that protect remaining fire-adapted native Central Highland species, such as Uapaca cf. bojeri⁵¹—which we recorded at multiple sites (Supplementary Fig. 2a). Alternatively, Cadia pubescens—a shrub native to the Central Highlands and recorded in one of our plots—could be prioritized for restoration based on its endangered status on the International Union for Conservation of Nature Red List, while including the ‘sacred’ Ficus lutea, widely planted by former regional kings for ritual and symbolic purposes and recorded in many of our plots, could provide important socio-cultural and spiritual functions²⁸. However, due to limited paleo-ecological and conservation-focused research in the region^51,55, evidence-based recommendations and guidance on species selection for long-term ecological restoration initiatives in and around Antananarivo remain challenging, stressing the need for targeted field studies and experiments. The dataset compiled in this study could support such efforts by providing spatially explicit information on where, and under what landscape conditions, particular species occur and regenerate naturally. In urban contexts specifically, institutions such as the Silo National des Graines Forestières (SNGF; https://www.sngf-silo.com/), under the Ministry of Environment and Sustainable Development (MEDD)—which produces forest seedlings for private individuals and has contributed substantially to introducing native species into the urban landscape—could play a key role in sourcing different species to support restoration initiatives.

Sacred sites as socio-cultural conservation strategy

Sacred sites on the surrounding hilltops, though often small, isolated and under increasing pressure, offer key conservation opportunities. This aligns with recent global meta-analyses indicating that conserving sacred forests is particularly effective for maintaining elevated plant diversity compared to other taxonomic groups⁵⁷. Because protection levels vary across sites, effective ecological conservation is inherently tied to maintaining and strengthening their historical and spiritual significance, ideally combined with active enrichment planting. Collaboration with local authorities, including traditional site guardians (known as Mpimasy), and with national—such as the National Office of Arts and Culture dedicated to the conservation and promotion of historical sites—and international institutions, is an important step in this direction. For instance, the royal Ambohimanga site contains one of the few remaining fragments of original highland forest, supporting many culturally important native species, such as ‘royal’ Ficus and Dracaena sp., and is protected through a combination of community-based practices, national heritage laws and designation by the United Nations Educational, Scientific and Cultural Organization²⁷. In Antsahadinta—another of the 12 sacred hill sites considered national heritage and the location with the most native and overall species in our study—forests have been protected from clearing and slash-and-burn agriculture through longstanding animistic-spiritual rules (fady)³⁰. This suggests that integrating socio-culturally driven land stewardship into national conservation planning could be key to sustaining native species diversity in and around Antananarivo.

This study provides a comprehensive, field-based assessment of woody plant diversity in natural regeneration across a rapidly transforming urban landscape in a lower-income tropical country. We found that urbanization further reduces native species richness across an already heavily degraded landscape, a pattern broadly reflected across gradients of different urban morphological and socio-economic variables. In contrast, exotic diversity showed no distinct pattern at either local or landscape scales. Although thematically diverse predictor variables captured direct and indirect dimensions of urbanization across different spatial scales, our results also underline the inherent challenges of studying plant communities in heavily human-modified urban landscapes. Biodiversity management strategies in such complex contexts must therefore account for both ecological constraints and socio-cultural realities. In Antananarivo, particularly within the sacred hill sites, ecological and historical conservation are closely interlinked. Integrating traditional and formal ecological perspectives and knowledge systems will be key to guiding long-term biodiversity conservation and restoration efforts in such landscapes.

Study area

The study was conducted in Antananarivo (18.9185° S, 47.5211° E), Madagascar’s capital and its peri-urban and rural surroundings (Fig. 1). Antananarivo’s population is about 3 million inhabitants⁵⁸ and increasing at a rate of about 5% annually⁵⁵. The climate is temperate subtropical with cool, dry winters, with mean temperatures ranging from ~11 to 27 °C (ref. ⁵⁹) and annual rainfall of ~1,000–1,500 mm, concentrated in the warm season (November–April)⁵⁵. Paleo-ecological and historical evidence indicate that the Central Highlands encompassing Antananarivo originally comprised a heterogeneous mosaic of forest, shrubby woodland and grasslands and that human settlement (approximately since 2 thousand years ago) and associated burning and grazing gradually transformed the region toward the largely treeless and heavily degraded savannah-like shrubby grasslands that dominate today^51,59,60. Vegetation around Antananarivo mostly consists of rice fields, short-rotation eucalyptus and pine plantations for charcoal production and timber and other agricultural land uses^{55,59,61,62,63}. Within urban areas, larger vegetation fragments consist of productive or ornamental vegetation in gardens and parks in different degrees of maintenance.

Vegetation data collection and processing

For sample site selection, we used a recent high-resolution land-cover map of Antananarivo^63,64 to classify all (n = 856) Fokontany (smallest administrative unit in Madagascar, similar to neighborhoods/villages) within a 40 × 40 km study area into four classes along an urban–rural gradient (QGIS: Jenks natural breaks for % impervious surface). We then randomly (using random draws) selected Fokontany from each class and inventoried 1–7 (mean 2.0 ± 0.1) vegetation patches per Fokontany, starting, where possible, with the largest one. If there were no (accessible) vegetation patches, we selected another Fokontany. Patches spanning multiple Fokontany were assigned to the one covering the largest area. Sample locations were selected using recent high-resolution satellite imagery (Google Earth Pro, v7.3), with vegetation patches initially defined as continuous areas of visible tree or shrub cover. During field visits, we recorded GPS coordinates and detailed notes on landscape context. On the basis of these records, we delineated patches (hand-drawing tool in Google Earth Engine) as (2) structurally functionally homogeneous units, in terms of vegetation type, abiotic environment and land management and (2) not separated by major roads or large buildings. In the special case of sacred patches, classification and delineation were based on a combination of (1) information gained from local landowners and/or land inhabitants, (2) Fokontany chiefs, (3) traditional site guardians, (4) expertise of the Malagasy co-authors and (5) external sources^{26,27,28,29,30,48,65,66} (Supplementary Table 5a).

All patches were ≥ 100 m² except one (65 m²) due to a mapping error. Accessibility was sometimes limited by land ownership (for example, private gardens and parks), physical barriers (for example, steep hill slopes) or safety considerations, especially in the most urbanized areas. Finding intact vegetation in the city center was difficult, whereas sampling in more rural areas was constrained by travel time and limited road access. Also, because Jenks natural breaks classification assigned more Fokontany to the most rural class 4, we adjusted sampling in this class to retain the study focus on more urbanized areas. Due to these logistical and methodological constraints, the proportional number of sampled Fokontany per class varied (Class 1: 42 of 105 Fokontany available within class; Class 2: 112/133; Class 3: 86/181; Class 4: 89/437). This initial ‘coarse’ urban–rural gradient served solely to randomize site selection and was not used in subsequent analyses.

Sampling was conducted between June 2022 and January 2024 in 40 m² transects (n = 604), each split into ten 2 × 2 m squares placed in the central/densest part of the patch, oriented east–west whenever possible. We placed at least one transect per patch, with more added (n = 1–6) based on patch size and (sometimes) time constraints. In each transect, we recorded all woody (trees and shrubs) plants > 5 cm tall and measured diameter at breast height (DBH) for individuals > 1.3 m tall. Some sampling locations contained several vegetation patches that were numbered accordingly (Fig. 5b). We collected photos, seeds, fruits and flowers of all species and deposited them in the University of Antananarivo herbarium database. We used the herbarium and the authors’ taxonomic expertise (Author contributions) in addition to global repositories (POWO: https://powo.science.kew.org/; GBIF: https://www.gbif.org/; Tropicos: https://www.tropicos.org) to identify and validate 99% of all recorded plant individuals (Supplementary Section 2). We classified species into native (= naturally occurring in Madagascar) and exotic using the mentioned repositories and reference literature^67,68,69. The data included 29,341 woody plant individuals in 329 vegetation patches across the study area (Fig. 1). To focus on natural regeneration, we restricted analyses to the 25,261 individuals with DBH ≤ 5 cm, an operational size threshold used in other studies as an approximate proxy for relatively recent or ongoing regeneration (that is, separating smaller regenerating individuals from established adult vegetation)^70,71. We estimate that >95% of these individuals represent unassisted recruits, established either from naturally dispersed seeds (germination) or via vegetative resprouting (from roots, stumps or cut stems). We distinguished spontaneous recruits from planted individuals using field observations and, where possible, statements from landowners or caretakers about land-use history (for example, Supplementary Table 5b). In unmanaged or abandoned sites, natural regeneration was usually evident from irregular growth form and heterogeneous size and age patterns, whereas planted vegetation growth was more uniform and usually limited to a few species. Additionally, during vegetation data sampling, we collected data on canopy opening and vegetation type.

Predictor variables

To examine if the diversity of natural regeneration varies predictably along direct and indirect urbanization gradients, we compiled an initial set of 30 potential predictor variables from a variety of sources (Supplementary Table 1a for more details), grouped in three thematic categories that either directly (urban form and configuration, or UFC, 12 variables) or indirectly (socio-economic and culture or SEC, 13 variables) reflect urbanization characteristics⁵ or environmental heterogeneity across the landscape (environment and landscape, EL; five variables).

All predictors of the UFC category were computed from an open buildings map^72,73. We included the interquartile range (IQR) and median for each of four building-level metrics (area, shape, compactness and elongation) as separate variables⁷³. Analysis of variance adjusted for multiple comparisons indicated IQR as a better predictor of species richness for compactness, elongation and shape index metrics (F ≈ 12.4–24.7; p_adj and p_Holm ≤ 0.003), whereas the median was slightly stronger for building area. For consistency across metrics, we retained IQR for all four variables (Supplementary Table 1f). We conducted principal component analysis (PCA) on standardized (z-scored) variables using Pearson correlation matrices to reduce the predictor sets within the UFC category to smaller, comparable sets of principal components (PCs). We visually inspected distributions to check for extreme outliers. PC1, PC2 and PC3 collectively explained >91% of total variation in the UFC category, with the four building-level variables loading highest on PC1, the two indicators of spatial compactness and aggregation—the area-weighted compactness index (AWCI)⁷⁴ and aggregation index (AI)⁷⁵)—on PC2 and building density and the nearest-neighbor index (NNI)⁷³ on PC3 (Supplementary Fig. 1c and Supplementary Table 1d,e). Accordingly, we retained PC1–PC3 for further analysis.

The SEC category consists of publicly available household-level data from Madagascar’s most recent (2018) national population and housing census conducted by the national statistics office (INSTAT⁷⁶) and data on population density (subnational population statistics⁷⁷), neighborhood age (variable historic records) and cultural–spiritual significance (Supplementary Table 1a). The latter three were retained individually as we considered them to be thematically distinct and/or were obtained from very different information sources. Combined use of household-level information on access to different services has been suggested in other African studies²² as a determinant of socio-economic status, particularly in cases where data on household income or education are lacking or are too coarse⁵. We assembled the ten housing-level census variables (items 13–22 in Supplementary Table 1a), standardized (z-scored) them and applied PCA to reduce dimensionality. PC1 explained 71.5% of the total variation and reflected a coherent socio-economic development gradient, with further PC axes contributing little additional thematic structure (Supplementary Fig. 1b and Supplementary Table 1a–c). Consequently, we retained only the PC1 score as a Socio-economic index, replacing the original variables. The five variables of the EL category (patch size, canopy openness, surrounding tree cover, topography and vegetation type) were retained individually as well, for the same reasons as outlined above.

Combined, this yielded 12 predictors: three UFC, four SEC and five EL variables. Because population density (Fokontany‑level) and socio‑economic index (commune‑level, each comprising on average 11 Fokontany) were only available at coarser scales, we assigned their respective values to all patches within those administrative units. To evaluate potential spatial clustering bias and pseudo-replication from such downscaling, we fitted negative-binomial generalized linear mixed models (GLMMs; lme4⁷⁸ in RStudio) for species richness with intercept-only random spatial effects (Fokontany or commune) and the two predictors as respective fixed effects. To reduce potential effects related to data skewness, we square-root transformed population density before model fitting. We then extracted Pearson residuals, linked them to patch centroids and computed Moran’s I (k = 8 nearest neighbors), which revealed no significant spatial autocorrelation for either (p > 0.05). Marginal R² (R²m) was modest and conditional R² (R²c) was higher (socio-economic index R²m ≈ 0.08–0.09, R²c ≈ 0.21–0.25; population density R²m ≈ 0.10–0.11, R²c ≈ 0.38–0.42), indicating any spatial structure was absorbed by random intercepts rather than remaining in residuals. These results justified aggregating these higher-order predictors to the 329 patches without further spatial correction.

All variables were calculated or compiled on either patch level, 100 m patch-surrounding buffer zones or neighborhood scale. The 100 m buffer was selected after initially testing 25, 50 and 100 m radii for two key base variables—percent impervious (urban) and forested (tree) land cover—which showed very high Pearson correlation across scales (r > 0.85) (Supplementary Fig. 1a). The 100 m radius provides a reasonable balance between capturing local landscape effects and avoiding oversensitivity to fine-scale noise and is commonly used in such type of research^79,80,81,82. Variance inflation tests indicated no problematic multicollinearity between the 12 final predictor variables (VIF for all variables < 1.78).

Modeling local-scale species richness variation

To test whether patch-level species richness of natural regeneration varied with the 12 predictors, we fitted an additive NB-GLM (using MASS⁸³, v7.3–65), separately for natives and exotics, on all 12 predictors (the global model) and all subsets of variables. Numeric predictors were z-standardized before modeling, and categorical predictors (vegetation type, neighborhood age, cultural–spiritual significance) were treated as factors with their respective levels (for example, vegetation type:ornamental; cultural–spiritual significance:sacred). We included a logarithmic offset for the number of transects per patch to account for unequal sampling intensity. We ranked candidate models by AICc (AICcmodavg⁸⁴, v2.3–4), performed model dredging using MuMIn (v1.48.11)⁸⁵, averaged the best-supported set (ΔAICc ≤ 2) and reported McFadden’s pseudo-R² to assess model fit. To quantify each predictor’s association with native and exotic diversity, we extracted model-averaged coefficients, summed Akaike weights across all top models and model inclusion counts. Top models showed no problematic multicollinearity (max VIF ≈ 3), and adding quadratic or cubic terms to z-scored predictors did mostly increase their AICc, so we retained linear specifications. To measure differences among multi-level (> 2) variables (vegetation type, neighborhood age), we obtained estimated marginal means from the main effect global NB-GLM, using emmeans (v2.0.0)⁸⁶ and then used Dunnett-adjusted contrasts⁸⁷ to compare each subcategory with a designated reference level through IRRs (Fig. 3). Potential zero-inflation (too many patches without either native or exotic species) was assessed by comparing the same negative binomial (NB) used for local analysis with zero-inflated NB-GLMs using AICc and parametric bootstrap LR tests (200 simulations) and simulation checks of zeros (n = 1,000)⁸⁸ (pscl⁸⁹, v1.5.9). NB was always favored over zero-inflated NB-GLMs (lower AICc for both native and exotic responses), and bootstrap LR tests did not support adding a zero-inflation component (for both p > 0.4). Observed zero counts (n = 112) fell within the NB-simulated 95% range for natives (103–132) and slightly below for exotics (12 versus 17–36). These diagnostics indicated that zero inflation was not a material issue for our data. To verify if including interaction effects would have substantially changed our inferences in terms of the directions and relative strengths of main effects, we ran the same procedure as described above but with a global model that included two-way interaction effects—which revealed nearly identical qualitative patterns (Supplementary Section 3 for data and specifications for the local species richness analysis).

To compare patch‑level diversity between sacred (n = 32) and non‑sacred sites (n = 297), we used two‑sided Brunner–Munzel rank‑order tests, which are robust to unequal group variances⁹⁰, and reported Cliff’s delta with 95% confidence intervals (CIs) (Fig. 5). Native (adj. p < 0.001) and exotic species richness (adj. p = 0.04) were tested in parallel using Holm adjustment to control for multiple comparisons. To test whether spatial clustering could explain the sacred–non-sacred contrast, we calculated Pearson residuals from non-spatial NB models (k = 8 nearest neighbors) based on patch centroids, whereas accounting for unequal patch sizes by including log(patch_size) as a covariate. Spatial autocorrelation was moderate but significant (p < 0.001) for both natives (Moran’s I = 0.20, z = 8.14) and exotics (I = 0.08, z = 3.14). Using mgcv (v1.9–3)⁹¹, we then fitted a generalized additive model (GAM) with a bivariate smooth of projected x and y coordinates to account for potential nonlinear spatial structure. For natives, the sacred effect remained significant (p < 0.001) and not attributable to spatial clustering; for exotics, it remained non-significant (p = 0.13). To assess if the cultural–spiritual significance predictor was confounded with other variables, we ran ‘drop-one’ analyses where each predictor was omitted in turn from the global NB-GLM used for dredging, then compared these reduced and full models using ΔAICc and LR statistics and examined changes in the sacred versus non-sacred coefficient⁹² (Supplementary Table 5c,d).

Estimating landscape-level species richness variation

We assessed variation in landscape-level species diversity of natural regeneration in relation to four variables from each of the three thematic categories (UFC, SEC, EL), all of which were included in at least 50% of the top models of the local analyses for either native or exotic species. For each variable, patches were assigned to 2–4 categories, and species richness was estimated separately for each of these categories. Categorical variables were kept as is, whereas numeric variables were categorized based on tertile thresholds (that is, 1/3, 2/3 and 3/3 quantiles; ~110 patches/tertile) (Fig. 4, Extended Data Fig. 4 and Supplementary Table 4a). PC2 and PC3 (UFC) were among those variables, but because these did not reflect an obvious thematic gradient, we decided to select the variables with strongest loading on these axes instead, namely building density, NNI and AI. AWCI was removed due to high correlation (Pearson r = 0.96) with the thematically more independent AI. PBefore species richness estimation, for AI and NNI only, we removed locations without or too little buildings required to calculate these indices. PC1 of the UFC category, reflecting closely related building-level metrics, was retained as morphological complexity. Pairwise Pearson correlations among retained predictors were modest, with some expected covariation among thematically close predictors (for example, population versus building density, r = 0.80). However, as multicollinearity was not problematic, with a maximum adjusted generalized variance inflation factor of 2.2 (ref. ⁹³), all predictors were retained for individual rarefaction/extrapolation analysis (Supplementary Section 4 provides data and specifications related to the landscape-level diversity analysis).

The number of patches and sampled individuals differed across categories. Typically, such differences are addressed by comparing species richness at equal sample sizes through rarefaction⁹⁴. However, consider two contrasting urbanization categories: the first with high local diversity, many rare species and pronounced differences in species composition among patches;and the second dominated by a limited number of widespread species. With equal sampling effort, these two categories would be compared at substantially different levels of sample coverage—that is, at different fractions of their total and actual landscape diversity^94,95. An increasingly common solution to reflect uneven abundance is to standardize species comparisons by equalizing sample coverage—the proportion of total expected individuals or occurrences captured—using a combination of rarefaction and extrapolation⁹⁶. In landscapes containing many rare species, this approach may yield more robust estimates of true diversity^94,95. The R package iNEXT⁹⁷ allows to generate coverage-based species accumulation curves, extrapolated up to twice the reference sample size and to estimate diversity at specified coverage levels. Here we estimated species richness using the estimateD function (iNEXT, v3.0.1, q = 0, Chao1-type estimator, abundance data)⁹⁸ at 95% sample completeness, which was the highest coverage achievable for each group within reasonable 95% bootstrap (n = 50) CIs (Supplementary Fig. 4b,c for a visual explanation of the method)⁹⁴. All categories reached 95% coverage through interpolation (rarefaction) only. Exploratory runs for Shannon (q = 1) and Simpson diversity (q = 2) showed consistent patterns across all diversity orders.

Data processing and analysis

Initial preparation of predictor variables was conducted in QGIS (v3.40.3) and Google Earth Engine⁹⁹, using the most recently available Google satellite imagery (2024–2025) and different open-access layers (administrative, roads, buildings, topography). In Python, we used GeoPandas (v0.8.1)¹⁰⁰, Shapely (v2.1.2)¹⁰¹, Rasterio (v1.4)¹⁰² and SciPy (v1.0, built on NumPy)¹⁰³ to calculate morphometrics—including density, AI, AWCI, NNI and building-level area, shape index, compactness, elongation and orientation—and Pearson correlation matrices for buffer size comparisons. Vector metrics used building polygons, while AI was computed from rasterized building footprints at 1-m resolution⁷⁵. All predictors were reprojected to a common coordinate reference system. In RStudio (v2024.12.1), tidyverse (v2.0.0, data wrangling)¹⁰⁴, vegan (v2.7–5, community analyses)¹⁰⁵, prcomp (v4.4.0, PCA; base R)¹⁰⁶, sp and spdep (v2.3.4, spatial autocorrelation)¹⁰⁷ and car (v3.1–5, VIF)¹⁰⁸ were used for a variety of data preparation and analysis steps. Visualizations were created with ggplot2¹⁰⁴ and corrplot (v0.95)¹⁰⁹. ChatGPT (OpenAI) was used to support limited parts of the coding, and all outputs were subsequently reviewed, tested and verified by the authors.

Inclusion and ethics statement

We acknowledge the importance of equitable representation of diverse authorship, perspectives and regions in ecological research. All components of this study were developed in close collaboration with local researchers, with Malagasy co-authors actively involved in every stage of the research process, including study design, fieldwork, data interpretation and manuscript preparation. We view this contribution as a foundation for future work in the field and hope it can support and inform studies led by scientists based in the regions represented.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.