The spatial analysis of biological interactions: morphological variation responding to the co‐occurrence of competitors and resources

By sharing geographic space, species are forced to interact with one another and the contribution of this process to evolutionary and ecological patterns of individual species is not fully understood. At the same time, species turnover makes that species composition varies from one area to another, so the analysis of biological interaction cannot be uncoupled from the spatial context. This is particularly important for clades that show high degree of specialization such as hummingbirds, where any variation in biotic pressures might lead to changes in morphology. Here, we describe the influence of biological interactions on the morphology of Hylocharis leucotis by simultaneously considering potential competition and diet resources. We characterized the extent of local potential competition and local available floral resources by correlating two measurements of hummingbird diversity, floral resources and the size of morphological space of H. leucotis along its geographic distribution. We found that H. leucotis shows an important morphological variability across its range and two groups can be recognized. Surprisingly, morphological variation is not always linked to local hummingbird richness or the phylogenetic similarity of. Only in the southern part of its distribution, H. leucotis is morphologically more variable in those communities where it coexist with closely related hummingbird species. We also found that morphological variation in H. leucotis is independent from the availability of floral resources. Our results suggest that abiotic factors might be responsible for morphological differences across populations in Hylocharis leucotis being biological interactions of minor importance.


Introduction
Species occurrence in local assemblages is ultimately determined by the joint influence of environment and biological interactions and their relative contribution is matter of debate (Diamond 1975, Sanderson et al. 2009, Harmon and Harrison 2015, Rabosky and Hurlbert 2015. The observed patterns of species co-occurrence were originally 2 interpreted in terms of competitive exclusion (Diamond 1975), however they might also be explained in terms of similarity and habitat preferences (Wiens 1989). Both ideas are not mutually exclusive and could act at different temporal or spatial scales in a nested fashion (Keddy 1992, Weiher and Keddy 1995, Weiher et al. 1998. At the same time, networks of biological interactions cannot be studied by using solely ecological methods because the composition of communities is not independent from the history of a region (Pigot and Etienne 2015). Moreover, dispersal might play an important role in community assembly as species distributions at continental scales are cohesive (but see Herrera-Alsina et al. 2018). From this perspective, the occurrence of species in a site is determined by the combined effect of evolutionary processes such as speciation and extinction, environmental conditions, interactions between species and rates of range expansion (Wiens andDonoghue 2004, Field et al. 2009).
An issue perhaps more intriguing and poorer understood is the feedback of species coexistence on evolutionary processes which gives spatial and temporal context paramount importance. The geographic mosaic theory of coevolution considers geographic heterogeneity as a key factor to describe the dynamics of the coevolutionary processes (Thompson 2009, Hembry et al. 2014). The coexistence of species can be described in two main clauses: 1) reciprocal natural selection pressure increases as more species locally co-occur and 2) the abiotic attributes of communities vary across the geographic range of the species (Thompson 2005, 2009, Benkman et al. 2010). According to this, variation in the strength of biological interactions produces coevolutionary hotspots or coldspots depending on whether there is a significant effect of reciprocal natural selection between the interacting species or not (Thompson 2005(Thompson , 2009). The intensity of biological interactions experienced by a given species could be different across its geographic range due to changes in local species composition which might lead species to undergo local specializations (Benkman et al. 2001). This is especially important because in turn, it could cause the ecological and evolutionary dynamics of species to follow different routes as their populations are spatially distributed under different conditions (Hanski 1998). Morphological variation across populations as a result of biological interactions (i.e. character displacement) has been suggested as a potential factor which triggers or completes speciation (Dayan and Simberloff 2005, Grant and Grant 2006, Pfennig and Pfennig 2009, 2010.
Some studies have addressed how morphological/functional traits could be influenced by interactions pairs or small groups of species. For instance, Herrera et al. (2006) found that regional variation in the local assemblages of pollinators is associated to variation of corolla traits of Lavandula latifolia (spike lavender). In the case of crossbills (Loxia curvirostra complex) and lodgepole pine (Pinus contorta spp. latifolia) it has been documented that an increase in seed defenses has favored an increase in bill size leading to strong divergent selection on crossbills which is mediated by the presence of a potential competitor (pine squirrel-Tamiasciurus hudsonicus) (Benkman et al. 2003(Benkman et al. , 2010. Nogueira et al. (2015) describes the complex relationship between extrafloral nectaries traits and functional properties of ant assemblage and how this interaction influences Anemopaegma album (Bignoniaceae). Stinchcombe and Rausher (2002) found that flowering plant ivyleaf morning glory Ipomoea hederacea shows evidence of selection on tolerance to deer damage which might depend on the presence of other natural enemies of this plant species. Although these studies contribute to the understanding of interspecific interactions and the potential feedback on a focal species, they all ignore the spatial context and the community structure as a whole (Althoff et al. 2014, Hembry et al. 2014. The high extent of ecological specialization in hummingbirds (McGuire et al. 2014, Sonne et al. 2019 suggests that any change in the pressure of biotic interactions can be reflected in morphological changes which makes them an ideal group to evaluate the influence of these interactions. Here, we explore whether morphological variation (i.e. morphological space) of Hylocharis leucotis (white-eared hummingbird) is influenced by the geographic mosaic of biological interactions in which it is embedded. In particular, we test whether Hylocharis leucotis has higher morphological variation in the localities where a) it experiences higher competition and b) there is a large diversity of floral resources.

Study system
We selected Hylocharis leucotis as focal species because its distribution in Mexico and Central American highlands is thoroughly documented and (different from other hummingbird species) its diet is very well described (Schuchmann 1999, Arizmendi andRodríguez-Flores 2012).
Hylocharis leucotis is a common species in highlands, allyear round resident inhabiting pine-oak, oak and pine-evergreen forest, clearings with flowers between 1200 and 3500 m a.s.l. It feeds and perches at low to mid-levels, often abundant along low banks of flowers. This species is distributed from the south of Arizona, northern Mexico, to Guatemala, El Salvador, Honduras and Nicaragua. Three subspecies have been documented: Hylocharis leucotis borealis, H. leucotis leucotis and H. leucotis pygmaea (Clements et al. 2018). We follow the nomenclature of Clements et al. (2018) who considered White-eared hummingbird part of the genus Hylocharis, although others authors considered that is part of the genus Basilinna (Schuchmann 1999, Gill andDonsker 2019).

Geographic data, floral resources and potential competition
We modelled the potential geographic range of H. leucotis as well as the distribution of all hummingbird species that inhabiting northern Central America and share at least partially the altitudinal interval of H. leucotis. We also modelled the potential distribution of 39 plant species that have been documented as dietary important for H. leucotis (Arizmendi and Rodríguez-Flores 2012). Hummingbirds and plant species included in the analysis are listed in Supplementary material Appendix 1. We obtained species' occurrence data from the global biodiversity information facility (GBIF) but we excluded records before 1950 as well as those records with geographic inconsistencies (i.e. records whose coordinates were evidently wrong, such as offshore records). Additionally, for hummingbird species, we compared many records to the known distributions of the species. In the case of plant species, we double-checked the Tropicos database (<www. tropicos.org/>) to look for any synonymy in the names of the plant species. Furthermore, we made sure that occurrence data coincided with the GBIF database. Finally, we consulted Dr Martínez-Gordillo, an expert in the Lamiaceae family (Martínez-Gordillo et al. 2013), to verify the authenticity of doubtful records (e.g. we found some records of Salvia elegans in GBIF database that were geographically inaccurate).
The climatic data were obtained from WorldClim 1.0 (Hijmans et al. 2005). Information on each set of variables in species modelling could be found in Supplementary material Appendix 1-2. We used genetic algorithm for rule-set production (GARP) model (Stockwell and Peters 1999) through the platform openModeller ver. 1.1.0 (Muñoz et al. 2011) to model the species' distribution. GARP generates and evaluates sets of rules representing nonrandom associations between climatic conditions in the localities known for a species, and those of the overall study region (Stockwell and Noble 1992, Stockwell and Peters 1999, Kobelkowsky-Vidrio et al. 2014). One of the appealing features of GARP algorithm is its high predictive accuracy Navarro 2009, Haverkost et al. 2010). Moreover, GARP is known to avoid the underestimation of presences (omission errors) but has a slight tendency to overestimation (i.e. commission errors; Peterson and Navarro 2009). We modeled only with spatially unique localities and with the 'best subsets' option implemented in openModeller to select the best supported model. The 'best subsets' are selected based on low omission and moderate commission errors (Anderson et al. 2003). For each species, we did 100 runs (i.e. generated models) with 70% of the points for training and 30% of the points for testing. Other GARP parameters for each species were: convergence limit of 0.01 or maximum number of iterations of 999; 20 models under omission threshold (soft omission, measured in percentage of predicted points); commission threshold of 50% (measured in percentage of pixels of the total area); and commission sample size of 999 (i.e. pseudo-absence data used to estimate commission). The GARP output was a consensus map from the 10 best models for each species (Kobelkowsky-Vidrio et al. 2014). open-Modeller provides a test (i.e. Accuracy Roc Score Kappa) to evaluate the quality of the generated models, the kappa score is a number between −1 and 1 where scores above 0.8 are generally considered good agreement, zero or lower means no agreement (see Results).
To obtain the final maps used further in our analysis, we reclassified and vectorised the raster outputs, considering the final vector of each species only the areas with more than seven consensus models (i.e. seven or more consensus models = 1 presence). In the final edition, we clipped each potential distribution map considering each species' accessible areas according to historic distribution and geographic barriers (Howell and Webb 1995, AOU 1998, Ridgely et al. 2003, Martínez-Gordillo et al. 2013. We used QGIS ver. 2.16.2 (<www.qgis.org/>) to visualize the entire set of distribution maps of hummingbird, we then overlaid grid of 552 equal-area cells (0.5 × 0.5° latitude and longitude near the equator) and only considered the area that matches the range of Hylocharis leucotis. Cells were regarded as local communities so we built a presence-absence matrix with species placed in rows and the grid cells in columns (55 rows and 552 columns). We repeated the same procedure for the plant species (39 rows and 552 columns). The geographic position where the museums' specimens of Hylocharis leucotis were collected was visualized in QGIS. Because sampled specimens were not evenly distributed in space, we consider the proximity of collection records to group individuals into operational geographical units (OGUs) to have a more homogeneous distribution. We defined 41 areas OGUs that included a minimum of 10 specimens (see morphological analysis section; Table 1, Fig. 1). Each OGU is integrated by a similar number of cells (ranging from 1 to 6) so the area is roughly similar. For the cells that included in each OGU's, we measured the morphological variation of H. leucotis and quantified the diversity of hummingbird species and floral resources.
The species richness of a given cell was calculated as the number of distributional maps that overlapped on this cell. In the case of plant species, species richness per cell is the abundance of floral resources within the range distribution of H. leucotis. By the same token, the total hummingbird richness per cell is the potential number of hummingbird species which H. leucotis coexists locally with, a figure that measures the potential competition experience by H. leucotis. If one OGU contains more than one cell, we calculate the mean of the diversity values across those cells in order to have a single value for the entire OGU. An alternative approach to assess the intensity of competition, is the calculation of the overall phylogenetic similarity (as a proxy for ecological similarity) between Hylocharis leucotis and the co-occurring hummingbird species for each OGU. Although exceptions have been documented (Uriarte et al. 2010), phylogenetic closeness is an accurate descriptor of competition intensity (Burns andStrauss 2011, Venail et al. 2014). We calculated the net relatedness index (NRI focal ) which, is the mean of phylogenetic distance between a given species and each one of the members of each OGU (Herrera-Alsina and Villegas-Patraca 2014). We used the phylogenetic tree from Bribiesca et al. (2018) which is in close agreement to the one proposed by McGuire et al. (2014) with the additional advantage that includes more species that coexist locally with H. leucotis.

Morphological analysis
We measured 693 bird skins, most specimens are housed at the Museum of Zoology 'Alfonso L. Herrera' (MZFC, UNAM) and the American Museum of Natural History (AMNH) but skins from other museums were included for a detailed coverage of the geographic distribution of Hylocharis leucotis (see Acknowledgments). We measured 11 linear variables that represent the size and shape of major functional modules of the bird external anatomy which are associated Table 1. ANOVA test for the morphological analysis of 11 linear variables that represent the size and shape of major functional modules of the bird external anatomy which are associated to the ecological traits (Claramunt 2010). The probabilities refer to one way ANOVAS, comparing the average of each OGU. P values: '***' = 0.001, '**' = 0.01, '*' = 0.05. 5 to ecological traits (Claramunt 2010): 1) bill length from the anterior border of the nostril to tip of the bill, 2) bill width at the level of the anterior border of nostrils, 3) bill depth (vertically) at the level of the anterior border of nostrils, 4) wing length to the longest primary, 5) wing length to the tenth primary, 6) length to the first secondary feather, 7) maximum tail length, 8) minimum tail length, 9) the width of the central rectrix at its mid-length, 10) tarsus length and 11) hallux length (including the claw).
We performed a global principal component analysis (PCA) with all the specimens and 11 variables and retained the seven variables that explain -globally -the most variation for each principal component. Those variables were bill length, wing length to the longest primary, length to the first secondary feather (wing width), maximum tail length, the width of the central rectrix at its mid-length (tail width), tarsus length and hallux length. We then performed a second PCA for each OGU with those seven variables. Using the eigenvalues of each PCA, the proper variance (Claramunt 2010) was computed for every OGU. The proper variance is an index for describing dispersion of morphological variation in a set of species, where the lower the value the less phenotypic variation in the morphological space.

