Genetic basis and expression of ventral colour in polymorphic common lizards

Colour is an important visual cue that can correlate with sex, behaviour, life history or ecological strategies, and has evolved divergently and convergently across animal lineages. Its genetic basis in non‐model organisms is rarely known, but such information is vital for determining the drivers and mechanisms of colour evolution. Leveraging genetic admixture in a rare contact zone between oviparous and viviparous common lizards (Zootoca vivipara), we show that females (N = 558) of the two otherwise morphologically indistinguishable reproductive modes differ in their ventral colouration (from pale to vibrant yellow) and intensity of melanic patterning. We find no association between female colouration and reproductive investment, and no evidence for selection on colour. Using a combination of genetic mapping and transcriptomic evidence, we identified two candidate genes associated with ventral colour differentiation, DGAT2 and PMEL. These are genes known to be involved in carotenoid metabolism and melanin synthesis respectively. Ventral melanic spots were associated with two genomic regions, including a SNP close to protein tyrosine phosphatase (PTP) genes. Using genome re‐sequencing data, our results show that fixed coding mutations in the candidate genes cannot account for differences in colouration. Taken together, our findings show that the evolution of ventral colouration and its associations across common lizard lineages is variable. A potential genetic mechanism explaining the flexibility of ventral colouration may be that colouration in common lizards, but also across squamates, is predominantly driven by regulatory genetic variation.

. Within species or groups of species, colour polymorphisms are more likely to lead to diversification and speciation compared to monomorphic colouration (Brock et al., 2021;Hugall & Stuart-Fox, 2012).Many species exhibit colour variants within populations, and these 'colour morphs' may be associated with differences in physiology (McKinnon & Pierotti, 2010), behaviour (Sinervo & Lively, 1996) and life-history traits (Roulin, 2004), often with a direct adaptive benefit (Maan & Sefc, 2013;Wollenberg Valero et al., 2019).However, associations between colour polymorphisms and other traits can be complex and even have detrimental effects.More specifically, theoretical work suggests that frequency-dependent polymorphisms (such as colour) can lead to a higher population extinction risk under environmental change (Bolton et al., 2015;Svensson & Connallon, 2019).Therefore, studying the mechanisms that give rise to colouration differences between individuals can provide crucial insights into understanding how diversity is generated and maintained, and may lead to adaptation and speciation.
Reptiles exhibit widespread colour variation (Olsson et al., 2013;Vitt & Caldwell, 2009), and several studies have identified colour polymorphism in reptiles as a heritable trait (Jin et al., 2022;Olsson et al., 2007;Sinervo et al., 2001;Vercken et al., 2007).Particularly within different groups of lizards (Andrade et al., 2019;Brock et al., 2021;Stuart-Fox et al., 2020), but also across the animal kingdom (Jamie & Meier, 2020), it has been suggested that ancestral variation in colouration genes may be responsible for the repeated evolution of colour polymorphisms.For example, repeated introgression across divergent lineages within a group leads to mosaic genomes (Yang et al., 2021), with potentially widely shared adaptive genomic regions, including colouration genes (Andrade et al., 2019).Among the lacertid lizards there is extraordinary intraspecific colour variation (Andrade et al., 2019;Brock et al., 2020;Stuart-Fox et al., 2020;Yang et al., 2021), and transitions from colour monomorphism to polymorphism in this group have been associated with a positive net diversification rate (Brock et al., 2021).
In reptiles, as across vertebrates, colouration is generated by pigment cells in the dermal layer of the skin, known as chromatophores.Reptiles have three types of chromatophores: melanophores contain the pigment eumelanin, responsible for black/ brown colouration; iridophores are pigment cells that contain light-reflecting structures, thus determining skin reflectance/ iridescence and xanthophores contain carotenoid and pteridine pigments, which impart yellow/orange or orange/red colouration respectively (Vitt & Caldwell, 2009).Synthesized pigment cells are stored in organelles such as carotenoid vesicles (carotenoids) and pterinosomes (pteridines).Each chromatophore has a specific structure, and its combined activity with other pigment cells results in distinct colouration patterns (Vitt & Caldwell, 2009).The synthesis of pigment cells is controlled, in part, by genetic mechanisms (Hubbard et al., 2010).However, few studies have attempted to identify the genetic mechanisms which allow variation in colour patterns to be expressed and transmitted to future generations (Olsson et al., 2013).Of these, most describe the mechanisms of melanin pigmentation (Olsson et al., 2013), with great focus on MC1R, a highly conserved gene across vertebrates that is responsible for the production of melanin (Rosenblum et al., 2004).
However, genes responsible for other pigmentation patterns remain largely uncharacterized.
The common lizard (Zootoca vivipara) is the most widespread terrestrial squamate, covering a large part of Eurasia spanning from west (Spain) to east (Japan) (Roitberg et al., 2013;Takeuchi et al., 2013).The species is split into six distinct evolutionary lineages that vary in geographic distribution and include two alternative reproductive modes, being either egg-laying (oviparous) or live-bearing (viviparous) (Recknagel et al., 2018;Surget-Groba et al., 2006).Common lizards vary considerably in their ventral colouration, from white to bright orange (Fitze et al., 2014;Vercken et al., 2007Vercken et al., , 2010)), and their colour polymorphism is associated with differences in behaviour (Vercken & Clobert, 2008), immunity (Vercken et al., 2008) and reproductive strategies (Vercken et al., 2007).For example, in one lineage (western viviparous), three female ventral colour morphs can be distinguished, each associated with a different level of aggression: the yellow morph being aggressive, the orange being submissive and the mixed morph (intermediate colour or both yellow and orange scales) displaying either of the two behaviours depending on the context (Vercken & Clobert, 2008).Maintenance of all three colour morphs was suggested by an overdominance effect of the mixed morph (Vercken & Clobert, 2008).In another lineage (western oviparous), males display three different ventral colour morphs, and these are associated with behaviours that follow a pattern of negative frequency-dependent selection (Fitze et al., 2014;Sinervo et al., 2007).Therefore, while the type of colour-behaviour associations are variable and both sex-and lineage-dependent (Fitze et al., 2014), ventral colour variation exists across the species.
Understanding the genetic mechanisms of colouration in the common lizard, thus, provides an excellent opportunity towards resolving a detailed understanding of colour polymorphism and its associations with other traits.
Here, we report on female ventral colour variation in two common lizard lineages and explore its associations with reproductive traits and molecular differences.We do so by leveraging an admixture zone where the alternative colour/parity lineages hybridize and by characterizing colour morph variation from a molecular to cellular to organismal level.First, we quantified and described the colour variation of oviparous and viviparous females, as well as hybrids.Using data on clutch and offspring size, we evaluated the extent to which colouration is associated with female reproductive investment, given that in this species colour has been found to be associated with reproductive behaviours.Second, we used genomic and transcriptomic approaches to explore the genetic basis of colour, both in regions of genomic association and in functional gene expression.With this information, we finally evaluated population genomic patterns and signals of selection in candidate colouration genes identified here and from the literature.

