Spatial Sorting Drives Morphological Variation in the Invasive Bird, Acridotheris tristis

The speed of range expansion in many invasive species is often accelerating because individuals with stronger dispersal abilities are more likely to be found at the range front. This ‘spatial sorting’ of strong dispersers will drive the acceleration of range expansion. In this study, we test whether the process of spatial sorting is at work in an invasive bird population (Common myna, Acridotheris tristis) in South Africa. Specifically, we sampled individuals across its invasive range and compared morphometric measurements relevant and non-relevant to the dispersal ability. Besides testing for signals of spatial sorting, we further examined the effect of environmental factors on morphological variations. Our results showed that dispersal-relevant traits are significantly correlated with distance from the range core, with strong sexual dimorphism, indicative of sex-biased dispersal. Morphological variations were significant in wing and head traits of females, suggesting females as the primary dispersing sex. In contrast, traits not related to dispersal such as those associated with foraging showed no signs of spatial sorting but were significantly affected by environmental variables such as the vegetation and the intensity of urbanisation. When taken together, our results support the role of spatial sorting in facilitating the expansion of Common myna in South Africa despite its low propensity to disperse in the native range.


Introduction
Ecologists often assume that species distributions are at their equilibriums ( [1] but see [2]) and treat species' geographical ranges as though they are static features [3]. However, changing environmental conditions can profoundly affect many aspects of species' geographic ranges and thus lead to dynamic distribution patterns [4][5]. It is increasingly recognized that populations can undergo rapid range shifts [2][3][6][7][8][9] and therefore dispersal ability at the edge of a species' geographic range becomes an important survival trait, particularly in the face of the negative effects of climate change [10][11]. In this regard, alien invasive species render excellent study systems to elucidate the evolutionary dynamics of dispersal traits in novel environments as they are often introduced to areas far removed from their historical native ranges and experience rapid range expansion.
Dispersal ability and survival in novel environments are paramount in shaping distribution ranges [12][13] and, not surprisingly, the evolution of dispersal in expanding populations has received much attention recently. For example, by simulating an expanding (invasive) population, Travis & Dytham [14] showed that individuals with a higher propensity to disperse were favoured at expanding edges of the distribution, resulting in increased dispersal rates at the range front. Phillips et al. [15] extended this notion by demonstrating that the dispersal kernel becomes less kurtotic and more skewed for individuals at the invading front during range expansion (invasion) events.
One potential mechanism by which dispersal is enhanced during range expansion has been termed spatial sorting, also known as spatial assortment by dispersal ability [16][17]. Spatial sorting can affect dispersal-relevant morphological traits such as wing size [18][19][20], leg size [21] or foot size [22] so that morphological variants associated with increased dispersal ability increases along the direction of range expansion.
Similarly, spatial sorting may also occur for cognitive-relevant traits that increase survival in novel environments. For instance, large brain size is typically associated with enhanced cognitive ability and has been linked to the establishment success and invasiveness in birds (see e.g. [23][24]). This is also true for behavioural traits such as aggression. For example, Duckworth [25] showed that more aggressive western bluebird individuals are more successful at expanding margins as they are more likely to occupy new habitats.
Evidently, spatial sorting is most relevant to traits that affect dispersal and fitness in novel environments. Once established, forces other than spatial sorting, such as local adaptation, competition and sexual selection [26][27], may act more intensely on other traits such as foraging traits (e.g. bill shape) that are mainly driven by available food resources [28][29]. It is therefore expected that traits associated with foraging to be not strongly affected by spatial sorting but, rather, by the characteristics of local environments. Thus, while spatial sorting will enhance dispersalrelevant traits towards a species' range margin, variation in traits such as foraging abilities will be more affected by local environmental conditions. It is possible to determine which processes shape morphological shifts by examining trait variation across entire expanding ranges of a given species. To this end, we studied morphological variations in a highly invasive bird species, the Common myna (Acridotheris tristis), across its entire distribution range in South Africa. This species was accidentally introduced to Durban on the Indian Ocean coast of South Africa in 1902, followed by a second introduction (putatively independent) inland in Johannesburg (557 km west from Durban) in 1938 [30]. The Common myna has a relatively low dispersal rate compared to other bird species (dispersal distances ,16 km; [31][32]) and is currently continually expanding its southern African range [33]. Specifically, we aim to determine whether any relationships exist between traits that are related to dispersal or survival in novel environments and the distance from the range core (i.e. the site of introduction). It is also expected that traits that are non-relevant to dispersal or survival in novel environments should not display similar relationships.

