Using crowdsourced photographic records to explore geographical variation in colour polymorphism

In wide‐ranging taxa occupying broad climate gradients, alternative colour pattern variants (i.e. morphs) may emerge as an adaptive response to local climatic regimes. We explore the patterns of geographical segregation among colour morphs of the Lace monitor and determine whether morphs occupy divergent climatic niches along a mesic–arid gradient.


| INTRODUC TI ON
Polymorphic species are those in which multiple distinct, genetically determined phenotypic variants (i.e. morphs) coexist in an interbreeding population, with the rarest morph occurring too frequently to be solely the result of recurrent mutation (Ford, 1945;Huxley, 1955). A ubiquitous form of polymorphism is colour polymorphism-the presence of distinct, heritable colour morphs that are fixed (i.e. they do not vary with ontogeny, physiology or sex; Gray & McKinnon, 2007;Mckinnon & Pierotti, 2010;Roulin, 2004;Roulin & Bize, 2007). Discrete, visually discernible polymorphisms provide clear phenotypic markers that are generally simple to quantify or categorise; hence, they have been used to investigate evolutionary processes long before the advent of modern molecular biology (Svensson, 2017). Colour polymorphisms can be maintained via multiple evolutionary processes acting within and between populations, including drift, gene flow and various forms of selection (both natural and sexual), such as disruptive selection or negative frequency-dependent selection (Gray & McKinnon, 2007;McLean & Stuart-Fox, 2014;Svensson, 2017).
In colour-polymorphic taxa, multiple morphs may coexist over much of a species' range, with geographical variation in the type, number and relative frequency of colour morphs (McLean et al., 2015;Pérez i de Lanuza et al., 2018;Runemark et al., 2010;Schilthuizen, 2013;Wunderle, 1981). In many instances, morphs are distributed non-randomly over an environmental cline, suggesting that each morph may have evolved to exploit alternative ecological niches (Da Silva et al., 2016;Fisher-Reid et al., 2013;Hoekstra et al., 2004;Ożgo, 2011;Skulason & Smith, 1995;Smith & Skúlason, 1996;Van Valen, 1965). This clinal variation in polymorphism is sometimes explained by ecogeographical rules and hypotheses, such as the tendency for heavily pigmented (darker) morphs to be found in more humid regions (Gloger's rule, reviewed in Delhey, 2019), or the tendency for darker morphs to occur in colder areas because low skin reflectance improves heat absorption (thermal melanism hypothesis, reviewed in Clusella Trullas et al., 2007).
Although the distribution, causes and consequences of colour polymorphism are well studied in some species, such as the brownlipped banded snail Cepaea nemoralis (Ożgo, 2011;Silvertown et al., 2011), primary data on the geographical distribution of morphs are scarce for most polymorphic species. If the geographical distribution of a polymorphism can be established, then such information can be used to test hypotheses about its evolution (McLean & Stuart-Fox, 2014). However, colour-polymorphic species often occupy many different habitats over broad geographical ranges (Delhey et al., 2013;Forsman & Åberg, 2008aForsman & Åberg, , 2008bGaleotti et al., 2003;Takahashi & Noriyuki, 2019); therefore, field studies which aim to document the geographical distribution of morphs in wide-ranging polymorphic species can be costly, time-consuming and hence impractical. This problem can be overcome by quantifying colour morphs depicted in georeferenced photographs of organisms, which are growing in availability through the advent of online, citizen science biodiversity databases (Mesaglio & Callaghan, 2021).
One of the most successful photograph-based citizen science initiatives is the iNaturalist organism occurrence recording tool (iNatu ralist.org), wherein georeferenced photographs of wild organisms are submitted by amateur naturalists, and other users verify the identification of each organism (Mesaglio & Callaghan, 2021). This adds a powerful new dimension to the modern practice of biogeographical research. Whereas conventional online biodiversity databases tell us the spatiotemporal occurrence of taxa only, iNaturalist tells us the spatiotemporal occurrence of taxa and their phenotypes.
The tremendous scientific potential of these data is becoming widely realised, and high-quality phenotypic datasets derived from iNaturalist photographs are increasingly being used to test ecogeographical hypotheses about the evolution of colour polymorphism (e.g. Drury et al., 2019;Lehtinen et al., 2020;Moore et al., 2021).
In Australia, the iNaturalist participation rate (number of observations, species and identifications) has increased exponentially since 2016, with particularly high record densities in eastern and south-eastern Australia, coinciding with major human population centres (Mesaglio & Callaghan, 2021). Accordingly, spatial variation in colour morphs may be well documented in polymorphic species that are common and whose distribution spans this region. One such widely distributed colour-polymorphic species is the Lace monitor Varanus varius (White, ex Shaw, 1790), a varanid (family Varanidae) distributed in eastern and south-eastern Australia and across a diverse environmental gradient (tropical, subtropical, cool temperate, semi-arid and arid environments). Individuals of V. varius exhibit one of two discrete colour morphs without any intermediates between them (Brown, 2014), which are easily visually discerned. The most commonly encountered phenotype is the 'Lace' morph, which has a dark grey to dull bluish black ground colour with pale, fine laceworklike patterning (Figure 1a). The second morph consists of few relatively broad, boldly contrasting yellow and black transverse bands ( Figure 1b). This morph is commonly referred to as the 'Bell's' morph, after London herpetologist Thomas Bell (Duméril & Bibron, 1836).
Based on extensive captive breeding, we know that both morphs can polymorphisms without needing to conduct expensive, time-consuming field studies on wild organisms.