| Sampling, husbandry and reproductive traits
Collection and husbandry followed Recknagel and Elmer (2019).

| Ventral colouration measurements and statistical analyses
Females were photographed using a Canon 70D and a 60-mm fixed lens under standard settings (fixed ISO: 160, aperture: 8, exposure compensation: 0) and with an X-rite ColorChecker to standardize colour between images.X-rite ColorChecker software (vers 1.0.1) was used to calibrate colour of the original images prior to analysis.
The colour profile was applied to the image and white balance was selected using the off-white square on the ColorChecker board.RGB (red, green, blue) values for the selected square were set to 230 nm using the eyedropper tool in Adobe Photoshop CC 2017.Five representative ventral scales were selected from the lower two-thirds of the belly using the lasso tool, excluding any melanic spots and glare.
The colour from the five selected areas was averaged with the RGB values recorded and these were converted to hue, saturation and lightness (HSL) values.
Repeatability of the colour measure over time was inferred from a subset of five females.Variation within the season was inferred from images taken once each week for 11 weeks, and variation within a day was inferred from images taken twice (morning and afternoon) on the day of measurement.
Absolute number of ventral melanic spots was counted on each female in Adobe Photoshop CC 2017.In addition, proportion of ventral area covered by melanic spots (proportion of spots) was calculated as follows: the area of melanic spots was selected using the magic wand tool, the selected area was extracted and divided by the total area of the ventral body (from vent to neck, excluding limbs).
Statistical analyses of colour were carried out in R (R Core team, 2015).A PCA was conducted in R with the 'prcomp' package to derive one major axis (PC1) summarizing HSL colour variation.For count data (number of spots), a generalized linear model (glm) with quasi-Poisson distribution was performed, followed by an analysis of variance (ANOVA).A beta regression was fitted to the proportional data (proportion of spots).Colour variables hue, saturation and lightness were fitted using a glm and ANOVA.
Model-based clustering using the R package 'mclust' was performed to classify colour PC scores, excluding hybrid individuals.We tested the most likely number of clusters in the data and if these were predicted by parity mode.
Generalized linear models (GLMs) were performed to test for an association between colour, parity mode and reproductive investment.HSL and the combinatorial effect of those variables in the PCA (PC1) were considered as the response variables, with parity mode (measured as number of days that it took offspring to emerge

| Admixture mapping
Genetic data corresponding to 496 of the 558 females phenotyped for ventral colouration, reproductive mode and reproductive investment were extracted from NCBI (PRJNA657575) as raw ddRADSeq data.These had been generated in a prior study on viviparous and oviparous common lizards (Recknagel et al., 2021).Data were first processed using STACKS (Catchen et al., 2013) to remove reads with uncalled bases and low quality (options -c and -q).Next, filtered reads of each individual were mapped to the common lizard reference genome (Yurchenko et al., 2020), deposited with NCBI (GCA_011800845.1,UG_Zviv_1, PRJNA638687).Genotypes were extracted from STACKS with the following conditions: each SNP had to be present in at least 66% of all individuals, an individual coverage of >5× and a minimum allele frequency of 10%.The resulting vcf file was phased and missing data were imputed using BEAGLE (Browning & Browning, 2007).Duplicate individuals (same individuals sampled across multiple years) were identified and removed using KING (Manichaikul et al., 2010); all phenotypic analyses were conducted excluding duplicate individuals (total N without duplicates: ventral colouration = 445, melanic spots = 414).
The proportion (Q) of each female's genome deriving from the oviparous or viviparous lineage was reconstructed based on all SNPs and k = 2 using ADMIXTURE (Alexander et al., 2009) For a single SNP showing the highest association with ventral colour (Log p-value: 7.78 on chromosome 5 at position 6,100,408 bp), we tested for the correlation between allelic state (homozygous A/A, homozygous G/G, heterozygous A/G) and hue.
PCAdapt (Luu et al., 2017) was used to infer signals of selection between the two parity modes from the genotypic data, which identifies divergently selected SNPs between two groups while correcting for background noise.The selection scores generated for each SNP (Benjamini-Hochberg corrected p-values) were then compared to the p-values obtained from admixture mapping to test whether any SNPs associated with colour were under divergent selection.
In addition, we tested if SNPs significantly associated with ventral colouration were linked to or within regions associated with parity mode in the common lizard.Genetic regions associated with parity were extracted from the literature (Recknagel et al., 2021).We considered colour-associated SNP as linked if they were within ~100 kb of a region associated with parity mode, following Recknagel et al. (2021).

| Differential expression analysis by RNA-seq
Seven ventral tissue samples were selected for RNA extraction: four pale oviparous (mean hue = 29.67 ± 0.88 std.error) and three yellow viviparous (mean hue = 42.00 ± 2 std.error) females.Lizards were sacrificed (as per the approved Schedule 1 method) and ventral scales from the lower two-thirds of the belly were excised.RNA was extracted using the PureLink RNA Mini kit (Life Technologies, Carlsbad, CA) following the manufacturer's instructions, with an additional homogenisation prior to extraction.After quality control, one sample from each colour was removed.RNAseq library preparation and sequencing were carried out by NovoGene on a NovaSeq using 150 bp paired-end chemistry.
After filtering of reads containing adapters and removing low-quality sequences with Trimmomatic version 0.39 (Bolger et al., 2014), all reads were mapped to the common lizard reference genome (Yurchenko et al., 2020) using STAR v2.7.2b_0928 (Dobin et al., 2013) with default parameters.Gene expression was quantified with htseq v0.11.2 (Anders et al., 2015) using the function htseqcount (parameters: --stranded = no, --order = pos, --type = CDS and idattr = gene_id).Genes with <20 reads across all samples were removed and counts normalized via log2 transformation with the R package DESeq2 v3.543 (Love et al., 2014).Differential expression analyses were carried out with DESeq2 and visualized using the R package EnhancedVolcano v1.0.1 (Blighe et al., 2020), applying a q-value (p-value adjusted for discovery of false positives) cut-off of <.05.A principal component analysis of the sequenced individuals was performed with the pca function and the method "svdImpute" in the R package pcaMethods (Stacklies et al., 2007).An analysis of gene ontology term enrichment was carried out using PANTHER v16.0 (Mi et al., 2021) GO-slim with the Anolis carolinensis reference set.
Through a review of the literature, genes with previous associations with melanin (N = 7) and pteridine (N = 8) pigment synthesis pathways (Braasch et al., 2007) and carotenoid metabolism (N = 25) (Waagmeester et al., 2009) were identified and exhaustively checked for differential expression between colour morphs in our dataset.Iridophore-related gene expression was investigated using a set of genes previously found to be enriched in iridophores (N = 409) (Higdon et al., 2013).

| Structural assessment of skin cell pigmentation
To identify the presence of cellular structures associated with colour, we conducted transmission electron microscopy (TEM) on ventral skin from one representative pale (oviparous) and yellow (viviparous) female used in the RNAseq.Skin samples were first fixed in 2.5% glutaraldehyde/0.1 M sodium cacodylate buffer, then dehydrated by ethanol series, stained and embedded in resin.Ultrathin sections (60-70 mm) were produced, additionally stained with 2% methanolic uranyl acetate and Reynolds lead citrate and then viewed on a JEOL1200EX transmission electron microscope at an accelerating voltage of 80 kV.Images were captured using an OLYMPUS ITEM soft imaging system.
Chromatophores were identified with the aid of reference images (Andrade et al., 2019).

| Identification of candidate genes
Candidate genes for ventral colouration and melanic spots were identified by first examining top SNPs (Benjamini-Hochberg corrected p-value <.05).These SNPs were assigned to candidate regions if at least two significantly associated SNPs were within a distance of 5 Mb.This conservative wider cut-off (e.g.Pallares et al., 2014;Pollinger et al., 2005) was chosen under the assumption that a relatively simple genetic basis underlies colouration with few loci of large effect, and because our marker density was only moderate, likely missing causal variants (Santure & Garant, 2018).The top SNP within a candidate region was first examined in more detail by identifying the nearest gene using bedtools (Quinlan, 2014).Next, the candidate gene was inspected for differential expression between yellow and pale females and its involvement in pigmentation pathways.If the nearest gene to the top SNP did not reveal any differentiation in expression and had no role in pigmentation pathways, we widened the search to the other genes within the 5 Mb candidate region.Any genes in a candidate region that were differentially expressed and showed an association with a pigmentation pathway were also considered candidate genes.

| Population genomics of colour candidates DGAT2 and PMEL
To investigate variation in colouration candidate genes, wholegenome sequence data of oviparous (N = 6) and viviparous (N = 16) females generated during a prior study (Recknagel et al., 2021) were extracted from NCBI (PRJNA657575).Candidate genes of interest were extracted from 5′ UTR to 3′ UTR using the common lizard reference genome annotation (Yurchenko et al., 2020).Allele frequency differences (AFD) and F ST (with 100 bp windows using a 50 bp step) were calculated using vcftools (Danecek et al., 2011).
Variants with F ST in the 97.5th percentile of the distribution were identified as outliers.Characterization of exons, introns and UTR regions was based on the common lizard reference genome annotation.S1).