Morphological measurements
Common myna (Acridotheres tristis) specimens were collected from 58 sampling sites spanning the current non-native range across South Africa (Fig. 1). The protocol for this study was approved by the animal care committees of Pretoria University (No EC010-10). All birds used in the study were invasive species shot by hunters. Birds were sampled from public sites. Collections were carried out prior to the breeding season (June-August 2010) to ensure that female weights were not biased due to reproductive function. Carcasses were frozen at 220uC until further use. Morphological measurements were captured by one investigator (CBS) for consistency. A total of 389 adult birds were weighed (0.001 g precision) and dissected to determine the sex (217 males and 172 females). A 0.1 mm-unit dial calliper was used for all morphological measurements unless otherwise stated. Seven measurements were taken namely (a) wing length (from the carpal joint to the tip of the longest primary), (b) tarsus length (from the tibio-tarsus joint to the distal end of the tarso-metatarsus), (c) bill length (from the bill tip to the skull at nasal-frontal hinge), (d) bill depth (at the proximate edge of the nostrils), (e) bill width (at the proximate edge of the nostrils), (f) head length (from the tip of bill to the back of head, minus bill length), and (g) tail length (from the uropygial gland to the tip of the longest rectrix using a 1 mm-unit flat ruler).
To correct for size, we performed a principal component analyses (PCA) from which PC1 and PC2 can be interpreted as the size and shape variables respectively [34]. All seven morphological measurements were then separately regressed against PC1 (i.e. size) as the sole explanatory variable in a general linear model. To avoid statistical dependency between the dependent variable and PC1, we ran separate PCAs with the dependent variable removed. Significant residuals from this general linear model would indicate a deviation from the expected model and therefore significant morphological variations. These residuals were then used in all subsequent analyses. In addition, we also calculated four morphological ratios, namely (a) bill ratio (BR: bill length-towidth), (b) head ratio (HR: head-to-body length), (c) tarsus ratio (TR: tarsus-to-body length) and (d) wing-to-tail ratio (WTR), as well as wing load (the residual from the regression of log 10 [wing length 2 ] by log 10 [body weight], [35]). In total, 14 morphological measurements were included in our analyses (Table 1; see also  Table S1 and Table S2).

Environmental variables
We compiled a set of 11 environmental variables, representing variation in climate, habitat and human impact, across all sampling sites (Table 2). Four bioclimatic variables namely mean winter precipitation, mean summer precipitation, mean summer maximum temperature and mean winter minimum temperature were obtained from the WorldClim database [36]. Habitat variables include altitude, the normalised difference vegetation index (NDVI, January mean from 1982 to 1999, [37]) and the core distance measured as distance from the sampling sites to the possible centre of introduction (Johannesburg and Durban). It is worth noting that distance from the sampling sites to Johannesburg is not inversely related to the distance from the same sampling sites to Durban. Instead, both reflect core distances of each introduction event, i.e. Johannesburg and Durban, and can be used as explanatory variables in a multivariate analysis. We performed a PCA on five transformed land covers variables, namely, cultivated land, degraded land, tree plantation, irrigation and urban build-up [37] calculated at the quarter degree resolution derived from 1994 and 1995 Landsat TM satellite images [38]. The two first axes of the PCA, named Urban 1 and 2, were used as variables indicative of the Human impact.

