Simultaneous effect of habitat remnancy, exotic species, and anthropogenic disturbance on orchid diversity in South Australia

Orchids are potentially useful as ecological indicators because of their sensitivity to habitat fragmentation and anthropogenic disturbance. While many studies explore the effect of single factors on orchid diversity, few investigate how the extent, configuration, and condition of surrounding habitat affect whole orchid communities. Here, we unravel the effect of biological invasions, anthropogenic disturbance (i.e., grazing pressure, ecological condition), and habitat fragmentation on an Australian orchid community. We sampled 39 plots across nine sites in the Mount Lofty Ranges, Australia. We recorded the number of orchid species and number of individuals per species in mid‐winter, early‐spring, and late‐spring to account for the effect of season on species visibility, with 115 surveys in total. We ranked grazing intensity and ecological condition, and estimated cover of exotic species. We analyzed the response of richness and diversity through generalized linear mixed models, and differences in species composition through non‐metric multidimensional scaling. Habitat configuration in the surrounding landscape had different effects at increasing radii, explaining 29%–87% of variance. Patch‐level orchid diversity was positively correlated with habitat edges in the immediate area, and with habitat cohesion at medium scales, whereas diversity was negatively correlated with increasing mean patch habitat area across larger surrounding areas. Orchids co‐existed with exotic species but were negatively affected once exotic cover exceeded 20%. Species composition was correlated with exotic cover. Our findings reveal a complex relationship between orchid communities and their surrounding environments suggesting orchids benefit from a somewhat disturbed environment at patch and landscape scales. These idiosyncratic responses suggest orchid diversity may be unreliable as early‐warning indicators of habitat disturbance.


| INTRODUCTION
Orchid species constitute the largest family of plants including over 700 genera and 27,000 species (Chase et al., 2015;Govaerts, 2016) occurring globally, except for the poles and extremely dry deserts (Jones, 2006;Tsiftsis, 2020). Due to the complexity of their symbiotic interactions with other species (e.g., dependence on specific associations with mycorrhizal fungi and pollinators ;Hutchings, 2010;Phillips et al., 2020) and their narrow niche breadth and high levels of geographic endemism (Dixon & Phillips, 2007;Orejuela-Gartner, 2012), the consequences of global change may be particularly severe for orchid communities (Gale et al., 2018).
Despite these threats, land-use change, including land clearing and habitat fragmentation, are recognized as the main threats to orchid diversity, having contributed to significant population declines in several parts of the world (Vogt-Schilb et al., 2015;Wraith & Pickering, 2019). Therefore, configuration of the surrounding landscape (e.g., vegetation heterogeneity, area of suitable habitat, edge length, and degree of cohesion) is likely to influence species assemblages (Fischer & Lindenmayer, 2007). As such, many orchid species are at a higher risk because they are confined to habitats that have been fragmented or undergone loss of area (Brummitt et al., 2015). Moreover, fragmented patches are more likely to have poor ecological condition and increased levels of disturbance (McIntyre & Hobbs, 1999). Habitat fragmentation, together with increasing ecological disturbance can lead to the interruption of interspecific interactions such as the ones established by orchids (Faast, 2010;Parra-Tabla et al., 2011;Pauw & Bond, 2011;Pauw & Hawkins, 2011;Phillips et al., 2015), pollination failure, and decreased reproductive success (Aguilar et al., 2006;Newman et al., 2013). Habitat fragmentation can indirectly foster weed invasion as well as enhance grazing intensity via affecting the proportion of edges around patches and accessibility for wild herbivores (Faast & Facelli, 2009). Additionally, proximity to roads and urban nuclei enhances the risk of orchids being collected and harvested illegally further threatening the conservation of orchid diversity in remnant habitats.
The rapid decline of orchid communities worldwide suggests urgent conservation measures are required (Phillips et al., 2020). Of the orchid species that have already been assessed for the IUCN Global Red List to date, more than half have been included under some form of threatened or near-threatened category (IUCN, 2021;Wraith & Pickering, 2018). Thus, from a conservation perspective, orchids constitute a key but complicated group given they display a great variety of environmental requirements and pollinator specificities (Reiter et al., 2016;Swarts & Dixon, 2017), for which process-based conservation actions or those preserving the habitat should be implemented alongside species-focused ones (Fay, 2018). Therefore, to design successful conservation strategies for this plant family, we need to fully understand orchid dynamics by studying simultaneously multiple variables affecting their survival (Fay, 2018).
Australian orchids provide a sound model for understanding the simultaneous effect of different threats to orchid diversity. Australia is a large and megadiverse continent that contains a high diversity of endemic orchid species (approximately 1800 species, 95% endemic; Backhouse, 2007;Wraith & Pickering, 2019). Orchids are overrepresented among Australia's threatened species, making up 17% of nationally threatened plants (Faast & Facelli, 2007). The susceptibility of orchids to environmental change is believed to be due to a complex set of threatening processes, including small population sizes, habitat modification and fragmentation, invasive species, grazing, and disruption of pollination interactions (Quarmby, 2010;Wraith & Pickering, 2019).
The aim of this study was to unravel the effect of the above factors on orchid community diversity in South Australia. The factors included: (i) habitat configuration in the surrounding landscape (including vegetation remnancy at increasing radii; adapting the method suggested by Parkes et al., 2003), and (ii) local disturbance, including grazing intensity and exotic species cover. Regarding habitat configuration in the surrounding landscape, we hypothesized that greater size of the patches of remnant habitat and connectivity of such patches would correlate with higher orchid richness, abundance, and diversity. Regarding local patch disturbance, we hypothesized that orchid diversity would be higher in plots with lower levels of disturbance, that is where ecological condition is good, and there is low intensity of grazing and low cover of exotic species.