| Ventral colour variation between lizards of different parity modes
Colour measurements were stable over time (Table S2).All three   S1).

| Reproductive investment is not associated with ventral colour
Reproductive investment was not associated with ventral colour variation after accounting for parity mode.More specifically, after model selection from generalized multiple regressions, only a significant association between parity mode and colour remained in the model (Hue: t-value = −14.90,p < .0001;saturation: t-value = −14.24,p < .0001;lightness: t-value = 12.73, p < .0001;PC1: t-value = −16.66,p < .0001;Table S3), while other variables such as clutch size and mass and offspring mass and survival were not associated with colour.

| Genetics of ventral colouration
Admixture analysis of 83,665 genome-referenced single nucleotide polymorphisms (SNPs) derived from 445 colour-phenotyped females revealed that 16.4% of individuals had genomes composed of at least 1% contribution of more than a single lineage (N oviparous = 169, N viviparous = 203, N hybrids = 73; Figure 2A).Using admixture mapping, we found 16 SNPs significantly associated with ventral colouration (Benjamini-Hochberg corrected p-value <.05) and most of these were within two genomic regions (Figure 2B).The most prominent candidate region (2.5 Mb) was on chromosome 5 and held six significantly associated SNPs, including the SNP with the strongest association (−log 10 p).The closest gene (~4 kb) to the SNP with the strongest association was DGAT2.Out of the top 28 associated SNPs genome wide (above −log 10 p > 4 cut-off), 11 (39%) were found to be associated with this candidate region (2.5 Mb) on chromosome 5 (Figure 2B).
Another candidate region was found on chromosome 2, spanning approximately 4 Mb (Figure 2C).This region had five SNPs showing an association with colour (p < .0001;Table S4) and included 86 genes.The top SNP fell within the COPRS gene, which is associated with gene regulation (Lacroix et al., 2008) and has no apparent role in pigmentation.
Finally, the second most significantly associated SNP overall was an isolated SNP located on chromosome 3 (Figure 2A), with the closest gene to this SNP being PLCB1.
We found two candidate regions (four SNPs) associated with the number and size of melanic spots on the female ventral region (Figure 2D).The most prominent region and the most strongly associated SNP were found on chromosome 1 (three significantly associated SNPs).The two closest genes downstream to the top SNP were the protein tyrosine phosphatases PTPRJ (~32 kb) and PTPN5 (~173 kb).An isolated SNP was found on chromosome 12, with the two closest candidate genes (LPP, ~19 kb; TPRG1, ~71 kb) not exhibiting an apparent role in pigmentation.

