Assemblage reorganization of South African dragonflies due to climate change

Climate change is expected to cause large shifts in species assemblages such as dragonflies and damselflies (Insecta: Odonata). Here, we assess the influence of environmental drivers of turnover on Odonata assemblages. Secondly, we map the predicted spatial variation in species composition, first as a gradient of assemblage similarity, and then as discrete bioregions delineating major areas of odonate endemism. Finally, we map the magnitude of expected change in species turnover in response to climate change under two emission scenarios.


| INTRODUC TI ON
Understanding the drivers of species distributions is essential for implementing effective conservation strategies for climate change adaptation, protecting regional diversity and ecosystem function, and preventing further species loss (Socolar et al., 2016). While considering policies aimed at biodiversity loss, there is a parallel need to better understand the candidate mechanisms linked to the drivers of ongoing and future biodiversity change. Global indicators continue to show rapid declines in the state of biodiversity worldwide, with shifts in species distributions and community structures already observed in a variety of taxa (Newbold et al., 2015). Species turnover may reflect the potential breakdown of local assemblagelevel community composition, with poorly known consequences.
For instance, many species are threatened by the loss of climatically suitable habitat (Bush, Nipperess, Duursma, et al., 2014;Krechemer & Marchioro, 2020;Niskanen et al., 2019). Changes in climate suitability is predicted to increase dispersion pressure in the future and may result in the substantial reorganization of many assemblages (Flenner & Sahlén, 2008;Wauchope et al., 2017).
In the light of these concerns, what approaches might best inform conservation planning and action? The collective diversity of a region cannot be adequately measured by site-level diversity indices, such as alpha and gamma diversity, as they do not explain patterns of spatial and temporal turnover (McGeoch et al., 2019). When assessing the impact of environmental change, knowing which species are present (i.e., species composition) is more informative than simply knowing how many species are present (i.e., species richness) (Hillebrand et al., 2018;Jones et al., 2017). Compositional turnover, or beta diversity, captures how fast local assemblages change their species composition over space and time (Socolar et al., 2016). The term beta diversity was first introduced by Whittaker (1960) and is defined as the "extent of change in community composition" among sites. Indices of beta diversity can therefore reveal how biodiversity is organized in space and over time (Socolar et al., 2016). While variation in local community composition is typically thought to be driven by species sorting along often spatially autocorrelated environmental gradients, other spatial processes such as dispersal limitation can further affect local community composition (Cottenie, 2005).
Consequently, when explored together with environmental and geographical gradients between sites, beta diversity can further reveal the underlying drivers of compositional turnover. A better understanding of large-scale patterns in community turnover and how these may shift with time can thus provide conservation-relevant insights about the maintenance and management of species over large areas.
Historical taxonomic biases in conservation focus have targeted charismatic large-bodied vertebrates, undermining the conservation of invertebrate taxa (Donaldson et al., 2016). These biases have resulted in invertebrates being inadequately managed and overlooked when facing challenges such as climate change. For instance, rising temperatures and a shift in rainfall are predicted for South Africa, which will likely threaten insects such as dragonflies (Odonata: Anisoptera) and damselflies (Odonata: Zygoptera), and for simplicity hereafter jointly termed dragonflies (Samways, 2010).
Dragonflies are widely used for freshwater assessments (Kietzka et al., 2017;Nasirian & Irvine, 2017;Schindler et al., 2003;Simaika & Samways, 2011). Just under 20% of the national fauna of South Africa's dragonfly species is endemic to the country (Samways et al., 2016). As a result of their narrow habitat tolerances and localized geographical ranges, endemic species are likely to be more susceptible to global changes than widespread species (Simaika & Samways, 2015;Vincent et al., 2020). In the face of increasing environmental change such as climate, the survival of dragonfly species will continue to be threatened unless intervention measures are put in place by conservation agencies to ensure the longevity of threatened species and cohesion of local assemblages (Nóbrega & De Marco Jr, 2011). As such, to assess how best to conserve biodiversity, and to predict how local assemblages may reshuffle due to environmental change, we first need to understand the processes that maintain diversity at regional spatial scales. With such knowledge on the drivers of regional diversity at hand (Deacon et al., 2020), we are equipped to explore anticipated community reorganization under different future climate change scenarios.
Changes in local temperature and precipitation regimes alter habitat suitability and consequently drive local and regional species turnover (Bush et al., 2013;Shah et al., 2012). Range-shifts in dragonflies are an adaptive response to this change in climate (Flenner & Sahlén, 2008;Hassall et al., 2007;Hickling et al., 2005). A biome covers regions with a similar climate, vegetation, disturbance regime and geography and will therefore also influence species composition as the disparity in resource and niche availability between biomes can drive selection (Stein et al., 2014). As all species within the order Odonata are dependent on water for their egg and larval phase, they are also particularly sensitive to water condition (Samways et al., 2016). This makes dragonfly assemblages excellent indicators of water quality change for aquatic ecosystem assessments (Kietzka et al., 2017;Schindler et al., 2003;Simaika & Samways, 2011).
Dragonflies have also been successfully used as model organisms to explore impacts of climate change (Bush et al., 2013;Hassall & Thompson, 2008;Simaika et al., 2013). Furthermore, research has indicated high congruence of odonate species richness and patterns of threat with that of other taxonomic groups, like mammals, birds and amphibians, which implies that odonates are an effective surrogate taxon (Darwall et al., 2011;Pinkert et al., 2020). Odonates could therefore act as a useful indicator for other taxonomic groups.
This study has three key objectives. Firstly, to assess the importance of environmental drivers of turnover on Odonata assemblages in South Africa. By relating species turnover of local assemblages to environmental gradients and geographical distance, we aim to uncover what processes are responsible for assemblage composition.
These relationships can also provide an indication of how the rate of compositional turnover varies along the gradients concerned. As the relationship between compositional turnover and environmental distance is not constant over the entire range of environmental gradients, we determine not only what covariates are important for maintaining turnover but also at which points along the gradients is turnover most sensitive. We achieve this using generalized dissimilarity modelling (GDM), which models spatial gradients of compositional turnover between pairs of geographical sites (Ferrier et al., 2007). Our second objective is to map predictions of compositional turnover across the entire region, first as a continuous gradient of assemblage similarity, and then as discrete bioregions depicting major areas of endemism based on clustering results under present climate conditions. For these present-day bioregions, we assess what species are representative of each bioregion by compiling a species list for conservation purposes, as well as a list of those narrow-ranged species that exist in only one of the identified bioregions. Thirdly, with space-for-time substitution (Blois et al., 2013), we estimate and map the magnitude of expected change in species turnover in response to climate change under two future climate conditions (i.e., RCP 4.5 and RCP 8.5) for years 2050 and 2070. With this, we also assess the stability of present-day bioregions under these future climate change scenarios and develop an index of bioregion instability. Results from this study carry important implications for conservation planning by providing information regarding biodiversity drivers, assemblage persistence and change as well as ecological integrity.