K E Y W O R D S
climatic niche divergence, colour pattern polymorphism, COUE, iNaturalist, Lace monitor, mesic-arid gradient, n-dimensional hypervolume, point-pattern process, varanid, Varanus varius occur in a single clutch; the Bell's morph appears to be inherited autosomally as a Mendelian dominant trait, with the Lace morph being recessive (Brown, 2014). The geographical context of this colour pattern polymorphism, and its ecological basis, has never been explored in detail. But notably, Horn (1980) collated 33 historical records of the Bell's morph and presented a simple map showing that virtually all localities of this morph existed inland of the Great Dividing Range (GDR)-a region which is comparatively more arid than Australia's mesic eastern and southern coastal areas.
In ectotherms, geographical variation in the composition and relative frequency of colour morphs has been linked to climate (Friedman et al., 2017;McLean et al., 2015;Pérez i de Lanuza et al., 2018), suggesting that climatically driven selection produces multiple phenotypes adapted to different climatic niches Lattanzio & Buontempo, 2021;Pérez i de Lanuza et al., 2018). For instance, climatic variables such as temperature, humidity and solar irradiation have been suggested to be major selective agents shaping ecogeographical patterns in the colour of body surfaces relevant for thermoregulation in ectotherms (Clusella Trullas et al., 2007Goldenberg et al., 2021;Kuyucu et al., 2018). Hence, if colouration is under locality-specific selection associated with some gradient in climate, it is likely that colour morphs will have divergent realised climatic niches along such gradients. This could explain why morphs of V. varius purportedly vary in relative frequencies across the species' distribution.
The aim of this study is to demonstrate the potential that iNaturalist photographs have for shaping the future of biogeographical research, while simultaneously answering a specific research question relating to a single species. We hypothesised that the colour pattern morphs of V. varius are generally segregated in geographical space, and that such segregation would be associated with realised climatic niche (herein referred to as "niche") divergence between morphs. Specifically, given the earlier observations of inland Bell's specimens (Horn, 1980), we predict that Bell's morph will occupy a more arid niche relative to the Lace morph, which may tend towards more mesic climates. To test this hypothesis, we first established the geographical context of the colour polymorphism by mapping occurrence records of each morph derived from crowdsourced photographs. We then tested for geographical divergence of morphs under an exploratory point-pattern process (PPP) analysis framework. We then analysed the topology of each morphs' niche in univariate and multivariate environmental space. And finally, we tested the equivalence and similarity of the morphs' niches. In doing so, we investigated whether there has been significant divergence among morphs' niches, and what climatic parameters are primarily associated with the niche divergence. We discuss the possibility that morphs are alternatively adapted to their respective climate niches as a result of selection on thermoregulation and crypsis.

