Elevation is a stronger predictor of morphological trait divergence than competition in a radiation of tropical lizards

Adaptations for efficient performance are expected to shape animal morphology based on selection for microhabitat use and ecological forces. The presence of competitor species is predicted to cause niches to contract and enhance trait divergence. Therefore, increased species richness is expected to lead to greater trait divergence, and to result in reduced overlap and similarity between morphologies of sympatric species. We examined patterns of morphospace occupancy and partitioning in the skink fauna of New Guinea, the world's largest tropical island. Because skink species richness is largely decoupled from elevation in New Guinea, we could examine the effects of both factors (as proxies for competition and abiotic conditions), on morphospace occupancy and partitioning. We measured 1860 specimens from 79 species of skinks throughout Papua New Guinea, and examined their morphospace occupancy in a spatial context. We calculated for each assemblage within equal-area cells, the volume of morphospace occupied by all skinks, the mean volume occupied per species, and the mean distance and overlap between all species pairs. We then examined if these metrics are related to species richness and elevation. Elevation is a stronger predictor of morphospace occupancy than species richness. As elevation increases intraspecific variation decreases and morphologies become more similar to each other, such that overall morphospace occupancy decreases. Highland skinks are, on average, smaller, thinner, and shorter limbed than lowland species. We hypothesize that harsh climates in the New Guinea highland habitats impose strong selection on skinks to occupy specific areas of morphospace that facilitate efficient thermoregulation in sub-optimal thermal conditions. We conclude that the effect of competition on trait divergence on a community and assemblage scale is eclipsed by abiotic selection pressures in these harsh environments.


| INTRODUCTION
Niche partitioning is among the most fundamental processes in ecology, and is a major force generating phenotypic diversity.
Ecological opportunity promotes diversification and partitioning of the niche space among species (Losos, 2010). Competing species are thought to be unable to co-occur if they are too ecologically similar (Hardin, 1960), with a limiting similarity to coexisting species (MacArthur & Levins, 1967). Competition between taxa occupying similar niches could be reduced via partitioning of microhabitat use, activity times and diet-all of which are thought to manifest in morphological adaptation (Brown & Wilson, 1956;Grant & Grant, 2006;Pianka, 1974;Slatkin, 1980). That said, such community-wide character displacement is not always manifested (Meiri et al., 2011;Simberloff & Boecklen, 1981;Stuart & Losos, 2013), the importance of competition in shaping community structure has long been contested (Connell, 1983), and the degree to which competition-driven processes are more important than abiotic selective pressures is still unclear.
Competition is thought to increase as richness increases (e.g. Soule, 1966;Mesquita et al., 2007;Losos, 2009;cf. Meiri et al., 2014) and the niche space becomes more saturated (Pianka, 1974), in a process termed niche packing (Pigot et al., 2016). The theoretical result of a fully packed niche space is one where distances and overlaps between species are constrained to a degree that still allows coexistence of all co-occurring and ecologically interacting species (Hutchinson, 1957). Theoretically then, as niche space becomes more and more packed, distances and overlap between species will decrease until the niche space reaches full saturation (Pianka, 1974).
Therefore, PNG skinks are an excellent model to examine the effects of interspecific competition and climatic conditions on morphological trait divergence. We examine these patterns using species richness and elevation as proxies for competition and habitat and climatic conditions, respectively.
We measured 1,860 specimens of PNG skinks (Figure 1a), comprising 79 species out of 100 currently recognised in 'mainland ' PNG, to test the following hypotheses: 1. As skink richness increases, niches will become narrower, leading to individual species occupying a smaller area of morphospace ( Figure 1b).
2. The overall 'volume' of the occupied morphospace will not vary greatly with richness, but the relative distances and overlap between various species will diminish with increasing richness. That is, the size of the proverbial morphological 'pie' remains the same, but the size of the 'slices' that the 'pie' is partitioned into will change as niche packing occurs (Pigot et al., 2016): a. As skink richness increases, the size of the morphospace will not vary ( Figure 1c).
b. As skink richness increases, distances between species will decrease as the morphospace becomes more saturated ( Figure 1d).
c. As skink richness increases, trait overlap between species will decrease due to diffuse competition reducing overlap between closely competing species in species-rich areas (Pianka, 1974;Figure 1e).
3. Competitive exclusion (Hardin, 1960) and limiting similarity (MacArthur & Levins, 1967) will prevent similar species from cooccurring. Therefore, the mean distances between morphology of a focal species and sympatric species will be longer than the distances between its morphology and those of allopatric species (Figure 1f), and species morphology will overlap more with that of allopatric relatives than with sympatric ones (Figure 1g).
Alternatively, habitat filtering will result in sympatric species being more similar to each other than to allopatric ones.

