Genetic , Demographic and Ecological Factors Contributing to the Evolution of Anolis Lizard Populations

Geographic and environmental factors can be catalysts of evolutionary processes and central drivers of genetic and phenotypic differentiation. Characteristic patterns of genetic and phenotypic variation elucidate whether predominately adaptive or nonadaptive processes are involved in shaping variation among populations. Such genetic and phenotypic signatures can provide insight in how specific extrinsic factors relate to adaptive and non-adaptive processes. This dissertation aims to investigate whether and how disturbance events, intraspecific competition and secondary contact after a species introduction contribute to genetic and phenotypic differentiation among populations of Anolis lizard. In chapter 1, I tested whether patterns of genetic diversity and differentiation of island populations in the Bahamas were consistent with disturbance-driven demographic fluctuations and strong genetic drift. I used ten microsatellite markers to measure genetic diversity and population differentiation. Population genetic structure and differentiation were consistent with expectations of strong genetic drift after founder events. Genetic bottlenecks were prevalent across the island populations and coalescent-based demographic models supported a colonization scenario for most islands, rather than scenarios for population bottlenecks or relative stability of populations. Low rates of gene flow among island populations have likely preserved the genetic signatures of colonization and extinction events. This study provides evidence that disturbance events such as hurricanes are catalysts of long-term, nonadaptive evolution when gene flow is limited. In chapter 2, I examined whether intraspecific competition is a potential driver of directional morphological change in male and female lizards (Anolis sagrei) among island populations in the Bahamas. I used measures of population density and injuries, as a proxy for aggressive encounters, to test whether body size and head size vary in the direction predicted by natural selection. I found that lizards in high-density populations had a higher proportion of injuries as compared to low-density populations, suggesting that the former experience higher levels of intraspecific aggression. Additionally, lizards in these higher density populations had larger heads (longer and wider, corrected for body size). In both sexes, relative head size and frequency of injury increased with population density, suggesting that females may play an active role in intraspecific competition. In chapter 3, I examined whether the introduction of the Cuban green anole (Anolis porcatus) to South Florida has resulted in local hybridization with its closely related sister species, the native green anole (Anolis carolinensis), or whether reproductive barriers prevent gene exchange across species boundaries. I used one mtDNA marker and 18 microsatellite loci to test for cyto-nuclear discordance, which is indicative of hybridization, and whether population-genetic patterns are consistent with recent or historic gene flow. I found mtDNA haplotypes of two species in South Miami, local A. carolinensis and A. porcatus from West Cuba. Contrary to expectations of recent hybridization, I did not find intermediate nuclear genotypes. Instead, the population in South Miami forms a separate, genetically homogeneous cluster, which is differentiated from both parental species. A historic gene flow analysis confirms that ~33% of the nuclear ancestry of the South Miami population is derived from West Cuban A. porcatus. Thus, reproductive isolation between A. porcatus and A. carolinensis is weak or absent, despite considerable divergence time in allopatry, which reinforces a proposal to revise the taxonomy of A. carolinensis and A. porcatus from West Cuba according to guidelines of the biological species concept. In summary, my dissertation evaluated how demographic and ecological factors relate to patterns of genetic and phenotypic variation in populations of Anolis lizards. My research shows that hurricanes can act as evolutionary catalysts by altering nonadaptive evolutionary processes when gene flow is limited. The same island populations vary morphologically in the direction predicted by natural selection as a consequence of intraspecific competition. Lastly, historic hybridization resulting from secondary contact between formerly allopatric introduced and native species can locally erode species boundaries and alter the genetic composition of native populations. By characterizing patterns of genetic and phenotypic variation, my research provides insight in to how extrinsic factors and long-term evolutionary processes are interconnected and potentially shape future trajectories of these

In chapter 2, I examined whether intraspecific competition is a potential driver of directional morphological change in male and female lizards (Anolis sagrei) among island populations in the Bahamas. I used measures of population density and injuries, as a proxy for aggressive encounters, to test whether body size and head size vary in the direction predicted by natural selection. I found that lizards in high-density populations had a higher proportion of injuries as compared to low-density populations, suggesting that the former experience higher levels of intraspecific aggression. Additionally, lizards in these higher density populations had larger heads (longer and wider, corrected for body size). In both sexes, relative head size and frequency of injury increased with population density, suggesting that females may play an active role in intraspecific competition.
In chapter 3, I examined whether the introduction of the Cuban green anole (Anolis porcatus) to South Florida has resulted in local hybridization with its closely related sister species, the native green anole (Anolis carolinensis), or whether reproductive barriers prevent gene exchange across species boundaries. I used one mtDNA marker and 18 microsatellite loci to test for cyto-nuclear discordance, which is indicative of hybridization, and whether population-genetic patterns are consistent with recent or historic gene flow. I found mtDNA haplotypes of two species in South Miami, local A. carolinensis and A. porcatus from West Cuba. Contrary to expectations of recent hybridization, I did not find intermediate nuclear genotypes.
Instead, the population in South Miami forms a separate, genetically homogeneous cluster, which is differentiated from both parental species. A historic gene flow analysis confirms that ~33% of the nuclear ancestry of the South Miami population is derived from West Cuban A. porcatus. Thus, reproductive isolation between A. porcatus and A. carolinensis is weak or absent, despite considerable divergence time in allopatry, which reinforces a proposal to revise the taxonomy of A. carolinensis and A. porcatus from West Cuba according to guidelines of the biological species concept.
In summary, my dissertation evaluated how demographic and ecological factors relate to patterns of genetic and phenotypic variation in populations of Anolis lizards.
My research shows that hurricanes can act as evolutionary catalysts by altering nonadaptive evolutionary processes when gene flow is limited. The same island populations vary morphologically in the direction predicted by natural selection as a consequence of intraspecific competition. Lastly, historic hybridization resulting from secondary contact between formerly allopatric introduced and native species can locally erode species boundaries and alter the genetic composition of native populations. By characterizing patterns of genetic and phenotypic variation, my research provides insight in to how extrinsic factors and long-term evolutionary processes are interconnected and potentially shape future trajectories of these populations.
x LIST OF TABLES