| Occurrence data collection
We obtained all research-grade records of V. varius on iNaturalist by downloading data through the Global Biodiversity Information Facility on 28 June 2021 (GBIF.org, 2021), resulting in an initial occurrence dataset of 1936 records. After removing records with coordinate uncertainty >5 km (n = 324) and records with photographs that did not legibly depict the animal's colour morph (n = 17), the photos uploaded with each record (n = 1595) were visually examined (by J.E.F. only, to avoid inter-observer bias) and individuals were scored as being either 'Lace' or 'Bell's' morph.
This resulted in 1536 records of the Lace morph and 59 records of Bell's morph.
To provide a more complete picture of the extent of occurrence of the Bell's morph, we supplemented our iNaturalist dataset with additional records from other sources. These additional records were obtained by contacting people with photographs of wild Bell's specimens on social media platforms (Facebook, Instagram and Flickr), and incorporating into our dataset those observations which the observer could provide confident coordinates (uncertainty <5 km) for the specimen. Notably, we impartially obtained all additional records irrespective of their location. We did not seek out additional records of the Lace morph from social media. This is because photographs of the Lace morph on social media significantly outnumbered those of the Bell's morph, making it infeasible to contact the owner of each photograph for meta data. Moreover, our iNaturalist dataset was already over-represented by Lace morph occurrences (96.3%); hence, these additional Bell's morph records reduce class imbalance. We argue that inconsistent sourcing of Bell's and Lace morphs records introduced negligible biases to posterior analysis, as iNaturalist and social media sources share a similar data collection mechanism (i.e. opportunistic sightings). The final occurrence dataset contained 1637 records, of which 1537 were Lace and 100 were Bell's morph (Supporting Information S1).
Occurrence data derived from iNaturalist and social media are typically recorded as unstructured opportunistic observations, which may bring about sampling biases in space and time. For instance, citizen science data tend to be over-sampled in areas closer to human population centres and transport infrastructure (Di Cecco et al., 2021;Hugo & Altwegg, 2017;Tiago et al., 2017). Moreover, species or phenotypes perceived as rare, colourful or otherwise interesting may be preferentially recorded (Caley et al., 2020;Isaac & Pocock, 2015). However, we argue that sampling biases concerning recorder behaviour are largely irrelevant or will pose negligible impact on our study. For example, iNaturalist observations are less dense in inland parts of Australia, owing to the lower activity of humans in these more remote areas (Mesaglio & Callaghan, 2021).
Though this may explain geographical biases in the number of records of V. varius across its range (e.g. more records in coastal areas), the propensity of observers to record a given morph should be mostly homogeneous across the species' range. It is in this respect that our measurement of spatial variation in colour polymorphism is expected to be negligibly biased for the purpose of our analyses.
Here, we removed records with identical coordinates, zero-value coordinates (records with absolute-zero values for latitude and longitude), equal coordinates (records with equal latitude and longitude values), coordinates falling in the sea and records from outside a temporal period of 1970 to present. Finally, to minimise the effects of sampling bias and spatial autocorrelation to subsequent analyses, we geographically filtered the quality checked occurrence datasets using 'spThin' package (Aiello-Lammens et al., 2015). Here, separately for each morph, a pairwise distance matrix of occurrence records is calculated, and records with the greatest number of neighbouring occurrences were removed. This process was iterated until there were no records of a nearest neighbour closer than 5 km. The final thinned occurrence dataset comprised 512 (87.5%) records for the Lace morph and 73 (12.5%) records for the Bell's morph. Bank, 2020) describing solar irradiance that directly reaches a surface for 1994-2019. We downloaded climatic variables at a resolution of 30 arcseconds (approximately 1 km 2 at the equator). From these climatic variables, we then selected an optimal subset of variables that align with varanid and reptile ecological theory ('stateof-the-art' predictor variable selection approach, see Petitpierre et al., 2017), are known to influence distributional patterns of colour polymorphism in other ectothermic species  and have an optimal variable correlation structure (Dormann et al., 2013). In reducing the dimensionality of environmental variables, we selected only the putatively ecologically most relevant variables (Lin & Wiens, 2017;Quintero & Wiens, 2013), while also minimising variable intercorrelation (Austin, 2002;Dormann et al., 2013). Variable correlation structure was quanti-