Data analysis
Exploratory analysis of the data showed that all the variables were roughly normally distributed with nine outliers (identified by using the function outlierTest 'Car' package in R) which were excluded in posterior analyses making a total of 240 females and 444 males. Then all the variables were rescaled and centered, using the generic function scale' in R. To recognize discontinuous subsets of the morphological space throughout the range of H. leucotis we performed a cluster analysis using different methods: single linkage agglomerative clustering (single), complete linkage agglomerative clustering (complete), average agglomerative clustering (UPGMA, WPGMA. UPGMC, WPGMC), Ward's minimum variance clustering (Ward, Ward2) (Borcard et al. 2011). We used a cophenetic correlation and Shepard-like diagrams to find the best clustering method; the optimal number of clusters was found according to silhouette widths (Rousseeuw quality index), and a silhouette plot of the final partition of clusters (Borcard et al. 2011). We also performed a simple variance analysis (ANOVA) to look for differences in the variation of each morphological variable of Hylocharis leucotis across the OGUs. Finally we performed a Pearson correlation test between the proper variance values (i.e. local morphological variation of H. leucotis) and 1) floral resources diversity, 2) hummingbird species richness and 3) NRI focal.
All analyses were run in R ver. 3.5.0 (R Core Team) using the packages 'car' (Fox and Weisberg 2011), 'vegan' (Oksanen et al. 2018) and 'cluster' (Maechler et al. 2018) for the main analyses and 'picante' (Kembel et al. 2010) for the calculation of net relatedness index.

