Phylogenetic and morphological influence on habitat choice in moisture‐harvesting horned lizards (Phrynosoma spp.)

Abstract In previous studies, the superhydrophilic skin of moisture‐harvesting lizards has been linked to the morphological traits of the lizards’ integument, that is, the occurrence of honeycomb‐shaped microstructures. Interestingly, these structures can also cover the skin of lizards inhabiting wet habitats. We therefore tested the influence of the microstructures’ main features on the habitat choice and wettability in the genus Phrynosoma. The genus Phrynosoma comprises moisture‐harvesting species as well as nonspecialists. Lizards of this genus inhabit large areas of North America with diverse climatic conditions. Remarkably, the differences in the manifestation of microstructures are just as versatile as their surroundings. The phylogeny of the lizards as well as the depth of their ventral microstructures, though independent of each other, correlated with the precipitation in their respective habitat. All other morphological traits, as well as the skin's wettability itself, could not predict the habitat of Phrynosoma species. Hence, it is unlikely that the microstructure influences the wettability, at least directly. Hence, we presume an indirect influence for the following reasons: (a) As the ventral side cannot get wet by rain, but the belly could easily interact with a wet surface, the microstructure might facilitate water absorption from wet soil following precipitation. (b) We found the number of dorsal microstructures to be linked to the occurrence of silt in the habitat. In our study, we observed scales being heavily contaminated, most likely with a mixture of dead skin (after shedding) and silt. As many lizards burrow themselves or even shovel sand onto their backs, deploying the substrate might be a mechanism to increase the skin's wettability.