Statistical analyses
Multi-scale pattern analyses. To investigate the spatial structure in our data, we used a PCA in which independent variables were regressed with predictors of spatial connectivity. This approach has been termed multi-scale pattern analysis (MSPA) [39] and is implemented in the package adegenet from the R statistical environment [40]. Specifically, we used Moran's eigenvector maps (MEMs), derived from Moran's index for spatial autocorrelation as spatial predictors [41]. MEMs have been shown to be good spatial predictors ( [39] and reference therein). MEMs are named in a decreasing order, where MEM 1 describes the spatial connectivity at the broadest spatial scale, with subsequent MEMs, describing spatial connectivity at reducing scales and the last MEM at the finest local spatial scale. The determination coefficient, R 2 , obtained from the MSPA was used to quantify the strength of association between the independent variables and MEMs [41]. Signals of spatial structure at particular scales in morphological measurements and environmental variables were detected for each sex separately (i.e. a total of four MSPAs). We then applied a redundancy analysis [42] to the MSPA for explaining the morphological variables by environmental variables at all spatial scales. This enabled us to combine spatial predictors and environmental variables simultaneously to describe the variation in morphological traits. A schematic description of the full procedure is summarised in Figure S1.
Linear mixed model. While the previous analysis was aimed at highlighting correlations between morphological traits and environmental variables across spatial scales, it will not provide the statistical significance of these interactions. Therefore, we performed a linear mixed model (LMM) for each morphological trait with the sampling site as a random factor and environmental variables as explanatory variables. The best model was identified by the Akaike information criterion (AIC, [43]) and the significance of each covariate tested using a Monte Carlo approach retaining only the significant covariates. When more than one variable was retained in the model, the potential collinearity between explanatory variables was also checked. If two explanatory variables exhibited a notable correlation (R 2 .0.1), the model was re-run for each variable separately to check whether significance was retained for each variable separately. An explanatory variable was removed from the final model if no significant effect was found. If both variables showed significant effects when run separately, both were retained in the model. We used a sequential regression to prioritise each explanatory variable by its significance and thus reduced the collinearity between covariates [44]. The LMM was first run for the entire dataset (females and males combined) and then retained models were run for each sex separately. This allowed us to examine whether environmental determinants, identified in the final model, are sexspecific.

Morphometric MSPA
The eigenvalues of the multi-scale pattern analysis (MSPA) for the morphological traits indicated that the first two principal component axes explained 40% and 42% of the variation in females and males respectively. For Moran's eigenvector maps (MEMs), we focused only on those variables that accounted for more than 3% of the morphological variation (see Fig. S1). Although this value is somewhat arbitrary, it nonetheless represents the point after which there was a strong decrease in the percentage of variation explained. Consequently, we identified five MEMs for females (MEM 1, 2, 6, 29 and 37 ) representing broadscale structures and two MEMs for males (MEM 1 at the broadest spatial scale and MEM 209 at a local scale). Spatial connectivity among samples failed to explain the observed variation in morphological traits as R 2 values were lower than 0.1 in both sexes. For females, two morphological traits (bill width and tail length) and bill ratio (BR) displayed the strongest, albeit still weak,  spatial structures with about 13 to 17% of their variation explained by the two first principal components of the MSPA. For males, a maximum of only 11% of variation explained for two traits (tail length and WTR).

Environmental MSPA
The spatial structure was strong in environmental variables and the eigenvalues of the MSPA demonstrated that the first two principal component axes explained more than 80% of the variation. Environmental spatial structure was consistently detected at broader scales (MEM 1-4 ). The broadest scale MEM 1 clearly resonated with the distance from Johannesburg (R 2 .0.4). The second and third MEMs (i.e. MEM 2 and MEM 3 ) were found mainly related to environmental factors that describe habitat quality, such as NDVI, altitude, temperature and precipitation (0.2,R 2 ,0.6).

Spatial and environmental MSPA
The MSPA redundancy analysis regressing morphological traits against environmental variables revealed some clear structure at different spatial scales (Fig. 2). For females, the first-group spatial predictors (MEM 1 and MEM 4 ) were related to traits influencing dispersal and cognitive capacity (wing length, wing loadings, tail length and head size) and explained 34 to 53% of the morphological variation ( Fig. 2A). The second-group spatial predictors (MEM 2 and MEM 3 ) were largely orthogonal to the first-group spatial predictors and concerned mainly traits related to foraging activities (bill size, bill ratio, bill depth and bill width), with 38 to 47.5% of the morphological variation explained ( Fig. 2A). Notably, the spatial predictors MEM 1 and MEM 4 were also associated with the distance from Johannesburg whilst the spatial predictors (MEM 2 and MEM 3 ) were related to the distance from Durban and other environmental factors of habitat quality. For males, the spatial predictor MEM 1 explained 58% of the variation in tail length and 74% of the variation in wing-tail ratio (WTR), whilst the second-group spatial predictors (MEM 2 and MEM 3 ) were weakly linked to nearly all the other traits, explaining from 7% of the variation in tail length to 45% of the variation in wing length, with no preferential traits highlighted (Fig. 2B).