FIGURE 1 (a)
A schematic representation of a skink with 14 morphometric measurements marked on it. The measurements are as follows: SVL: snout-vent length; HL: head length (from anterior edge of ear opening to snout); HW: head width (at the widest point); HD: head depth (at the deepst point); PeW: pectoral girdle width; PeD: pectoral girdle depth; PlW: pelvic girdle width; PlD: pelvic girdle depth; AL: arm length (shoulder to elbow); FaL: forearm length (elbow to wrist); ML: manus length (wrist to tip of longest digit); FL: femur length (from pelvis to knee); CL: crus length (from knee to heel); PL: pes length (from heel to tip of longest digit). Two additional measurements not displayed here are as follows: FrL: length of front limb (sum of AL, FaL and ML); HiL: length of hind limb (sum of FL, CL and PL). (b-e) Schematic representations of the morphospace of a hypothetical assemblage of three skink species, denoted by blue squares, red circles and yellow stars, each representing a different species. (b) Mean volume is calculated as the mean of all individual species' volumes of morphospace in the assemblage, represented here by coloured ellipses. (c) Total volume is calculated as the volume of morphospace occupied by individuals from all species in the assemblage, represented here by the grey shape. (d) Mean pairwise distance is calculated as the mean of all distances between centroids of all pairs of species in the assemblage, the centroids represented here by hexagons. (e) Mean overlap is calculated as the mean of all overlaps between all pairs of species in the assemblage, represented here as purple shaded polygons with dark outlines. (f-g) Hypotheses tested for the comparison of mean distance from and mean overlap with the focal species within and outside their ranges. The dashed line represents the null regression with an intercept of 0 and a slope of 1. (f) The blue-shaded area represents mean distances from the focal species that are higher within the focal species' range, that is, regression slope shallower than 1.
(g) The purple-shaded area represents mean overlaps with the focal species that are lower within the focal species' range, that is, regression slope steeper than 1. Thus, these shaded areas represent the expectations of the ecological competition paradigm   (Slavenko, Tamar, et al., 2020) and were excluded. The remaining 17 species were omitted due to scarcity of samples in the collections.
All have narrower distributions, on average, than included species (Supporting Information S1). We note that these 17 species are similar in body size and the proportions of lowland and highlands species to the pool of included species (Supporting Information S1).
Therefore, their exclusion from this study is unlikely to cause systematic bias in our results.
The border between PNG and Indonesia is only political and has no biogeographical or biological meaning. However, the Indonesian half of New Guinea is very poorly studied in comparison to PNG, leading to an artificial drop in species richness compared to PNG (see e.g. richness maps in Roll et al., 2017). We therefore limit ourselves to only examining the better-known eastern half of New Guinea, comprising the nation of PNG, without its outlying islands.
From each specimen, we measured 16 traits to quantify variation in body size and shape ( Figure 1a). We selected measurements relating to relative limb length, head size and shape, and body length and width. Relative limb length and body width are related to locomotor performance and clinging capability, and therefore to microhabitat use (Garland Jr. & Losos, 1994;Melville & Swain, 2000;Pianka, 1969). Head size is a strong proxy for bite force, and is often indicative of dietary preferences in lizards (Herrel et al., 2001;Van Damme et al., 1995;Verwaijen et al., 2002). Competition for microhabitats and food resources could both lead to trait divergence (Grant & Grant, 2006;Losos, 2009).
Limb measurements were all taken from the left side of the body, barring cases where the left limbs were damaged (e.g. missing digits), in which case all measurements for that specimen were taken from the right side of the body (186 specimens, 10% of total).
All measurements were taken using a digital caliper to the nearest 0.1 mm. We obtained distribution maps from Roll et al. (2017) for all the species included in the study. We overlayed the maps onto a 25 × 25 km equal-area Berhmann projection grid of PNG (sensu Tallowin et al., 2017), comprising 620 cells, excluding cells with <50% land cover to prevent edge effects. For each cell, we tallied the total skink species richness, and recorded which species occur in it.
We obtained elevation data from the Papua New Guinea Resource Information System (3rd ed.; Bryan & Shearman, 2008).
We then calculated mean elevation (in m a.s.l.) for each cell and assigned cells to either 'highland' or 'lowland', based on habitat categorisations in Bryan and Shearman (2015), who define lower montane forest from c. 1,000 m and above. This agrees with a breakpoint we uncovered in the relationship between skink richness and mean elevation at ~1,076 m a.s.l. (Supporting Information S1; Figure S1). Therefore, we used 1,076 m as the cut-off point for our categorisation, with cells below 1,076 m defined as 'lowland', and cells above 1,076 m defined as 'highland'. Similarly, we averaged mean elevation per cell (in m a.s.l.) across the distribution of each species and assigned species to either 'highland' or 'lowland' using the same categorisation scheme.
All the analyses were performed in R v4.0.2 (R Core Team, 2020) and the data and code are available in Dryad Digital Repository https://doi.org/10.5061/dryad.hdr7s qvgs .

