Urban aliens and threatened near-naturals: Land-cover affects the species richness of alien- and threatened species in an urban-rural setting

Urbanisation has strong effects on biodiversity patterns, but impacts vary among species groups and across spatial scales. From a local biodiversity management perspective, a more general understanding of species richness across taxonomic groups is required. This study aims to investigate how fine-scale land-cover variables influence species richness patterns of locally threatened and alien species. The study was performed in Trondheim, Norway, covering a steep urbanisation gradient. Spatially correlated Generalised Linear Mixed Effects Models predicting the number of all-, threatened-and alien species by taxon, habitat, habitat heterogeneity and mean aspect within 500 m×500 m grid cells were constructed. The habitat categories were based on detailed land-cover maps. The highest number of threatened species was found in habitats relatively less affected by humans, whereas the number of alien species were only dependent on taxonomic group and spatial correlation. It is shown that land-cover variables within an administrative border can be used to make predictions on species richness within overarching species groups. Recommendations to biodiversity management agencies are to ensure protection of natural habitats to favour locally threatened species, and closely monitor urban areas to mitigate the introduction and spread of alien species.


Materials and Methods
Study area. The study was carried out within Trondheim Municipality (Norway) administrative borders, around 63.42°N, 10.38°E (Fig. 1a,b). It is a southern-boreal 43 , coastal municipality with an area of 342 km 2 , a population of approximately 190,000 people 44 , and annual mean temperature and precipitation are approximately 5 °C and 887 mm 45 . The municipality holds a steep urbanisation gradient; from the city centre and industrial areas, through rural areas including agricultural areas and commercial forests, and to nature reserves and areas managed for biodiversity conservation. The municipality covers highly different nature types, including both coastline, subalpine areas and limnic systems, and thus has a high potential for varied biological communities and high levels of biodiversity 44 . Trondheim municipality is fairly well-sampled with regards to species occurrence records, e.g. due to the presence and activity of the University Museum.
Data retrieval and data cleaning. Land cover data. Land cover was based on the Norwegian AR5 maps (Land Resource map 1:5000) from NIBIO 46 . Shapefiles of the land cover maps were provided by the Trondheim Municipality in April 2018. The AR5 maps are both continually and periodically updated, and provides the most complete data on national land resources 47 . Land cover is categorised based on land cover type, tree cover type, timber productivity and soil condition, giving 66 functionally unique categories within Trondheim municipality (hereafter called "land cover types") (Supplementary Material 1, Table S.1). The map was overlaid by a 500 m × 500 m grid.
Updates of the AR5 maps are mainly done if the categorical classification of an area changes, and the responsible authorities are notified of this change 46 . Consequently, "unannounced" changes are not reflected in the data. As the land cover data was matched with GBIF records from 2013-2018, changes within this period are not taken into account.
Aspect of the terrain was retrieved from a Digital Terrain Model raster with a resolution of 25 m × 25 m. The circular aspect (unit: degrees) was transformed to a "northness"-measure by =°Northness cosin Aspect ( ( ) ) , hence fitting a scale of −1 to 1 (in this definition: −1 = south-facing, 1 = north-facing). The values were rescaled to a gradient from 0 to 1 to match the scale of other included variables. All flat areas were given NA-values. For each grid cell, northness was calculated as the mean of all raster cells within the overlaid grid cell.
GBIF occurrence records. Large amounts of data on species occurrences are available from online databases, such as the Global Biodiversity Information Facility (GBIF) 48 , and the Norwegian Biodiversity Information Centre 49 .
All occurrence records from a bounding box around Trondheim Municipality (the exact municipality border was too detailed to include in the process) were downloaded from GBIF on 06/03-2018 50 (864,715 records in 1. Records containing a full species name for comparability with the threatened-and alien species lists. 2. Coordinate uncertainty of ≤354 m, (1/2 length of the diagonal of 500 m ×500 m grid cells). 3. Records made between January 1 st 2013 and March 6 th 2018 to ensure compatibility with the used land cover maps, and a negligible amount of land cover change.
Of these records, 94.9% were within the kingdom Animalia (91.4% of the total data set were birds), 3.7% within Plantae, and 1.3% within Fungi. 0.65‰ (163 records) were from outside these kingdoms (Supplementary  material 2, Table S.2). The data set was divided into threatened-and alien species (only including animals, plants and fungi).
The "threatened" status was defined based on one or more assessments of the national Norwegian Red List from 2006, 2010 and 2015, provided by The Norwegian Biodiversity Information Centre. See Supplementary material 3 for detailed description of inclusion details.
The "alien" status was based on the Alien Species List (v. 2012 51 ) from The Norwegian Biodiversity Information Centre. Only species alien in mainland Norway were retained (excluding species alien only to Svalbard). All alien species were included, regardless of risk category. Discrepancies in nomenclature between GBIF records and species lists were resolved using the Taxonomic Name Resolver (function "tnrs" from the taxize-package 52 ). Only terrestrial and limnic species were included in the data sets. All species classified as marine by The Norwegian Biodiversity Information Centre were manually excluded from the lists (excluding birds; all bird species in the data set were regarded as terrestrial). 32 Table S.2). The risk of species mis-identification is considered negligible, as the majority of records are associated with organisations deemed reliable regarding species identification (e.g. the Norwegian Ornithological Society, the Norwegian Botanical Society and the NTNU University Museum herbarium). Furthermore, as individual species are not analysed, it is unlikely that single erroneous records will affect the aggregated species pool.
The number of threatened-and alien species, and the overall species richness, registered within each grid cell was calculated, and divided into five taxonomic groups: birds, non-avian animals, plants, fungi and other taxa. "Other taxa" was excluded from further analyses due to a low amount of data. preparation of variables. Land cover variables. To reduce the number of land cover types while avoiding subjectively defining categories, hierarchical cluster analysis was used to identify grid cells with similar composition, creating a limited number of clustered land cover type categories. All grid cells within the administrative border of the Trondheim municipality were used for the cluster analysis, including cells only partially within the www.nature.com/scientificreports www.nature.com/scientificreports/ municipality border, including only the within-municipality area. Marine grid cells (entirely covered by ocean) were not included, resulting in 1509 grid cells in total.
The cluster analysis was done using the function "hclust" on a dissimilarity matrix based on the AR5 land cover in each grid cell, using "Complete linkage" as the clustering method, and a Bray-Curtis dissimilarity matrix of the individual grid cells (function "vegdist", package vegan 53 ). Cut-off value was set at height=0.99 (referring to the height of the cluster-tree, where height=1 indicates no clustering, and height=0 each individual branch (grid cell) being an autonomous cluster), resulting in 17 clusters in total, of which 6 included ≤3 grid cells. The clusters including ≤3 grid cells were mainly found on the municipality border. These were excluded from further analysis. Each cluster will hereafter be referred to as a "habitat".
The habitats were named according to the (on average) dominating land cover types within the cells (  Table 1). The most frequent habitat within the municipality was Cultivated, followed by Coniferous forest and Urban/developed areas.
Outliers in the number of records or number of species (evaluated separately for each taxon level) were excluded based on Tukey's method (0.75 quantile + 1.5*IQR). Subsequently, the habitats Open firm ground and forest and Open firm ground and cultivated land were excluded from the analyses, as only one and two grid cells remained, respectively. 485 grid cells were included in the subsequent analyses (Fig. 2).
Habitat heterogeneity was calculated for each grid cell as the Simpson's Diversity Index. The index is calculated , and n is the total area of a particular land cover type, N is the total area of the grid cell. The index ranges between zero (completely homogeneous land cover) and one (infinite heterogeneity in land cover; a hypothetical value). The grid cells included in the analyses ranged between 0.012 and 0.884.