Results
Within the distribution of H. leucotis, the number of hummingbird distributions overlapping at local scale varied from 6 to 35, peaking at the highlands of south of Mexico and Central America and declining northwards (Fig. 2a, Supplementary  material Appendix 3). According to the built-in test of open-Modeller to evaluate the quality of models generated with GARP (i.e. Accuracy Roc Score Kappa), all the generated models have scores around 0.8-0.9 which means good reliability (Supplementary material Appendix 1). The frequency distribution of hummingbird richness had one single mode, local mean species richness was s = 18.05 (SD = 52.08) and showed a positive skewness (g1 = 0.32), indicating that in most cells the number of species was lower than average. Geographic ranges of the hummingbirds varied in size from 2 to 421 cells (mean = 148.36 cells). The number of overlapping ranges of floral resources varied from 1 to 39 cells, the highest number of ranges overlapping is located in central and southern Mexico and in the highlands of central America, declining towards the northern and southern extremes of the distribution of the H. leucotis (Fig. 2b, Supplementary material Appendix 3). Mean species richness was 27.90 (SD = 98.4) and showed right-skewed distribution (g1 = −1.070, i.e. most cells have a higher number of species than the average).
The cluster analysis shows that populations of H. leucotis are divided morphologically in two groups regardless of sex. However, in the case of females we found support for a gradual separation: one group is distributed from the Isthmus of Tehuantepec northwards whereas the other one stretches from the Isthmus of Tehuantepec to Nicaragua (Fig. 3). We found that males follow a similar pattern but the separation is more abrupt, starting from Guatemala and central America. The first group covers all the range of H. leucotis, whereas the second group is mainly concentrated in central America (Fig. 3). We found with the ANOVA analysis, that eight variables for males and seven variables for females have significant differences across the 41 OGUs (Table 1, Fig. 4a-b), suggesting that there is an important variation in the morphological space of H. leucotis and pointing towards a larger size in specimens collected in the northern part of its geographical distribution.
We found a positive but not statistically significant relationship (R = 0.174, p-value = 0.331) between the hummingbird local richness and morphological variation (i.e. proper variance) of H. leucotis (Fig. 5b) whereas the relationship between phylogenetic similarity and morphological variation was positive and not significant (R = 0.043, p-value = 0.809) (Fig. 5c). Moreover we found a positive relationship between plant richness and morphological variation (R = 0.232, p-value = 0.193) but it is not significant (Fig. 5a).
Because the cluster analysis suggested two main groups, we divided the OGUs into two sub-regions: north and other south of the Isthmus of Tehuantepec and repeated the analysis in each sub-region. We found that when analyzing the northern subregion separately, the conclusions are similar to those reached in the overall regional analysis (Supplementary material Appendix 4). In contrast, the southern sub-region showed different results than the regional analysis: we found a positive but not statistically significant relationship (R = 0.240, p-value = 0.503) between the hummingbird local richness and morphological variation of H. leucotis (Fig. 5e) whereas the relationship between phylogenetic similarity and morphological variation was positive and significant (R = 0.705, p-value = 0.022) (Fig. 5f). Finally we found a positive relationship between plant richness and morphological variation (R = 0.258, p-value = 0.470) but it is not significant (Fig. 5d).