| Morphometric analyses
We size adjusted all morphometric measurements apart from body length (Snout-Vent Length; SVL) by dividing them by SVL. Therefore, we had one non-transformed measurement, SVL (in mm), and 15 size-adjusted measurements (as proportion from each individual's SVL). We reduced dimensionality in the morphometric data by performing a PCA using the 'prcomp' function from the StAtS package in R, and we used the 'broken stick' method (Jackson, 1993) to select only those PC axes which explained a higher than expected portion of variance for downstream analyses. These PCs were used to define the morphospace of PNG skinks. We used the hypErVOLumE package (Blonder & Harris, 2018) to calculate the volume of morphospace occupied by each species (using the 'hypervolume_gaussian' function), measured as a Gaussian 95% kernel density (a method to estimate the probability density function) estimate around the coordinates of individuals of each species in morphospace.

| Spatial analyses of morphospace
We calculated four morphospace metrics for each cell ( Skink richness in PNG peaks along the northern versant and is largely decoupled from elevational gradients in parts of PNG below ~1,076 m, but is negatively correlated with elevation in parts of PNG above ~1,076 m (Supporting Information S1; Figure S1). We therefore examined the effect of both species richness (as a proxy for competition) and elevation (as a proxy for habitat and climatic conditions) on our four morphospace metrics.
Species richness is thought to reflect interspecific competition, under the assumption that as species richness increases, the number of competitors increases (see above). We therefore use the species richness of skinks as a proxy for the strength of interspecific competition, although we are aware that species identity and abundance may be more important determinants of the strength of competition.
Environmental conditions in PNG change predictably with elevation. Temperatures drop by roughly 0.5-0.6°C for every 100 m gained in altitude (Bourke, 2010). Largely due to the effects of changing temperature, vegetation structure, and therefore available microhabitat for skinks, also changes with elevation.
With increasing elevation, coastal mangroves give way to tropical broadleaf forests of diminishing canopy height, then giving way to shrubs, grasslands, and even tundra-like vegetation and bare rock at the top of the highest mountains (Allison, 2009;Bryan & Shearman, 2015). We therefore use elevation as a proxy for environmental conditions.
We ran Generalised Additive Models (GAM) for each metric using the 'gam' function from the mgcV package (Wood, 2011), with mean elevation and species richness as linear predictor variables, as well as spatial coordinates (latitude and longitude in decimal degrees) modelled as a thin plate spline to account for spatial autocorrelation (Wood, 2003 (Grömping, 2006).
As an added sensitivity analysis, we divided our dataset into two subsets, one of only lowland cells (<1,076 m; n = 408) and one of only highland cells (>1,076 m; n = 212). We repeated the above analyses for each subset of the data.
We calculated the centroids of PC scores for each cell, which we used to quantify the commonest area of morphospace occupied in each cell. We compared highland to lowland cells by performing a PERMANOVA analysis with the 'adonis' function from the VEgAN package (Oksanen et al., 2019), to test if there is a shift in PC centroids between lowland and highland cells. We generated histograms of the PC scores on each PC for each cell and examined changes in the distributions of PC scores with elevation by fitting linear regression models of mean elevation against the standard deviation, skewness and kurtosis of PC scores in each cell. Similar to the previously described analyses, we examined if the overall morphospace occupied per cell shifted between species in the highlands and lowlands. After assigning each species to either 'highland' or 'lowland' (as described in the 'data collection' section), we used a PERMANOVA analysis to compare morphospace occupancy between highland and lowland species.
Finally, we tested for signals of competition between species sharing similar areas of morphospace in their distributions (Figure 1fg). We did this by tallying which cells are within or outside the range of each species. We then calculated the mean distance and overlap from the focal species for each cell (i.e. the average of all distances or overlaps between each species in the cell and the focal species), and calculated the mean distances from the focal species across all cells within its range, and outside its range. So if, for example, species a occurs in cells i and j, but not in cells m and n, then we would calculate the mean distance and overlaps from species a for all other species occurring within cells i and j as the mean distance or overlap within the species' range, and the mean distance and overlaps from species a for all other species occurring within cells m and n as the mean distance or overlap outside the species' range.
We then performed two tests: (a) t tests comparing, for each species, the mean distances and overlaps between cells within the focal species' range and outside the focal species' range, and (b) a linear regression through the origin with mean distance/overlap within the focal species' range on the x-axis, and mean distance/overlap outside the focal species' range on the y-axis. In the linear regression, we tested whether the slope differed significantly from 1 (a slope of 1 meaning no difference in distances or overlaps within and outside focal species' ranges), our null hypothesis with no signal of competitive exclusion or limiting similarity. The alternative hypotheses we examined are a slope shallower than 1 for mean distances (meaning distances between species are higher within focal species' ranges), and a slope steeper than 1 for mean overlap (meaning overlap between species is lower within focal species' ranges), both of which we interpret as indicative of competition driving species to be more dissimilar in sympatry.

| Randomisation tests
Correlations between our morphospace metrics and species richness may exist due to random community assembly processes, even without the effect of competition. For instance, we might expect the total size of the occupied morphospace to increase with species richness even if species positions in morphospace are random in respect to one another, simply because there is a higher chance to 'sample' a species with an unusual morphology when more species are tallied.
To account for this, we generated null expectations of the relationship between morphospace occupancy and species richness based on a randomisation process.
For each value of species richness in our dataset, we generated 100 randomly sampled assemblages of species, in which species are sampled from the PNG pool with a sampling probability weighted by range size. Then, for each such random assemblage, we calculated the four morphospace metrics as described above-to represent morphospaces populated without competition or habitat filtering.
We then generated, for each random sample and each of the four morphospace metrics, full GAM models as described above (species richness, mean elevation, latitude and longitude). We extracted the regression coefficients with richness and elevation from each model, and thus generated, for each morphospace metric, a null distribution of the regression coefficients expected under a random process of morphospace occupation. We calculated p values for each metric as the proportion of the distribution as extreme as, or more extreme than the observed coefficient.

| RESULTS
The broken stick test determined that only the first two PCs explained a substantial portion of variance ( Figure S2). Cumulatively, these PCs explain 78.1% of the variance in morphology (PC1 61%, PC2 17.1%; Figure S3) Both skink richness and mean elevation were retained in the best models for three out of four morphospace metrics-all apart from mean morphospace volume (Table 1). Model R 2 values range between 78% (for mean pairwise overlap) and 96% (for mean morphospace volume). Mean morphospace volume, total morphospace volume and mean pairwise distance decrease with mean elevation. Mean pairwise overlap increases with mean elevation.
Mean pairwise overlap decreases with skink richness, whereas total morphospace volume and mean pairwise distance increase with increasing richness. Spatial patterns in morphospace metrics vary in importance-they explain between 11% of the variation (for total morphospace volume) to 78% of the variation (for mean morphospace volume). Mean pairwise distance and mean pairwise overlap vary more with longitude, whereas latitude explains most of the variation in mean morphospace volume. Mean elevation explains between 4.6 and 5.8 times more variation in total morphospace volume and mean pairwise distance than skink richness (Table 1). Only for mean pairwise overlap does species richness In panels g-h, the blue line represents the observed correlation coefficients from GAM models including both mean elevation and species richness as predictors, with spatial coordinates modelled as a thin plate spline. The blue shading around the line represents confidence intervals from the same models. The grey lines represent 100 regression lines from similar GAM models using randomised assemblages as the underlying data. Full summaries for the plotted models are listed in Table 1 Mean pairwise distance Mean pairwise jaccard similarity coefficient Mean volume Total volume 0 1,000 2,000 3,000 0 1,000 2,000 3,000 01 ,000 2,000 3,000 01 ,000 2,000 3,000    Figure S4; Table 1). The decrease in mean morphospace volume with mean elevation is likely not different from null expectations. However, the decreases in total morphospace volume and mean pairwise distance with mean elevation, and the increase in mean pairwise overlap with mean elevation, all differ significantly from null expectations (all p ≤ 0.02; Figure 3h; Table 1). We ran linear regressions through the origin, with mean distance/overlap within the focal species' range on the x-axis, and mean distance/overlap outside the focal species' range on the y axis  Figure   S6; Table S3). There is also a shift in the modal values, with the centroids of PC1 increasing and the centroids of PC2 decreasing with an increase in elevation (Figure 5b,c; Figure S6; Table S3).
This effect is much more pronounced for PC2 centroids, in which mean elevation explains 39% of the variance (Table S3). Similarly, species from lowland and highland cells occupy slightly distinct areas of morphospace (PERMANOVA based on 999 permutations, p = 0.001, R 2 = 0.012; Figure S4), but the difference is much smaller than the assemblage-level comparison, and the model has very low explanatory power.

| Little evidence of competition-driven trait divergence
Classic ecological theory predicts that, as the number of competitors increases, species morphology will adapt to reduce similarity due to the effects of diffuse competition (Hutchinson, 1957;Pianka, 1974).
However, we found only little evidence of competition-driven trait divergence in PNG skinks in traits linked to microhabitat use and dietary preferences. The morphospace does appear to be more tightly packed in species-rich areas (Table 1)  Comparisons of (a) mean distance from and (b) mean overlap with the focal species within and outside their ranges. The dashed line represents the null regression slope of 1. The blue line represents a linear regression through the origin between observed values within (0.99 ± 0.01, p = 0.26, R 2 = 0.99) and outside (0.98 ± 0.02, p = 0.34, R 2 = 0.98) focal species' ranges, with confidence intervals represented by the shaded area.
p values correspond to the null hypothesis that the slope = 1. Each dot is coloured based on the significance of t tests comparing the metric within and outside the focal species' range-turquoise for significant difference between ranges (p < 0.05), and pink for no significant difference between ranges result may suggest that some elements of niche packing occur (Pigot et al., 2016), but our randomisation models have uncovered that the patterns of richness-related morphospace occupancy do not differ from null expectations (Figure 3g; Figure S5). Since our null models explicitly fill the morphospace randomly with regard to species morphology, this would suggest that the observed patterns of the morphospace becoming larger, and more tightly packed, in species-rich areas, need not reflect any competition-driven process.
Morphospace occupancy appears to be more strongly influenced by elevation, a proxy for environmental conditions and habitat type, than it is by species richness, a metric for interspecific competition ( Figure 3; Table 1). Furthermore, there is little evidence of competitive exclusion or limiting similarity having driven assemblage structure ( Figure 4; Table S2). Therefore, in PNG skinks, similar to what has been described in several other cases (e.g. Meiri et al., 2011;Stuart & Losos, 2013), the effects of competition on shaping trait divergence appear weak. There are several possible explanations for the weak effect of competition: First, there may be other axes of niche space that are not captured by the morphometric measurements used here, and competition may affect some axes of trait space more strongly than others.
For instance, several species of skinks from the genera Prasinohaema and Lipinia have developed subdigital lamellae, likely in association with their arboreal habits (Greer, 1974). Traits such as these, which have well-documented impacts on locomotory performance (Irschick et al., 1996;Macrini et al., 2003), are likely important in microhabitat niche partitioning but were not accounted for in our analyses.
Second, the method we employed here assumes that all species that occur within a 25 × 25 km cell are sympatric and syntopic, an assumption that is not here tested due to the resolution of the data used, and very likely is not true for all cells in PNG. For instance, skink species on Mt Kiandi are separated by microhabitat and elevation (Allison, 1982), meaning that despite all co-occurring in the same cell, few of them interact ecologically with one another. Furthermore, microhabitat and temporal partitioning (assuming these do not by themselves influence morphology, which is logical at least for the temporal niche) could allow even syntopic species with a similar morphology to cohabit. Species richness may therefore be an inadequate measure of competition in this system. The extent to which species in the same cell interact with each other across PNG, and therefore experience competition, and whether patterns of co-occurrence change spatially and with elevation, need to be determined.
Third, competition may simply be less important in shaping trait divergence in our examined traits compared to abiotic selective pressures. In fact, if strong selection exists, co-occurring species may be even more similar than expected by chance due to environmental filtering. Our models instead strongly support elevation, a proxy for habitat and climatic differences, as a stronger driver of trait divergence.

| The size of occupied morphospace changes with elevation
We found that highland species occupy a smaller portion of morphospace than do lowland species (Figure 3d,h; Figure S4). Indeed, many of the most morphologically distinct skinks in PNG, which occupy the margins of the morphospace, are lowland species-these include the semi-aquatic and unusually long-limbed Fojia bumui, the extremely large Tiliqua gigas, or several species of large-headed Tribolonotus, noted for their strongly keeled scales and bony protrusions. However, despite more morphologies being represented there, the morphospace in the lowlands does not appear to be more tightly packed compared to the highlands, as overlap increases with increasing elevation and decreases with species richness (Table 1).
An important caveat is that the species in this study are not phylogenetically independent, and failing to account for phylogeny may impact results of morphological analyses, since both time from divergence and the rate of trait evolution would affect the accumulation of trait differences. If phylogenetic niche conservatism (Blomberg et al., 2003;Wiens & Graham, 2005)

| Lowland and highland species occupy distinct areas of morphospace
As well as occupying a smaller portion of morphospace in general, highland skinks also occupy a somewhat different area of morphospace compared to skinks from lowland cells ( Figure 5). Such differences in morphospace could arise from several distinct mechanisms: for instance, environmental differences could lead to different selection regimes driving morphological evolution (Pinto et al., 2008).
Similarly, ecological opportunity and the evolution of key innovations might lead to colonisations of novel areas of the adaptive landscape (Losos, 2010). However, highland species in PNG do not so much represent a distinct area of morphospace, but rather a diminished subset of the morphologies found in the lowlands ( Figure S4). Therefore, the differences between lowland and highland assemblages are derived from a shift in the distributions and modal values of morphologies in the assemblages (Figure 5b-c), rather than the presence of unique morphologies in highland assemblages.
As elevation increases, the habitat changes drastically, particularly with reduced stratification as canopy height decreases and eventually disappears entirely above the tree line (Bryan & Shearman, 2015). Skinks from highland cells are, on average, smaller, thinner and with shorter limbs compared to skinks from lowland cells (Figures 2 and 5). Relative limb length is related to locomotor performance, and therefore to microhabitat use (Garland Jr. & Losos, 1994). Longer limbs are typically associated with faster sprinting speeds and enhanced jumping capability (Irschick et al., 2005;Melville & Swain, 2000;Pianka, 1969). Therefore, species occupying open habitats tend to have longer limbs (Melville & Swain, 2000;Goodman et al., 2008;cf. Schulte II et al., 2004).
Conversely, arboreal and ground-dwelling species in leaf litter and dense vegetation tend to be smaller and shorter limbed (Melville & Swain, 2000). Many of the highland PNG skinks occur in alpine grasslands, shrublands and tree-fall gaps (Allison & Greer, 1986;Greer et al., 2005), and locomotion in dense grasslands might provide selective pressure for short relative limb length.
Another possible avenue for selection on body size and shape in highland skinks is climate. The lapse rate in PNG is roughly 0.5-0.6°C for ever y 10 0 m gained in altitude (Bourke, 2010), meaning that montane regions of New Guinea are typified by extremely cold temperatures with occasional frosts, with mean temperatures dropping as low as 5-6°C (Allison, 2009), and nightly temperatures reaching freezing levels. These harsh climates are far from ideal for ectotherms such as lizards. Many montane species exhibit adaptations for cold climate such as ovoviviparity (Greer et al., 2005). A relationship between body size and climate has long been contended, although evidence for such a relationship in reptiles is weak (Slavenko et al., 2019). The driving mechanism behind this change is purported to be thermoregulation-small body sizes lead to increased surface area-to-volume ratios, which in ectotherms such as lizards may improve heat gain in cold climates (Carothers et al., 1997). Thermoregulation may be extremely difficult in montane habitats (Monasterio et al., 2009), and therefore needs to be especially efficient (Ortega et al., 2016). In high elevations, low ambient temperatures lead to short windows available for activity. Therefore, ectotherms require sunbathing to achieve optimal temperatures, and the more efficient and faster heat gain is, the more effective the use of this window for activity.
We suggest that the harsh climatic conditions of montane habitats in PNG drive morphological evolution. Morphologies that are suboptimal for thermoregulation may be adequate in the mild climates of the lowlands, but lead to reduced performance in harsh highland climates, particularly in intraspecific interactions such as competition for limited food resources or mating opportunities (Irschick et al., 1996;Macrini et al., 2003;Pinto et al., 2008). This hypothesis could be tested by comparing the thermoregulatory efficiency of species with different morphologies, or even using inert models to compare passive rates of heat acquisition and retention.

| CONCLUSIONS
Competition is classically considered to be a foremost driver of niche partitioning (Hutchinson, 1957;Pianka, 1974). Ecological theory predicts that competitors cannot coexist in the same niche (Brown & Wilson, 1956;Grant & Grant, 2006;Hardin, 1960;Pianka, 1974;Slatkin, 1980), leading in many cases to morphological trait divergence (Grant & Grant, 2006). However, evidence suggests that competition is not always the strongest driver of niche partitioning or trait divergence (e.g. Connell, 1983;Meiri et al., 2011;Stuart & Losos, 2013). Our results support the latter conclusion, and that abiotic selective pressures in harsh environmental conditions may be more important and overshadow the effects of biotic interactions in determining trait evolution and divergence. In particular, we think that shifting climatic regimes with elevation drive adaptation to montane habitats in our study system of PNG skinks, and are a strong influencer of morphological divergence. We strongly recommend that competition not be taken to be a strong driver of community assembly and trait divergence a priori, and that its effects should always be quantified and clearly tested and contrasted with the predicted and modelled effects of abiotic pressures. Guinea who provided access to their ancestral lands.

AUTHORS' CONTRIBUTIONS
A.S. collected data and designed and performed the statistical analyses; A.A. and S.M. both contributed to supervision of the project and conceptual design; A.S. wrote the first draft of the manuscript, and all authors contributed substantially to revising the writing and study design. All authors gave final approval for publication.

DATA AVAIL ABILIT Y STATEMENT
All data used in the study and code to run the analyses are available from Dryad Digital Repository: https://doi.org/10.5061/dryad.hdr7s qvgs .