| Environmental background
For each morph, we defined the spatial extent of its environmental background as all terrestrial biomes it has presumably sampled -that is, the terrestrial ecoregions (Dinerstein et al., 2017) intersected by the morph's geographical range.
We defined the geographical range of each morph based on a 100-km-radius buffer of its geographically thinned occurrence records (described above). We then randomly sampled approximately 10,000 background points across the spatial extent of environmental backgrounds (Barbet-Massin et al., 2012).

| PPP analysis
To assess spatial segregation among morphs, we conducted an exploratory spatial data analysis of the distribution of the Lace and Bell's morph under a PPP analytical framework using 'spatstat' package (v.2.2; Baddeley & Turner, 2005). Here, we defined each morphs' PPP using the geographically thinned occurrence records as points, the associated colour morph identifier (i.e. Lace or Bell's) of each record as marks, and the species' geographical range as the sampling window. In contrast to each morph's geographical range, the species' geographical range was based on a 100-km-radius buffer of pooled (Lace and Bell's morph) geographically thinned occurrence records.
We first tested for complete spatial randomness of each morphs' PPP using the chi-square test (two-sided Pearson's χ 2 goodnessof-fit test) based on quadrat counts. Here, the sampling window is divided into 10 × 10 quadrats with equal sizes. The Pearson's χ 2 statistic is then computed based on the intensity of points in each quadrat versus random expectations.
We then tested for independence of both morphs' PPPs using a random labelling test (H 0 = the distribution of each colour morph's points are conditionally independent and identically distributed).
Here, 999 simulated PPPs were generated from both morphs' pooled occurrence records, holding their locations fixed and randomly resampling their marks with replacement. For both observed and simulated PPPs, we computed for G ij -function (aka G 'cross-type' or 'i-to-j' function) of the points of Bell's morph relative to that of the Lace morph, and vice versa. The Finally, we analysed spatial segregation of morphs across the species' range. We first spatially smoothed the geographically thinned occurrence records through Gaussian kernel density smoothing (bandwidth = 100 km). We then developed bias maps by computing the relative proportions of intensity of points in space based on the ratios of spatially smoothed occurrences of the Bell's versus Lace morph. We visually analysed negative association of the intensity and thereby spatial segregation of morphs across the species' range.
To quantitatively support this, we tested for covariation in the distribution of the Bell's and Lace morph using Kolmogorov-Smirnov test.

| Climate niche analysis
Inspired by the studies of Pili et al. (2020), Rödder et al. (2017) and Tingley et al. (2014), we analysed similarities and differences between the Lace and Bell's morphs realised climatic niches through univariate and multivariate analysis. We conducted a univariate analysis by computing the density profiles of climatic variables sampled in the locations of geographically thinned occurrence records. We then conducted a multivariate analysis using the COUE (Centroid shift, Overlap, Unfilling and Expansion) framework . Because different statistical frameworks can yield variable conclusions, we complemented the findings of COUE by conducting a multivariate analysis using a competing framework-the n-dimensional hypervolume (Blonder et al., 2014(Blonder et al., , 2018).

| Defining the topology of niches in principal component environmental space
We analysed both morphs' niches in environmental space defined by the first two (COUE framework) or four (n-dimensional hypervolume framework) principal components of a PCA. We calibrated the PCA using the seven climatic variables (described above) sampled at the pooled geographically thinned occurrence records and environmental background points of both morphs. These first two or four principal components, respectively, capture 88.95% or 97.88% of the variation in the climatic data. The correlation structure of climatic variables with the principal component axes is shown in Figure S2.1.