| Study area
The study was conducted in the Mount Lofty Ranges, which are located to the east of Adelaide, the capital of South Australia, comprising 11,208 km 2 (Daniels & Good, 2015), and extending north-south for approximately 300 km (Robinson et al., 2018). The Mount Lofty Ranges are a climatic refugium at the continental scale due to their heterogeneity of climatic and topographic conditions (e.g., they display a wide elevation range within a relatively small area) (Byrne, 2008;Crisp et al., 2001;Guerin et al., 2016;Guerin & Lowe, 2013). The area has a Mediterranean climate. Mean annual precipitation ranges from 460 to 990 mm, whereas minimum and maximum annual mean temperatures range from 7.8 to 10.3 C and from 18.1 to 21.7 C, respectively (Harwood et al., 2016), and the maximum altitude is 936 m (Robinson et al., 2018). Geologically, the Mount Lofty ranges have been formed by repeated faulting from the late Cenozoic until the late Pleistocene (Preiss, 2019). The soil is formed by metamorphosed sedimentary rocks from the Proterozoic, with quartzite predominating in the most elevated points (Cochrane, 1963). The Mount Lofty Ranges, like other Mediterranean-climate ecosystems around the world, have been subjected to intensive disturbance, especially since European settlement, including land clearing for agrarian purposes, habitat fragmentation, species introductions, and regular wildfire regimes (Casado et al., 2018;Martín-Forés, 2017;Martín-Forés et al., 2017;Westphal et al., 2007). Originally, the vegetation in the Mount Lofty ranges mainly consisted of dry sclerophyll forests containing several species of the genus Eucalyptus, and eucalypt-dominated grassy woodlands in areas with lower elevation (Cochrane, 1963). At present, only about 8% of the original vegetation remains, with patches of remnant native vegetation mainly located either in private lands, or confined to conservation parks and roadside vegetation (Daniels & Good, 2015).

| Data collection
We selected nine sites covering different ecological conditions and level of disturbance from existing monitoring networks (Bond et al., 2019;Guerin et al., 2014). The sites corresponded to remnant vegetation patches, mainly woodlands dominated by species from the genus Eucalyptus and surrounded with agrarian landscapes. Sampling took place between winter and spring 2019. Depending on the size of the patch, between three and five 30 Â 30 m plots were established at each site, with a total of 39 plots surveyed (Table 1). Plots comprised a gradient of land-use from heavily disturbed to relatively undisturbed areas. We surveyed each plot three times, in mid-winter, earlyspring, and late-spring, as some of the orchid species flower or emerge only in winter. These three surveys ensured the most complete composition sampling possible across the main orchid flowering season, during which many species flower for short, specific periods. Due to logistical constraints, six plots were only visited twice and three were only visited once; thus, a total of 115 surveys comprised the final dataset.
In each plot, pairs of observers searched for orchids along 5-m-wide panels. We recorded and photo-vouchered each orchid species and counted the number of individuals per species per plot. Records were uploaded to the Wild Orchid Watch iNaturalist project for confirmation of species identifications. For each plot, we ranked grazing intensity from 0 (no signs of grazing) to 3 (severe grazing intensity) according to Croft et al. (2009) (Supporting Information Material S1), and we estimated the percentage of cover occupied by exotic species. Exotic species refer to all alien plant species that are naturalized, regardless of the invasive status. We also scored overall ecological condition at plot level based on Keighery Condition Scale (Keighery, 1994; Supporting Information Material S1) from 0 (completely degraded) to 1 (very good), to explore its effect on orchid species richness and abundance.