| Differential gene expression of ventral colouration
From RNAseq of ventral skin, a total of 1798 genes were differentially expressed between pale and yellow females (Figure 3A).

Gene ontology (GO) analysis revealed 54 significantly enriched GO
terms in the ventral skin (Table S5), and these were mostly related to metabolic and cellular processes.Background differences in ventral skin gene expression were expected to be high and potentially confounded due to the ca. 4 million year evolutionary divergence between the two colour/parity lineages (Cornetti et al., 2014); therefore, we focused our analysis on the set of known colour genes we had extracted from the literature (Table S6).
We found that all seven of the candidate genes involved in the melanin synthesis pathway (Braasch et al., 2007;McLean et al., 2017) were expressed in the ventral skin tissue.Three of these were significantly differentially expressed -TYR, TYRP1 and SLC24A5 -and upregulated in the yellow (viviparous) individuals (Figure 3b; Table S6).
Another melanin synthesis gene, PMEL, was marginally significant (q = 0.056).This gene was located within the candidate region for colouration on chromosome 2 that was identified in the genetic mapping.
All 25 of the known carotenoid metabolism genes were expressed in the ventral tissue (Table S6).Five of these genes (ALDH18A1, DGAT1, DGAT2, DHRS3, SCARB1) were found to be significantly differentially expressed between colour morphs (Figure 3C); four genes were up-regulated in yellow relative to pale females, while DHRS3 was up-regulated in the pale skin.The differentially expressed carotenoid genes included DGAT2 (q < 0.05), which was the gene closest to the most significantly associated SNP on chromosome 5 identified in the genetic mapping.
Six known pteridine synthesis genes (from a total set of eight genes) were expressed in the ventral skin but none were differentially expressed between yellow and pale (Figure 3d; Table S6).
Additionally, 115 of the genes known to be associated with iridophores were found to be expressed in ventral tissue, including 15 genes that were differentially expressed between the two colour morphs, approximately evenly divided between up-regulated in pale and up-regulated in yellow (Table S7).