| COUE framework
We primarily analysed similarities and differences in both morphs' niches under the COUE framework , using the 'ecospat' package (v3.2. Di Cola et al., 2017). Here, we first projected the PC scores of occurrence records and environmental backgrounds onto a two-dimensional grid consisting of 100 × 100 cells and defined by the first two principal components (herein 'PCA grid'). We then smoothed the densities in each grid cell using a Gaussian kernel density function with a standard bandwidth (Silverman, 1986). Based on the topology of both morphs' climatic niches as depicted in the PCA grid, we decomposed niche stability (i.e. proportion of the Bell's morph's niche overlapping with that of the Lace morph's), expansion (i.e. proportion of the Bell's morph's niche non-overlapping with that of the Lace morph's) and unfilling (i.e. proportion of the Lace morph's niche non-overlapping with that of the Bell's morph's) between the niches of the Bell's versus Lace morph in analogue environmental space. We calculated these COUE indices in the intersection of the 95th percentile of environmental backgrounds of each morph to remove marginal climates resulting from kernel smoothing. As such, these indices are not confounded by non-analogue climates, making measured expansion indicative of difference in the realised niches .
We then quantified the overlap of the climatic niche of the Bell's versus Lace morph using Schoener's statistic D and Hellinger distance I for niche overlap (Schoener, 1970;Warren et al., 2008), which estimates the overall similarity between the two niches over the whole environmental space (Warren et al., 2008). Schoener's D and Hellinger I ranges from zero (no overlap) to one (complete overlap).
We then quantified niche conservatism by testing estimates of D and I under two alternative hypotheses, representing two extremes across a spectrum of niche conservatism: niche equivalency, which determines whether realised climatic niches of both morphs are conserved in the strictest sense (i.e. whether observed niche overlap is effectively indistinguishable when randomly reallocating the occurrences of both entities among the two ranges), and niche similarity, which estimates whether the realised climatic niche occupied by one colour morph is more similar to the niche occupied by the other than would be expected by chance Warren et al., 2008).
We tested the hypothesis of niche equivalency by randomly splitting pooled occurrence records of both morphs into two datasets, maintaining the same number of occurrences as observed in the original datasets, and estimating D and I. We iterated this process 999 times to generate a null distribution of simulated D and I values. We then compared the observed D and I values with the distribution of that of simulated values; we assumed niche equivalency (i.e. the realised climatic niches of both colour morphs are more similar than two random niches) if the observed overlap value fell within the density of 95% of that of null distribution.
Meanwhile, in testing for niche similarity, we first simulated a random niche by randomly shifting the entire observed occurrence density of one morph (e.g. Lace). We then calculated the overlap of the simulated niche to the observed niche of the other morph (e.g. Bell's). We tested niche similarity in three ways: (1)
Here, we quantified the Lace and Bell's morphs' climatic niches based on multidimensional hypervolumes of random records.
These random records were derived from the Gaussian kernel density estimates of the distributions of each morph's occurrence records in multidimensional environmental space defined by the first four components of PCA described above. The bandwidth for Gaussian kernel density estimates was estimated using the Silverman bandwidth estimator. We decomposed niche differences and similarity using hypervolume distance and intersection metrics (minimum distance, centroid distance, volume of the intersection, the unique fraction of hypervolumes). We also computed the similarity of niches in multidimensional environmental space using Jaccard and Sorensen similarity indices. Jaccard and Sorensen similarity indices range from zero (not similar) to one (similar). We calculated the hypervolume metrics and similarity indices at the 95th percentile probability boundary of hypervolumes (to remove marginal climates resulting from kernel smoothing).

| Geographical distribution of colour morphs
Visual inspection of occurrence records of colour morphs in geographical space ( Figure 2) indicates that, while the Bell's morph can be found across much of the species' geographical range, it is relatively more widespread and frequent west of the GDR, whereas records of the Lace morph are more concentrated along vast stretches of the east coast, such as that from Sydney to Melbourne. Notably, there were several observations in our dataset in which both Lace and Bell's specimens were captured in a single photograph (e.g. at Shepparton, Victoria; Frogmore, New South Wales; Woodgate, Queensland), demonstrating not only that the two morphs overlap geographically, but that they are sympatric in some locations, as per the definition of a true polymorphism (Ford, 1945;Huxley, 1955).

| Spatial point-pattern analysis
Chi-square test indicated that both morphs' PPPs are non-homogenous and cluster in space (Table 1). Similarly, random-labelling randomisation test showed the Bell's and Lace morphs' PPPs are conditionally nonindependent and non-identically distributed (Figure 3a). Finally, spatial bias maps (Figure 3b,c) and Kolmogorov-Smirnov test (Table 1) showed spatial segregation of morphs across the species' range.

| Univariate climate niche analysis
Comparison of density profiles of climatic variables at the occurrence locations of Bell's versus Lace morphs revealed that the Bell's morph is nested within the variable range of the Lace morph in all variables ( Figure S2.2). However, despite the availability and access to similar environmental conditions (similar environmental background), unimodal peaks of density profiles indicate that the Bell's morph tends to occupy areas with a broader temperature range, with lower precipitation, that are more arid, with higher potential evapotranspiration, and that receive more solar irradiation than that of areas occupied by the Lace morph.

| COUE framework
The topology of occurrence densities in PC environmental space ( Figure 4a) showed partial overlap between the niches of the Bell's versus the Lace morph. The niches of both morphs showed low to moderate overlap (D = 0.29; I = 0.51) yet were more similar (p > 0.05) than expected by chance, but not equivalent (p < 0.001; see Table S2.3). Furthermore, despite the availability and access to similar environmental conditions, the Bell's morph's niche slightly expanded (4.45%) and its centroid shifted towards environments that are more arid and receive more solar irradiation. Conversely, relative to the Lace morph's niche, the Bell's morph's niche unfilled (28.63%) from wetter and colder areas.

| n-Dimensional hypervolume framework
In corroboration with the findings using the COUE framework, the topology of the Bell's and Lace morph in n-dimensional hypervolume environmental space (Figure 4b) showed that the Bell's morph's niche is partially overlapping with that of the Lace morph's (75.22%).

| DISCUSS ION
Climate has been implicated as a primary or contributing factor influencing variation in the type and frequency of colour morphs across the geographical range of many ectothermic species (Friedman et al., 2017;McLean et al., 2015;Pérez i de Lanuza et al., 2018). The potential drivers of such variation, as we have here shown, can be meaningfully explored using the phenotypic data captured in georeferenced photographs obtained from citizen science platforms. In this case study of V. varius, we demonstrate that, albeit found across most of the species' geographical range, the Bell's morph occurs more frequently inland of the GDR, whereas the Lace morph occurs more frequently in relatively coastal regions. PPP analyses indicate that the distributions of both morphs are distinct and spatially segregated. The findings of univariate analysis of realised climatic niches, the two competing statistical frameworks used for multivariate analysis of realised climatic niches (i.e. COUE and n-dimensional hypervolume), and tests of niche conservatism (i.e. niche similarity and equivalency), all corroborated and showed similar yet distinct niches of the Bell's and Lace morph. The Bell's morph occupies environments that are more arid, have higher solar irradiation, a broader temperature range, higher potential evapotranspiration and lower precipitation compared to environments occupied by the F I G U R E 2 Distribution of two discrete colour morphs in the Lace monitor Varanus varius. Blue and red points indicate occurrence records of the Lace morph (n = 1537) and Bell's morph (n = 100), respectively. The digital elevation model illustrates the general coastal-inland stratification of the two morphs across high-elevation areas, particularly in South-Eastern Australia where the Great Dividing Range is broadest and highest. Photographs by Kristian Bell (Bell's morph) and Leo Berzins (Lace morph).
Lace morph. Overall, these findings support the hypothesised association between the Bell's morph's geographical distribution and climatic aridity, and supports previous research indicating that climatic variables such as temperature, precipitation and solar irradiation are correlated with geographical variation in colour polymorphism in ectothermic species Clusella-Trullas et al., 2008;Lattanzio & Buontempo, 2021;McLean et al., 2015;Pérez i de Lanuza et al., 2018).
Although our analytical framework is correlative in nature, our findings provide compelling insights into the climatic selective processes that might have led to such geographical variation in colour polymorphism, and the plausible ecological functions that ensue as a result of each morph. Colouration often has direct effects on thermoregulation and camouflage, both of which influence survival and hence evolutionary fitness (Cheng et al., 2018;Clusella-Trullas et al., 2008;Matthews et al., 2018;Olsson et al., 2013;Xing et al., 2016). Herein, we discuss the potential for climatic selection on thermoregulation and camouflage to generate the observed geographical distribution of colour morphs in V. varius.