| Data analyses
We combined observations of orchid community composition, including richness, abundance, and diversity metrics with spatial data. We calculated fragstats from spatial data with the function ClassStat from the SDMTools R package (VanDerWal et al., 2015) and the package raster (Hijmans, 2020). We calculated habitat configuration (Jaeger, 2000), and assessed the effect of the mean patch area of the remnant native vegetation patches, total edge length, and T A B L E 1 Surveyed sites, and geographic coordinates and altitude of the plots within each site habitat cohesion on orchids at increasing spatial scales (i.e., increasing radii around the plots; Hirao et al., 2008).
Habitat cohesion refers to the patch cohesion index proposed by Schumaker (1996) to quantify the connectivity of habitat and it is an area-weighted measure that is proportional to the mean perimeter-area ratio divided by the mean patch shape index or standardized perimeter-area ratio (Mcgarigal, 2015). We performed Kruskal-Wallis tests to look for differences associated with sampling season (i.e., the time during the main orchid flowering season in which the survey took place: winter, early spring, and late spring). Results of differences among the three survey times are reported and discussed in the Supporting Information Material S2.
For orchid species richness and abundance, we conducted a suite of generalized linear mixed models (GLMM) fitted to a Poisson error distribution, whereas for Shannon diversity index we used linear mixed effects models (LMM) fitted to a Gaussian error distribution. To test for the effect of habitat configuration within the surrounding landscape on species richness, abundance and Shannon diversity index, we conducted comparisons of sets of mixed-effects models at four different landscape scales, including circular areas around the plot of 1-, 2-, 5-, and 10-km radii. Fixed effects included mean patch area of the remnant vegetation, edge length, and habitat cohesion (VanDerWal et al., 2015). To test for the effect of disturbance on species richness, abundance and Shannon diversity index, we conducted the same suite of models but including grazing, cover of exotic species, and the interaction between both as fixed effects. For all models, we included plot nested within site and season nested within plot as random effects, to account for both potential spatial autocorrelation and the variation among revisits associated with different surveys conducted in different seasons.
In every case, we extracted the best subset of models based on ΔAICc ≤ 2, and we selected the most parsimonious one (i.e., the one with the lowest number of variables). Afterward, we calculated marginal (i.e., the proportion of variance explained by fixed effects) and conditional (i.e., the proportion of variance explained by fixed and random effects) R 2 for the LMMs (i.e., Shannon diversity models) with the MuMIn R package (Barton, 2020). We estimated R 2 for the species richness and abundance GLMMs with the delta method from a model with a Poisson distribution also with the MuMIn R package. For the best subsets of models, we computed ANOVAlike tables testing the inclusion of random-effect terms in the models; we used the function ranova from the lmerTest R package for the LMMs. For the GLMMs, we compared the models with and without the inclusion on the random effect with an anova function.
The response of orchid species composition to disturbance (including overall ecological condition at plot level, grazing intensity, and percentage of cover of exotic species) was explored by conducting non-metric multidimensional scaling. We used as input a matrix that combined the species abundances records of the three surveys, by keeping the maximum abundance found across the different seasons per species. We calculated Bray-Curtis dissimilarities among plots with the vegdist function of the vegan package (Oksanen et al., 2019). We explored the ideal number of dimensions to consider with the NMDS.scree function from vegan and afterwards set the number of dimensions to three (displaying stress < 0.1).
We conducted NMDS with the metaMDS function setting the maximum random starts to 100. We fitted environmental factors (i.e., grazing intensity and percentage of vegetation cover occupied by exotic species) as correlation vectors onto the obtained ordination with 999 permutations to test significance of associations. We obtained the projections of points onto vectors that had maximum correlation with corresponding environmental variables, the squared correlation coefficient of these vectors with the obtained ordination, and the associated p values.
All statistical analysis and calculations were performed using R (R Core Team, 2020) employing the packages vegan (

| Effect of spatial configuration
We recorded a total of 68 orchid species, all of them terrestrial geophytes. Within-site variability of habitat parameters was in general very low, whereas among sites, the values of some habitat parameters, such as mean patch area of the remnant habitat, differed considerably (Supporting Information Material S3). The effect of habitat configuration in the surrounding landscape differed across the spatial scales considered. The best models for species richness and abundance regarding the effect of habitat configuration in an area of 1-km radius included edge length, having a significant positive effect, whereas for orchid Shannon diversity, the best model was the null one (i.e., no fixed effects included). The best models for orchid species richness and abundance regarding habitat configuration in an area of 2-km radius included a significant positive effect of habitat cohesion. In a surrounding area of 5-km radius, habitat cohesion had a significant and positive effect for the three diversity metrics, mean patch area was negatively related with species richness and abundance, and edge length had a positive effect of orchid abundance. Finally, at greater spatial scale (i.e., 10-km radius) the best models for the diversity metrics included mean patch area, which displayed a significant and marginally significant negative effect on species richness and Shannon diversity index, respectively (Figure 1; Supporting Information Materials S4-S6). The variance explained by the fixed factors of the best models for species richness and abundance increased with increasing radii of the surrounding area, varying from 29% to 87%, and from 34% to 77%, respectively. The inclusion of the random effects improved the model in every case (see details in Supporting Information Material S6).

| Effect of ecological disturbance
The best model for orchid species richness, abundance and Shannon index with regards to the effect of ecological disturbance included exotic species, having a significant negative effect (Supporting Information Materials S4 and S5). However, the variance explained by the cover of exotic species only accounted for 14%-18% for all diversity metrics, with the inclusion of the random effects significantly improving the models (Supporting Information Material S6). Although the effect of exotic species cover was negative overall, the most diverse orchid communities coexisted with intermediate levels of exotic cover (i.e., 15%-20%; Figure 2). When the cover of exotic species surpassed 20%, the negative association with orchids became stronger, with highly invaded plots occurring in three sites displaying few or no orchids (Figure 2).
Similarly, there was one site with very low values of overall ecological condition at plot level (i.e., ranked condition <0.4) that did not host any orchids (Figure 2).

| Floristic differences among sites
NMDS showed that orchid species composition differed among sites but also among replicate plots within them. Within some sites (such as CRO and to a lesser extent BLA; see Table 1 for sites names), there was a lot of variability in orchid species composition among plots, whereas plots within some other sites such as CHA or SAN, were quite similar ( Figure 3). Additional overlap in species composition was associated with season (Supporting Information Material S2).

| Effect of spatial configuration
Contrary to our expectations, the mean patch area of remnant vegetation in the surrounding landscape did not have a positive effect on orchid diversity. The effect of habitat configuration in the surrounding landscape at increasing radii varied depending on spatial scale. At the local scale (i.e., 1 km radius), orchid species richness and abundance were positively correlated with the total edge F I G U R E 1 Representation of the effect of landscape habitat parameters on local patch orchid diversity at different spatial scales and the sign of the effect. The effect of mean area on Shannon diversity at 10-km surrounding area is only marginally significant. See Supporting Information Material S4 for the best subsets of models, and S5 for estimates, statistics, and significance for fixed effects of the best models F I G U R E 2 Legend on next page. length of the remnant habitat. When considering medium-sized surrounding areas (2-5 km radii) orchids were favored mainly by habitat cohesion. At larger surrounding areas (5-10 km radii), we observed a negative effect of the mean patch area of the remnant vegetation on the patch-level diversity.
The positive correlation with edge length at the local scale (1 km radius) can be related to the fact that within a forest matrix, the edges can sustain equal or even more orchid species than the forest core (Parra S anchez et al., 2016). Edges are in general related to ecotone habitats, supporting greater environmental heterogeneity and that in turn is reflected in higher diversity (Senft, 2009). Ecotones have considerably lower levels of competition for edaphic and light resources (i.e., low abundance of dominant plant species and low cover of trees and shrubs); these open niches can therefore be colonized by pioneer species with weak competitive ability or high light requirements such as orchids (Arditti & Ghani, 2000;Fekete et al., 2020).
Additionally, the surveyed areas are close to urban nuclei, within anthropogenic habitats, such as cemeteries, orchards and roadsides; thus, the edges can provide refugial conditions for orchid species (Fekete et al., 2017) and act like ecological corridors by facilitating pollination and seed dispersal between suitable habitats and therefore spread of the orchid taxa (Van Rossum & Triest, 2012).
The positive effect of habitat cohesion at intermediate spatial scales (i.e., surrounding areas of radii 2 km and 5 km) is supported by previous findings. Greater habitat cohesion implies reduced fragmentation and consequently associated disturbance. It is well known that the dynamics of pollinating insects, and therefore the orchidpollinator relationship, as well as the fruit set are positively associated with both the density and the condition of the surrounding vegetation matrix (Akhalkatsi et al., 2014;Parra-Tabla et al., 2011).
At larger spatial scales (i.e., surrounding landscape of 10 km radii), we found an a priori unexpected negative F I G U R E 2 Orchid species richness and individual abundance recorded in 115 seasonal surveys of 39 plots at different levels of disturbance, with regards to (a) percentage of cover occupied by exotic species, (b) overall ecological condition, and (c) grazing intensity. Boxplots and whiskers graphics show the median value (black line in the middle of the box), first and third quartiles (lines limiting the box), and minimum and maximum values (extremes of the whiskers) F I G U R E 3 Non-metric multidimensional scaling, in which the differences in orchid community composition among sites can be differentiated by the colors of the polygons, and differences within sites by the size of the polygon. Vectors show the disturbance variables that were significantly correlated to the ordination space effect of the mean patch area of the remnant vegetation on the diversity of orchid communities at patch-scale. It has been previously stated that peripheral orchid populations can perform better than populations occurring in the center of vegetation patches, despite being under more stressful conditions (García et al., 2010). Therefore, larger patches of remnant habitat could diminish orchids' potential to find edges or semi-disturbed conditions in which they might perform better. In this sense, it has been observed that habitat fragments of modest size can be ideal to maintain functional pollinator populations (Cane, 2001) which might indirectly lead to greater orchid abundance. Thus, maintaining modest size patches of remnant habitat ensuring good connection among them (e.g., close to each other to promote high cohesion while maintaining high values of total edge length), could be a good strategy to support more diverse orchid communities.

| Effect of ecological disturbance
As expected, high disturbance negatively affected orchid diversity, with exotic species cover >20% being associated with decreased orchid diversity, although the plots with the highest cover of exotics were found within a small number of sites. However, orchid communities appeared to be most diverse in patches that had some level of disturbance.
High cover of exotic species had a negative effect on all orchid diversity metrics, whereas low cover seemed to not negatively affect orchid communities (Figure 2), despite exotic plant species ostensibly posing a threat to 65% of the orchid species present in Australia (Nevill, 2010;Wraith & Pickering, 2019). In addition, cover of exotic plant species had a significant relationship to the orchid community composition, possibly related to the ability of the different orchid species to compete for resources (e.g., light, water, pollinators) with generalist exotic plants, or some additional factor of disturbance that affects both. Exotic plants are often pollination generalists (Bartomeus et al., 2008), dominating the diet of pollinators (Morales & Traveset, 2009;Pyšek et al., 2011), and achieving high reproductive success (Thompson & Knight, 2018). In addition, in cases where exotic plants also form mycorrhizal associations, native orchids can become displaced due to overlapping webs of fungi (Bonnardeaux et al., 2007). Moreover, generally weeds are indicators of a history of high anthropogenic disturbance (Bonilla-Valencia et al., 2021), with greater presence of exotic species being often associated with greater past or present grazing intensity (Wraith & Pickering, 2019). Therefore, the orchid species that are able to survive in more invaded places could be benefiting from grazing being mainly directed to palatable exotic grasses.
Regarding the effect of grazing, our results only showed a weak positive trend between orchid richness and abundance and grazing intensity (attributable mainly to rabbits at highly disturbed sites and Western Gray Kangaroos or Euros elsewhere; Figure 2). Apart from that, we found no significant evidence that orchid communities in South Australia are negatively affected by level of grazing intensity, despite herbivores having an increasing grazing impact on protected areas of native vegetation in the study region (Prowse et al., 2019).
Regarding overall ecological condition at plot level, we observed a minimum threshold of ecological condition below which orchid species were not observed ( Figure 2). Previous studies showed that poor overall ecological condition is frequently associated with habitat modification and fragmentation, high levels of solar irradiance, and increased levels of grazing and exotic cover (Akhalkatsi et al., 2014;Wraith & Pickering, 2019). Despite this, some of the orchid species restricted to Mediterranean climate ecosystems are known to be able to grow in habitats that appear disturbed to a certain degree (Martin-Fores, 2017;Tsiftsis & Antonopoulos, 2017) or even to benefit from post-fire conditions (Bradshaw et al., 2011;Coates et al., 2006;Lamont & Downes, 2011). Nevertheless, good ecological condition (e.g., areas with greater native vegetation and that have not been affected by fire, dieback, or intensive grazing) is generally related to a diminished proportion of orchids at risk of extinction (Wraith & Pickering, 2019).

| Orchids as ecological indicators
The combined effects of biological invasions, anthropogenic disturbance, and habitat fragmentation on orchid communities are complex to interpret and our results suggest that orchids as a family are not ideal as a general or early warning ecological indicators. This can be related to the disparity of traits that species within this family display (Taylor et al., 2021). For example, orchids can have different environmental and habitat preferences (e.g., orchids that prefer open areas typical from grasslands versus orchids that prefer forested areas), tolerance to disturbance (e.g., disturbance opportunistic versus disturbance-intolerants; Collins & Brundrett, 2015), or reproduction strategy (e.g., some orchids species are highly clonal).
Preventing degradation of ecosystems and promoting good ecological condition and low exotic species cover is crucial to preserve orchids. However, the recommendations toward habitat configuration appear unclear, and need to be carefully considered at different spatial scales in order to inform management strategies. The orchids recorded in this study appear to be favored by edges at local scale, and small-sized, well-connected habitat patches at medium and large spatial scales. Previous literature suggested that the effect of habitat configuration in the surrounding landscape can differentially influence orchid population viability depending on the pollination syndrome displayed by the different species and on their specificity with the pollinator involved (Newman et al., 2013). Further studies should aim to elucidate the effect of habitat configuration on orchid diversity controlling for the species' pollination syndrome, clonal ability, environmental requirements, and tolerance to disturbance. Unraveling responses of different functional groups within the orchid family to habitat fragmentation, biological invasions, and disturbance would be necessary as a preliminary step to be able to use orchids as reliable indicators of the suitability of a certain habitat. This will also help to design improved conservation strategies adjusted to the configuration of the area to be managed and the specificities and preferences of the orchid community.

ACKNOWLEDGMENTS
This work was supported by the grant 327-2018 -Australian Orchid Foundation-Climate and habitat condition controls on orchid populations-research outcomes associated with a citizen science program, and the grant from the Department of Industry, Innovation and Science, Australian Government-Inspiring Australia, Citizen Science Grants-Wild Orchid Watch. The authors thank the TERN Ecosystem Surveillance team, Rosalie Lawrence, Robert Lawrence, Penny McLachlan, Katie Irvine, Michael Starkey, Sally O'Neill, Candy Guerin, Anthelia Bond, the District Councils of Mount Barker, and Alexandrina and the South Australian Department for Environment and Water. They would also like to thank two anonymous reviewers for their insightful comments.

CONFLICT OF INTERESTS
The authors declare that there are no conflict of interests.

AUTHOR CONTRIBUTIONS
Irene Martín-Forés and Greg R. Guerin conceived the ideas. Greg R. Guerin and Samantha L. Bywaters collected the data. Irene Martín-Forés analyzed the data and drafted the paper. All authors reviewed the final version of the manuscript and approved its submission.

DATA AVAILABILITY STATEMENT
The dataset will be published in a public repository upon the acceptance of the manuscript. The dataset associated with this publication is publicly available on Figshare https://doi.org/10.25909/19443821.v1 ORCID Irene Martín-Forés https://orcid.org/0000-0003-3627-0347