| Molecular variation in candidate colouration genes
Using whole-genome sequences, we examined the two candidate genes that were both differentially expressed in RNAseq and strongly associated with colouration in the genetic mapping, DGAT2 and PMEL.We found no fixed differences either in the coding sequence, intronic regions or UTR regions of the genes (Figure 4A-C).
DGAT2 and PMEL allele frequency differentiation between colour morphs was relatively low compared to other differentially expressed pigmentation genes (Figure 4A).Close investigation of DGAT2 and PMEL candidate gene sequences revealed that variants in exonic regions had generally lower differentiation between colour morphs and were never F ST outliers (Figure 4B,C).Further, we found no evidence that the genomic regions containing the colourassociated loci were under selection (Figure S3), as inferred from PC analysis of outlier loci.Rather, the majority of SNPs significantly different between colour morphs (13 out of 15 in DGAT2, and all 7 in PMEL) were located in introns (Figure 4B,C).Two of the outliers in the DGAT2 sequence were located in the 3′ UTR region.
We found that colour candidate genes DGAT2 and PMEL were not more closely linked to previously identified parity mode genes (Recknagel et al., 2021) than expected by chance (across all SNPs, 57.4% are closer to a parity mode region than DGAT2, and 53.6% are closer than PMEL; Figure S4).Similarly, candidate genes for melanic spots are not linked to previously identified parity mode candidate genes (64.3% of SNPs are closer to a parity mode region than PTPRJ) (Figure S4).
From the RADseq dataset, analysis of the most highly associated SNP close to the DGAT2 gene showed that females with pale ventral colour were always homozygous for one allele (genotype G/G), while hybrid and yellow females had all three genotypes (G/G, G/A, A/A) (Figure 4D).Nonetheless, individuals that had the A allele had a higher hue value, meaning they tended to be more yellow (F (2,442) = 36.0,p < .0001).This was the case for both heterozygous and homozygous individuals, which did not differ significantly (Figure 4D).

F I G U R E 3
Differential expression in melanin-and carotenoid-based pigmentation between colour morphs.(A) Volcano plot depicting differential gene expression between pale (oviparous) and yellow (viviparous) ventral skin.One gene, MATN1, exhibited an extreme log 2 FC value (−25.8) and was excluded for graphical reasons.Labels indicate pigment pathway genes that show significant differential expression between colour morphs.Horizontal dashed line: q-value = 0.05.Vertical dashed lines: log 2 fold-change = 2 or −2.Light grey: Non-significant data points; Orange: higher expression in yellow viviparous individuals; Dark grey: higher expression in pale oviparous individuals.(B) Mean log 2 -transformed counts of genes involved in melanin synthesis.(C) Mean log 2 -transformed counts of genes involved in carotenoid metabolism.(D) Mean log 2 -transformed counts of genes involved in pteridine synthesis.Significant differences in gene expression are denoted by asterisks (*q < 0.05; **q < 0.01; ***q < 0.001).

