The global loss of floristic uniqueness

Quantification of changes in floristic similarity

To quantify changes in floristic similarity by naturalized flowering plant species, we extracted regional lists of alien species from the Global Naturalized Alien Flora (GloNAF) database⁴⁵ and regional lists of native species from the Global Inventory of Floras and Traits (GIFT) database⁴⁶. The GloNAF database contains lists of naturalized vascular plant taxa for 861 regions (countries or subnational administrative units), ranging in size from 0.03 to 6,864,961 km² (median size is 15,152 km²) and covering >80% of the terrestrial ice-free surface globally⁴⁷. GloNAF includes 13,803 plant taxa that, according to the original data sources, are alien plants and have established self-sustaining wild populations in the respective regions (i.e., are naturalized⁵). The GIFT database is a compilation of floras and checklists of predominantly native vascular plant species with an indication of their floristic status for more than 300,000 species across nearly 3000 regions with near global coverage⁴⁶. We first selected regions that matched perfectly between GloNAF and GIFT. Additionally, we merged some GloNAF regions to match a larger GIFT region, and vice versa, by comparing the overlapping area of nested regions using the R package ‘sf’ (version 0.8-0)⁴⁸.

To ensure the highest data quality, and to be on the conservative side, we restricted our analysis to regions with complete or nearly complete checklists of both native and naturalized alien species. For GloNAF, we only included regions for which there was at least one species list judged to include more than 50% of the naturalized taxa for that region⁴⁵. Although the judgment of species-list completeness is coarse and for most lists made by the GloNAF curators, it allows the exclusion of regions for which the data are obviously poor. For GIFT, we included a region only if at least one species list aimed to represent its entire native angiosperm flora. Our strict selection criteria resulted in a dataset including native and naturalized species for 658 non-overlapping regions, including 154 island regions, 503 mainland regions and one region including both islands and mainland areas (Chile). These regions covered all continents, except Antarctica, but there was low coverage for parts of Africa and Asia (Fig. 4).