Discussion
We found that there is an important morphological variation of H. leucotis across its geographic distribution and, intriguingly, this variation is not associated in general with neither The spatial distribution of hummingbird species within the geographic range of H. leucotis follows the general pattern reported for a variety of animal and plant groups where richness peaks in the tropics and decreases towards temperate zones (McGuire et al. 2014). Floral resources showed a similar pattern of distribution. Such overlap might suggest facilitation, however finding evidence of this interaction is difficult to prove (Graham et al. 2009, Lessard et al. 2016. According to the results, individuals in the northern part of the geographical distribution show a tendency to have a larger body size. This finding is in line with Schuchmann (1999)     . Morphological variation of Hylocharis leucotis and its association with local diet availability and potential competition. Phenotypic variation inside the morphological space of H. leucotis is summarized using proper variance (y-axis). Potential competition is measured using either number of coexisting species i.e. local richness (panels b and e) or phylogenetic similarity between H. leucotis and its coexisting species (c and f ). Diversity of floral resources corresponds (a and d) to the local richness of plant species that are documented to be part of the diet of H. leucotis. Top panels (a, b and c) show these relationships in the regional analysis whereas the panels in the bottom (d, e and f ) show these relationships in the sub-region south of the Isthmus of Tehuantepec. The analysis of sub-region north of the Isthmus of Tehuantepec is similar to the regional one (top panels) and can be found in the Supplementary material Appendix 4. Only panel (f ) shows an association with p < 0.05. 9 the north is larger, whereas H. pygmaea is smaller and thrives in the southern part of the range. From the cluster analysis, we found support for a gradual separation of two main groups of females and males that meet around the Isthmus of Tehuantepec, such gradual change agrees with the genetic differences in this species reported by Zamudio-Beltrán (2011).
Our study covered the vast majority of the geographical distribution of H. leucotis but unfortunately, there is one area with no specimens represented in museums and collections. This area is between OGU 5 and OGU 6, located on the Sierra Madre Occidental more specifically in the state of Sinaloa. We believe that this lack of sampling in this area is related to poor accessibility in combination with potential high rates of organized crime. According to proper variance values of OGU 5 (49.32) and OGU 6 (71.11), there might be an intermediate morphology that the analysis is missing (Fig. 1, Table 2). Furthermore, the diversity values of OGU 5 and OGU 6 are very similar which suggests that the contribution to this unsampled area to the overall pattern could be minor. We therefore expect that this missing information does not compromise our conclusions.
In general, we found that H. leucotis has a (non-significant) tendency to be morphologically more variable in local communities with higher hummingbird richness and also when it co-occurs with species that are phylogenetically more related, in other words the morphological space seems larger when interspecific competition is more intense. In those sites where H. leucotis coexists with phylogenetically closer species, the overlap in resource use should be high which leads to the focal species to develop a less specialized use of the resources. This suggests that competition plays a major role in community assembly by fostering the expansion of species' morphological space to make use of a wider range of resources (perhaps unexploited resources; Malpica et al. 2017). The response of species to intense competition could be the reduction of the morphological space associated with the specialization on a given resource. Alternatively, the increase in the amplitude of morphological space in sites with lower competitive pressure suggests competitive release (Moulton and Pimm 1986). However studies on vertebrate assemblages point to the expansion of morphological space (volume-increasing mechanism) as a more common process in natural communities (Ricklefs and Schluter 1993, Ricklefs and Miles 1994, Moreno et al. 2006. Because the analysis showed a lack of statistical support in the most of cases, those conclusions should be taken with caution. In the southern sub-region, phylogenetic similarity was found to be positively related to morphological variation which suggest that within a large region, there are areas subject to higher ecological pressures which is ultimately associated with biotic interactions. However, the difference in mechanisms operating at sub-regional level makes that the signal of the relationship morphology and phylogenetic similarity be lost at the regional analysis, which highlights the importance of studying ecological processes at different scales (García 2006). This suggests that biotic interactions play a minor role in the morphological variation of H. leucotis, and therefore its separation into subspecies is associated with other factors, at least in most of its geographical distribution. Perhaps the morphological variation we found could be the result of adaptation of the focal species to changing environmental conditions across its geographic distribution (Brown and Gibson 1983). Recently, Soteras et al. (2018) reported that morphological variation in sword-billed hummingbird is related to the geographic distribution of long-flowered plant species which contrasts to our results, this difference is likely to be due to the extent of specialization of these two hummingbird species. H. leucotis might have a less specialized bill morphology that enables it to feed on many different plant species, this makes its dependence on one single type of flower less strong when compared with the extremely specialized and long bill of sword-billed hummingbird. Our results suggest that character displacement is not a key mechanism for H. leucotis as local variation in morphology is independent from community structure. This finding is important because it contributes to the debate on diversity fostering/hindering rates of speciation. However, it can be the case that the morphology of a focal species is not only a function of hummingbird community composition but it also depends on the relative abundances of species. This could also be true for floral resources where the local richness of plant did not explain variation in morphology. It is important to note that we only assessed the effect of hummingbird species on H. leucotis morphology, however its foraging specialization could also be constrained or driven by the presence of non-avian animals such as bees, moths and butterflies. It is likely that the consumption rate of floral resources by insect species is lower than the consumption by hummingbirds because of the difference in size and energy requirement of these groups. This could suggests a minor contribution of insect species to hummingbird evolution, however, very high abundances of insects can have a much larger impact, and unfortunately sufficient data is not available to test the possible effect of competition across taxa (Sonne et al. 2019).
Here we found that 'antagonistic interspecific biotic interactions', or the 'availability of floral resources' are not associated with the morphological variation of the focal species. However it is possible that the greatest pressure for biotic interactions to which H. leucotis is subject is the 'intraspecific competition', so it would be interesting to know the local variation in the abundance of the same species throughout its distribution area.
Finally, the variability in the morphological space of H. leucotis that we reported could then be explained by other factors such as phenotypic plasticity ) and environmental differences along its geographical distribution (Hawkins et al. 2007). Because its large range size, H. leucotis is subject to a great variation in environmental conditions, this variability could influence the importance of biological interaction thus the signal of competition could be eroded and untractable at the spatial scale of our analysis. Some authors reported that at regional scale, the spatial distribution of body size in flycatcher assemblages shows a scaling pattern coincident with macroecological rules of Bergmann. In highland forests they found that morphological traits increase northeast to southwest and suggested that morphological variation is explained in part by the climatic gradients (Cortés-Ramírez et al. 2019). So is possible that morphological variation of H. leucotis is associated with Bergmann rule too as it shows an increase in body size from south to north. Further exploration on climatic gradient along the range size of focal species is needed.