| INTRODUC TI ON
High temperatures as well as rare and irregular precipitation characterize arid regions across the world. To survive in such a harsh environment, animals require specialized adaptations, for example, avoiding of overheating and gathering enough water. Efficient water collection from rain, moist air, or soil has been described as "moisture harvesting" and evolved independently in several species (Cardwell, 2006;Lillywhite & Stein, 1987;Norgaard & Dacke, 2010;Repp & Schuett, 2008;Sherbrooke et al., 2007;Yenmiş et al., 2015).
Moisture-harvesting lizards, as one example for these species, show a stereotypical behavior: In several species of the genera Phrynocephalus, Trapelus, and Phrynosoma, this behavior includes (a) lowering of head and tail, (b) raising of abdomen in an arch, (c) splaying of legs, and (d) flattening of the body (Schwenk & Greene, 1987;Sherbrooke, 1990Sherbrooke, , 2004Vesely & Modry, 2002). On the other hand, the integument (skin and appendages like scales) of such lizards shows specific adaptations as well. Whereas in wet habitats, like rainforests, the lizards' scales exhibit superhydrophobic and selfcleaning properties to prevent infections (Joel & Weissbach, 2019;Watson et al., 2015), in arid regions they have to prevent water evaporation from their body and simultaneously absorb, retain, and/or facilitate drinking of water (Comanns, 2018;Joel et al., 2018).
The reptilian epidermis in general consists of an inner α-and an outer β-layer, mainly composed of α-and β-keratins, respectively. These two layers are divided by the mesos layer (Alibardi & Maderson, 2003;Maderson et al., 1998). Because of its composition of lipids, it is responsible for the water impermeability of the lizard's integument and hence protects the lizard from desiccation. The outer β-layer is covered by the corneous protective Oberhäutchen (Nick, 2014). For harvesting water, the lizards' scales overlap, forming an interconnected channel network (Comanns et al., 2015;Gans et al., 1982;Sherbrooke et al., 2007). Because of capillary forces within these channels, water transportation takes place passively and sometimes even with a pre-determined flow direction toward the mouth (Comanns et al., 2011(Comanns et al., , 2015(Comanns et al., , 2016Yenmiş et al., 2015). In addition to this channel network, the lizards must have a hydrophilic integument. The hydrophilicity was so far linked to a honeycomb-shaped microstructure (or microornamentation) on the Oberhäutchen, keeping a film of water after pre-wetting (Berthé et al., 2009;Comanns et al., 2011;Hofling & Renous, 2009;Riedel et al., 2015;Sherbrooke et al., 2007;Vesely & Modry, 2002;Yenmiş et al., 2015). This microstructure, however, is also found on nonmoisture-harvesting reptiles, such as the modest forest dragon Hypsilurus modestus (own observation) or the smooth helmeted iguana Corytophanes cristatus (Lang, 1989).
Our study focuses on the genus Phrynosoma (horned lizards), which comprises at least three moisture-harvesting species (P. cornutum, P. modestum, and P. platyrhinos) as well as at least one nonspecialist (P. hernandesi) (Dodge, 1938;Peterson, 1998;Sherbrooke, 1990Sherbrooke, , 2002. Additionally, members are endemic to southwestern regions of North America and inhabit mainly, but not only, arid habitats (Baur & Montanucci, 1998;Bryson et al., 2012;Lara-Resendiz et al., 2014, 2015Newbold & MacMahon, 2014;Sherbrooke, 2002;Zamudio & Parra-Olea, 2000). In case the microstructure does indeed entail increased wettability, we assume that Phrynosoma species inhabiting drier areas exhibit more and/or deeper honeycomb-shaped structures on their Oberhäutchen. Additionally, we assume that closely related species would inhabit similar habitats, as their morphological features would be beneficial for colonization. Hence, in this study we analyze the linkage between microstructure, phylogeny, habitat choice as well as wettability of Phrynosoma spp.

| Study animals
Six captive-bred Phrynosoma platyrhinos were kept in an arid terrarium (146 cm × 55 cm × 80 cm). A soil mixture consisting of sand and clay was used as substrate. Below the heat lamp, temperature was 45-50°C, whereas the cooler regions of the terrarium were 28°C.
The lizards were fed with ants, cowpea weevils, and microcrickets.
Please refer to Baur and Montanucci (1998) for further information about keeping Phrynosoma.
Additional to samples of our individuals, shed skin of captivebred P. asio, P. cerroense, P. coronatum, P. goodei, P. hernandesi, P. orbiculare, P. modestum, and P. taurus was kindly provided by private breeders. Please note that we were not able to match these samples to individuals within one species, and thus, there might be pseudoreplicates included in the shed skin data.
Museum specimens of wild-caught P. asio, P. cerroense, P. cornutum, P. hernandesi, P. orbiculare, P. platyrhinos, P. solare, P. taurus, and Moloch horridus were obtained from the collection of the Zoological Research Museum Alexander Koenig, Bonn, Germany. The used biological material is also listed in the supplement (Table S1).

| Wetting experiments
We refrained from using museum specimens for wettability tests, as we observed ruptures of the Oberhäutchen in several specimens, leading to an absorption of the applied fluid by the subjacent skin layers. Hence, we chose to use only shed skin for wetting experiments. Shed skin samples exhibiting a rupture within the region of interest were discarded.
Three pieces of shed skin, ventrally and dorsally, were taken for examination. Scale structure differs between ventral and dorsal sides in Phrynosomatidae, making discrimination easy. The pieces of shed skin were placed onto a glass slide, and 2 µl of distilled water with 0.6% blue food dye (Queen Fine Foods PTY Ltd., Alderley, Australia) for coloration was applied. Spreading of the liquid was observed via a high-speed recording microscope (250 fps, VW-9000C; Keyence Cooperation, Osaka, Japan). The differences (%) of the area covered by the droplet (0 to 60 s) were evaluated with the help of a customized Python (3.5) script. For calculation of mean and statistics, we replaced data for perfect wetting (i.e., droplet wets complete skin within 60 s and thus a rim is unidentifiable) with the largest value we were still able to measure in our data set (raw data: Table S7).
We validated our data with the wettability of living animals (P. asio, P. modestum, P. platyrhinos, P. solare, P. taurus). For this, a droplet of tap water was applied to the dorsal side of the animal to visually check whether spreading occurs ("yes" vs. "no"). For all but P. platyrhinos, these data were provided by private breeders.

| Cleaning of samples
All skin samples were contaminated to varying degrees. This was not only true for museum specimens but also for shed skin samples.
Hence, the contamination is not (only) due to storage and/or age of the museum samples. Different cleaning methods were exploited to remove the contamination on the shed skin of P. platyrhinos, without any success: We first used an ultrasonic bath with different cleaning agents (distilled water, 70% ethanol, and chloroform) for 30 s and 60 s. In a second approach, shed skin samples were incubated and swirled for 90 min in chloroform and subsequently treated for 30 s in chloroform in the ultrasonic bath. The same treatment was carried out with soapy water and the final cleaning with distilled water. Alternatively, the shed skin pieces were either fixed with double-sided sticky tape and cleaned with a fine brush and/or compressed air of 10 bar, or a second stripe of adhesive tape was slightly pressed onto the distal side of the shed skin and then slowly pulled off. In a last trial, the shed skin was fixed upside down to a centrifugation tube, either empty or filled with distilled water. The samples were then centrifuged with 2,000 g for 2 min.

| Morphological examination of microstructures
Skin samples (approx. 5 × 10 mm, ventrally and dorsally) were taken with a scalpel from preserved museum specimens of P. asio, P. cerroense, P. cornutum, P. hernandesi, P. orbiculare, P. platyrhinos, P. solare, and P. taurus. Samples were dehydrated by an ascending ethanol series. Remaining alcohol was washed out with hexamethyldisilazane and dried overnight. Afterward, one side of the sample was clean cut with a razor blade. Slices were either placed flat or with the cutting edge facing up on SEM plates. Afterward, samples were sputtered with gold and observed with a scanning electron microscope (Philips AG, Amsterdam, the Netherlands). Shed skin samples of P. asio, P. cerroense, P. coronatum, P. goodei, P. hernandesi, P. modestum, P. orbiculare, P. platyrhinos, and P. taurus were prepared as described above without the dehydration step.
Data are presented as mean ± SD, where SD is the standard deviation (raw data: Table S7, including biological material source). SEM samples were analyzed using Keyence VW-9000 Software (version 1.4.0.0). The clean cut samples were used to measure the depth of the microstructure (with n as high as possible for each sample), whereas the flat samples were used to calculate the microstructures per mm 2 (with three spots per sample).

| Histology
Histological examinations were performed with skin samples of museum specimens of P. solare, P. cornutum, and M. horridus. As museum specimens, they were wild-caught animals and not bred for several generations on a non-natural substrate. Thus, the contamination should appropriately reflect any contamination occurring in the wild.
Samples were prepared and dehydrated as described before and then cut with razor blades into slices of ca. 1 mm width. These slices were embedded in epoxy resin. Semi-thin sections of ~0.7 µm thickness were cut with an OM U3 microtome (Reichert, Vienna, Austria).
For the subsequent staining, the epoxy resin was removed as described by Mayor et al. (1961). Prepared sections of each lizard and body side were stained with either methylene blue staining, alizarin red staining, periodic acid-Schiff staining, or gram staining (Mulisch & Welsch, 2015). Cover-slips with sections, stained with methylene blue, alizarin red, or periodic acid-Schiff stain, were glued to an object slide using Merckoglas. For gram-stained sections, Canada balsam was used for final fixation. All sections were examined with a BA310 digital microscope (Motic Deutschland GmbH, Wetzlar, Germany).

| Phylogenetic tree
Depending on the taxonomic concept, the genus Phrynosoma currently comprises 17-21 species (Montanucci, 2015;Uetz et al., 2020;Zamudio et al., 1997). However, by testing several species delimitation methods, Blair and Bryson (2017) found no support for the classification proposed by Montanucci (2015), who resurrected the names P. brevirostris and P. ornatissimum and used morphological information to describe two new taxa (P. bauri and P. diminutum).
Therefore, we follow the taxonomy proposed by Blair and Bryson (2017) within the present study.
Phylogenetic relationships within Phrynosoma were inferred with BEAST version 1.8.3 (Drummond et al., 2012), which estimates phylogeny using Bayesian MCMC. We used previously published data (Hodges & Zamudio, 2004;Leaché & McGuire, 2006;Nieto-Montes de Oca et al., 2014) of six nuclear loci (BDNF, EXPH5, NKTR, R35, RAG1, and SOCS5, Tables S2 and S3). Substitution models were calculated using ModelTest (Posada & Crandall, 1998) as implemented in the package phangorn for Cran R. Model parameters were estimated separately for each locus, but loci were not partitioned by codon positions in order to avoid overparameterization that impeded the MCMC mixing. The Yule speciation process was set as tree prior to using a random starting tree. Results of the BEAST analyses were summarized from four independent runs with 2 × 10 7 generations each, sampling every 2 × 10 3 trees and omitting the initial 10% as burn-in after checking for convergence and sufficient effective sample sizes with Tracer v1.7 (Rambaut et al., 2018).

| Hydrological and pedological niche assessment
With a spatial resolution of 30 arc seconds, climatic information characterizing humidity and precipitation conditions within the range of the studied Phrynosoma species was obtained from the CHELSA database representing an average between 1979 and 2013 (Karger et al., 2017). As we were predominantly interested in the hydric niche of the species, we downloaded monthly mean, minimum, and maximum temperature as well as precipitation and computed bioclimatic variables as listed in Table 1 using the ENVIREM package for R (Title & Bemmeis, 2017). These variables describe annual averages, extremes, and different aspects of the conditions the species experience within their habitats. Therefore, these parameters are likely to represent the selective landscape influencing the morphological adaptations associated with moisture harvesting. Soil data for the average conditions within a depth range of 0-5 cm and with a spatial resolution of 250 m were obtained from the SoilGrids database (https://www.isric.org/explo re/soilg rids).
As geographic range surrogates, georeferenced records for each species were obtained from the global biodiversity information facility (www.gbif.org) and the VertNet database ( (Pagel, 1999) and Blomberg´s K (Blomberg et al., 2003) using the phylosignal package for R (Keck et al., 2016).
Based on each set of records, the median position of the species in climate space was computed and median traits of all three morphological variables were tested in a phylogenetic general least-squares (PGLS) framework in R using a Brownian motion model as phylogenetic constrain (Orme et al., 2018). We acknowledge the possibility of alpha error accumulation using false discovery rate (fdr) adjustment (p.adjust function in R). We performed a principal component analysis using 10,000 random samples throughout the study area considering those variables, where we detected a significant relationship between morphological features and environmental parameters irrespective of potential alpha errors (see results). Based on the median positions of each species in the resulting principal components, we reconstructed ancestral conditions using the "fastBM" function and trait diagrams using the "phenogram" function of the phytools package for R (Revell, 2012). These trait diagrams illustrate a projection of the phylogenetic tree in a space defined by the climatic niche and time (Evans et al., 2009). The niche breadth of each species was illustrated using density profiles (function "ggdensity," "ggpubr," and "ggplot2" packages for R; Wickham (2016) and Kassambra (2020)). The multidimensional niche position and volume were quantified by computing hypervolumes based on support vector machines (hypervolume package for R; Blonder, 2015;Blonder et al., 2014).

| Phylogenetic tree and divergence times
Our dated phylogeny of Phrynosoma obtained from the BEAST analyses is generally well supported, with most nodes having posterior probabilities >.95 (Figure 1). The inferred relationships largely reflect the subgeneric divisions of Phrynosoma and the topology is in concordance to previously published phylogenies of the genus (Leaché & Linkem, 2015;Leaché & McGuire, 2006;Wiens et al., 2013).

| Phylogenetic signal in hydric niche dimensions
Based on phylogeny, we analyzed the dispersion of Phrynosoma across America. We indeed found strong evidence that the habitat choice of Phrynosoma is in dependence on its phylogenetic origin.
There were significant phylogenetic signals analyzing hydric as well as the pedological niche dimensions (Table 2; i.e., Thornthwaite aridity index, annual mean precipitation (bio_12), minimum (bio_14) and maximum (bio_13) precipitation of the driest/wettest month, precipitation of the wettest quarter (bio_16), the climatic moisture index, and the proportion of soil organic carbon compounds (soc) and clay).
As suggested by both Pagel's lambda and Blomberg's K, a Brownian motion model of evolution may explain best the relationships in hydric niche dimensions. Phylogenetic signals in pedological niche dimensions were less pronounced.

| Characterization of the morphological traits in correlation to phylogeny
According to the phylogenetic signal when analyzing the niche dimensions, we assumed that likewise morphological traits correlate to habitat choice, as morphological traits are likely to be similar in closely related species. We therefore tested the linkage between phylogeny and microstructure characteristics. Phrynosoma's Oberhäutchen are covered with pentagonal or hexagonal microstructures ( Figure 2). The surface of the scales was not uniformly structured, though, but microstructures were pronounced in the periphery projecting into the channels between the scales. A microstructure was often missing in the middle of the scales, especially on scales of P. cornutum and on ventral scales of almost all species (to a different degree) except for P. taurus.
On the scales, ~ 1,300 to ~ 2,800 structures per mm 2 were counted, whereas the mean depth varied from 0 to 7.5 µm (Figure 1, Table S5).

| Correlation of morphology and habitat
Though morphological traits are independent of phylogeny, habitat choice was linked to morphological traits additionally to being linked to phylogenetic relations. By testing if morphological features can be explained by climatic factors, the phylogenetic generalized least squares regressions suggested a significant correlation in 9 of 175 tests when treating tests separately (Table 3, Appendix S1, but see potential caveats below). Annual mean precipitation (bio_12), precipitation of the wettest month as well as quarter (bio_13, bio_16), precipitation seasonality (bio_15), potential evapotranspiration (PET) of the coldest quarter, and PET seasonality were significant predictors of ventral microstructure depth (six cases, Figure 3) and the interaction of ventral and dorsal microstructure depth (two cases). Only one pedological parameter was significantly correlated with morphology (silt ~number of dorsal microstructures).
As our data analyses were explorative, we need to acknowledge the potential of alpha error accumulation, and indeed, accounting for false discovery rates (fdr), none of our correlations remained significant. Testing 175 correlations at an alpha of 0.05 results in an adjusted significance level of 0.00028, which is very unlikely to be met by any correlation considering our sample size and natural variation in both morphological and ecological space. Therefore, our results need to be interpreted with caution.

| Comparisons of hydric niches among taxa
A principal components (PC) analysis based on those habitat variables significantly correlated with morphological features suggested two PCs, which in total represent 84.7% of the total variation. PC1 (63.1%, eigenvalue = 3.8) is strongly negatively correlated with PET seasonality and positively correlated with annual mean precipitation (bio_12), precipitation of the wettest month (bio_13), and quarter (bio_16). PC2 (21.7%, eigenvalue = 1.3) is mainly driven by precipitation seasonality (bio_15) and PET of the coldest quarter (cf. correlation circle in Figure 4).
Comparing the distribution of species in the environmental space revealed strongly varying degrees of niche position and breadth ( Figure 5). The niche space as suggested by the hypervolume F I G U R E 1 Dated phylogeny of the genus Phrynosoma (left), traits measured (middle) as well as distribution records used in this study (right). The tree was obtained from a Bayesian relaxed-clock analysis of six nuclear loci. The outgroup (Sceloporus bicanthalis) is not shown for clarity. Gray bars show 95% highest posterior density (HPD) of divergence times, numbers on nodes refer to Table S4, where the detailed values are given. Solid nodes have posterior probabilities >.95 and the topology presented here is identical to the fully supported phylogeny based on phylogenomic data of Leaché and Linkem (2015). Next to the phylogenetic tree are the mean values for wettability (wet) as well as microstructure number (MSn)

| Wettability of Phrynosoma
Foreshadowing a link between morphology of the skin and habitat choice, we finally tested whether surface structures indeed increase wettability. We detected large differences in wettability of Phrynosoma species: Neither all species nor both body sides of one species (dorsal or ventral) were equally well wettable (Figures 1 and 6, Table S6). Only the shed skin of four of the nine tested species (P. asio, P. cerroense, P. hernandesi, and P. modestum) showed a prompt spreading of water on their skin (i.e., perfect wetting), and only for P. modestum, this was true for both body sides ( Figure 6, Table S6). This pattern of wettability could be replicated in live animals (only dorsal side): P. asio, P. modestum, and P. solare (no shed skin samples available for the latter species) proved to be perfectly wettable, whereas water did not spread well on the skin of P. platyrhinos and P. taurus. The differences in wettability were neither following an obvious phylogenetic nor morphological pattern. We additionally observed that the stereotypic moisture-harvesting behavior does not have to indicate a good wettability of the skin: Even though P. platyrhinos performed the stereotypical behavior (Figure 6a), the wettability of its skin was rather low compared to other Phrynosoma species.

| Contamination as an influencing factor?
The "non-predictability" of wettability by neither morphological trait nor phylogeny would stay in contrast to the described correlation between some climatic factors and morphological traits. However, applying fdr corrections advise us to interpret some of our data with caution. During SEM analysis, we encountered problems with contaminants on the scales of almost all Phrynosoma species, varying in shape and size as well as quantity even within one individual (Figures 2 and 7). None of our cleaning approaches resulted in successful removal of these contaminants. We hypothesized contamination might be a concomitant feature of the microstructure and could be the link between morphology, phylogeny, and hydric and pedological niches.
For a first insight into such an interplay, we analyzed the origin of the contaminations on museum skin samples of P. cornutum, for which moisture harvesting has been described in literature, as well as P. solare, which has been observed to be perfectly wettable (tests in live animals). Additionally, we included the moisture-harvesting lizard Moloch horridus in our study. M. horridus is often used as a comparison species for Phrynosoma, as it inhabits a similar habitat in Australia, but developed moisture-harvesting independently of Phrynosoma (Meyers & Herrel, 2005;Pianka & Parker, 1975;Pianka & Pianka, 1970;Sherbrooke, 1990Sherbrooke, , 1993Sherbrooke et al., 2007).
All three species exhibited a high amount of contamination on their skin.
Semi-thin sections of resin-embedded skin of P. cornutum, P. solare, and M. horridus were stained with methylene blue. We observed the contamination to be stained more strongly than the skin and purple instead of blue (Figure 8). Alizarin red dyes inorganic particles, while, for example, skin layers stay unstained. We observed stained particles within the contamination, often measuring only a TA B L E 3 Summary statistics of significant phylogenetic generalized least-square regressions. P-values were adjusted using the FDR method of p.adjust in R. Significance is marked in red. For abbreviations, see as also observed in SEM images (Figure 7), were heavily stained pink.

| D ISCUSS I ON
Previous studies linked the superhydrophilic skin of moistureharvesting lizards to the occurrence of honeycomb-shaped microstructures. Our results show a correlation of precipitation in the habitat with the ventral microstructure depth as well as a correlation of hydric as well as the pedological niche dimensions with F I G U R E 3 Significant correlations between hydric niche dimensions and ventral microstructure depth as revealed by pgls analyses. None of these correlations remained significant when correcting for alpha error accumulation using false discovery rates. Abbreviations: bio_12 annual mean temperature; bio_13 precipitation of the wettest month; bio_15 precipitation seasonality; bio_16 precipitation of the wettest quarter; PET potential evapotranspiration. Species are abbreviated with the first three letters of their specific epithet to control for potential alpha errors revealed that none of these correlations remained significant. We conclude that it is unlikely that the microstructures influence the wettability, at least directly. We speculate that Phrynosoma and other moisture-harvesting lizards could deploy the soil of their habitat to influence the skin's wettability. Future studies have to validate this hypothesis. Wiens et al. (2013) suggested that Phrynosoma originated in comparatively arid environments and mesic habitats were occupied only later. Adaptation to arid environments may therefore be the ancestral condition. According to our dated phylogeny, diversification within the genus mainly happened 15-46 Ma ago ( Figure 1) and within this period, climatic conditions in Northern and Central

| Evolution of environmental niches
America have been subject to numerous fluctuations and changes (Chapin, 2008). Interestingly, we found a phylogenetic signal for the choice of habitat, regarding precipitation-related variables and two of three pedological factors. This is partly in concert with previous findings: Analyzing both temperature and precipitation-related variables, Luxbacher and Knouft (2009)  In contrast, P. goodei and P. platyrhinos are closely related species with limited habitat overlap. In accordance with the prediction that phylogeny influences niche choice, their habitats showed, for example, similar annual precipitation ( Figure 3). Nevertheless, aside from the dorsal depth of microstructures, both are morphologically distinct and equally poorly wettable. On the contrary, the next related species, P. modestum, has no habitat overlap with both and is superhydrophilic.
It has equally flat ventral microstructures compared to P. platyrhinos and similar quantity of dorsal microstructures compared to P. goodei.
We found correlations for both these features with hydric (ventral depth), respectively, pedological (dorsal quantity) niche dimensions. Though we found correlations between morphological traits and niche, applying fdr corrections to control for potential alpha errors erased all correlations. Hence, our data have to be interpreted with caution and we assume it is rather unlikely that the microstructures (directly) influence wettability. Skin structures have been shown to influence tribological properties (Baum et al., 2014), which could be reflected in worn down microstructures at exposed parts of the scales in Phrynosoma.

| Moisture harvesting in Phrynosoma spp
Phrynosoma is a typical exemplary genus for moisture (or rain) harvesting lizards, and for three species, the stereotypical behavior has already been described. For P. modestum, our study could confirm that moisture harvesting implies a superhydrophilic skin. However, for P. platyrhinos, we found the contrary: the stereotypical moistureharvesting behavior was accompanied by a badly wettable skin with only few and flat microstructures. Interestingly, P. hernandesi was formerly observed to conventionally drink when thirsty (Dodge, 1938), but we characterized their skin as superhydrophilic. Future studies hence should also focus on the important feature of the wettability of the skin when describing moisture harvesting in reptiles.
Four of the five species with the largest niche space characterization show a superhydrophilic skin (P. cornutum is described as moisture harvesting having a predominant flow direction of spreading water (Comanns et al., 2015)

| Origin of sample contamination
It was surprising to find an indicated link between deeper ventral microstructures and higher precipitation, as the ventral side of Phrynosoma is not exposed to rain. On the other hand, higher precipitation also leads to wetter soil, which has regular contact to the ventral integument. Analyzing the relation between Phrynosoma's morphological traits and soil parameters revealed that the occurrence of silt (particle size: ≥0.002 and ≤0.05 mm) in the habitat correlates with the number of dorsally located microstructures. Particles of the contamination on our samples were stained by alizarin red and thus contain calcium (Mulisch & Welsch, 2015 Pianka & Pianka, 1970). There was no phylogenetic signal for silt, but other soil parameters in the habitat (clay and organic particles; Table 2). Hence, we speculate a link between silt and microstructures could have developed independently.
Additionally, there was a strong signal for shed skin, especially remnants of the lower α-layer. An accumulation of fragments of old skin seems plausible due to the lizards' innate shedding behavior and these skin remnants would be stained by gram staining due to the amount of cystine in α-keratin (Alibardi, 1998;Fischer, 1953). The organic contamination on the P. solare scales (both sides) also consisted of fungi, seen in the SEM as well as in our histological examinations. Fungi are found on several healthy reptiles as biofilm (Paré et al., 2003). We found fungi only on P. solare, and hence, assume this contamination was due to a captured unhealthy animal or incorrect storage of samples and has no functionality.

| Influence of the honeycomb-shaped microstructure on wettability
If the microstructure of moisture-harvesting lizards' scales would be the key factor for the efficient spreading of water on their skin, being clean would be essential. However, Phrynosoma were never directly observed nor described in literature to perform grooming or similar behavior. On the contrary, Phrynosoma as well as M. horridus evolved habits that increase their level of contamination: Depending on the species, behaviors like rubbing their venters on the ground, shoveling sand on their dorsal surface, or burrowing themselves in a sandy and dusty substrate were described (Comanns et al., 2016;Pianka et al., 1998;Sherbrooke, 1993;Weese, 1919). All of these can contribute to silt accumulating in the honeycomb-shaped microstructures of the lizards' scales. Additionally, the lizards' skin often sheds in pieces, especially the α-layer in direct contact with the new forming Oberhäutchen (Maderson et al., 1998). These residual skin pieces complement the contamination. Therefore, we hypothesize that the scales' microstructure does not contribute to the lizards' wettability directly, but can serve as container for contamination.
Smaller particles, like shed skin and especially silt, would influence the wettability positively due to (a) smaller capillaries forming between the grains and particles and (b) in case of silt: reduced contact angle compared to the skin's keratin due to different material chemistry. Further investigations should evaluate this by, for example, comparing the wettability of freshly shed animals (no silt accumulation) to animals exposed to silt in a controlled manner.

ACK N OWLED G M ENTS
We

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The supplement includes the following tables: S1 Samples used for investigation. S2 Taxon sampling and GenBank accession numbers.