We restricted our analyses to flowering plants (angiosperms), which had the most complete species lists, and to species with accepted names in The Plant List²⁴ (http://www.theplantlist.org/). We excluded species with an uncertain native/alien status or with a conflicting status, i.e., being native to a region according to GIFT but being alien to the same region according to GloNAF. Furthermore, since the native/alien status of many infraspecific taxa and hybrid taxa are less clear, we restricted our analyses to the species level (i.e., infraspecific taxa were assigned to the binomial species name), and we excluded hybrids. Our final dataset included 1,139,254 native species-by-region records for 189,110 species and 141,762 naturalized species-by-region records for 10,130 species.

For all 216,153 possible pairwise combinations of the 658 regions, we quantified the taxonomic and phylogenetic similarities between their native floras (SimTax_native, SimPhyl_native), and between their floras including both native and naturalized alien species (SimTax_{native+naturalized}, SimPhyl_{native+naturalized}). As the regions vary largely in species richness (ranging from 11 to 13,720 species with a median of 1704), we used the Simpson similarity index for taxonomic similarity (Eq. 1)⁴⁹, which is largely insensitive to species richness:⁵⁰

Here a is the number of species common to both regions, b is the number of species that occur in the first region but not in the second and c is the number of species that occur in the second region but not in the first⁵¹. Likewise, we calculated the Simpson phylogenetic similarity index as phylogenetic similarity (Eq. 2) as implemented in the R package ‘betapart’ (version 1.5.1)⁵²:

Here A is the total length of the phylogenetic branches in the phylogenetic tree that are shared by the species of both regions, B is the total length of the phylogenetic branches that are shared only by the first region and C is the total length of the phylogenetic branches that are shared only by the second region⁵¹. To quantify changes in similarity due to naturalization of alien species, we calculated the degree of homogenization H (or differentiation, see below) for each pair of regions as

A small value of 0.001 was added to both similarities to avoid infinite values. A positive log-response ratio indicates homogenization (i.e., increased floristic similarity between two regions), and a negative one indicates differentiation (i.e., decreased floristic similarity). As an alternative to the Simpson similarity index, we also calculate the Sørensen similarity index, which additionally takes into consideration the nestedness of the floras in the paired regions⁵¹. As the results were not sensitive to the choice of similarity indices (Supplementary Fig. 14), we focused our analyses on the Simpson similarity index.

To quantify phylogenetic similarity, we used a phylogenetic tree including all angiosperms with accepted names in The Plant List (Supplementary Fig. 2). The tree was developed based on the mega phylogeny of Smith and Brown⁵³. We added missing species (n = 71,124, of which 733 are naturalized in other regions) with their accepted names in The Plant List to the root of their genus or families. For details on the development of the phylogenetic tree, see ref. ⁴⁷.

Quantification of geographic distances and climatic distances

We calculated the pairwise geographic distance between regions as the distance between their geographic centroids using the R package ‘geosphere’ (version 1.5-10)⁵⁴. We also calculated the nearest distance between the geographic borders of regions. However, since the distances between geographic centroids are highly correlated with distances between region borders (n = 216,153, r = 0.996, P < 0.001), we only used distance between region centroids in our analysis.

We quantified the pairwise climatic distance between regions as the distance between their positions in multidimensional climate space. We extracted all 19 biologically relevant variables of temperature, precipitation and seasonality from the WorldClim database⁵⁵ (Supplementary Table 1) at a resolution of 2.5 arc-min. As some of these bioclimatic variables are highly correlated, we first conducted a PCA on them, and used the first three principal axes, which are orthogonal to each other, as new climatic descriptors. Many of the bioclimatic variables had a skewed distribution and varied greatly in magnitude. Therefore, before including them in the PCA, we transformed each variable to be as approximate to a normal distribution as possible (Supplementary Table 3) using the R package ‘normalizer’ (version 0.1.0)⁵⁶. We then scaled them to have a mean of zero and a standard deviation of one. For each 2.5 arc-min grid cell, we obtained its scores along the first three principal component axes, which in total explained 85.9% of the variance in the original bioclimatic data (Supplementary Fig. 15a). We then calculated the climatic centroid of cells extracted for each region within the three-dimensional climatic space, and quantified climatic distance between two regions as the Euclidean distance between their climatic centroids. In addition, as the first principal component axis (PC1) was mainly related to temperature variables, and PC2 mainly to precipitation variables (Supplementary Fig. 15b), we also calculated distances for PC1 and PC2 separately.

Compiling data on administrative relationships

To assess current and past administrative relationships among the 658 regions, we first assigned each region to its current country, and then compiled data on the past colonial relationships between the 110 countries represented by our regions. For colonial relationships, we only considered the period after 1492 (i.e., the year of discovery of the Americas by Columbus), which is the time period during which the major introductions and naturalizations of alien species happened⁵⁷. As basis for our current and past administrative relationship dataset, we used the TRADHIST dataset⁵⁸. Since TRADHIST is not comprehensive and describes only the administrative relationship between countries or states in the last two centuries, we added data on colonial relationships for the main colonial empires as listed in the Wikipedia article ‘Colonial empire’ (https://en.wikipedia.org/wiki/Colonial_empire). We next extracted data on relationships of dependent territories, which are defined as territories that do not possess full political independence as a sovereign state but remain politically outside the controlling state’s integral area⁵⁹ (e.g., Guam [USA]) from the Wikipedia article ‘Dependent territory’ (https://en.wikipedia.org/wiki/Dependent_territory).

We were able to group administrative relationships of each pairwise combination of the 658 studied regions into three main categories: (i) same country: the two regions are currently part of the same country (e.g., California and Texas, both part of the USA); (ii) dependency relation: the two regions are either currently dependent territories or past colonies of the same country (e.g., Hong Kong and Ireland, both were historically colonies of the UK), or the two regions belong to two different countries of which one is or was the dependent territory or colony of the other (e.g., Guam [USA] and England; California [USA] and England); (iii) no administrative relation: the two regions have no current or past administrative relationships.

Assessing the association of the change in floristic similarity with geographic distance, climatic distance and administrative relationship

To describe the non-linear relationships of taxonomic and phylogenetic floristic similarities (SimTax_native, SimTax_{native+naturalized}, SimPhyl_native, SimPhyl_{native+naturalized}) along the gradients of geographic distance and climatic distance, we fitted single-predictor log-binomial generalized linear models (GLMs) following ref. ¹⁷. The intercept of the model with geographic distance as the predictor was fixed at 1, assuming complete similarity at a distance of 0 km. Following ref. ^16,17, we calculated and compared the halving distance (i.e., the distance at which a given similarity value is predicted to have decreased by 50%) of each of the four similarity indices.

To statistically test how changes in taxonomic and phylogenetic similarities (i.e., the degree of homogenization or differentiation) between two regions vary with geographic distance, climatic distance, administrative relationship, and their interactions, we used multiple regression on distance matrices (MRM)⁶⁰. To account for non-independency of data points, caused by use of each region in multiple region pairs, the statistical significance of the MRM-model-coefficient estimates was assessed by comparing them to the null distribution of the coefficient estimates⁶⁰. The latter was produced by simultaneously shuffling the rows and columns of the response matrix (i.e., degree of homogenization between regions). This permutation was done 999 times (see Ref. ⁶⁰ for an example of the permutation process). As the MRM models assume linear relationships, we additionally used Generalized Additive Models (GAMs) to check whether any possible nonlinear effect of geographic distance or climatic distance would change our conclusions from the linear MRM model. Since the results of both models were qualitatively consistent with regard to the main effects, we present only the results of the linear MRM model in the main text, and the GAM results and some deviations with regard to the interaction effects in Supplementary Fig. 5. To assess whether the effects of climatic distance were mainly due to temperature or precipitation variables, we ran additional MRM analyses in which we replaced the single climatic distance measure with distances based on PC1 and PC2, respectively (Supplementary Fig. 15b).

Assessing homogenization hotspots and associations with characteristics of the regions and their floras

To assess where the global homogenization hotspots are, and how homogenization values depend on characteristics of the regions and their floras, we first averaged for each region all of its pairwise homogenization values H. We then extracted a set of variables characterizing the regions and their floras, including the species richness and phylogenetic diversity of native and naturalized floras, the proportion of endemic species, the donor score of the region (calculated as the average number of non-native regions each native species is naturalized in), the size of the region, and whether the region is an island or a mainland region. We then did a linear multiple regression of the average degree of homogenization of the regions on the variables characterizing the regions and floras. Of these predictors, species richness and phylogenetic diversity, the donor score, and the region size were log-transformed to increase their linear association with homogenization. To compare their relative importance, predictors were scaled to have a mean of zero and a standard deviation of one. We reduced spatial autocorrelation in the residuals of the linear model by adding a spatial autocovariate that incorporates a matrix of longitude and latitude coordinates of the centroids of the regions with the R package ‘spdep’ (version 1.1-3)⁶¹ (Supplementary Fig. 16). All data extraction, statistical analyses and figures were done using R version 4.1.0⁶².

Reporting summary

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

Source: Ecology - nature.com