Chromosomal inversions and common quail diversification in Macaronesian archipelagos and the mainland

Recommended citation: Ravagni, S. (2023). Chromosomal inversions and common quail diversification in Macaronesian archipelagos and the mainland. PhD Thesis. Universidad de Sevilla, Seville, Spain. Cover illustrations are original creations of Paula Martin. Doctoral dissertation CHROMOSOMAL INVERSIONS AND COMMON QUAIL DIVERSIFICATION IN MACARONESIAN ARCHIPELAGOS AND THE MAINLAND SARA RAVAGNI Advisors: Carles Vilà Arbones and Inés Sánchez Donoso Programa de Doctorado en Biología Integrada A dissertation submitted in fulfillment of the requirements for the degree of Doctor of Philosophy at University of Seville, 2023 APPROVED BY Dr. Carles Vilà Arbones Research Professor at the Doñana Biological Station (EBD-CSIC) Dra. Inés Sánchez Donoso Researcher at the Doñana Biological Station (EBD-CSIC) This thesis was funded by projects CGL2016-75227-P and PID2019-108163GB-I00 and a PhD fellowship (BES-2017- 081291) from the Spanish Ministry of Economy and Competitiveness. A ARAY A CARLOS ALLA MIA FAMIGLIA “Nothing in biology makes sense except in the light of evolution” Theodosius Dobzhansky, 1973 “If I have seen further, it is by standing on the shoulders of giants” Isaac Newton, 1676 INTRODUCTION INTRODUCTION 16 | PA G E Islands offer a unique opportunity to investigate mechanisms leading to population differentiation and adaptation. The replicated nature of the processes of differentiation, occasionally resulting in convergent evolution, offers the possibility to investigate molecular mechanisms involved in diversification (Figure 1). Common patterns of molecular change appearing independently over different islands are a strong indication of repeated selective forces. The search of patterns of selection over the genome traditionally focuses on protein coding sequences. These searches normally target genes subject to strong selective forces. However, recent research and advances in genome sequencing approaches have paved the way for the study of chromosomal arrangements that can have an important effect on the evolutionary history and adaptation of populations. This thesis focuses on the study of the molecular mechanisms associated with the differentiation of island and mainland lineages of common quails (Coturnix coturnix), with an emphasis on the differentiation due to chromosomal rearrangements. FIGURE 1. Islands are ideal systems to study many evolutionary processes, due to their small size and isolated nature. They are affected by geological changes that convert them into dynamic systems, where the interplay of colonization, evolution and extinction determines the community composition. From Whittaker et al. (2017). CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 17 | PA G E ISLANDS AS STUDY SYSTEM FOR EVOLUTIONARY PROCESSES Islands have been fascinating biologists since the first expeditions around the globe in the 19th century and they inspired Darwin and Wallace in the development of the theory of evolution through natural selection. Later, the seminal work on the theory of island biogeography by Wilson and MacArthur (1963, 1967) opened the path to the development of many theories that were initially based on island biodiversity but that applied to more general questions and contexts in ecology, conservation and biogeography (Santos et al., 2016). As Mayr stated in 1967: “An island may demonstrate certain biological phenomena almost with the clarity of a test-tube experiment; indeed, every island biota is an experiment of its own”. Islands offer unique opportunities to study ecological and evolutionary processes in discrete and isolated areas that differ in geological age, area, elevation and habitat diversity and present simplified biotas. This allows the simultaneous study of the interaction of multiple factors that affect speciation and community composition in a simplified scenario. Islands can be considered as replicates of ecological and evolutionary processes under comparable conditions and therefore they have influenced our knowledge of lineage differentiation, species distribution and community ecology (Losos & Ricklefs, 2009; Warren et al., 2015; Santos et al., 2016; Patiño et al., 2017; Whittaker et al., 2017). Islands are dynamic systems that change over time and present peculiarities due to their unique geological and ecological characteristics. Because of their smaller sizes, islands host a reduced number of species compared to mainland systems (e.g., Gonzalez-Castro et al., 2012; Traveset et al., 2016). The isolation and the hard boundaries limit immigration from the mainland and favor the establishment of species that have long dispersal ability (Wright et al., 2016; Shaw et al., 2016; Heads et al., 2017; but see Heleno & Vargas, 2015; Garcia-Olivares et al., 2017). Therefore, islands present an unbalanced set of taxa compared to the mainland, which is known as disharmony, and many ecological niches are left vacant. Those are occupied by species that broaden their niche width through “ecological release” (Diamond, 1970), leading to peculiar lineages with unusual adaptations (e.g., woodpecker finches and blood-sucking finches on the Galapagos Islands, dragon’s blood trees on Socotra and Macaronesia, pollinating reptiles on the Canary Islands). The ecological conditions can facilitate the diversification of some lineages into a broad array of species with different ecological, phenotypic and behavioral traits, leading to spectacular adaptive radiations such as Darwin’s finches in the Galapagos, drosophilid flies in Hawaii and Anolis lizards and Eleutherodactylus frogs in Caribbean Islands (Losos et al., 1998; INTRODUCTION 18 | PA G E Losos & Schneider, 2009; Lapoint et al., 2013; Lamichhaney et al., 2015; Dugo-Cota et al., 2019; Gillespie et al., 2020). In some cases, the species divergence can occur exceptionally fast, such as in Darwin's finches where a new reproductively isolated lineage evolved in just a few generations (Lamichhaney et al., 2018). In other cases, differentiation can be the result of a very long process. The mechanisms that lead to the differentiation of insular forms can be diverse, including genetic drift, natural selection and phenotypic plasticity, and they are often unknown (Baeckens & Van Damme, 2020). However, the current use of next generation sequencing technologies makes the identification of large numbers of genomic markers available and affordable even for non-model species. These genomic data can shed light on the role played by each factor on the origin of the diversity in islands (e.g. Lamichhaney et al., 2015; Armstrong et al., 2018; Salces-Castellano et al., 2019; Martin et al., 2021). One important question still in debate is related to the relative importance of isolation and gene flow in the speciation process, and molecular data from archipelagos, comparing within- and among-island diversity with mainland areas, can provide further understanding on the interplay between both drivers of speciation (Patiño et al., 2016). Although different evolutionary processes and mechanisms drive differentiation and speciation, there are some morphological, ecological and behavioral traits that appear repeatedly in many taxa and different archipelagos, suggesting some kind of deterministic evolutionary forces acting on islands. New large phenotypic datasets have confirmed that these traits, that are defined as the “island syndrome” (Baeckens & Van Damme, 2020), have undergone a repeated convergent evolution across multiple lineages and therefore represent a frequent phenomenon in islands (Losos & Ricklefs, 2009; Lens et al., 2013; Wright et al., 2016; Whittaker et al., 2017; Benítez- López et al., 2021; but see Lomolino et al., 2013; Itescu et al., 2014). These traits include the tendency towards dwarfism or gigantism in mammals and in plants, loss of dispersal ability in birds and insects, tameness in many vertebrates and reduced clutch size in birds. Rails provide a striking example of parallel evolution of insular flightlessness in birds, with 31 out of 150 species (including extinct ones) independently evolving flightlessness from a volant ancestor (McCall et al., 1998; Kirchman, 2012; Hume et al., 2019). There is increasing evidence that this is not limited to a few taxa and has been observed in four orders of birds, including actively volant species (Wright et al., 2016). The loss of dispersal capability seems to be related to changes in flight muscles and an increased investment from forelimbs to hindlimbs in response to a reduced number of predators (Wright et al., 2016). CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 19 | PA G E Islands are interesting and unique systems, often constituting biodiversity hotspots. They can hold high proportions of endemic species that generally present small effective population sizes and geographically limited distribution ranges, and sometimes they host relictual or paleoendemic lineages. The island syndrome traits, in particular the loss of dispersal ability and the naïveté towards predators, make insular forms more vulnerable and prone to human-driven extinction. Since human pressure is increasing, reducing habitat availability and leading to alien species introduction, we need to broaden our knowledge on island species and communities to increase our efforts in the conservation of island ecosystems (Wood et al., 2017). MACARONESIAN ARCHIPELAGOS The Macaronesian region includes 40 main volcanic islands that form the archipelagos of Azores, Madeira, Selvagens, Canary Islands and Cape Verde, located in the Northeastern part of the Atlantic Ocean, ranging from 39°N 31°W to 15°N 23°W (Figure 2). The Macaronesian archipelagos show a high variation in ecological conditions and habitat diversity due to their intrinsic characteristics, such as isolation from the mainland (from 96 km for the Canary Islands to 1365 km for Azores), latitude (affecting the climate, from oceanic temperate in the Azores, to Mediterranean in the more central archipelagos and warm arid in Cape Verde), geological ages (35,000 years for Alegranza in the Canary Islands to 25.7 million years for Selvagem Grande) and altitude (maximum height from 130 to 3715 m, above sea level; Fernandez- Palacios, 2010; Illera et al., 2012; Florencio et al., 2021). Macaronesia is recognized as one of the world's biodiversity hotspots due to their high proportion of endemisms: around 25% of the terrestrial taxa are endemic to the region (Izquierdo et al., 2004; Arechavaleta et al., 2005; Borges et al., 2008, 2010) and it is home to an increasing number of threatened species. The unique combination of isolation, evolutionary history, and diverse environmental and geographic conditions, as well as high species diversity, have made Macaronesia an ideal field of study for biogeographers, evolutionary biologists, and conservationists, with an increasing number of investigations focused on this region (Florencio et al., 2021). The distance from the continent limits the distribution of reptiles, which are absent from the more distant Azorean archipelago, and non-volant mammals, found only in the closest Canary Islands. For example, Mateo et al. (2022) list 31 species of reptiles for Cape Verde (8 of them probably extinct), 17 for the Canary Islands (2 probably extinct), 2 for Selvagens, 1 for Madeira and 0 for Azores. INTRODUCTION 20 | PA G E FIGURE 2. Map of Macaronesia, showing the different archipelagos and their maximum age. CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 21 | PA G E The Macaronesian avifauna is also characterized by a high number of endemic species and subspecies (Stattersfield et al., 1998; Martín & Lorenzo, 2001). The presence of several bird subspecies differentiated among archipelagos, combined with geographical features such as the varying distance from the mainland, provides the ideal environment to study colonization patterns, genetic diversity, and past and present gene flow between insular forms and their mainland counterparts (e.g., Armstrong et al., 2018; Martin et al., 2021; Recuerda et al., 2021). CHROMOSOMAL INVERSIONS AS MECHANISM OF POLYMORPHISM MAINTENANCE AND DIFFERENTIATION With the reduction in the prices of next generation sequencing, more and more species and individuals can be sequenced and the variation that can be detected is not only limited to a reduced number of microsatellites or single nucleotide polymorphisms (SNPs) and can span entire genomes. Chromosomal rearrangements, including inversions, copy number variants and insertions-deletions, could be the molecular basis of adaptation and speciation in many cases and therefore should be integrated in the study of ecological and evolutionary processes (Wellenreuther et al., 2019; Merot et al., 2020; Weissensteiner et al., 2020). Among them, inversions are an important source of genetic differentiation and are considered as an evolutionary mechanism that can be involved in local adaptation and, eventually, in speciation (Wellenreuther & Bernatchez, 2018). Chromosomal inversions were first described based on genetic maps of different species of Drosophila by Sturtevant in 1921 and they were detected, a few years later, in the polytene chromosomes in D. melanogaster (Dobzhansky & Sturtevant, 1938). It was demonstrated that they have a huge effect on transmission because they cause suppression of recombination when found in a heterozygous state (Kirkpatrick, 2010). Chromosomal inversions take place when a segment of a chromosome is reversed within the same chromosomal location and therefore the same alleles are present, but they are found in a different order compared to the ancestral state (Figure 3). In heterozygous individuals, carrying one chromosome with the inversion and the homologous one without the inversion, the two chromosomes cannot properly pair during meiosis and this makes recombination very difficult in this region. Thus, the sequence within the inversion evolves independently from the homologous sequence in chromosomes that maintain the original gene order (Faria et al., 2019a). The suppression of recombination in the inverted sequence can play an important role in INTRODUCTION 22 | PA G E local adaptation and diversification because it allows maintaining favorable combinations of alleles on the same chromosome despite gene flow. Therefore, the presence of an inversion can lead to the formation of a supergene (Wellenreuther & Bernatchez, 2018; Schwander et al., 2014). Supergenes are defined as sets of tightly linked co-adapted alleles that control complex traits and are consistently inherited together as a single Mendelian unit (Thompson & Jiggins, 2014). The interest in chromosomal inversions and their related supergenes has recently increased and has revealed their widespread presence in multiple taxa. They are linked to the maintenance of complex morphological and behavioral adaptive phenotypes segregating within species and are often associated with the origin and evolution of sex chromosomes. They also favor reproductive isolation and speciation (Schwander et al., 2014; Thompson & Jiggins, 2014; Wellenreuther & Bernatchez, 2018). FIGURE 3. Chromosomal inversions are structural variations in which a segment of a chromosome breaks off, rotates and reattaches to the same chromosome in the reverse orientation. This results in a change in the linear order of genes within the inverted segment and can affect the expression of genes by altering their regulatory regions or disrupting the reading frame. Many inversions affect behavioral and ecological traits. The Asian swallowtail butterfly (Papilio polytes) exhibits female-limited mimicry that is controlled by a single gene within an inversion (Kunte et al., 2014; Nishikawa et al., 2015). In red fire ants (Solenopsis invicta), a 13-Mb genomic region, harboring three large chromosomal inversions, controls colony social organization, determining if a colony has multiple or a single queen (Wang et al., 2013; Helleu et al., 2022). Interestingly, other four CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 23 | PA G E supergenes have been described in other ant species and all of them govern social organization in colonies (Kay et al., 2022). In other species inversions are directly linked to different migratory strategies, such as in the Atlantic cod (Gadus morhua), where four inversions maintain the divergence between stationary and migratory ecotypes (Berg et al., 2017; Matschiner et al., 2022), and in willow warblers (Phylloscopus trochilus), in which two different migratory patterns are determined by two inversions in different chromosomes (Lundberg et al., 2023). In two bird species inversions show dramatic effects on the mating systems (Figure 4): in the ruff Calidris pugnax three extremely different alternative male morphs, differing in mating behavior, size, plumage and other characteristics, are defined by the ancestral order and two different inversion haplotypes (Lamichhaney et al., 2016). The territorial and parental behavior is also affected by inversions in white-throated sparrows (Zonotrichia albicollis), but in this case both males and females show differences and mate exclusively with the same morph, leading to an unusual situation that results in four effective ‘sexes’ (Thomas et al., 2008; Zinzow-Kramer et al., 2015; Campagna, 2016). Inversions can span large portions of the genome, being 8.4 megabases (Mbp) the average size, and encompass hundreds of genes (Wellenreuther & Bernatchez, 2018). The largest inversion -before this thesiswas found in the white-throated sparrow and it involves 100 Mbp, including a large portion of the genome, around 10% (Thomas et al., 2008; Campagna, 2016). Several species present multiple inversions, for example in the common fruit fly Drosophila melanogaster at least eight large inversions have been described (Corbett-Detig & Hartl, 2012; Kapun et al., 2016) and in the marine snail Littorina saxatilis 17 putative inversions were found in different frequencies according to the habitat (Faria et al., 2019b; Koch et al., 2022). In some cases, inversions are found in adjacent positions and affect the same traits, such as the case of three collinear rearrangements that govern the social organization of colonies in the red fire ant Solenopsis invicta (Wang et al., 2013; Helleu et al., 2022) and the mimetic coloration pattern in Heliconius numata (Jay et al. 2018, 2021). Selection could favor the persistence of larger inversions because they maintain linkage among a larger number of distant loci. It has been shown that a limited level of recombination can occur between inverted and non-inverted regions due to double crossover and gene conversion and it can increase the genetic variation within the inversions (Kirkpatrick and Burton, 2006; Kirkpatrick, 2010; Wellenreuther & Bernatchez, 2018; Korunes & Noor, 2019). INTRODUCTION 24 | PA G E FIGURE 4. Two examples showing the dramatic impact that inversions can have on both phenotype and mating behavior. Above, the white-striped and tan-striped forms of the whitethroated sparrow that differ in territorial and parental behavior as a result of a chromosomal inversion (photo credit: Simon Pierre Barrette). Below, a female and two alternative male morphs in ruffs that differ in the presence of an inversion (photo credit: Winfried Wisniewski, from Taylor & Campagna, 2016). CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 31 | PA G E So far, very little is known about the genetic structure and the genomic basis of the morphological polymorphism that characterizes this bird species. Some preliminary analyses carried out with microsatellites (I. Sanchez-Donoso and C. Vilà, personal communication) did not show any genetic differentiation among quail populations. However, the coexistence of two morphotypes of quails with differences in pigmentation and perhaps migratory behavior in the same localities seems difficult to explain in genetically homogenous populations and should be studied in more detail. A proper understanding of the sources of phenotypic variation in common quails requires the genome-wide characterization of quails from diverse origins, with differentiated phenotypes, and an investigation of the degree of admixture between the different genetic lineages. OBJECTIVES OBJECTIVES 34 | PA G E The general objective of this thesis is to study the genomic basis of the mechanisms that facilitate the maintenance of polymorphisms within a species and its diversification in insular environments, focusing mainly on the role of chromosomal rearrangements in common quails. By studying the structural variability of the common quail genome, we aim to elucidate the processes that enable the persistence of polymorphisms in small island populations and facilitate sympatric divergence. To reach the main goal of the thesis, I focused on different specific questions: i) Determination of the biases in the use of a well-known software, STRUCTURE, to identify hybrids and introgressed individuals (Chapter 1). This knowledge is central to properly understanding the estimates of admixture between populations. ii) Identification of the genomic basis of morphological polymorphism in populations of common quails and assessment of the role of a chromosomal inversion (Chapter 2). iii) Reconstruction of the evolutionary history of an island endemic: divergence and genetic differentiation in the Azorean common quail and the role of a chromosomal inversion (Chapter 3). iv) Evaluation of the role of chromosomal inversions on the genome evolution and diversification of common quails (Chapter 4). MATERIAL AND METHODS MATERIAL AND METHODS 36 | PA G E FORWARD SIMULATIONS AND POPULATION GENETICS ANALYSES The program STRUCTURE (Pritchard et al., 2000) is often used to assess the degree of admixture between two species or separated populations. However, previous analyses have shown that the selection of samples can affect the inferences obtained by the program. Thus, we set out to investigate the reliability of the admixture (level of introgression) estimates generated by this program and that often are the basis for conservation and management programs. We simulated genetic data for two populations and carried out forward simulations to test different scenarios with varying levels of admixture. First, we generated the initial populations with different genetic differentiation and varying numbers of two types of genetic markers. In the second step we used a custom Python script to simulate asymmetric gene flow from one population to the other under varying conditions. This resulted in more than 5,000 different datasets that were analyzed with STRUCTURE and later compared with the real ancestry values calculated for each individual during the simulations (Chapter 1). The use of simulations for testing software reliability offers numerous advantages, especially when combined with empirical data. Simulations enable the generation of a large number of conditions and replicates, which would otherwise be impossible to test. However, it is crucial to critically assess the assumptions made during simulations, as they could affect the conclusions and the application to natural populations. Thus, simulations allow generating datasets with a structure similar to empirical datasets, while controlling diverse parameters. This makes it possible to assess how diverse methodological and biological factors could impact the accuracy of ancestry estimates. FIELD SAMPLING Male common quails were captured during the breeding season of several years (from 2007 to 2019) in different localities along their distribution range: Spain, Portugal, Morocco, Italy, Netherlands and the Macaronesian archipelagos of Azores, Madeira, Canary Islands and Cape Verde. The breeding season occurs from March to June in most localities except Cape Verde, where breeding takes place in November. Most quails were captured using a horizontal net. We first determined the presence of males in specific locations using recorded female calls and we extended the net CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 37 | PA G E over the cereal crop where one male quail was heard responding to the female calls (Figure 7). We tried to attract it using the female call recordings and once the individual approached the decoy and was under the net, we scared it, so the bird tried to fly and was then trapped in the net. Although common quails are difficult to capture, this method allows sampling male individuals during the breeding season. However, it is exceedingly difficult to capture females and most studies are limited to the data provided by one sex. This can be a strong limitation in some analyses and should be taken into account. FIGURE 7. Once a male is detected by its response to female calls from an electronic decoy, an horizontal net is extended on the cereal field to capture it. A small number of quails, 5 from Italy and 7 from southern Spain, were captured during migration using vertical mist nets and playing recordings of both female and male calls. For the Azores, the sampling was complemented with birds provided by hunters. From each captured quail we collected a high number of phenotypic data, such as weight, tarsus length, pectoral lipid band and length of different feathers. We also took pictures to study throat pigmentation and design patterns (Figure 8). The age was inferred from the pattern of molt and wear of the primaries (Saint Jalme & Guyomarc’h, 1995). MATERIAL AND METHODS 38 | PA G E FIGURE 8. Collection of phenotypic (morphology and pigmentation data) from a captured common quail male. CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 39 | PA G E We took two primary feathers from captured males: one grown during the previous breeding season and the other one grown during the wintering season. The isotopic signature of these feathers corresponds to that of the environment in the region where the feathers were formed. This information can provide insights into the approximate areas where the quails spent the previous Summer and Winter, allowing us to infer their migratory behavior (Rubenstein & Hobson, 2004). For genomic analyses we extracted 100 µl of blood from the jugular vein. The blood samples were stored in a buffer for the preservation of DNA and were frozen immediately upon arrival to the lab. Birds were released in the same sampling locality as soon as all data and samples were collected. We also sampled the testicles of 15 individuals from Portugal and Catalonia (NE Spain) that were used for immunofluorescence essays (Chapter 2). In this case, we dissected the birds in the field and testicles were immediately frozen in a portable tank with liquid nitrogen to ensure immediate fast freezing without any cell degradation. PHENOTYPIC ANALYSES We analyzed different morphological parameters and we statistically compared oneyear old male common quails of different genotypes in Chapter 2 (homokaryotypes with and without the inversions and heterokaryotypes). Among all the measures taken, we chose to represent phenotypic variability with weight, width of the lipid band, tarsus length and throat pigmentation, which was classified into six levels, and we calculated a modified Holynski Index (Holynski, 1965; Lockwood et al., 1998; Fiedler, 2005) based on the lengths of the primary feathers, the second secondary feather and the folded wing. This index provides an estimate of the wing pointedness, which is directly related to the flight efficiency and therefore to the migratory capability (Figure 9). In Chapter 3, the objective was to test if common quails in the Azores were morphologically different from quails from the continent and other Macaronesian islands. As we did not know the genotype of all individuals and we demonstrated in the previous chapter that different genotypes coexist in some localities, we opted for another approach, a Gaussian Mixture Model (GMM) cluster analysis (Fraley & Raftery, 1998; Fraley & Raftery, 2002), that does not assume prior differences based on the geographic location or the genotype. GMM is a probabilistic model that MATERIAL AND METHODS 40 | PA G E assumes that data are generated under different Gaussian distributions and therefore it estimates various parameters that are used by a clustering algorithm to classify the samples into clusters. FIGURE 9. The wing pointedness, calculated through the Holynski Index, is related to the flight efficiency: the more rounded, the less efficient. IMMUNOFLUORESCENCE ASSAYS We used fluorescence in situ hybridization (FISH) as a complementary method to confirm the presence of chromosomal inversions (Chapter 2). This technique relies on the use of fluorescently labelled antibodies that selectively bind to target antigens, enabling their visualization under a fluorescence microscope. We identified the position, on the Japanese quail genome, of probes developed for chicken (Zlotina et al., 2012) and we selected two probes, one chosen to hybridise inside the presumed inverted region in chromosome 1 (K19) and the other outside this region (J23). The two probes were marked with different fluorochromes, and their relative position was visualized using a fluorescence microscope, allowing the measurement of the distance along the chromosome of the two fluorochromes. If a large region of the chromosome was inverted, the two fluorescent marks would appear at different distances in individuals homozygotes with and without the inversion (homokaryotypes, with two copies of the same chromosomal types), as illustrated by individuals AA and BB in Figure 10. In individuals carrying one chromosome with the inversion and the other one without the inversion (heterokaryotypes), the CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 47 | PA G E Azorean quail population is genetically and morphologically differentiated from other populations. Azorean quails are characterized by a smaller size, dark throat pigmentation and the absence of migratory behavior. Contrary to the assumptions that quails were introduced by Portuguese settlers, our findings point towards a long evolution in isolation, with a divergence from other existing populations dating back to at least 0.8 million years ago. We discovered that some individuals in this population carry the inversion in chromosome 1 that we previously identified in quails from the mainland and other Macaronesian archipelagos. Despite only half of the Azorean quails carrying this inversion, none exhibit migratory behavior. This implies that the inversion could play a different role in the Azorean archipelago and suggests that the loss of migration could be associated with different genomic mechanisms. Furthermore, the two haplotypes of the inversion are very different from those in other localities and have coexisted in the Azores for a long time, supporting the role of balancing selection in maintaining this polymorphism. As Azorean quails constitute a separate lineage, conservation practices should take into account their unique genetic composition and avoid translocations from other archipelagos or the continent. Finally, the last chapter focuses on the genomic effects of chromosomal inversions in the evolution of common quail populations. In Chapter 4 the de novo assembly of a common quail genome allowed a better characterization of the inversion in chromosome 1 and the detection of two additional inversions in chromosome 2, which span 7 and 25 Mbp. These regions segregate independently from the one located on chromosome 1, are not restricted geographically and show no clear association with morphological traits. The two inversions are strongly linked, as the genotypes across all quails remain consistent in both regions. This could suggest the persistence of errors in the assembly. However, the differences in heterozygosity, divergence time and ratio of nonsynonymous to synonymous mutations (N/S) and previous cytogenetic results support the hypothesis of their origin from a double inversion event. The suppressed recombination in the inverted regions reduces selection efficiency, leading to ineffective purging and increase in deleterious mutations. This effect is evident in the inversion in chromosome 1, where Ne is reduced due to its restricted geographical distribution, and in the largest inversion in chromosome 2. The N/S values vary along the genome, indicating the uneven action of selective forces in different segments of the genome. Chromosomal inversions evolve independently from the collinear genome and can accelerate the divergence in this bird lineage. CHAPTER 1 BIASED ASSESSMENT OF ONGOING ADMIXTURE USING STRUCTURE IN THE ABSENCE OF REFERENCE SAMPLES ADAPTED FROM: Sara Ravagni, Ines Sanchez-Donoso & Carles Vilà (2020). Biased assessment of ongoing admixture using STRUCTURE in the absence of reference samples. Molecular Ecology Resources, 21: 677–689. https://doi.org/10.1111/1755-0998.13286 CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 51 | PA G E ABSTRACT Detection of hybridization and introgression is important in ecological research as in conservation and evolutionary biology. STRUCTURE is one of the most popular software to study introgression and allows estimating what proportion of the genome of each individual belongs to each ancestral population, even in cases where no reference sample from the ancestral nonadmixed populations is previously identified. In spite of its frequent use, some studies have indicated that ancestry estimates may not always be reliable. We simulated population data under different conditions with regard to the genetic differentiation between ancestral populations, number of loci considered, number of alleles per marker and hybridization rate, and analyzed data with STRUCTURE. When reference samples were not included, the comparison of the known degree of admixture for each simulated individual and the value estimated with STRUCTURE revealed a strong underestimation of the level of introgression, classifying many admixed individuals as nonadmixed. This derives from an inaccurate estimation of the ancestral allele frequencies. When samples from the nonadmixed ancestral population were included as reference in the analyses, the bias in the estimations was reduced. The most accurate estimates were obtained when potentially admixed samples were few in relation to reference samples. Thus, whenever possible, a very large proportion of nonadmixed reference samples should be included in admixture assessments and different approaches should be combined. The misestimate of the amount of introgression can impair our understanding of the evolutionary history of species and misguide conservation efforts. KEYWORDS Bayesian clustering, hybrid zone, hybridization, introgression, population admixture, simulation CHAPTER 1 52 | PA G E INTRODUCTION Hybridization is a major concern in conservation biology (Rhymer & Simberloff, 1996) but also a source of evolutionary innovation (Abbott et al., 2013; Arnold, 2015) and adaptive introgression (the introduction of genes from one evolutionary lineage into the gene pool of another) has been shown to play an important role in the evolution and diversification of different clades (Burgarella et al., 2019; Hamilton & Miller, 2016; Oziolor et al., 2019; Suarez-Gonzalez et al., 2018). Thus, the identification of admixed individuals and the assessment of the levels of introgression are fundamental steps to study the effect of hybridization on populations. However, the detection of these admixed individuals can be problematic when the parental species show limited phenotypic differentiation (Randi, 2008). The ineffectiveness of morphological criteria in differentiating cryptic hybrids or admixed individuals (Oliveira et al., 2008) promoted the use of highly polymorphic genetic markers, such as microsatellites or SNPs. The development of software programs based on the analysis of panels of genetic markers and the implementation of Bayesian clustering models has facilitated the detection of hybrids and the assessment of introgression of genes of one species into the gene pool of another. One of the most widely used programs for these analyses is STRUCTURE (Pritchard et al., 2000) which allows identification of the origin of the genome of each individual. After choosing a K value, i.e., the number of populations, STRUCTURE subdivides the sample into K different clusters - trying to minimize departures from Hardy- Weinberg equilibrium and linkage disequilibrium - and estimates for each individual the proportion of the genome that could originate from each cluster (Barilani et al., 2007). In this way, STRUCTURE is very efficient at separating populations and identifying admixed individuals. A Scopus search (on 26 June 2019) for articles citing the original paper of Pritchard et al. (2000) and containing the words ‘hybrid*’ or ‘introgression’ in the title, abstract or keywords returned a total of 3,446 articles, highlighting the popularity of the software to detect hybridization and population admixture. However, the power of STRUCTURE to accurately estimate the amount of introgression has been questioned as it could underestimate the number of admixed individuals (Randi, 2008; Sanchez-Donoso et al., 2014; Sanz et al., 2009). The main limitation seems to be related to the detection of old admixture as opposed to F1 hybrids and first generations backcrosses (Oliveira et al., 2008). It has been suggested that estimates of the amount of introgression could be improved by increasing the number of loci (Pritchard et al., 2000; Randi, 2008; Vähä & Primmer, 2006) and using both linked and unlinked markers (Lecis et al., 2006), as this could enable a better CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 53 | PA G E assignment of admixed individuals in separate genotypic classes and the identification of past events of hybridization. The use of high number of loci derived from genome-wide data increases resolution in the detection of varying levels of introgression by identifying regions of the genome of different origin (Gómez- Sánchez et al., 2018). Unfortunately, this is not always feasible due to the lack of a reference genome, or due to a trade-off between costs and number of specimens to analyse. In order to characterize population variability, large numbers of samples may need to be studied, making genomic approaches unaffordable and, despite technical advances, studies with reduced number of loci continue to be frequent in day-to-day studies of hybridization and admixture (for example, see Alacs et al., 2010; Arias et al., 2019; De Barba et al., 2017; Sujii et al., 2019). In addition, large numbers of unmapped markers are not necessarily an advantage for the study of ongoing hybridization because they could yield redundant information if they are in complete linkage. Another factor affecting the performance of STRUCTURE is the sampling scheme, i.e. whether or not samples of each parental population or species are similar in size, and the phylogenetic relationships between the two species (Neophytou, 2014; Puechmaille, 2016). A previous study using empirical data of populations with known ancestry showed that STRUCTURE outperforms other common approaches in the identification of admixed individuals (Bohling et al., 2013). However, this study emphasizes the importance of including in the analyses a portion of individuals that can be a priori diagnosed as nonadmixed for the two parental classes and this may not always be feasible. In many studies genetic analyses are carried out to identify potentially admixed individuals without previously defining reference samples (for example, see Godinho et al., 2011; Muñoz-Fuentes et al., 2007; Oliveira et al., 2008; Ortego et al., 2017; Trigo et al., 2013). This can be particularly important in cases with high levels of hybridization and persistent introgression throughout the distribution range, or in cases where introgression takes place in geographically structured populations for which reference samples from another population may not be appropriate (for example, see Glover et al., 2017; Lavretsky et al., 2019; Sullivan et al., 2016). In this study we used simulations to determine the accuracy of STRUCTURE in the estimation of the individual level of introgression when nonadmixed reference individuals were not available. We assessed the importance of the number and polymorphism of the loci used, the divergence between the ancestral populations and different rates of hybridization between them. We also studied how adding appropriate reference samples impacts the reliability of the estimates. Our goal was CHAPTER 1 54 | PA G E to identify ways to improve the accuracy of the estimates of the degree of introgression. MATERIALS AND METHODS We carried out simulations of asymmetric gene flow from one population to another to reduce complexity and to facilitate analyses because the allele frequencies for one of the populations would not change over time. However, this is not an uncommon situation. Some examples of this kind of gene flow are the hybridization between domestic and vulnerable wild species (Godinho et al., 2011), admixture between a rare species and an abundant one (An et al., 2017), restocking of game species with farmed animals of alien origin (Sanchez-Donoso et al., 2014), or directional gene flow as a result of the biology of the species hybridizing (Muñoz-Fuentes et al., 2007). SIMULATION OF ANCESTRAL POPULATIONS We used the software EASYPOP v. 2.0.1 (Balloux, 2001) to simulate pairs of diploid random mating populations, genotyped for 200 unlinked loci assuming different levels of polymorphism, with (a maximum of) either two, five or 10 alleles per locus. Mutation rate for markers with five or 10 alleles was set at 10–3 and for markers with two alleles at 10–8, as commonly assumed for microsatellites and SNPs (Drake et al., 1998; Ellegren, 2004; Payseur & Nachman, 2000). We assumed a mutation model (K- allele model: KAM) for markers with two alleles: each allele had the same probability to mutate to the other allelic state. We used a mixed model including single step mutation model (SSM) for markers with up to five or 10 alleles, with a proportion of 0.3 of KAM events (Ellegren, 2004). We generated between 700 and 5,000 individuals per generation and ran the simulations for at least 1,000 generations to assume a long time of separation between the populations and approach mutation-drift equilibrium and stable differentiation. We generated 100 pairs of populations for each kind of marker and with genetic differentiation, measured as FST, around 0.05, 0.1 and 0.2. Thus, we simulated 900 pairs of populations genotyped for 200 loci each (3 kinds of markers × 3 levels of differentiation × 100 replicates = 900 pairs of populations). Parameters used for the runs are reported in Table S1. For details of the simulations see Supporting Information, Figures S1 and S2. CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 55 | PA G E SIMULATION OF ALLELE INTROGRESSION We wrote a script in Python 2.7 to simulate different hybridization rates and to subsample a random subset of individuals from each population to be analysed with STRUCTURE (see Figure 1 for a schematic representation of the simulations). Population A was where admixture took place due to some individuals arriving from population B. The goal of the analyses was to assess if it was possible to correctly estimate the degree of introgression in individuals sampled from the population A using STRUCTURE. Hybridization rate, that is the proportion of breeders originating from population B that contributed to the offspring of population A every generation, was around 1% (0.01) or 5% (0.05). Although high, rates of hybridization this high have been described for diverse taxa (for example, see Lavretsky et al., 2019; Muñoz- Fuentes et al., 2007; Nussberger et al., 2014). A key factor here is that we assume recurrent hybridization every generation and that admixed individuals do not have a reduced fecundity so that they are able to freely interbreed with other individuals in population A. To initiate the simulations (at generation 0), the program randomly selected 1,000 individuals from populations A and B from one of the 100 initial pairs of populations generated with EASYPOP for a given combination of FST and type of marker. For 10 generations, the programme generated the same number of individuals (1,000) by selecting two parents from the previous generation. Each parent was randomly selected from population A with a probability of (1 − m), where m was the hybridization rate (0.01 or 0.05); otherwise, the parent was randomly chosen from population B (geneflow was unidirectional from B to A). The genotype of the offspring was obtained by randomly selecting one of the two alleles from each parent at each locus with equal probability; loci were independent. For each individual, we calculated the proportion of the genome belonging to population A (qreal): for individuals originating from population B, this proportion was qreal = 0.0; for individuals from population A at generation 0, this proportion was qreal = 1.0; for subsequent generations, the proportion was the average of the values of the two parents. To confirm that the script was working as expected, we calculated the level of introgression expected after 10 generations in the different scenarios as in Verdu and Rosenberg (2011). Indeed, the value corresponded to the average level of introgression per individual in population A. At the tenth generation, a random subsample of 100 individual genotypes was taken both from population B and from the introgressed population A to be analysed with STRUCTURE. In these subsamples, the number of loci was reduced to 10, 30 or 100 loci chosen at random. The data were used to generate an input file for STRUCTURE. CHAPTER 1 56 | PA G E FIGURE 1. Pipeline to simulate introgression and subsampling for STRUCTURE analyses. First, 100 pairs of populations for each set of cases (pairwise FST = 0.05, 0.1 or 0.2, typed at 200 marker loci with two, five or 10 alleles) were simulated with EASYPOP. For 10 generations in population A, we randomly selected pairs of genotypes from the previous generation in population A or from population B according to the hybridization rate, and generated genotypes of 1,000 offspring. Finally, we randomly subsampled 100 individuals from each population A and B (with genotypes for 10, 30 or 100 randomly chosen loci) and genotypes were analysed in STRUCTURE. From the output of this program, we extracted estimates of ancestry and compared to the real values derived from the proportion of ancestry from population A in the simulations (see text). CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 63 | PA G E We used generalized linear models to assess both the effect of the degree of differentiation between the ancestral populations (FST) and the number of alleles per marker on the absolute difference between qSTRUCTURE and qreal (the bias in the inference of q) for the two hybridization rates and considering 100 loci. Both variables proved to have a highly significant effect (p < 10–16). The increase in genetic differentiation showed the stronger effect in reducing the difference between the two q values (Table 1), therefore improving STRUCTURE estimates, but having markers with more alleles also helped. TABLE 1. Generalized linear models testing the effect of FST and number of alleles over the absolute difference between qreal and qSTRUCTURE for samples from the admixed population A. Response variable Hybridization rate Explanatory variables Estimates z p |qSTRUCTURE - qreal| 1% (intercept) -2.541 -300.71 <2e-16 *** FST -1.007 -26.93 <2e-16 *** Number of alleles -0.014 -14.72 <2e-16 *** |qSTRUCTURE - qreal| 5% (intercept) -0.838 -109.26 <2e-16 *** FST -0.943 -27.98 <2e-16 *** Number of alleles -0.0129 -15.4 <2e-16 *** Note: The two factors have a significant effect. We used 100 markers with five or 10 alleles and the two rates of hybridization. The degree of overestimation of q by STRUCTURE can be better appreciated in the distribution of qSTRUCTURE − qreal for population A (Figure 3, S6 and S7). While the overestimation was limited when the hybridization rate was 1%, it dramatically increased when the rate was 5%. We compared the distribution of qreal and qSTRUCTURE for an example case (FST = 0.1, 30 markers with 10 alleles) for both hybridization rates (Figure 4). Although the two hybridization rates resulted in very different populations with regard to the level of introgression (see qreal in Figure 4), the results obtained with STRUCTURE were almost identical (see qSTRUCTURE), suggestive of a relatively low introgression, and showing that STRUCTURE was unable to differentiate the two cases. CHAPTER 1 64 | PA G E FIGURE 3. Density plots of qSTRUCTURE − qreal for individuals from population A. Simulations were carried for a differentiation of FST = 0.1 between the ancestral populations. (a) hybridization rate of 1%; (b) 5%. Data should ideally distribute around value 0 (marked with a vertical line). However, STRUCTURE analyses tended to overestimate the proportion of genome belonging to population A (qSTRUCTURE − qreal), being the overestimation bigger with higher hybridization rate. FIGURE 4. Density plots for qreal (a) and qSTRUCTURE (b) resulting from hybridization rates of 1% and 5%). Simulations were carried for a degree of differentiation between the ancestral populations (FST) of 0.1 and 30 loci with 10 alleles. The values around q = 0 corresponded to individuals from population B, while the rest reflect the admixed population A. Although qreal indicated different biological situations for the two hybridization rates (a), with very different number of nonadmixed A individuals (q close to 1), qSTRUCTURE values suggested that the two situations resulted in similar introgression (b). CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 65 | PA G E The allele frequencies estimated by STRUCTURE for population A were progressively diverging from the ancestral as the number of generations of introgression increased (Figure 5). At generation 0, before any admixture, the estimates of allele frequencies for populations A and B were similar to the frequencies calculated from the ancestral populations. The distances in this case were not 0 because the estimates obtained by STRUCTURE were based on subsamples of the populations, resulting in sampling errors. As introgression increased the number of alleles from population B into population A in subsequent generations, STRUCTURE estimates of the allele frequencies in population A diverged more strongly from the ancestral ones, while estimates of the allele frequencies for population B remained unaffected. As estimations of q are associated to the inferred allele frequencies, the inaccuracy in the ancestral allele frequency estimates could lead to the observed overestimations in qSTRUCTURE for individuals originating from population A. FIGURE 5. Distance between estimates of ancestral allele frequencies obtained with STRUCTURE and their real values for populations A and B as introgression of alleles from B to A advanced. A total of 100 STRUCTURE estimates of ancestral allele frequencies were obtained at generations 0, three, six and 10. In population B, estimates corresponded closely to the real values in all generations. For the admixed population A, divergence from the ancestral allele frequencies increased with introgression showing that STRUCTURE did not correctly estimate ancestral allele frequencies. CHAPTER 1 66 | PA G E A possible solution to improve accuracy in the estimates of the allele frequencies could be the inclusion of reference individuals from the ancestral population A. Generalized linear models showed that (with a hybridization rate of 5%) the accuracy of qSTRUCTURE estimates improved by increasing the proportion of reference individuals and the number of alleles per marker, as well as activating the POPALPHAS option and marking reference individuals with USEPOPINFO (Table 2; Figure 6 and S8). However, qreal and qSTRUCTURE were still very different (Figure 6). When the hybridization rate was lower (1%) the effect of including reference samples was not obvious, and the generalized linear model could not be fitted because the model assumptions were not fulfilled. TABLE 2. Generalized linear models explaining the differences between qreal and qSTRUCTURE for samples from the admixed population A. Response variable Hybridization rate Explanatory variables Estimates z p |qSTRUCTURE - qreal| 5% (intercept) -0.863 -140.89 <2e-16 *** Reference -3.432 -205.56 <2e-16 *** Number of alleles -0.014 -19.32 <2e-16 *** USEPOPINFO (activated) -0.133 -33.92 <2e-16 *** POPALPHAS (activated) -0.499 -140.03 <2e-16 *** Note: The explanatory variables were the proportion of samples used as a reference (0%, 10% or 30%), number of alleles of the markers (five or 10) and use of the USEPOPINFO and POPALPHAS options in STRUCTURE for simulations using 30 loci and a hybridization rate of 5%. The strongest effect is associated to the proportion of individuals used as reference (an increase in the proportion of individuals used as reference leads to a decrease in the difference between the estimates). Despite the improvement in the estimates with the inclusion of reference individuals, the biases were still very apparent even in the case when 30% of the individuals from A corresponded to reference samples (Figure 6) and qSTRUCTURE values were still significantly higher than the corresponding qreal (p < 10–16). This could be due to inherent biases in STRUCTURE or difficulties in the inference of the ancestral frequencies when a large proportion of the samples derived from the admixed population. To investigate which was the case, we compared qreal and qSTRUCTURE obtained when analysing small sets of 5 admixed individuals with 100 samples from A before admixture and 100 samples from B. The results (Figure 7) show that the biases in the estimates of qSTRUCTURE practically disappeared. CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 67 | PA G E FIGURE 6. Comparison of qSTRUCTURE and qreal when 30% of the individuals from the target population are sampled from the ancestral population and are used as reference. Individuals used as reference and those belonging to population B were excluded from the plots. Simulations were carried out with FST = 0.1 and 30 loci, varying the number of alleles per marker and the hybridization rate. The accuracy of qSTRUCTURE estimates improved notably with the inclusion of reference individuals (compare to Figure 2b) and when activating the POPALPHAS option, especially for the cases with higher hybridization rate, where no pure individuals remained in the admixed population. Variance of qSTRUCTURE estimates was quite large when using 30 markers but centred around the corresponding qreal values (along the diagonal in the figures), but still qSTRUCTURE values were significantly higher than qreal (p < .003). As expected, using 100 markers greatly reduced this variance and with 10 alleles qSTRUCTURE and qreal values were not significantly different (p-value = .510). These results imply that STRUCTURE estimates were not intrinsically biased and including a very large number of reference samples and high number of markers helped to reduce biases in the estimates. CHAPTER 1 68 | PA G E FIGURE 7. Comparison of qSTRUCTURE and qreal when most of the analysed individuals were used as reference and admixture was evaluated in just a few individuals. Simulations were carried out for a differentiation of FST = 0.1 between the ancestral populations. Each panel represents 10 independent qSTRUCTURE estimates for 100 admixed individuals representing the full range of qreal values. In each STRUCTURE run, 100 individuals from the ancestral population A and 100 from population B were analysed together with five admixed samples. The points lie around the diagonal showing that qreal and qSTRUCTURE tended to be similar. The analysis of a data set from a natural population of common quails that experienced introgression from farm quails showed the same pattern when we added reference samples belonging to the same population but from before most of the restocking campaigns (Figure S9). The q values estimated after adding reference showed more admixture than those estimated without adding a proper reference for the wild population (Wilcoxon signed rank test, p < 10–6). This result confirmed the pattern observed in the simulations. Given the biases observed with STRUCTURE, we carried out additional analyses with simulated data corresponding to ancestral populations differing by FST = 0.1 and 100 biallelic loci with ADMIXTURE, Ohana and sNMF to assess if the same biases were present in all cases. We graphically compared results from the four programs and STRUCTURE exhibited the worst performance under the set of conditions that we evaluated (Figure S10). The Wilcoxon signed rank test confirmed that q values were not overestimated with ADMIXTURE, Ohana and sNMF when the hybridization rate was 1%. CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 69 | PA G E DISCUSSION Our results show that when appropriate reference samples are not included in the analyses, ancestry estimates provided by STRUCTURE can be very biased. In fact, the results provided by this software can be very similar even when comparing populations with completely different levels of introgression (Figure 4). When populations experience hybridization and admixture during multiple generations, the proportion of the genome deriving from the ancestral local population tends to be overestimated by STRUCTURE, leading to an underestimation of the real degree of introgression and of the number of admixed individuals. Although the precision of the ancestry estimates provided by STRUCTURE tends to improve with higher number of markers and alleles (as suggested by previous studies; McFarlane & Pemberton, 2019; Vähä & Primmer, 2006), this increase in precision does not correspond to an increase in accuracy and similar biases are observed using a small or a larger number of markers. The overestimate is particularly extreme in scenarios of hybrid swarms, as simulated with a hybridization rate of 5%: while none of the individuals was free of introgression, the estimation offered by the program suggested that purebred individuals were majority in the sample. The poor performance of STRUCTURE under the simulated scenarios could impact the interpretation of admixture patterns in deeply introgressed populations, affected by various generations of hybridization, as in contact zones (Baldassarre et al., 2014; Johnson et al., 2015; Ortego et al., 2018). Extensive admixture could also occur in other cases, such as in the intercrossing between wild and domestic species, which have been coexisting for centuries or decades, or between different wild species that become in contact after a long time of evolution in isolation (Beaumont et al., 2001; Burgarella et al., 2018; Glover et al., 2017; Mckelvey et al., 2016; Scarcelli et al., 2017). It may seem that a hybridization rate of 5%, as in our models, is unlikely to exist in nature. However, similar values have been reported in the literature based on the identification of F1 hybrids (Lorenzini et al., 2014; Muñoz-Fuentes et al., 2007; Pacheco et al., 2017; Sullivan et al., 2016). The limitations and the poor performance of STRUCTURE in admixed populations were already in part highlighted in the recommendations of Pritchard et al. (2000) about a proper utilization of the software to obtain reliable ancestry estimates. These authors indicate that, in cases of extensive admixture, STRUCTURE cannot estimate ancestral allele frequencies and it cannot give accurate estimates of q because of the high variance in how many of the individual's alleles derive from one or the other CHAPTER 1 70 | PA G E population. Our results confirm this and show that introgression leads to poor estimates of the ancestral allele frequencies, as we observe how the allele frequencies estimated for the ancestral population were increasingly different from the real ancestral frequencies (Figure 5), reflecting the allele frequencies for an already admixed population. The biases in the estimates of ancestry persist even in the cases with 100 loci of high polymorphism (Figure 2), suggesting that the estimation of ancestral allele frequencies may not be corrected just by increasing the number of unlinked loci. After repeated introgression during several generations, STRUCTURE may be unable to reliably reconstruct ancestral allele frequencies. Including purebred reference individuals in the data set resulted in a larger improvement in the estimates of q than just increasing the number of markers. However, even replacing 30% of the individuals from the admixed population A by reference individuals was not enough to completely remove the bias (Figure 6) and the best estimates were obtained when analysing just a few target individuals together with many reference individuals (Figure 7). In order to develop cost-effective genetic tools for the assessment of introgression, efforts have been concentrating in the identification of the most suitable combination of markers. A reduced number of informative loci with high diagnostic power could be as effective as a high number of less informative loci, indicating that the discriminating power could be more important than their number (Oliveira et al., 2015; Randi et al., 2014). However, it is not completely clear what would be the best strategy to identify the most informative loci without previous data on marker variability across populations. On the other hand, in our simulations we considered all loci to be unlinked. The use of linked markers could also improve the identification of older admixture between populations (Falush et al., 2003; Lecis et al., 2006) because introgression can differently affect regions of the genome (for example, see Anderson et al., 2009), but may be less suitable to study ongoing hybridization and introgression. Our study highlights the importance of carrying out simulations in each study case to assess the reliability of estimates. The level of introgression can be assessed combining different approaches and simulations trying to properly reflect the functioning of the study system (McFarlane & Pemberton, 2019) should be carried out to test the power of the analysis, determine accuracy and assess possible biases (for example, see Oliveira et al., 2008; Randi et al., 2014; Sanchez-Donoso et al., 2014 Sanz et al., 2009). The estimates by STRUCTURE show a remarkable increase in accuracy when many nonadmixed reference individuals are included in the analysis (Figure 6). Therefore, we stress the importance of including samples from reference nonadmixed CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 71 | PA G E populations (for example, museum specimens dating from before an admixture event; see also Sanchez-Donoso et al., 2014) to increase the reliability of STRUCTURE analyses. The availability of these samples can be limited, and historical DNA extraction and amplification can be costly and effort-demanding, but the accuracy in the analysis improves notably, increasing the reliability of the results. When only a limited number of reference samples is available, it could be useful to carry on multiple STRUCTURE analyses including only a small proportion of those samples whose ancestry is unknown (Figure 7), using USEPOPINFO to define the reference samples and comparing the results obtained with different test sample sets. Also, the use of the option POPALPHAS can improve the analyses when source populations are unequally represented and if there is unbalanced sampling (see Wang, 2017). The use of a small number of markers is also known to influence the results of STRUCTURE (Toyama et al., 2020). Nevertheless, Lawson et al. (2018) have shown that different demographic histories can lead to identical results suggesting admixture in STRUCTURE and emphasize the importance of combining analytical approaches to obtain a more robust analysis of recent demographic history. Our comparison of the results provided by different programs showed that not all of them suffer the same biases or to the same degree. Consequently, we strongly suggest to combine STRUCTURE with other approaches and simulations, and evaluate the consistency in the results, especially when it is not possible to include suitable reference samples in the analyses. This is especially important, for example, in the cases where hybridization and introgression can be relevant in the design of management and conservation plans. In ecology, conservation and evolutionary biology, it is important to efficiently identify hybrids and admixed individuals, as well as to determine gene flow among different populations. STRUCTURE analyses showed a tendency to classify admixed individuals as nonadmixed when reference samples were not included. The misidentification of the degree of introgression can impact our understanding about the evolutionary history of species or the risks of genetic homogenization and extinction, and therefore its implications for management and conservation plans should be carefully considered. ACKNOWLEDGEMENTS We thank Carlos Rafii for providing access to servers to run simulations and analyses. Part of the computer work was carried out in Genomics servers of the Doñana's Singular Scientic-Technical Infrastructure (ICTS-RBD). The Conservation and Evolutionary Genetics Group at the Estación Biológica de Doñana provided valuable CHAPTER 1 72 | PA G E support and comments on the manuscript and Dr Jennifer Leonard also helped us to improve the English. The study is supported by project CGL2016-75227-P to CV from the Spanish Government and an FPI (Formación de Personal Investigador) fellowship (BES-2017-081291) to SR. We are also grateful to the anonymous referees and associate editors for their constructive comments. AUTHOR CONTRIBUTION C.V., and I.S.-D. designed the study, S.R. carried out the project and wrote the scripts, S.R., and I.S.-D. analysed the data, S.R. wrote the first draft of the manuscript and all authors contributed to the text. This study was initiated as part of the Master in Biodiversity and Conservation Biology of S.R. at the University Pablo de Olavide (Seville, Spain). DATA AVAILABILITY STATEMENT Genotypes corresponding to simulated starting populations, scripts used to simulate introgression and sampling, and output of these scripts that were used for subsequent analyses (as well as R scripts used to analyse and plot data) are available at https://github.com/sararvg/introgression_structure. Quail data is available at the original publication (Sanchez-Donoso et al., 2012). REFERENCES Abbott, R., Albach, D., Ansell, S., Arntzen, J. W., Baird, S. J. E., Bierne, N., Boughman, J., Brelsford, A., Buerkle, C. A., Buggs, R., Butlin, R. K., Dieckmann, U., Eroukhmanoff, F., Grill, A., Cahan, S. H., Hermansen, J. S., Hewitt, G., Hudson, A. G., Jiggins, C., … Zinner, D. (2013). Hybridization and speciation. Journal of Evolutionary Biology, 26(2), 229–246. https://doi.org/10.1111/j.1420- 9101.2012.02599.x Alacs, E. A., Georges, A., FitzSimmons, N. N., & Robertson, J. (2010). DNA detective: A review of molecular approaches to wildlife forensics. Forensic Science, Medicine, and Pathology, 6, 180–194. https://doi.org/10.1007/s1202 4- 009-9131-7 Alexander, D. H., & Lange, K. (2011). Enhancements to the ADMIXTURE algorithm for individual ancestry estimation. BMC Bioinformatics, 12(1), 246. https://doi.org/10.1186/1471-2105-12-246 An, M., Deng, M., Zheng, S.-S., Jiang, X.-L., & Song, Y.-G. (2017). Introgression threatens the genetic diversity of Quercus austrocochinchinensis (Fagaceae), an endangered oak: A case inferred by molecular markers. Frontiers in Plant Science, 8, 229. https://doi.org/10.3389/fpls.2017.00229 CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 79 | PA G E SUPPORTING INFORMATION Supplementary text. Simulation of ancestral populations. Supplementary table S1. Input parameters used in EASYPOP. Figure S1. FST obtained in the population simulations. Figure S2. Final average number of alleles per locus in the simulated. populations versus the maximum number of alleles and in relation to level of differentiation (FST). Figure S3. Individual proportion of genome belonging to population A estimated with STRUCTURE (qSTRUCTURE) compared to the real value (qreal) calculated during the simulations, for FST = 0.05. Figure S4. Individual proportion of genome belonging to population A estimated with STRUCTURE (qSTRUCTURE) compared to the real value (qreal) calculated during the simulations, for FST = 0.2. Figure S5. Linear regressions for qSTRUCTURE in relation to qreal for individuals from population A. Figure S6. Density plots of qSTRUCTURE - qreal for individuals from population A resulting from the simulations with FST = 0.05 between the ancestral populations. Figure S7. Density plots of qSTRUCTURE - qreal for individuals from population A resulting from the simulations with FST = 0.2 between the ancestral populations. Figure S8. Comparison of qSTRUCTURE and qreal when 10% of the target sample are reference individuals from the ancestral population. Figure S9. Comparison between qreal and q estimates obtained with different software (for simulations with FST=0.1, and 100 loci with 2 alleles): ADMIXTURE, Ohana, sNMF and STRUCTURE. CHAPTER 1 | APPENDIX 80 | PA G E SUPPLEMENTARY TEXT Simulation of ancestral populations A high number of simulations for markers with different number of alleles were carried out in EASYPOP v. 2.0.1 (Balloux, 2001) in order to obtain pairs of populations with the desired levels of genetic differentiation. The parameters used in the simulations to reach the desired level of differentiation are indicated in Table S1. A total of 100 pairs of populations were simulated with each set of parameters. Although each simulation was different, the FST values observed in all cases were very close to the intended values and no outliers were detected (Figure S1). Consequently, all simulated datasets were appropriate to represent the targeted FST values. Similarly, since the EASYPOP simulations were run for thousands of generations, there was a risk that the final mean number of alleles per locus could be very much lower than the initial value due to drift. However, we simulated large populations and the average number of alleles remained very close to the initial number in all cases (Figure S2). Only in the case of simulations for FST=0.2 for markers with a maximum of 10 alleles, the final average number of alleles was clearly lower (8.04). However, this average remains very much higher than 5 alleles and the simulations remain relevant to represent sets of markers with high polymorphism. We are aware that our simulations may imply a high level of homoplasy due to the combination of a high mutation rate (especially for microsatellite-like markers, for which it was set at 10-3), elevated effective population size, and long simulations. Therefore, we carried out tests with different running parameters to test the potential effect of homoplasy on downstream analysis. On one hand, we simulated populations with lower effective size (N=400) and higher migration rate (0.001) and we run these simulations during a very high number of generations (50000) to reach driftmigration-mutation equilibrium (as in Evanno, Regnaut, & Goudet, 2005; Puechmaille, 2016). On the other hand, to eliminate completely the effect of homoplasy we assumed a mutation rate of 0 and a lower number of generations (600). In both cases, for a given value of population differentiation we obtained practically identical results, the same as for the analyses shown in the main text of the manuscript (data not shown), which shows that homoplasy did not have an effect on the analyses. CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 81 | PA G E REFERENCES Balloux, F. (2001). EASYPOP (Version 1.7): A Computer Program for Population Genetics Simulations. Journal of Heredity, 92(3), 301–302. https://doi.org/10.1093/jhered/92.3.301 Evanno, G., Regnaut, S., & Goudet, J. (2005). Detecting the number of clusters of individuals using the software STRUCTURE: A simulation study. Molecular Ecology. https://doi.org/10.1111/j.1365-294X.2005.02553.x Puechmaille, S. J. (2016). The program STRUCTURE does not reliably recover the correct population structure when sampling is uneven: Subsampling and new estimators alleviate the problem. Molecular Ecology Resources, 16(3), 608–627. https://doi.org/10.1111/1755-0998.12512 CHAPTER 1 | APPENDIX 82 | PA G E SUPPLEMENTARY TABLE Table S1. Input parameters used in EASYPOP to simulate populations with the desired level of FST for the different number of alleles. Mutation model: 1, KAM (each allele has the same probability to mutate to the other allelic state); 2, Mixed model of SSM (single step mutation model) with a 0.3 proportion of KAM mutation events. FST Number of alleles Number of generations Individuals per population Mutation rate Mutation model 0.005 2 1000 5000 10-8 1 0.1 2 1300 3000 10-8 1 0.2 2 1300 1400 10-8 1 0.005 5 1500 3000 10-3 2 0.1 5 3000 1500 10-3 2 0.2 5 4000 700 10-3 2 0.005 10 1500 3000 10-3 2 0.1 10 3000 1500 10-3 2 0.2 10 4000 700 10-3 2 CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 83 | PA G E SUPPLEMENTARY FIGURES Figure S1. FST obtained in the population simulations. Input parameters in Table S1 were chosen to target the expected FST values. Those parameters were used to obtain 100 pairs of populations for each combination of expected FST and number of alleles. The small variance in all cases indicates that all the simulated populations indeed represented the intended level of differentiation and were suitable for subsequent analyses. Figure S2. Average number of alleles per locus in the simulated populations. This was evaluated to confirm that the number of alleles in the simulated EASYPOP populations were close to the maximum number of alleles set in the simulations and had not declined excessively due to drift. Except for markers with 10 alleles and FST= 0.2, the effective number of alleles remained very close to the initial number set in the EASYPOP simulations. CHAPTER 1 | APPENDIX 84 | PA G E Figure S3. Individual proportion of genome belonging to population A estimated with STRUCTURE (qSTRUCTURE) compared to the real value (qreal) calculated during the simulations, for FST = 0.05. (a) Hybridization rate of 1%, (b) Hybridization rate of 5%. Compared to Figure 2, precision and accuracy in estimates decrease, showing that a lower genetic differentiation between the ancestral populations affects negatively the accuracy in the estimates of ancestry. (A) CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 85 | PA G E (B) CHAPTER 1 | APPENDIX 86 | PA G E Figure S4. Individual proportion of genome belonging to population A estimated with STRUCTURE (qSTRUCTURE) compared to the real value (qreal) calculated during the simulations, for FST = 0.2. (a) Hybridization rate of 1%, (b) Hybridization rate of 5%. Compared to Figures 2 and S3, this shows an improvement in accuracy and precision with increasing level of genetic differentiation between the ancestral populations. (A) CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 87 | PA G E (B) CHAPTER 1 | APPENDIX 88 | PA G E Figure S5. Linear regressions for qSTRUCTURE in relation to qreal for individuals from population A. We only show the cases where regression lines could be adjusted (residuals followed normal distribution). In all cases the intercept was significantly higher than 0, showing that STRUCTURE tended to overestimate the proportion of genome from population A. a, intercept; b, slope. All coefficients had a significant p value (<10-16). CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 95 | PA G E MATERIALS AND METHODS EXPERIMENTAL MODEL AND SUBJECT DETAILS The common quail is a Galliform with a wide distribution in the Palearctic region and in southern Africa. Eurasian populations winter in the Sahel and India (Guyomarc’h et al., 1998), reaching breeding grounds in northern Africa and Eurasia between February and March. It is a highly mobile species that combines these trans-Saharan migratory movements with nomadic movements during the breeding season in search of suitable but ephemeral habitats of mainly winter cereal crops, such as wheat and barley (Rodríguez-Teijeiro et al., 2009). As cereal maturation depends on latitude and altitude, destruction of this habitat due to harvesting forces quails to move in an asynchronous manner, following a ‘‘green wave.’’ After harvest, only females with broods remain in the breeding areas to raise their chicks (Rodríguez-Teijeiro et al., 2009). Chicks remain with their mothers for about 4 weeks (Mur, 1994) and females can lay several clutches, at very different locations, during the breeding season (Puigcerver, 1997). As sex ratio is unbalanced toward males, unpaired males move extensively during the breeding season in search of females (Rodríguez-Teijeiro et al., 2006). The common quail is a popular game species in many countries in Western Europe and ring recoveries confirm that philopatry is small and that quails move extensively across Western Europe (Rodríguez-Teijeiro et al., 1992). The population size of this species is estimated at 15-35 million individuals (BirdLife International, 2021). Sample information and all morphological and stable isotope data are included as supplemental datasets (Data S1A and S1B, respectively). Sampling and animal handling procedures conformed to relevant institutional and national guidelines and regulations in bioethics. All our experimental work was carried out in accordance with the ASAB/ABS Guidelines for the treatment of animals. Fieldwork in Spain and Morocco was conducted under the supervision of J.D. Rodríguez-Teijeiro and M. Puigcerver following the guidelines of the Federation of European Laboratory Animal Science Associations (FELASA), under sampling permits awarded by the regional governments. Portuguese samples were obtained under sampling permits 634/2017/CAPT and 635/2017/CAPT from ICNF (‘‘Instituto da Conservacão da Natureza e das Florestas’’) to David Gonçalves and Pedro Andrade. Sampling permit for Madeira was granted by ‘‘Secretaria Regional do Ambiente e Recursos Naturais,’’ 8/PNM/2015-MAD. Cape Verde samples were obtained under sampling permit 08/2016 from ‘‘Ministerio da Agricultura e Ambiente’’ to J.D. Rodríguez-Teijeiro. The work fulfilled the ethical recommendations of the European Union and the Spanish CHAPTER 2 96 | PA G E legislation (Spanish Law 5/1995 and Decree 214/1997) and was approved by the Ethics Committee on Animal Experimentation from the University of Barcelona. FIELD SAMPLING We captured common quail males during the breeding season in Spain, Portugal, Morocco, Italy, Netherlands and the Atlantic archipelagos of Madeira, Canary Islands and Cape Verde (Data S1A). We collected phenotypic data from each captured quail. We weighed the captured quails with a precision scale to the nearest ± 0.5 g and measured the width of their pectoral lipid band (viewed through the skin) with a caliper (±0.01 mm). We measured maximum folded wing length with a zero-stop ruler (±0.5 mm) by flattening the camber and straightening the lateral curvature of the wing, and the length of each primary feather and the second secondary from one wing with a pin ruler (±0.5 mm; Jenni & Winkler, 1989). We also took pictures of the throat to study pigmentation. In common quails, molt is split between the breeding area and the wintering area (Netwon, 2010), and some feathers are grown in each area. For each individual we inferred age from the patterns of molt and wear of the primaries (Saint Jalme & Guyomarc’h, 1995) and took feather samples for isotope analyses. The feathers were chosen taking into account the age category, and we sampled one feather grown in the breeding area and another one grown in the wintering grounds. About 100 µl of blood was extracted from each individual from the jugular vein for genomic analyses. After sampling, birds were released in the place of capture except for 15 individuals for which we sampled the testes for the immunofluorescence assays (Table S1). In this case, birds were immediately dissected and testicles were frozen in liquid nitrogen or dry ice. MOLECULAR GENETICS Genotyping by sequencing DNA was extracted by using the DNeasy Blood & Tissue Kit of QIAGEN. We genotyped quail samples using a reduced representation genotype-by-sequencing (GBS) approach, which consists of next generation sequencing short fragments of DNA (64 bp) starting at the cutting point of some restriction enzymes (Elshire et al., 2011). We built a shotgun library from DNA of 80 quails from 9 localities (plus a control) after digestion with the restriction enzyme EcoT22I. This library was sequenced using Illumina HiSeq 2500 (1x100bp) at Cornell University, Biotechnology Resource Center. We called SNPs using the TASSEL-3-GBS pipeline (Glaubitz et al., CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 97 | PA G E 2014) and mapped the obtained sequences against the Japanese quail genome (Coturnix japonica, RefSeq assembly accession number: GCF_001577835.1) using BWA-MEM v.0.7.8-r455 (Li & Durbin, 2009). We filtered for biallelic SNPs with minor allele frequency higher than 0.03, missing data less than 25%, minimum quality of 20 and depth higher than 5 and lower than 100 using VCFTOOLS v.0.1.13 (Danecek et al., 2011). After applying these filters, we retained 60,024 SNPs. We assessed population structure with a principal component analysis (PCA) with the function snpgdsPCA (Patterson et al., 2006; Price et al., 2006) in the R package SNPRelate (Zheng et al., 2012). To assess regions of the genome responsible for population differentiation, we calculated Weir and Cockerham’s FST (Weir & Cockerham, 1984) among the 9 localities for each SNP along the genome estimator implemented in VCFTOOLS. We plotted FST values with the package qqman (Turner, 2018) in R 4.0.3 (R Development Core Team, 2020), using RStudio v1.3.959 (RStudio Team, 2020). Analyses of population structure were carried out with ADMIXTURE v1.3.0 (Alexander et al., 2009). To determine the number of clusters (K) that best explain the data for different SNP datasets, we ran 20 independent analyses using the flag ‘‘K’’ from K = 1 to K = 5 with default parameters and we compared the cross-validation errors to find the most likely K. We combined the results from the different runs using the CLUMPAK server (Kopelman et al., 2015) with default settings. The initial dataset was complemented with 39 additional samples that were also typed using a genotyping-by-sequencing approach but utilizing the restriction enzyme Nsil, which is an isoschizomer of EcoT22I, and using an Illumina NovaSeq 2x150 bp sequencing at the University of Wisconsin-Madison Biotechnology Center. We optimized the filtering parameters, mainly the maximum proportion of missing data, to obtain a subset of SNPs that produced genotypes consistent with those obtained when analyzing the original dataset, which were used to genotype the new individuals. We used this smaller set of markers to increase sample size for the distribution of the karyotype variants. Whole genome sequencing To better characterize the chromosomal limits of the inversion, we sequenced the complete genome of 16 homokaryotype males at a coverage of about 10x. We selected 10 homokaryotypes from the southwestern type (BB) from separate locations on the continent and Atlantic islands, and 6 from northern and eastern quails (AA; Table S2). DNA was extracted by phenol chloroform followed by an CHAPTER 2 98 | PA G E ethanol precipitation. Whole genome sequencing libraries were constructed following standard protocols (Camacho-Sanchez et al., 2018), and sequenced at the SNP&SEQ Technology Platform in Uppsala (Sweden) using Illumina HiSeqX (2x150bp). IMMUNOFLUORESCENCE ASSAYS Spermatocyte nuclei surface spread To obtain spread nuclear chromatin, we performed surface spreading of spermatocyte nuclei on glass slides from frozen quail testis (Roig et al., 2010). Briefly, a small portion of frozen tissue was minced in cold PBS to obtain cell suspensions, cells were then treated with 1% lipsol for 40 min, fixed with fresh 1% paraformaldehyde containing 0.15% Triton X-100 for 2h in a humid chamber, airdried and washed with 0.4% Photo-Flo. Immunofluorescence staining To identify the chromosomes, surface spreads were stained against the axial core marker, SYCP3, using a previously described method (Pacheco et al., 2015): slides were blocked for 15 min with blocking buffer (0.2% BSA, 0.2% gelatin, 0.05% Tween- 20 in PBS) and then incubated with primary antibody antirabbit SYCP3 (1:400, Abcam Cat# ab15093, RRID:AB_301639), overnight at 4°C; after washes, slides were incubated with a goat anti-rabbit DyLight 649 (1:400, Jackson ImmunoResearch Lab Cat# 108-495-003, RRID:AB_2337515) for 1 h at 37°C. Finally, slides were mounted with Vectashield mounting medium containing DAPI. Slides were examined on a Zeiss Axiophot epifluorescence microscope. We captured 20-40 well-spread pachytene spermatocyte cells from each slide and saved their locations on each slide. Fluorescence in situ hybridization (FISH) assay Using the Japanese quail genome (C. japonica) as a reference (GENBANK ID GCA_001577835.1) and probes developed for chicken (Zlotina et al., 2012), we synthesized probes for FISH assays. One of the probes, K19 (NCBI code WAG- 22K19), was designed to hybridize to a DNA sequence inside the divergent region of chromosome 1 and the other one was outside this region (J23: WAG-22J23). Probes were labeled by nick-translation with digoxigenin-11-dUTP (Roche, for probe K19) or biotin-16-dUTP (Roche, J23). We performed an adapted FISH assay to detect the locations of the two probes on spermatocyte chromosomes (Pacheco et al., 2018). CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 99 | PA G E Slides were denatured in 70% formamide in 2% SSC for 20 min at 77°C, incubated with 1M sodium thiocyanate for 3h, denatured again with denaturation solution for 20 min at 77°C, and dehydrated in a series of ethanol solutions. At the same time, 200-400ng DIG or Biotin-labeled probes were precipitated with Cot-1 DNA (Invitrogen) and salmon sperm DNA (Stratagene). The precipitated probe DNAs were resuspended in hybridization buffer (50% formamide, 2X SSC, and 10% dextran sulfate) under constant agitation at 37°C, denatured for 10min at 79°C, and reannealed for at least 30min at 37°C. Finally, reannealed probe DNAs were incubated on denatured slides for 72h at 37°C. After hybridization, slides were washed with 50% formamide, 2 x SSC at 45°C, rinsed with 4 x SSC containing 0.1% Tween 20, and blocked with blocking buffer (4 x SCC, 4mg/ml BSA, 0.1% Tween20) for 30min at 37°C; finally, slides were incubated with Alexa Fluor 555 streptavidin (Invitrogen, 1:100) and Anti-Digoxigenin-Fluorescein (Roche, 1:10) for 1h at 37°C, washed and air-dried. Since SYCP3 signals were lost during the FISH assay, a second round of immunostaining against SYCP3 was performed after FISH. Then both probe and SYCP3 signals were examined and captured in pre-selected pachytene spermatocyte cells. Slides were examined using a Zeiss Axioskop epifluorescence microscope equipped with a Point Gray Research camera controlled by the ACO XY Software (A.COLOMA Open microscopy). Images were processed by Adobe Photoshop to match the intensity observed under the microscope. Measurement of distances between probes was performed using ImageJ (Rasband, 2018). STABLE ISOTOPE ANALYSES We washed the feathers with absolute ethanol, rinsed them with distilled water and put them in a stove to dry. We cut into small pieces two portions (0.3 mg each approximately) of the apical end of the feathers, and placed them into silver capsules for the δ2H analyses and into tin capsules for the δ13C analyses. δ2H quantification Isotope measurements were performed on H2 derived from high-temperature flash pyrolysis at 1450°C by means of Flash HT Plus elemental analyzer coupled to a Delta- V Advantage isotope ratio mass spectrometer via a CONFLO IV interface (Thermo Fisher Scientific, Bremen, Germany). Stable isotope ratio was expressed in the standard δ-notation (‰) relative to the Vienna Standard Mean Ocean Water (VSMOW). Based on laboratory standards, the measurement error was ± 3‰. The CHAPTER 2 100 | PA G E standards used were: CBS, KHS, USGS-42 (keratin standards supplied by Environment Canada) and LIE-PA2 (feathers of Razorbill, internal standard). The δ2H analyses were carried out using the comparative equilibration approach described in Wassenaar and Hobson (2003), and by using calibrated keratin isotope reference materials to avoid effects of H exchange with ambient water vapor. δ13C quantification Samples were combusted at 1020°C using a continuous flow isotope-ratio mass spectrometry system by means of Flash HT Plus elemental analyzer coupled to a Delta-V Advantage isotope ratio mass spectrometer via a CONFLO IV interface (Thermo Fisher Scientific, Bremen, Germany). Stable isotope ratios were expressed in the standard δ-notation (‰) relative to Vienna Pee Dee Belemnite. Based on laboratory standards, the measurement error was ± 0.1‰. The standards used were: EBD-23 (cow horn, internal standard), LIE-BB (whale baleen, internal standard) and LIE-PA (feathers of Razorbill, internal standard). These laboratory standards were previously calibrated with international standards supplied by the International Atomic Energy Agency (IAEA, Vienna). QUANTIFICATION AND STATISTICAL ANALYSIS Phenotypic data We calculated a modified Holynski Index (Holynski, 1965; Lockwood et al., 1998; Fielder, 2005) to measure wing pointedness based on the lengths of the primary feathers, the second secondary feather and the folded wing: 𝑀𝑜𝑑𝑖𝑓𝑖𝑒𝑑 𝐻𝑜𝑙𝑦𝑛𝑠𝑘𝑖 𝐼𝑛𝑑𝑒𝑥 =∑𝑃𝑟𝑜𝑥−∑𝐷𝑖𝑠𝑡 𝑊𝑖𝑛𝑔 ∗100 Being: ∑𝑃𝑟𝑜𝑥 =∑(𝑃𝑚−𝑃𝑟𝑜𝑥𝑃𝑚−𝑖)+(𝑃𝑚−𝑆2) 𝑦 𝑖=1 ∑𝐷𝑖𝑠𝑡 =∑(𝑃𝑚−𝐷𝑖𝑠𝑡𝑃𝑚+𝑗) 𝑧 𝑗=1 CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 101 | PA G E Pm: length of the largest primary feather (m) y: number of primaries proximal to m z: number of primaries distal to m ProxPm-i: length of the ith proximal primary feather in relation to m DistPm+j: length of the jth distal primary feather in relation to m S2: length of the second secondary feather Wing: length of the folded wing We statistically compared weight, the width of the lipid band and the modified Holynski Index among karyotype groups (AA, AB and BB) in one year old males (Euring 5) by fitting a linear model for each one of the response variables studied. Karyotype group was the explanatory variable in the three models. Homoscedasticity and normality of residuals were checked by visual exam of scatterplots (both were fulfilled). Significance was evaluated with an F-test with the function Anova() available in the R car package (Fox & Weisberg, 2019). Pairwise comparisons between genomic clusters were done with Tukey’s posthoc tests using the function emmeans() from the emmeans package (Lenth, 2020). Throat pigmentation was categorized into six levels, from pale (level 1) to completely pigmented (level 6, Figure S5), by visual examination of throat pictures of one year old males (Euring 5) and comparing to a pre-defined scale. In any case, to avoid subjective appreciation, the characterization of quails was always done by the same researchers (JDR-T and MP) that have been working together in the study of quails for decades. We performed a Kruskal-Wallis test to test the significance of the differences in throat pigmentation among groups using the function kruskal.test(), and used Wilcoxon-Mann-Whitney tests to assess pairwise differences between karyotype groups with the function wilcox.test(). Statistical analyses were done in RStudio. GENOMIC DATA SNP calling and filtering We checked sequencing data quality by using FASTQC v.0.11.5 (Andrews, 2010) and we removed adaptor sequences and quality trimmed with CUTADAPT v.1.17 (Martin, 2011). We mapped reads against the Japanese quail genome (RefSeq assembly accession number: GCF_001577835.2) using BWA-MEM v.0.7.15 (Li & Durbin, 2009) and we performed read group tagging and duplicate marking with CHAPTER 2 102 | PA G E PICARD v.2.4.1 (https://broadinstitute.github.io/picard/). We used GATK v3.6 (McKenna et al., 2010) for joint indel realignments. We used the depth option implemented in SAMtools v. 1.6.0 (Li et al., 2009) to calculate the mean depth of coverage. We called variants across the 16 genomes using FREEBAYES v.1.3.1 (Garrison, 2012), with the flags ‘‘-X’’ and ‘‘-u’’ to ignore complex events and multi-nucleotide polymorphisms. Variants were filtered from the initial variant call set via hard-quality filters in vcffilter (from vcflib) as follows: mapping quality across samples higher than 40, additional contribution of each observation at least 10 log units, reads on both strands and at least two reads balanced to each side of the site. To remove false variants, we excluded sites with quality score smaller than two times the depth when the depth was greater than the mean depth plus the squared root of the mean depth using VCFTOOLS v.0.1.15 (Puritz et al., 2014). Afterward, we removed loci with a mean depth higher than 3 times the mean depth and those in repetitive regions, based on the soft-masked Japanese quail genome. We kept only biallelic variants located on assembled chromosomes and with less than 50% missing data and minor allele count greater than 2. After applying quality and depth filters, we kept 17,836,922 biallelic SNPs. Population structure and differentiation We assessed population structure with a principal component analysis (PCA) using PLINK v1.90b6.10 (Purcell et al., 2007) and with ADMIXTURE and CLUMPAK. These analyses were performed with SPNs from whole-genome sequences, from the inverted region in chromosome 1 and also excluding the inverted region. After confirming the existence of a chromosomal inversion using immunofluorescence, in order to better identify the region of high differentiation, we calculated average FST between homokaryotype individuals (AA versus BB) across non overlapping 1 Kbp and 1,000 Kbp windows using the FST implemented in VCFTOOLS. We plotted FST values along the genome with the package ggplot2 (Wickham, 2009) in RStudio. To confirm the suppression of recombination in the inverted region, we thinned the SNP dataset for sites at a minimum distance of 10 Kbp and we calculated R2 values in chromosome 1 with the LDHeatMap package (Shin et al., 2006) in R for all samples and only for inverted and non-inverted samples. The putative position of the breakpoints was estimated from the location of the transition in FST values in the inverted region of chromosome 1. However, we considered the physical location to be slightly within the divergent region and not to CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 103 | PA G E correspond exactly with the position of the outlier SNPs, as the increased divergence can extend beyond the breakpoints due to genetic hitchhiking. To determine the position of the breakpoints, we visually checked the read mappings around the expected position of breakpoints with IGV v.2.5.239 (Thorvaldsdóttir, 2013), to identify differential patterns of pair orientations compatible with the presence of an inversion between inverted and non-inverted samples. We carried out pairwise alignments between regions including the putative location of breakpoints using the MAFFT aligner (Katoh & Standley, 2013; https://mafft.cbrc.jp/alignment/server/). The visual inspection of the read mapping and the pairwise alignments did not show clear breakpoints around the region of the transition of FST. This could be due to the use of a reference genome belonging to a different species or to the absence of the breakpoint region in the reference genome. These analyses suggested that the breakpoints could be located around positions 52,100,000 and 168,610,000 on chromosome 1. Based on a literature search, we identified 5 genes located in the inverted region that could be potential candidate for pigmentation: TYR, RAB-38, OCA2, DCT and GPR143 (Loftus et al., 2002; Protas & Patel, 2008; Bourgeois et al., 2016; El Sabry et al., 2017; Liu et al., 2019; Mastrangelo et al., 2020; Xi et al., 2020), and one related to migratory behavior or restlessness: NPAS2 (Steinmeyer et al., 2009; Mettler et al., 2015; Bazzi et al., 2017). To assess if the gene function could be different in quails with and without the inversion, we annotated the SNP dataset with SnpEff v.4.3t (Cingolani et al., 2012) using the GTF annotation file from the Japanese quail genome. We identified nonsynonymous mutations from the annotated SNP set to check if the gene function could have been affected by the inversion and we found one fixed non synonymous change between inverted and non-inverted chromosomal regions in RAB-38 and other three in TYR (Table S3). Time of divergence of the sequences in the inversion We estimated absolute divergence (dXY) between the inverted and non-inverted haplotypes, which is a measure of the fraction of nucleotide differences between the two sequences, using a custom perl script (https://github.com/MattChristmas/ Quail_inversion_scripts). As divergence is not consistent across the entire length of the inversion, being highest at the ends and decreasing in the middle likely due to double crossover and gene conversion events, we estimated dXY for the first and last 2 Mbp of the inversion. This gave us the best estimate of accumulated divergence since the origin of the inversion event. We then converted this divergence estimate to an estimate of divergence time using the equation: CHAPTER 2 104 | PA G E T = dXY/(2*µ) Where T is the divergence time in generations and µ is the mutation rate in mutations per site per generation. The mutation rate has not been estimated for the common quail so we estimated a range for T using published mutation rates for other bird species, including the ancestral bird lineage (1.23 x 10-9 site-1 generation-1), chicken (1.91 x 10-9 site-1 generation-1), zebra finch (2.21 x 10-9 site-1 generation-1) from Nam et al. (2010) and collared flycatcher (4.6 x 10-9 site-1 generation-1) from Smeds et al. (2016). With a generation time of ~1 year (Puigcerver et al., 1992), the value of T is directly translatable to an estimate of the number of years of divergence. STABLE ISOTOPE DATA We calculated the difference between the isotopic value of the primary feather grown during the previous breeding season minus the one of the primary feather grown during the last wintering season for each isotope. We fitted two linear models (one per isotope) to test whether the differences observed among genetic groups (fixed factor) were statistically significant. Homoscedasticity and normality of residuals were checked by visual exam of scatterplots (both were fulfilled). Pairwise comparisons between genomic clusters were done with the Tukey’s posthoc test using the function emmeans() from the emmeans package. RESULTS AND DISCUSSION To assess the genetic structure of common quails that differ in phenotype in the western limit of the distribution (Figures 1A and S1), we sampled 80 male quails from Italy, Spain, Portugal, Morocco, and the eastern Atlantic archipelagos of Canary Islands and Madeira, including localities where different color morphs were dominant. We characterized genome-wide variation across these individuals using a reduced representation approach (supplemental information) and identified a large block of SNPs along chromosome 1 that was associated with very high differentiation (FST) between localities (Figure 1B). We assessed genetic structure for this region of chromosome 1 by using ADMIXTURE assuming 2 clusters. While the genetic composition of many quails had 100% assignment to one cluster or the other, 30 had ~50:50 assignment to each cluster (Figure 1C), suggesting the existence of two chromosomal types. A principal-component analysis (PCA) of this region of CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 111 | PA G E We next investigated correlations between phenotype and karyotype. We used a subset of SNP markers to assess chromosomal composition of a larger sample of sexually active male quails from different locations and compiled morphological data. The three karyotype groups of quails differed in weight, with AA quails being the lightest (AA: mean ± SD = 94.8 g ± 6.9, n = 43; AB: 97.3 ± 6.8, n = 42; BB: 99.8 ± 6.4, n = 18; F2,100 = 3.73; p = 0.027; Figure 4A; see supplemental information). They also showed significant differences in pigmentation: BB quails had darker throat coloration, while it was intermediate for AB (Figures 4B and S5; median color category for BB was 5, 4 for AB, and 2 for AA; H(2) = 50.96; p < 0.001). Interestingly, the genes RAB-38 and TYR, which have previously been found to affect pigmentation in birds (Abolins-Abols et al., 2018; Mastrangelo et al., 2020), are very close to the inversion breakpoint and the genome sequences for AA and BB quails differ by several nonsynonymous substitutions in these genes (Table S3), suggesting that they could be associated with differences in pigmentation in this system too. We also observed differences in the width of the pectoral lipid band as viewed through the skin, which is positively related to motor restlessness (Saint-Jame, 1990). This band was significantly wider for AA quails than for AB and BB quails (AA: 6.58 mm ± 3.58, n = 43; AB: 2.98 ± 2.05, n = 42; BB: 2.44 ± 1.87, n = 18; F2,100 = 23.32; p < 0.001; Figure 4C). We used wing measurements to calculate a modified Holynski index (Fiedler, 2005) for each individual. Higher values for this index indicate that the wings are relatively narrower and pointier, which is expected to indicate higher flight efficiency (Mönkkönen, 1995). The comparison of this index among the three karyotypes showed that AA quails had significantly higher values, while BB had lower values and AB had intermediate values (AA: 116.0 ± 10.6, n = 21; AB: 104.9 ± 8.2, n = 38; BB: 98.4 ± 11.0, n = 16; F2,72 = 16.78; p < 0.001; Figure 4D), suggesting that the quails without the inversion (AA) had a more efficient flight apparatus. The limited flight efficiency of the quails with the inversion (AB and BB) could be associated with lack of or reduced migration. These phenotype patterns were not just an effect of latitudinal differences. Excluding AA quails from northern latitudes in the analyses gave the same results because AA quails did not differ in northern and southern localities (Figure S6). CHAPTER 2 112 | PA G E FIGURE 4. Morphology and isotopic signature in relation to the inversion (A) Body weight of AA, AB, and BB male common quails. (B) Proportion of individuals with different degrees of throat pigmentation. (C) Width of the pectoral lipid band. (D) Variation in the Holynski index, which measures wing pointedness and is related to flight efficiency. These morphological differences are not associated to latitude; see Figure S6. (E and F) Proportion of deuterium δ2H (E) and carbon δ13C isotopes (F) in two primary feathers from 1-year-old males, one grown during the previous breeding season and the other during the wintering season. AB and BB quails show small differences in the isotopic signature when comparing the two feathers, suggesting very reduced migratory movements compared to AA birds. ***p ≤ 0.001; **p ≤ 0.01; *p ≤ 0.05; ns, non-significant. Figure related to Data S1A and S1B. CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 113 | PA G E To study the migratory behavior, we analyzed the isotopic signature of primary feathers from captured individuals. The isotopic signature for carbon (δ13C) and deuterium (δ2H) of feathers is greatly determined by the frequency of these isotopes at the location where they were formed (Rubenstein & Hobson, 2004). We used feather wear and molting patterns to identify the age of quails. From 83 1-year-old males (Euring 5, the majority of the captured individuals), we took one primary feather that had been grown during the previous breeding season. This feather informed about the isotopic signature at the place of hatching. A second feather, grown in the wintering ground, allowed the comparison of the two locations. For the two isotopes, AA quails showed large differences between feathers grown in breeding and wintering grounds, while individuals that carried the inversion, either as homokaryotype BB or as heterokaryotype AB, showed very small differences (Figures 4E and 4F; supplemental information). This indicates that AB and BB birds did not migrate or that their migratory movements were much reduced as compared to quails without the inversion. Reduced dispersal associated with the inverted sequence is further supported by a neighbor-joining tree representing the evolutionary relationships between the 16 genomes considering only the region within the inversion. The 10 BB quails sequenced formed groups corresponding to different geographic origins: Madeira, Canary Islands, and Iberian Peninsula plus Morocco (Figure 3D; see also Figure 1D). This suggests reduced mobility for quails carrying B haplotypes that results in population structure. On the other hand, A haplotypes did not show population differentiation. At the same time, sequences outside the inversion do not show any population structure (Figures 3C and S2). These results imply that quails with and without the inversion can freely interbreed, and their genomes recombine over most of their length. However, those that carry one or two copies of the inversion (AB and BB) have reduced mobility, leading to population structure for the inverted sequence among geographic regions. Chromosomal inversions are known to play an important role in evolution (Kirkpatrick & Barton, 2006; Wellenreuther & Bernatchez, 2018) and can facilitate the coexistence of well-differentiated morphologies and behaviors (Küpper et al., 2016; Huang et al., 2018), including differences in migration patterns (Kirubakaran et al., 2016). In the case of the common quail, we have demonstrated that an inversion is associated with the coexistence of two complex sets of traits over part of the distribution range. The inversion is exceptionally large (Wellenreuther & Bernatchez, 2018), resulting in a supergene that represents a considerable portion of the quail genome, more than 12%. Several quail species inhabited the Macaronesian CHAPTER 2 114 | PA G E archipelagos and went extinct during the Late Holocene as a result of human impacts. Some of these species had larger body sizes and reduced flying capacities (Rando et al., 2020) and could have been the source of the inversion. The introgression of genes from a related species can facilitate adaptation (Kirkpatrick & Barton, 2006; Jay et al., 2018; Christmas et al., 2019), as in the case of the white-throated sparrow (Tuttle et al., 2016), where the introgression of an inversion has led to the coexistence of morphological and behavioral phenotypes (Merritt et al., 2020). However, we cannot exclude the possibility that the inversion in quails represents an ancestral polymorphism. The presence of a large number of heterokaryotypes (AB) implies the absence of complete pre-zygotic isolation between the two chromosomal types. Moreover, the lack of population structure excluding the inversion suggests the absence of postzygotic genetic incompatibilities. The restricted distribution of the inversion implies that heterozygotes also have reduced dispersal capacity and inter-island dispersal may have been the result of sweepstakes dispersal events (de Queiroz, 2005). This has been rare enough to result in genetic differentiation for the inversion, between islands and the mainland and between archipelagos. Patterns of genomic differentiation between AA and BB quails tend to decrease toward the center of the inversion (Figures 3A and S4), probably caused by occasional non-homologous pairing, which could allow some transfer of genetic material between the two chromosomal variants by double crossover events or by gene conversion (Navarro et al., 1997). These double crossovers are occasionally observed (Stevison et al., 2011) and may have contributed to an increase in genetic diversity in the region of the inversion in AA quails (Figure 3B) compared to the rest of the genome. We do not have evidence of differences in fitness between the two chromosomal types. However, the distribution of the darker-colored birds carrying the inversion is restricted to warmer areas, as expected under Gloger’s rule (Delhey et al., 2019). This could be associated with camouflage capacity or higher resistance of melanized feathers against keratinolytic bacteria, particularly in humid environments (Burtt & Ichida, 2004; Delhey et al., 2019). The habitats occupied by AB and BB quails at the southern limit of the breeding range may also offer suitable habitats for the species all year round, allowing permanent habitation. In fact, AB and BB individuals were found in areas where the common quail populations have been suggested to be resident (Figure1E; BirdLife International, 2021). The reduced recombination associated with the inversion facilitates the divergence between genes in this CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 115 | PA G E chromosomal region. This generates sufficient opportunity for the accumulation of substitutions that could facilitate local adaptation or even incompatibilities that could facilitate speciation (Faria & Navarro, 2010). Alternatively, frequency-dependent or temporally fluctuating selection due to changes in land uses induced by humans or climate (Kirkpatrick & Barton, 2006; Lamichhaney et al., 2016) could lead to stable polymorphisms over thousands of years. However, we know too little on the population dynamics of the two chromosomal types to understand the mechanisms of coexistence. ACKNOWLEDGEMENTS The authors thank all the people who helped with the fieldwork and especially Francesc Sardà-Palomera, Jorge Falagán, and Fernando Gavilán (Spain); Mohamed Maghnouj (Morocco); Jan Staal (the Netherlands); and Pierfrancesco Micheloni (Italy). Joan Navarro and Francisco Ramírez provided support for isotopic analyses. Irene Quintanilla, Anna Cornellas, Sarai López, and the Laboratories of Molecular Ecology (LEM-EBD) and Stable Isotopes (LIE-EBD) at the Doñana Biological Station (EBD-CSIC) helped with the laboratory work. Genome sequencing was performed by the SNP&SEQ Technology Platform in Uppsala. The facility is part of the National Genomics Infrastructure (NGI) Sweden and Science for Life Laboratory. The SNP&SEQ Platform is also supported by the Swedish Research Council and the Knut and Alice Wallenberg Foundation. Informatic facilities were provided by Doñana ICTS-RBD. Microproject (C.V.) was funded by Severo Ochoa award from the Spanish Government to EBD-CSIC; Spanish Government grants PID2019-108163GB-100 (C.V. and J.D.R.-T.), BFU2016-80370-P, and PID2019-107082RB-I00 (I.R.); Spanish Government fellowship BES-2017- 081291 (S.R.); China Scholarship Council fellowship (Y.H.); and Fundação para a Ciência e Tecnología Research Fellowship PD/BD/114028/2015 (P.A.). Federación de Caza de Euskadi provided financial support for part of the field sampling. AUTHOR CONTRIBUTION C.V. and I.S.-D. conceived the research. J.D.R.-T., M.P., I.J.-B., I.S.-D., C.V., P.A., and D.G. led the fieldwork, including sampling and phenotypic data collection. I.S.-D. and J.A.L. led stable isotope and molecular genetic lab work. I.R., Y.H., A.M.-L., and S.R. carried out immunofluorescence assays. S.R., I.S.-D., M.J.C., M.T.W., G.F., J.A.L., and C.V. carried out molecular genetic data analyses. I.S.-D., I.J.-B., and J.D.R.-T. carried CHAPTER 2 116 | PA G E out phenotype analyses. C.V., I.S.-D., and S.R. wrote the first version of the manuscript. All authors discussed and commented on the manuscript. DATA AND CODE AVAILABILITY Genomic data are deposited in GenBank (GenBank: PRJNA730394). Genotyping-by- sequencing data (non filtered SNPs) are accessible in https://digital. csic.es/handle/10261/251279 (https://dx.doi.org/10.20350/digitalCSIC/13989). Results of the FISH assays are included in Table S1 in this publication. Original script for dxy calculations is provided in https://github.com/MattChristmas/ Quail_inversion_scripts. Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request. REFERENCES Abolins-Abols, M., Kornobis, E., Ribeca, P., Wakamatsu, K., Peterson, M.P., Ketterson, E.D., & Milá, B. (2018). Differential gene regulation underlies variation in melanic plumage coloration in the dark-eyed junco (Junco hyemalis). Mol. Ecol. 27, 4501–4515. Alexander, D. H., Novembre, J. & Lange, K. (2009). Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 19, 1655–1664. Andrews, S. (2010). FastQC: a quality control tool for high throughput sequence data. https://www.bioinformatics.babraham.ac.uk/projects/fastqc/. Bazzi, G., Podofillini, S., Gatti, E., Gianfranceschi, L., Cecere, J.G., Spina, F., Saino, N., & Rubolini, D. (2017). Candidate genes have sex-specific effects on timing of spring migration and moult speed in a long-distance migratory bird. Curr. Zool. 63, 479–486. BirdLife International (2021). "Species factsheet: Coturnix coturnix" (2021, downloaded from http://www.birdlife.org on 03/02/2021). Bourgeois, Y.X.C., Bertrand, J.A., Delahaie, B., Cornuault, J., Duval, T., Milá, B., & Thébaud, C. (2016). Candidate gene analysis suggests untapped genetic complexity in melanin-based pigmentation in birds. J. Hered. 107, 327–335. Burtt, E. H. & Ichida, J. M. (2004). Gloger’s rule, feather-degrading bacteria, and color variation among Song Sparrows. Condor 106, 681–686. Camacho-Sanchez, M., Quintanilla, I., Hawkins, M.T.R., Tuh, F.Y.Y., Wells, K., Maldonado, J.E., & Leonard, J.A. (2018). Interglacial refugia on tropical mountains: Novel insights from the summit rat (Rattus baluensis), a Borneo mountain endemic. Divers. Distrib. 24, 1252–1266. Christmas, M. J. Wallberg, A., Bunikis, I., Olsson, A., Wallerman, O., & Webster, M.T. (2019). Chromosomal inversions associated with environmental adaptation in honeybees. Mol. Ecol. 28, 1358–1374. CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 117 | PA G E Cingolani, P., Platts, A., Wang, L., Coon, M., Nguyen, T., Wang, L., Land, S.J., Lu, X., & Ruden, D.M. (2012). A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff. Fly 6, 80–92. Danecek, P., Auton, A., Abecasis, G., Albers, C.A., Banks, E., DePristo, M.A., Handsaker, R.E., Lunter, G., Marth, G.T., Sherry, S.T., et al.; 1000 Genomes Project Analysis Group (2011). The variant call format and VCFtools. Bioinformatics 27, 2156–2158. de Queiroz, A. (2005). The resurrection of oceanic dispersal in historical biogeography. Trends Ecol. Evol. 20, 68–73. Delhey, K. (2019). A review of Gloger’s rule, an ecogeographical rule of colour: definitions, interpretations and evidence. Biol. Rev. 94, 1294–1316. El Sabry, M. I., Eshak, M. G., Stino, F. K. R. & Bondioli, K. R. (2017). Comparing growth, immune and pigmentation related gene expression in three lines of Japanese and wild European quail. Anim. Sci. Pap. Reports 35, 407–418. Elshire, R. J., Glaubitz, J.C., Sun, Q., Poland, J.A., Kawamoto, K., Buckler, E.S., & Mitchell, S.E. (2011). A Robust, Simple genotyping-by-sequencing (GBS) approach for high diversity species. PLoS One. 6, e19379. Faria, R. & Navarro, A. (2010). Chromosomal speciation revisited: Rearranging theory with pieces of evidence. Trends Ecol. Evol. 25, pp. 660–669. Fiedler, W. (2005). Ecomorphology of the external flight apparatus of blackcaps (Sylvia atricapilla) with different migration behavior. Ann. N. Y. Acad. Sci. 1046, 253–263. Fox, J. & Weisberg, S. An R Companion to Applied Regression. Sage Publications, Thousand Oaks, California, ed. 3 (2019). Garrison, E. & Marth, G. (2012). Haplotype-based variant detection from short-read sequencing. http://arxiv.org/abs/1207.3907. Glaubitz, J. C., Casstevens, T.M., Lu, F., Harriman, J., Elshire, R.J., Sun, Q., and Buckler, E.S. (2014). TASSEL-GBS: A high capacity genotyping by sequencing analysis pipeline. PLoS One. 9, e90346. Guyomarc’h, J. C., Combreau, O., Puigcerver, M., Fontoura, P.A. & Aebischer, N.J. (1998). Coturnix coturnix Quail. Birds of the Western Palearctic Update 2, 27-46. Holynski, R. (1965). The methods of analysis of wing-formula variability (English summary). Notatki Orn. 6, 21-25. Huang, Y. C., Dang, V. D., Chang, N. C. & Wang, J. (2018). Multiple large inversions and breakpoint rewiring of gene expression in the evolution of the fire ant social supergene. Proc. R. Soc. B Biol. Sci. 285. Jay, P., Whibley, A., Frézal, L., Rodríguez de Cara, M.Á., Nowell, R.W., Mallet, J., Dasmahapatra, K.K., & Joron, M. (2018). Supergene evolution triggered by the introgression of a chromosomal inversion. Curr. Biol. 28, 1839-1845.e3. Jenni, L. & Winkler, R. (1989). The feather-length of small passerines: A measurement for wing-length in live birds and museum skins. Bird Study 36, 1– 15. Jiménez-Blasco, I., Puigcerver, M., Gonçalves, D., García-Galea, E. & Rodríguez- Teijeiro, J. D. (2015). Morphological differentiation in individuals with different migratory strategies: the case of the common quail (Coturnix coturnix). Presented CHAPTER 2 118 | PA G E at the XXXII Congress International Union Game Biologist (IUGB), Puebla, Mexico, 24-28 August 2015. Katoh, K. & Standley, D. M. (2013). MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol. Biol. Evol. 30, 772– 780. Kirkpatrick, M. & Barton, N. (2006). Chromosome inversions, local adaptation and speciation, Genetics 173, 419–434. Kirubakaran, T. G., Grove, H., Kent, M.P., Sandve, S.R., Baranski, M., Nome, T., De Rosa, M.C., Righino, B., Johansen, T., & Ottera, H. (2016). Two adjacent inversions maintain genomic differentiation between migratory and stationary ecotypes of Atlantic cod. Mol. Ecol. 25, 2130–2143. Kopelman, N. M. Mayzel, J. Jakobsson, M., Rosenberg, N. A. & Mayrose, I. (2015). Clumpak: A program for identifying clustering modes and packaging population structure inferences across K. Mol. Ecol. Resour. 15, 1179–1191. Küpper, C., Stocks, M., Risse, J.E., Dos Remedios, N., Farrell, L.L., McRae, S.B., Morgan, T.C., Karlionova, N., Pinchuk, P., & Verkuil, Y.I. (2016). A supergene determines highly divergent male reproductive morphs in the ruff. Nat. Genet. 48, 79–83. Lamichhaney, S., Fan, G., Widemo, F., Gunnarsson, U., Thalmann, D.S., Hoeppner, M.P., Kerje, S., Gustafson, U., Shi, C., Zhang, H., et al. (2016). Structural genomic changes underlie alternative reproductive strategies in the ruff (Philomachus pugnax). Nat. Genet. 48, 84–88. Lenth, R. (2020). emmeans: estimated marginal means, aka least-squares means. https://cran.r-project.org/web/packages/emmeans/index.html. Li, H. & Durbin, R. (2009). Fast and accurate short read alignment with Burrows- Wheeler transform. Bioinformatics 25, 1754–1760. Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., Marth, G., Abecasis, G., & Durbin, R.; 1000 Genome Project Data Processing Subgroup (2009). The sequence alignment/map format and SAMtools. Bioinformatics 25, 2078–2079. Liu, X., Zhou, R., Peng, Y., Zhang, C., Li, L., Lu, C., & Li, X. (2019). Hair follicles transcriptome profiles in Bashang long-tailed chickens with different plumage colors. Genes Genomics 41, 1357–1367. Lockwood, R., Swaddle, J. P. & Rayner, J. M. (1998). Avian wingtip shape reconsidered: wingtip shape indices and morphological adaptations to migration. J. Avian Biol. 29, 273–292. Loftus, S.K., Larson, D.M., Baxter, L.L., Antonellis, A., Chen, Y., Wu, X., Jiang, Y., Bittner, M., Hammer, J.A. III, & Pavan, W.J. (2002). Mutation of melanosome protein RAB38 in chocolate mice. Proc. Natl. Acad. Sci. U. S. A. 99, 4471–4476. Madge, S. & McGowan, P. (2002). Pheasants, partridges and grouse (Princeton University). Martin, M. (2011). Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal 17, 10. CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 119 | PA G E Mastrangelo, S., Cendron, F., Sottile, G., Niero, G., Portolano, B., Biscarini, F., & Cassandro, M. (2020). Genome-wide analyses identifies known and new markers responsible of chicken plumage color. Animals 10, E493. McKenna, A., Hanna, M., Banks, E., Sivachenko, A., Cibulskis, K., Kernytsky, A., Garimella, K., Altshuler, D., Gabriel, S., Daly, M., & DePristo, M.A. (2010). The genome analysis toolkit: A MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 20, 1297–1303. Merritt, J.R., Grogan, K.E., Zinzow-Kramer, W.M., Sun, D., Ortlund, E.A., Yi, S.V. & Maney, D.L. (2020). A supergene-linked estrogen receptor drives alternative phenotypes in a polymorphic songbird. PNAS 117, 21673–21680. Mettler, R., Segelbacher, G. & Schaefer, H. M. (2015). Interactions between a candidate gene for migration (ADCYAP1), morphology and sex predict spring arrival in blackcap populations. PLoS One 10, e0144587. Mönkkönen, M. (1995). Do migrant birds have more pointed wings?: A comparative study. Evol. Ecol. 9, 520–528. Mur, P. (1994). Contribution à la gestion des populations paléarctiques de caille des blés dans la phase européenne de son cycle annuel. Recherches méthodologiques sur la cinétique démographique et appréciation des facteurs de fluctuations. PhD thesis (University of Rennes I). Nam, K. et al. Molecular evolution of genes in avian genomes. Genome Biol. 11, R68 (2010). Navarro, A., Betrán, E., Barbadilla, A., & Ruiz, A. (1997). Recombination and gene flux caused by gene conversion and crossing over in inversion heterokaryotypes. Genetics 146, 695–709. Newton, I. (2010). Bird migration (William Collins). Pacheco, S., Maldonado-Linares, A., Marcet-Ortega, M., Rojas, C., Martínez- Marchal, A., Fuentes-Lazaro, J., Lange, J., Jasin, M., Keeney, S., Fernaández- Capetillo, O., et al. (2018). ATR is required to complete meiotic recombination in mice. Nat. Commun. 9, 1–14. Pacheco, S., Marcet-Ortega, M., Lange, J., Jasin, M., Keeney, S., & Roig, I. (2015). The ATM Signaling cascade promotes recombination-dependent pachytene arrest in mouse spermatocytes. PLOS Genet. 11, e1005017. Patterson, N. Price, A. L. & Reich, D. (2006). Population structure and eigenanalysis. PLoS Genet. 2, 2074–2093. Price, A.L., Patterson, N.J., Plenge, R.M., Weinblatt, M.E., Shadick, N.A., & Reich, D. (2006). Principal components analysis corrects for stratification in genome-wide association studies. Nat. Genet. 38, 904–909. Protas, M. E. & Patel, N. H. (2008). Evolution of coloration patterns. Annu. Rev. Cell Dev. Biol. 24, 425–446. Puigcerver, M., Gallego, S., Rodríguez-Teijeiro, J. D., & Senar, J. C. (1992). Survival and mean life span of the quail Coturnix c. coturnix. Bird Study 39, 120–123. Puigcerver, M., Rodrigo-Rueda, F. J., Rodríguez-Teijeiro, J. D. & Gallego, S. (1997). On the second clutches in the common quail (Coturnix coturnix). Gibier Faune Sauvage 14(4), 617–622. CHAPTER 2 120 | PA G E Purcell, S., Neale, B., Todd-Brown, K., Thomas, L., Ferreira, M.A., Bender, D., Maller, J., Sklar, P., de Bakker, P.I., Daly, M.J., & Sham, P.C. (2007). PLINK: A tool set for whole-genome association and population-based linkage analyses. Am. J. Hum. Genet. 81, 559–575. Puritz, J.B., Hollenbeck, C.M. & Gold, J.R. (2014). dDocent: a RADseq, variantcalling pipeline designed for population genomics of non-model organisms. PeerJ 2:e43. R Development Core Team (2020). R: A language and environment for statistical computing (R Foundation for Statistical Computing). Rando, J. C., Alcover, J.A., Pieper, H., Olson, S.L., Hernández, C.N., and López- Jurado, L.F. (2020). Unforeseen diversity of quails (Galliformes: Phasianidae: Coturnix) in oceanic islands provided by the fossil record of Macaronesia. Zool. J. Linn. Soc. 188, 1296–1317. Rasband, W.S. (2018). ImageJ (US National Institutes of Health). https://imagej.nih.gov/ij/. Robinson, J.T., Thorvaldsdóttir, H., Winckler, W., Guttman, M., Lander, E.S., Getz, G., & Mesirov, J.P. (2011). Integrative genomics viewer. Nat. Biotechnol. 29, 24– 26. Rodríguez-Teijeiro, J. D., Barroso, A. Gallego, S. Puigcerver, M. & Vinyoles, D. (2006). Orientation-cage experiments with the European Quail during the breeding season and autumn migration. Can. J. Zool. 84, 887–894. Rodríguez-Teijeiro, J. D., Puigcerver, M. & Gallego, S. (1992). Mating strategy in the European quail (Coturnix c. coturnix) revealed by male population density and sex ratio in Catalonia (Spain). Gibier Faune Sauvage 9, 377-386. Rodríguez-Teijeiro, J. D., Sardà-Palomera, F., Nadal, J., Ferrer, X., Ponz, C., & Puigcerver, M. (2009). The effects of mowing and agricultural landscape management on population movements of the common quail. J. Biogeogr. 36, 1891–1898. Roig, I., Dowdle, J.A., Toth, A., de Rooij, D.G., Jasin, M., & Keeney, S. (2010). Mouse TRIP13/PCH2 is required for recombination and normal higher-order chromosome structure during meiosis. PLoS Genet. 6, e1001062. RStudio Team (2020). RStudio: integrated development for R (R Foundation for Statistical Computing). Rubenstein, D. R. & Hobson, K. A. (2004). From birds to butterflies: Animal movement patterns and stable isotopes. Trends Ecol. Evol. 19, pp. 256–263. Saint-Jalme, M. & Guyomarc’h, J. C. (1995). Plumage development and moult in the European Quail Coturnix c. coturnix: criteria for age determination. Ibis 137 (4), 570-581. Saint-Jame, M. (1990). La reproduction chez la caille des blés (Coturnix c. coturnix); études expérimentales des cycles saisonniers et de la variabilité interindividuelle. PhD thesis (Université de Rennes I). Schwander, T., Libbrecht, R. & Keller, L. Supergenes and complex phenotypes (2014). Curr. Biol. 24, R288–R294. CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 127 | PA G E Figure S4. Linkage disequilibrium within chromosome 1. Related to Figure 3. Weir and Cockerham's FST along chromosome 1 in 1Kbp windows comparing homokaryotype AA (N= 6) and BB (N= 10) quails. A 115 Mbp region of chromosome 1 shows very high levels of genetic differentiation due to an inversion. The orange line represents FST values in 1 Mbp windows. Below, linkage disequilibrium (LD) heatmap. Blue regions indicate higher LD, corresponding with the inverted region. CHAPTER 2 | APPENDIX 128 | PA G E Figure S5. Variability of throat pigmentation. Related to STAR Methods and Data S1A. Throat pigmentation was categorized into six levels, from pale (level 1) to completely pigmented (level 6), by visual examination of throat pictures and comparing them to a predefined scale. Figure S6. Morphology of AA common quail males in relation to their sampling locality. Related to Figure 4 and Data S1A. Northern localities (N): Netherlands, Italy, NW Spain, NE Spain and CE Spain. Southern localities (S): S Spain, Morocco, Canary Islands and Cape Verde. (A) Body weight. (B) Proportion of individuals with different degrees of throat pigmentation. (C) Width of the pectoral lipid band. (D) Variation in the Holynski index, which measures wing pointedness and is related to flight efficiency. CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 129 | PA G E SUPPLEMENTARY TABLES Table S1. Immunofluorescence essays on spermatocytes in meiotic prophase. Related to Fig 2A and STAR Methods. Average distance between two fluorochromes (K19 and J23) along chromosome 1 of quail spermatocytes. Karyotype composition indicates individuals identified as homokaryotypes without inversion (AA), heterokaryotypes (AB) or homokaryotypes with the inversion (BB) as represented in Fig. 2A. Kar., karyotype group. Sample ID Sampling locality 1 Sampling locality 2 Kar. Average distance (K19-J23), μm 17061401 NE Spain St. Boi de Lluçanès AA 24.79 17061402 NE Spain St. Boi de Lluçanès AA 20.11 17061600 NE Spain Alp AA 20.03 17061609 NE Spain Alp AA 22.55 17061614 NE Spain Alp AA 31.54 17061615 NE Spain Alp AA 23.67 17062601 S Portugal Aljezur BB 11.43 17062602 S Portugal Tavira BB 13.50 17062603 S Portugal Tavira BB 13.53 17062609 S Portugal Aljezur BB 10.53 17062613 S Portugal Aljezur BB 10.85 17062604 S Portugal Tavira AB 28.37 17062606 S Portugal Aljezur AB 20.71 17062607 S Portugal Aljezur AB 31.02 17062610 S Portugal Aljezur AB 20.95 CHAPTER 2 | APPENDIX 130 | PA G E Table S2. Whole genome resequencing of common quails. Homokaryotype AA and BB quails, as identified by genotyping-by-sequencing, resequenced. Kar., karyotype group. Individual ID in Fig. 3 Sampling locality 1 Sampling locality 2 Kar. Number of mapped reads Coverage 09051003 Italy Italy Monte Brisighella AA 82,693,872 9.6 08030502 Morocco 1 Morocco El Gara AA 70,519,520 8.5 08030503 Morocco 2 Morocco El Gara BB 67,394,401 8.2 080410XX Morocco 3 Morocco Beni Mellal BB 90,962,857 10.6 14052112 S Portugal 1 S Portugal Tavira BB 85,400,857 10.1 14052312 S Portugal 2 S Portugal Odiaxere BB 92,760,611 10.7 14042507 NE Spain NE Spain Figuerola AA 78,214,047 9.6 12052702 NW Spain NW Spain Valdesogo AA 97,266,382 11.6 15051702 S Spain S Spain Coronil AA 114,679,132 13.5 14041210 El Hierro 1 Canary Islands El Hierro BB 94,729,231 11.4 14041313 El Hierro 2 Canary Islands El Hierro AA 79,754,285 9.4 15041807 Gran Canaria Canary Islands Gran Canaria BB 75,682,609 8.7 14040907 Tenerife 1 Canary Islands Tenerife BB 77,717,165 9.3 14040910 Tenerife 2 Canary Islands Tenerife BB 91,084,220 10.6 15061003 Madeira 1 Madeira Madeira BB 76,739,862 9.3 15061103 Madeira 2 Madeira Madeira BB 87,358,186 10.0 CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 131 | PA G E CHAPTER 3 EVOLUTIONARY HISTORY OF AN ISLAND ENDEMIC, THE AZOREAN COMMON QUAIL ADAPTED FROM: Sara Ravagni, Ines Sanchez-Donoso, Irene Jiménez-Blasco, Pedro Andrade, Manel Puigcerver, Ana Chorão Guedes, Raquel Godinho, David Gonçalves, Manuel Leitão, Jennifer A. Leonard, José Domingo Rodríguez-Teijeiro & Carles Vilà (2023). Evolutionary history of an island endemic, the Azorean common quail. Molecular Ecology. https://doi.org/10.1111/mec.16997 CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 135 | PA G E ABSTRACT Oceanic islands are characterised by conditions that favour diversification into endemic lineages that can be very different from their mainland counterparts. This can be the result of fast phenotypic divergence due to drift or the result of slower adaptation to local conditions. This uniqueness can obscure their evolutionary history. Here we used morphological, stable isotope, genetic and genomic data to characterize common quails (Coturnix coturnix) in the Azores archipelago and assess the divergence from neighbouring common quail populations. Historical documents suggested that these quails could have a recent origin associated with the arrival of humans in the last centuries. Our results show that Azorean quails constitute a welldifferentiated lineage with small size and dark throat pigmentation that has lost the migratory ability and that diverged from mainland quail lineages more than 0.8 mya, contrary to the notion of a recent human-mediated arrival. Even though some Azorean quails carry an inversion that affects 115 Mbp of chromosome 1 and that has been associated with the loss of the migratory behaviour in other common quail populations, half of the analysed individuals do not have that inversion and still do not migrate. The long coexistence and evolution in isolation in the Azores of two chromosomal variants (with and without the inversion) is best explained by balancing selection. Thus, a unique and long evolutionary history led to the island endemic that we know today, C. c. conturbans. KEYWORDS Biogeography, chromosomal rearrangements, differentiation, genotyping by sequencing, Macaronesian islands, subspecies CHAPTER 3 136 | PA G E INTRODUCTION Oceanic islands are unique systems to study evolution due to their intrinsic characteristics (Losos & Ricklefs, 2009; Warren et al., 2015; Whittaker et al., 2017). They represent independent replicates of evolutionary processes with often known geological ages and hard boundaries that limit immigration from other islands or the mainland. Fragmentation results in differential selective forces and accelerated genetic drift, due to small population size, facilitates increased divergence from the ancestral forms (Leroy et al., 2021). At the same time, taxa can evolve to use empty ecological niches, leading to community composition well differentiated from the closest mainland. The isolation and simpler biotic community structure in islands also offer unique opportunities for adaptive radiations as in Anolis lizards or Eleutherodactylus frogs in Caribbean islands (Dugo-Cota et al., 2019; Losos et al., 1998; Losos & Schneider, 2009). Moreover, there are hints of deterministic patterns and some repeated trends are found in separate islands ("the island syndrome", see Baeckens & Van Damme, 2020). This is exemplified by the tendency towards dwarfism and gigantism in mammals (Athanassiou et al., 2019; Benítez-López et al., 2021; Bromham & Cardillo, 2007) or flightlessness in birds and insects (Losos & Ricklefs, 2009; Whittaker et al., 2017). Under these conditions, divergence can proceed fast (Knope et al., 2012; Lamichhaney et al., 2018) and the degree of phenotypic differentiation between islands and continental forms may not be an indication of the time of divergence. This has sometimes resulted in problems assessing if divergent insular forms are the result of an ancient evolutionarily independent lineage or represent a plastic response to the insular conditions (for example, Aleixandre et al., 2013; Yahnke et al., 1996). One example of island diversification is that of the common quail (Coturnix coturnix) in the Macaronesian archipelagos (including Azores, Madeira, Canary Islands and Cape Verde), in the eastern Atlantic Ocean. The common quail is a small migratory bird with a wide distribution range in the Palearctic region, from the Atlantic islands to Mongolia, and extending to southern Africa and southern India (Figure 1, Guyomarc’h et al., 1998). Six geographic subspecies have been described based on plumage patterns and wing length (del Hoyo et al., 1994; Guyomarc'h & Perennou, 2009). However, half of the subspecies correspond to the Macaronesian archipelagos, highlighting the relevance of oceanic islands in differentiation: Coturnix coturnix coturnix, found in Europe, North Africa and Asia; C. c. erlangeri, in East Africa (Ethiopia); C.c. africana, in Southern Africa; C.c. conturbans, in the Azores; C.c. confisa, in Canary Islands and Madeira; C. c. inopinata, in Cape Verde. CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 143 | PA G E K. We combined results from the different runs using the CLUMPAK server (Kopelman et al., 2015) with default settings. We calculated individual heterozygosity using VCFTOOLS for the entire SNP dataset, excluding SNPs from the inverted region to remove the main source of differentiation between AA and BB samples, and only with SNPs from the inverted region. To estimate the effective population size of the Azorean quails, we used the software GONE (Santiago et al., 2020). It estimates population Ne values in the last 150-200 generations based on linkage-disequilibrium patterns. We generated input files using PLINK selecting the 10 Azorean quails and excluding the inverted region in chromosome 1. We ran 50 independent estimates of the GONE analysis using default parameters except for the recombination rate. The recombination rate was calculated based on the Japanese quail genome as (genetic map length in cM)/(genome size in Mb), in our case 3057/927.66 = 3.295 cM/Mb (Recoquillay et al., 2015). We assumed the recombination rate to be similar in the common quail. As the estimates for the most recent 5 to 10 generations are less reliable because of the fewer recombination events, we estimated the current effective population size as the average for 10 to 20 generations ago. To estimate the absolute divergence (dXY) between Azorean and AA quails, we ran the STACKS Version 2.62 pipeline (Catchen et al., 2011; Catchen et al., 2013). We mapped reads against the Japanese quail genome using BWA-MEM v. (Li & Durbin, 2009) and we used the ref_map.pl script implemented in Stacks to call variants and calculate population genetic statistics based on a population map file separating samples according to their genotypes or sampling locality (Azores, AA, AB and BB). We first selected only the inverted region in chromosome 1 and then the rest of the genome, excluded heterozygous samples from the analysis and included the flag -X “populations:--min-maf 0.03 -p 2 -R 0.75 --fstats –smooth” to enable haplotype-based statistics. We estimated 95% confidence intervals for dXY using 1000 bootstrap replicates with the package boot (Canty & Ripley, 2022) in R. We calibrated the divergence time between Azores and AA samples based on the results obtained in Sanchez-Donoso et al. (2022), where we estimated that inverted BB and not inverted AA sequences diverged at least 0.9 mya (mean divergence time ranging from 0.9 to 3.2 mya, depending on the mutation rate used for the calculations) based on published mutation rates for several bird groups (Nam et al., 2010; Smeds et al., 2016). Population structure and PCA analyses revealed that some individuals from the Azores carried the inversion in chromosome 1 (haplotype B, see below). This inverted sequence was found to be divergent between Madeira, the Canary Islands and the CHAPTER 3 144 | PA G E mainland (Sanchez-Donoso et al., 2022). The study of the homologous inverted sequence (haplotype BAz) in the Azores could reveal if it was likely to derive from the (natural or artificial) recent arrival of quails carrying the inversion from one of those areas. However, we did not find homozygous individuals for the inverted sequence among the Azorean quails (see Results) and phasing the long inversion from so few individuals turned out to be a complex problem. To shed light on the origin of the inverted chromosome BAz in the Azores, we used SNP data from homozygous quails without the inversion from the Azores (we will call them AAzAAz, since the sequence is clearly divergent from AA reported by Sanchez-Donoso et al., 2022) and BB quails (homozygous with the inversion) from Canary Islands, Madeira and the Iberian Peninsula. We created synthetic AAzB genotypes from each location that could be compared with the heterozygous AAzBAz individuals from the Azores. If the Azorean BAz derived from a recent introduction (or natural arrival) from one of the other populations, we would expect AAzBAz genotypes to be similar to the synthetic AAzB generated for that locality. If they were different from all, the results could suggest that the BAz haplotype in the Azores had been evolving separately from the others. We first converted the vcf file generated by TASSEL into the genepop format using PGDSpider v.2.2.2.5 (Lischer & Excoffier, 2012). We then generated 20 synthetic hybrids for each combination of Azores AAzAAz (n = 4 individuals) with BB from Portugal (n = 5), BB from the Canary Islands (n = 4) and BB from Madeira (n = 2) with the software HybridLab version 1.1 (Nielsen et al., 2006). We performed a PCA including the simulated hybrids and the heterozygous Azorean quails with the function dupi.pca implemented in the R package adegenet (Jombart & Bateman, 2008). RESULTS MORPHOLOGICAL ANALYSES A PCA of morphological and throat pigmentation data for 437 quails including 49 individuals from the Azores, 84 from the Canary Islands, 24 from Madeira and 280 from the continental mainland showed that quails from the Azores were differentiated from those from the other localities along the first component (Figure 2a), with smaller weight, tarsus, and wing length. Azorean quails had dark throat colouration along the second component. On the other hand, quails from the other locations showed extensive overlap and were not clearly separated. CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 145 | PA G E FIGURE 2. (a) Principal component analysis of 437 common quail males analysed for size (weight, wing length, tarsus length) and throat colouration. Colours represent localities, and larger dots represent the centroid for each population. The centroid for Azores is very separated from the centroids for the other populations. (b) Proportion of individuals in each locality belonging to the three clusters resulting from a GMM applied to morphological data from 437 common quail males. (c) Morphological characterization of the three clusters (mean and 95% confidence intervals). CHAPTER 3 146 | PA G E We used these data to investigate if distinct clusters of phenotypes could be identified that grouped the quails. The Hopkins index had a value of 0.716, which is higher than the 95% quantile of the Beta distribution for 10% of the sample (n = 43; β(43,43) = 0.588), confirming that data are clusterable (p < .001). The Bayesian Information Criterion suggested that data were best explained by the presence of three clusters. Despite the overlap between the clusters in the PCA (Figure S1), cluster 2 showed geographic segregation from the other two and was greatly restricted to the Azores (Figure 2b). Although there was variation in the phenotype of Azorean quails, over 60% of the individuals were included in this cluster while less than 5% of the individuals from the other Macaronesian archipelagos or the mainland were assigned to this group. This cluster included small birds (small weight, tarsus length and wing length) with dark throats (Figure 2c). Although the overall body size of the birds in clusters 1 and 3 was similar, birds from Cluster 1 were slightly heavier, with longer tarsus, shorter wings, and darker throat colouration (Figure 2c). This is the phenotype previously described for quails carrying the inversion in chromosome 1 (haplotype B), which is associated with reduced migration and is more frequent in the Canary Islands and Madeira. It is also present but at a lower frequency on the mainland, concentrating in the southwest of the Iberian Peninsula and northern Morocco (Sanchez-Donoso et al., 2022). Cluster 3 birds had lighter throats, smaller bodies and longer wings, and correspond to individuals with the previously described and most common mainland chromosome A. ISOTOPE ANALYSES The ratios of δ13C and δ2H isotopes in feathers are determined by the environmental conditions and diet in the place where the feather was grown. Therefore, differences in the isotopic composition between breeding and wintering feathers are informative regarding the use of different breeding and wintering grounds. In Sanchez-Donoso et al. (2022), quails carrying one or two copies of the inversion (homozygotes BB and heterozygotes AB) showed very little differences between breeding and wintering feathers compared to AA homozygotes, suggesting very limited migratory movements or lack of them. The analysis of the isotopic composition of 10 Azorean quails showed that seasonal differences were significantly smaller than for AA quails, but similar to BB quails from the mainland and other Macaronesian islands with reduced migration (δ13C model: F(2,53 )= 6.377, p = .003; Azores vs. AA, p = .020; Azores vs. BB, p = .573; δ2H model: F(2,57 )= 4.841, p = .011; Azores vs. AA, p = .007; Azores vs. BB, p = .379; Figure 3). This suggests that quails in the Azores do not migrate to wintering grounds outside the archipelago. CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 147 | PA G E FIGURE 3. Proportion of deuterium δ2H (a) and carbon δ13C isotopes (b) in two primary feathers, one grown during the previous breeding season and the other during the wintering season, in 1-year-old Azorean male quails as well as male quails with BB and AA genotypes from the mainland, Canary Islands and Madeira from Sanchez-Donoso et al. (2022). GENETIC SUBSTRUCTURE WITHIN THE AZORES No apparent population structure was found among the islands with a principal coordinate analysis of microsatellite data (PCoA, Figure S2). Nevertheless, an AMOVA analysis suggested some degree of differentiation, which was significant both when comparing all islands separately (p = .001) or Central and Eastern groups (p < .001). The differences were driven by a slight divergence of São Miguel quails compared to the Central group quails, as suggested by the positions along the first axis of the PCoA. However, only 2% of the total genetic variance could be associated with population differentiation (Table S2). CHAPTER 3 148 | PA G E Based on the log probability for each value of K (number of clusters) in STRUCTURE, the value that best explained the data was K = 1, which was further confirmed by the admixture observed in all individuals when we considered two ancestry clusters (K = 2, see Figure S3). A slight degree of population substructure within the archipelago was only noted when adding sampling locality information. In this case, the most likely value of K was K = 4, but the groups that were generated did not have a clear biological meaning other than a slightly reduced gene flow between islands. Given this very small differentiation, all Azorean quails were considered as a single group in genomic analyses. GENOMIC ANALYSES Genotyping-by-sequencing analyses revealed 56,039 SNPs that were used to compare quails from different origins. A PCA including Japanese quails and game farm hybrids (used for restocking in the Iberian Peninsula; Sanchez-Donoso et al., 2012) showed that common quail populations -including the Azoresdid not have evidence of introgression of Japanese quail alleles (Figure S4). Consequently, introgression of Japanese quail alleles was not likely to contribute to population differentiation and Japanese and admixed quails were excluded from subsequent analyses. The exploratory PCA of all common quails showed that Azorean quails were genetically divergent, mainly along the second axis (Figure 4a). The quails sampled on the mainland and other Macaronesian archipelagos grouped according to the presence of the chromosomal inversion in chromosome 1 as described in Sanchez-Donoso et al. (2022), while the Azorean samples were differentiated, also forming clusters. The study of Sanchez-Donoso et al. (2022) showed that the inversion in chromosome 1 played a key role in the differentiation of common quail populations. To investigate if quails from the Azores carried the inversion, we did a PCA including only SNPs located within the inversion. The PCA showed that five Azorean quails clustered close to AB quails from other localities, closer to the quails from other Macaronesian archipelagos, confirming the presence of heterozygous individuals for the inversion FIGURE 4. Principal-component analysis (PCA) for genomic data of 90 male quails separated by location and genotype. Azorean quails (in green) appear differentiated from the other common quails, which group according to the presence or absence of a very large inversion in chromosome 1 (A, chromosome 1 without inversion; B, chromosome 1 with the inversion). (a) Using SNPs from the whole genome (56,039 SNPs); (b) analysing only the inverted region (7,578 SNPs) and (c) excluding the inverted region in chromosome 1 (48,461 SNPs). Data from locations other than the Azores come from Sanchez-Donoso et al. (2022). CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 149 | PA G E CHAPTER 3 150 | PA G E (Figure 4b). The other five individuals from Azores clustered close to AA quails. In both cases, Azorean quails remained somewhat separated from all others. The presence of the inversion in some birds from the Azores was supported by analyses with the software ADMIXTURE. While analyses with the complete SNP dataset showed similarity between Azorean and AA quails (Figure S5a), analysis using SNPs from the region of the inversion and K = 2 showed that five Azorean quails clustered together with the AA quails while the remaining five individuals were a half-and-half mix of the two clusters, the same as AB quails (Figure S5b). For K = 3 the Azorean quails without the inversion separated from the other AA and AB, indicating sequence divergence. As this implies that the non-inverted haplotype in chromosome 1 in the Azores was genetically differentiated from the A chromosome in other localities, we refer to it as AAz. Since the inversion affects a large portion of the genome (12%), it can have a strong effect on the clustering analyses when assessing the differences over the entire genome. Thus, we repeated the previous analyses excluding the inverted region in chromosome 1. The PCA showed that Azorean quails grouped separately from all other quails (Figure 4c), and the same pattern was observed in the ADMIXTURE analysis at K = 2 (Figure S5c). Thus, while all other quails formed a single group when the region of the inversion was excluded, suggesting extensive admixture between populations (see also Sanchez-Donoso et al. 2022), Azorean quails remained very differentiated. Among these, two individuals appeared also differentiated, mainly due to divergences along chromosome 1 but outside the inversion. Even though ADMIXTURE analyses showed that the quails from the Azores were well differentiated from other quails, the PCA (Figure 4a,c) suggested that they contain high intrapopulation variability (even after excluding the two particularly divergent individuals), which was surprising considering the limited geographic distribution and potentially small effective population size. The estimated effective population size of the Azorean quails was large, around 4600 (50 replicates, mean= 4602, SD= 228). The individual heterozygosity calculated for the whole genome and excluding the inverted region was similar for AA, BB and AB quails, confirming the lack of a strong effect due to reduced population size (Figure S6). This genetic diversity was associated to phenotypic diversity as well. Although Azorean quails are morphologically differentiated from other populations, a PCA of the morphological characteristics of the 49 quails sampled in the Azores showed intrapopulation diversity (Figure 5). Among these, the quails whose genome was analysed and that were carrying the inversion presented a darker throat pigmentation compared to those without the inversion. However, they did not differ in body size. CHROMOSOMAL INVERSION AND COMMON QUAIL DIVERSIFICATION 151 | PA G E FIGURE 5. PCA of 49 common quail males from the Azores analysed for size (weight, wing length, tarsus length) and throat colouration. The 10 genotyped individuals are coloured according to their haplotype for the inversion in chromosome 1. Despite the small sample size, the inversion seems to have an impact on the pigmentation, similar to the effect shown in other archipelagos and the mainland (Sanchez-Donoso et al. 2022). However, when comparing the differences in isotopic signatures between breeding and wintering feathers, quails not carrying the inversion (AAzAAz) did not show larger differences, suggesting lack of migratory behaviour. Consequently, this behaviour was not associated with the inverted haplotype BAz (Figure S7). We used a Manhattan plot to evaluate the differentiation between Azorean and AA quails along the genome. The mean weighted FST was 0.044, and highly differentiated SNPs were distributed across the whole genome rather than forming well-defined clusters (Figure S8). This confirmed that Azorean quails were very divergent from the other quails across the inversion as well as the entire genome. This sequence difference suggests ancient divergence between Azorean AAz and other A haplotypes. The estimated time of divergence between A and B is at least 0.9 mya, based on dXY = 0.010 (Sanchez-Donoso et al. 2022). Applying this calibration to the divergence between Azorean AAz and the widespread A (dXY = 0.009, 95% CI: 0.0089-0.0092), it results in a minimum divergence time of 0.814 mya. Haplotypes B CHAPTER 3 152 | PA G E and AAz are even more divergent (dXY = 0.011, 95% CI: 0.0109-0.0114) and so likely diverged even earlier, at least 1 mya. As we are using the fastest of several mutation rates based on other species, the divergence times are rough estimates and are likely to be a minimum estimate of the divergence time. Also, these estimates assume no gene flow. If there had been occasional gene flow, the real divergence times could be even older. Overall, our results suggest that the divergence between A and AAz is more recent than between A and B, and that these divergences are very old events. The divergences estimated for the rest of the genome are similarly large, showing that the estimates are not just a consequence of the inversion but reflect the evolutionary history of the population: dXY = 0.009 when comparing Azorean quails not carrying the inversion with both AA and BB (95% CI: 0.0086-0.0088 and 0.0094- 0.0095, respectively). The presence of haplotype BAz in the Azores could be explained by natural dispersal or release of quails originating from other areas of the distribution, a common practice associated with restocking for hunting (Caravaca et al., 2022). The individuals carrying haplotype B have low dispersal capacity (Sanchez-Donoso et al., 2022) and this has resulted in different variants of the haplotype in the continental landmasses and each one of the Atlantic archipelagos studied (Canary Islands and Madeira). Therefore, if the inversion had arrived to the Azores as a result of an ancient natural dispersal and evolved independently, we would expect it to be genetically differentiated from the other localities. Alternatively, if it was introduced by humans, we expect to find the same lineage present in other localities, probably mainland Portugal. Since no BAzBAz homozygous individuals were found in our sample from the Azores, we tested this hypothesis by generating synthetic heterozygous AAzB genotypes between Azorean AAzAAz quails and BB homozygotes from Portugal, Madeira and the Canary Islands. These synthetic genotypes were then compared to Azorean heterozygous quails (AAzBAz). A PCA showed that they formed very differentiated groups (Figure 6). Since all of them included haplotype AAz from the Azores, the result implies that the B haplotype in the Azores is differentiated from those in the Canary Islands, Madeira, and the mainland. Consequently, this inverted sequence BAz was not likely to derive from recent restocking with quails from these other populations. Instead, the BAz haplotype in the Azores may be the result of evolution in isolation from other populations. Although we were not able to assess the divergence time between the inverted sequence BAz in the Azores and the other localities directly, we estimated the divergence of the inversion between homozygous BB from the Canary Islands, Madeira and the mainland as dXY. Using the mutation rate from Sanchez-Donoso et al. (2022), we inferred a minimum divergence time between the B haplotype in the Canary Islands and Madeira of 104,500 years. The