Linear regression
When both sexes were included, the LMM analysis revealed significant sexual dimorphism in the Common myna, with females being smaller for all traits except for tail length, shape (PC2), bill ratio (BR) and bill width (Table 3). Sex explained 70 to 97% of the variation in morphological traits. Dispersal traits (wing length, wing loading and tail length) increased with altitude and the distance to Johannesburg, whereas birds collected further away from Johannesburg had bigger heads.
When the LMM was run for males and females separately, some of the environmental factors that have significant effects on morphological traits were found to be sex-specific (Table 3). For females, wing length and wing loading were influenced by altitude and the distance to Johannesburg, whereas only altitude had a significant effect on these traits in males. Similarly, only female head size increased with the distance to Johannesburg, and only the tail of male birds increased with both altitude and the distance to Johannesburg. Female tarsus was found shorter at the high altitude, whereas male tarsus dwindled with both altitude and winter precipitation. The environmental variable Urban 2, corresponding to degraded areas and commercial tree plantations, had a significant impact on male birds to have elongated and flattened bills (i.e. increasing bill length but reducing bill depth), whereas the bill width in male birds increased with the maximum temperature. The vegetation index NDVI was found mainly affecting head and bill ratios for both sexes. Value distributions of traits that have been subjected to spatial sorting are expected to shift between individuals at the core and at the edge of the distribution. We plotted value distributions for two traits that were only influenced by one factor to ensure the observation of only one effect. Head size in females (cognitive traits influence by the distance to Johannesburg) showed a clear shift due to spatial sorting between core and marginal individuals (Fig. 3A), whereas bill width in females (foraging traits) showed no obvious shift and thus a weak effect of spatial sorting (Fig. 3B). Discussion Spatial sorting theory predicts that dispersal ability should be enhanced at the range front of an expanding population into novel habitats. In other words, a significant positive correlation between dispersal ability and distance from the core of the distribution should be observed. This was supported by our results where a significant correlation was found between one of the core distances (Johannesburg) and traits related to dispersal ability (i.e. flying). The spatial assortment for stronger dispersers at the expanding edge creates a positive feedback that can accelerate the speed of range expansions. Such accelerating range expansion has been commonly identified in many highly invasive species, such as bush crickets [45] and cane toads [21]. Evidently, ever-increasing dispersal ability at the range margin, manifested through dispersal and cognitive traits, could play an important role in range dynamics and invasion success [21].
If sex-biased dispersal occurs in a particular species, one would expect a stronger effect from spatial sorting on one of the sexes. This is indeed the case for Common mynas in South Africa. Dispersal-relevant traits in females, specifically wing measurements (and also head size), showed a strong correlation with the distance to one of the introduction sites (Johannesburg) while no such a correlation was detected for males. Although no studies on the dispersal of Common mynas are available, the species does show signs of female-biased dispersal, consistent with Greenwood's [46][47] proposal that female-biased dispersal and male-biased philopatry are common in birds. Specifically, dispersal strategies are often linked to mating systems, resulting in resource defence in monogamy (typical in birds) where males take the lead role in acquisition and defence of resources and receive considerable benefits by remaining philopatric. As males are often more subjected to predation and aggression [48][49], aggression-related traits (such as tail length) may be more affected in males and therefore could explain the high morphological variation in tails for male Common mynas. In addition, sex biased dispersal could potentially lead to different sex ratios in expanding populations (e.g. a slightly lower, yet insignificant, sex ratio ( = 0.45) was found for birds within 250 km radius to Johannesburg than beyond (sex ratio = 0.49)).
An increase in cognitive ability would strongly favour individuals by allowing them to better cope under novel environmental conditions [23][24]25,50]. Head size is a qualitative proxy for brain size [24] and, subsequently, cognitive abilities [23]. For instance, increased head size has been linked to increased offspring defence and higher success at expanding range edges [25]. Head size in female Common mynas, the dispersing sex, increased significantly with distance away from one of the points of introduction (Johannesburg).
Why a significant correlation was found between many morphological traits and the distance to Johannesburg but not to Durban remains unknown. Two possible explanations warrant further investigation. First, dispersal-related traits often become homogenized once the range expansion stops [15] so that while the spatial sorting influences morphological variation in expanding populations, its effect will be diluted once populations reach their equilibriums. Specifically, when populations no longer expand, the backward migration of individuals from the margins to the core will homogenize morphological variations across the distribution range. Since the introduction to Durban predates the introduction to Johannesburg by nearly thirty years [33], it may be that the Durban expansion has potentially filled up most suitable habitats and reached the distributional equilibrium. Second, distinct environmental characteristics of these two introduction points could have differentially influenced populations' expansion. Johannesburg is located within the grassland biome of South Africa, bordering the savannah biome [51], whereas Durban is located within a subtropical thicket that extends along the east coast of the country. While the open grassland savannah may be more conducive to dispersal, the thicket and coastal forests surrounding Durban but also the Drakensberg mountain ridge seems impenetrable and may have contributed to prevent high levels of dispersal from this coastal introduction point.
Spatial sorting favours dispersal-relevant traits and those that provide a benefit at the range front during range expansion. Traits that do not provide any benefits in novel environments should therefore not show strong spatial sorting effect. Here we found that foraging traits in Common myna show no significant correlation with distance from the core. Conversely, factors of habitat quality could affect these non-dispersal-related traits. This is exactly what we found between the morphological variation in foragingrelevant traits and environmental variables (e.g. NDVI and urbanization). Specifically, urbanization can modify the quality and type of food resources and therefore influence bill shape (bill length and depth), consistent with previous studies [28][29]. As an indicator of primary productivity (and thus the habitat quality and food resources), NDVI was found to significantly influence the head ratio (HR) and bill ratio (BR) in both sexes. Besides food Table 3. Summary of environmental variables with significant effect on traits measured using linear regression models.  (2) Wing length -+|m(+)|f(+) +|f(+) Wing loading -+|m(+)|f(+) +|f(+) (2) TR -2|m (2) 2|m (2) BR -2|m (2) resources, other environmental factors are also expected to influence morphology and we found winter precipitation to negatively influence tarsus length. Similarly altitude appears to affect tarsus length negatively but wing and tail length positively. The latter findings are consistent with the evidence from Redbilled chough (Pyrrhocorax phyrrhocorax) and Alpine chough (Pyrrhocorax graculus) [52]. Spatial sorting has been one of the main hypotheses put forward to explain the increase in dispersal abilities in alien invasive species during range expansion. Spatial sorting allows dispersal-enhancing traits to accumulate at the expanding range front in only a few generations through phenotypic plasticity. Moreover, since dispersal can affect the processes of intraspecific competition and kin (inbreeding) avoidance, it is also under classical natural selection [14]. High dispersal rates can enable individuals to reduce the intensity of kin competition [53] and therefore could have an evolutionary advantage at the range front [54][55]. However, high dispersal rates can also act to constrain adaptation at the periphery and slow down range expansion [56,57]. The potential trade-off between dispersal and reproduction can further complicate any simple deductions from classical natural selection [58]. Therefore, the roles of spatial sorting (which does not rely upon fitness differentials) vs. natural selection (which assumes fitness differentials) during range expansion, while not mutually exclusive, remains a matter of debate [17,[59][60]. Spatial sorting generally acts at a faster pace than natural selection and thus suffices to enhance dispersal traits at the range front of a younger population; however, genetic analyses and common garden experiments are needed to quantify the relative contribution from each of these processes during range expansion.
In conclusion, we found that Common mynas in South Africa seemed to be subjected to spatial sorting and consequently dispersal-relevant traits showed a shift towards higher dispersal ability at the range of the species' current distribution (e.g., see [17]). Our findings confirm that this invasive species has enhanced its dispersal ability during the range expansion under novel environmental conditions despite being a relatively poor disperser in its native area. This has important implications for invasion biology in general as life-history traits from a species' native range are often used to inform ecologists about the potential to disperse and spread in non-native ranges [27].