| DISCUSS ION
Common lizards represent a well-known example for the role of ventral colour polymorphisms in ecological and evolutionary strategies (Fitze et al., 2014;San-Jose et al., 2014;Sinervo et al., 2007;Vercken et al., 2007Vercken et al., , 2008Vercken et al., , 2010)).Here, we show that oviparous and viviparous common lizards differ significantly in their ventral colouration within a contact zone: oviparous females have pale bellies covered with melanic spots, while viviparous females have yellow bellies with fewer or no spots (Figure 1).Hybrid individuals were intermediate in both colouration and the number of melanic spots.
While ventral colour polymorphisms have been described in other common lizard lineages, ours is a first evaluation of two cooccurring but evolutionarily distinct viviparous and oviparous F I G U R E 4 Evolutionary history and divergence of pigmentation genes.(A) Allele frequency difference (AFD) of genetic mapping candidate genes and differentially expressed genes between oviparous and viviparous individuals derived from whole genome sequences.Shown are differences in coding regions, introns and UTRs.The first two genes were the candidates derived from genetic mapping (DGAT2 and PMEL), and the other eight genes were differentially expressed.Box limits show the first quartile (25th percentile) and third quartile (75th percentile).The whiskers extend to 1.5× inter-quartile range from the box limit.(B) Mean F ST between oviparous and viviparous individuals within the DGAT2 sequence.(C) Mean F ST between oviparous and viviparous individuals within the PMEL sequence.The outlier threshold was calculated as the 97.5th percentile of the distribution separately for each gene.Fifteen outliers were identified in DGAT2 and seven outliers were identified in PMEL.Gene sequences are displayed above the graphs, with exons, introns and UTR regions shown by chromosomal position.Exons are coloured in orange, introns are coloured in blue and UTR regions are coloured in purple.(D) Frequency of genotypes of the SNP showing the highest association with colouration and close to the DGAT2 gene.Females with the A allele had a higher hue value and were more yellowish than individuals homozygous for the G allele.Violin plot and box plots with the median, lower and upper quartiles and outliers are shown.This suggests that common lizard colour variation and its genetic links with other reproductive traits such as strategies and behaviour are evolutionarily unstable.
A key to understanding the evolution of colour polymorphisms, their adaptive significance and correlations with other traits is to identify the genetic mechanism and basis of colouration.Our admixture mapping approach and differential expression data identify a strong colour candidate gene: DGAT2.The SNP showing the highest association with yellow ventral colouration is closest to this gene, and it is also differentially expressed in ventral scales of the two parity modes differing in colouration.Vertebrates cannot synthesize carotenoids endogenously and must acquire them through their diet (von Lintig, 2010).During the metabolism of dietary carotenoids, lipids are involved in their transport and delivery to the chromatophores.Our gene ontology (GO) analysis revealed 54 significantly enriched GO terms in the ventral skin (Table S5), mostly related to metabolic and cellular processes, including lipid biosynthesis, an important process for carotenoid metabolism (Toews et al., 2017).
The DGAT2 gene encodes a transmembrane protein that promotes the synthesis of triacylglycerol droplets (Bhatt-Wessel et al., 2018;Chitraju et al., 2019), where dietary carotenoids are stored.The protein is involved in the cellular uptake of carotenoids and is associated with carotenoid-based colouration in lizards (de Mello et al., 2021;McLean et al., 2017), birds (Toomey et al., 2017) and fish (Ahi et al., 2020).A direct link between expression differences in DGAT2 and presence and absence of carotenoid esters in yellow and white cichlid skin (Ahi et al., 2020) suggests that DGAT2 might be involved in the esterification of carotenoids, a process that may affect the deposition and distribution of these pigments in the skin, and thereby influence colouration.The most significantly associated SNP from our admixture mapping that lies close to the DGAT2 gene suggests a dominant allele effect for yellow colouration.We found no fixed differences between pale oviparous and yellow viviparous individuals in the coding region, and together with the differential expression of DGAT2, this may suggest that regulatory variation rather than fixed coding differences play a role in regulating the colour polymorphism.
A second candidate gene for ventral colouration that we identified is PMEL, located within a region of 4 Mb associated with ventral colour in our admixture analysis.Further, PMEL was marginally significant in the differential expression analysis (Figure 3b).PMEL is a transmembrane glycoprotein involved in early melanogenesis and production of eumelanin in melanocytes (Hellström et al., 2011;Watt et al., 2013).More specifically, PMEL regulates the formation of fibrils in the melanosome (Watt et al., 2013).In later stages of melanogenesis, these fibrillar sheets are packed with melanin, leading to black or brown colouration.In vertebrates, PMEL is associated with pigmentation-related mutants that can result in yellow colouration, particularly in epidermal tissue such as skin (Burgon et al., 2020), scales (Wang et al., 2022), and feathers (Ishishita et al., 2018;Zheng et al., 2020).For example, in golden tilapia, CRISPR knockouts of PMEL genes lead to weak vertical stripes and a yellowish colouration (Wang et al., 2022).
Finally, the SNP with the second largest −log 10 p score was nearest to PLCB1.This gene is also related to colour, as it is part of the protein kinase C pathway and is involved in melanocyte differentiation (Vidács et al., 2022).A study on common North American woodpeckers found that PLCB1 and PLCB4 were associated with yellow and orange plumage colouration (Aguillon et al., 2021).Overall, support for PMEL and PLCB1 genes as candidates is less clear than for DGAT2 but highly suggestive and worth further research.
In addition to identifying candidate genes for ventral colouration, our admixture mapping of ventral melanic spots revealed another interesting candidate gene class: protein tyrosine phosphatases (PTPs).These proteins have been previously shown to be involved in pigmentation phenotypes (Baxter et al., 2019).More specifically, the gene PTPN11 is responsible for the LEOPARD syndrome in humans, a phenotype that includes the formation of multiple lentigines (melanic spots) (Digilio et al., 2002;Legius et al., 2002;Motegi et al., 2015), and mutations in the PTPN6 gene are involved in a range of abnormal skin pigmentation phenotypes in mice (Green & Shultz, 1975;Nesterovitch et al., 2011).In common lizards, the extent of black spotting on the ventral scales is sexually dimorphic, with females having fewer (Vroonen et al., 2013).In males, it has been associated with fitness traits (San-Jose et al., 2017), while its role in females remains currently unclear.As far as we are aware, this is the first genetic mapping to identify candidate genes for melanic spots in a squamate.Given its potential covariation with fitness, it is ripe for further exploration within and across species.
Across common lizards and through extensive behavioural ecology research, ventral colour has been thought to vary with reproductive strategies (Fitze et al., 2014;Vercken et al., 2007Vercken et al., , 2010;;Vercken & Clobert, 2008).Our findings suggest that colouration is not associated with, nor genetically linked to, parity mode or reproductive investment.Considering that linkage is already weak (below an R 2 of .15) at a distance of around 100 kb (Recknagel et al., 2021), it is unlikely that proximity to parity mode genes allowed for genetic hitchhiking of DGAT2 (closest to parity mode gene at ~11.3 Mb; Figure S4a) or PMEL (closest to parity mode at ~7.3 Mb; Figure S4b).Combined evidence from the absence of a signal of selection and a lack of association with variation in reproductive investment suggests that ventral colouration is not a trait under strong selection in females in these common lizard lineages.Instead, we suggest that colour is genetically independent of parity mode and reproductive investment.Future behavioural experiments co-varying colour and parity mode would be valuable to infer if there is lineage-based mate preference for ventral colour.
Colour polymorphism has frequently been implicated as a trait facilitating species diversification (Brock et al., 2021;Hugall & Stuart-Fox, 2012).For example, colouration is extremely important in mate choice as a visual cue and can quickly lead to assortative mating and divergent adaptation (McKinnon & Pierotti, 2010;Seehausen & Van Alphen, 1998).Our results and published data on colouration of other squamate species suggest that colouration evolves flexibly and convergently.In the wall lizard, a lacertid lizard species closely related to the common lizard, it was found that entirely different genes (BCO2, SPR) regulate differences in orange/ yellow ventral colouration driven by pteridines and/or carotenoids (Andrade et al., 2019).In a broader evolutionary context, we suggest that colouration polymorphisms evolve quickly and through diverse genetic mechanisms (Hofreiter & Schöneberg, 2010;Keene et al., 2015;Kratochwil, 2019).
Colouration genes are usually best known for their involvement in biological pathways leading to the production of pigments.
However, most of these genes have pleiotropic roles and are crucial components of metabolic pathways during ontogeny and the maintenance of organisms (Hoekstra, 2006;Hubbard et al., 2010), and loss-of-function mutations can have lethal consequences (Nadeau et al., 2008).Therefore, it is likely that non-coding, regulatory regions play an important role in colour variation (Galbraith & Hayward, 2023;Kratochwil et al., 2018;Orteu & Jiggins, 2020;Toomey et al., 2017;van der Burg et al., 2020).In common lizards, we showed that there are no fixed differences in protein-coding regions in the set of colour candidate genes between the pale-bellied oviparous and yellow-bellied viviparous lineages that we identified.Indeed, this would be compatible with research showing that squamate reptiles have the highest proportion of intraspecific cisregulatory mutations regulating pigmentation (Elkin et al., 2023).