| Potential adaptive functions of the Lace and Bell's morphs
Light skin colouration has higher thermal reflectance and hence slower rates of heat absorption (Geen & Johnston, 2014;Porter & Gates, 1969;Watt, 1968), thereby reducing the risk of overheating and enabling individuals to remain active under relatively hot ambient temperatures (Gibert et al., 1998;Kuyucu et al., 2018;Norris, 1967). Consequently, light skin colouration is often associated with areas experiencing high temperatures or high solar irradiation (Clusella-Trullas et al., 2008, 2009Goldenberg et al., 2021;Kuyucu et al., 2018). Similarly, our results show that the lighter Bell's morph occurs within an arid subset of the species' niche that experiences more incident solar irradiation. The Bell's morph may therefore constitute a thermoregulatory adaptation which is selected for in arid, highly radiative environments because it more effectively dissipates heat from direct and indirect solar irradiation, thereby preventing overheating (Gibert et al., 1998;Goldenberg et al., 2021).
The converse may be true of the Lace morph; its overall more melanic appearance, according to the thermal melanism hypothesis, would confer a faster heating rate in relatively coastal, shaded forest habitats, which are wetter and have lower solar irradiation than relatively open arid areas (Clusella-Trullas et al., 2008;Xing et al., 2016).
Thus, the divergence of colour morphs along a mesic-arid gradient may be driven by selection on thermal reflectance under different climates.
A complementary explanation to the evolution of these morphs being a direct physiological response to climate is an indirect response to climate-mediated differences in habitats across the species' range. Colour polymorphisms can emerge via selection on various forms of crypsis, such as background matching or disruptive colouration (Hall et al., 2013;Matthews et al., 2018;Price et al., 2019). In the case of the Bell's morph, its large and homogeneous areas of pale yellow on the body likely function in F I G U R E 4 Multivariate analysis of the realised climatic niche of the Bell's (red) and Lace (blue) colour morphs in the Lace monitor Varanus varius using (a) the Centroid shift, Overlap, Unfilling, and Expansion framework and (b) the n-dimensional hypervolume framework. The niches of both morphs are visualised in environmental space defined by seven climatic variables recorded in occurrence records and environmental background points. Principal component 1 (PC1) is negatively correlated (r > 0.9) to maximum temperature of warmest month and annual potential evapotranspiration and positively correlated to precipitation of driest quarter; PC2 is positively correlated (r > 0.9) to precipitation of wettest quarter; PC3 is primarily positively correlated (r > 0.4) to aridity; PC4 is primarily negatively correlated (r > 0.25) to precipitation of wettest quarter. See Figure S2.1 for more details on the correlation of the seven environmental variables to the first four principle component axes.
background matching, concealing the animal against the predominantly pale yellow/orange substrate colours that typify the visual environment in relatively arid climatic regions west of the GDR.
Similarly, selection for light body colouration in populations inhabiting pale sandy substrates has been demonstrated in numerous polymorphic species (e.g. Hoekstra et al., 2004;Rosenblum, 2006).
Moreover, the Bell's morph consists of bold, alternating light and dark bands, rather than being uniformly light. Such boldly contrasting bands can disguise body outlines, reducing predator attack success more effectively than colour static morphs (Cuthill et al., 2005;Hall et al., 2013;Price et al., 2019;Webster et al., 2013). The lighter Bell's morph may be poorly camouflaged and hence maladaptive in mesic (coastal) habitats (e.g. Cheng et al., 2018), which could additionally explain its absence/rarity in such environments. Mesic forest environments tend to offer darker backgrounds because they harbour dark soils rich in organic matter, and dense vegetation blocks sunlight, creating a well-shaded forest floor with dappled sunlight patches (Delhey, 2019;Endler, 1993). Darker colouration tends to prevail in these humid environments (i.e. Glodger's rule); hence, the Lace morph may be well concealed in mesic forest habitats, with its overall dark dorsal colouration and finely scattered light spots functioning to match the forests' dark background and obscure the body outline against the fine heterogeneity of sun and shade patches (Cheng et al., 2018;Delhey, 2019

| Geographical cline in morph frequency
Colour morphs can be mutually exclusive across a species' geographical range, but they often vary in local frequencies over environmental gradients (Fisher-Reid et al., 2013;Lotter & Scott, 1977;McLean et al., 2015;Pérez i de Lanuza et al., 2018;Schilthuizen, 2013 Figure 4). Indeed, differences in the topology of both morphs' realised climatic niches in environmental space support the observed geographical segregation; the Bell's morph is more prevalent in arid inland areas than mesic coastal areas, whereas the Lace morph exhibits the opposite occurrence pattern (Figures 2 and   3c). Field observations also support this geographical cline in relative morph frequency; for example, only the Lace morph occurs in mesic coastal regions near Sydney (Figure 2), whereas in sandy, semi-arid ecosystems at Nombinnie Nature Reserve (500 km inland of Sydney along a similar latitude), the ratio of Lace to Bell's specimens is approximately 6:4 (K. Bell, pers. comm.). Similarly, monomorphic Lace populations occur in mesic forests of far eastern Victoria (Figure 2), whereas the morph ratio is approximately 1:1 in populations on the inland side of the GDR in NE Victoria near Wangaratta (Heard & Black, 2003). Together, these findings highlight that mesic-arid transitions can be primary axes of environmental variation over which geographical variation in colour polymorphism emerges. In agreement, aridity has previously been linked to spatial variation in colour morph composition and frequency in the lizard Ctenophorus modestus (McLean et al., 2015).
Further field measurements of morph frequency at a range of sites spanning coastal-inland (i.e. mesic-arid) gradients are needed to elucidate the specific evolutionary process driving the geographical cline in morph frequency in V. varius.

| The utility of iNaturalist for studying spatial phenotypic variation
Georeferenced photographs, if used cautiously, can be profitably exploited by researchers wanting to explore the spatial occurrence of phenotypes in ways that traditional techniques cannot facilitate.
The potential for this technological advancement to further the field of biogeography becomes obvious when we consider whether the present study would be possible without such data; V. varius has a distribution spanning the entire eastern seaboard of Australia, and can be difficult to detect in an area without considerable search effort. One would likely need to implement a field survey program for a duration much longer than one's own lifetime to achieve the same coverage and intensity of phenotypic sampling as that of our iNaturalist data. Nonetheless, citizen science data are not without their biases (Dickinson et al., 2010;Isaac & Pocock, 2015), and investigators must assess the demands that their research question places on these citizen science data before ignoring data obtained by more traditional means (e.g. structured surveys with fixed sampling protocol).

| CON CLUS ION
Our study provides a clearer understanding of the geographical distribution of the colour polymorphism in V. varius. It is a true polymorphism in the sense that both morphs occur in sympatry across much of the species' range. However, there is also clear geographical variation in the polymorphism that corresponds to a mesic-arid climatic gradient, with the Bell's morph occupying a more arid climate niche than the Lace morph. There also appears to be monomorphic (Lace morph only) populations, particularly in coastal regions of the southeast. Based on these findings, we hypothesise that the Bell's morph may function as an adaptation to the arid environments with which it appears to be predominantly associated, possibly emerging as a direct evolutionary response to climate itself (e.g. to enhance thermoregulatory efficiency), or the need for individuals to match the visual environments found in different climatic regions (e.g. to enhance crypsis and avoid predation). Future research should directly measure changes in the relative frequency of morphs over coastal-inland (i.e. mesic-arid) gradients, and test experimentally how each morph influences fitness under varying ecological contexts. As for most taxa, we still lack basic information on phenotypic variation (sensu Hortal et al., 2015). We have demonstrated here how georeferenced photographs from citizen science databases can be used not only to overcome this knowledge shortfall, but also enable us to explore the geographical distributions of phenotypes without needing to conduct expensive, time-consuming field studies on wild organisms.

ACK N OWLED G EM ENTS
We thank those, too numerous to list, who submitted their photo-

CO N FLI C T O F I NTE R E S T
The authors declare that they do not have any conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The occurrence data, environmental variables and the R script of the