| Study region
Our study is focused in South Africa at the southernmost tip of the African continent. South African climate ranges from a Mediterranean type in the southwest, temperate in the interior, and subtropical in the northeast. Average annual rainfall is c. 450 mm but varies considerably from west (100 mm) to east (900 mm) (see Figure S1-BIO12 in Appendix S1), with some desert regions along the western edge of the country. The weather is largely influenced by ENSO (El Niño-Southern Oscillation) events, which result in cyclical drought and flood periods (Dilley & Heyman, 1995). Published predictions of global warming are cause for concern as even a small temperature increase will place considerable stress on South African insect populations (Samways, 2010).

| Odonata dataset
The spatial dataset of South African Odonata, described in detail by Samways (2009) andSamways et al. (2016), is based on a compilation of monitoring and museum collections and sightings data from the 1900s to present, as well as records published in Pinhey (1984Pinhey ( , 1985. Records are also continually updated with the most recent maps published in the handbook on freshwater assessment for South Africa (Samways et al., 2016). At the time of writing, the dataset comprised 164 species and over 25,000 presence records throughout South Africa, around 17,000 of which were collected after year 2000. South Africa is one of the most significant centres of dragonfly endemism in the world as it contains most relictual species on the African continent (Clausnitzer et al., 2012;Simaika et al., 2013). The dataset's large sample size and expert verification of all insects collected makes this an ideal dataset to focus efforts to predict assemblage reorganization of South African dragonflies induced by climate change.

| Environmental covariates
The choice of climatic variables should aim for ecological relevance (Araújo et al., 2019). Environmental variables were therefore selected based on the large body of odonate literature focused on understanding dragonfly assemblage structure. For example, climate can be considered a dominant factor driving insect distributions (Bush et al., 2013;Hassall & Thompson, 2008). The life cycle of dragonflies is primarily governed by temperature and water availability that are essential for larval development and breeding. Annual temperatures as well as temperature range are also known to affect dragonfly physiology, such as phenology, seasonal regulation, immune function and pigment production used for thermoregulation (Hassall & Thompson, 2008). Precipitation of the driest quarter was selected as prolonged periods of drought are known to drastically impact larvae survival (Pernecker et al., 2020). Thus, we extracted four present-day bioclimatic variables, with two describing air temperature and two describing precipitation (www.world clim. org; Hijmans et al., 2005), over the period from 1970 to 2000. These climatic variables comprise annual mean temperature -"BIO1," temperature annual range -"BIO7," annual precipitation -"BIO12" and precipitation of the driest quarter -"BIO17." Secondly, data characterizing river conditions were obtained from the National Freshwater Ecosystem Priority Areas (NFEPA) database available through SANBI's Biodiversity GIS portal (bgis.sanbi.org/ SpatialDataset; Council for Scientific & Industrial Research, 2011).
Because dragonflies are susceptible to changes in the quality of aquatic systems, dragonfly response has been used extensively as an indicator of ecological integrity (Simaika & Samways, 2011). River condition was categorised into five classes: unmodified/largely natural with few modifications -"AB," moderately modified -"C," largely modified -"D," extremely modified -"EF" or tributary condition modelled as no longer intact -"Z" (see Table S2 in Appendix S1).
Length of each river condition was calculated per grid cell (cells described in detail below under generalised dissimilarity modelling).
Thirdly, we included human population density data for 2010, obtained from NASA's open data portal (beta.sedac.ciesin.columbia.edu/data; Ciesin, 2016). Population density was used as human presence (and density thereof) is largely associated with biodiversity loss. We explored a long list of other anthropogenic factors, such as human footprint and human influence index, but found they were largely correlated with human population density (Variance Inflation Factor >7) and yet less intuitive to interpret and were therefore not included.
Lastly, we sourced data on vegetation types, including the percentage cover of each biome (n = 9, see Appendix S1) within each geographical grid cell for the year 2011 (bgis.sanbi.org/SpatialDataset; South African National Parks, 2011). We specifically used vegetation types as we expect species-specific habitat requirements to drive dragonfly diversity and distribution patterns (Schindler et al., 2003).
Indeed, vegetation gradients have been found to be an important driver of assemblage patterns (Osborn & Samways, 1996). org). In total, four sets of variables (two RCPs ×two periods) were included to predict dragonfly assemblage reorganization under future climates in South Africa. Climate predictions from an additional General Circulation Model (GCM), namely CNRM-CM6-1, were also considered, although analyses using this GCM produced rather similar results (see Appendix S3 and Figure S10) and therefore were not discussed in the main text.

| Generalized dissimilarity modelling
Generalized dissimilarity modelling (GDM) is formulated as a nonlinear extension of the more traditional distance approach of linear matrix regression (Ferrier et al., 2007). It is a statistical method for analysing and predicting patterns of compositional turnover in species assemblage, usually in response to environmental gradients that vary in space. GDM uses the compositional dissimilarity between all site pairs and fits these dependent variables as a function of the difference in environmental and geographical covariates between these site pairs (Ferrier et al., 2007). Unlike in the case of species distribution models (SDMs), GDM does not target the presence and absence of individual species but rather species turnover of the entire assemblage between sites. For this reason, SDMs typically require a comprehensive geographical coverage of records for estimating the representative niche of the focal species; in contrast, GDM requires a representative list of species for sites of interest to compute realistic level of turnover in the model.
All analyses were carried out in r (version 3.6.1; R Development Core Team, 2019). South Africa was divided into 1,978 quarter degree cells in accordance with national atlas data resolution standards of previous projects (Harrison & Cherry, 1997). Each cell was considered a site and characterized by its central geographical coordinates; species composition (presence and pseudo-absence) and associated environmental covariates were summarized using the "LetsR" r package (Vilela & Villalobos, 2015). Since species presence/absence can be sensitive to sampling effort, we removed cells with low sampling counts (less than ten records/observations, which is the median number of records per quarter degree cell across South Africa) ( Figure 1) leaving 351 quarter degree cells containing 22,552 dragonfly records to model representative dragonfly assemblages (Appendix S2). To further assess the role of sampling effort on observed local species composition, we included the total number of Odonata records in each cell as a surrogate for sampling effort in our analyses; this means that a higher number of Odonata records in a grid cell is likely to reflect a higher survey effort and thus yield a higher number of discovered species. Increased sampling efforts will yield overall a more representative species list of local assemblages (Hortal & Lobo, 2005) (Appendix S2; Figure S2) and will ideally reduce issues related to sampling inadequacy (e.g., pseudoabsences; Hui et al., 2011). Geographical distance was calculated as orthodromic distance between pairs of cells to account for changes in turnover resulting from the decay of compositional similarity with distance (Morlon et al., 2008;Soininen et al., 2007). This distance decay can be a product of multiple mechanisms, such as ecological filtering through spatially autocorrelated environmental gradients (the Moran effect; Hansen et al., 2020) and dispersal limitation .
All covariates showing excessive skewness were transformed and only those where this transformation increased data symmetry were retained (See Appendix S1 for detail; Figure S1 for transformed variables). Variance Inflation Factor (VIF) measures how much the variance of a regression coefficient is inflated due to multicollinearity within the model. All covariates were tested for multicollinearity at a quarter degree scale by computing the VIF using the "car" r package (Fox et al., 2012). In general, a VIF value that exceeds five or in some cases ten indicates a problem with multicollinearity (Kutner et al., 2005). To be conservative we removed all variables with a VIF exceeding five, which was only the percentage cover of grasslands within the biome data. The final selection comprised 20 explanatory variables, including geographical distance and sampling effort, which were used as environmental predictors in the GDM (see Appendix S1 for data layer sources, projections and transformations as well as Table S1 for a summary of the twenty covariates used in the model).
The differences in species composition between two sites were expressed as a single ecological distance, the Jaccard distance metric. To explain compositional dissimilarity between pairs of cells, we transformed each covariate into five monotonic I-splines with three internal knots using generalized linear models embedded in the GDM (using the "gdm" R package; Manion et al., 2018). The number of Isplines determines the precision of pinpointing the response along the environmental gradient represented by a covariate, while the number of knots determines the number of interlinked spline pieces and was kept low to avoid overfitting (Perperoglou et al., 2019). The percentage contribution of each explanatory variable to the variation of Jaccard dissimilarity was then estimated from 20 separate models (removing one predictor at a time). The variance explained by the full GDM model was then compared to each of the partial models to evaluate the importance of each explanatory variable (Table 1). The uncertainty in the fitted I-splines was assessed using bootstrapping of 70% of the sample site pairs with 50 bootstrap iterations (Appendix S3; Figure S5). The absolute model error of the observed versus predicted partial ecological distance was calculated and mapped for each cell (Appendix S3; Figure S6). A set of Bayesian Bootstrap Generalized Dissimilarity Models (BBGDMs) were also run to assess the robustness of the response curves for key model predictors (see Figure S9 in Appendix S3) (Woolley et al., 2017).
Results from the GDM were presented as response curves, where firstly the height of the response curve of each predictor provides an indication of the total amount of compositional turnover associated with a specific environmental gradient. Secondly, the slope of the response curve at a specific point provides an indication of how the rate of compositional turnover varies along the environmental gradient concerned. Additionally, to assess whether estimated turnover reflects true turnover (i.e., species gain and loss) or nested structure (i.e., narrow-ranged species are geographically nested within the range of wide-ranged species), we partitioned dissimilarity between all site pairs into total beta diversity, true turnover and nestedness components, using the Sørensen dissimilarity index, the Simpson dissimilarity index and the difference of the two (sensu Soininen et al., 2018; see Figure S8 in Appendix S3).

| Assemblage composition mapping
Once the GDM was fitted to observed data, we used the fitted model to predict dissimilarity for unobserved locations in South Africa (cells without dragonfly records, or fewer than 10). For this we used all data layers within the GDM as well as a simulated constant layer representing "adequate sampling" (75 records per quarter degree grid cell) to replace that of real sampling effort. The level of adequate sampling was obtained from the I-spline analysis of the response curve (the number of records when the dissimilarity prediction levels off). This allows us to predict "true" turnover especially in under-sampled regions. We used the "gdm.transform" function, which uses the observed dissimilarity matrix and predictors as an F I G U R E 1 South African map of dragonfly species records at a quarter degree cell level. Data were obtained from the spatial dataset of South African Odonata, described in detail by Simaika and Samways (2009) and Samways et al., 2016. Number of Odonata records in each grid cell is used as a surrogate for sampling effort in our analyses, with a higher number of records in a grid cell representing a higher survey effort and thus a higher number of discovered species. Cells with less than ten records were removed in the analyses input to produce a predicted dissimilarity matrix (N-by-N in dimension), with its element depicting pairwise dissimilarity between sites. For dimensionality reduction and visualization, we used the Principal Component Analysis (PCA) of the predicted dissimilarity matrix and assigned the first three principal components (PCs) to a Red-Green-Blue (R-G-B) colour palette for mapping (Koubbi et al., 2011;Ramette, 2007). Principal Coordinates Analysis (e.g., Picazo et al., 2020) and PCA are both popular methods; however, as we are mapping the predicted dissimilarity matrix, not the observed dissimilarity matrix from co-occurrence matrix, we opted to use the PCA. When mapping the PCs, we weighted each PC to reflect the different proportions of variance explained by each PC (see the subsection below titled "Current and possible shifts in dragonfly assemblage composition relative to bioregion"). Two areas with similar colours show similar species compositions, whereas two areas with different colours represent different species compositions.
Discrete bioregion boundaries were delineated by classifying the first three PC values (represented as three bands: R-G-B) using the r packages "vegclust" (De Cáceres et al., 2010) and "diceR" (Chiu & Talhouk, 2019). Specifically, we first identified eight cluster prototypes from the present-day R-G-B scenario using a simple K-means clustering algorithm (De Cáceres et al., 2010). The significance of these cluster results was then assessed using the SigClust K-Means algorithm (Chiu & Talhouk, 2019). The choice of eight clusters/bioregions was based on exploratory analyses that identified the clearest clusters. For practical application, we further determined what species exist within each of these eight present-day bioregions and selected the top ten species occurring within a particular bioregion as representative species for each (Table 2). In contrast, a species occurring in only one bioregion was considered a narrow-ranged species (see Table S3 in Appendix S4).
To explore anticipated spatial turnover in future climates, we assessed how dragonfly assemblages would potentially respond to future climate change scenarios by re-running the fitted GDM with all previously included environmental layers but with climate layers from RCP 4.5 and RCP 8.5 scenarios for 2050 and 2070 instead of present-day climate layers. The R-G-B of the principal components of these pairwise dissimilarities were then classified using a hard noise clustering (HNC) algorithm where sites were assigned from one to eight prototype clusters (those from the present-day K-means clustering) or a new "noise" cluster based on a threshold distance (0.8) to all cluster centroids (De Cáceres et al., 2010). Each of the outputs was then plotted to produce five maps representing bioregions of the present-day, RCP 4.5 for 2050, RCP 8.5 for 2050, RCP 4.5 for 2070 and RCP 8.5 for 2070. Note, as no new "noise" cluster was found, each future bioregion was linked precisely to one specific present-day bioregion (i.e., bioregions of the same colour in the maps represent the same bioregion shifted in future climate scenarios).
Through this technique, we can view the reorganization of bioregions under different emission scenarios and years. A summary of the above methodology is provided in Figure S7  Lastly, to explore temporal turnover in future climates we utilized the GDM "predict" function, under space-for-time substitution (available in the "GDM" package). In this case, raster predictions were provided for the two time periods of interest (2050 and 2070) and the two climate change scenarios (RCP 4.5 and RCP 8.5) to estimate the magnitude of expected temporal turnover in dragonfly communities within a grid cell in response to environmental change associated with the four selected climate change scenarios. Outputs were then plotted to produce four additional maps.

| Drivers of dragonfly assemblage composition turnover
The observed beta turnover, measured by the Sørensen index, reflects largely true turnover of species gain and loss between sites TA B L E 1 The percentage contribution (R 2 ) of each predictor variable in explaining compositional turnover in the generalized dissimilarly model (the Simpson index) and only slightly nestedness of species distributions ( Figure S8 in Appendix S3). The percentage deviance explained by the full GDM (Pearson's R 2 ) was 35.5% (p < .05). The relationship between observed dragonfly compositional dissimilarity of each site pair and the linear prediction of the regression equation from the GDM (otherwise known as the predicted ecological distance between site pairs) is shown in Figure S3 (Appendix S3). Sampling effort (31%) and annual mean temperature (21%) contributed the most to explain species turnover of South African dragonflies (Table 1).
According to the response curves, annual mean temperature evokes the largest turnover in community composition, followed by sampling effort and then geographical distance between sites measured in quarter degrees ( Figure 2). We find annual mean temperature exceeds all other environmental explanatory variables at explaining turnover in community composition by at least four times, meaning it is the principal abiotic variable driving turnover between sites ( Figure 2a). The slope of this relationship is roughly linear, suggesting compositional turnover increased with increasing annual mean temperature at a rate of c. 8.6% per °C (Figure 2a). For sampling effort, turnover in community composition stabilizes once the effort of 75 records per quarter degree cell is reached (Figure 2b). For geographical distance between sites, at around seven quarter degrees (c. 175 km) compositional turnover also stabilizes ( Figure 2c). This indicates that beyond c. 175 km, complete turnover is driven by environmental gradients and not geographical distance. Compositional turnover is also sensitive to small, as well as large, changes in percentage savanna coverages (Figure 2d). When comparing the various levels of water condition, class D (the largely modified class) had the greatest amount of associated compositional dissimilarity as well as the most rapid rate of change (Figure 2e). Compositional turnover is only affected by temperature annual range above 24℃ and increases further above 30℃ range (Figure 2f). The response of compositional turnover to other predictors is shown in Figure S4 (Appendix S3). In addition, the BBGDM confirmed these response curves but suggested a smaller effect size for sampling effort (Appendix S3; Figure S9).

| Current and possible shifts in dragonfly assemblage composition relative to bioregion
The R-G-B plot of the principal components of the predicted beta dissimilarity matrix illustrates the spatial variability of the dissimilarity of assemblage composition under present climate conditions ( Figure 3b). From the PCA, PC1 explained 58% of the variance and PC2 and PC3 explained 21% and 6%, respectively (total 85%). The compositional gradient represented by PC1 is dominant in the high elevation interior, while PC2 is dominant along the East and PC3 in the mountainous areas all across the escarpment of the country ( Figure 3b). K-means clustering of these three PCs produced eight distinct bioregions (p < .001) according to the SigClust K-means test of significance (Chiu & Talhouk, 2019). Figure 3a shows the nine South African vegetation biomes according to Mucina and Rutherford (2006). Figure 3c shows the identified present-day dragonfly bioregions derived from eight distinct clusters of compositionally similar dragonfly assemblages. Areas of similar assemblage composition are distinctly visible as unique colours (Figure 4c-g). Each bioregion differs in size and shape as it corresponds to a clear turnover/transition of dragonfly assemblages. We can consider the regions of transition between the eight dragonfly assemblages as ecotones or boundaries between unique regional dragonfly assemblages.

| Major drivers of compositional turnover
The trends visible in the I-spline response curves in Figure 2 give us insight into the manner in which turnover varies along the gradient concerned. Our results show that turnover in dragonfly community composition increased with increasing annual mean temperature at a rate of c. 8.6% per °C. As temperature governs development rates and is known to affect dragonfly physiology (Hassall & Thompson, 2008), the behaviour of dragonflies and their survival depend largely on temperature (Gilbert & Raworth, 1996). As such, we see a steady effect of temperature on compositional turnover.
Therefore, it is likely that even small changes in annual mean temperature will cause notable responses in dragonfly communities, such as range shifts and consequently species turnover.
The other major predictor of observed species turnover was sampling effort, making it a confounding factor that could strongly hinder our ability to accurately predict spatial patterns. By including sampling effort in our model, we were able to compensate for areas with lower sampling effort (sampling inadequacy) and determine whether sampling had been done adequately enough to predict Of all considered variables annual mean temperature and sampling effort had the greatest relative contribution to variance explained in compositional turnover. Water condition, only contributed marginally to turnover (Table 1; Appendix S3, Figure S4) In the response curve for geographical distance, the slope is steepest at lower distances but reaches a plateau around seven quarter degree cells (c. 175 km) (Figure 2c). This suggests that turnover beyond this distance is restrained, and notable distance decay of similarity only occurs within 175 km. This is most likely a reflection of disparity in niche requirements and the mobility of Odonata species. Sites that are closer to each other often share similar environmental characteristics and thus harbour a similar species composition, whereas sites that are further apart with distinct environmental characteristics will harbour distinct assemblages. The plateauing of geographical distance also indicates that beyond this geographical range of 175 km, the establishment and composition of dragonfly assemblages cannot be modified by distance-related processes, such as the dispersal capacity of dragonflies within the regional environmental context. Under the space-for-time substitution (Blois et al., 2013), this range thus reflects the limit that a dragonfly community can shift to retain its compositional similarity in response to climate change.
Using historical data of British Odonata, Hickling et al. (2005) demonstrated a shift of 74 km over 36 years between 1960 and 1995, or an annual range shift of merely 2.1 km. This corresponds to the upper dispersal ability of some extensively studied Zygoptera (Conrad et al., 1999;Watts et al., 2004). Some Anisoptera in Namibia, however, have shown the potential for longer dispersal distances (Suhling et al., 2017), such as Phyllomacromia contumax, Phyllomacromia picta and Olpogastra lugubris, which are also present in our dataset. Collectively, 175 km could mark the upper limit of dispersal distance, in the context of regional ecological and environmental heterogeneity, to contribute to compositional turnover; beyond this distance, environmental differences could be solely responsible for observed community dissimilarity. Dispersal limit is a complex phenomenon dependant on species traits, population structure and habitat complexity (e.g., river and habitat configuration), as well as the time-scale represented by records within grid cells, so this result of the distance decay of similarity reflects the emerged interplay of multiple factors and must be interpreted in this context, not necessarily a reflection of a species' inner dispersal ability.

| The future of dragonfly assemblages under climate change
We found mean annual temperature to be the greatest environmen- is expected to exceed 4℃ in the same period. Considering Figure 4, our GDM predicts a magnitude of 64% change in species turnover in coastal regions and up to 78% turnover inland for the RCP 4.5 scenario. For RCP 8.5, we see even greater magnitudes of species turnover, that is between 65% and 80%, with the greatest magnitude shown in the year 2070 for RCP 8.5. The specific hotspots of concern are in the arid centre of the Northern Cape, the border between the North West and Limpopo as well as the far north of the Free State. Interestingly, these hotspots of concern mirror the findings of Coetzee et al. (2009)

on the turnover of birds across Southern
Africa. They also identified higher rates of species turnover in the central arid areas of South Africa. This is also consistent with other findings from the same region on a broad range of taxa (plants, mammals, reptiles, invertebrates;Broennimann et al., 2006;Erasmus et al., 2002;Foden et al., 2007). These hotspots of species turnover are typically described as having a unique climate compared to the rest of South Africa, with no or few species endemic to these areas will likely be retained should temperatures continue to rise (Coetzee et al., 2009). Dornelas et al. (2014) showed that on a global scale, many species assemblages are rapidly changing, with an estimated mean of 10% turnover of montane species per decade. Gibson-Reinemer et al. (2015) predicted a similar level of 12% species turnover per decade in Asia, Europe, North America, South America and the Indo-Pacific as climates continue to warm. Species turnover is likely to be even greater for ectotherms, such as insects. This presents a great challenge for conservation planning as species assemblages that previously occurred in protected areas will not stay intact as species move in response to climate change, in different directions and at different rates. This will in turn also create new species interaction networks that may result in highly complex, nonlinear and sometimes extreme responses (Walther, 2010). This further emphasizes the demand for adaptive management strategies, as many species that occur in areas today will emigrate or worse, become extinct under future predicted temperature increases. Indeed, according to  and also our results, major shifts in the distribution of Odonata are expected to occur because of climate change, which presents major challenges not only for Odonata but for conservation of freshwater biodiversity in general.
This complex, nonlinear response is echoed in Figure 3. These shifts indicate that even at the level of bioregions, 36% of South Africa will likely undergo a complete bioregion shift in the near future, signalling a highly unstable future for dragonfly assemblages.
As bioregions represent areas of similar species compositions, we should expect major future alterations to present assemblages. Overall, our results across all scenarios indicate that a rapid shift in dragonfly assemblages is likely to occur between now and 2050 in South Africa. Such assemblage reorganization is predicted due to either rapid species turnover (Flenner & Sahlén, 2008) or shifts in range margins (Hickling et al., 2005) from ongoing climate change.

| Practical application of the bioregions
One may want to focus on species that are representative of a bioregion for tracking shifts in assemblage structure. These species that show the greatest sampling counts in each of their bioregions should be regularly assessed as early indicators of potential shifts in larger assemblage structures (Table 2). Narrow-ranged species, on the other hand, should be targeted for conservation efforts as they lack redundancy between bioregions. For example, cluster 4 holds five narrow-ranged species that occur only within this bioregion (Appendix S4; Table S3). This cluster therefore maintains a highly unique assemblage of species and should be targeted for conservation and reassessment efforts. Not all the species investigated are native however and should be tagged as an additional point of concern for conservationists. Aciagrion dondense is one such species in cluster 4. A recent colonizer into South Africa, possibly driven by floods farther North, and a possible facilitator of the local demise of Agriocnemis ruberrima (Samways, 2010), also sampled within this cluster.
The findings presented here have direct implications for the selection of protected areas for conservation purposes in future. A principal goal of conservation planning is to establish a network of protected areas that represent the largest variety of species (Linke et al., 2007). In this case, bioregions ( Figure 3 and Table 2) can be used to either establish appropriate reference sites or to define areas in which randomized sites for sampling should be selected (Hawkins et al., 2000). With the dragonfly's ability to reflect changes in temperature, these monitoring programmes can further be used to track forefronts of climate driven changes to other species distributions as well. where the environment may be suitable (Bush, Nipperess, Duursma, et al., 2014). Furthermore, alternate clustering techniques as well as different RCP scenarios could lead to different cluster outcomes.

| Modelling limitations
If the distance dragonfly species need to move in response to climate change exceeds their ability to disperse, many species may need to adapt or face extinction. Our model also assumed that species responses in the future are the same as they are now, which when coupled with adaptation may not always be the case (Stoks & Cordoba-Aguilar, 2012). Furthermore, we used a GDM which applies pairwise beta diversity as a measure of species turnover and thus only captures turnover between pairs of sites. Future research could extend this pairwise analysis to a multi-site biodiversity assessment, such as zeta diversity (Hui & McGeoch, 2014), which may provide a more informative approach to explaining the response of rare versus common species. Lastly our model assumptions were that climate will shift while our other variables remain the same in future years.
To improve future models, we can (with the availability of such data) also include predicted change of these other variables such as river condition and human density. This will likely increase the robustness of the model prediction.

| CON CLUS ION
Our results show that the turnover in dragonfly community composition increased with rising annual temperature at a rate of 8.6% per °C. This indicates that South African dragonflies are vulnerable to the effects of climate change, specifically based on their sensitivity to changes in temperature. We also found that with sufficient sampling these insects may be valuable large-scale biological indicators of climate change.
When we mapped the predictions of compositional turnover across the region and delineated eight discrete bioregions based on compositional similarity of predicted dragonfly assemblages, there was clear turnover/transition of dragonfly assemblages expected across South Africa. These spatial clusters of similar assemblages are predicted to wax and wane in future climates, and similar patterns at a global scale are yet to be explored. Across these bioregions we found 12 species out of the 164 total were narrow ranging species, existing in only one of the eight present-day bioregions (Table 2).
This list of narrow ranging species will be valuable for targeted conservation efforts as it highlights the lack of redundancy in these 12 species across the region. So too, can the representative species list of each bioregion be used to establish long-term monitoring sites, select new reserve boundaries and help track potential shifts in the distributional limits of some species (Appendix S4, Table S3).
When assessing temporal turnover, we detected the largest changes to be expected between now and the year 2050, with predictions reaching up to 80% turnover, suggesting a need for the strategies of rapid response to be developed through coordinated conservation efforts. Indeed, when assessing how present-day bioregions may dissolve and shift, we observed major shifts between now and 2050, and less between 2050 and 2070. Present-day bioregions are highly unstable and may show nonlinear and rapid changes in assemblage compositions. Under all emission scenarios, we can expect far-reaching consequences from the large-scale reorganization of dragonfly assemblage due to climate change. Facing a rapidly changing climate, conservation planning that aims to mitigate climate change impacts and species loss, needs to be mindful of these anticipated shifts and reshuffles as currently demarcated protected areas may not be adequate as climate changes.

ACK N OWLED G EM ENTS
This work was supported by the National Research Foundation of South Africa (grant 89967) and the Global Insect Threat-Response

CO N FLI C T O F I NTE R E S T
The author(s) declare(s) that there is no conflict of interest.

PEER R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13422.

DATA AVA I L A B I L I T Y S TAT E M E N T
Original and generated data, as well as R scripts for data analyses, can be found at https://doi.org/10.5061/dryad.4mw6m 90b9.