| CON CLUS ION
Here, via genetic mapping and differential expression analyses, we identified genomic regions and candidate genes that show signals consistent with regulating ventral colour variation between syntopic oviparous and viviparous common lizards of pale versus yellow ventral colouration: the DGAT2 gene, and a larger region containing the PMEL gene as a potential candidate.Moreover, we identified a genomic region associated with ventral melanic spots, which is potentially under control of a PTP candidate gene.The ventral colour and melanic spot genes warrant further exploration to validate their function and infer their frequencies and phenotype associations more broadly in these and other common lizard lineages.
Contrary to other common lizard lineages, reproductive investment is not associated with ventral colouration in the oviparous-viviparous contact zone studied here.The degree of colour variation in common lizards and its association with behavioural syndromes has been questioned and remains contentious in the scientific community (Cote et al., 2008;Vercken et al., 2008).Therefore, the consequences of ventral colouration appear to vary across lineages and must be further addressed to get a clearer understanding of genetic correlations and their evolutionary significance in common lizards.
Our findings provide valuable candidate gene information and genomic patterns that can be explored in other populations and settings.Finally, our findings suggest that the evolution of differential colouration by regulatory changes may be a mechanism by which behavioural syndromes and other colouration-life-history trait correlations can co-evolve quickly, but this warrants further investigation.

