This is the accepted version of the journal article: Cardador, Laura. «Incorporating phylogeographic information in alien bird distribution models increases geographic extent but not accuracy of predic- tions». Biological invasions, Vol. 24, issue 3 (March 2022), p. 683–695. DOI 10.1007/s10530-021-02673-7 This version is available at https://ddd.uab.cat/record/257159 under the terms of the license 1 Incorporating phylogeographic information in alien bird distribution models 2 increases geographic extent but not accuracy of predictions 3 4 5 Laura Cardador1,2, Pedro Abellán3, Tim M. Blackburn1,4 6 7 1 8 Centre for Biodiversity and Environment Research, Department of Genetics, Evolution, 9 and Environment, University College London, London WC1E 6BT, United Kingdom 10 2 11 Ecological and Forestry Applications Research Centre, Campus de Bellaterra 12 (UAB) Edifici C, 08193 Cerdanyola del Vallès, Spain 13 3 14 Department of Zoology, Universidad de Sevilla, Facultad de Biología, 41012 Seville, 15 Spain 16 17 Institute of Zoology, Zoological Society of London, Regent’s Park, London, NW1 4RY, 4 18 UK 19 20 21 *Corresponding author. 22 Email: l.cardador@creaf.uab.cat; lcardador81@gmail.com 23 24 25 ACKNOWLEDGEMENTS 26 27 We thank all authors and observers who contributed to the datasets on species 28 distributions and phylogeographic patterns, without which this study could not have 29 been possible. This project has received funding from the European Union’s Horizon 30 2020 research and innovation programme under the Marie Skłodowska-Curie Grant 31 Agreement No. 752149. L.C. also received funding from the Beatriu de Pinós 32 fellowship program (funded by the Catalan Government and EU COFUND program nº 33 801370). PA is funded by ‘V Plan Propio de Investigación’ of the Universidad de 34 Sevilla (Spain). The authors acknowledge the use of the UCL Myriad High Performance 35 Computing Facility (Myriad@UCL), and associated support services, in the completion 36 of this work. 37 38 39 40 41 42 43 44 45 46 47 48 49 50 1 51 ABSTRACT 52 53 Species distribution models (SDM) have been proposed as valuable first screening tools 54 for predicting species responses to new environmental conditions. SDMs are usually 55 conducted at the species level, assuming that species-environment relationships are a 56 species-specific feature that do not evolve and show no variability across a species’ 57 range. However, broad environmental tolerances at the species level can encompass 58 narrower and different environmental tolerances for specific lineages or populations. In 59 this study, we evaluate whether SDMs that account for within-taxon niche variation in 60 climate and human-habitat associations provide better fits between projected 61 distributions and real occurrence data for alien bird species than species-level SDMs. 62 Our study focuses on eight alien bird species with established alien populations for 63 which detailed phylogeographic information was available. Similarity in climates and 64 human disturbance conditions occupied by different phylogenetic groups within species 65 was low and not greater than random expectations. Accounting for intraspecific niche 66 variation in SDMs modified the distribution and extent of suitable habitat predicted as 67 susceptible to invasion, but did not result in more accurate model predictions in alien 68 ranges. Until more accurate information on intraspecific variability is available, species- 69 level models can be reasonable candidates. When phylogeographic information is 70 available, the use of the most conservative criterion (i.e. to model both species and 71 lineages on the basis of the actual range) is recommended. 72 73 74 Keywords: alien species, birds, climate, human disturbance, intraspecific niche 75 variation, invasion risks, species distribution models 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 2 101 INTRODUCTION 102 103 Human activities are accelerating the rate and magnitude of changes in species 104 geographic ranges worldwide. Climate change, land-use transformations, habitat 105 fragmentation or environmental pollution, among others, threaten the persistence of 106 several species in formerly suitable locations, leading to species range contractions or 107 extinctions, or pushing species to track new suitable environmental conditions beyond 108 former range limits (Pimm et al. 1995; Vitousek 1997; Parmesan et al. 1999). The 109 increase and development of trade and transport infrastructures, in turn, have 110 accelerated the dispersal, and subsequent establishment, of alien species in places far 111 away from their native range (Blackburn et al. 2009; Hulme 2009). Such changes have 112 notably altered the distribution of species worldwide, homogenizing species 113 assemblages (McKinney and Lockwood 1999; Capinha et al. 2015; Sayol et al. 2021), 114 with detrimental effects on biodiversity and ecosystems. 115 Although biotic interactions and dispersal are important in constraining species 116 ranges, environmental factors can exert a primary role (Gaston 2003; Huntley et al. 117 2007). Species distribution models (SDMs; models that statistically relate observed 118 species occurrences to environmental variables) have thus been proposed as valuable 119 first screening tools for predicting species responses to new environmental conditions in 120 new geographic areas (e.g. invasion risk assessments) or under future environmental 121 scenarios (e.g. global climate-change) based on current species occurrence-environment 122 relationships (Guisan and Thuiller 2005; Araújo and Peterson 2012). SDMs rely on 123 ecological niche theory, which predicts that for relatively recent events such as 124 biological invasions, the environmental niche is expected to be conserved (Peterson 125 2011). SDMs often focus on macroclimatic variables. Additionally, accounting for 126 environmental factors other than climate, such as habitat characteristics and human 127 disturbance, can substantially improve model predictions (Strubbe et al. 2015; Cardador 128 and Blackburn 2020). 129 SDMs are usually conducted at the species level, assuming that species- 130 environment relationships are a species-specific feature that does not evolve and shows 131 no variability across a species’ range. However, widely distributed species often 132 encompass different taxonomic or evolutionary units, which can reflect the existence of 133 ecotypes and locally adapted populations (Smith et al. 2019). Spatial heterogeneity in 134 environments coupled with reduced gene flow can encourage local adaptation and 135 functional differences, leading to divergence in niches among closely related lineages. 136 In fact, recent work has suggested that broad environmental tolerances at the species 137 level usually encompass narrower and different environmental tolerances for specific 138 lineages or populations within the species (Peterson and Holt 2003; Pearman et al. 139 2010). Hence, modelling a species as a single undifferentiated entity may obscure the 140 possibility that these lineages occupy distinct niches and, as a consequence, miss the 141 idiosyncratic response of intraspecific lineages to changing environmental conditions 142 (Pearman et al. 2010; Lecocq et al. 2019). Accounting for intraspecific niche variation 143 in SDMs has thus been highlighted as important for forecasting species range shifts 144 under changing environmental conditions, particularly under global climate change 145 (Pearman et al. 2010; Peterson et al. 2019). In particular, phylogeographic structures 146 have been proposed as a useful proxy to incorporate intraspecific differentiation in 147 SDMs. However, the efficiency of using these proxies in SDMs remains largely 148 unknown, in part because spatiotemporally independent data to test the accuracy of 149 model predictions are often unavailable (Peterson et al. 2019). 3 150 Biological invasions represent unique, large‐scale biogeographical experiments 151 for evaluating model transferability (Liu et al. 2020). Large numbers of alien species 152 have been introduced well outside their native ranges, resulting in geographically 153 independent datasets. However, the relevance of including intraspecific niche variation 154 in invasion risk assessments has not often been addressed (but see exceptions, Strubbe 155 et al. 2015; Godefroid et al. 2016). Many alien species currently arrive at new areas 156 because they are imported as trade commodities (Hulme et al. 2008; Abellán et al. 2016; 157 Cardador et al. 2017, 2019). If individuals coming from different geographic origins 158 belong to different lineages that have particular ecological niches (in terms of climate or 159 human tolerances), their invasion success is also likely to differ across different 160 recipient environments. The omission of intraspecific niche structure from niche 161 modelling exercises may lead to some lineages having little representation in the 162 resulting species models (Pearman et al. 2010; D’Amen et al. 2013). This can lead to 163 underestimation of the climate tolerances of alien species and, as a consequence, their 164 potential for establishment and spread in new environments. 165 Here, we assess the key assumption of distribution modelling theory – that the 166 environmental niche remains conserved across species native ranges – as applied to the 167 bird invasion process. For this, we focused on eight alien bird species with established 168 alien populations, selected because of the availability of robust phylogeographic 169 information, and for which different phylogenetic lineages (i.e., phylogroups) have been 170 identified in previous studies. We explored niche variation in the climatic and human 171 disturbance spaces occupied by the different lineages in the native range in order to 172 assess whether the different phylogeographic lineages or genetic units occupy different 173 niches. It should be noted that the niches considered here relate to the realized niche 174 (occupied niche) and the Grinnellian niche concept: that is, the response of species to a 175 set of non-consumable environmental variables that influence their large-scale 176 geographical distribution (Soberón 2007). As evidence of niche conservatism was not 177 found, we assessed how accounting for intraspecific niche variation in SDMs influences 178 predictions about potential distributions in adventive regions for the whole species. 179 Accounting for intraspecific variation is expected to improve the representation of 180 different phylogroups in models (particularly the scarcest and narrowly distributed in 181 the native range), and thus to produce better fits between projected distributions and real 182 occurrence data for alien bird species than species-level models not considering within- 183 taxon niche structure. 184 185 186 METHODS 187 188 Phylogenetic and occurrence data 189 190 Our study focuses on 8 alien bird species that have established alien populations in 191 different regions of the world and for which detailed phylogeographic information 192 derived from analyses of mitochondrial and/or nuclear DNA sequences was available 193 from the literature (Table 1). We only considered studies reporting a clear definition of 194 intra-specific phylogenetic divisions and covering substantial parts of species’ native 195 ranges. Thiessen polygons were applied to locations with available phylogenetic data to 196 delimit geographic boundaries of different phylogroups (Figs. 1a and S1 in Supp. Info.) 197 (Strubbe et al. 2015). Thiessen polygons define an area of influence around each 198 sampled point (points with genetic data in our study), where every location of the study 199 area (species ranges in our case) is nearer to this point than to all the others. Thiessen 4 200 polygons derived from sampled points where a given phylogroup was identified were 201 merged together to obtain the phylogroup geographic boundaries. We used ArcGIS 10.5 202 for those analyses. Note that for most species, different phylogroups were allopatric 203 (Fig. S1), although cases of partial overlap in geographic ranges also occur. For each 204 species, occurrence data for niche and modelling analyses were compiled from the 205 Global Biodiversity Information Facility (GBIF, GBIF.org, 2017,Table S1). GBIF 206 records spanned the years 1744 to 2017 (Table S2). Compiled records were classified as 207 pertaining to the native breeding range or alien established range according to range 208 maps provided by the BirdLife International & NatureServe (2014) and the Global 209 Avian Invasions Atlas (Dyer et al. 2017), respectively. Note that in the case of 210 Copsychus saularis, BirdLife International & NatureServe (2014) provided separate 211 maps for Copsychus saularis and Copsychus mindanensis. Based on the phylogenetic 212 data available, we considered these two taxa as the same species and considered their 213 range maps jointly (Sheldon et al. 2009). Species occurrences in the native range were 214 then assigned to different phylogroups according to geographic boundaries defined by 215 Thiessen Polygons (Figs. 1a and S1 in Supp. Info.). Occurrence data were aggregated at 216 5-arcminute resolution, which corresponds approximately to 10 × 10 km. This 217 resolution was considered to be representative of the size of the smallest cities, and thus 218 adequate to capture the main responses of bird species to humanized environments with 219 acceptable computing time. Samples with reported geographical issues, location 220 uncertainty above 5 km, or with central grid coordinates of atlases of >10 km resolution, 221 were removed from analyses. Duplicate samples at the 5-arcminute resolution were 222 handled as single observations. Final sample sizes ranged from 3,156 to 32,052 for 223 different species in the native range, and from 2 to 2,368 in the alien range (Table 1). 224 225 Environmental variables 226 227 We considered eight bioclimatic variables (obtained from WorldClim 1, 228 http://www.worldclim.org/) (Hijmans et al. 2005), which are known to affect bird 229 distributions (Strubbe et al. 2015; Cardador et al. 2016): annual mean temperature, 230 temperature seasonality (standard deviation ×100), maximum temperature of the 231 warmest month, minimum temperature of the coldest month, annual precipitation, 232 precipitation of the driest month, precipitation of the wettest month and precipitation 233 seasonality (coefficient of variation). However, annual mean temperature, temperature 234 seasonality and annual precipitation were highly correlated with other climate variables 235 across the world (r ≥0.90) and thus removed from analyses (Cardador and Blackburn 236 2019). We considered two variables as descriptors of human transformed environments: 237 i) the Global Human Influence Index, which provides a weighed composite map of 238 anthropogenic impacts including urban extent, population density, land cover, night 239 lights and distance to roads, railways, navigable rivers and coastlines (Sanderson et al. 240 2002), and ii) the percentage of urban habitats, as a more specific descriptor of 241 urbanization. The percentage of urban habitats at the 5-arcminute resolution was derived 242 from MODIS-based global land cover climatology data at 500m resolution (Broxton et 243 al. 2014). 244 245 Niche analyses 246 247 We compared the climatic and human disturbance niches of different phylogroups of a 248 given species in its native range using the framework proposed by Broennimann et al. 249 (2012). These analyses involved four steps: (1) definition of a two-dimensional gridded 5 250 environmental space, (2) calculation of the density of occurrences along the two- 251 dimensional environmental space using kernel smoothers (3) measurement of niche 252 overlap between occurrence densities of two given groups along the environmental 253 space and (4) statistical tests of niche similarity using a randomization approach. For 254 climatic niche analyses the gridded environmental space was that formed by the first 255 two axes of a PCA on the five climatic variables considered (see above) in 5-arcminute 256 pixels across the world. These axes explained 78% of the inertia. The first PCA axis 257 (45%) predominantly represented temperature gradients (with higher values 258 representing warmer climates) while the second axis (33%) represented precipitation 259 gradients (with higher values representing drier climates with higher precipitation 260 seasonality, Table S3). For the human space, the two axes represented the two human- 261 related variables considered (i.e., the Global Human Influence Index and the percentage 262 of urban environments - see above). Niche similarity was assessed using the Schoener's 263 D metric, calculated from the occupancies in the environmental space depicted by the 264 two first axes of the PCA. This metric indicates the overall match between two niches 265 over the whole climatic or human spaces and ranges from 0 (no overlap) to 1 (complete 266 overlap). We calculated niche similarity between each phylogroup of a species and all 267 the other phylogroups using one-to-one comparisons. 268 We assessed niche conservatism by conducting niche similarity tests, whereby 269 each obtained value of niche similarity was compared against a null distribution of 100 270 simulated similarity values (obtained when comparing the observed niche of one 271 phylogroup with niches obtained by drawing occurrences at random within available 272 habitats and vice versa) (Warren et al. 2008; Broennimann et al. 2012). Note that the 273 niche similarity test is thus bidirectional, and two tests were conducted for each 274 comparison between two phylogroups. As available habitat for each phylogroup in 275 species’ native ranges for analyses, we considered two alternative approaches. First we 276 considered all ecoregions (Olson et al. 2001) occupied by each species in its native 277 range (Figs. 1b and S2), as this might represent the complete gradient of climatic and 278 human conditions that the study species could have reasonably encountered considering 279 that dispersal is expected to be mainly limited by major biogeographical barriers in 280 native ranges (Soberon and Peterson 2005). Second, we accounted for potential spatial 281 constraints limiting the access of different phylogroups to available habitat for the 282 whole species. For this, we repeated the niche similarity tests considering as available 283 habitat for each phylogroup that present within the geographic range actually occupied 284 by that particular phylogroup: i.e., that limited by the geographic boundaries derived 285 from Thiessen polygons (Fig.S1). The results of both analyses were highly consistent 286 (Table S4), and for simplicity we thus only provide those for the former approach in the 287 main manuscript. All analyses were conducted using the ‘ecospat’ library in R software 288 (Broenniman et al. 2014). The minimum sample size for analyses was five occurrences 289 (Broennimann et al. 2012). 290 291 Species distribution models 292 293 We fitted SDMs calibrated on occurrences of different phylogroups (single-phylogroup 294 models) and for the species as a whole (species-level model) to generate global 295 predictions of species potential distributions outside their native ranges. For robustness 296 of analyses only phylogroups with more than 50 occurrence locations (exceptions: 1 297 phylogroup for the northern cardinal Cardinalis cardinalis, 1 for the oriental magpie- 298 robin Copsychus saularis and 9 for the ring-necked parakeet Psittacula krameri) were 299 retained for single-phylogroup models (Stockwell and Peterson 2002). Data from 6 300 phylogroups with fewer than 50 occurrence locations were also disregarded from 301 species-level models, to avoid biases in model comparisons related to differences in the 302 data used. Predictions were derived from an ensemble model of three techniques – 303 generalized linear models, MAXENT and random forest – using R library ‘biomod2’. 304 Both the linear and quadratic terms of the climate and human predictors were 305 considered, to account for positive or negative responses to intermediate values of the 306 variables. All models were run with a single set of a maximum of 10,000 pseudo- 307 absences randomly drawn from all ecoregions occupied by each species across its native 308 range as with niche analyses (Figs. 1b and S2). Presences and pseudoabsences were 309 weighted as such to ensure neutral (0.5) prevalence. 310 To reduce the potential effect of sampling biases in the data, a bias file was 311 created by retrieving from GBIF occurrence data at the family level for each species 312 (Elith et al. 2010; Cardador and Blackburn 2019). We derived a kernel density map of 313 sampling bias at a 5-arcminute resolution using ArcMap 10.5 to be included as a fixed 314 effect in model training. Occurrence data from species in the same taxonomic family are 315 expected to suffer from the same detection limitations, reducing the effect of sampling 316 biases in observed distribution patterns. To further account for potential effects of the 317 data selection, we conducted 10 replicates for each model by using random samples 318 (70%) of the complete datasets. Final ensemble model predictions for each species and 319 phylogroup were generated as averaged means of all model replicates conducted (Fig. 320 1c-e). Sampling bias was set to its maximum value for model predictions. For each 321 species, a composite model prediction integrating information on all single-phylogroup 322 models was then developed (phylogroup-composite model) (Fig. 1f). For this, single- 323 phylogroup predictions were first standardized to a maximum value of 1 to make them 324 comparable. We then calculated the mean probability of occurrence of at least one of 325 the related phylogroups using the multiplicative probability method described in 326 Pearman et al. (2010). We converted continuous model predictions into binary 327 presence–absence maps by implementing a threshold for species presence that 328 maximized sensitivity plus specificity (Liu et al. 2005) in the training region (Fig. 1g). 329 In the case of composite models, binary maps were obtained by assigning species 330 presence to each cell that was predicted suitable for at least one single-phylogroup 331 model (Fig. 1h). To reduce problems related to model extrapolation, model projections 332 were adjusted using multivariate environmental similarity surfaces (MESSs) (Mateo et 333 al. 2014) (Fig. S3). Environmental suitability in dissimilar areas (MESS <0) was 334 considered to be zero. However, analyses using non-adjusted model projections were 335 highly concordant (see results). 336 Model accuracy of phylogroup-composite and species-level models in predicting 337 species occurrences in the native range were evaluated using the Boyce index (which 338 ranges from -1 to 1, with higher values indicating higher match (Hirzel et al. 2006)) and 339 AUC (which ranges from 0 to 1, with values up to 0.5 representing models not better 340 than random (Phillips et al. 2006)), using the libraries ‘pROC’ and ‘ecospat’ in R. 341 Sensitivity (i.e., the proportion of correctly classified presences) was also computed 342 using the binary maps derived from continuous predictions. Potential differences in 343 accuracy metrics among species-level and phylogroup-composite models in native 344 ranges were evaluated using non-parametric paired Wilcoxon signed rank tests. 345 346 Geographic extent of SDM predictions 347 348 We calculated the extent of predicted suitable habitat in km2 for each species at a global 349 scale according to binary map projections of single-phylogroup, phylogroup-composite 7 350 and species-level models, using the library ‘raster’ from R. Accounting for intraspecific 351 variation is expected to improve the representation of different phylogroups in models 352 (particularly the most scarce and narrowly distributed in the native range), and thus to 353 increase the geographical distribution of habitat conditions predicted as suitable for 354 alien birds: environments occupied by rarer phylogroups are expected to be marginal, 355 and thus result in low probabilities of occurrence in those environments when 356 distribution is modelled at the species level. Differences in the extent of phylogroup- 357 composite and species-level model predictions in alien ranges were evaluated using 358 one-tailed paired t-tests. We also tested the hypothesis that the niches occupied by the 359 most narrowly distributed phylogroups in the native range are little represented in 360 species-level models. For this, we assessed the relationship between the percentage of 361 habitat predicted as suitable by single-phylogroup models also predicted as suitable by 362 species-level models and phylogroup relative range size (i.e., the proportion of species 363 native range covered by a given phylogroup). We used general linear models (GLM) for 364 those analyses. 365 366 Accuracy of SDM predictions 367 368 SDM-based predictions were tested against real occurrence data in alien ranges using 369 the Boyce index, AUC and sensitivity using the libraries ‘pROC’ and ‘ecospat’ in R as 370 in native ranges. As background for AUC and Boyce calculations, we followed the 371 framework proposed by (Strubbe et al. 2013, 2015) to estimate the area that could have 372 been effectively accessible to introduced birds. We buffered each 5-arcminute alien 373 locality with a distance equal to the minimum invasion speed recorded for birds (i.e. 374 4.59 km year, derived from (Blackburn et al. 2009)) multiplied by the number of years 375 since introduction (Fig. S4). For localities with duplicate records the oldest year was 376 used in analyses. When the specific year of introduction was not provided (≤ 4% of total 377 alien localities for each species), we were conservative and only considered that locality 378 (not a buffer around it) for background calculations. Buffers were only allowed to cover 379 areas effectively outside species native breeding ranges. Both occurrence localities in 380 regions where the species has effectively established (see ‘Phylogenetic and occurrence 381 data’ section) and known introduction localities (compiled from (Redding et al. 2019)) 382 were used. Differences in the accuracy of phylogroup-composite and species-level 383 model predictions in alien ranges were evaluated using one-tailed paired t-tests (AUC) 384 and one-tailed paired Wilcoxon signed rank tests (Boyce and sensitivity) according to 385 fit of normality and homoscedasticity assumptions. 386 387 RESULTS 388 389 Intraspecific climatic and human niche variation 390 391 Phylogenetic groups within species occupied partially overlapping portions of the 392 climate and human spaces available in the native range. However, climatic niche 393 similarity among phylogroups of a species was low (mean ± SD, D = 0.07 ± 0.09, N = 394 8) and not more similar than expected by chance for the vast majority of phylogroups 395 within species (95% of 101 reciprocal similarity tests involving 32 phylogroups from 396 eight species had P > 0.05, Table S4). Human niche similarity was higher than climatic 397 niche similarity (D = 0.32 ± 0.13), but generally not more similar (97% of tests had P > 398 0.05) than expected by chance for most comparisons. 8 399 400 Species distribution models in native ranges 401 402 Species-level and phylogroup-composite model predictions showed a good agreement 403 with species occurrences in the native ranges (mean ± SD, species-level: Boyce = 0.88 ± 404 0.24, AUC = 0.86 ± 0.11, sensitivity = 0.78 ± 0.19; phylogroup-composite: Boyce = 405 0.91 ± 0.12, AUC = 0.84 ± 0.09, sensitivity = 0.89 ± 0.10). No significant differences in 406 model performance between species-level and phylogroup-composite models were 407 observed (one-tailed paired Wilcoxon rank test: Boyce, V = 21.5, P = 0.71; AUC, V = 408 28, P = 0.93; sensitivity, V = 6, P = 0.054). 409 410 Geographic extent and accuracy of model projections in alien ranges 411 412 Species-level and phylogroup-composite models provided similar but not equal 413 predictions of habitat susceptible to invasion (Figs. 1e-h, S5 and S6). Notably, 414 phylogroup-composite models projected significantly larger habitat suitability at a 415 global scale than species-level models (Fig. 2a, t = -2.0, df = 7, P = 0.04). On average, 416 only 30 ± 16% of total pixels predicted as suitable by binary maps derived from both 417 species-level and phylogroup-composite models were coincident between both types of 418 models, while 52 ± 30% of pixels were predicted as suitable only by phylogroup- 419 composite models and 18 ± 19% only by species-level models. The capacity of species- 420 level models to predict the occurrence of individual phylogroups was significantly 421 related to phylogroup range size (estimate: 0.68 ± 0.13, P < 0.001, R2 = 0.55; Fig. 3). 422 When used to predict occurrences in alien ranges, model accuracy was on average 423 good according to different metrics considered (Figs. 2b-d), but high variability was 424 observed among species (Figs. 2b-d; see also Fig S7 for comparisons of models not 425 adjusted by MESS analyses). Accounting for intraspecific niche variation did not 426 improve model accuracy (Boyce, V = 21, P = 0.69; AUC, t = -0.09, df = 7, P= 0.47; 427 sensitivity, V = 9, P= 0.22). These result hold when omitting Platycercus elegans from 428 analyses (Boyce, V = 14, P = 0.53; AUC, t = -0.29, df = 6, P= 0.39; sensitivity, V = 9, 429 P= 0.22), for which sample size in the alien range was very low (Table 1). 430 431 DISCUSSION 432 433 There is a limited number of studies considering intraspecific niche variation in models 434 assessing environmental susceptibility to the colonization by alien birds, and they have 435 rarely compared modelled predictions against independent sets of occurrence data 436 (Peterson et al., 2019). Our results agree with recent evidence suggesting that 437 conspecific phylogenetic lineages of a species can differ in the climates and human 438 disturbance conditions they experience in native ranges (Peterson and Holt 2003; 439 Pearman et al. 2010; D’Amen et al. 2013): we find little climatic niche similarity among 440 phylogroups of the 8 species in our study. Accounting for intraspecific variation in 441 SDMs modified the distribution and the extent of potential suitable habitat for the whole 442 species. However, contrary to our expectations, accounting for intraspecific niche 443 variation did not result in more accurate model predictions, according to current 444 distributions of established alien species. 445 The increased extent of phylogroup-composite model predictions, when compared 446 to classical species-level models, suggests that the omission of intraspecific niche 447 structure from species distribution models underestimates intraspecific realized niche 448 variation, and thus species-level prediction of habitat susceptible to invasion (Peterson 9 449 and Holt 2003; D’Amen et al. 2013; Godefroid et al. 2016). The extents of geographic 450 areas occupied by different phylogroups of a species in our study were not equivalent 451 (Fig. S1), and this may lead to some phylogroups (particularly the most scarce and 452 narrowly distributed in the native range) having little representation in model 453 predictions obtained by classical species level models (Fig. 3). 454 However, contrary to our expectations, accounting for intraspecific niche variation 455 did not improve model accuracy when predicting occurrences in alien ranges. These 456 results could be related to different, not mutually exclusive processes. First, it should be 457 noted that phylogroup separation – in this study and in general – relies on neutral 458 molecular markers. Strong structuring in these markers across populations indicates 459 reduced dispersal and thus increased potential for local genetic adaptation to emerge 460 (Lenormand 2002). Yet, we do not know whether genetic structuring in neutral markers 461 really reflects local genetic adaptation to climatic and human environmental conditions 462 experienced (Holderegger et al. 2006), and thus in these species tolerances, or just 463 reflects differences in the realized (i.e., occupied) niche of different phylogroups 464 (Guisan et al. 2014; Peterson et al. 2019). In this sense, since most phylogroups within 465 species considered have allopatric distributions, observed realized niche divergence may 466 have been driven by different environmental conditions in the range of each 467 phylogroup, rather than by adaptation to different conditions within a shared spatial 468 distribution (Maia-Carvalho et al. 2018). If so, dividing a species’ range into several 469 groups might have yielded different climate–occupancy relationships for each group 470 even in the absence of local adaptation. In our study, this could explain the larger 471 species range predictions of models considering phylogroup information, but the quite 472 similar prediction performances. For some species it is also possible that local 473 adaptation occurs at higher or lower taxonomic levels (Peterson et al. 2019; Smith et al. 474 2019) or geographic scales (Cardador et al. 2016). 475 Second, lineages little represented in classical species-level models are also those 476 less likely to be introduced and subsequently established in new areas, given the 477 positive effects of abundance and geographic range size on introduction and 478 establishment success in alien species (Blackburn and Duncan 2001b, a). If this is the 479 case, accounting for intraspecific niche variation would result in small differences in 480 model prediction for the most commonly translocated taxa. At the same time, 481 accounting for intraspecific niche variation might overestimate the potential alien range 482 size, as the environmental tolerances of alien individuals would actually be narrower 483 than that of the species as a whole. Under this hypothesis (i.e., the more common 484 phylogroups in native areas are more often introduced and established in alien ranges), 485 higher model accuracy for the more common phylogroups should be expected. 486 However, post-hoc analyses assessing the relationship between accuracy of single- 487 phylogroup model predictions and phylogroup relative range size (i.e., the proportion of 488 species native range covered by a given phylogroup) offer little support for this 489 hypothesis (Pearson correlation coefficients between accuracy metrics and phylogroup 490 range size, AUC: r = -0.20, P =0.33; Boyce: r = -0.01, P = 0.97, N = 25). 491 Third, as most bird introductions are relatively recent, species might not occupy 492 all of the potential suitable environments available in the invaded range, due to dispersal 493 limitations (Blackburn et al. 2009; Ascensão et al. 2020). For some species, current 494 alien distributions may thus reflect the characteristics of the new introduction localities 495 rather than optimum environmental conditions, which might be more likely to be moved 496 into during spread (Abellán et al. 2017). Furthermore, while climate and human 497 variables appeared to be major factors shaping alien species distributions (Cardador and 498 Blackburn 2019), omission of other important drivers of bird distributions, such as 10 499 interspecific interactions (Blackburn et al. 2009; Redding et al. 2019), might also 500 produce an overestimate of the projected suitable area for species. Some species might 501 even undergo niche shifts in alien ranges, although this seems to be less frequent 502 (Broennimann et al. 2007; Strubbe et al. 2013; Cardador and Blackburn 2020). All of 503 these issues may have contributed to reduce model accuracy for both species-level and 504 phylogroup-composite model predictions in alien ranges, masking potential differences. 505 Overall, our results show variability in model predictions linked to taxonomic 506 level considered. However, comparisons of model predictions with current available 507 distribution data in alien ranges do not provide evidence of an improvement in 508 prediction accuracy for models accounting for intraspecific niche variation. This result 509 is relevant given the lack of information about the presence and geographical 510 distribution of phylogeographic lineages for many species: according to our results, 511 until more accurate information on intraspecific variability is available, species-level 512 models can be reasonable candidates. However, when phylogeographic information is 513 available, the use of the most conservative criterion (i.e. to model both species and 514 lineages on the basis of the actual range, e.g. Mori et al. 2019) is recommended, given 515 that the ultimate purpose of such modelling exercises is to reduce invasion risks and 516 their consequences on biodiversity conservation. 517 518 519 DECLARATIONS 520 521 Funding 522 This project has received funding from the European Union’s Horizon 2020 research and 523 innovation programme under the Marie Skłodowska-Curie Grant Agreement No. 752149. 524 L.C. also received funding from the Beatriu de Pinós fellowship program (funded by the 525 Catalan Government and EU COFUND program nº 801370). PA is funded by ‘V Plan 526 Propio de Investigación’ of the Universidad de Sevilla (Spain). 527 528 Conflicts of interest/Competing interests 529 The authors have no conflict of interest to disclose. 530 531 Availability of data and material 532 533 All data used in this manuscript have already been published or archived and can be 534 accessible via Figshare http://dx.doi.org/10.6084/m9.figshare.4234850 (Dyer et al. 535 2016), GBIF.org (2018 - for a complete list of GBIF occurrence downloads see Table 536 S1), WorldClim database (Hijmans et al. 2005) available for download from 537 http://www.worldclim.org, BirdLife International & NatureServe, (2014) available for 538 download from http://datazone.birdlife.org/species/requestdis and terrestrial 539 ecoregions of the world (Olson et al. 2001) available for download from 540 http://www.worldwildlife.org/publications/terrestrial-ecoregions-of-the-world. Data on 541 geographic extent and accuracy of model projections will be accessible via Figshare 542 should the manuscript be accepted. 543 544 545 Code availability 546 547 Code associated with this manuscript will be accessible via Figshare should the manuscript 548 be accepted. 549 11 550 Authors' contributions 551 ‘Not applicable’ 552 553 Ethics approval 554 ‘Not applicable’ 555 556 Consent to participate 557 ‘Not applicable’ 558 559 Consent for publication 544 560 ‘Not applicable’ 561 562 REFERENCES 563 564 Abellán P, Carrete M, Anadón JD, et al (2016) Non-random patterns and temporal 565 trends (1912-2012) in the transport, introduction and establishment of exotic birds 566 in Spain and Portugal. Divers Distrib 22:263–273. 567 https://doi.org/10.1111/ddi.12403 568 Abellán P, Tella JL, Carrete M, et al (2017) Climate matching drives spread rate but not 569 establishment success in recent unintentional bird introductions. Proc Natl Acad 570 Sci 114:9385–9390. https://doi.org/10.1073/pnas.1704815114 571 Araújo MB, Peterson AT (2012) Uses and misuses of bioclimatic envelope modeling. 572 Ecology 93:1527–1539. https://doi.org/10.1890/11-1930.1 573 Ascensão F, Latombe G, Anadón JD, et al (2020) Drivers of compositional dissimilarity 574 for native and alien birds: the relative roles of human activity and environmental 575 suitability. Biol Invasions. https://doi.org/10.1007/s10530-020-02196-7 576 BirdLife International, NatureServe (2014) Bird Species Distribution Maps of the 577 World 578 Blackburn TM, Duncan RP (2001a) Determinants of establishment success in 579 introduced birds. Nature 414:195–197. https://doi.org/10.1038/35102557 580 Blackburn TM, Duncan RP (2001b) Establishment patterns of exotic birds are 581 constrained by non-random patterns in introduction. 927–939 582 Blackburn TM, Lockwood JL, Cassey P (2009) Avian invasions. The ecology and 583 evolution of exotic birds. Oxford University Press, Oxford 584 Broenniman O, Petitpierre B, Randin C, et al (2014) Package ‘ ecospat ’ 585 Broennimann O, Fitzpatrick MC, Pearman PB, et al (2012) Measuring ecological niche 586 overlap from occurrence and spatial environmental data. Glob Ecol Biogeogr 587 21:481–497. https://doi.org/10.1111/j.1466-8238.2011.00698.x 588 Broennimann O, Treier UA, Müller-Schärer H, et al (2007) Evidence of climatic niche 589 shift during biological invasion. Ecol Lett 10:701–709. 590 https://doi.org/10.1111/j.1461-0248.2007.01060.x 591 Broxton PD, Zeng X, Sulla-Menashe D, Troch PA (2014) A global land cover 592 climatology using MODIS data. J Appl Meteorol Climatol 53:1593–1605. 593 https://doi.org/10.1175/JAMC-D-13-0270.1 594 Capinha C, Essl F, Seebens H, et al (2015) The dispersal of alien species redefines 595 biogeography in the Anthropocene. Science (80- ) 348:1248–1251. 596 https://doi.org/10.1126/science.aaa8913 597 Cardador L, Blackburn TM (2020) A global assessment of human influence on niche 598 shifts and risk predictions of bird invasions. Glob Ecol Biogeogr 599 Cardador L, Blackburn TM (2019) Human-habitat associations in the native 600 distributions of alien bird species. J Appl Ecol 56:1189–1199. 12 601 https://doi.org/10.1111/1365-2664.13351 602 Cardador L, Carrete M, Gallardo B, Tella JL (2016) Combining trade data and niche 603 modelling improves predictions of the origin and distribution of non-native 604 European populations of a globally invasive species. J Biogeogr 43:967–978. 605 https://doi.org/10.1111/jbi.12694 606 Cardador L, Lattuada M, Strubbe D, et al (2017) Regional Bans on Wild-Bird Trade 607 Modify Invasion Risks at a Global Scale. Conserv Lett 10:717–725. 608 https://doi.org/10.1111/conl.12361 609 Cardador L, Tella J, Anadón J, et al (2019) The European trade ban on wild birds 610 reduced invasion risks. Conserv Lett 12:e12631 611 D’Amen M, Zimmermann NE, Pearman PB (2013) Conservation of phylogeographic 612 lineages under climate change. Glob Ecol Biogeogr 22:93–104. 613 https://doi.org/10.1111/j.1466-8238.2012.00774.x 614 Dyer EE, Redding DW, Blackburn TM (2016) Data from: The global avian invasions 615 atlas, a database of alien bird distributions worldwide. Figshare, 616 https://doi.org/106084/m9.figshare4234850 617 Dyer EE, Redding DW, Blackburn TM (2017) The global avian invasions atlas, a 618 database of alien bird distributions worldwide. Sci Data 4:170041. 619 https://doi.org/10.1038/sdata.2017.41 620 Elith J, Kearney M, Phillips S (2010) The art of modelling range-shifting species. 621 Methods Ecol Evol 1:330–342. https://doi.org/10.1111/j.2041-210X.2010.00036.x 622 Gaston KJ (2003) The structure and dynamics of geographic ranges. Oxford University 623 Press, New York 624 GBIF.org (2017) GBIF Home Page 625 Godefroid M, Rasplus J-Y, Rossi J-P (2016) Is phylogeography helpful for invasive 626 species risk assessment? The case study of the bark beetle genus Dendroctonus. 627 Ecography (Cop) 39:1197–1209. https://doi.org/10.1111/ecog.01474 628 Guisan A, Petitpierre B, Broennimann O, et al (2014) Unifying niche shift studies: 629 Insights from biological invasions. Trends Ecol Evol 29:260–269. 630 https://doi.org/10.1016/j.tree.2014.02.009 631 Guisan A, Thuiller W (2005) Predicting species distributions: offering more than simple 632 habitat models. Ecol Lett 8:993–1009 633 Haring E, Gamauf A, Kryukov A (2007) Phylogeographic patterns in widespread corvid 634 birds. Mol Phylogenet Evol 45:840–862. 635 https://doi.org/10.1016/j.ympev.2007.06.016 636 Hijmans RJ, Cameron SE, Parra JL, et al (2005) Very high resolution interpolated 637 climate surfaces for global land areas. Int J Climatol 25:1965–1978 638 Hirzel AH, Gwenaëlle LL, Helfer V, et al (2006) Evaluating the ability of habitat 639 suitability models to predict species presences. Ecol Modell 1999:142–152. 640 https://doi.org/https://doi.org/10.1016/j.ecolmodel.2006.05.017 641 Holderegger R, Kamm U, Gugerli F (2006) Adaptive vs. neutral genetic diversity: 642 implications for landscape genetics. Landsc Ecol 21:797–807. 643 https://doi.org/10.1007/s10980-005-5245-9 644 Hulme PE (2009) Trade, transport and trouble: Managing invasive species pathways in 645 an era of globalization. J Appl Ecol 46:10–18. https://doi.org/10.1111/j.1365- 646 2664.2008.01600.x 647 Hulme PE, Bacher S, Kenis M, et al (2008) Grasping at the routes of biological 648 invasions: A framework for integrating pathways into policy. J Appl Ecol 45:403– 649 414. https://doi.org/10.1111/j.1365-2664.2007.01442.x 650 Huntley B, Green RE, Collingham YC, Willis SG (2007) A climatic atlas of European 13 651 breeding birds. Durham University, The RSPB and Lynx Edicions, Barcelona 652 Joseph L, Dolman G, Donnellan S, et al (2008) Where and when does a ring start and 653 end? Testing the ring-species hypothesis in a species complex of Australian 654 parrots. Proc Biol Sci 275:2431–40. https://doi.org/10.1098/rspb.2008.0765 655 Lecocq T, Harpke A, Rasmont P, Schweiger O (2019) Integrating intraspecific 656 differentiation in species distribution models: Consequences on projections of 657 current and future climatically suitable areas of species. Divers Distrib 25:1088– 658 1100. https://doi.org/10.1111/ddi.12916 659 Lenormand T (2002) Gene flow and the limits to natural selection. Trends Ecol Evol 660 17:183–189. https://doi.org/10.1016/S0169-5347(02)02497-7 661 Liu C, Berry PM, Dawson TP, Pearson RG (2005) Selecting thresholds of occurrence in 662 the prediction of species distributions. Ecography (Cop) 3:385–393 663 Liu C, Wolter C, Xian W, Jeschke JM (2020) Species distribution models have limited 664 spatial transferability for invasive species. Ecol Lett ele.13577. 665 https://doi.org/10.1111/ele.13577 666 Liukkonen-Anttila T, Uimaniemi L, Orell M, Lumme J (2002) Mitochondrial DNA 667 variation and the phylogeography of the grey partridge (Perdix perdix) in Europe: 668 From pleistocene history to present day populations. J Evol Biol 15:971–982. 669 https://doi.org/10.1046/j.1420-9101.2002.00460.x 670 Maia-Carvalho B, Vale CG, Sequeira F, et al (2018) The roles of allopatric 671 fragmentation and niche divergence in intraspecific lineage diversification in the 672 common midwife toad (Alytes obstetricans). J Biogeogr 45:2146–2158. 673 https://doi.org/10.1111/jbi.13405 674 Mateo RG, Broennimann O, Petitpierre B, et al (2014) What is the potential of spread in 675 invasive bryophytes? Ecography (Cop) n/a-n/a. https://doi.org/10.1111/ecog.01014 676 McKinney M, Lockwood J (1999) Biotic homogenization: a few winners replacing 677 many losers in the next mass extinction. Trends Ecol Evol 14:450–453 678 Mori E, Menchetti M, Zozzoli R, Milanesi P (2019) The importance of taxonomy in 679 species distribution models at a global scale: the case of an overlooked alien 680 squirrel facing taxonomic revision. J Zool 307:43–52. 681 https://doi.org/10.1111/jzo.12616 682 Olson DM, Dinerstein E, Wikramanayake ED, et al (2001) Terrestrial ecoregions of the 683 world: a new map of life on Earth. Bioscience 51:933–938 684 Parmesan C, Ryrholm N, Stefanescu C, et al (1999) Poleward shifts in geographical 685 ranges of butterfly species associated with regional warming. Nature 399:579–583. 686 https://doi.org/doi.org/10.1038/21181 687 Pearman PB, Amen MD, Graham CH, et al (2010) Within-taxon niche structure : niche 688 conservatism , divergence and predicted effects of climate change. 689 https://doi.org/10.1111/j.1600-0587.2010.06443.x 690 Peterson AT (2011) Ecological niche conservatism: a time-structured review of 691 evidence. J Biogeogr 38:817–827. https://doi.org/10.1111/j.1365- 692 2699.2010.02456.x 693 Peterson AT, Holt RD (2003) Niche differentiation in Mexican birds: Using point 694 occurrences to detect ecological innovation. Ecol Lett 6:774–782. 695 https://doi.org/10.1046/j.1461-0248.2003.00502.x 696 Peterson ML, Doak DF, Morris WF (2019) Incorporating local adaptation into forecasts 697 of species ’ distribution and abundance under climate change. Glob Chang Biol 698 775–793. https://doi.org/10.1111/gcb.14562 699 Phillips SJ, Anderson RP, Schapire RE (2006) Maximum entropy modeling of species 700 geographic distributions. Ecol Modell 190:231–259. 14 701 https://doi.org/10.1016/j.ecolmodel.2005.03.026 702 Pimm SL, Russell GJ, Gittleman JL, Brooks TM (1995) The future of biodiversity. 703 Science (80- ) 269:347–350 704 Redding DW, Pigot AL, Dyer EE, et al (2019) Location-level processes drive the 705 establishment of alien bird populations worldwide. Nature 571:103–106. 706 https://doi.org/10.1038/s41586-019-1292-2 707 Sanderson EW, Malanding J, Levy MA, et al (2002) The human footprint and the last of 708 the wild. Bioscience 52:891–904. https://doi.org//10.1641/0006- 709 3568(2002)052[0891:THFATL]2.0.CO;2 710 Sayol F, Cooke RSC, Pigot AL, Blackburn TM, Tobias JA, Steinbauer MJ, Antonelli A, 711 Faurby, S. (2021). Loss of functional diversity through anthropogenic extinctions 712 of island birds is not offset by biotic invasions. Science Advances, in press. 713 Sheldon FH, Lohman DJ, Lim HC, et al (2009) Phylogeography of the magpie-robin 714 species complex (Aves: Turdidae: Copsychus) reveals a Philippine species, an 715 interesting isolating barrier and unusual dispersal patterns in the Indian Ocean and 716 Southeast Asia. J Biogeogr 36:1070–1083. https://doi.org/10.1111/j.1365- 717 2699.2009.02087.x 718 Smith AB, Godsoe W, Rodríguez-Sánchez F, et al (2019) Niche Estimation Above and 719 Below the Species Level. Trends Ecol Evol 34:260–273. 720 https://doi.org/10.1016/j.tree.2018.10.012 721 Smith BT, Escalante P, Hernández Bãos BE, et al (2011) The role of historical and 722 contemporary processes on phylogeographic structure and genetic diversity in the 723 Northern Cardinal, Cardinalis cardinalis. BMC Evol Biol 11:. 724 https://doi.org/10.1186/1471-2148-11-136 725 Soberón J (2007) Grinnellian and Eltonian niches and geographic distributions of 726 species. Ecol Lett 10:1115–1123 727 Soberon J, Peterson AT (2005) Interpretation of Models of Fundamental Ecological 728 Niches and Species’ Distributional Areas. Biodivers Informatics 2:1–10. 729 https://doi.org/10.17161/bi.v2i0.4 730 Stockwell DR. B, Peterson AT (2002) Effects of sample size on accuracy of species 731 distribution models. Ecol Modell 148:1–13. https://doi.org/10.1016/S0304- 732 3800(01)00388-X 733 Strubbe D, Broennimann O, Chiron F, Matthysen E (2013) Niche conservatism in non- 734 native birds in Europe: niche unfilling rather than niche expansion. Glob Ecol 735 Biogeogr 22:962–970. https://doi.org/10.1111/geb.12050 736 Strubbe D, Jackson H, Groombridge J, Matthysen E (2015) Invasion success of a global 737 avian invader is explained by within-taxon niche structure and association with 738 humans in the native range. Divers Distrib 21:675–685. 739 https://doi.org/10.1111/ddi.12325 740 Vitousek PM (1997) Human Domination of Earth’s Ecosystems. Science (80- ) 741 277:494–499. https://doi.org/10.1126/science.277.5325.494 742 Warren DL, Glor RE, Turelli M (2008) Environmental niche equivalency versus 743 conservatism: quantitative approaches to niche evolution. Evolution 62:2868–83. 744 https://doi.org/10.1111/j.1558-5646.2008.00482.x 745 Zink RM, Pavlova A, Drovetski S, Rohwer S (2008) Mitochondrial phylogeographies 746 of five widespread Eurasian bird species. J Ornithol 149:399–413. 747 https://doi.org/10.1007/s10336-008-0276-z 748 749 15 Table 1. Bird species selected on the basis of the existing literature on phylogeographic structure. Species name, number of phylogroups per species, source for phylogenetic data and the number of occurrences available in the native and alien ranges according to GBIF (www.gbif.org) are provided. For native occurrences, the minimum and maximum occurrences available for different phylogroups are shown in parentheses. Species Phylogroups References Native occurrences Alien occurrences Alauda arvensis 2 Zink et al. 2008 19789 (232- 19557) 2378 Cardinalis cardinalis 6 Smith et al. 2011 41372 (9- 40079) 200 Copsychus saularis 3 Sheldon et al. 2009 3893 (20 - 3726) 36 Corvus frugilegus 2 Haring et al. 2007 7445 (52- 7393) 53 Perdix perdix 2 Liukkonen-Anttila et al. 2002 6444 (1804 - 6291) 2316 Pica pica 2 Haring et al. 2007 32052 (782 - 31270) 115 Platycercus elegans 2 Joseph et al. 2008 5046 (79- 4967) 2 Psittacula krameri 17 Strubbe et al. 2015 3156 (1-1834) 701 750 16 751 Figure Legends 752 753 Figure 1. Phylogroup distribution, native occurrences and model predictions for the 754 Eurasian skylark Alauda arvensis. (a) Distribution of two different phylogroups (red 755 and orange) as derived from Thiessen polygons and locations with genetic data (points); 756 (b) native occurrence data (black dots) and background (i.e., all ecoregions occupied by 757 each species in its native range, grey) used for species distribution modelling; (c-d) 758 continuous predictions of single phylogroup models, polygons show the distribution of 759 the phylogroup used in models; (e) species-level model continuous predictions; (f) 760 phylogroup-composite model continuous predictions; (g) species-level model binary 761 predictions; (h) phylogroup-composite model binary predictions. In (g) and (h), green 762 represents habitat predicted as suitable and grey as non-suitable. 763 764 Figure 2. Comparisons of (a) extent of suitable habitat and (b – d) model accuracy in 765 alien ranges between species distribution models accounting or not for phylogroup 766 niche variation. Extent of suitable habitat (a) is derived from binary maps using the 767 maximum sensitivity plus specificity threshold. Model accuracy in alien ranges is 768 assessed by the AUC (b), TSS (c) and sensitivity (d). N = 8. 769 770 Figure 3. Similarity between single-phylogroup and species-level model predictions in 771 relation to phylogroup relative range size. Similarity refers to the percentage of pixels 772 predicted as suitable by a single-phylogroup model also predicted as suitable by 773 species-level models at a global scale. 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 17 801 Figure 1. 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 18 827 Figure 2. 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 19 850 Figure 3. 851 852 853 20