Statistical analyses.
Generalised Linear Mixed Effects Models (GLMM) were constructed, predicting the threatened-(Poisson error distribution), alien-(Poisson error distribution), and overall (negative binomial error distribution due to overdispersion) number of species in each grid cell by habitat, habitat heterogeneity, northness, including an interaction with taxon (birds, non-avian animals, plants, fungi and other taxa) for all variables. Total number of records within grid cells were used as offset to account for differences in sampling effort. To account for spatial autocorrelation (Moran's I > 1 in exploratory Generalised Linear Models), a Matérn correlation function was used as a random effect (package spaMM 54 ). The models were fitted using Maximum Likelihood. Model selection was performed as stepwise backwards selection, based on AIC on the full models of the form:

No species habitat taxon habitat heterogeneity taxon northness taxon
Matern longitude latitude ( 1 ) . = * + * + * + + . The models were subsequently used to predict species richness across all grid cells within the Trondheim municipality, using 100 records as an offset.
All data preparation, analyses and figures were made in R, version 3.6.1 55 .

Results
Different models proved to be optimal for the three species groups (all-, threatened-and alien species). For overall species richness, all predictors and interaction terms were retained, whereas threatened species richness was predicted by habitat, northness and taxon. Alien species richness was only predicted by taxon (Tables 2-4  total species richness. For overall species richness, northness had a negative effect on species richness, the magnitude varying with taxon ( Table 2, Fig. 4): non-avian animals responded most strongly to northness, followed by fungi, plants and birds. The response to habitat heterogeneity varied by taxon: plants and birds responding positively to increasing levels of habitat heterogeneity, fungi and non-avian animals having a negative response (Fig. 5). Similarly, the response to habitat differed among taxa, all other variables being held constant at mean values: fewest birds are predicted in Open marsh and coniferous forest followed by Coastal, Freshwater, Coniferous forest; high production and Urban/developed. The highest number was predicted for Cultivated, followed by Urban/vegetated/riparian, Coniferous forest; low production and -medium production. However, 0.95 C.I. overlapped for all habitats. For non-avian animals, Cultivated and Urban/developed had lower predicted species richness compared to Open marsh and coniferous forest, and Cultivated was lower than Coniferous forest; low production and -medium production as well. All other 0.95 C.I. overlapped. The highest number of fungi species was predicted for Open marsh and coniferous forest, 0.95 C.I. only overlapping with Coniferous forest; high production. The lowest number was predicted for Coastal, 0.95. C.I. overlapping with Freshwater, Cultivated, Urban/vegetated/ riparian and Urban/developed. The lowest number of plants was predicted for Freshwater, followed by Urban/ developed (0.95 overlapping with Urban/vegetated/riparian, Cultivated and Coniferous forest; high production). The highest number was predicted for Open marsh and coniferous forest, 0.95 C.I. overlapping with Coastal, Coniferous forest; medium production, -low production and Coniferous forest; high production ( Fig. 6). threatened species richness. For threatened species, increasing values of northness increase the predicted number of species (Table 3, Fig. 4). The highest species richness values are found for birds, followed by non-avian animals, fungi and plants. However, 0.95 C.I. overlap for all taxa except for birds and plants in Urban/developedand Cultivated areas. The highest numbers of species are found in Open marsh and coniferous forest, followed by Cultivated, Coastal, Freshwater, Urban/developed, Coniferous forest; low production, -high production, Urban/ vegetated/riparian, and Coniferous forest; medium production. However, all 0.95 C.I. overlap (Fig. 6).
Alien species richness. For alien species, only taxon was retained as a predictor; the highest number of species predicted for plants, followed by non-avian animals, fungi and birds. However, the 0.95 C.I. overlapped for all taxa except birds and plants (Table 4, Fig. 6).

Discussion
Urban areas are often found to have high levels of biodiversity, but little is known on how fine-scale land use is structuring species diversity in cities. We used species occurrence records from GBIF and official land cover classifications to determine how habitat affects total species richness, and the number of threatened and alien species. We did so by constructing spatially correlated Generalised Linear Mixed Effects Models based on habitat, habitat heterogeneity, aspect and taxonomic group within 500 m × 500 m grid cells across the Trondheim municipality, selecting the best models based on ∆AIC. The best models varied for overall-, threatened and alien www.nature.com/scientificreports www.nature.com/scientificreports/ species richness, with total species richness depending on all predictors and their interaction with taxon, whereas threatened species richness depended on habitat, aspect and taxon, and alien species richness only depended on taxon. The relationship between species richness in general are highly complex and dependent on multiple factors and interactions (Table 2, Figs. 3-6). Threatened, native species are associated with non-anthropogenic habitats ( Table 3, Figs. 4 and 6), whereas alien species are mainly affected by spatial correlations on the investigated spatial  www.nature.com/scientificreports www.nature.com/scientificreports/ scale ( Table 4). The key findings of this study advance our understanding of the field by confirming the association of threatened, native species with more natural habitats, and the potential for establishment of alien species across all habitats on a management-relevant spatial scale.
The retention of all predictors and interactions in the model of overall species richness illustrate the complex relationships between environmental variables and different taxonomic groups. Nevertheless, the overall negative effect of northness reflect the low species richness of north-facing slopes, compared to south-facing ones 42 (Fig. 4). The different taxa responded differently to increasing habitat heterogeneity, the only unidirectional response being for plants (positive) (Fig. 5) 56 , in which respectively habitat heterogeneity and habitat richness were positively associated with species richness in urban areas. However, other studies have found a positive correlation for restricted taxonomic groups, such as arthropods 29,40 , birds and mammals 30 , which was not observed here.
Non-surprisingly, the different taxa responded differently to various habitats (Fig. 6). Interestingly, whereas plants, fungi and non-avian animals generally responded negatively to urban areas (differences not necessarily significantly different from other habitats however), the effect was less pronounced for birds. This could reflect their high mobility, and potential for an "urban adapter/exploiter"-status of some bird species 13 . In contrast, the habitat with the highest predicted number of both plant-, non-avian animal-, and fungi species, had the lowest predicted number of bird species.
Threatened species richness generally responded positively to increasing northness, in contrast to what would be expected (Fig. 4). This could potentially be an artefact of the habitat associations; Coastal areas had higher northness values (mean = 0.758, S.E. = 0.022 compared to the overall mean = 5.35, S.E. = 0.005). The effect of taxon reflects the differences in the number of species within each taxonomic group being classified as threatened; (50 bird species, 26 plant species, 12 non-avian animal species, 33 fungi species included in the study). For all taxa, the lowest species numbers are predicted for all variants of coniferous forest, contrary to the initial expectations, and urban areas (Fig. 6). The negative effect of urban areas on threatened species richness mirrors the findings of Aronson et al. (2014) 8 , and emphasises how vulnerable native species are not pre-adapted to the changed environments of the city. Contrary to expectations, the effect of the various forest habitats on threatened species is lower than for most other habitats ( Table 3). The low number of threatened species in forests can be due to the lack of sampling, showing a spatial bias in the data rather than an effect. This should however be accounted for by using the number of records as an offset in the models. Rather, large parts of the forested areas in Trondheim are srongly affected by previous afforestation for timber production, where mainly coniferous species (both native and alien) were planted 33 . These forests might not provide the needed conditions for native species 57 . Plantations and secondary vegetation have been shown to harbour fewer species than primary forests 58,59 . The lack of association between forested areas and threatened species calls for a nuanced perspective on what forest types constitute suitable habitats for species of interest, as indicated by Ingram et al. (2015) 58 and Horák et al. (2019) 59 . The highest species numbers are predicted for Open marsh and coniferous forest and Coastal areas (Fig. 6); the former is likely the habitat category reflecting the lowest human impact. The high number of threatened species in coastal habitats can likely be ascribed to these habitats being ecotones, providing varied habitat conditions. Ecotones have    38 found ecotonal species to mainly be natives, which is supported by the findings here.
Interestingly, in the model of alien species richness, only taxon was retained as a significant predictor, reflecting the differences in the number of species within each taxonomic group being classified as alien (5 bird species, 156 plant species, 10 non-avian animal species, 6 fungi species included in the study). The lack of response to either of the other investigated variables stands in stark contrast to the expectations and previous findings, but can be attributed to alien species often being generalist opportunists; the spatial scale investigated does not reflect  Fig. S.2), it is evident that founder events and subsequent spread of alien species are of crucial importance: on the investigated scale, even more important than the configuration of environment. As many alien species are introduced through urban areas mainly due trade and traffic 12,15,31 , emphasis must be put on the www.nature.com/scientificreports www.nature.com/scientificreports/ importance of avoiding unintentional introduction of potential invasive species. As an example, the review by Kowarik (2011) 5 found cities to be hotspots of alien plant species. In addition, port cities have been suggested as even greater hotspots of introductions, leaving Trondheim even more vulnerable 60,61 . www.nature.com/scientificreports www.nature.com/scientificreports/ As the explanatory variables used in these models are "indirect" (sensu Guisan and Zimmermann (2000) 62 ), the habitats are proxies for underlying environmental (direct) drivers. Therefore, a direct extrapolation to other geographical areas should be cautious 62 . However, the general methods are highly applicable elsewhere.
Of the 1,509 grid terrestrial cells, 485 qualified for analyses; species occurrence data was sparse in the rest. Those used in the analyses were biased towards urban areas (Table 1), supporting the general trend in citizen science data; concentrated around inhabited or areas otherwise accessible to the public 63,64 . For example, areas within Trondheim municipality relatively far from human activity are under-sampled, with two habitats not being represented in the analyses at all (Open firm ground and forest -and cultivated land). This bias was accounted for in the models, but the differentiated sampling effort nevertheless leaves varying degrees of uncertainty for each habitat and taxon. The sample sizes differed among species groups, with many more observations of threatened than alien species. The differences thus might reflect sampling strategy rather than reality.
As the models are by nature rather crude, they inevitably lack predictor variables, which could have increased model accuracy. However, including highly detailed variables was not the aim of this study. Since the data set included a wide array of species, these will not respond in similar ways to variation in the included variables, or to missing variables 65 . The more species included in the models, the more opposing mechanisms are attempted to be fitted within a single modelling framework, giving a poorer fit, compared to models with a narrower scope.
The number of GBIF records have increased in recent years (see Speed et al. (2018) 64 ). Of all species recorded in Trondheim, approximately 1/3 have been recorded within the municipality from year 2013 to 2018. Of the 6,020 species from the downloaded data set not included in the analyses, 33.9% (2,039) have only been recorded once, and 85.5% (5,150) have been recorded <10 times. Most of these infrequent species are insects. This taxonomic skew is likely due to this species group being poorly sampled or requiring expert knowledge to identify to species level.
Different correlations with environmental variables are expected at different spatial scales for different organisms 19,23,66 . Taxa and species with opposing responses to the included variables could mask each other, thus not revealing the underlying mechanisms 56 . Simultaneously, the mechanisms underlying species distributions vary with spatial scale, not necessarily in the same direction for different taxa 19,23,67 . As multiple different taxonomic groups were included in this study, the used spatial scale is potentially inappropriate for all taxa.
According to Pautasso (2007) 19 , a negative correlation between urbanisation and species richness is expected when the study grain is smaller than 1 km, as in this study, but positive at larger scales. This is ascribed to the larger scale reflecting human settlement in productive areas, competing for space with other species, whereas the small scales reflect more detailed environmental-and land cover effects.
Our results indicate that if the Trondheim municipality is to be managed to favour biodiversity, favouring threatened species and excluding alien species, the following actions can be recommended: Habitat heterogeneity on a relatively small spatial scale should be ensured, favouring overall species richness. This should however not be confused with fragmentation of natural habitat.
To favour threatened species, non-anthropogenic-and coastal areas should be monitored and protected, potentially expanding the actions to ecotones in general.
To limit the spread of alien species, initial introduction and establishment should be avoided. Thus, urbanand other anthropogenic areas should be closely monitored and managed 12,68 .
Protection of important and heterogeneous habitat should be accounted for in unison with ensuring large habitat patches, rather than multiple smaller ones; a metastudy by Beninde et al. (2015) 56 showed patch area to have the largest positive effect on urban biodiversity.

conclusions
Overall-, threatened-and alien species richness are not determined by the same land-cover variables. Totaland threatened species richness responds to both habitat and aspect, whereas alien species richness does not respond to any of the variables included in this study. The highest numbers of threatened species are associated with non-anthropogenic habitats, but in contrast to expectations, not more positively associated with forested areas than other habitats, calling for detailed investigations of the importance of forest environments for threatened species. Our finding that alien species do not respond to land-cover variables, but only spatial correlations, confirms the importance of founder events, and highlights the status of cities as gateways for alien species in general.
To mitigate the knowledge gaps from under-sampled habitats, we urge for sampling outside inhabited areas and for less investigated taxa. Using models build on administrative land cover maps and open database occurrence records can be a useful tool for local biodiversity management, by providing guidelines regarding where to aim future efforts, both regarding future conservation efforts and future investigations. Further work is however needed in dealing with the inherent biases of such databases.
In the case of Trondheim, an averaged sized Northern European city, the recommendations for biodiversity management are to ensure protection of natural habitats within a spatial resolution of 250,000 m 2 , and to closely monitor and manage urban areas to mitigate the introduction and spread of alien species.   Table 4. Model output from the spatially correlated GLMM of alien species richness. The model was constructed with a Poisson error structure. The factor level Aves is used as intercept, thus categorical predictor values are relative to this.