[
FPPM]), snout-vent length (SVL) was measured, a tail-tip tissue sample was taken and females were released at site of capture.Their clutches were counted (clutch size [CS]) and weighed (clutch mass [CM]) using a digital balance on the morning after parition and incubated at 24°C until hatching.Each hatchling was weighed (offspring mass [OM]) before being released to the site of their mother's capture.From these data, we derived the following measures to describe reproductive investment: average egg mass (CM/CS = EM), relative clutch mass (CM/FPPM = RCM), relative offspring mass ([CS*HM]/ FPPM = ROM) and hatchling survival (N hatched offspring/CS = OS).
from eggs) and reproductive investment (clutch size [CS], offspring mass [OM], offspring survival [OS], relative clutch mass [RCM], relative offspring mass [ROM] and average egg mas [EM]) as explanatory variables and size (SVL) as a covariate.Model simplification was used, removing the least likely explanatory variable (based on p-value) successively until only significant results remained.
. To estimate standard errors in Q, 200 bootstrap replicates were performed.Individuals in which the standard error did not cross 1 or 0 were counted as a hybrid (McFarlane & Pemberton, 2019), while the remaining individuals were counted as either pure oviparous or viviparous.Multivariate linear mixed models (MVLMMs) in GEMMA (Zhou & Stephens, 2012) were used to map (i) colouration (using the HSL measurements; N = 445) and (ii) melanic spots (measured as number of spots and proportion of spotted area; N = 414).A centred relationship matrix was estimated to correct for individual relatedness and population stratification.Statistical significance was assessed by adjusting p-values with the Benjamini-Hochberg correction.

F
I G U R E 1 Phenotypic divergence in colouration between oviparous, hybrid and viviparous females.(A) Representative pale oviparous (up), hybrids (middle) and brightly yellow-coloured viviparous female (down).(B) Boxplots illustrating differences in ventral colouration, measured as (B) hue, (C) saturation and (D) lightness and the (E) number of ventral melanic spots and (F) proportion of ventral melanic spots.Hybrids are intermediate between the two parity modes; all differences are significant between parity modes.Letters a-c, denote significant differences.Box plots show the median, lower and upper quartile and outliers.(G) Violin plot showing differences between parity modes in PC2 of the ventral colouration PCA.(H) Ventral colouration PCA (derived from hue, saturation, lightness and number of spots) showing differences and overlap between the parity modes and hybrids.
the chromatophore layers (xanthophores, iridophores and melanophores) were observed at the cellular level in both parity modes (Figure S2).Oviparous females had more melanic spots and these occupied a larger area on their ventral region relative to viviparous females (number of spots (ANOVA): F (2,411) = 60.01,p < .0001;proportion of spots (beta regression): z = 11.61,R 2 = .14,N = 413, p < .0001)(Figure 1E,F).Hybrid individuals (N = 73, confirmed with genomic admixture analysis) were intermediate both in colouration and in number of melanic spots.Differences between hybrid and oviparous or viviparous colouration are significant in all but one comparison (Figure 1B-F; Table

F
Genetics of ventral colouration in oviparous and viviparous common lizards.(A) Admixture graph based on 445 female individuals from the contact zone of hybridizing pale oviparous and yellow-bellied viviparous common lizards (to facilitate visualization, standard error bars are not shown).(B) Genetic mapping results of the ventral colouration phenotype show a peak on the start of chromosome 5, with the strongest SNP closest to the gene DGAT2.(C) Closer view of chromosome 5 and the candidate region with the colour locus.Regions associated with parity mode are shown in red.The closest distance between colour loci and parity mode loci is not less than 10 Mb, which is beyond the distance of LD decay(Recknagel et al., 2021) and therefore not expected to be in strong linkage.The candidate region (~5 Mbp) is indicated by a yellow vertical bar.Density of SNPs is shown in blue along the top.(D) The candidate region on chromosome 2 contains the colouration candidate gene PMEL and is more than 7 Mb apart from the closest parity mode genes and is also not expected to be linked.The candidate region (~5 Mbp) is indicated by a yellow vertical bar.Density of SNPs is shown in blue along the top.(E) Genetic mapping of the ventral melanic spots illustrating a peak on chromosome 1.The closest genes to the top SNP are PTPRJ and PTPN5, two different protein tyrosine phosphatases and candidate genes for pigmentation phenotypes.
colour morph.In contrast to previous research within the Western Viviparous lineage(Vercken et al., 2007(Vercken et al., , 2010)), we found no association between female ventral colouration and female reproductive investment (based on measures of clutch size, weight and viability; N = 558) once accounting for differences in parity mode.This suggests that female ventral colour is not linked to reproductive investment.We conclude that while parity mode (being oviparous or viviparous) has been shown to be associated with reproductive investment(McLennan et al., 2019; Recknagel &   Elmer, 2019), ventral colour does not play a role in these differences.