Introduction
The frequency and intensity of catastrophic disturbance events such as fires, floods and hurricanes is increasing as a result of global change (Aponte et al. 2016;Arnell & Gosling 2016;Bender et al. 2010;Goldenberg et al. 2001;Westerling et al. 2006). Disturbance events can catalyze ecological change by rapidly restructuring ecosystems and species communities (Turner 2010). Such disturbance-driven ecological changes are relevant for evolution when they alter environmental conditions, meta-community structure, and population sizes (Banks et al. 2013;Hanski 2012;Pelletier et al. 2009;Schoener 2011). Disturbance can shift the direction and strength of natural selection and consequently drive adaptive evolution (Grant & Grant 2002;Grant et al. 2017). Alternatively, changes in genetic diversity following population size reductions and altered gene-flow patterns (Apodaca et al. 2013;Spear et al. 2012) can increase the relative importance of genetic drift as a driver of nonadaptive evolution (Banks et al. 2013;Davies et al. 2016). However, empirical studies that highlight the relative importance of disturbance events as non-adaptive driver are rare, because genetic data are often limited to a small number of populations or lack undisturbed reference populations (Apodaca et al. 2013;Fleming & Murray 2009), making it difficult to distinguish the effects of population size reductions from preexisting population structure and post-disturbance gene flow (Banks et al. 2013).
We took advantage of a replicated small-island system in the Bahamas with a wellunderstood disturbance regime and used genetic data to study the non-adaptive evolutionary consequences of reoccurring colonization and extinction events. We include a subset of islands with known colonization history as reference populations and compared patterns of genetic diversity and differentiation of experimental islands to islands with natural colonization histories.
The Bahamas are located in a major hurricane trackway and previous work has shown that hurricanes are an important driver of extinction-recolonization dynamics in populations of the brown anole, Anolis sagrei (Losos et al. 2003;Schoener et al. 2001aSchoener et al. , 2004Spiller et al. 1998;Spiller & Schoener 2007). Hurricane-induced extinction and subsequent recolonization follows a complex pattern in which island size, elevation, and timing of reproduction determine population survival and natural recolonization (Schoener et al. 2001a(Schoener et al. , 2004 In this study, we systematically investigated whether and how repeated population-size reductions, extinctions, and recolonizations have affected the genetic diversity and differentiation of 17 island populations of brown anoles near Staniel Cay. We hypothesized that hurricane-induced extinction and natural recolonization is frequent in this system and that strong genetic drift is a central driver of population diversity and differentiation. We include four reference islands, that were experimentally colonized in 1977 with 5-10 lizards from Staniel Cay, a larger source island for an unrelated study (Schoener & Schoener 1983); the remainder are all thought to have been occupied continuously since before 1977. Because the timing of colonization and the number of colonizers is known, these islands serve as valuable comparisons for our analysis decreasing the ambiguity associated with inferring demographic processes from patterns. We tested whether patterns of genetic diversity and differentiation across these 17 islands were consistent with disturbance-driven population demography. First, we used overall genetic diversity, population structure, and differentiation to test for signatures of strong genetic drift. Second, we used demographic models to reconstruct the potential processes leading to the observed patterns of genetic diversity. Finally, we compared population-genetic parameters and demographic models of the reference populations to naturally established populations with unknown demographic histories to gain insight in whether similar demographic processes gave rise to similar patterns of genetic diversity and differentiation.
In the case of strong drift, we predict low genetic diversity, random genetic differentiation and high inbreeding coefficients. In contrast, if an alternative process, such as vicariance due to rising sea levels has occurred, we predict genetic differentiation to be related to geographic distance among islands (i.e., isolation by distance) and genetic diversity to be a function of island size. If hurricanes are major drivers of recolonization and extinction dynamics, we predict widespread signatures of genetic bottlenecks for populations on islands susceptible to hurricanes, but not on larger, less-susceptible islands; we used estimated elevation above sea level and island size as a proxy for hurricane susceptibility (Schoener et al. 2001a(Schoener et al. , 2004. Finally, if natural colonization is prevalent in our system, we predict that patterns of genetic diversity and differentiation of the experimentally colonized reference islands will mimic patterns of similar, but naturally colonized islands.

Study system
Anolis sagrei occurs naturally on islands throughout the Bahamas and more widely throughout the Caribbean and beyond (Losos 2009). It is a trunk-ground habitat specialist, occupying a structural habitat including the ground as well as tree trunks and branches close to the ground. Female lizards become reproductively mature within 6-9 months and lay several single-egg clutches during the breeding season.
Females are capable of storing sperm from multiple males for up to two month (Calsbeek et al. 2007), which might facilitate natural recolonization when mates are absent. Previous work on populations in the Bahamas has shown that viable populations can be founded with as few as two individuals (Kolbe et al. 2012).

Study site and sampling
In 2011, we sampled populations of A. sagrei on 16 small islands ranging in area from ~500 to 3320 m 2 and Staniel Cay, which is considerably larger at ~1.4 km 2 (Fig 1).
The entire study area spans ~14 km (SE to NW) and is ~4.5 km wide (SW to NE) and the two closest neighboring islands (311 and 305) are ~150 m apart from each other ( Fig 1A). We sampled a total of 469 adult individuals (12 -35 per population; mean = 27.6 ± 6.0), including both males and females. Lizards were captured using a noose or by hand. Tail tips were collected and preserved in ethanol for DNA extraction.
Population sizes were estimated using log-linear capture-recapture methods based on a three-day census on each island (Heckel & Roughgarden 1979).

DNA extraction and microsatellite genotyping
Genomic DNA was extracted from tail tips using a standard isopropanol precipitation protocol. We amplified 10 microsatellite loci using PCR and fluorescently labeled primers. We combined primers from the literature ( Bardeleben et al. 2004;Wordley et al. 2011) and newly designed primers ( Genetic variables were allelic richness, expected heterozygosity, and observed heterozygosity.

Population structure
We used two approaches to determine population structure among islands. First, we performed a Bayesian cluster analysis with STRUCTURE v2.3.4 (Rosenberg 2004), using an admixture model and correlated allele frequencies. We assumed gene flow among island populations and sequentially increased the number of clusters (K = 3 -18). Iterations for each run were set to 10, with 10 6 generations and a burn-in of 500,000 steps. We used STRUCTURE HARVESTER (Earl & Vonholdt 2012) to determine the most likely number of clusters based on deltaK. To visualize the result, we used CLUMPP (Jakobsson & Rosenberg 2007) and DISTRUCT (Rosenberg 2004). Second, we used discriminant analysis of principal components (DAPC), a multivariate-distance based approach, to characterize population structure without the assumption of HWE using the R package ADEGENET (Jombart 2008). We retained the first 80 principal-component axes for the discriminant analysis and estimated genetic clusters rather than using population priors.

Contemporary gene flow
To test for contemporary gene flow between islands, we used a Bayesian method implemented in the software BayesAss v3 (Wilson & Rannala 2003). This method estimates the proportion of immigrants for each population from the past two generations and the present generation, without assuming HWE (Wilson & Rannala 2003). The program assumes well differentiated populations (FST ≥ 0.1) and a low migration rate (m = 0.01) among populations (Faubet et al. 2007;Meirmans 2014).
We conducted 10 independent runs starting from random seeds, 20*10 6 iterations, a burn-in of 10 6 generations sampling every 2000 th generation. We optimized the mixing parameters deltaA, deltaM and deltaF according to the user manual. To assure convergence and to select the best run, we examined the tracer file for high effective

Genetic bottlenecks and demographic history
We used two complementary approaches to test for genetic bottlenecks and alternative demographic histories. First, we used Approximate Bayesian Computation to test for the presence of past genetic bottlenecks using the program DIYABC v2.0 (Cornuet et al. 2014). We constructed three historic demographic scenarios (Fig. 2): colonization (scenario 1), population bottleneck (scenario 2), and stable population size (scenario 3). We used uniform prior distributions for all variable parameters (Table S2). The census population size of each population was used as an upper limit for the effective populations size (NP). Time of the demographic event was allowed to vary between 10 and 100 generations and extended in 100-generation intervals up to 1000 generations (Table S2). Because Anolis sagrei becomes sexually mature in 6-9 months (Losos 2009), we assumed one generation per year to be a good estimate.
Further, because lizard populations can grow rapidly after colonization (Schoener & Schoener 1983), the duration of the bottleneck (db) was allowed to vary between 1 and 5 generations. The colonization scenario (Scenario 1) models a demographic scenario, in which a small number of individuals (NC) from an unknown source population (i.e., a ghost population) colonizes the island population. After 1-5 generations, the population reaches the current effective populations size (NP). All study islands are relatively close (≤ 2 km) to considerably larger islands in the east (Fig 1), which represent potential sources of founder individuals and were not sampled. We thus allowed the effective population size of the ghost population to vary between 1,000 and 10,000, which prior work suggests is a realistic range for large islands of these sizes. In the bottleneck scenario (Scenario 2), the island population prior to the bottleneck (NB) is of similar size as the present population (NP). During the bottleneck, NB is reduced to half its size (NH) for one to five generations. Scenario 3 models a population of size NP that has remained constant or larger (NU) for the past 100 generations. Some islands can potentially support larger lizard populations and we thus set the upper limit of the past effective population size to twice the size of the

Genetic diversity and population differentiation
The small island populations had low measures of genetic diversity as compared to Staniel Cay (Table 1). Average allelic richness was lower in all 16 small-island populations (mean AR = 4.2 ± 1.14; range 1.9 -5.9) than in the largest island population, Staniel Cay (AR = 7.5). Observed heterozygosity was lower than expected heterozygosity in all populations, suggesting deviations from HWE that might be the result of inbreeding or small population size. Eleven island populations had relatively high inbreeding coefficients (FIS > 0.25), five islands had moderate inbreeding coefficients (0.25 > FIS > 0.13), and only island 931 had a low inbreeding coefficient of 0.08. All but two populations (i.e., 305 and WBC) had five-or-more loci deviating significantly from HWE, suggesting widespread non-equilibrium.
We did not find evidence for isolation by distance: genetic differentiation between island populations was not correlated with geographic distance between islands (Mantel test: r = 0.03; P = 0.37; dbRDA: DF = 10; F = 1.46; P = 0.27). In addition, we detected no relationship between island characteristics and measures of genetic diversity: island area was not related to genetic diversity (Pearson's correlation: AR r = -0.04; P = 0.90; He r = 0.24; P = 0.46; Ho r = -0.01; P = 0.97), and we found no support for a relationship between island elevation and genetic diversity  (Table S4).
Contemporary gene flow among islands was low overall, with an average migration rate among the 17 islands of m = 0.01 ± 0.0001 (combining the three best runs of BayesAss; Table S5). Migration rates were slightly but significantly higher when restricting the analyses to smaller groups of islands based on geography (North: m=0.03 ± 0.03; Central: m = 0.04 ± 0.04; South: m=0.03 ± 0.03; Table S6). Overall, migration rates were lower than 0.1 and these populations can accordingly be considered independently evolving units (Wilson & Rannala 2003).

Demographic history and population bottlenecks
Evidence for genetic bottlenecks was prevalent across island populations. For the ABC analysis, we distinguished between two types of genetic bottlenecks: populationsize reduction and founder event. Twelve island populations had high posterior probabilities for the founder event scenario, including the four experimentally colonized populations (Fig. 3). Four other island populations (931, 311, 930 and 936) showed support for the founder-event scenario and a scenario with stable population size (absence of genetic bottleneck), albeit with lower posterior probabilities. The average number of founder individuals was 6.21 ± 0.87 (Table 2), based on posterior parameter estimation. For Staniel Cay, the bottleneck scenario was supported with relatively high posterior probabilities, followed by the scenario with stable population size ( Fig. 3).
Results of the ABC analysis were largely consistent with genetic bottlenecks identified with the M-ratio. Fifteen populations had an M-ratio < 0.68, suggesting genetic bottlenecks (Table 1). The only discrepancy was island 332, which showed evidence for a genetic bottleneck in the ABC analysis, but not in the M-ratio test.
Island 931 showed evidence for a genetic bottleneck in the ABC analysis, but not in the M-ratio test. was excluded from the regression, due to the exceptionally large confidence interval.
We found no correlation between island elevation and timing of colonization events (Pearson's correlation: r = 0.09, P = 0.78).

Discussion
Multiple lines of evidence support the hypothesis that hurricanes shape genetic diversity and population structure by driving extinction-recolonization dynamics of Anolis lizard populations on small islands in the Bahamas. We found that smaller, lower-lying islands (which should be most vulnerable to hurricane disturbance; Schoener et al. 2001aSchoener et al. , 2004 have undergone population size reductions as a consequence of founder events during recolonization. Timing of the founder events suggests that populations on larger islands have persisted for longer timespans than populations on smaller islands. Low genetic diversity, population differentiation not related to geographic distance, and high inbreeding coefficients are characteristic signatures of founder effects and provide support for the prediction that strong genetic drift is a major driver of non-adaptive evolution in this system. Principal factors determining presence or absence and timing of genetic bottlenecks are likely related to island characteristics. Population persistence and recolonization depends, among other factors, on island size and elevation relative to storm surges (Schoener et al. 2001a(Schoener et al. , 2004. During natural restoration in times between storms, larger islands are likelier targets for immigration and are less susceptible to extinction when hurricanes reoccur (Losos & Ricklefs 2009;Schoener et al. 2001a). Islands smaller than 1860 m 2 show evidence for complete extinction and recolonization, whereas islands larger than 2000 m 2 (311, 930, 936, Staniel Cay) support multiple demographic scenarios including relative stability as well as colonization ( Fig. 2, Table 2). The only exceptions to this pattern are islands 931 (1070 m 2 ) and 926 (3320 m 2 ). Island 931 has the highest elevation (9 m), the lowest inbreeding coefficient, and the oldest colonization time, and its population has likely persisted during past hurricanes when more low-lying islands were submerged by storm surges. The relatively low elevation of island 926 (1.5 m, among the three lowest values among all islands), might explain the presence of a genetic bottleneck despite its large size. Alternative causes of population size reduction, such as predation, disease and resource depletion, are unlikely explanations for the widespread occurrence of genetic bottlenecks in this system. Once colonized, lizard populations can grow rapidly and persist for at least several decades, and seem to be relatively unaffected by avian predators, resource depletion and disease (Schoener & Schoener 1983). Rodent predators occur occasionally on offshore islands in the Bahamas, but only one study was able to link rat predation to population size reduction (Gasc et al. 2010). Reduced population size prior to hurricanes, however, does not explain population survival or lack of it alone (Schoener et al. 2001b). Predation may have a synergistic effect on survival during hurricanes by making populations more susceptible to extinction and by reducing the likelihood of subsequent population recovery (Schoener et al. 2001b).However, curly tail lizards (Leiocephalus carinatus) were not present on any of the 16 small islands at the time of sampling.
The average number of colonizers from the ABC analysis (6.21 ± 0.87) is strikingly similar to the number of lizards that were used to found the experimentally colonized reference islands (5-10 lizards) in 1977 (Schoener & Schoener 1983).
Lizards are capable of overwater dispersal and females can store sperm from multiple males (Calsbeek et al. 2007), which might facilitate establishment after colonization even in the absence of mates (Schoener & Schoener 1983). Anolis sagrei can survive up to 24 h floating in sea water, but evidence for the actual occurrence and frequency of overwater dispersal is scarce (Logan et al. 2016;Schoener & Schoener 1984).
Migration rates in our study were lower (m = 0.01 ± 0.0001) than previously reported rates of A. sagrei on a similar archipelago in the Bahamas. Logan et al. (2016) attribute gene flow to thermal similarity among islands and proposes that physiological preadaptation is a more important factor than geographic distance. Thus, factors unrelated to geographic distance might promote or constrain rates of gene flow

Conclusion
Our study highlights the evolutionary significance of disturbance events as driver of non-adaptive evolution (Grant et al. 2017). We show that disturbance events can catalyze evolution by driving extinction and recolonization dynamics. In the absence of gene flow, random genetic differentiation persists as a consequence of small population size during recolonization. The finding that timing of the colonization events is related to island size is consistent with the prediction that larger islands are less susceptible to disturbance-induced extinction and better targets for recolonization provide further support that hurricanes are major drivers of extinction: the highest island at 9 m above sea level shows no genetic evidence of recolonization and likely has not experienced complete extinction since low sea levels connected the island populations ~1000 years ago. The second largest island 926 with an elevation of only 1.5 m shows evidence for recent colonization despite its large size. As hurricanes become more frequent, the time intervals between extinction events likely become shorter, leaving less time for recolonization, population recovery and adaptive evolution. Thus, random, non-adaptive processes will increasingly play an important role in the evolution of island populations and coastal ecosystems that are susceptible to hurricane disturbance.

Acknowledgments
We thank Dave Spiller, Todd Palmer, Liam Revell, and Travis Ingram for help during data collection in the field. We also thank Tim Thurman, Rowan Barrett and Clara Cooper-Mullin for valuable comments on the manuscript. KPM was supported by an Erasmus Mundus Scholarship in Evolutionary Biology (MEME).   Table S1. Primers and amplification conditions for microsatellite loci.

Abstract
In polygynous lizards, male-male competition is an important driver of morphological traits associated with dominance during intraspecific encounters. The extent to which females engage in aggressive behavior and thus contribute to phenotypic evolution is not well studied. Most previous studies rely on population density as a proxy for intraspecific competition and lack measures that may better reflect aggressive behavior. We used injury frequency of male and female brown anoles (Anolis sagrei) from 16 island populations to test whether aggressive encounters increase with population density. We further asked whether intraspecific competition is a potential driver of phenotypic traits related to dominance. We found that populations with higher densities had higher proportions of injuries, consistent with increased intraspecific competition. Relative head size (both length and width) of male and female lizards increased with population density, suggesting that larger heads might be advantageous when competition is higher. Morphological changes and injuries were detected in both males and females and show that females play a more active role in intraspecific competition than previously appreciated. More studies are needed to determine whether female aggressive encounters are restricted to intrasexual competition and how morphological traits of females are related to dominance and reproductive success.

Introduction
In many vertebrate species, males compete for mating opportunities and engage regularly in aggressive male-male interactions (Knell 2009;Kokko & Rankin 2006).
In polygynous species, dominant males mate with the majority of females in their territory, whereas subordinate ones reproduce substantially less, if at all (Kokko & Rankin 2006). Thus, the ability to obtain and defend territories disproportionally increases reproductive success for males. One key prediction is that sexual selection shapes phenotypic traits associated with dominance and performance during aggressive male-male interactions.
In lizards, morphological traits associated with dominance and performance in male-male competition are well characterized. Larger males hold larger territories in the wild (Lappin & Husak 2005) and tend to "win" over smaller ones in staged encounters (Jenssen et al. 2005;Lailvaux et al. 2004). When equally sized males are matched, head size and bite force are key predictors of dominance (Gvozdík & Damme 2003;Lailvaux et al. 2004;Perry et al. 2004). In the majority of species, body size is correlated with head size and bite force and thus is the primary trait linked to reproductive success of males (Herrel et al. 2001;Herrel et al. 2007;Herrel et al. 2010;Verwaijen et al. 2002;Wittorski et al. 2016; but see Lappin & Husak 2005).
While the idea that sexual selection shapes body size and head size has been supported empirically in several lizard species (Donihue et al. 2015;Jenssen et al. 2000;Lailvaux et al. 2004;Perry et al. 2004), most studies rely on population density as a proxy for intraspecific competition and lack more direct measures for aggressive behavior. Additionally, the extent to which females engage in aggressive behavior and contribute to phenotypic evolution is not well studied because past studies mostly focus on males.
Females are often described as less aggressive (Claessen et al. 2000) and less territorial (Nunez et al. 1997;Schoener & Schoener 1982), devoting more time to resource acquisition than competition with conspecifics. In several species, females have smaller body and head size and disproportionally weaker bite force than males, even when corrected for body size (Lappin et al. 2006;Wittorski et al. 2016). The observed sexual dimorphism led to the interpretation that intraspecific competition plays a minor role in shaping morphological characters of female lizards (Herrel et al. 2010;Lappin et al. 2006;Wittorski et al. 2016). However, the few available studies that include females suggest that aggressive behavior might be a relevant driver of morphological characters (Calsbeek & Smith 2007;Comendant et al. 2003).
Injuries are a frequent consequence of intraspecific interactions and can be used as a proxy for aggressive behavior (Donihue et al. 2015;Gvozdik 2000). Biting and jaw locking during agonistic encounters leave bite scars and amputated body appendages, such as claws, digits, limbs and tails Donihue et al. 2015;Vervust et al. 2009). While limb related injuries are associated with intraspecific aggression, tail amputations are more ambiguous Gvozdik 2000;Schoener & Schoener 1980). The frequency of tail amputations can vary with species richness and type of predators as well as intraspecific aggression Cooper et al. 2004;Pafilis et al. 2009). Tails are shed more easily (i.e., requiring lower bite force) in populations with predators as compared to predator-free populations, making interpretations of tail loss more complex Cooper et al. 2004;Pafilis et al. 2009).
In this study, we use brown anoles (Anolis sagrei) from 16 island populations to test whether intraspecific competition is a potential driver of phenotypic traits in males and females. We use the frequency of injuries as a proxy for intraspecific aggressive encounters. If intraspecific competition is higher in populations with higher densities, we expect the proportion of injuries to increase with population density. Larger body and head size enhance dominance and performance during male-male competition (Gvozdík & Damme 2003;Jenssen et al. 2005;Lailvaux et al. 2004;Perry et al. 2004), but little is known about female-female interactions. Therefore, if dominance related traits are advantageous during aggressive encounters, we expect larger lizards with larger heads in higher density populations.

Sampling and morphological measurements
Male and female brown anoles were collected in 2011 on 16 small islands in the Bahamas near Staniel Cay where they were the only lizard species present (N = 419; Nfemale = 177; Nmale = 242, Table 1). We took x-ray and scanner images of live lizards.
We scored injuries as the number of missing claws, digits and limbs from the x-ray and scanner images. We measured body size (snout-vent length, or SVL), head length, and head width from x-ray images using the plug-in ObjectJ for the software package Image J (Abràmoff et al. 2004). We used Mosimann's geometric mean (Mosimann 1970) to correct for the effect of body size following Butler & Losos (2002). We calculated the geometric mean of all measured traits as an overall measure of size (SIZE) and used the log transformed ratio of each trait to SIZE (log(trait/SIZE)).
Population density was calculated using the estimated number of individuals per vegetated island area. We used vegetated area, rather than island size, to better represent the actual usable habitat of the lizards. We estimated population size with a log-linear capture-recapture method (Heckel & Roughgarden 1979) based on a threeday census on each island.

Injury and intraspecific competition
As a proxy for intraspecific aggressive encounters, which we predict will increase with population density, we calculated the proportion of injured lizards on each island.
We used a linear model to test whether the percentage of injured lizards is related to population density. We then compared injury percentages between males and females using Welch's t-test to test for sex differences in aggressive encounters. Percentages were arcsine transformed prior to the analysis. All statistical analyses were conducted in R v3.3.2 (Team 2016).

Morphological differentiation and intraspecific competition
To test whether variation in body size and head shape is related to population density, we used separate linear mixed models for each trait. We included sex and population density as fixed effects and island as a random effect (Pinheiro et al. 2015).
Since the morphology of males and females might be affected differently by population density, we included sex*density as interaction term. For instance, if selection on phenotypic traits is stronger for males, we would expect a steeper slope for the relationship between a given trait and population density. If morphological traits of the sexes are similarly affected, we expect no differences in slope between males and females, resulting in a non-significant interaction term. Non-significant interaction terms were excluded from the analysis.

Results
We tested whether the percentage of lizards injured, a proxy for aggressive encounters, increases with population density. Injury percentage was higher in populations with higher densities (P = 0.05; t = 2.15; DF = 14; R 2 = 0.25; Fig 1), suggesting an increase in intraspecific aggression. Injury percentage for females was more variable among populations and on average higher than for males (Females = 0.29% ± 0.20%; Males = 0.24% ± 0.12%), however, the difference was not significant (P = 0.62; t = 0.49; DF = 26.12; Fig 2).
We tested whether body size and head shape varies with population density in the direction predicted by natural selection. Overall, males were larger than females (P < 0.001; t = 41.03; DF = 402; Table 2, Fig 3). Body size did not vary with population density (P = 0.31; t = -1.04; DF = 14; Table 2) and the interaction term (sex*density) was not significant. The random effect of island accounted for 13.8 % of the total variance in body size. Females had relatively longer heads than males, but relative head width did not differ between the sexes (Fig 3). Relative head length and head width increased with population density (head length: P = 0.003; t = 3.68; DF = 14; head width: P = 0.05; t = 2.19; DF = 14; Table 2; Fig 4), but the interaction terms were not significant. The random effect of island explained 12.2% of the total variance for relative head length and 20.9% of variance for relative head width.

Discussion
How morphological traits of males and females relate to intraspecific competition is critical for advancing our understanding of non-random evolutionary processes in polygynous mating systems (Kokko & Rankin 2006). Using morphological data from 16 island populations, we found evidence that intraspecific competition plays a role in shaping head proportions of male and female A. sagrei. Lizards in high-density populations had a higher proportion of injuries as compared to low-density populations, suggesting the former experience higher levels of intraspecific aggression (Donihue et al. 2015). Lizards in these high-density populations had larger heads (longer and wider, corrected for body size), a trait associated with dominance during intraspecific interactions (Gvozdík & Damme 2003;Lailvaux et al. 2004;Perry et al. 2004). Relative head size increased for both sexes and injuries were detected in males and females, suggesting that females play an active role in intraspecific interactions (Calsbeek 2009;Calsbeek & Smith 2007;Donihue et al. 2015). Thus, dominance related traits are likely advantageous for both males and females when intraspecific competition is high.
The relationship between head proportions and dominance during agonistic encounters in males is well documented (Gvozdík & Damme 2003;Lailvaux et al. 2004;Perry et al. 2004). In several lizard species, larger head proportions increase bite force, which is associated with larger territory size and reproductive success (Donihue et al. 2015;Herrel et al. 2010). Our data are in agreement with this prediction and show that both sexes have larger heads (longer and wider heads relative to their body size) in populations with higher densities. In males, natural selection might drive head shape variation due to the beneficial effect of bite force, similar to patterns observed in other lizard species (Donihue et al. 2015;Herrel et al. 2010). Whether head size in females increases social dominance and reproductive success remains to be determined in future studies.
In agreement with previous studies, our data show that the number of aggressive encounters leading to injuries increases with population density (Donihue et al. 2015;Vervust et al. 2009). However, since populations with higher densities have larger heads, and head size typically increases bite force (Herrel et al. 2001;Herrel et al. 2007;Herrel et al. 2010;Verwaijen et al. 2002;Wittorski et al. 2016), we acknowledge that injuries might be more frequent in high-density populations because of greater bite force rather than increased aggressive encounters. Future studies could measure bite force in these populations to evaluate this possibility.
The overall proportion of injuries in our study did not differ between males and females, suggesting that both sexes engage in aggressive behavior. Our findings are inconsistent with the assumption that males engage substantially more in aggressive behaviors than females, in this polygynous mating systems (Kokko & Rankin 2006).
Moreover, empirical studies fail to establish any clear pattern related to sex and injury, with some studies showing no difference between the sexes (Donihue et al. 2015;Vervust et al. 2009), others showing a higher proportion of injuries in males (Gvozdik 2000) or a higher proportion in females (Passos et al. 2013). Whether, and if so how, the sexes differ in aggressive behavior needs to be determined in future studies by including behavioral observations of natural populations or behavioral experiments including females.

Conclusion
Our findings highlight that the current understanding of polygynous mating systems in Anolis lizards is incomplete with respect to females (Kamath & Losos 2017) and their role in intraspecific competition. Sexual dimorphism and head shape variation have been attributed to sexual selection, enhancing male performance during agonistic encounters (Donihue et al. 2015;Jenssen et al. 2005;Lailvaux et al. 2004;Perry et al. 2004). This conclusion is mainly based on the assumption that females engage substantially less in aggressive and territorial behavior. Our data show that relative head size, a trait related to dominance, increased in both males and females with increasing population density. Injuries were detected in both sexes, suggesting that females may engage in aggressive intraspecific encounters as well. Previous associations of head size and bite force indicate that both males and females might benefit from enhanced performance when levels of competition are high (Calsbeek 2009;Calsbeek & Smith 2007;Donihue et al. 2015;Herrel et al. 2010). Our data suggest that females play a more active role in intraspecific competition than previously assumed. More studies are needed to examine whether female aggressive encounters are restricted to intrasexual interactions or involve both sexes, and whether female dominance increases reproductive success, similar to the pattern observed in males.      and females (gray circles). Head proportions are relative head length (left) and relative head width (right).

Abstract
In allopatric species, reproductive isolation evolves through the accumulation of genetic incompatibilities. The degree of divergence required for complete reproductive isolation is highly variable across taxa, which makes the consequences of secondary contact between allopatric species unpredictable. Since the Pleistocene, two species of Anolis lizards, A. carolinensis and A. porcatus, have been allopatric, yet this period of independent evolution has not led to substantial species-specific morphological differentiation and therefore they might not be reproductively isolated. In this study,

Introduction
In allopatric species, reproductive isolation evolves through the accumulation of genetic incompatibilities in geographically separated lineages (Dobzhansky and Dobzhansky 1937;Orr 1995  leading to admixture among genetically distinct source populations in several cases (Kolbe et al. 2007). Despite widespread intraspecific admixture, hybridization between recognized species is considered rare among anoles (Losos 2009 allisoni in Central Cuba ) and A. pulchellus x A. krugi in Puerto Rico (Jezkova, Leal, and Rodríguez-Robles 2013). Hybridization between A. carolinensis and A. porcatus has been suggested repeatedly, mainly because the two species have no species-specific morphological characters despite considerable divergence time (Kolbe et al. 2007;Camposano 2011;Tollis 2013). Sufficient evidence for reproductive isolation or lack thereof has not been shown.
In this study, we examine whether A. porcatus and A. carolinensis are reproductively isolated species, and characterize the genetic consequences of secondary contact in South Miami. We used one mtDNA marker and 18 nuclear microsatellite loci to test whether hybridization has occurred between the two species.
We distinguished between contemporary and historic gene flow and estimated the timing of the admixture event. Discordance between nuclear and cytoplasmic markers is characteristic of hybridization and commonly used to identify hybrid individuals (Toews and Brelsford 2012). If the two species interbreed in South Miami, we expect a high frequency of individuals with mismatches between nuclear genotypes and mtDNA haplotypes.

Sample collection
We sampled 32 Anolis carolinensis individuals from the J.W. Corbett Wildlife Management in South Florida ~200 km north of Miami, 92 green anole individuals from the putative hybrid population in South Miami and 54 A. porcatus individuals from Western Cuba (Table S1). Genomic DNA was extracted from tail tips and liver tissue using a modified ethanol precipitation protocol.

Molecular methods
We amplified a region of 343-571bp of the mtDNA NADH dehydrogenase We amplified 18 microsatellite markers using PCR with fluorescently labeled primers. We used seven newly designed primers (Table S1)  Hill at Yale University. Markers for all samples were analyzed with the software GeneMapper® v4.1 and visually checked to ensure accurate peak calling.

Phylogenetic analysis and haplotype divergence
To determine the species identity of mtDNA haplotypes for individuals sampled in South Miami and the geographic origin of the introduction, we constructed a maximum likelihood phylogeny including samples from the geographic range of both species ( Fig. 1; Table S3). We included 111 individuals of A. porcatus from East and West Cuba spanning the entire native range, 83 individuals of A. carolinensis sampled throughout Florida and 86 individuals from South Miami (Fig. 1A). Anolis loysiana was used as outgroup taxon. Sequences were aligned and visually inspected for accuracy using the MUSCLE plugin in Geneious v7.1.9 (Kearse et al. 2012). We collapsed individual sequences into distinct haplotypes using DNAcollapser implemented in FaBox v1.41 (Villesen 2007). To retain individuals with short mtDNA sequences, we generated two separate alignments. One alignment consisted of 571bp for 200 individuals, resulting in 156 haplotypes. The second alignment was 343bp long and included all 280 samples resulting in 182 haplotypes. We used RAxML v8.0 (Stamatakis 2006) implemented in the CIPRES Science Gateway v3.3 (Miller, Pfeiffer, and Schwartz 2011) to generate maximum likelihood phylogenies. Bootstrap values were obtained from 1000 iterations using rapid bootstrapping.
We used pairwise sequence divergence to determine the degree of nucleotide divergence between native and introduced A. porcatus haplotypes. We identified the genetically most similar individuals between introduced-range and Cuban haplotypes based on the fewest number of pairwise nucleotide differences. Pairwise sequence divergence was calculated as the number of nucleotide differences divided by the sequence length.

Population structure and differentiation
First, we performed a Bayesian cluster analysis with STRUCTURE v2.3.4 (Rosenberg 2004), using the admixture model and correlated allele frequencies. We allowed for gene flow among populations and modeled six different clustering scenarios, sequentially increasing the number of clusters K (1-6). We conducted 10 independent runs for each scenario with a burn-in of 500 000 and 1000 000 MCMC iterations. We used delta K to determine the most likely number of clusters following the Evanno method (Evanno, Regnaut, and Goudet 2005) implemented in STRUCTURE HARVESTER v0.6.94 (Earl and Vonholdt 2012). We combined the genotype proportions of each cluster (q-matrix) from 10 independent runs with CLUMPP (Jakobsson and Rosenberg 2007) and visualized the results with the R package ggplot2 v2.1.0 (Wickham 2011). We repeated the Bayesian cluster analysis with population pairs (SFL-MIA and WCU-MIA) and WCU separately to identify potential population substructure. Model parameters were used as described above.
Second, we used discriminant analysis of principal components (DAPC) to determine the degree of differentiation between clusters using the R package adegenet (Jombart 2008). To characterize and find genetic clusters, DAPC uses a multivariate approach and PCA transformed allele frequencies. In contrast to the Bayesian clustering approach, DAPC does not rely on specific population model assumptions, such as HWE. The number of clusters (K) was sequentially increased starting with one cluster.
The model fit for K clusters was determined with the Bayesian Information Criterion (BIC).

Maximum likelihood and ABC modeling of historic admixture
To detect historic gene flow, we used a tree-based maximum likelihood approach with the program TreeMix (Pickrell and Pritchard 2012). This approach uses allele frequencies to model relatedness among populations as a non-bifurcating tree.
Migration edges are added as additional branches to a bifurcating tree allowing for population ancestry from more than one parental populations. Migration edges are added stepwise to the tree model until the covariance of the model best matches the covariance of the data. Residual matrices were used to determine the model fit.
Positive residuals indicate greater genetic variation in the population than explained by the simple tree model suggesting admixture (Pickrell and Pritchard 2012). The model assumes migration within a single generation. The fraction of alleles derived from migration is represented as weight of the migration edge.
To infer the timing of the admixture event, we used approximate Bayesian computation to model the demographic history of the three populations using DIY-ABC v2.0 (Cornuet et al. 2014). A set of summary statistics was used to assess the fit between simulated datasets and empirical data. We used mean number of alleles, mean genetic diversity, pairwise FST-values and the maximum likelihood coefficient of admixture (Choisy, Franck, and Cornuet 2004). The demographic scenario simulates divergence between SFL and WCU and a more recent admixture event that gave rise to the MIA population. The prior for the divergence between SFL and WCU was set between [6000 000-13 000 000] generations, based on previous divergence time estimates (Campbell-Staton et al. 2012;Tollis and Boissinot 2014). We set the prior for the effective population size as [100-10 000] using a uniform prior distribution. To estimate timing of the admixture event, we used a prior of [1-4000] generations assuming one generation per year. We simulated 1000 000 datasets and used the 1000 dataset with the smallest Euclidean distance to the empirical data for parameter estimation.

Phylogenetic analysis and mtDNA haplotype divergence
We constructed a maximum likelihood phylogeny from mtDNA haplotypes to determine the haplotype identity of individuals sampled in South Miami. Individuals sampled in South Miami (N=86) were not monophyletic. Thirty samples (35%) representing six haplotypes were nested within a well-supported clade of A. porcatus from West Cuba. Fifty-six samples (65%) clustered with A. carolinensis from South Florida, representing 20 haplotypes (Fig. 1C). Haplotypes from A. carolinensis and A.
porcatus were co-distributed across the study area in South Miami (Fig. 1B).
The 571bp alignment resulted in a total of 155 haplotypes, 41 from South Florida, 43 from Central and north Florida, 52 from West Cuba and 19 from East Cuba. Anolis porcatus from West Cuba is sister to the monophyletic A. carolinensis (Fig S1). Since, the 343bp alignment resulted in an overall similar tree topology (Fig. S2), we focus on the 571bp in the main text. Individual haplotypes and sampling locations for the 343bp alignment can be accessed in the supplementary material.
Introduced A. porcatus haplotypes were nested in a well-supported clade of A.
porcatus from seven sampling locations near Havana in West Cuba ( Figure 1A).
Branches within the clade were not well supported (bootstrap < 95), which limits more specific identification of potential source locations of the introduction. Average sequence divergence between introduced A. porcatus haplotypes and the genetically most similar ones from Cuba ranged from 0.0-1.75% divergence (mean = 1.14 ± 0.68%; Table 1). One individual from South Miami (MIA640) shared the same haplotype with one individual from Havana (JJK2796).

Genetic diversity and differentiation using microsatellite loci
Genetic structure and diversity was assessed for populations from WCU, SFL and MIA using nuclear microsatellite markers. Three microsatellite loci (Ac2, F06, g01) deviated significantly (P < 0.05) from HWE. Excluding those loci from the analysis did not affect the results and we thus included them in subsequent analyses. Allelic richness was similar across populations (mean AR = 10.62 ± 0.55; Table 2). Observed heterozygosity was lower than expected heterozygosity in all populations (mean Ho = 0.70 ± 0.03; mean He = 0.81 ± 0.03). FST-values showed similar degrees of differentiation between populations (mean pairwise FST = 0.08 ± 0.01, Table 2).
Individual allele frequencies of all markers are shown in Fig. S2.

Population structure and differentiation
The genetic cluster analysis recovered three distinct genetic clusters (

Maximum likelihood and ABC modelling of historic admixture
Tree-based maximum likelihood analysis of microsatellite markers supported one migration event between WCU and MIA ( Fig. 2A). The weighted migration edge suggested that ~33% of the nuclear genetic ancestry in the MIA population is derived from WCU. Including the migration edge significantly improved the fit of the model as compared to a strictly bifurcating tree model (P < 0.001, Table S4). The migration model explained 99% of the total variance in the data, whereas the strictly bifurcating the tree model accounted for 80% (Fig. S8).
Time estimates from the ABC analysis suggest that the admixture event between MIA and WCU occurred within the last 245 -2670 generations (median TA = 887, mode TA = 554, 95% CI 245-2670; Table 3, Fig. 3). The median rate of admixture was RA = 0.24 (95% CI 0.14-0.35), which is consistent with the maximum likelihood coefficient of admixture = 0.31 obtained from the summary statistics (Table S5).
Estimates for the remaining parameters used in the model are shown in Table S5.
Summary statistics generated from the posterior probability distribution show similar values compared to the observed data and are largely non-significant, suggesting that modeled parameters provide a good fit for the data (Table S5). Performance measures for parameter estimates were consistently low and indicate accurate estimates (Table   S6).

Discussion
In an effort to characterize the genetic consequences of secondary contact between the native A. carolinensis and the closely related introduced A. porcatus, our data show evidence for past hybridization followed by differentiation of the hybrid population. We found discordance between nuclear microsatellite markers and mtDNA haplotypes in the South Miami population, which is indicative for hybridization (Hailer et al. 2012;Miller et al. 2012;Roy et al. 2015). Thirty-five percent of mitochondrial haplotypes in the South Miami population are derived from A. porcatus from West Cuba and 65% from the native A. carolinensis in south Florida.
Genetic cluster analyses of nuclear markers show that the South Miami population is homogeneous and genetically distinct from populations of both parental species, which is characteristic of hybrid ancestry rather than ongoing hybridization (James and Abbott 2005;Kronforst et al. 2006;Mims et al. 2010;Keller et al. 2014;Lavretsky et al. 2015). Tree-based maximum likelihood analysis confirms that ~33% of the nuclear genetic ancestry is derived from West Cuba. Strikingly, the proportion of nuclear ancestry derived from West Cuba (~33%) was remarkably similar to the proportion of A. porcatus mtDNA haplotypes in South Miami (~35%). Thus, reproductive barriers between A. porcatus and A. carolinensis are weak or absent despite considerable divergence in allopatry. Thus, secondary contact after species introduction has led to a genetically distinct hybrid population.
Time estimates from ABC analyses suggest that hybridization occurred between 245 -2670 generations ago with a skewed distribution towards the present (Fig. 3 (Powell 1992). However, a thorough evaluation of morphological differences between the species concluded that morphological characters are inadequate for species delimitation (Camposano 2011). Genetic evaluation of the A. carolinensis species complex revealed paraphyly of A. porcatus, dividing this species into an eastern and western clades in Cuba, with the western clade being sister to A. carolinensis (Glor et al. 2004;Glor, Losos, and Larson 2005). Our study provides a thorough genetic evaluation of species boundaries between A. carolinensis and A. porcatus. According to the biological species concept, populations of distinct species are incapable of effectively interbreeding with one another (Mayr 1982), which is inconsistent with the findings of our study. Thus, genetic evidence for successful hybridization, as well as morphological similarities between the two species (Camposano 2011) reinforces a proposal to revise the taxonomy of A. carolinensis and A. porcatus from West Cuba, according to guidelines of the biological species concept (Mayr 1982).
Several Anolis species have been introduced to Florida and some from multiple native-range source populations in Cuba (Kolbe et al. 2007). In agreement with previously collected A. porcatus haplotypes from Miami (Kolbe et al. 2007), the phylogenetic analysis identified sampling sites located near Havana in West Cuba as potential source of the introduction. Haplotypes from these locations did not show evidence for further spatial structure (Fig. 1A), which is consistent with a single West Cuban populations based on microsatellite data ( Fig S6). Thus, the source locations likely resemble a single panmictic population. However, a more comprehensive sampling approach of West Cuban populations is needed to clarify whether population structure exists and whether the introduction of A. porcatus involves a single or multiple independent introductions.

Conclusion
The aim of this study was to characterize the genetic consequences of secondary    Table S2: Primer sequences and annealing temperatures (Tm) for the 18 microsatellite loci and partial mtDNA region of the NADH dehydrogenase subunit 2.   Genetic clusters from the Bayesian cluster analysis for K = 2 and K = 3. The most likely number of cluster was K = 3. C) DAPC analysis with K = 3 clusters, axis 1 accounting for 51% variation and axis 2 for 49%. Bottom figures show the density of each cluster for axis 1 (left) and axis 2 (right).