Genomics of divergence, gene flow and selection in the genus Lynx

Genomics of divergence, gene flow and selection in the genus Lynx Enrico Bazzicalupo Supervisor: José Antonio Godoy Tutor: Pedro Abellán ABSTRACT............................................................................................................................ 5 INTRODUCTION....................................................................................................................7 Evolutionary Forces Driving Genetic Differentiation....................................................8 Effective Population Size And Population Connectivity............................................... 8 Divergent Selection And Local Adaptation................................................................. 10 Gene Flow And Introgression...................................................................................... 11 Importance Of Conservation Genomics.......................................................................12 The Eurasian And Iberian Lynx................................................................................... 13 Thesis Summary And Goals.........................................................................................15 REFERENCES...................................................................................................................16 CHAPTER 1: History, demography and genetic status of Balkan and Caucasian Lynx lynx (Linnaeus, 1758) populations revealed by genome-wide variation.......................................................33 ABSTRACT.......................................................................................................................34 INTRODUCTION..............................................................................................................34 METHODS.........................................................................................................................37 Sampling and DNA extraction.....................................................................................37 Sequence data generation and processing....................................................................38 Mitogenomic phylogeography..................................................................................... 40 Nuclear demographic and divergence reconstruction using PSMC.............................40 Nuclear divergence and admixture reconstruction using TreeMix.............................. 41 Nuclear genomic structure........................................................................................... 41 Recent demographic history.........................................................................................42 Nuclear and mitogenomic diversity............................................................................. 42 Runs of Homozygosity (ROH) and inbreeding............................................................43 RESULTS...........................................................................................................................43 Mitogenomic phylogeography..................................................................................... 43 Nuclear divergence and admixture reconstruction using TreeMix.............................. 44 Nuclear demographic and divergence reconstruction using PSMC.............................46 Nuclear genomic structure........................................................................................... 48 Recent demographic history.........................................................................................50 Nuclear and mitogenomic diversity............................................................................. 51 Runs of Homozygosity (ROH) and inbreeding............................................................51 DISCUSSION.................................................................................................................... 53 Taxonomic implications............................................................................................... 54 Conservation implications............................................................................................55 REFERENCES...................................................................................................................56 APPENDIX........................................................................................................................62 CHAPTER 2: Genome-environment association analyses reveal geographically restricted adaptive divergence across the range of the widespread Eurasian carnivore Lynx lynx (Linnaeus, 1758).........................................................................................................................................83 2 ABSTRACT.......................................................................................................................84 INTRODUCTION..............................................................................................................84 METHODS.........................................................................................................................87 Sampling, DNA extraction, sequencing and read alignment....................................... 87 Variant calling and filtering..........................................................................................88 Environmental predictors.............................................................................................89 Identifying candidate regions of selection....................................................................90 Functional enrichment of candidate genes...................................................................92 Adaptive population structure......................................................................................92 Analysing spatial gradients of adaptive genetic differentiation...................................92 RESULTS...........................................................................................................................93 Variable selection and variance partitioning................................................................ 93 Genomic scans of selection..........................................................................................94 Functional enrichment of candidate genes – panther overrepresentation test..............94 Adaptive population structure......................................................................................95 Spatial gradients of adaptive genetic differentiation – GDM.......................................97 DISCUSSION.................................................................................................................... 98 Evidence of local adaptation........................................................................................ 99 Drivers of local adaptation and genes involved......................................................... 101 Implications for conservation and taxonomy.............................................................102 REFERENCES.................................................................................................................104 APPENDIX...................................................................................................................... 111 CHAPTER 3: Convolutional neural networks help unravel the history of introgression between two lynx species in South-Western Europe............................................................................... 123 ABSTRACT.....................................................................................................................124 INTRODUCTION............................................................................................................124 METHODS.......................................................................................................................128 Sampling, DNA extraction and sequencing...............................................................128 Sequencing reads alignment, variant calling and initial filtering...............................129 Demographic inference.............................................................................................. 130 Phasing variants..........................................................................................................131 Introgression scans.....................................................................................................131 RESULTS.........................................................................................................................134 Demographic inference.............................................................................................. 134 Introgression scans.....................................................................................................138 DISCUSSION.................................................................................................................. 141 REFERENCES.................................................................................................................145 APPENDIX......................................................................................................................152 DISCUSSION....................................................................................................................... 156 Evolutionary Forces Behind Lynx Divergence.......................................................... 157 3 Guiding Lynx Conservation....................................................................................... 158 Expanding Eurasian Lynx Sampling..........................................................................159 Additional Insight From Structural Variation.............................................................160 REFERENCES.................................................................................................................160 CONCLUSIONS...................................................................................................................166 4 ABSTRACT La diferenciación genética entre especies y poblaciones resulta de la interacción de fuerzas evolutivas como la deriva genética y la selección divergente, que promueven la divergencia, y el flujo génico, que la contrarresta. En esta tesis, analizamos secuencias de genomas completos de varios individuos de lince Eurasiático e Ibérico para entender cómo estas fuerzas modelan la historia evolutiva de estas poblaciones silvestres, examinando las implicaciones de nuestros resultados para su taxonomía y conservación. El primer capítulo presenta el primer análisis de genoma completo de dos poblaciones de lince Eurasiático: el lince Balcánico, en peligro crítico de extinción, y el lince Caucásico, previamente no muestreado, que reside en un posible refugio glacial. Se evalúan tanto las divergencias mitocondriales como autosómicas para entender sus relaciones con otras poblaciones. Evaluamos el estado genético de estas poblaciones en términos de tamaño efectivo de la población, diversidad genética y consanguinidad, y sacamos conclusiones sobre su estatus taxonómico y de conservación. El segundo capítulo investiga las adaptaciones locales a las condiciones ambientales en el rango de distribución del lince Eurasiático. Al escanear los genomas de múltiples individuos de diferentes poblaciones, identificamos ventanas genómicas que responden a las presiones selectivas impuestas por el ambiente, y destacamos qué condiciones ambientales específicas están impulsando la adaptación, junto con los genes asociados. Se comparan las divergencias genéticas adaptativas y neutrales, y se reconstruye la composición genética adaptativa en respuesta a variables ambientales clave a lo largo del rango de la especie. Estos hallazgos se discuten en el contexto de la taxonomía y biología de la conservación del lince Eurasiático. El tercer capítulo explora el flujo génico interespecífico mediante el análisis de patrones de introgresión entre el lince Ibérico y el Eurasiático. Usando un enfoque novedoso basado en redes neuronales convolucionales profundas, identificamos ventanas genómicas que muestran señales de introgresión en tres poblaciones distintas de linces. Basándonos en los patrones de introgresión detectados, discutimos la historia demográfica de las poblaciones remanentes junto con su relación con la población localmente extinta que habitaba la zona híbrida entre las dos especies. Además, se identifican las áreas del genoma más propensas a exhibir introgresión y se analizan las consecuencias de la introgresión sobre la diversidad genética de las poblaciones receptoras. En conjunto, estos capítulos ofrecen un análisis completo de las fuerzas genéticas que modelan las poblaciones de lince Eurasiático e Ibérico, proporcionando valiosas perspectivas sobre sus historias evolutivas, procesos adaptativos y necesidades de conservación. Genetic differentiation among species and populations results from the interplay of evolutionary forces such as genetic drift and divergent selection, which promote divergence, and gene flow, which counteracts it. In this thesis, we analyze whole genome sequences of various Eurasian and Iberian lynx individuals to understand how these forces shape the evolutionary history of these wild populations, examining the implications of our results for their taxonomy and conservation. The first chapter presents the first whole genome analysis of two Eurasian lynx populations: the critically endangered Balkan lynx and the previously unsampled Caucasian lynx, which resides in a potential glacial refugium. Both mitochondrial 5 and autosomal divergences are assessed to understand their relationships with other populations. We evaluate the genetic status of these populations in terms of effective population size, genetic diversity, and inbreeding, and draw conclusions regarding their taxonomic and conservation status. The second chapter investigates local adaptations to environmental conditions across the Eurasian lynx distributional range. By scanning the genomes of multiple individuals across different populations we identify genomic windows responding to selective pressures imposed by local environments, and highlight what specific environmental conditions are driving adaptation, together with the associated genes. Adaptive and neutral genetic divergences are compared, and the adaptive genetic composition in response to key environmental variables is reconstructed across the species’ range. These findings are discussed in the context of the Eurasian lynx taxonomy and conservation biology. The third chapter explores interspecific gene flow by analyzing introgression patterns between the Iberian and Eurasian lynx. Using a novel approach based on deep convolutional neural networks, we identify genomic windows showing signals of introgression in three distinct lynx populations. Based on the introgression patterns detected, we discuss the demographic history of remnant populations together with their relationship with the locally extinct population which inhabited the hybrid zone between the two species. Additionally, the areas of the genome most likely to exhibit introgression are identified, and the consequences of introgression on the genetic diversity of recipient populations are analyzed. Together, these chapters provide a comprehensive analysis of the genetic forces shaping the Eurasian and Iberian lynx populations, offering valuable insights into their evolutionary history, adaptive processes, and conservation needs. 6 INTRODUCTION 7 Evolutionary Forces Driving Genetic Differentiation Genetic differentiation can accumulate among lineages through changes in allele frequencies caused by the joint action of four main agents of evolutionary change: mutation, genetic drift, selection and migration. Mutation, by creating new alleles, is the ultimate source of genetic diversity. However, it typically occurs at a much lower rate compared to other evolutionary forces, which makes mutation alone cause only very gradual changes in allele frequencies over many generations (Kimura, 1964). Genetic drift is the product of random sampling of reproducing individuals in finite populations, and causes the population to experience allele frequency changes across generations that are essentially random (Buri, 1956). Because of the randomness of the process, the strength of genetic drift, which is the relative change in frequency that an allele can undergo in one generation, is dictated by the effective number of individuals in the population (Charlesworth, 2009). Neutral alleles in relatively small populations will have relatively large changes in frequencies across generations, while the opposite is true in relatively large populations. Selection interferes with the random process of drift, dictating the direction and strength of the changes in frequencies across generations. Alleles can be favoured or disfavoured by selective forces, increasing and decreasing their frequency respectively, while the selection coefficient of a particular allele will affect how large the change is from one generation to the next (Coop, 2020). These effects of selection are only observable when the change is big enough to surpass the variance in the change introduced by genetic drift (Charlesworth, 2009), that is if the selection coefficient or the population size are sufficiently high. Lastly, gene flow will shift allele frequencies in a population by introducing new sets of alleles from other populations. This can counteract the effects of both selection and drift by lowering the adaptive and neutral genetic differentiation accumulated between the two populations (Lenormand, 2002; Yeaman & Otto, 2011). Overall, diverging lineages go through a dynamic tug of war between the differentiating forces of genetic drift and divergent selection against the homogenizing action of gene flow (Comeault & Matute, 2016; Roux et al., 2016). By studying how these forces act in natural populations, we can unravel the mechanisms behind the maintenance of genetic differentiation and speciation despite ongoing gene flow. Effective Population Size And Population Connectivity The demographic history of a population is characterized by changes in its size over time. In the context of evolution and conservation, population size can be described in two ways: the census size (Nc), which represents the total number of individuals, and the effective population size (Ne). The latter, introduced in early studies on genetic drift (Wright, 1931), is defined as the size of an idealized population that would exhibit the same level of genetic drift as the observed population. Although various methods for defining and calculating Ne exist (Caballero, 1994), it generally reflects the number of individuals in a population that will effectively contribute to the gene pool of the next generation. This means that Ne is usually much lower than Nc in natural populations (Frankham, 1995) and directly influences the rate of change in genetic composition and diversity of a population through genetic drift (Charlesworth, 2009). Effective population size is positively correlated with genetic diversity 8 and negatively correlated with inbreeding and the accumulation of deleterious mutations, three critical factors influencing a population's risk of extinction (Frankham, 1995). Population structure is another genetic pattern influenced by demographic history, as it affects genetic differentiation and reflects the connectivity among subpopulations (Gagnaire, 2020). Population subdivisions are ubiquitous in nature and can be induced by behavioural segregation (Nanda & Singh, 2012; Selz et al., 2014; Uy et al., 2018), by habitat fragmentation (Dias et al., 2013; Templeton et al., 1990) and simply by other factors reducing dispersal potential, like geographic distance (isolation by distance, Wright, 1943). By exchanging individuals between the subpopulations through migration, population structure directly impacts Ne. Factors such as the balance between migration and drift, asymmetries in subpopulation sizes, and unequal migration rates, can have drastic effects on the overall Ne of a metapopulation, which can either exceed or fall below the sum of the Ne values of individual subpopulations depending on these dynamics (Nunney, 1999; J. Wang & Caballero, 1999; Whitlock & Barton, 1997). Through the reconstruction of the demographic history of a species we can thus gain insights on how genetic drift and migration have shaped the genetic diversity within subpopulations and the genetic differentiation between them. At the same time, the current patterns of genetic variation can be used to infer the demographic history of the species and its populations (Loog, 2021; Lowe & Allendorf, 2010). Distinct aspects of genomic data can be harnessed to study the magnitude and timing of population size changes by recovering the demographic trajectories of a population at different timescales. Coalescent theory has been extensively used to reconstruct the demographic history of a population from just one (H. Li & Durbin, 2011) or multiple whole genome sequences (Palacios et al., 2015; Schiffels & Durbin, 2014; Terhorst et al., 2017). These methods provide fairly accurate estimates of effective population sizes over thousands of generations in the distant past and even estimate divergence times among lineages (Cahill et al., 2016). When the goal is the reconstruction of the most recent demographic trajectory, methods based leveraging the information of linkage disequilibrium patterns along the genomes of multiple individuals can be more useful (Corbin et al., 2012; Ragsdale & Gutenkunst, 2017; Santiago et al., 2020). The consequences of a population’s demographic history on the genetic health of a population can be directly measured in terms of its genetic diversity and inbreeding levels. Individual and population measures of genetic diversity can be used to assess the level of genetic erosion a population is suffering (Leroy et al., 2018). Homozygosity patterns along the genomes of multiple individuals can be investigated to estimate how inbred a particular population is at a given time (Ceballos et al., 2018) and even characterize the timing of inbreeding events (Bertrand et al., 2019; Druet & Gautier, 2022). From genomic data we can also infer ancestry composition and historical relationships among populations with the aid of a plethora of different methods (Lawson et al., 2012; Meisner & Albrechtsen, 2018; Pickrell & Pritchard, 2012; Pritchard et al., 2000), although some of these results should be interpreted with caution (Lawson et al., 2018). Inferring a population's demographic trajectory, its current genetic status, and its connectivity to other populations can provide essential information to understand its evolutionary history and evaluate its risk of extinction. 9 evolutionary and conservation dynamics. The Balkan lynx, representing the most endangered remnant population of Eurasian lynx, has experienced dramatic population declines, making genomic insights particularly valuable for guiding conservation efforts (Melovski et al., 2013, 2020, 2022). In contrast, the Caucasian lynx descends from a population that likely occupied a distinct glacial refugium (Tarkhnishvili et al., 2012; Vereshchagin, 1959), a region not included in prior studies. Both populations are recognized as unique, which led to their description as separate subspecies (Kitchener et al., 2017). By analyzing whole-genome sequences from these two key populations, we aim to understand their relationships with other Eurasian lynx populations, reconstruct their demographic histories, and evaluate their genetic health. We aim to understand the implications of our findings for the taxonomic classification of these populations and their importance for conservation management strategies. In the second chapter, we employ a combination of genome-environment association approaches to uncover genomic variation with potential adaptive significance in the Eurasian lynx. We identify genomic regions associated with local adaptation in different populations, determine the environmental variables and the specific genes involved in these processes, and characterize patterns of adaptive genetic divergence among populations. Again, we aim to discuss the implications of adaptive population structure for defining taxonomic and conservation units in the Eurasian lynx, providing insights to further guide assisted gene flow and genetic rescue initiatives. In the third chapter, we investigate genomic regions showing signals of interspecific introgression between Iberian and Eurasian lynx populations. Using demographic modelling informed by newly sequenced individuals and high-quality reference genomes, we reconstruct the recent and historical demographic histories of these populations. This modelling also provides realistic simulated data for training our introgression detection methods. We then identify and compare the extent of introgression across populations and genomic regions in both species. Additionally, we explore which Eurasian lynx population may have been most closely related to the extinct population that once coexisted with the Iberian lynx in Southwestern Europe. By examining what genomic regions are more permeable to introgression, we aim to understand the impact of introgression on local genetic diversity and its role in shaping evolutionary dynamics. REFERENCES Abascal, F., Corvelo, A., Cruz, F., Villanueva-Cañas, J. L., Vlasova, A., Marcet-Houben, M., Martínez-Cruz, B., Cheng, J. Y., Prieto, P., Quesada, V., Quilez, J., Li, G., García, F., Rubio-Camarillo, M., Frias, L., Ribeca, P., Capella-Gutiérrez, S., Rodríguez, J. M., Câmara, F., … Godoy, J. A. (2016). Extreme genomic erosion after recurrent demographic bottlenecks in the highly endangered Iberian lynx. Genome Biology, 17(1), 251. https://doi.org/10.1186/s13059-016-1090-1 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 16 Abdul-Muneer, P. M. (2014). Application of Microsatellite Markers in Conservation Genetics and Fisheries Management: Recent Advances in Population Structure Analysis and Conservation Strategies. Genetics Research International, 2014, 1–11. https://doi.org/10.1155/2014/691759 Agrawal, A. F., & Whitlock, M. C. (2012). Mutation Load: The Fitness of Individuals in Populations Where Deleterious Alleles Are Abundant. Annual Review of Ecology, Evolution, and Systematics, 43(1), 115–135. https://doi.org/10.1146/annurev-ecolsys-110411-160257 Aitken, S. N., Jordan, R., & Tumas, H. R. (2024). Conserving Evolutionary Potential: Combining Landscape Genomics with Established Methods to Inform Plant Conservation. Annual Review of Plant Biology, 75(1), 707–736. https://doi.org/10.1146/annurev-arplant-070523-044239 Aitken, S. N., & Whitlock, M. C. (2013). Assisted Gene Flow to Facilitate Local Adaptation to Climate Change. Annual Review of Ecology, Evolution, and Systematics, 44(1), 367–388. https://doi.org/10.1146/annurev-ecolsys-110512-135747 Alexander, D. H., Novembre, J., & Lange, K. (2009). Fast model-based estimation of ancestry in unrelated individuals. Genome Research, 19(9), 1655–1664. https://doi.org/10.1101/gr.094052.109 Anderson, E. C., & Thompson, E. A. (2002). A Model-Based Method for Identifying Species Hybrids Using Multilocus Genetic Data. Genetics, 160(3), 1217–1229. https://doi.org/10.1093/genetics/160.3.1217 Arnegard, M. E., McGee, M. D., Matthews, B., Marchinko, K. B., Conte, G. L., Kabir, S., Bedford, N., Bergek, S., Chan, Y. F., Jones, F. C., Kingsley, D. M., Peichel, C. L., & Schluter, D. (2014). Genetics of ecological divergence during speciation. Nature, 511(7509), 307–311. https://doi.org/10.1038/nature13301 Arnold, M. L., & Martin, N. H. (2009). Adaptation by introgression. Journal of Biology, 8(9), 82. https://doi.org/10.1186/jbiol176 Balkenhol, N., Dudaniec, R. Y., Krutovsky, K. V., Johnson, J. S., Cairns, D. M., Segelbacher, G., Selkoe, K. A., Von Der Heyden, S., Wang, I. J., Selmoni, O., & Joost, S. (2017). Landscape Genomics: Understanding Relationships Between Environmental Heterogeneity and Genomic Characteristics of Populations. In O. P. Rajora (Ed.), Population Genomics (pp. 261–322). Springer International Publishing. https://doi.org/10.1007/13836_2017_2 Barrett, R., & Schluter, D. (2008). Adaptation from standing genetic variation. Trends in Ecology & Evolution, 23(1), 38–44. https://doi.org/10.1016/j.tree.2007.09.008 Behzadi, F., Malekian, M., Fadakar, D., Adibi, M. A., & Bärmann, E. V. (2022). Phylogenetic analyses of Eurasian lynx (Lynx lynx Linnaeus, 1758) including new mitochondrial DNA sequences from Iran. Scientific Reports, 12(1), 3293. https://doi.org/10.1038/s41598-022-07369-z Bell, D. A., Robinson, Z. L., Funk, W. C., Fitzpatrick, S. W., Allendorf, F. W., Tallmon, D. A., & Whiteley, A. R. (2019). The Exciting Potential and Remaining Uncertainties of Genetic Rescue. Trends in Ecology & Evolution, 34(12), 1070–1079. 17 https://doi.org/10.1016/j.tree.2019.06.006 Bertrand, A. R., Kadri, N. K., Flori, L., Gautier, M., & Druet, T. (2019). RZooRoH: An R package to characterize individual genomic autozygosity and identify homozygous‐by‐descent segments. Methods in Ecology and Evolution, 10(6), 860–866. https://doi.org/10.1111/2041-210X.13167 Bierne, N., Lenormand, T., Bonhomme, F., & David, P. (2002). Deleterious mutations in a hybrid zone: Can mutational load decrease the barrier to gene flow? Genetical Research, 80(3), 197–204. https://doi.org/10.1017/S001667230200592X Breitenmoser, U., & Breitenmoser-Würsten, C. (2008). Der Luchs – Ein Grossraubtier in der Kulturlandschaft. Salm Verlag. Breitenmoser, U., Breitenmoser-Würsten, C., Lanz, T., von Arx, M., & Antonevich, A. (2015). Lynx lynx (errata version published in 2017). IUCN Red List of Threatened Species. https://www.iucnredlist.org/en Bürger, R., & Lynch, M. (1995). EVOLUTION AND EXTINCTION IN A CHANGING ENVIRONMENT: A QUANTITATIVE‐GENETIC ANALYSIS. Evolution, 49(1), 151–163. https://doi.org/10.1111/j.1558-5646.1995.tb05967.x Buri, P. (1956). Gene Frequency in Small Populations of Mutant Drosophila. Evolution, 10(4), 367. https://doi.org/10.2307/2406998 Bush, W. S., & Moore, J. H. (2012). Chapter 11: Genome-Wide Association Studies. PLoS Computational Biology, 8(12), e1002822. https://doi.org/10.1371/journal.pcbi.1002822 Caballero, A. (1994). Developments in the prediction of effective population size. Heredity, 73(6), 657–679. https://doi.org/10.1038/hdy.1994.174 Caballero‐Gómez, J., Rivero‐Juarez, A., Zorrilla, I., López, G., Nájera, F., Ulrich, R. G., Ruiz‐Rubio, C., Salcedo, J., Rivero, A., Paniagua, J., & García‐Bocanegra, I. (2022). Hepatitis E virus in the endangered Iberian lynx ( Lynx pardinus ). Transboundary and Emerging Diseases, 69(5). https://doi.org/10.1111/tbed.14624 Cahill, J. A., Soares, A. E. R., Green, R. E., & Shapiro, B. (2016). Inferring species divergence times using pairwise sequential Markovian coalescent modelling and low-coverage genomic data. Philosophical Transactions of the Royal Society B: Biological Sciences, 371(1699), 20150138. https://doi.org/10.1098/rstb.2015.0138 Capblancq, T., & Forester, B. R. (2021). Redundancy analysis: A Swiss Army Knife for landscape genomics. Methods in Ecology and Evolution, 12(12), 2298–2309. https://doi.org/10.1111/2041-210X.13722 Casas-Marce, M., Marmesat, E., Soriano, L., Martínez-Cruz, B., Lucena-Perez, M., Nocete, F., Rodríguez-Hidalgo, A., Canals, A., Nadal, J., Detry, C., Bernáldez-Sánchez, E., Fernández-Rodríguez, C., Pérez-Ripoll, M., Stiller, M., Hofreiter, M., Rodríguez, A., Revilla, E., Delibes, M., & Godoy, J. A. (2017). Spatiotemporal Dynamics of Genetic Variation in the Iberian Lynx along Its Path to Extinction Reconstructed with Ancient DNA. Molecular Biology and Evolution, 34(11), 2893–2907. https://doi.org/10.1093/molbev/msx222 18 Ceballos, F. C., Joshi, P. K., Clark, D. W., Ramsay, M., & Wilson, J. F. (2018). Runs of homozygosity: Windows into population history and trait architecture. Nature Reviews Genetics, 19(4), 220–234. https://doi.org/10.1038/nrg.2017.109 Chapron, G., Kaczensky, P., Linnell, J. D. C., von Arx, M., Huber, D., Andrén, H., López-Bao, J. V., Adamec, M., Álvares, F., Anders, O., Balčiauskas, L., Balys, V., Bedő, P., Bego, F., Blanco, J. C., Breitenmoser, U., Brøseth, H., Bufka, L., Bunikyte, R., … Boitani, L. (2014). Recovery of large carnivores in Europe’s modern human-dominated landscapes. Science, 346(6216), 1517–1519. https://doi.org/10.1126/science.1257553 Charlesworth, B. (2009). Effective population size and patterns of molecular evolution and variation. Nature Reviews Genetics, 10(3), 195–205. https://doi.org/10.1038/nrg2526 Clarke, J. G., Smith, A. C., & Cullingham, C. I. (2024). Genetic rescue often leads to higher fitness as a result of increased heterozygosity across animal taxa. Molecular Ecology, 33(19), e17532. https://doi.org/10.1111/mec.17532 Clavero, M., & Delibes, M. (2013). Using historical accounts to set conservation baselines: The case of Lynx species in Spain. Biodiversity and Conservation, 22(8), 1691–1702. https://doi.org/10.1007/s10531-013-0506-4 Comeault, A. A., & Matute, D. R. (2016). Reinforcement’s incidental effects on reproductive isolation between conspecifics. Current Zoology, 62(2), 135–143. https://doi.org/10.1093/cz/zow002 Coop, G. (2020). Population and Quantitative Genetics (Third Edition). The Coop Lab, University of California, Davis. https://github.com/cooplab/popgen-notes Corbin, L. J., Liu, A. Y. H., Bishop, S. C., & Woolliams, J. A. (2012). Estimation of historical effective population size using linkage disequilibria with marker data. Journal of Animal Breeding and Genetics, 129(4), 257–270. https://doi.org/10.1111/j.1439-0388.2012.01003.x Crandall, K. A., Bininda-Emonds, O. R. P., Mace, G. M., & Wayne, R. K. (2000). Considering evolutionary processes in conservation biology. Trends in Ecology & Evolution, 15(7), 290–295. https://doi.org/10.1016/S0169-5347(00)01876-0 Crow, J. F. (1948). ALTERNATIVE HYPOTHESES OF HYBRID VIGOR. Genetics, 33(5), 477–487. https://doi.org/10.1093/genetics/33.5.477 Darul, R., Gavashelishvili, A., Saveljev, A. P., Seryodkin, I. V., Linnell, J. D. C., Okarma, H., Bagrade, G., Ornicans, A., Ozolins, J., Männil, P., Khorozyan, I., Melovski, D., Stojanov, A., Trajçe, A., Hoxha, B., Dvornikov, M. G., Galsandorj, N., Okhlopkov, I., Mamuchadze, J., … Schmidt, K. (2022). Coat Polymorphism in Eurasian Lynx: Adaptation to Environment or Phylogeographic Legacy? Journal of Mammalian Evolution, 29(1), 51–62. https://doi.org/10.1007/s10914-021-09580-7 DeGiorgio, M., Huber, C. D., Hubisz, M. J., Hellmann, I., & Nielsen, R. (2016). S weep F inder 2: Increased sensitivity, robustness and flexibility. Bioinformatics, 32(12), 1895–1897. https://doi.org/10.1093/bioinformatics/btw051 DeGiorgio, M., & Szpiech, Z. A. (2022). A spatially aware likelihood test to detect sweeps 19 from haplotype distributions. PLOS Genetics, 18(4), e1010134. https://doi.org/10.1371/journal.pgen.1010134 Delibes, M. (1980). Feeding Ecology of the Spanish Lynx in the Coto Doñana. Acta Theriologica, 25(24), 309–324. Dias, M. S., Cornu, J., Oberdorff, T., Lasso, C. A., & Tedesco, P. A. (2013). Natural fragmentation in river networks as a driver of speciation for freshwater fishes. Ecography, 36(6), 683–689. https://doi.org/10.1111/j.1600-0587.2012.07724.x Druet, T., & Gautier, M. (2022). A hidden Markov model to estimate homozygous-by-descent probabilities associated with nested layers of ancestors. Theoretical Population Biology, 145, 38–51. https://doi.org/10.1016/j.tpb.2022.03.001 Duforet-Frebourg, N., Bazin, E., & Blum, M. G. B. (2014). Genome Scans for Detecting Footprints of Local Adaptation Using a Bayesian Factor Model. Molecular Biology and Evolution, 31(9), 2483–2495. https://doi.org/10.1093/molbev/msu182 Durvasula, A., & Sankararaman, S. (2019). A statistical model for reference-free inference of archaic local ancestry. PLOS Genetics, 15(5), e1008175. https://doi.org/10.1371/journal.pgen.1008175 Fan, R., Lange, K., Zhang, Y., & Zhong, M. (2011). A cross-population extended haplotype-based homozygosity score test to detect positive selection in genome-wide scans. Statistics and Its Interface, 4(1), 51–63. https://doi.org/10.4310/SII.2011.v4.n1.a6 Fariello, M. I., Boitard, S., Naya, H., SanCristobal, M., & Servin, B. (2013). Detecting Signatures of Selection Through Haplotype Differentiation Among Hierarchically Structured Populations. Genetics, 193(3), 929–941. https://doi.org/10.1534/genetics.112.147231 Fay, J. C., & Wu, C.-I. (2000). Hitchhiking Under Positive Darwinian Selection. Genetics, 155(3), 1405–1413. https://doi.org/10.1093/genetics/155.3.1405 Fitzpatrick, M. C., & Keller, S. R. (2015). Ecological genomics meets community-level modelling of biodiversity: Mapping the genomic landscape of current and future environmental adaptation. Ecology Letters, 18(1), 1–16. https://doi.org/10.1111/ele.12376 Forester, B. R., Lasky, J. R., Wagner, H. H., & Urban, D. L. (2018). Comparing methods for detecting multilocus adaptation with multivariate genotype-environment associations. Molecular Ecology, 27(9), 2215–2233. https://doi.org/10.1111/mec.14584 Frankham, R. (1995). Effective population size/adult population size ratios in wildlife: A review. Genetical Research, 66(2), 95–107. https://doi.org/10.1017/S0016672300034455 Frankham, R. (2005). Genetics and extinction. Biological Conservation, 126(2), 131–140. https://doi.org/10.1016/j.biocon.2005.05.002 Frankham, R., Ballou, J. D., Eldridge, M. D. B., Lacy, R. C., Ralls, K., Dudash, M. R., & Fenster, C. B. (2011). Predicting the Probability of Outbreeding Depression: Predicting Outbreeding Depression. Conservation Biology, 25(3), 465–475. 20 https://doi.org/10.1111/j.1523-1739.2011.01662.x Fraser, D. J., & Bernatchez, L. (2001). Adaptive evolutionary conservation: Towards a unified concept for defining conservation units. Molecular Ecology, 10(12), 2741–2752. https://doi.org/10.1046/j.1365-294X.2001.t01-1-01411.x Funk, W. C., McKay, J. K., Hohenlohe, P. A., & Allendorf, F. W. (2012). Harnessing genomics for delineating conservation units. Trends in Ecology & Evolution, 27(9), 489–496. https://doi.org/10.1016/j.tree.2012.05.012 Gagnaire, P. (2020). Comparative genomics approach to evolutionary process connectivity. Evolutionary Applications, 13(6), 1320–1334. https://doi.org/10.1111/eva.12978 Galtier, N. (2024). An approximate likelihood method reveals ancient gene flow between human, chimpanzee and gorilla. Peer Community Journal, 4, e3. https://doi.org/10.24072/pcjournal.359 Gautier, M. (2015). Genome-Wide Scan for Adaptive Divergence and Association with Population-Specific Covariates. Genetics, 201(4), 1555–1579. https://doi.org/10.1534/genetics.115.181453 Geneva, A. J., Muirhead, C. A., Kingan, S. B., & Garrigan, D. (2015). A New Method to Scan Genomes for Introgression in a Secondary Contact Model. PLOS ONE, 10(4), e0118621. https://doi.org/10.1371/journal.pone.0118621 Gil-Sánchez, J. M., Ballesteros-Duperón, E., & Bueno-Segura, J. F. (2006). Feeding ecology of the Iberian lynxLynx pardinus in eastern Sierra Morena (Southern Spain). Acta Theriologica, 51(1), 85–90. https://doi.org/10.1007/BF03192659 Green, R. E., Krause, J., Briggs, A. W., Maricic, T., Stenzel, U., Kircher, M., Patterson, N., Li, H., Zhai, W., Fritz, M. H.-Y., Hansen, N. F., Durand, E. Y., Malaspinas, A.-S., Jensen, J. D., Marques-Bonet, T., Alkan, C., Prüfer, K., Meyer, M., Burbano, H. A., … Pääbo, S. (2010). A Draft Sequence of the Neandertal Genome. Science, 328(5979), 710–722. https://doi.org/10.1126/science.1188021 Grummer, J. A., Booker, T. R., Matthey‐Doret, R., Nietlisbach, P., Thomaz, A. T., & Whitlock, M. C. (2022). The immediate costs and long‐term benefits of assisted gene flow in large populations. Conservation Biology. https://doi.org/10.1111/cobi.13911 Guan, Y. (2014). Detecting Structure of Haplotypes and Local Ancestry. Genetics, 196(3), 625–642. https://doi.org/10.1534/genetics.113.160697 Hahn, M. W., & Hibbins, M. S. (2019). A Three-Sample Test for Introgression. Molecular Biology and Evolution, 36(12), 2878–2882. https://doi.org/10.1093/molbev/msz178 Harris, K., & Nielsen, R. (2016). The Genetic Cost of Neanderthal Introgression. Genetics, 203(2), 881–891. https://doi.org/10.1534/genetics.116.186890 Hedrick, P. W. (2013). Adaptive introgression in animals: Examples and comparison to new mutation and standing variation as sources of adaptive variation. Molecular Ecology, 22(18), 4606–4618. https://doi.org/10.1111/mec.12415 Hedrick, P. W., & Fredrickson, R. (2010). Genetic rescue guidelines with examples from Mexican wolves and Florida panthers. Conservation Genetics, 11(2), 615–626. 21 https://doi.org/10.1007/s10592-009-9999-5 Hedrick, P. W., & Garcia-Dorado, A. (2016). Understanding Inbreeding Depression, Purging, and Genetic Rescue. Trends in Ecology & Evolution, 31(12), 940–952. https://doi.org/10.1016/j.tree.2016.09.005 Hibbins, M. S., & Hahn, M. W. (2022). Phylogenomic approaches to detecting and characterizing introgression. Genetics, 220(2), iyab173. https://doi.org/10.1093/genetics/iyab173 Hoban, S., Kelley, J. L., Lotterhos, K. E., Antolin, M. F., Bradburd, G., Lowry, D. B., Poss, M. L., Reed, L. K., Storfer, A., & Whitlock, M. C. (2016). Finding the Genomic Basis of Local Adaptation: Pitfalls, Practical Solutions, and Future Directions. The American Naturalist, 188(4), 379–397. https://doi.org/10.1086/688018 Hohenlohe, P. A., Funk, W. C., & Rajora, O. P. (2021). Population genomics for wildlife conservation and management. Molecular Ecology, 30(1), 62–82. https://doi.org/10.1111/mec.15720 Huerta-Sánchez, E., Jin, X., Asan, Bianba, Z., Peter, B. M., Vinckenbosch, N., Liang, Y., Yi, X., He, M., Somel, M., Ni, P., Wang, B., Ou, X., Huasang, Luosang, J., Cuo, Z. X. P., Li, K., Gao, G., Yin, Y., … Nielsen, R. (2014). Altitude adaptation in Tibetans caused by introgression of Denisovan-like DNA. Nature, 512(7513), 194–197. https://doi.org/10.1038/nature13408 Ingvarsson, P. K., & Whitlock, M. C. (2000). Heterosis increases the effective migration rate. Proceedings of the Royal Society of London. Series B: Biological Sciences, 267(1450), 1321–1326. https://doi.org/10.1098/rspb.2000.1145 Johnson, W. E., Eizirik, E., Pecon-Slattery, J., Murphy, W. J., Antunes, A., Teeling, E., & O’Brien, S. J. (2006). The Late Miocene Radiation of Modern Felidae: A Genetic Assessment. Science, 311(5757), 73–77. https://doi.org/10.1126/science.1122277 Joly, S., McLenachan, P. A., & Lockhart, P. J. (2009). A Statistical Approach for Distinguishing Hybridization and Incomplete Lineage Sorting. The American Naturalist, 174(2), E54–E70. https://doi.org/10.1086/600082 Jones, M. R., Mills, L. S., Alves, P. C., Callahan, C. M., Alves, J. M., Lafferty, D. J. R., Jiggins, F. M., Jensen, J. D., Melo-Ferreira, J., & Good, J. M. (2018). Adaptive introgression underlies polymorphic seasonal camouflage in snowshoe hares. Science, 360(6395), 1355–1358. https://doi.org/10.1126/science.aar5273 Kaczensky, P., Chapron, G., Von Arx, M., Huber, D., Andrén, H., & Linnell, J. (2013). Status, management and distribution of large carnivores-bear, lynx, wolf & wolverine-in Europe. Verlag nicht ermittelbar. Kaplan, N. L., Hudson, R. R., & Langley, C. H. (1989). The “hitchhiking effect” revisited. Genetics, 123(4), 887–899. https://doi.org/10.1093/genetics/123.4.887 Kim, B. Y., Huber, C. D., & Lohmueller, K. E. (2018). Deleterious variation shapes the genomic landscape of introgression. PLOS Genetics, 14(10), e1007741. https://doi.org/10.1371/journal.pgen.1007741 Kim, Y., & Nielsen, R. (2004). Linkage Disequilibrium as a Signature of Selective Sweeps. 22 Genetics, 167(3), 1513–1524. https://doi.org/10.1534/genetics.103.025387 Kitchener, A. C., Breitenmoser-Würsten, C., Eizirik, E., Gentry, A., Werdelin, L., Wilting, A., Yamaguchi, N., Abramov, A. V., Christiansen, P., Driscoll, C., Duckworth, J. W., Johnson, W. E., Luo, S. J., Meijaard, E., O’Donoghue, P., Sanderson, J., Seymour, K., Bruford, M., Groves, C., … Tobe, S. (2017). A revised taxonomy of the Felidae: The final report of the Cat Classification Task Force of the IUCN Cat Specialist Group. Cat News Special Issue, 11, 80. Klassmann, A., & Gautier, M. (2022). Detecting selection using extended haplotype homozygosity (EHH)-based statistics in unphased or unpolarized data. PLOS ONE, 17(1), e0262024. https://doi.org/10.1371/journal.pone.0262024 Kleinman-Ruiz, D., Lucena-Perez, M., Villanueva, B., Fernández, J., Saveljev, A. P., Ratkiewicz, M., Schmidt, K., Galtier, N., García-Dorado, A., & Godoy, J. A. (2022). Purging of deleterious burden in the endangered Iberian lynx. Proceedings of the National Academy of Sciences, 119(11), e2110614119. https://doi.org/10.1073/pnas.2110614119 Kratochvil, J. (1968). Survey of the distribution of populations of the genus Lynx in Europe. Acta Sc. Nat. Brno, 4, 5–12. Kurtén, B., & Granqvist, E. (1987). Fossil pardel lynx (Lynx pardina spelaea Boule) from a cave in southern France. Annales Zoologici Fennici, 24(1), 39–43. Lande, R. (1994). RISK OF POPULATION EXTINCTION FROM FIXATION OF NEW DELETERIOUS MUTATIONS. Evolution, 48(5), 1460–1469. https://doi.org/10.1111/j.1558-5646.1994.tb02188.x Lawson, D. J., Hellenthal, G., Myers, S., & Falush, D. (2012). Inference of Population Structure using Dense Haplotype Data. PLoS Genetics, 8(1), e1002453. https://doi.org/10.1371/journal.pgen.1002453 Lawson, D. J., Van Dorp, L., & Falush, D. (2018). A tutorial on how not to over-interpret STRUCTURE and ADMIXTURE bar plots. Nature Communications, 9(1), 3258. https://doi.org/10.1038/s41467-018-05257-7 Layton, K. K. S., Brieuc, M. S. O., Castilho, R., Diaz‐Arce, N., Estévez‐Barcia, D., Fonseca, V. G., Fuentes‐Pardo, A. P., Jeffery, N. W., Jiménez‐Mena, B., Junge, C., Kaufmann, J., Leinonen, T., Maes, S. M., McGinnity, P., Reed, T. E., Reisser, C. M. O., Silva, G., Vasemägi, A., & Bradbury, I. R. (2024). Predicting the future of our oceans—Evaluating genomic forecasting approaches in marine species. Global Change Biology, 30(3), e17236. https://doi.org/10.1111/gcb.17236 Lenormand, T. (2002). Gene flow and the limits to natural selection. Trends in Ecology & Evolution, 17(4), 183–189. https://doi.org/10.1016/S0169-5347(02)02497-7 Leroy, G., Carroll, E. L., Bruford, M. W., DeWoody, J. A., Strand, A., Waits, L., & Wang, J. (2018). Next‐generation metrics for monitoring genetic erosion within populations of conservation concern. Evolutionary Applications, 11(7), 1066–1083. https://doi.org/10.1111/eva.12564 Li, G., Davis, B. W., Eizirik, E., & Murphy, W. J. (2016). Phylogenomic evidence for ancient 23 hybridization in the genomes of living cats (Felidae). Genome Research, 26(1), 1–11. https://doi.org/10.1101/gr.186668.114 Li, G., Figueiró, H. V., Eizirik, E., & Murphy, W. J. (2019). Recombination-Aware Phylogenomics Reveals the Structured Genomic Landscape of Hybridizing Cat Species. Molecular Biology and Evolution, 36(10), 2111–2126. https://doi.org/10.1093/molbev/msz139 Li, H., & Durbin, R. (2011). Inference of human population history from individual whole-genome sequences. Nature, 475(7357), 493–496. https://doi.org/10.1038/nature10231 Lindtke, D., & Buerkle, C. A. (2015). The genetic architecture of hybrid incompatibilities and their effect on barriers to introgression in secondary contact: THE GENETIC ARCHITECTURE OF HYBRID INCOMPATIBILITIES. Evolution, 69(8), 1987–2004. https://doi.org/10.1111/evo.12725 Loog, L. (2021). Sometimes hidden but always there: The assumptions underlying genetic inference of demographic histories. Philosophical Transactions of the Royal Society B: Biological Sciences, 376(1816), 20190719. https://doi.org/10.1098/rstb.2019.0719 Lowe, W. H., & Allendorf, F. W. (2010). What can genetics tell us about population connectivity?: GENETIC AND DEMOGRAPHIC CONNECTIVITY. Molecular Ecology, 19(15), 3038–3051. https://doi.org/10.1111/j.1365-294X.2010.04688.x Lucena-Perez, M., Bazzicalupo, E., Paijmans, J., Kleinman-Ruiz, D., Dalén, L., Hofreiter, M., Delibes, M., Clavero, M., & Godoy, J. A. (2022). Ancient genome provides insights into the history of Eurasian lynx in Iberia and Western Europe. Quaternary Science Reviews, 285, 107518. https://doi.org/10.1016/j.quascirev.2022.107518 Lucena‐Perez, M., Kleinman‐Ruiz, D., Marmesat, E., Saveljev, A. P., Schmidt, K., & Godoy, J. A. (2021). Bottleneck‐associated changes in the genomic landscape of genetic diversity in wild lynx populations. Evolutionary Applications, 14(11), 2664–2679. https://doi.org/10.1111/eva.13302 Lucena‐Perez, M., Marmesat, E., Kleinman‐Ruiz, D., Martínez‐Cruz, B., Węcek, K., Saveljev, A. P., Seryodkin, I. V., Okhlopkov, I., Dvornikov, M. G., Ozolins, J., Galsandorj, N., Paunovic, M., Ratkiewicz, M., Schmidt, K., & Godoy, J. A. (2020). Genomic patterns in the widespread Eurasian lynx shaped by Late Quaternary climatic fluctuations and anthropogenic impacts. Molecular Ecology, 29(4), 812–828. https://doi.org/10.1111/mec.15366 Lucena-Perez, M., Paijmans, J. L. A., Nocete, F., Nadal, J., Detry, C., Dalén, L., Hofreiter, M., Barlow, A., & Godoy, J. A. (2024). Recent increase in species-wide diversity after interspecies introgression in the highly endangered Iberian lynx. Nature Ecology & Evolution, 8(2), 282–292. https://doi.org/10.1038/s41559-023-02267-7 Mable, B. K. (2019). Conservation of adaptive potential and functional diversity: Integrating old and new approaches. Conservation Genetics, 20(1), 89–100. https://doi.org/10.1007/s10592-018-1129-9 Mallet, J. (2005). Hybridization as an invasion of the genome. Trends in Ecology & Evolution, 20(5), 229–237. https://doi.org/10.1016/j.tree.2005.02.010 24 Mallet, J., Besansky, N., & Hahn, M. W. (2016). How reticulated are species? BioEssays, 38(2), 140–149. https://doi.org/10.1002/bies.201500149 Manel, S., Joost, S., Epperson, B. K., Holderegger, R., Storfer, A., Rosenberg, M. S., Scribner, K. T., Bonin, A., & Fortin, M. (2010). Perspectives on the use of landscape genetics to detect genetic adaptive variation in the field. Molecular Ecology, 19(17), 3760–3772. https://doi.org/10.1111/j.1365-294X.2010.04717.x Martínez, F., Manteca, X., & Pastor, J. (2013). RETROSPECTIVE STUDY OF MORBIDITY AND MORTALITY OF CAPTIVE IBERIAN LYNX ( LYNX PARDINUS ) IN THE EX SITU CONSERVATION PROGRAMME (2004–JUNE 2010). Journal of Zoo and Wildlife Medicine, 44(4), 845–852. https://doi.org/10.1638/2011-0165R4.1 Matute, D. R., Butler, I. A., Turissini, D. A., & Coyne, J. A. (2010). A Test of the Snowball Theory for the Rate of Evolution of Hybrid Incompatibilities. Science, 329(5998), 1518–1521. https://doi.org/10.1126/science.1193440 Matyushkin, E. N., & Vaisfeld, M. A. (Eds.). (2003). The Lynx – Regional Features of Ecology, Use and Protection. Nauka. Mecozzi, B., Sardella, R., Boscaini, A., Cherin, M., Costeur, L., Madurell-Malapeira, J., Pavia, M., Profico, A., & Iurino, D. A. (2021). The tale of a short-tailed cat: New outstanding Late Pleistocene fossils of Lynx pardinus from southern Italy. Quaternary Science Reviews, 262, 106840. https://doi.org/10.1016/j.quascirev.2021.106840 Meisner, J., & Albrechtsen, A. (2018). Inferring Population Structure and Admixture Proportions in Low-Depth NGS Data. Genetics, 210(2), 719–731. https://doi.org/10.1534/genetics.118.301336 Meli, M. L., Cattori, V., Martínez, F., López, G., Vargas, A., Palomares, F., López-Bao, J. V., Hofmann-Lehmann, R., & Lutz, H. (2010). Feline leukemia virus infection: A threat for the survival of the critically endangered Iberian lynx (Lynx pardinus). Veterinary Immunology and Immunopathology, 134(1–2), 61–67. https://doi.org/10.1016/j.vetimm.2009.10.010 Melovski, D., Ivanov, G., Stojanov, A., Avukatov, V., Trajce, A., Hoxha, B., von Arx, M., Breitenmoser, C., Hristovski, S., Shumka, S., & Breitenmoser, U. (2013). Distribution and conservation status of the Balkan lynx (Lynx lynx balcanicus Bureš, 1941). Proceedings of the 4th Congress of Ecologists of Macedonia with International Participation, 50–60. Melovski, D., Krofel, M., Avukatov, V., Fležar, U., Gonev, A., Hočevar, L., Ivanov, G., Leschinski, L., Pavlov, A., Stojanov, A., Veapi, E., & Mengüllüoğlu, D. (2022). Diverging ecological traits between the Balkan lynx and neighbouring populations as a basis for planning its genetic rescue. Mammalian Biology. https://doi.org/10.1007/s42991-022-00268-w Melovski, D., von Arx, M., Avukatov, V., Breitenmoser-Würsten, C., Đurović, M., Elezi, R., Gimenez, O., Hoxha, B., Hristovski, S., Ivanov, G., Karamanlidis, A. A., Lanz, T., Mersini, K., Perović, A., Ramadani, A., Sanaja, B., Sanaja, P., Schwaderer, G., Spangenberg, A., … Breitenmoser, U. (2020). Using questionnaire surveys and occupancy modelling to identify conservation priorities for the Critically Endangered 25 https://doi.org/10.1111/j.0014-3820.2000.tb01232.x Whitlock, M. C., & Barton, N. H. (1997). The Effective Size of a Subdivided Population. Genetics, 146(1), 427–441. https://doi.org/10.1093/genetics/146.1.427 Whitlock, M. C., Ingvarsson, P. K., & Hatfield, T. (2000). Local drift load and the heterosis of interconnected populations. Heredity, 84(4), 452–457. https://doi.org/10.1046/j.1365-2540.2000.00693.x Whitney, K. D., Randell, R. A., & Rieseberg, L. H. (2006). Adaptive Introgression of Herbivore Resistance Traits in the Weedy Sunflower Helianthus annuus. The American Naturalist, 167(6), 794–807. https://doi.org/10.1086/504606 Willi, Y., & Hoffmann, A. A. (2009). Demographic factors and genetic variation influence population persistence under environmental change. Journal of Evolutionary Biology, 22(1), 124–133. https://doi.org/10.1111/j.1420-9101.2008.01631.x Wright, S. (1931). EVOLUTION IN MENDELIAN POPULATIONS. Genetics, 16(2), 97–159. https://doi.org/10.1093/genetics/16.2.97 Wright, S. (1943). Isolation by Distance. Genetics, 28(2), 114–138. https://doi.org/10.1093/genetics/28.2.114 Wright, S. (1949). THE GENETICAL STRUCTURE OF POPULATIONS. Annals of Eugenics, 15(1), 323–354. https://doi.org/10.1111/j.1469-1809.1949.tb02451.x Yeaman, S., & Otto, S. P. (2011). ESTABLISHMENT AND MAINTENANCE OF ADAPTIVE GENETIC DIVERGENCE UNDER MIGRATION, SELECTION, AND DRIFT. Evolution, 65(7), 2123–2129. https://doi.org/10.1111/j.1558-5646.2011.01277.x Zhang, X., Kim, B., Lohmueller, K. E., & Huerta-Sánchez, E. (2020). The Impact of Recessive Deleterious Variation on Signals of Adaptive Introgression in Human Populations. Genetics, 215(3), 799–812. https://doi.org/10.1534/genetics.120.303081 32 CHAPTER 1: History, demography and genetic status of Balkan and Caucasian Lynx lynx (Linnaeus, 1758) populations revealed by genome-wide variation Enrico Bazzicalupo, Maria Lucena-Perez, Daniel Kleinman-Ruiz, Aleksandar Pavlov, Aleksandër Trajçe, Bledi Hoxha, Bardh Sanaja, Zurab Gurielidze, Niko Kerdikoshvili, Jimsher Mamuchadze, Yuriy A. Yarovenko, Muzigit I. Akkiev, Mirosław Ratkiewicz, Alexander P. Saveljev, Dime Melovski, Alexander Gavashelishvili, Krzysztof Schmidt, José A. Godoy 33 ABSTRACT Genome-wide genetic data can provide key input for both taxonomy and conservation, but its use in this context remains limited. In this study, we performed the first genome-wide assessment of genetic variation in two populations of the Eurasian lynx, the Balkan population, the most threatened, and the Caucasian population, a possible glacial refugium, with the aim to place them in the context of the species, investigate their demographic history and evaluate their genetic status. We obtained whole genome resequencing data from seven Balkan and 12 Caucasian lynx, and analysed them along with novel and existing data from other populations. Based on a total 105 whole genome and 114 mitogenome sequences, we reconstructed phylogenetic and historical relationships, ancient and recent demography, and patterns of genetic diversity and inbreeding. Both the Balkan and the Caucasian lynx appear as distinct mitochondrial lineages that diverged from the rest of the Eurasian lynx lineages ca. 92.6 kya, and from each other ca. 46.4 kya. Autosomal data suggest, however, that the Balkan lynx is closely related with the Carpathian population, revealing alarmingly low genetic diversity and high inbreeding. In contrast, the Caucasian lynx shows a longer history of relative isolation from the rest of lynx populations and high genetic diversity, consistent with its large long-term effective population size. The taxonomic status of the Balkan lynx remains unresolved due to the evidence of long-term isolation in the mitogenome, contrasting with extensive autosomal admixture and intense recent genetic drift in the nuclear genome. Our results alert on genetic risks and call for the consideration of genetic rescue from closely related Carpathian lynxes. In contrast, substantial mitogenomic and autosomal divergence with no signs of genetic drift supports the identification of the Caucasian lynx as a separate subspecies with good genetic health. INTRODUCTION Widespread mammalian carnivores face particular conservation and taxonomic challenges. For millennia they have been susceptible to human impacts through direct persecution and habitat destruction and fragmentation, with these pressures peaking in intensity over the last few centuries (Ripple et al., 2014), showing signs of recovery only in the last few decades (Chapron et al., 2014). They often occupy diverse habitats and broad climatic ranges, and present subtle morphological and ecological differences, which have led to the recognition of multiple subspecies in some cases (e.g. Uphyrkina et al., 2001; Wilting et al., 2015). One notable example of a widespread carnivore is the Eurasian lynx (Lynx lynx Linnaeus, 1758), whose range extends across most of Eurasia, encompasses a wide range of habitats and climates, and includes populations with contrasting past and recent demography (Lucena-Perez et al., 2020; Matyushkin et al., 2003). Our recent genome-wide analysis of samples distributed across most of the species’ range (Lucena-Perez et al., 2020) revealed several mitogenomic lineages of L. lynx with divergence dates ranging from ca. 100 to 17 kya, likely originating during periods of isolation in separate glacial refugia. However, nuclear data suggests that most of these lineages were subsequently admixed through male-biased gene flow, resulting in a shallower nuclear structure, where 34 Asia-Europe appears as the deeper separation. Asian and European L. lynx populations started to diverge around 100 kya, prompted by the onset of a period of continuous and widespread decline in population size which intensified during the last millennia, leading to the increasing fragmentation and isolation of populations. Several remnant populations in Europe (i.e. Norway, NE Poland and the Carpathians) were identified as genetically eroded due to sustained decline and isolation in the last few centuries. The Balkan lynx (Lynx lynx balcanicus Bureš, 1941) has been considered the most threatened autochthonous subspecies of the Eurasian lynx. Recently classified as critically endangered by IUCN, its numbers are estimated to be a mere 50 or less mature individuals (Melovski et al., 2015). Its current distribution range is confined to the south-west Balkan Peninsula, with core areas along the Macedonian-Albanian border (Melovski et al., 2020; Figure 1.1). The evolutionary history of the Balkan lynx since the post-glacial period was first postulated by Hemmer (1993), who suggested a post-glacial re-colonization from the Balkans to southern Scandinavia, thus suggesting a shared history of Balkan, Carpathian and Scandinavian populations. Based on the mitochondrial DNA of specimens from Skopje, Macedonia dated in the 20th century, Gugolz et al. (2008) found two Balkan lynx haplotypes, which are not known to the other European populations, one unique to the Balkans, and the other shared with the Caucasus. More recently, whole mitogenome sequences revealed a highly divergent mitogenomic lineage exclusive to the Balkans, reinforcing the idea of a separate southern refugium (Lucena-Perez et al., 2020). However, the nuclear genetic status of the relict Balkan lynx and its relationships with other populations or subspecies were left unassessed due to insufficient sampling and low-quality sequencing, allowing only mitogenome reconstructions. The Caucasian lynx (Lynx lynx dinniki, Satunin, 1915) is another putative subspecies of Eurasian lynx that occurs in the Caucasus region and extends into Turkey, Iraq and Iran (Breitenmoser et al., 2015; Figure 1.1), and for which high mitochondrial haplotype diversity typical of glacial refugia was previously observed (Rueness et al., 2014). Because the Caucasus region is well known as a biodiversity hotspot and late Pleistocene refugium for a number of taxa (Tarkhnishvili et al., 2012; Vereshchagin, 1959), it is likely that the Caucasian lynx has a long independent evolutionary history. However, previous studies using mitochondrial and microsatellite markers have left the taxonomic status of Caucasian lynx lineage and its relationships to the rest of the species unresolved (Gugolz et al., 2008; Rueness et al., 2014). This subspecies was recently found to be paraphyletic at the mitogenome with respect to L. l. balcanicus and divergent with respect to other Eurasian lynx (Mengüllüoğlu et al., 2021). Lynx abundance in the region used to be high and the species was treated as a natural resource, as in other parts of the Soviet Union, with kill numbers showing a decreasing trend during most of the 20th century. A general recovery since the 1980s has been reported, with numbers being still low in the Russian and Turkish Caucasus, but lynx being quite widespread across protected areas in Armenia, Azerbaijan, Georgia and Iran (Askerov et al., 2020). 35 Figure 1.1 Sampling and distribution map of Balkan and Caucasian lynx. Extant and possibly extant distributions of L. l. balcanicus and L. l. dinniki are colour-coded. Coordinates of sampled individuals are marked by yellow (Balkans) and blue (Caucasian) dots. Spatial geographic range downloaded from IUCN (Breitenmoser et al., 2015). Both the Balkan and the Caucasian lynx are provisionally listed as separate subspecies L. l. balcanicus and L. l. dinniki, respectively, in the latest taxonomic review of the Felidae, based on morphological, biogeographical and limited genetic information, and the same status applies to the Carpathian lynx, L. l. carpathicus (Kitchener et al., 2017). Some doubts remain on whether balcanicus should become part of dinniki, as they shared mitochondrial haplotypes (Gugolz et al., 2008), and the balcanicus mitogenomes are descendant of the common ancestor of dinniki mitogenomes (Mengüllüoğlu et al., 2021), and since the two southern refugia were likely connected by land bridges during the last Ice Age. The delimitation of carpathicus is supported by the possibility of an extra-Mediterranean forest refugium in the region, although its limited mitogenomic variation, well embedded within that of Polish, Baltic and western Russia lynxes, argues otherwise (Lucena-Perez et al., 2020). A deeper understanding of the evolutionary and demographic history of lynx populations at the southwestern part of the species range is thus needed to resolve their taxonomic and conservation status, and to gain a more complete understanding of the evolutionary history of the species as a whole. Genome-wide genetic data provide unique insights into these issues with unprecedented power and resolution, allowing for the reconstruction of events such as the onset of population isolation and structure, demographic declines/expansions, admixture and selection. This information is particularly relevant for taxonomy and conservation, as long-term isolation, genetic distinctiveness and local adaptation can support the delineation of subspecies and conservation units. However, the application of genomics in this context is not 36 without hurdles, challenges and controversies, including the relative weight to give to historical versus adaptive processes, the lack of objective and general divergence thresholds, and the peril of granting evolutionary and taxonomic significance to divergence due to genetic drift caused by recent bottlenecks (Coates et al., 2018). Finally, with regards to conservation, genomic data offer a snapshot of the genetic status of the populations by providing reliable and comparable estimates of genetic diversity and inbreeding, which are essential for genetic risk assessment. When genetic erosion is confirmed, it can be alleviated by facilitating gene flow with other conspecific populations. An appropriate source populations must be identified which, ideally, would restore genetic diversity while preserving possible local adaptations and minimizing the probability of reducing further the fitness through outbreeding depression, issues which can also be informed by genomic data (Ralls et al., 2018). In this study, we produced whole-genome resequencing data from representative samples of the Balkan and the North-Eastern portion of the Caucasian lynx populations (Figure 1.1) and analysed them in combination with available whole-genome data from across the species distributional range (Lucena-Perez et al., 2020). We aimed to provide additional evidence in the assessment of their taxonomic status as distinct subspecies and conservation units, and to evaluate their genetic status. We have done this by reconstructing their history and demography, assessing their genetic diversity and inbreeding levels, and analysing their relationship between each other and with the rest of Eurasian lynx populations genomically assessed so far. Our results add support to the distinction between Balkan and Caucasian lynx, confirm the anticipated depauperate genetic status of Balkan lynx and help to identify a suitable population to serve as a source for genetic rescue for these critically endangered felids. METHODS Sampling and DNA extraction We sampled blood from seven Lynx lynx individuals from the Balkan population and tissue from 12 Caucasian individuals. Samples from the Balkans cover the small extant range of the subspecies, including Macedonia, Montenegro, Kosovo and Albania, whereas samples from the Caucasus cover the NE portion of the Caucasus lynx distribution, including Georgia, Armenia and the Republics of Dagestan and Kabardino-Balkaria, Russian Federation (Figure 1.1). We also sampled tissue of one additional individual from each of the populations previously sampled in Lucena-Perez et al. (2020), Carpathians, Latvia, Norway, NE-Poland, Tuva, Urals and Yakutia, two from the Primorsky Krai population, although no new individuals were sequenced from the Mongolia population. Together with the individuals analysed in Lucena-Perez et al. (2020), these sum up to 124 L. lynx individuals which yielded mitogenomic sequences (114 individuals) and/or nuclear sequences (105 individuals; Table S1.1). We used data of a Lynx rufus (bobcat) individual from Jerez Zoo (Spain) as an outgroup in some of the analyses (see below). The license for lynx live-trapping and blood sampling in the Balkans as well as the license for tissue sampling in the rest of the 37 populations was obtained from the Macedonian Ministry of Environment and Physical Planning and the National Ethics Committee for Animal Experiments (11-2186/2; 11-546/2; 11-1006/10). This license adheres to the Directive 2010/63/EU on the protection of animals used for scientific purposes. No animals were harmed during live-trapping and handling. For DNA extractions, we performed different protocols. Most samples were processed by overnight digestion using proteinase K and extracted using silica-coated paramagnetic beads (NucleoMag® Tissue, MACHEREY-NAGEL GmbH & Co. KG). Samples yielding DNA concentration too low to be sequenced were re-extracted using a classical phenol-chloroform protocol. One Caucasian sample was a wet claw, and DNA extraction was performed using a protocol meant for bone (Rohland & Hofreiter, 2007) in a sterile lab. Sequence data generation and processing Library preparation was carried out in Centro Nacional de Análisis Genómico (CNAG) facilities. Briefly, gDNA was sheared, size selected, end-repaired and adenylated following the appropriate Illumina protocol. After indexed paired-end adapter ligation, quantity, quality and size of the libraries were assessed. Finally, libraries were sequenced using Illumina HiSeq2000 in CNAG facilities. Primary data analysis was carried out using the standard Illumina pipeline. We performed a quality control of our data using fastqc (https://www.bioinformatics.babraham.ac.uk/projects/fastqc), and sequencing reads were trimmed using the software Trimmomatic 0.38 (Bolger et al., 2014). Mitogenomic haplotypes for Balkan and Caucasian samples were reconstructed as described by Lucena-Perez et al. (2020). Briefly, random subsamples of ca. 5 M reads were mapped to the L. lynx mitochondrial reference genome generated by Abascal et al. (2016) using bwa-mem (Li, 2013) with default parameters, then SNPs were called using Freebayes 1.0.2-29-g41c1313 in haploid mode (minimum read mapping quality 30, minimum base quality 20, a minimum of two bases supporting an alternative allele call, and a maximum read mismatch fraction of 0.2) (Garrison & Marth, 2012), and consensus sequences were generated from called SNPs and the reference genome with FastaAlternateReferenceMaker command in gatk (McKenna et al., 2010), with sites not covered by any read or with low-quality genotypes coded as N. Repetitive regions RS2 and RS3 (positions 16096-16382 and 16908-174, respectively) (Sindičić et al., 2012), were excluded from the analysis. The calling procedure was done conjunctly with the previously analysed samples. In our previous mitogenome reconstructions, Multiple Nucleotide Polymorphisms (MNPs) and indels were filtered out but were not masked (Lucena-Perez et al., 2020), so some mitogenomes erroneously included the reference base at these positions. These sites are now called or left as N depending on coverage, so the sequences of the mitogenome haplotypes analysed here slightly differ from those previously reported (GenBank acc. Numbers: MK229199.2-MK229293.2 for previously reported now updated sequences; OK539600-OK539619, for novel sequences). 38 Additionally, and in order to expand the geographical coverage of the study, we re-analysed cytochrome b (cytb; 1140 pb; S = 7) and cytochrome oxidase I (COI; 630 pb; S = 5) sequences of Turkish lynx (İbİş et al., 2019), along with overlapping sequences in our mitogenome dataset. The resulting alignment of the concatenated sequences (cytb + COI 1770 bp; S = 19) contained our previously reported mitogenomic haplotypes, and the four haplotypes reported by İbİş et al. (2019). Sequencing reads were subsequently mapped to two nuclear reference genomes separately. Samples not previously analysed in Lucena-Perez et al. (2020), were mapped to the same 2.4 Gb Lynx pardinus reference genome LYPA1.0; (GCA_900661375.1; Abascal et al., 2016), using bwa-mem (Li, 2013), following the same workflow of Lucena-Perez et al. (2020). This was done in order to complete that dataset and compare the results for newly sequenced populations (i.e. the Balkans and the Caucasus), with previously published L. lynx genomes (Lucena-Perez et al., 2020). In order to take advantage of a more contiguous reference genome for some analyses (see below), all newly and previously sequenced whole-genome samples were also mapped to a 2.4 Gb Lynx canadensis nuclear reference genome (mLynCan4_v1.p; GCA_007474595.1; Rhie et al., 2020), using the same workflow. After alignment, duplicate reads were marked with the MarkDuplicates function of Picard Tools (http://broadinstitute.github.io/picard), and INDELs were realigned using the gatk (McKenna et al., 2010) commands RealignerTargetCreator and IndelRealigner. Average sequencing depth for each sample was calculated using SAMtools depth (Li et al., 2009), and can be found in Table S1.1. Variant calling was performed on the genome alignments of all the L. lynx samples to the L. canadensis reference genome using the gatk (4.1.4.1) HaplotypeCaller command. We obtained a first dataset of unfiltered variants consisting of 25.52 million SNPs. We then proceeded to filter out variants found in low mappability and repetitive regions, as well as all INDELs and non-biallelic variants. Low quality variants were also filtered out following gatk’s suggested thresholds of variant quality (i.e., QUAL < = 30, QD < = 2.0, FS > = 60.0, MQ < = 40.0) after inspecting the genome-wide distribution of values of these statistics. We maintained only the variants that were represented in at least four samples of each population and with a minimum depth of 3× per sample. Variants with a depth exceeding the average depth plus 1.5 times the standard deviation were also filtered out, in order to remove variants in regions where multiple paralogs are collapsed in the reference genome. Mean depth and its standard deviation were calculated using the software SAMtools depth. The final filtered VCF dataset consisted of ~5 million SNPs. We repeated the same steps with the single genome alignment of the L. rufus sample, in order to include it as an outgroup in the TreeMix analysis (see below). For Principal Component Analysis (PCA), NGSadmix, diversity analysis, genetic distances and FST analyses (see below), bam files were subsampled to a depth within the range of the lower depth data (5–6×) using SAMtools view -s (Li et al., 2009) to avoid biases associated with differences in sequence depth. Also, for these analyses, we only considered intergenic regions (61% of the nuclear genome, ~1.5 Gb) to account for neutral evolution. Sex 39 chromosomes were identified and excluded from diversity analysis as in Lucena-Perez et al. (2020). Mitogenomic phylogeography Consensus mitogenome sequences were aligned and collapsed into distinct haplotypes using the pegas R package (Paradis, 2010). Phylogenetic relationships among haplotypes were inferred by constructing a median-joining haplotype network (Bandelt et al., 1999) in PopART (Leigh & Bryant, 2015) and by inferring Maximum-likelihood phylogenetic trees with RAxML 8.2.11 (Stamatakis, 2014), as implemented in the Geneious Prime suite, under a GTR model with a gamma distribution of rates across sites and a proportion of invariant sites (GTRGAMMAI), with node support estimated with 100 nonparametric bootstrap replicates. The cytb + COI alignment was analysed in a similar way. We used beast 2.5.2 (Bouckaert et al., 2014) to construct a dated phylogeny of lynx mitogenomic haplotypes. We excluded the control region and ND6 gene, leaving a total of 14,923 sites with quite homogeneous base composition and evolutionary parameters that can be analysed as a single partition, and applied a HKY + G + I model, a strict clock calibrated at 1.54 × 10–8 mutations/site/year, and selected a Coalescent Bayesian Skyline tree prior (Ritchie et al., 2016). The MCMC chain was run for 10 M steps and was sampled every 1,000 steps, and the first 10% sampled steps were discarded as burn-in. All parameters showed an ESS >500. TreeAnnotator 2.4.8 was used to obtain the Maximum Clade Credibility tree, which was visualized with FigTree 1.4.3 (https://github.com/rambaut/figtree). Nuclear demographic and divergence reconstruction using PSMC Changes in effective population size (Ne) between 2 Mya and 10 Kya were estimated on the basis of the nuclear genome, using a Pairwise Sequentially Markovian Coalescent model (PSMC; Li & Durbin, 2011). This model infers population size history from the distribution of the local density of heterozygous sites in a single diploid sequence. For this analysis, we used autosomal whole genome data of 15 L. lynx samples sequenced at higher depth (>19×), and aligned to the L. canadensis reference genome. All of our sampled populations were represented with at least one individual, with the exception of the Mongolia population for which we did not have any high-depth sample. For each sample, we generated a diploid consensus sequence using SAMtools mpileup as suggested by Li and Durbin (2011). Minimum read depth was set to 5× and maximum read depth to 50× for all samples. Non-autosomal scaffolds were discarded and the PSMC analysis was run with the following flags: “-N25 -t15 -r5 –p” and "4 + 25 × 2 + 4 + 6”. The PSMC model can also be used to infer divergence time between populations by creating a diploid sample combining two-phased haplotypes from two samples of different populations into a pseudo-diploid sequence (Cahill et al., 2016; Chikhi et al., 2018). The point in time at which the estimated Ne for the pseudo-diploid sample rises to infinitely large values can be interpreted as the moment when gene flow ceases between the two populations. For our phased haplotypes, we used the X chromosome sequences of male samples. To 40 identify male and female samples, a custom script was run to calculate the average depth across the autosomes, and the X and Y chromosomes using SAMtools depth. Samples were considered to be males when the average X chromosome depth was around half of the autosomal average, and the average Y chromosome depth was larger than 1. We investigated split times for Balkan and Caucasian populations with the rest of the available populations. As we did not have a high depth male sample for all populations, we only analysed all possible pseudo-diploid combinations between one male from the Caucasian and Balkan populations, and a male from each of the NE-Poland, Carpathian and Yakutia populations, as well as between each other. X chromosome consensus sequences were extracted for each of these samples and the pseudo-diploids were created using seqtk mergefa (https://github.com/lh3/seqtk). PSMC analyses were then conducted on the pseudo-diploids as described above for true diploids. We performed 100 bootstraps (Li & Durbin, 2011) for both the original and the pseudo-diploids analyses. Outputs were plotted, and scaled to time and population sizes assuming a mean generation time of 5 years (Lucena-Perez et al., 2018) and a mutation rate per site per generation of 6 × 10−9 (Abascal et al., 2016) for the autosome. For the X chromosome, we calculated a X chromosome to autosome mutation rate ratio of 0.944, based on a male mutation rate bias of 1.4 (Sayres et al., 2011) and the contribution of two thirds of the total mutation rate in the X chromosome in females. This translates into a mutation rate for the X chromosome of 5.664 × 10−9 per site per generation. Nuclear divergence and admixture reconstruction using TreeMix To infer patterns of splits and admixture between L. lynx populations we used TreeMix 1.12 (Pickrell & Pritchard, 2012), using the SNPs data from the mapping to L. canadensis as input. Following Lucena-Perez et al. (2020), allele counts from autosomal intergenic regions were extracted from the L. lynx populations VCF file as well as from a second VCF file which included the one L. rufus sample to be used as outgroup; then both resulting allele count files were merged under the assumption that any SNP absent in one of the species (but present in the other) would be fixed for the reference allele. We ran TreeMix 1.12 (Pickrell & Pritchard, 2012) setting L. rufus as the outgroup and the block size to 100. We modelled between zero and six migration events (0 ≤ m ≤ 6) and calculated the proportion of variance in relatedness between populations explained by each model. To assess the consistency of migration edges, we performed 10 additional runs for m = 2 and another ten for m = 3 using different random seeds. All results were plotted using the R script included in TreeMix. Nuclear genomic structure To assess the genetic relationships among the two populations of interest and all previously analysed populations, we performed PCA and NGSadmix analyses with the available mappings to the L. pardinus reference together with the newly sequenced individuals, and using ANGSD as in (Lucena-Perez et al., 2020), that is – uniqueOnly 1 – remove_bads 1 – only_proper_pairs 1 – baq 1 – C 50 – minMapQ 30 –minQ 20 –doCounts 1 –minInd (number of individuals in the population/2) – setMaxDepth (average (AVR) depth for the population + 41 Nuclear genomic structure Caucasian and Balkans populations are homogeneous in terms of ancestry. The Evanno et al. (2005) method gave overwhelming support to K = 2, followed by K = 3 (Figure S1.6). Consistently with previous results from Lucena-Perez et al. (2020), the uppermost level of structure (K = 2) splits samples in an East-West divide, with the Balkan population showing 100% western ancestry, like all other European populations west of Urals (Urals included), and the Caucasian population showing mostly European but also some Asian ancestry (Figure 1.4a). The complete western ancestry in the Balkans contrasts with previous results showing partial Asian ancestry (Lucena-Perez et al., 2020), which we provisionally attribute to the larger and higher quality samples used in this study. Congruent patterns are recovered by PCA analyses, with Balkans aligning with other European populations and Caucasus occupying an intermediate position but closer to Western populations in PC1 (Figure 1.4b). At higher Ks, lynxes from the Caucasus and the Balkans separate as an additional cluster at K = 3, and as independent ones starting at K > 4 (Figure 1.4); they also occupy opposite positions along PC3 (Figure 1.4). Other populations are assigned entirely to either the Western or Eastern cluster at K = 2 and K = 3, with few exceptions. Moderate and low amounts of western ancestry are recovered both in the Tuva and Yakutia populations, respectively, as previously reported (Lucena-Perez et al., 2020). Carpathian lynx show about half of its ancestry in Western and half in the Balkan-Caucasian clusters at K = 3, and are placed between these two populations in PC2 and closer to the Balkan population in PC3 (Figure 1.4). Increasing number of clusters resulted in further splits of populations into separate clusters, with all populations separating into different clusters at K = 8, except for Kirov and Urals which group into one cluster, while Latvia, Tuva and Mongolia show mixed ancestry from different clusters. Groupings at K = 8 are consistent with those recovered at K = 6 in Lucena-Perez et al. (2020), where the Caucasian population was missing and the Balkan population was represented by only three low-quality samples, what might have reduced the ability to recover it as an independent cluster (Puechmaille, 2016). At K = 9–13, convergence was much weaker, with many different reconstructions and many of the runs producing implausible clusters splitting one or a few individuals of a single population (data not shown), similarly to what was observed in Lucena-Perez et al. (2020) for K > 6. 48 Figure 1.4 Relationship among individuals based on nuclear autosomal genotypes. (a) Individual ancestry inferred in the NGSadmix analysis assuming from two to eight genetic clusters (K2–K8). The result of the run with the highest likelihood is displayed for each K. Populations are sorted from west to east. Partial Asian ancestry is inferred for Caucasus at K2 but at K3 Caucasus groups with Balkan in a “Southern” cluster and Carpathians appear with mixed western and southern ancestry. Balkan and Caucasian lynx tend to form separate and exclusive clusters in runs at K ≥ 4. (b) PCA separates eastern and western individuals in the first axis (6.43% of the variance explained), and westernmost populations in the second axis (3.41% of the variance explained). Balkans aligns with the other European populations in the first axis, whereas Caucasian lynx are placed in between eastern and western populations but closer to western 49 Recent demographic history The software GONE infers generally declining numbers for all of the L. lynx populations during the last 750 years (Figure 1.5a). Most populations show relatively steady numbers followed by a sharp decline from around 50 to 70 generations ago, until around 15 to 25 generations ago. A less sharp but still declining trend is shown by Kirov, Latvia, Urals and Balkan populations. The populations with the highest recent effective population sizes are Kirov and Tuva, with a mean of around 550 in the last 10 generations (Figure 1.5b). Recent effective sizes range from ~150 to ~400 in the rest of the populations, except for bottlenecked European populations (Carpathians, NE-Poland, Norway and Balkans) which are all below 80. The Balkans remains constantly among the smallest populations throughout the reconstructed time period, and the smallest in the 10 most recent generations (Ne =30.6). Some populations (e.g. NE-Poland and Mongolia) show high variation in the effective population size reconstructed across repetitions of the analysis for less recent generations (>40), so effective population size estimations should be interpreted with caution (Figure S1.7). Figure 1.5 Changes in effective population size (Ne) in the last 150 generations (a), and in the last 10 only (b), for the different populations sampled, inferred using the software GONE. The mean value of 20 independent replicates is represented. See Figure S1.7 for variation among replicates 50 Nuclear and mitogenomic diversity The Balkan population shows the lowest nuclear θW diversity levels amongst all Eurasian lynx populations analysed so far, and its lowest θW X/A ratio suggests a contribution of recent bottlenecks (Figure 1.6a,b; Figure S1.8; Table S1.4). Contrastingly, the diversity of Caucasian lynx is within the range of Asian populations, and its θW X/A is closer to (although lower than) that expected under equilibrium (i.e. 0.75) (Pool & Nielsen, 2007), indicating larger and more stable recent population sizes (Figure 1.6a,b; Figure S1.8; Table S1.4). Patterns in π diversity are similar, with Caucasus showing values closer to those from Kirov and Urals. It must be noted that θW reported in Figure 1.6a for Caucasus (4.41 × 10–4) corresponds to a subset of nine individuals, because we observed extremely high θW diversity when considering all individuals (n = 12; 6.7 × 10–4). This was apparently due to three low-quality and possibly slightly contaminated DNA samples, which showed a high percentage of reads that did not pass the quality filter during trimming, a high percentage of reads with clipped data, and a high divergence compared with other Caucasian individuals based on individual genetic distances and PCA (Figure S1.9; samples c_ll_ca_0245, c_ll_ca_0248 and c_ll_ca_0254). The contrasting recent demography of the Balkan and Caucasian populations is further reflected in the site frequency spectrum (SFS), which is almost flat in the former, but shows an excess of low frequency alleles in the latter (Figure S1.10). In agreement with the autosomal data, mitochondrial diversity in the Balkan population is among the lowest, with only two haplotypes differing at a single site (Hd = 0.36 ± 0.16; π = 2.16 × 10–5 ± 2.58 × 10–5), while Caucasian population, with Hd = 0.78 ± 0.10 and π = 2.18 × 10–4 ± 1.33 × 10–4, is among the highest (Figure S1.11). Both Balkan and Caucasian lynx are the most differentiated when nucleotide distances are taken into account (Table S1.3). Runs of Homozygosity (ROH) and inbreeding The different populations show different levels of average population inbreeding coefficients (avg FROH) (Figure 1.6c). All the populations, except for Balkans, have an average FROH of less than 0.04, with Mongolia, Tuva, Latvia and the Carpathians being the least inbred, with an average FROH of less than 0.01. Some individuals appear as outliers within their population, showing levels of inbreeding much higher than the population average, with one individual from the Caucasian population (avg FROH = 0.025) showing strikingly high levels of inbreeding (FROH = 0.198). The most inbred population is the Balkan population by a large margin, with an average FROH of 0.158, although one individual shows very low inbreeding (FROH < 0.01). Cumulative ROH lengths follow a linear relationship with the total number of ROH blocks, an indication of relative uniformity of ROH block length across samples (Figure S1.12). The most notable exception is the single inbred Caucasian individual. This individual shows a level of inbreeding (i.e. cumulative ROH length) similar to Balkan individuals, but with fewer – and thus larger – ROHs, an indication of a recent inbreeding rather that small historical population sizes (Ceballos et al., 2018). 51 Figure 1.6 Genetic diversity and inbreeding coefficient in the different lynx populations. Populations are sorted from west to east. (a) Autosomal diversity is estimated as Watersson's theta (θW). (b) Ratio of X chromosome to autosomal Watersson's theta diversity. (c) Individual and population average inbreeding coefficients (FROH), represented by dots and lines, respectively, estimated as the cumulative length of autosomal ROHs of >2Mb divided by the total length of the autosome 52 DISCUSSION In our study, we report the first comprehensive genomic analyses of the Balkan and Caucasian populations of the Eurasian lynx, exploring both mitogenomic and nuclear genomic patterns. These analyses provide valuable insights into their history, their phylogenetic relationships with the rest of the Eurasian lynx populations and their current genetic status. The expanded mitogenomic data confirm the existence of a divergent mitogenomic lineage exclusive to Balkan lineage (HG 1) and confirm the existence of a divergent lineage exclusive to Caucasian lynx (HG 6; lineage B1 in Mengüllüoğlu et al. (2021)). These two lineages diverged from other lynx mitogenomic lineages around 92 kya and between each other around 46 kya, most likely due to isolation in separate glacial refugia (Figure S1.13). We did not observe the second dinniki sublineage that is more closely related to Balkan mitogenomes than to the more divergent Caucasian sublineage and was sampled by Mengüllüoğlu et al. (2021) in Southwestern Anatolia, a region that was not covered in this study. Past episodes of isolation reflected in deep mitogenomic divergences were, however, followed by extensive admixture driven by male-biased dispersal, causing homogenization of the nuclear genome and overall discordance of population relationships derived from mitogenomic and nuclear data (Lucena-Perez et al., 2020). Genetic relationships based on the nuclear genome are generally closest among neighbouring populations despite some pairs representing reciprocally monophyletic groups with deep mitogenomic divergences. The most striking case is that of the Carpathians and Balkans populations, which represent two haplogroups diverging around 100 Kya, but appearing as closely-related in Treemix, admixture and PCA analyses of nuclear data. Furthermore, it must be noted that intense genetic drift in recent times in these and other bottlenecked European populations (i.e. Norway and NE-Poland) tends to increase nuclear genetic divergence between populations, so the genetic relationship of Balkans and Carpathians might have been even closer in pre-bottleneck times. In contrast, the Caucasus shows a relatively large divergence with other western populations at both the mitochondrial and the nuclear genomes, which, along with the absence of signals of intense genetic drift, is suggestive of long-term isolation from other populations. This interpretation is further supported by the older dates of cessation of gene flow between Caucasus and other populations (around 20 and 30 kya), whereas the Balkan-Carpathian isolation is the most recent (around 15 kya). It must be noted, though, that a more extensive sampling is needed to integrate currently unexplored areas, delimit the distributional range of each haplogroup and to test the monophyly of the mitogenome of the different subspecies. Indeed, our reanalysis of short cytb + COI sequences indicate the occurrence of haplotypes external to HG 6 and HG 1 in Turkey, within the proposed range of L. l. dinniki (Figure 1.1), and Gugolz et al. (2008) reported a control region haplotype shared by Balkans and Caucasus. More recently, a mitogenome analyses revealed a mitogenomic sublineage related to the Balkan HG 1 in Anatolia, and a basal divergent mitogenomic haplotype occurring in one lynx from the Khingan Mountains, Northeast China, within the range of L. l. isabellinus (Mengüllüoğlu et al., 2021). 53 The results of our study add the relatively deep divergence of the Caucasian population to the East–West genomic divergence observed in the previous study (Lucena-Perez et al., 2020). These divergences started around 100 Kya but increased during a generalized post-glacial demographic decline that intensified during the Holocene (Lucena-Perez et al., 2020). These processes became more intense in western Europe in historical times, resulting in the extirpation of the species from most of the area and leaving a few small, isolated and genetically differentiated populations in the westernmost limit of the current range (Lucena-Perez et al., 2020). Our results confirm this generalized decline and show moderate to small historical effective sizes and drastic declines occurring in most eastern and western populations during the last several centuries. Climate cooling that occurred between 1300 and 1850 – the so called “Little Ice Age” had negative impacts on abundance of deer species sensitive to harsh winters and not specifically adapted to deep snow, including the roe deer, a main prey of the Eurasian lynx in many areas of its distribution (Ossi et al., 2015; Telfer & Kelsall, 1984). These likely negative effects of climate in prey species, in conjunction with the increase in direct and indirect human impacts, could be behind the inferred drastic lynx declines during the last few centuries. Current effective sizes are generally modest, but are alarmingly small in the bottlenecked western European populations, including the Carpathians, NE-Poland, Norway and, especially, the Balkans. This generalized declining trend highlights the need for the implementation or revision of active conservation plans. Where conservation plans are already in place, these must be revised following our results and eventually intensified. Signs of recent populations increases are observed for Latvia and Tuva, but not for Scandinavia or the Carpathian Mts., despite reports suggesting recent recovery (Linnell et al., 2009); these latter recoveries may be too recent and/or weak to be recorded as increases in effective size with our methodology (Santiago et al., 2020). Population levels of inbreeding do not fully correlate with what might be expected looking at population level genomic diversity. Moderate inbreeding is found for relatively highly diverse populations, such as Kirov, Urals, Yakutia and Caucasus, although the average of these last three populations might be influenced by a few highly inbred individuals. On the other hand, relatively low levels of inbreeding are found in some of the populations that show low genomic diversity, such as Carpathians, Poland and to some extent Norway, which has many non-inbred individuals. This suggests that in these cases, population reductions have not been drastic or lasting enough, and have not been accompanied by levels of fragmentation high enough, to result in frequent mating between close relatives. The population with the highest inbreeding coefficient, together with the lowest genomic diversity, is the Balkans, indicating a concerning level of genetic erosion. Taxonomic implications Our genomic analyses of lynxes representing five of the six currently recognized subspecies (Kitchener et al., 2017) provide substantial evidence relevant for the taxonomic status of L. l. balcanicus and L. l. dinniki, and for addressing the possibility of L. l. balcanicus being part of 54 L. l. dinniki (von Arx et al., 2004; Kitchener et al., 2017). The combination of deep mitogenomic and autosomal divergence in the absence of signals of genetic drift supports the long-term genetic distinctiveness of the Caucasian lynx and its taxonomic identification as a separate subspecies, L. l. dinniki. Although sharing a common ancestor that diverged ~90 kya from the rest of Eurasian lynx, Balkan and Caucasian lynx mitogenomic lineages diverged ~40 kya, and they currently share no mitogenomic haplotypes, at least in our set of samples. The observation by Mengüllüoğlu et al. (2021) of a second mitogenomic sublineage in Southwestern Anatolia, more closely related to Balkan mitogenomes than to the other Caucasian sublineage, makes dinniki paraphyletic with respect to balcanicus, suggesting a more recent mitogenomic divergence and supporting the inclusion of the Balkan lynx as part as dinniki. However, the two subspecies appear to have had different nuclear genetic histories, with gene flow between the two becoming negligible around the Last Glacial Maximum, decreasing support to these two subspecies being synonymous. The Balkan lynx mitogenomic distinctness appears to have been countered by intense interglacial male-biased gene-flow with neighbouring European populations, especially the Carpathians, and a large part of the nuclear divergence observed now may have been caused by intense genetic drift, driven by low effective population size and increasing isolation during the Holocene. Therefore, the status of L. l. balcanicus as a distinct subspecies is still debatable. With respect to the other subspecies, the status of L. l. carpathicus may also be questioned, based on their lack of mitogenomic divergence, its single haplotype being part of HG 2 distributed from Poland, across the Baltic states, to Urals, and its nuclear relationships with Balkans, possibly as the consequence of admixture, as suggested by clustering analyses (K = 3–5; Figure 1.4). Our updated genomic survey reinforces the support for two widely distributed subspecies: L. l. lynx (Europe) and L. l. wrangeli (Asia), in line with Kitchener et al. (2017), with a recent suggestion that two samples identified as L. l. isabellinus in Lucena-Perez et al. (2020), should instead be ascribed to L. l. wrangeli based on mitogenomic data (Mengüllüoğlu et al., 2021). In any case, a genome-wide analyses of samples across southern Asia is necessary to clarify the status of L. l. isabellinus and to obtain a more complete picture of intraspecific L. lynx taxonomy. Conservation implications The Balkan lynx is the most genetically eroded and inbred of the Eurasian lynx populations analysed so far, which alerts of a high risk of extinction due to inbreeding depression and reduced adaptive potential. Even though the genomic signatures of possible local adaptations were not assessed in this study, this possibility may be considered unlikely given modest environmental differentiation, extensive gene flow in the past and relatively short isolation times; furthermore, any possible past adaptation may have already been lost, while maladaptive alleles may have accumulated due to intense genetic drift. We therefore encourage the consideration of genetic rescue through assisted gene flow, as the risks of possible outbreeding depression are much lower than the benefits brought by the entrance of new alleles into the population (Ralls et al., 2018). Its closest genomic relatedness to Carpathian lynx identifies this population as an optimal source of animals for genetic 55 restoration, pending the genomic analyses of Anatolian lynxes, a part of the distribution that was not covered by our analysis. Using males would be advisable to mimic natural male-biased gene flow in this large carnivore (Samelius et al., 2012) and to preserve the uniqueness of the Balkan lynx mitochondrial genome. In conclusion, we have shown the potential of using genomic data with the aim of addressing both intraspecific taxonomy and conservation questions regarding wild non-model species. Patterns of genomic and mitogenomic differentiation, and the reconstruction of past and recent demography provided necessary lines of evidence for the delineation of taxonomic and conservation units, allowed the identification of populations in need of genetic restoration, and identified most suitable sources. Taxonomic and conservation issues are not independent, though, as conservation planning for many endangered species is often guided by subspecies taxonomy. Additional input into these issues will however require the explicit analysis of adaptive variation and the investigation of possible adaptive divergences. REFERENCES Abascal, F., Corvelo, A., Cruz, F., Villanueva-Cañas, J. L., Vlasova, A., Marcet-Houben, M., Martínez-Cruz, B., Cheng, J. Y., Prieto, P., Quesada, V., Quilez, J., Li, G., García, F., Rubio-Camarillo, M., Frias, L., Ribeca, P., Capella-Gutiérrez, S., Rodríguez, J. M., Câmara, F., … Godoy, J. A. (2016). Extreme genomic erosion after recurrent demographic bottlenecks in the highly endangered Iberian lynx. Genome Biology, 17(1), 251. https://doi.org/10.1186/s13059-016-1090-1 Askerov, E., Karen, M., Gurielidze, Z., Mousavi, M., Shmunk, V., Trepet, S., Kutukcu, A. E., Heidelberg, A., & Zazanashvili, N. (2020). Status of large carnivores in the Caucasus. In Ecoregional Conservation Plan (ECP) For The Caucasus 2020 Edition (Supplementary Reports, pp. 37–47). WWF, KfW. Bandelt, H. J., Forster, P., & Rohl, A. (1999). Median-joining networks for inferring intraspecific phylogenies. Molecular Biology and Evolution, 16(1), 37–48. https://doi.org/10.1093/oxfordjournals.molbev.a026036 Bolger, A. M., Lohse, M., & Usadel, B. (2014). Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics, 30(15), 2114–2120. https://doi.org/10.1093/bioinformatics/btu170 Bouckaert, R., Heled, J., Kühnert, D., Vaughan, T., Wu, C.-H., Xie, D., Suchard, M. A., Rambaut, A., & Drummond, A. J. (2014). BEAST 2: A Software Platform for Bayesian Evolutionary Analysis. PLoS Computational Biology, 10(4), e1003537. https://doi.org/10.1371/journal.pcbi.1003537 Breitenmoser, U., Breitenmoser-Würsten, C., Lanz, T., von Arx, M., & Antonevich, A. (2015). Lynx lynx (errata version published in 2017). IUCN Red List of Threatened Species. https://www.iucnredlist.org/en Cahill, J. A., Soares, A. E. R., Green, R. E., & Shapiro, B. (2016). Inferring species divergence times using pairwise sequential Markovian coalescent modelling and low-coverage genomic data. Philosophical Transactions of the Royal Society B: 56 Biological Sciences, 371(1699), 20150138. https://doi.org/10.1098/rstb.2015.0138 Canty, A., & Ripley, B. D. (2017). boot: Bootstrap R (S‐Plus) Functions (Version 1.3) [Computer software]. https://cran.r-project.org/web/packages/boot Ceballos, F. C., Joshi, P. K., Clark, D. W., Ramsay, M., & Wilson, J. F. (2018). Runs of homozygosity: Windows into population history and trait architecture. Nature Reviews Genetics, 19(4), 220–234. https://doi.org/10.1038/nrg.2017.109 Chapron, G., Kaczensky, P., Linnell, J. D. C., von Arx, M., Huber, D., Andrén, H., López-Bao, J. V., Adamec, M., Álvares, F., Anders, O., Balčiauskas, L., Balys, V., Bedő, P., Bego, F., Blanco, J. C., Breitenmoser, U., Brøseth, H., Bufka, L., Bunikyte, R., … Boitani, L. (2014). Recovery of large carnivores in Europe’s modern human-dominated landscapes. Science, 346(6216), 1517–1519. https://doi.org/10.1126/science.1257553 Chikhi, L., Rodríguez, W., Grusea, S., Santos, P., Boitard, S., & Mazet, O. (2018). The IICR (inverse instantaneous coalescence rate) as a summary of genomic diversity: Insights into demographic inference and model choice. Heredity, 120(1), 13–24. https://doi.org/10.1038/s41437-017-0005-6 Coates, D. J., Byrne, M., & Moritz, C. (2018). Genetic Diversity and Conservation Units: Dealing With the Species-Population Continuum in the Age of Genomics. Frontiers in Ecology and Evolution, 6, 165. https://doi.org/10.3389/fevo.2018.00165 Evanno, G., Regnaut, S., & Goudet, J. (2005). Detecting the number of clusters of individuals using the software structure: A simulation study. Molecular Ecology, 14(8), 2611–2620. https://doi.org/10.1111/j.1365-294X.2005.02553.x Fumagalli, M., Vieira, F. G., Linderoth, T., & Nielsen, R. (2014). ngsTools: Methods for population genetics analyses from next-generation sequencing data. Bioinformatics, 30(10), 1486–1487. https://doi.org/10.1093/bioinformatics/btu041 Garrison, E., & Marth, G. (2012). Haplotype-based variant detection from short-read sequencing. arXiv:1207.3907 [q-Bio]. http://arxiv.org/abs/1207.3907 Gugolz, D., Bernasconi, M. V., Breitenmoser‐Würsten, C., & Wandeler, P. (2008). Historical DNA reveals the phylogenetic position of the extinct Alpine lynx. Journal of Zoology, 275(2), 201–208. https://doi.org/10.1111/j.1469-7998.2008.00428.x Hemmer, H. (1993). Felis (Lynx) lynx Linnaeus, 1758-Luchs, Nordluchs. Handbuch Der Säugetiere Europas, 5, 1119–1167. İbİş, O., Özcan, S., Kırmanoğlu, C., Keten, A., & Tez, C. (2019). Genetic Analysis of Turkish lynx (Lynx lynx) Based on Mitochondrial DNA Sequences. Russian Journal of Genetics, 55(11), 1426–1437. https://doi.org/10.1134/S1022795419110061 Kitchener, A. C., Breitenmoser-Würsten, C., Eizirik, E., Gentry, A., Werdelin, L., Wilting, A., Yamaguchi, N., Abramov, A. V., Christiansen, P., Driscoll, C., Duckworth, J. W., Johnson, W. E., Luo, S. J., Meijaard, E., O’Donoghue, P., Sanderson, J., Seymour, K., Bruford, M., Groves, C., … Tobe, S. (2017). A revised taxonomy of the Felidae: The final report of the Cat Classification Task Force of the IUCN Cat Specialist Group. Cat News Special Issue, 11, 80. 57 c_ll_ur_0200 Urals Russia 59.00 55.00 m y/y lucena2020 13.33 H_07 c_ll_ur_0202 Urals Russia 60.89 54.38 f y/n chapter 1 23.47 n.a. c_ll_ur_0203 Urals Russia 57.36 55.00 f y/y lucena2020 13.79 H_11 c_ll_vl_0107 Primosky Krai Russia 136.47 44.94 f y/y lucena2020 6.54 H_20 c_ll_vl_0108 Primosky Krai Russia 135.58 45.29 f y/y lucena2020 11.58 H_20 c_ll_vl_0109 Primosky Krai Russia 136.73 45.67 m y/y lucena2020 5.94 H_05 c_ll_vl_0110 Primosky Krai Russia 132.91 49.01 m y/y lucena2020 7.08 H_04 c_ll_vl_0112 Primosky Krai Russia 137.01 45.86 f y/y lucena2020 29.33 H_05 c_ll_vl_0113 Primosky Krai Russia 137.01 45.86 m y/y lucena2020 8.75 H_05 c_ll_vl_0114 Primosky Krai Russia 137.61 47.29 f y/n this study 23.4 n.a. c_ll_vl_0128 Primosky Krai Russia 137.43 45.91 m y/y lucena2020 8.29 H_20 c_ll_vl_0132 Primosky Krai Russia 137.43 45.91 m y/y lucena2020 8.66 H_20 c_ll_vl_0137 Primosky Krai Russia 133.46 51.39 f y/n chapter 1 25.47 n.a. c_ll_ya_0138 Yakutia Russia 132.07 59.90 m y/y lucena2020 8.21 H_21 c_ll_ya_0139 Yakutia Russia 129.20 61.19 m y/y lucena2020 8.28 H_22 c_ll_ya_0140 Yakutia Russia 129.16 61.16 m y/y lucena2020 8.48 H_23 c_ll_ya_0141 Yakutia Russia 136.18 66.89 f y/n chapter 1 23.6 n.a. c_ll_ya_0142 Yakutia Russia 136.18 66.89 m y/y lucena2020 8.63 H_07 c_ll_ya_0143 Yakutia Russia 136.18 66.89 m y/y lucena2020 8.14 H_07 c_ll_ya_0145 Yakutia Russia 129.45 61.89 m y/y lucena2020 8.59 H_24 c_ll_ya_0146 Yakutia Russia 130.68 60.78 m y/y lucena2020 23.58 H_23 c_ll_ya_0147 Yakutia Russia 127.30 60.75 f y/y lucena2020 8.26 H_23 64 Table S1.2. Distribution of mitochondrial haplotypes in populations. The number indicates the number of observations of each haplotype in each population. 65 Table S1.3. Mitogenomic diversity in Eurasian lynx populations. N, samples size; # Hap, number of haplotypes; S, number of segregating sites; Hd, haplotype diversity; Hd SD, standard deviation of Hd; π, nucleotide diversity; π SD, standard deviation of π; K, average number of nucleotide differences; D and F are neutrality indices of Tajima, (1989) and Fu and Li (1993), respectively; Fst (nuc), population Fst based on nucleotide distances. 66 Population N # Hap S Hd Hd sd Pi Pi sd K D D p-valu e F Fst (nuc) Balkans 10 2 1 0.36 0.16 2.16E-05 2.58E-05 0.36 0.01 0.99 0.68 0.95 Caucasus 13 6 11 0.78 0.10 2.18E-04 1.33E-04 3.97 0.46 0.65 0.78 0.92 Carpathians 6 1 0 0.00 0.00 0.00E+00 0.00E+00 0.00 NA NA NA 0.91 Kirov 14 5 23 0.70 0.11 5.39E-04 2.95E-04 9.73 1.28 0.20 1.61 0.69 Latvia 6 3 9 0.73 0.16 2.80E-04 1.84E-04 4.60 0.99 0.32 1.36 0.83 Mongolia 8 5 6 0.82 0.10 1.17E-04 8.45E-05 2.61 1.09 0.28 0.95 0.86 NE-Poland 12 2 7 0.15 0.13 6.57E-05 5.19E-05 1.17 -1.98 0.05 -2.3 2 0.89 Norway 17 2 1 0.12 0.10 7.15E-06 1.33E-05 0.12 -1.16 0.24 -1.4 8 0.91 Primorski Krai 8 3 12 0.68 0.12 3.41E-04 2.08E-04 6.14 1.61 0.11 1.30 0.79 Tuva 6 2 2 0.33 0.22 2.03E-05 2.67E-05 0.67 -0.93 0.35 -0.8 6 0.89 Urals 6 3 23 0.60 0.22 4.46E-04 2.81E-04 7.67 -1.50 0.13 -1.2 1 0.76 Yakutia 8 5 19 0.86 0.11 4.62E-04 2.75E-04 7.61 0.20 0.84 0.55 0.76 Table S1.4. Autosomal diversity and X/A diversity ratios in Eurasian lynx populations. N, samples size; θW, Watterson’s theta estimator; π, nucleotide diversity; X/A, ratio of chromosome X to autosome diversity; SE, standard error; SD, standard deviation. 67 Population N θw θw SE π π SE θw X/A θw X/A SD π X/A π X/A SD Balkans 7 1.82E-04 1.22E-06 2.19E-04 1.68E-06 0.521 0.022 0.424 0.022 Caucasus 7 4.12E-04 1.52E-06 4.31E-04 1.99E-06 0.685 0.013 0.549 0.014 Caucasus 9 4.41E-04 1.52E-06 4.49E-04 1.99E-06 0.727 0.013 0.559 0.014 Caucasus 12 6.70E-04 1.51E-06 5.19E-04 2.02E-06 0.880 0.008 0.636 0.012 Carpathians 6 3.14E-04 1.62E-06 3.72E-04 2.15E-06 0.551 0.019 0.509 0.021 Kirov 13 3.57E-04 1.58E-06 4.62E-04 2.27E-06 0.599 0.017 0.536 0.018 Latvia 6 3.72E-04 1.75E-06 4.36E-04 2.17E-06 0.600 0.017 0.536 0.016 Mongolia 8 4.62E-04 1.80E-06 5.34E-04 2.40E-06 0.683 0.013 0.588 0.014 Norway 8 2.65E-04 1.45E-06 3.39E-04 2.06E-06 0.615 0.020 0.520 0.020 NE-Poland 8 2.85E-04 1.55E-06 3.65E-04 2.12E-06 0.530 0.018 0.457 0.016 Primorski Krai 8 4.24E-04 1.73E-06 4.99E-04 2.27E-06 0.601 0.015 0.523 0.015 Tuva 6 4.87E-04 2.03E-06 5.43E-04 2.48E-06 0.658 0.015 0.593 0.015 Urals 6 3.93E-04 1.94E-06 4.50E-04 2.35E-06 0.607 0.017 0.520 0.019 Yakutia 8 4.49E-04 1.75E-06 5.25E-04 2.29E-06 0.677 0.016 0.547 0.015 Figure S1.1 Phylogenetic tree of Eurasian lynx mitogenomic haplotypes. The tree was inferred with RAxML 8.2.11, as implemented in the Geneious Prime suite, under a GTR model with a gamma distribution of rates across sites and a proportion of invariant sites (GTRGAMMAI). Numbers above branches indicate node support estimated with 100 non-parametric bootstrap replicates. Haplogroups are numbered following Lucena‐ Perez et al. (2020) and adding the novel HG 6 exclusive to the Caucasus. 68 Figure S1.2 Haplotype network of Eurasian lynx mitogenome haplotypes. Phylogenetic relationships among haplotypes were inferred by constructing a medianjoining haplotype network in PopART. Pie charts at nodes represent the haplotype frequency in different populations. Solid back dots represent not-observed inferred ancestral haplotypes. Bars on branches connecting the haplotypes indicate mutations. 69 Figure S1.3 Phylogenetic tree of concatenated cytb and COI sequences. The alignment includes the Eurasian lynx mitogenome haplotypes and the four Turkish haplotypes previously reported by (İbİş et al., 2019). The limited resolution of the analyzed fragment reduces the number of observed haplotypes to 11, grouped in only three supported clades: Balkan, Caucasian and the rest, with most non-southern samples sharing a single widespread haplotype. Three of the four sequences from Turkey grouped with Caucasian samples, whereas a fourth sequence belonged to the widespread non-southern haplotype. 70 Figure S1.4. Additional population trees inferred with TreeMix. (a) Tree with no migration events (m=0). (b) Most frequently observed tree (9/10 replicates) at m=2. 71 Figure S1.5 PSMC analysis bootstrap replicates. Main run (thick line) and bootstrap replicates (thin lines) of PSMC analysis are represented for the Caucasian individual (a) and the Balkan individual (b), for autosomal and all pseudodiplod pairs. Estimations are made assuming a per site per generation mutation rate of 6 × 10−9 (Abascal et al., 2016) for the autosome, a rescaled per site per generation mutation rate for the X-chromosome of 5.66 × 10−9, and a 5-year generation time (Lucena-Perez et al., 2018). 72 Figure S1.6 Delta K plot of Evanno’s test derived from CLUMPAK. Results for all Ks are displayed on the right. Large support for K=3 compared to K≥4 is more evident with a close up excluding K=2 (left panel). 73 Figure S1.13. Variation in temperature (oC) recorded in the Vostok ice core over the last 160,000 years (Petit et al. 1999; EPICA 2004) and average date of main mtDNA divergence events and their 95% confidence intervals. Indeed, the 95% interval limits of the estimated divergence times overlaps the last glacial period (LGP). The most ancient divergence, coincides with an earlier milder period within LGP, while the latter interval coincides with a subsequent harsher period preceding the last glacial maximum (LGM) (Figure S1.8). Assuming that temperate and boreal wooded areas comprise the optimal habitat for the Eurasian lynx (Breitenmoser et al., 2015; Breitenmoser & Breitenmoser-Würsten, 2008; Guggisberg, 1975; Matyushkin & Vaisfeld, 2003; Nowell et al., 1996; Schmidt et al., 2011), and based on the LGM biome model of Gavashelishvili & Tarkhnishvili (2016) and global biome distributions reconstructed over the last 120 kyr by Hoogakker et al. (2016) , the first milder event of the LGP fragmented the optimal habitat more intensely in northern latitudes, by replacing boreal forests (i.e. taiga) with steppe, desert, tundra and glaciers. The second harsher event of the LGP was associated to an extension of boreal forest in northern latitudes and increased fragmentation in Southern latitudes, likely impeding gene flow between Caucasian and Balkan populations and also isolating the Carpathian-Baltic populations from the remaining clades. This second cold period coincides well with the first divergence event (~47 kya) of the haplotype 3a1 of brown bears, when individuals from present-day Romania and Bulgaria became isolated from other populations of the species (Anijalg et al., 2018). 80 Abascal, F., Corvelo, A., Cruz, F., Villanueva-Cañas, J. L., Vlasova, A., MarcetHouben, M., Martínez-Cruz, B., Cheng, J. Y., Prieto, P., Quesada, V., Quilez, J., Li, G., García, F., Rubio-Camarillo, M., Frias, L., Ribeca, P., Capella-Gutiérrez, S., Rodríguez, J. M., Câmara, F., ... Godoy, J. A. (2016). Extreme genomic erosion after recurrent demographic bottlenecks in the highly endangered Iberian lynx. Genome Biology, 17(1), 251. Anijalg, P., Ho, S. Y. W., Davison, J., Keis, M., Tammeleht, E., Bobowik, K., Tumanov, I. L., Saveljev, A. P., Lyapunova, E. A., Vorobiev, A. A., Markov, N. I., Kryukov, A. P., Kojola, I., Swenson, J. E., Hagen, S. B., Eiken, H. G., Paule, L., & Saarma, U. (2018). Large-scale migrations of brown bears in Eurasia and to North America during the Late Pleistocene. Journal of Biogeography, 45(2), 394–405. https://doi.org/10.1111/jbi.13126 Breitenmoser, U., & Breitenmoser-Würsten, C. (2008). Der Luchs – Ein Grossraubtier in der Kulturlandschaft. Salm Verlag. Breitenmoser, U., Breitenmoser-Würsten, C., Lanz, T., von Arx, M., & Antonevich, A. (2015). Lynx lynx (errata version published in 2017). IUCN Red List of Threatened Species. https://www.iucnredlist.org/en Fu, Y. X., & Li, W. H. (1993). Statistical tests of neutrality of mutations. Genetics, 133(3), 693–709. https://doi.org/10.1093/genetics/133.3.693 Gavashelishvili, A., & Tarkhnishvili, D. (2016). Biomes and human distribution during the last ice age. Global Ecology and Biogeography, 25(5), 563–574. https://doi.org/10.1111/geb.12437 Guggisberg, C. A. W. (1975). Wild cats of the world. In Wild cats of the world (pp. 49– 58). David & Charles (Holdings) Limited. Hoogakker, B. A. A., Smith, R. S., Singarayer, J. S., Marchant, R., Prentice, I. C., Allen, J. R. M., Anderson, R. S., Bhagwat, S. A., Behling, H., Borisova, O., Bush, M., Correa-Metrio, A., de Vernal, A., Finch, J. M., Fréchette, B., Lozano-Garcia, S., Gosling, W. D., Granoszewski, W., Grimm, E. C., ... Tzedakis, C. (2016). Terrestrial biosphere changes over the last 120 kyr. Climate of the Past, 12(1), 51–73. https://doi.org/10.5194/cp-12-51-2016 İbİş, O., Özcan, S., Kırmanoğlu, C., Keten, A., & Tez, C. (2019). Genetic Analysis of Turkish lynx (Lynx lynx) Based on Mitochondrial DNA Sequences. Russian Journal of Genetics, 55(11), 1426–1437. https://doi.org/10.1134/S1022795419110061 Lucena-Perez, M., Marmesat, E., Kleinman-Ruiz, D., Martínez-Cruz, B., Węcek, K., Saveljev, A. P., Seryodkin, I. V., Okhlopkov, I., Dvornikov, M. G., Ozolins, J., Galsandorj, N., Paunovic, M., Ratkiewicz, M., Schmidt, K., & Godoy, J. A. (2020). Genomic patterns in the widespread Eurasian lynx shaped by Late Quaternary climatic fluctuations and anthropogenic impacts. Molecular Ecology, 29(4), 812–828. https://doi.org/10.1111/mec.15366 Lucena-Perez, M., Soriano, L., López-Bao, J. V., Marmesat, E., Fernández, L., Palomares, F., & Godoy, J. A. (2018). Reproductive biology and genealogy in the endangered Iberian lynx: Implications for conservation. Mammalian Biology, 89, 7–13. https://doi.org/10.1016/j.mambio.2017.11.006 81 Matyushkin, E. N., & Vaisfeld, M. A. (Eds.). (2003). The lynx – regional features of ecology, use and protection. Nauka. Nowell, K., Jackson, P., & IUCN/SSC Cat Specialist Group (Eds.). (1996). Wild cats: Status survey and conservation action plan. IUCN. Schmidt, K., Ratkiewicz, M., & Konopiński, M. K. (2011). The importance of genetic variability and population differentiation in the Eurasian lynx Lynx lynx for conservation, in the context of habitat and climate change. Mammal Review, 41(2), 112–124. https://doi.org/10.1111/j.1365-2907.2010.00180.x Tajima, F. (1989). Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics, 123(3), 585. 82 CHAPTER 2: Genome-environment association analyses reveal geographically restricted adaptive divergence across the range of the widespread Eurasian carnivore Lynx lynx (Linnaeus, 1758) Enrico Bazzicalupo, Mirosław Ratkiewicz, Ivan V. Seryodkin, Innokentiy Okhlopkov, Naranbaatar Galsandorj, Yuriy A. Yarovenko, Janis Ozolins, Alexander P. Saveljev, Dime Melovski, Alexander Gavashelishvili, Krzysztof Schmidt, José A. Godoy 83 ABSTRACT Local adaptations to the environment are an important aspect of the diversity of a species and their discovery, description and quantification has important implications for the fields of taxonomy, evolutionary and conservation biology. In this study, we scan genomes from several populations across the distributional range of the Eurasian lynx, with the objective of finding genomic windows under positive selection which may underlie local adaptations to different environments. A total of 394 genomic windows are found to be associated with local environmental conditions, and they are enriched for genes involved in metabolism, behaviour, synaptic organization and neural development. Adaptive genetic structure, reconstructed from SNPs in candidate windows, is considerably different from the neutral genetic structure of the species. A widespread adaptively homogeneous group is recovered occupying areas of harsher snow and temperature climatic conditions in the north-western, central and eastern parts of the distribution. Adaptively divergent populations are recovered in the westernmost part of the range, especially within the Baltic population, but also predicted for different patches in the western and southern part of the range, associated with different snow and temperature regimes. Adaptive differentiation driven by climate does not correlate much with the subspecies taxonomic delimitations, suggesting that subspecific divergences are mostly driven by neutral processes of genetic drift and gene flow. Our results will aid the selection of source populations for assisted gene flow or genetic rescue programs by identifying what climatic patterns to look for as predictors of pre-adaptation of individuals. Particularly, the Carpathian population is confirmed as the best source of individuals for the genetic rescue of the endangered, isolated and genetically eroded Balkan population. Additionally, reintroductions in central and western Europe, currently based mostly on Carpathian lynxes, could consider the Baltic population as an additional source to increase adaptive variation and likely improve adaptation to their milder climate. INTRODUCTION Local adaptations are expected to shape species diversity and contribute to intraspecific divergence and may eventually lead to speciation (de Queiroz, 2007). The identification of locally adapted populations should thus play a key role in the delimitation of evolutionary significant units (ESUs) for species conservation (Funk et al., 2012; Ryder, 1986), designation of ecotypes (Stronen et al., 2022) and subspecies (Schmidt et al., 2019). While there are still some aspects of the significance of adaptive genetic variation that remain not fully grasped, it is widely recognized that it plays an important role in speciation pathways and the global distribution of organisms (Doebeli & Dieckmann, 2004; Weissing et al., 2011). Understanding and describing species-wide patterns of functional genetic diversity is therefore key to prioritizing conservation efforts and guiding taxonomic delimitation (Schmidt et al., 2019; Stronen et al., 2022; Teixeira & Huber, 2021). One approach to predict and quantify adaptive population structure is trying to locate the areas of the genome which have been subjected to differential selective pressures among populations. Different approaches have been developed to achieve this goal (Hoban et al., 84 2016). Many studies have focused on the identification of genomic regions maximally differentiated among populations through summary statistics of genetic differentiation measured along the genome (Duforet-Frebourg et al., 2014; Gautier, 2015). One limitation of this type of inference is the inability to identify the underlying environmental or ecological factors driving the observed differentiation, beyond postulations based on the function of the genes involved. Other approaches have tried to work around this limitation by assessing directly the correlation between potential predictors (e.g. environmental variables) and genetic variation. These methods are commonly known as genome-environment association (GEA) analyses, and may be based on univariate and multivariate approaches, or a combination of both, which provides a higher power of detection (Forester et al., 2018). Loci identified as associated to local adaptation by any method will likely show a pattern of genetic structure different to that observed at neutral loci. Whereas the former will reflect the gradients of adaptive pressures (Nosil et al., 2008; Wang & Bradburd, 2014), the latter will be determined by the interaction of genetic drift and gene flow across space and time, and thus is usually correlated with geographic distance or past isolation in glacial refugia (Sobel et al., 2010; Wright, 1943). Although long-range dispersal was previously believed to counteract the population structuring, processes like genotype driven dispersal (Bolnick & Otto, 2013) or isolation by environment (Wang & Bradburd, 2014) may contribute to population differentiation and local adaptations. Large carnivores, in particular, often show genetic differentiation between geographically close populations, despite their ability to move long distances within their individual life-span. For example, climate and habitat variability were found to affect both genetic and phenotypic dissimilarities among grey wolf populations (Geffen et al., 2004; Musiani et al., 2007). Schweizer et al. (2016) found correlations of functional gene variants with temperature, precipitation, and vegetation variables among wolf ecotypes in genes related to morphology, vision, metabolism, and thermoregulation. Similarly, population structure of Canadian lynx has been associated with differences in snow conditions across its distribution (Stenseth, Shabbar, et al., 2004). The Eurasian lynx (Lynx lynx Linnaeus, 1758) is a large carnivore playing an essential role as a top-chain predator which inhabits most of the Eurasian continent, with populations having contrasting past and recent demographic histories. Although according to the IUCN the global status of this felid is Least Concern and most populations are stable, the status and trends vary greatly within its European range, with some populations being classified as endangered, while the status of Asian populations remain poorly known (Breitenmoser et al., 2015). Large and relatively diverse populations inhabit the Russian part of north-eastern Europe and most of the Asian continent, while highly fragmented populations that suffered recent demographic bottlenecks are found in Scandinavia and the Baltic states (Bazzicalupo et al., 2022; Lucena-Perez et al., 2020). Two other autochthonous isolated populations inhabit the Carpathian and Balkan regions (Bazzicalupo et al., 2022; Lucena-Perez et al., 2020), whereas few small reintroduced populations have been established in several areas in Central Europe using the Carpathians as their main source population, most of them exhibiting genome-wide diversity loss (Mueller et al., 2022). Throughout its wide distribution, the 85 Eurasian lynx occurs in a variety of ecological and climatic conditions. They inhabit a variety of climates, with extreme differences of minimum temperatures, levels of precipitation and snow cover, and habitats including mostly evergreen boreal forest, but also deciduous woodlands, alpine zones, grasslands, xeric shrublands and semi-deserts. The species also displays a broad phenotypic variability affecting size, morphology and pelage patterns (Darul et al., 2022; Matyushkin & Vaisfeld, 2003; Schmidt et al., 2011). Based mostly on this variability, a number of distinct Eurasian lynx subspecies have been described in different parts of the species range: L. l. lynx in central and eastern Europe, L. l. balcanicus in the south-western Balkans, L. l. carpathicus in the Carpathian mountains, L. l. isabellinus in central Asia, L. l. wrangeli in central and eastern Asia and L. l. dinniki in the Caucasian mountains and Anatolia (Kitchener et al., 2017). Nevertheless, the genetic support based on intergenic, purportedly neutral variation has been lacking or inconclusive for some of these intraspecific subunits (Bazzicalupo et al., 2022; Lucena-Perez et al., 2020). Genetic structure based on neutral markers, although not strong, was observed between Eurasian lynx populations, both at local and continental scales (Bazzicalupo et al., 2022; Behzadi et al., 2022; Lucena-Perez et al., 2022, 2020; Mengüllüoğlu et al., 2021; Rueness et al., 2003, 2014; Schmidt et al., 2009). The overall consensus of these studies, backed by both nuclear and mitochondrial markers, pointed to the existence of at least three major lineages within the species, one in the East (L. l. wrangeli), one in the Northwest (L. l. lynx) and one in the Southern part of the distribution (L. l. dinniki). Another divergent mitogenomic lineage was suggested to occur in China (L. l. isabellinus) pending on nuclear corroboration (Mengüllüoğlu et al., 2021). While the genetic differentiation observed within the Western part of the Eurasian lynx range has largely been attributed to recent human-related population declines, the genetic differentiation observed between the Eastern, Western and Southern lineages is caused by genetic isolation initiated in late Pleistocene, and intensified throughout the Holocene (Bazzicalupo et al., 2022; Lucena-Perez et al., 2020). In a first attempt to find ecological and environmental drivers of genetic differentiation, microsatellite and mtDNA differentiation was found to be related to snow cover, NAO index and diet in the western part of the Eurasian lynx range (Ratkiewicz et al., 2014), but no studies to date have tried to describe patterns of adaptive population structure at a range-wide scale. In this study, we apply a combination of GEA approaches with the objective of revealing genomic variation carrying possible adaptive significance for the Eurasian lynx. We aimed at identifying areas of the genome responsible for local adaptation in different populations, identifying what genes are involved in these processes, and characterizing patterns of adaptive genetic divergence. We finally discuss the relevance of adaptive population structure in the delimitation of Eurasian lynx taxonomic and conservation units, to inform assisted gene flow and genetic rescue programs. 86 METHODS Sampling, DNA extraction, sequencing and read alignment We analyzed whole genome sequences of a total of 103 Lynx lynx individuals, comprising 12 extant populations distributed across the global range of the species and representing all proposed subspecies, except L. l. isabellinus (Figure 2.1). Most samples were tissues collected from legally hunted individuals (Norway, Latvia, Russia) or animals opportunistically found dead in the field (NE-Poland, Carpathian Mountains, Caucasus, and Mongolia). Blood samples of live-trapped lynx (Białowieża Primeval Forest in NE Poland and Carpathian Mountains, in Poland, and the western part of North Macedonia) were also obtained. The licenses for lynx live-trapping and blood sampling adhere to the Directive 2010/63/EU on the protection of animals used for scientific purposes (see Bazzicalupo et al. (2022) and Lucena-Perez et al. (2020) for details on specific license numbers). No animals were harmed during live-trapping and handling. DNA extractions and sequencing were carried out during previous studies by our research group (Bazzicalupo et al., 2022; Lucena-Perez et al., 2020). Samples were processed by overnight digestion using proteinase K and extracted using silica-coated paramagnetic beads (NucleoMag® Tissue, MACHEREY-NAGEL GmbH & Co. KG). Samples yielding DNA concentration too low to be sequenced were re-extracted either using a classical phenol-chloroform protocol or, in case of one Caucasian sample that was a wet claw, using a protocol meant for bone (Rohland & Hofreiter, 2007) in a sterile lab. Subsequently, gDNA was sheared, size-selected, end-repaired and adenylated following the appropriate Illumina protocol. After ligating indexed paired-end adapters, DNA fragments were amplified via PCR (polymerase chain reaction) if required, and the quantity, quality and size of the libraries were assessed. Finally, libraries were sequenced using Illumina HiSeq2000 or Illumina HiSeq X-10, in Centro Nacional de Análisis Genómico (CNAG) or Macrogen facilities, respectively. In all cases, samples were sequenced using Illumina protocols, and primary data analysis was carried out with the standard Illumina pipeline. We performed a quality control of our data using fastqc (https://www.bioinformatics.babraham.ac.uk/projects/fastqc), and sequencing reads were trimmed using the software trimmomatic 0.38 (Bolger et al., 2014). Sequenced reads were aligned to a 2.4 Gb Lynx canadensis nuclear reference (mLynCan4_v1.p; GCA_007474595.1; Rhie et al., 2021), following the same workflow as our previous analyses (Bazzicalupo et al., 2022; Lucena-Perez et al., 2020). Sample location information, as well as sequencing depth and European Nucleotide Archive (ENA) accession numbers for alignment files, can be found in Table S2.1. 87 Figure 2.1 Map highlighting the distributional range of the Eurasian lynx, delimited by a green line, and locations of samples used in this study (circles, color coded by population). Variant calling and filtering Variant calling was performed on the genome alignments of all the samples using the GATK v4.1.4.1 (McKenna et al., 2010) HaplotypeCaller command. A first dataset of unfiltered variants consisting of 25.52 million SNPs was then filtered, removing variants found in low mappability and repetitive regions, as well as all INDELs and non-biallelic variants. Low quality variants were also filtered out following GATK's suggested thresholds of variant quality (i.e., QUAL ≤ 30, QD ≤ 2.0, FS ≥ 60.0, MQ ≤ 40.0), after inspecting the genome-wide distribution of values of these statistics. Variants that were not represented in at least four samples of each population and with a minimum depth of 3× per sample were filtered out. In order to remove variants in regions where multiple paralogs are collapsed in the reference genome, we removed variants with a read depth exceeding the average read depth plus 1.5 times the standard deviation, calculated using SAMtools depth (Li et al., 2009). The filtered VCF dataset containing the genotype information of all the individuals (whole-set VCF) consisted of 4,983,054 SNPs. From this whole-set VCF we then proceeded to exclude some of the populations defined in Lucena-Perez et al. (2020) and Bazzicalupo et al. (2022) in order to carry out Genome Environment Association (GEA) analyses with confidence (see “Identifying candidate regions of selection” below). We excluded the Balkans, Carpathians, North-Eastern Poland, and Norway populations, which reduced the sample size of the selection scans to 71, with the remaining 8 populations represented by a minimum of 7 and a maximum of 13 individuals. The excluded populations were characterized by high differentiation, lowest levels of genetic diversity and small recent effective population sizes, caused by intense genetic drift (Bazzicalupo et al., 2022; Lucena-Perez et al., 2020), a process that may confound the genomic signals of positive selection (Huber et al., 2016; Weigand & Leese, 2018). We finally proceeded to remove all the variants that were below a 88 minimum allele frequency (MAF) of 0.05, yielding a final VCF consisting of 2,100,553 SNPs (final-set VCF). Environmental predictors Extracting environmental data In order to analyze the association between the genotype and the environment, we extracted environmental data from the geographic area around the coordinates of each of the sequenced individuals using the R package “raster” (Hijmans, 2021). As possible environmental predictors of genetic variation, we used the set of 19 WorldClim bioclimatic variables (Fick & Hijmans, 2017; see Table S2.2 for the complete list of WorldClim variables used, together with their definition and abbreviation adopted in this work). Additionally, we included two variables related to snow precipitation, the mean number of days with snow cover throughout the year (Mean_snow_days) and the average snow depth during the month of January (Jan_mean_depth). These variables were associated with genetic structure in a previous study (Ratkiewicz et al., 2014) and their values for the years between 1999 and 2018 were extracted from the Northern Hemisphere subset of the Canadian Meteorological Centre operational global daily snow depth analysis (Brown & Brasnett, 2010). Environmental predictor values for each individual were calculated by averaging the values within a circle with a diameter of 100 km around the coordinate location of the individual using the R package “geobuffer” (Ștefan, 2019), a diameter that corresponds to the minimal home range described for the species (Breitenmoser & Breitenmoser-Würsten, 2008). Variable selection In order to select which of the WorldClim environmental predictors had the strongest effect on the genetic variance observed in lynx populations, we implemented a predictive approach based on forward selection using redundancy analysis (RDA), as described by Capblancq and Forester (2021). Briefly, environmental predictors are added to an empty RDA model one by one selecting the one that increases the variance explained by the model the most (in terms of adjusted r2 values), until the variance explained by the full model including all the environmental predictors is reached or a permutation-based significance test does not reach the desired significance threshold. We implemented this using the functions rda and ordiR2step of the package “vegan” (Oksanen et al., 2020), using a p-value cut-off of 0.01 and 1000 permutations. Additionally, given their reported importance in determining genetic structure (Ratkiewicz et al., 2014), we decided to add to our final set of environmental predictors the snow-related variable which would increase the variance explained by the forward-selected RDA model the most. To reduce collinearity among our environmental predictors, we excluded from the analysis any variable whose variance inflation factor (VIF) exceeded 7 and Pearson's correlation coefficient (r) exceeded 0.85 with any other variable, keeping in the analysis the variable with highest adjusted r2 contribution to the RDA model. 89 first principal component, as well as occupying the negative extreme of the second principal component. It is also, following Latvia, the second most adaptively distant population to the rest of Eurasian lynx populations. Overall, as the first principal component explained around 30% of the total variation in the adaptive genetic dataset, which represented more than five times the amount explained by the second principal component (<6%), Latvia and NE-Poland are the most and second most adaptively divergent populations, respectively. Other than these two, most populations (Caucasus, Kirov, Mongolia, Primorsky Krai, Tuva, Urals, Yakutia) belong to an adaptively homogeneous group, and the three remaining populations (Balkans, Carpathians, Norway) are slightly differentiated from each other and from the homogeneous group, but to a much smaller degree than NE-Poland and Latvia. 96 Figure 2.2 Adaptive genetic structure as recovered by PCA (top) and population frequency distributions of differences in pairwise Euclidean distances between neutral and adaptive loci (bottom). In the PCA biplots relationships between samples, color coded by population, are represented for the first two principal components (top left) and the third and fourth principal components (top right). (bottom). A frequency distribution graph of the difference between adaptive and neutral distances between the members of the population and the individuals of each of the other populations (color-coded). A vertical dashed line at the zero separates pairs of populations that are adaptively similar but neutrally divergent, whose differential distances are zero or negative, and adaptively divergent pairs, whose differential distances are positive. Spatial gradients of adaptive genetic differentiation – GDM Results from the GDM analysis allow us to visualize the expected variation in genetic composition in geographic space (turnover) as a response to geographic location and different environmental predictors (Figure 2.3 and Figure S2.7). The percentage of null deviance explained by the fitted GDM model was 57.25% for the neutral dataset and 62.14% for the adaptive dataset. Biplots of the first two principal components with labelled vectors are also included in Figure 2.3 and Figure S2.7, giving an idea of the magnitude and direction of the most correlated environmental predictors. Relatively homogeneous turnover is observed for neutral genetic composition (Figure S2.7). Overall, changes in neutral genetic compositions are predicted to be mainly a consequence of geographic location (Figure S2.7), affected by both latitude and longitude (biplot in Figure S2.7 and Table S2.7). Other environmental predictors have much lower relative importance, with the highest non-geographic predictors being isothermality (Iso_T), with less than half the relative importance of geography (Figure S2.7). In sharp contrast to neutral variation, the turnover of adaptive genetic composition is greatly influenced by two environmental predictors, T_range_day and Mean_snow_days, with 60.35% and 32.77% of the overall importance respectively (Figure 2.3 and Table S2.8). Throughout most of the Northern part of the species range, adaptive genetic composition is determined by the effect of high values of Mean_snow_days, with predicted composition shifts between areas affected by lower or higher values of T_range_day (range of darker to lighter green in Figure 2.3, see biplot). Most of the populations sampled in this study fall within these regions of high values of both these variables, including all the Asian populations (Tuva, Mongolia, Yakutia, Primorsky Krai), the populations from Western Russia (Kirov, Urals) and Scandinavia (Norway). Patches of similar adaptive compositions are predicted to be found in highland regions of the Southern part of the range. These are the mountainous ranges of Europe (occupied by the Carpathian and Balkan populations), as well as the Caucasian and Himalayan mountains. The rest of the Southern range is characterized by either lower Mean_snow_days but high T_range_day, as in the lowland areas of the European range (where we sampled the NE-Poland population and part of the Latvia population), Anatolia, and Iran (pink areas in Figure 2.3, see biplot) or lower values of both, as in the Latvian Courland Peninsula (where the rest of the Latvia population is sampled) and the southern tip of Scandinavia (purple and blue areas in Figure 2.3, see biplot). The differences in turnover patterns between the neutral and adaptive datasets are highlighted by 97 the geographical representation of Procrustes superimposition residuals (Figure S2.8). Higher residual values are observed as we move along the longitudinal range, as the main difference is represented by the effect of geography, which strongly influences neutral genetic composition, but does not affect the adaptive one. Figure 2.3 Predicted turnover in genetic composition as calculated by GDM on adaptive genetic markers. Locations with similar colors are expected to have similar adaptive genetic composition, based on the effects of the different environmental predictors. The contribution of different variables to genetic composition is indicated in the associated biplot, with vectors pointing to the direction of the contributing variable and indicating the magnitude of its contribution. The bar plot indicates the relative importance of each environmental variable to the overall GDM model calculated as the rescaled maximum value of the model's fitted I-Splines. Circles represent sampled individuals, and a green line delimits Eurasian lynx current distributional range (Breitenmoser et al., 2015). A map zoomed on the range occupied by the adaptively divergent Latvia and NE-Poland populations is presented to appreciate their distinct expected genetic compositions. DISCUSSION In this study, we examined range-wide patterns of adaptive differentiation driven by environmental conditions across the distributional range of the Eurasian lynx. Loci potentially involved in local adaptation were identified through their association with environmental predictors. We assessed how local adaptation in Eurasian lynx populations affects the global genetic structure of this large carnivore species by describing the way neutral and adaptive genetic differentiation patterns differ, and by estimating both the magnitude and the direction of these changes. We also described how genetic composition is expected to change throughout the species' range in response to geographic distance and environmental factors. By analyzing the changes in adaptive genetic compositions across the distributional range of the species we discovered areas of expected strong adaptive divergence driven by specific environmental conditions. We also identified the biological processes that were the main targets of this adaptive differentiation as those overrepresented among the genes with variation associated with environmental variables. 98 Evidence of local adaptation The genetic structure reconstructed with candidate adaptive loci was substantially different to that reconstructed with neutral ones. The three main clades recovered in the neutral genome, that divide Eurasian lynx populations in three main lineages (Bazzicalupo et al., 2022; Lucena-Perez et al., 2020) are not observed when looking at the areas of the genome identified as responding adaptively to environmental conditions. From the adaptive PCA a highly homogeneous group is instead recovered, which includes representatives from three distinct neutral clades, the Caucasus, the Western lynx clade (Kirov, Urals) and the Eastern lynx clade (Mongolia, Primorsky Krai, Tuva, Yakutia), with adaptive genetic distances being lower than neutral genetic distances in all pairwise comparisons. A similar pattern is also observed when mapping the turnover in expected genetic composition in response to environmental conditions. Across the north-western, central, and eastern part of the distribution, areas occupied by lynxes of this adaptively homogenous group, genetic composition is modelled to be mostly homogeneous and largely driven by a high mean number of days with snow cover (Mean_snow_days) and high mean diurnal temperature ranges (T_range_day). Given the low adaptive distance among the populations and the similar environmental pressures they seem to face, it appears as the neutral differentiation observed among these populations is consistent with a scenario of pure neutral evolution driven by historical isolation and limitations in gene flow caused by distance, as previously discussed (Bazzicalupo et al., 2022; Lucena-Perez et al., 2020). On the other hand, the variants involved in the adaptation to local environmental conditions maintain more similar allelic frequencies among neutrally differentiated populations and may be thus more resistant to genetic drift. The lack of signals of differential adaptation within this adaptively homogeneous populations contrasts with the extensive variation in phenotype (e.g. pelage pattern and size) and in habitat and prey across the species' wide distributional range (Breitenmoser & Breitenmoser-Würsten, 2008; Darul et al., 2022; Kitchener et al., 2017). This suggests that either the observed phenotypic variation is evolving as neutral, or that it is a plastic response to the environment. In fact, pelage pattern variation was found to be poorly associated to local habitat and most closely related to least-cost distances from inferred glacial refugia, suggesting that the current patterns of pelage variation are more of a legacy of past phylogeographic history since the last glacial maximum than the outcome of local adaptation to the environment (Darul et al., 2022). On the other hand, epigenetic differentiation was found to be responsible for the variation in size between Alaskan and insular Newfoundland populations of the closely related Canada lynx, in a background of neutral homogeneity (Meröndun et al., 2019). This suggests that even though we are not observing adaptive divergences among Eurasian lynx in the north-west, central, and eastern parts of the distribution, environmental conditions could still be potentially inducing adaptive plastic responses in phenotypic traits through epigenetic changes. Additionally, environmental variables not included in this study, such as soil type or land cover, could also be driving local adaptations within this homogeneous group that we have not been able to detect in this study. Overall, it seems that ecological flexibility might have allowed this homogeneous group to colonize widespread areas with diverse climatic and ecological conditions. It will be interesting to address whether that is also the case for other felids, such 99 as jaguars (Eizirik et al., 2001) and pumas (Culver, 2000), widespread species which may also base their evolutionary success as apex predators on phenotypic plasticity rather than natural selection on genomic variants. Most of the divergence observed in adaptive loci is represented by the differentiation of the samples belonging to the westernmost region of the Baltic population of Eurasian lynx, represented in our study by the Latvia and NE-Poland populations. Samples from both populations are differentiated in the adaptive PCA and show higher adaptive than neutral distances from the rest of samples. Genetic composition in the regions occupied by these populations is again mostly driven by two environmental predictors, the mean number of days with snow cover throughout the year (Mean_snow_days) and the ranges in mean diurnal temperature (T_range_day). Low values of both variables are expected to determine the genetic composition of lynx inhabiting the Latvian Courland Peninsula, while higher Mean_snow_days characterize the rest of the Latvia population and higher T_range_day influences the genetic composition of NE-Poland. The high adaptive differentiation observed in these populations and the distinctive environmental conditions they are subjected to suggest an ongoing process of adaptive differentiation in the regions occupied by these populations. This contrasts with what is observed when looking at the neutral genome, where individuals from Latvia and NE-Poland are relatively closer to other members of what is known as the European lowland or Baltic lynx populations, which also includes members of the Western-Russia populations of Kirov and Urals (Lucena-Perez et al., 2020; Schmidt et al., 2021). Some evidence of genetic structuring between the westernmost Baltic lynx individuals and the Scandinavian and West-Russian was already detected when reconstructing their mitochondrial phylogenomic relationships, which revealed the presence of at least two distinct haplogroups (Lucena-Perez et al., 2020; Mengüllüoğlu et al., 2021), one of which is in fact more related to other Eurasian lynx from the Asian part of the distributional range. It could be possible that adaptive divergence to different environmental conditions occurred during the isolation of these more anciently divergent mitochondrial clades in separate refugia with different environments. Given that the adaptive genetic composition of the previously inhabited lowland areas of the European landscape is predicted to be similar to that of these adaptively divergent populations, these adaptations may have been widespread across Europe in the past. In that case, we would now be observing the remnants of a differentially adapted lowland lynx following the extensive habitat reduction and fragmentation that occurred in the area. Gene-flow between these adaptively diverging populations is however still ongoing, as demonstrated by their sharing of mitochondrial haplotypes with other western lynx and their relative homogeneity in the neutral genome (Lucena-Perez et al., 2020; Mengüllüoğlu et al., 2021). Weaker signals of differentiation in the adaptive PCA are observed for individuals of the Balkans, Carpathians, and Norway populations, whose adaptive distances to members of other populations than their own do not differ from the corresponding neutral distances or are only slightly higher. These populations occupy regions with environmental conditions that more closely resemble the ones observed in the northern-most part of the distribution, characterized by high Mean_snow_days and T_range_day, as the consequence of either their 100 altitudinal position, in the mountainous ranges of the Carpathians and the Balkans, or their northward latitudinal position, in the case of Norway. As their neutral and adaptive distances to the rest of populations are of similar magnitude, and they do not appear to occupy environments with unique climatic conditions, their overall genomic differentiation is most parsimoniously interpreted as the consequence of genetic drift augmented by isolation and low effective population sizes. It must be noted, however, that these populations were excluded from our search for adaptive variation specifically because of intense recent genetic drift, so we may have missed potential variants underlying adaptations exclusive to these populations. Drivers of local adaptation and genes involved Variance partitioning analysis of the RDA model, although limited by the number of observations we have compared to the number of variables considered, allowed us to disentangle the effects of demographic history, geographic location, and environmental conditions on genetic variation in Eurasian lynx, and revealed that part of the genetic variation among individuals could be attributed to local climatic conditions. Gradients in annual mean diurnal ranges in temperature (T_range_day) and the mean number of days with snow cover throughout the year (Mean_snow_days) appear to have the strongest effects on adaptive genetic variation of Eurasian lynx populations, indicating that these environmental variables – or others closely correlated with these (see below) – might be important factors for the biology and evolution of this species. Apart from its previously reported correlations with skull morphology and diet (Yom-Tov et al., 2011), snow precipitation was also already suggested as a possible predictor of genetic differentiation in the western range of the Eurasian lynx, based on both mtDNA and microsatellite markers (Ratkiewicz et al., 2014; Schmidt et al., 2011). It was suggested that familiarity and adaptation to local snow conditions would increase efficiency in prey capture and discourage dispersal to other climatic regions (Nilsen et al., 2009), similarly to how snow conditions are predicted to affect population differentiation and dynamics in Canada lynx (Stenseth, Ehrich, et al., 2004; Stenseth, Shabbar, et al., 2004). It is hard to predict the selective pressures that diurnal ranges in temperature generate on lynx populations. The identification of Mean_snow_days and T_range_day as the most important variables in shaping adaptive genetic composition should be taken with a grain of salt, as it might be the result of an indirect correlation. Other climatic or ecological variables correlated to these two might be the actual causal factors, and the effect we observe here could be a byproduct of their correlation. In particular, local adaptations in response to T_range_day might be associated with specific mechanisms that determine how Eurasian lynx directly cope with temperature fluctuations. On the other hand, T_range_day might represent a proxy of ecological environments, as diurnal ranges in temperature have been linked to forest structure heterogeneity (Ehbrecht et al., 2019), which ultimately affects the ecological community composition of an area (Levine & HilleRisLambers, 2009). Indeed, a significant effect of the habitat structure, mediated by stalking cover availability, on habitat suitability and on the 101 distribution of Eurasian lynx was recently reported in central Europe (Schmidt et al., 2023). Specialization on the hunting of particular prey, mainly larger ungulates, is being observed for lynx inhabiting the Baltic region (Valdmann et al., 2005), which are here recovered as the most adaptively divergent. Local climatic conditions might be facilitating the acquisition of these larger prey, allowing locally adapted lynx to thrive, and even colonize new territories. With the data at hand, precise descriptions of the mechanisms by which the genes involved in local adaptation directly affect the fitness of Eurasian lynx cannot be made. Our results suggest that four main categories of biological functions are carried out by the genes involved in local adaptation to the environment: metabolic processing of proteins and sugars by O-glycosylation, locomotory and general behavior, synaptic regulation and organization, and nervous system development. O-linked glycosylation is a type of protein glycosylation that adds sugar molecules to specific serine or threonine residues on proteins. This modification can have a very wide range of effects on a number of different biological processes, by influencing protein conformation and expression, by regulating processes mediated by the organisms molecular interaction events, such as fertilization or immune responses, and by modulating the activity of signaling molecules such as hormones (Hounsell et al., 1996; Van den Steen et al., 1998). Strong adaptive pressures on behavior and neural regulation, organization and development might be common among hypercarnivores, as enrichment in genes with similar functions has also been observed when analyzing other carnivores' evolution (Mittal et al., 2019; Scholl et al., 2021). Broadly speaking, local adaptation might be affecting the way lynx are experiencing, navigating, and regulating their physiology in the different ecological conditions imposed by snow, temperature, and other possibly associated ecological variation. Implications for conservation and taxonomy The neutral variation observed in the Eurasian lynx and the adaptive variation described in this study are not similar in the way they differentiate among populations and subspecific clades. The divergent groups corresponding to the subspecies of L. l. lynx in the west, L. l. wrangeli in the east, and L. l. dinniki in the Caucasian mountains, revealed by the neutral genome (Bazzicalupo et al., 2022; Lucena-Perez et al., 2020), are not observed when looking at adaptive variation. Across all of the sampled range, except for the westernmost part of the distribution, these three divergent lineages appear as adaptively homogeneous. Their divergence is therefore most probably a consequence of the varying interaction of genetic drift and gene flow across space and time, it is strongly correlated with geographic location and with isolation in separate glacial refugia, and does not appear to be driven by adaptive processes. As for the remaining subspecies we have sampled, L. l. balcanicus from the Balkans and L. l. carpathicus from the Carpathian Mountains, we detect some differentiation at adaptive variants, which occurs however along less important axis of adaptive differentiation and is much weaker than that observed with neutral variants. It is thus likely that genetic drift might have affected allele frequencies in these otherwise selected variants (Bazzicalupo et al., 2022; Lucena-Perez et al., 2020). 102 The strongest signals of adaptive divergence driven by the environment were observed in the populations in Latvia and NE-Poland, which seem to be locally adapted to different snow conditions and diurnal temperature fluctuations, experiencing a relatively mild climate when compared, for example, to northern Scandinavia or the mountainous ridges of the Carpathians and the Caucasus. However, without additional support from phenotypic and ecological data, we advise against using this evidence to classify Baltic populations as a novel subspecies separate from L. l. lynx. Apart from taxonomic inferences, understanding how adaptive variation shapes the relationships between different populations of Eurasian lynx can assist the identification of evolutionary significant units (ESUs; Funk et al., 2012; Ryder, 1986) and ecotypes (Stronen et al., 2022), with the aim of guiding conservation management decisions (Teixeira & Huber, 2021). Our results suggest that differences in local environmental conditions, especially linked to snow precipitation and daily temperature fluctuations, should be an important factor to consider when delimiting Eurasian lynx ESUs. Within the distributional range of the L. l. lynx subspecies, we were able to identify at least two possibly distinct adaptive units. One of them is adapted to harsher conditions of high snow and temperature fluctuations, which is widely distributed in areas of Western-Russia and Scandinavia, and the other is adapted to the milder climatic conditions of European lowlands. Conservation and restoration programs might want to take these two adaptive units into consideration when selecting potential sources for reintroduction or reinforcement. In particular, reintroductions in western Europe, currently based mostly on Carpathian lynx (Mueller et al., 2022), might consider the Baltic population as an additional and valuable source. Indeed, the reintroduced lynx population in western Poland is assumed to be founded by individuals that, based on their genotypes, correspond to the Baltic population (Tracz et al., 2021). A more challenging case is represented by the Balkan lynx, L. l. balcanicus, whose adaptive differentiation observed in this study is hard to disentangle from the neutral divergence caused by intense genetic drift and prolonged low effective population sizes. Our previous assessment of the Carpathian population as representing the best source of individuals for genetic rescue of this highly endangered population, based on its closest relationship in neutral autosomal variation (Bazzicalupo et al., 2022) and also supported by their ecological similarity (Melovski et al., 2022), could be additionally confirmed by their relative overall adaptive similarity and by their similar environment defined by the most relevant environmental variables. Our results are also relevant for evaluating the long-term effects of climate change on the Eurasian lynx distribution. The two environmental variables identified as drivers of adaptive divergence are likely to alter their geographical distribution in the near future due to climate change, so associated variants may also expand and contribute to species resilience and persistence. In other species, the occurrence of adaptive responses to environmental pressures has been shown to reduce range loss projections under future climatic conditions (Razgour et al., 2019). It is thus critically important to conserve the range of adaptive variation within species by protecting populations harboring locally adapted variants if we aim to assure their long-term viability. 103 REFERENCES Arnegard, M. E., McGee, M. D., Matthews, B., Marchinko, K. B., Conte, G. L., Kabir, S., Bedford, N., Bergek, S., Chan, Y. F., Jones, F. C., Kingsley, D. M., Peichel, C. L., & Schluter, D. (2014). Genetics of ecological divergence during speciation. Nature, 511(7509), 307–311. https://doi.org/10.1038/nature13301 Ashburner, M., Ball, C. A., Blake, J. A., Botstein, D., Butler, H., Cherry, J. M., Davis, A. P., Dolinski, K., Dwight, S. S., Eppig, J. T., Harris, M. A., Hill, D. P., Issel-Tarver, L., Kasarskis, A., Lewis, S., Matese, J. C., Richardson, J. E., Ringwald, M., Rubin, G. M., & Sherlock, G. (2000). Gene Ontology: Tool for the unification of biology. Nature Genetics, 25(1), 25–29. https://doi.org/10.1038/75556 Bazzicalupo, E., Lucena‐Perez, M., Kleinman‐Ruiz, D., Pavlov, A., Trajçe, A., Hoxha, B., Sanaja, B., Gurielidze, Z., Kerdikoshvili, N., Mamuchadze, J., Yarovenko, Y. A., Akkiev, M. I., Ratkiewicz, M., Saveljev, A. P., Melovski, D., Gavashelishvili, A., Schmidt, K., & Godoy, J. A. (2022). History, demography and genetic status of Balkan and Caucasian Lynx lynx (Linnaeus, 1758) populations revealed by genome‐wide variation. Diversity and Distributions, 28(1), 65–82. https://doi.org/10.1111/ddi.13439 Behzadi, F., Malekian, M., Fadakar, D., Adibi, M. A., & Bärmann, E. V. (2022). Phylogenetic analyses of Eurasian lynx (Lynx lynx Linnaeus, 1758) including new mitochondrial DNA sequences from Iran. Scientific Reports, 12(1), 3293. https://doi.org/10.1038/s41598-022-07369-z Beissinger, T. M., Rosa, G. J., Kaeppler, S. M., Gianola, D., & de Leon, N. (2015). Defining window-boundaries for genomic analyses using smoothing spline techniques. Genetics Selection Evolution, 47(1), 30. https://doi.org/10.1186/s12711-015-0105-9 Bolger, A. M., Lohse, M., & Usadel, B. (2014). Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics, 30(15), 2114–2120. https://doi.org/10.1093/bioinformatics/btu170 Bolnick, D. I., & Otto, S. P. (2013). The magnitude of local adaptation under genotype‐dependent dispersal. Ecology and Evolution, 3(14), 4722–4735. https://doi.org/10.1002/ece3.850 Breitenmoser, U., & Breitenmoser-Würsten, C. (2008). Der Luchs – Ein Grossraubtier in der Kulturlandschaft. Salm Verlag. Breitenmoser, U., Breitenmoser-Würsten, C., Lanz, T., von Arx, M., & Antonevich, A. (2015). Lynx lynx (errata version published in 2017). IUCN Red List of Threatened Species. https://www.iucnredlist.org/en Brown, Ross, & Brasnett, Bruce. (2010). Canadian Meteorological Centre (CMC) Daily Snow Depth Analysis Data, Version 1 [Dataset]. NASA National Snow and Ice Data Center DAAC. https://doi.org/10.5067/W9FOYWH0EQZ3 Capblancq, T., & Forester, B. R. (2021). Redundancy analysis: A Swiss Army Knife for landscape genomics. Methods in Ecology and Evolution, 12(12), 2298–2309. https://doi.org/10.1111/2041-210X.13722 104 Culver, M. (2000). Genomic ancestry of the American puma (Puma concolor). Journal of Heredity, 91(3), 186–197. https://doi.org/10.1093/jhered/91.3.186 Darul, R., Gavashelishvili, A., Saveljev, A. P., Seryodkin, I. V., Linnell, J. D. C., Okarma, H., Bagrade, G., Ornicans, A., Ozolins, J., Männil, P., Khorozyan, I., Melovski, D., Stojanov, A., Trajçe, A., Hoxha, B., Dvornikov, M. G., Galsandorj, N., Okhlopkov, I., Mamuchadze, J., … Schmidt, K. (2022). Coat Polymorphism in Eurasian Lynx: Adaptation to Environment or Phylogeographic Legacy? Journal of Mammalian Evolution, 29(1), 51–62. https://doi.org/10.1007/s10914-021-09580-7 de Queiroz, K. (2007). Species Concepts and Species Delimitation. Systematic Biology, 56(6), 879–886. https://doi.org/10.1080/10635150701701083 Doebeli, M., & Dieckmann, U. (2004). Adaptive Dynamics of Speciation: Spatial Structure. In U. Dieckmann, M. Doebeli, J. A. J. Metz, & D. Tautz (Eds.), Adaptive Speciation (1st ed., pp. 140–168). Cambridge University Press. https://doi.org/10.1017/CBO9781139342179.008 Duforet-Frebourg, N., Bazin, E., & Blum, M. G. B. (2014). Genome Scans for Detecting Footprints of Local Adaptation Using a Bayesian Factor Model. Molecular Biology and Evolution, 31(9), 2483–2495. https://doi.org/10.1093/molbev/msu182 Ehbrecht, M., Schall, P., Ammer, C., Fischer, M., & Seidel, D. (2019). Effects of structural heterogeneity on the diurnal temperature range in temperate forest ecosystems. Forest Ecology and Management, 432, 860–867. https://doi.org/10.1016/j.foreco.2018.10.008 Eizirik, E., Kim, J.-H., Menotti-Raymond, M., Crawshaw JR., P. G., O’Brien, S. J., & Johnson, W. E. (2001). Phylogeography, population history and conservation genetics of jaguars (Panthera onca, Mammalia, Felidae). Molecular Ecology, 10(1), 65–79. https://doi.org/10.1046/j.1365-294X.2001.01144.x Fick, S. E., & Hijmans, R. J. (2017). WorldClim 2: New 1‐km spatial resolution climate surfaces for global land areas. International Journal of Climatology, 37(12), 4302–4315. https://doi.org/10.1002/joc.5086 Fitzpatrick, M. C., & Keller, S. R. (2015). Ecological genomics meets community-level modelling of biodiversity: Mapping the genomic landscape of current and future environmental adaptation. Ecology Letters, 18(1), 1–16. https://doi.org/10.1111/ele.12376 Fitzpatrick, M., Mokany, K., Manion, G., Nieto-Lugilde, D., & Ferrier, S. (2021). gdm: Generalized Dissimilarity Modeling (Version 1.5.0-1) [R]. https://CRAN.R-project.org/package=gdm Forester, B. R., Lasky, J. R., Wagner, H. H., & Urban, D. L. (2018). Comparing methods for detecting multilocus adaptation with multivariate genotype-environment associations. Molecular Ecology, 27(9), 2215–2233. https://doi.org/10.1111/mec.14584 Funk, W. C., McKay, J. K., Hohenlohe, P. A., & Allendorf, F. W. (2012). Harnessing genomics for delineating conservation units. Trends in Ecology & Evolution, 27(9), 489–496. https://doi.org/10.1016/j.tree.2012.05.012 105 c_ll_mo_0187 Mongolia L. l. wrangeli Mongolia 103.78 48.64 f 5.85 yes c_ll_mo_0190 Mongolia L. l. wrangeli Mongolia 108.27 47.40 f 7.58 yes c_ll_mo_0191 Mongolia L. l. wrangeli Mongolia 108.68 47.63 m 8.05 yes c_ll_ki_0090 Kirov L. l. lynx Russia 50.40 59.81 f 21.68 yes c_ll_ki_0091 Kirov L. l. lynx Russia 50.40 59.81 m 6.31 yes c_ll_ki_0092 Kirov L. l. lynx Russia 50.40 59.81 f 5.37 yes c_ll_ki_0093 Kirov L. l. lynx Russia 50.40 59.81 m 5.46 yes c_ll_ki_0094 Kirov L. l. lynx Russia 50.40 59.81 m 6.06 yes c_ll_ki_0095 Kirov L. l. lynx Russia 48.42 60.81 m 5.52 yes c_ll_ki_0096 Kirov L. l. lynx Russia 48.42 60.81 m 6.37 yes c_ll_ki_0097 Kirov L. l. lynx Russia 48.42 60.81 f 6.12 yes c_ll_ki_0098 Kirov L. l. lynx Russia 48.42 60.81 f 5.88 yes c_ll_ki_0099 Kirov L. l. lynx Russia 48.42 60.81 f 6.1 yes c_ll_ki_0100 Kirov L. l. lynx Russia 46.73 61.19 m 6.42 yes c_ll_ki_0101 Kirov L. l. lynx Russia 46.73 61.19 f 6.38 yes c_ll_ki_0102 Kirov L. l. lynx Russia 46.73 61.19 m 6.34 yes c_ll_la_0044 Latvia L. l. lynx Latvia 22.53 56.64 m 7.76 yes c_ll_la_0045 Latvia L. l. lynx Latvia 27.58 56.31 f 23.97 yes c_ll_la_0047 Latvia L. l. lynx Latvia 28.15 56.38 f 7.02 yes c_ll_la_0048 Latvia L. l. lynx Latvia 22.82 57.36 f 5.2 yes c_ll_la_0052 Latvia L. l. lynx Latvia 22.60 57.35 m 6.74 yes c_ll_la_0053 Latvia L. l. lynx Latvia 22.82 57.05 f 7.67 yes c_ll_la_0054 Latvia L. l. lynx Latvia 25.57 57.67 m 8.47 yes c_ll_no_0065 Norway L. l. lynx Norway 21.34 69.51 f 22.98 no c_ll_no_0075 Norway L. l. lynx Norway 12.36 64.51 m 5.59 no c_ll_no_0076 Norway L. l. lynx Norway 10.68 63.82 m 6 no c_ll_no_0077 Norway L. l. lynx Norway 10.95 63.57 m 5.1 no c_ll_no_0078 Norway L. l. lynx Norway 12.15 61.38 m 5.39 no c_ll_no_0079 Norway L. l. lynx Norway 8.41 59.40 m 5.11 no c_ll_no_0080 Norway L. l. lynx Norway 8.73 60.66 f 5.3 no c_ll_no_0081 Norway L. l. lynx Norway 8.97 60.27 f 5.59 no c_ll_no_0082 Norway L. l. lynx Norway 9.03 60.14 m 5.51 no c_ll_po_0001 North-Eastern Poland L. l. lynx Poland 23.60 52.95 m 6.73 no c_ll_po_0002 North-Eastern Poland L. l. lynx Poland 23.84 52.82 f 6.67 no c_ll_po_0003 North-Eastern Poland L. l. lynx Poland 23.72 52.80 f 6.49 no c_ll_po_0011 North-Eastern Poland L. l. lynx Poland 23.51 52.89 m 6.61 no c_ll_po_0014 North-Eastern Poland L. l. lynx Poland 23.95 52.77 m 6.61 no c_ll_po_0019 North-Eastern Poland L. l. lynx Poland 23.74 52.68 m 6.59 no 112 c_ll_po_0105 North-Eastern Poland L. l. lynx Poland 23.47 53.09 m 5.44 no c_ll_po_0106 North-Eastern Poland L. l. lynx Poland 23.66 53.10 f 5.76 no c_ll_po_0150 North-Eastern Poland L. l. lynx Poland 23.60 52.76 m 23.99 no c_ll_tu_0153 Tuva L. l. wrangeli Russia 96.53 51.48 f 8.25 yes c_ll_tu_0154 Tuva L. l. wrangeli Russia 96.53 51.48 f 22.02 yes c_ll_tu_0157 Tuva L. l. wrangeli Russia 96.53 51.48 m 8.16 yes c_ll_tu_0158 Tuva L. l. wrangeli Russia 96.53 51.48 f 8.15 yes c_ll_tu_0159 Tuva L. l. wrangeli Russia 96.53 51.48 f 7.91 yes c_ll_tu_0165 Tuva L. l. wrangeli Russia 96.53 51.48 f 7.88 yes c_ll_tu_0166 Tuva L. l. wrangeli Russia 96.53 51.48 m 8.38 yes c_ll_ur_0194 Urals L. l. lynx Russia 59.64 56.05 f 11.29 yes c_ll_ur_0195 Urals L. l. lynx Russia 60.13 55.15 m 12.27 yes c_ll_ur_0196 Urals L. l. lynx Russia 60.13 55.15 m 12.91 yes c_ll_ur_0199 Urals L. l. lynx Russia 59.79 55.18 m 13.83 yes c_ll_ur_0200 Urals L. l. lynx Russia 59.00 55.00 m 13.33 yes c_ll_ur_0202 Urals L. l. lynx Russia 60.89 54.38 f 23.47 yes c_ll_ur_0203 Urals L. l. lynx Russia 57.36 55.00 f 13.79 yes c_ll_vl_0107 Primosky Krai L. l. lynx Russia 136.47 44.94 f 6.54 yes c_ll_vl_0108 Primosky Krai L. l. lynx Russia 135.58 45.29 f 11.58 yes c_ll_vl_0109 Primosky Krai L. l. wrangeli Russia 136.73 45.67 m 5.94 yes c_ll_vl_0110 Primosky Krai L. l. wrangeli Russia 132.91 49.01 m 7.08 yes c_ll_vl_0112 Primosky Krai L. l. wrangeli Russia 137.01 45.86 f 29.33 yes c_ll_vl_0113 Primosky Krai L. l. wrangeli Russia 137.01 45.86 m 8.75 yes c_ll_vl_0114 Primosky Krai L. l. wrangeli Russia 137.61 47.29 f 23.4 yes c_ll_vl_0128 Primosky Krai L. l. wrangeli Russia 137.43 45.91 m 8.29 yes c_ll_vl_0132 Primosky Krai L. l. wrangeli Russia 137.43 45.91 m 8.66 yes c_ll_vl_0137 Primosky Krai L. l. wrangeli Russia 133.46 51.39 f 25.47 yes c_ll_ya_0138 Yakutia L. l. wrangeli Russia 132.07 59.90 m 8.21 yes c_ll_ya_0139 Yakutia L. l. wrangeli Russia 129.20 61.19 m 8.28 yes c_ll_ya_0140 Yakutia L. l. wrangeli Russia 129.16 61.16 m 8.48 yes c_ll_ya_0141 Yakutia L. l. wrangeli Russia 136.18 66.89 f 23.6 yes c_ll_ya_0142 Yakutia L. l. wrangeli Russia 136.18 66.89 m 8.63 yes c_ll_ya_0143 Yakutia L. l. wrangeli Russia 136.18 66.89 m 8.14 yes c_ll_ya_0145 Yakutia L. l. wrangeli Russia 129.45 61.89 m 8.59 yes c_ll_ya_0146 Yakutia L. l. wrangeli Russia 130.68 60.78 m 23.58 yes c_ll_ya_0147 Yakutia L. l. wrangeli Russia 127.30 60.75 f 8.26 yes 113 Table S2.2. List of WorldClim bioclimatic variables used in this study, including the description provided by the WorldClim website and the abbreviation adopted in this study. WorldClim variable Meaning Abbreviation BIO1 Annual Mean Temperature T_mean_year BIO2 Mean Diurnal Range (Mean of monthly (max temp – min temp)) T_range_day BIO3 Isothermality (BIO2/BIO7) (×100) Iso_T BIO4 Temperature Seasonality (standard deviation ×100) T_seasonality BIO5 Max Temperature of Warmest Month T_max_warm BIO6 Min Temperature of Coldest Month T_min_cold BIO7 Temperature Annual Range (BIO5-BIO6) T_range_year BIO8 Mean Temperature of Wettest Quarter T_wet_quart BIO9 Mean Temperature of Driest Quarter T_dry_quart BIO10 Mean Temperature of Warmest Quarter T_warm_quart BIO11 Mean Temperature of Coldest Quarter T_cold_quart BIO12 Annual Precipitation P_annual BIO13 Precipitation of Wettest Month P_wet_month BIO14 Precipitation of Driest Month P_dry_month BIO15 Precipitation Seasonality (Coefficient of Variation) P_seasonality BIO16 Precipitation of Wettest Quarter P_wet_quart BIO17 Precipitation of Driest Quarter P_dry_quart BIO18 Precipitation of Warmest Quarter P_warm_quart BIO19 Precipitation of Coldest Quarter P_cold_quart 114 Table S2.3. Summary of the forward model building, reporting variables, their adjusted r2, degrees of freedom (Df), AIC value, F value and p-value, in the order the were added to the null model. Variable Adjusted r2 Df AIC F p-value + T_dry_quart 0.027322 1 950.94 2.9101 0.002 + T_range_day 0.059092 1 949.62 3.2622 0.002 + P_seasonality 0.088930 1 948.34 3.1616 0.002 + P_wet_quart 0.098198 1 948.56 1.6680 0.004 + P_wet_month 0.108735 1 948.67 1.7566 0.002 All variables 0.154015 Table S2.4. Summary of the analysis of variance partitioning. Inertia, r2, p-value and proportion of both explainable and total variance are reported for the full model, and the three distinct partial models. Confounded variance is the amount of the full model inertia that is not explainable by partial models. Partial RDA models Inertia r-squared p-value Proportion of explainable Variance Proportion of total Variance Full model: F ~ clim. + geog. + struct. 252685 0.2615709 0.001*** 1 0.26 Pure climate: F ~ clim. | (geog. + struct.) 83295 0.086 0.001*** 0.33 0.09 Pure structure: F ~ struct. | (clim. + geog.) 55197 0.057 0.001*** 0.22 0.06 Pure geography: F ~ geog. | (clim. + struct.) 42160 0.044 0.001*** 0.17 0.04 Confounded climate/structure/geography 72033 0.29 0.07 Total unexplained 713344 0.74 Total inertia 966029 1.00 115 Table S2.5. Summary of GenWin analysis for window boundary definition based on Bayes Factor values calculated by BayPass for each predictor environmental variable, and their overlap with candidate SNPs from RDA. Variable Number of Windows Mean Length (bp) Minimum Length (bp) Maximum Length (bp) SD Length (bp) Windows with RDA SNPs Max RDA SNPs in Window T_dry_quart 1335 20951 10000 70000 10789 96 7 T_range_day 1342 20447 10000 70000 10295 108 8 P_seasonality 1329 21128 10000 70000 10567 54 5 P_wet_quart 1338 20440 10000 80000 10458 68 4 Mean_snow_days 1325 20943 10000 90000 10712 279 9 Table S2.7. Rotation values of the different environmental predictors in the PCA conducted on the GDM transformed raster layers using neutral loci. Predictor PC1 PC2 PC3 x-coord 0.999 -0.003 -0.03 y-coord -0.023 -0.544 -0.498 T_mean_year 0.001 0.141 0.048 T_range_day 0.003 0.062 -0.008 Iso_T -0.02 0.799 -0.392 T_max_warm 0.012 0.114 0.22 T_warm_quart 0.002 0.059 0.035 P_seasonality 0.026 0.158 0.039 P_warm_quart 0.009 0.004 0.737 Table S2.8. Rotation values of the different environmental predictors in the PCA conducted on the GDM transformed raster layers using candidate loci. Predictor PC1 PC2 PC3 T_range_day 0.916 0.401 0.005 Jan_mean_depth -0.082 0.175 0.981 Mean_snow_days -0.393 0.899 -0.193 116 Figure S2.1. Correlogram of the scores of the first two principal components of a neutral PCA, environmental predictors and geographic locations. 117 Figure S2.2. Comparison of eigenvalues percentages of the neutral PCA to a broken stick model. Figure S2.3. Scree plot representing the amount of variance explained by each axis of the redundancy analysis, including the 5 selected environmental predictors and the pruned SNPs dataset. 118 Figure S2.4. SNP and sample loadings on RDA axes 1 and 2. Non-significant SNPs are represented by a dark grey cross. Significant SNPs (small circles) are color coded to reflect which predictor variable (vector lines) they most strongly correlate with. Samples (big circles) are color-coded by their population. 119 Figure S2.5. Manhattan plots showing the results of the intersection of multivariate (RDA) and genomic windows delimited by GenWin, based on individual (BayPass) analyses of each of the five environmental predictors. Non-candidate windows are represented by gray crosses. Outlier windows (>99 percentile of Wstat) are represented by yellow crosses. Outlier windows that overlap with at least one candidate SNP from the multivariate analysis are represented by a yellow circle. 120 Figure S2.6. Representation of neutral PCA, with PC1 and PC2 on the left and PC1 and PC3 on the right. Samples are color-coded by population. Figure S2.7 Predicted turnover in neutral genetic composition as calculated by GDM (right) as consequence of gradients of the different environmental predictors (left). Locations with similar colors are expected to have similar genetic composition, based on the effects of the first three principal component gradients (biplot in middle). Circles represent sampled individuals and green line delimits Eurasian lynx current distributional range (Breitenmoser et al., 2015). 121 Lucena‐Perez et al., 2021). Individuals belonging to the Western and Southern clades of Eurasian lynx were selected, as the genetic distinction of these clades was corroborated by genomic data (Bazzicalupo et al., 2022; Lucena‐Perez et al., 2020) and they represent the two clades geographically closest to the possible contact point between the two lynx species (Sommer & Benecke, 2006). As the methodological approach relies on training a model based on simulated data under different introgression scenarios, we first performed demographic modeling between the two population pairs, with the aim of both unraveling the recent and past demographic histories of these populations with greater detail allowed by newly sequenced individuals and high-quality reference genomes, and creating a training dataset that was most similar to real data as possible. Then, with the genomic regions with signals of introgression identified in each species, we estimated the amount of introgression observable in each species’ genome, across different populations and genomic regions. We address the question of which of the remnant Eurasian lynx lineages is more closely related to the extinct population that coexisted with the Iberian lynx in Northern Iberia and Southern Europe. Finally, by looking at the genomic regions where introgression is happening, we assess what areas of the genome might be more permeable to introgression, to which extent introgression has affected genic and intergenic regions, and how introgression is affecting local patterns of genetic diversity in the different populations. METHODS Sampling, DNA extraction and sequencing Whole genome sequences of a total of 54 lynx individuals were analyzed in this study. Individuals from the Andujar population (LPA, n=22) of Iberian lynx (Lynx pardinus) were chosen to be analyzed alongside individuals belonging to the Western (WEL, n=20) and Southern (SEL, n=12) clades of Eurasian lynx (Figure 3.1). Some of the sequences used in this analysis were generated during previous genomic studies of both Iberian (Lucena‐Perez et al., 2021) and Eurasian lynx (Bazzicalupo et al., 2022; Lucena‐Perez et al., 2020). In this study we generated resequencing data for 10 additional individuals of Iberian lynx. DNA was extracted from tissue and blood samples using a combination of silica-coated paramagnetic beads (NucleoMag® Tissue, MACHEREY-NAGEL GmbH & Co. KG) and a classical phenol-chloroform protocol in LEM-EBD facilities (Seville, Spain). gDNA samples were then sent for paired-end library preparation and sequencing on Illumina NovaSeq X Plus platform at Novogene Ltd facilities (Cambridge, UK). In table S1 there are details regarding each sample’s origin, original sequencing study, sequencing technology used and average sequencing depth. In all cases, primary data analysis was carried out with the standard Illumina pipeline. Sequencing reads were processed using the software fastp (Chen et al., 2018) and quality control of processed reads was carried out using the software fastqc (https://www.bioinformatics.babraham.ac.uk/projects/fastqc). We enabled base correction 128 using overlapping paired sequences, trimmed adapters and poly-G sequences and selected a minimum length of 30 base pairs per read after trimming. Sequencing reads alignment, variant calling and initial filtering To avoid reference bias in our analyses, sequencing reads were aligned to the Lynx rufus reference genome (mLynRuf2.2, https://denovo.cnag.cat/lynx_rufus_data), using BWA-MEM (H. Li, 2013). After alignment, duplicate reads were marked with the MarkDuplicates function of Picard Tools (http://broadinstitute.github.io/picard), and INDELs were realigned using GATK RealignerTargetCreator and IndelRealigner (McKenna et al., 2010). We used the software Qualimap 2 (Okonechnikov et al., 2016) to assess the alignment quality and calculate average sequencing depth for each sample (Table S1). Variant detection and calling were performed using the software GATK HaplotypeCaller v4.2 (McKenna et al., 2010), which generated a Genome Variant Calling Format (GVCF) file for each of our individuals. These were then combined into a single GVCF of all of the individuals using the CombineGVCFs tool of GATK. Finally, the GenotypeGVCFs tool was used to call variants and genotypes of all of the individuals in the unified GVCF. We limited our analysis to the autosomal regions of the reference genome, which includes 18 distinct chromosomes, excluding any variants not falling in those. We identified low complexity regions and interspersed repeats in the reference genome using the combination of RepeatModeler and RepeatMasker softwares (http://www.repeatmasker.org). Any variant found within these regions, was removed from the VCF together with any INDEL and non-biallelic variant. Low quality variants were also filtered out following GATK's suggested thresholds of variant quality (i.e., FS ≥ 60, MQ ≤ 40), using a slightly stricter threshold for the Quality by Depth parameter (QD < 6 instead of the suggested QD < 2), after inspecting the genome-wide distribution of values of these statistics. We then separated our dataset into three distinct VCFs, one for each population, to perform population specific variant filtering. From each of the three population VCFs we filtered out any SNP with 15% or more missing data. Additionally, to avoid regions with multiple paralogs collapsed in the reference genome, mean read depth in consecutive 10 kbp windows along the genome of each individual was calculated using the software mosdepth (Pedersen & Quinlan, 2018). For each window we then calculated the sum of sequencing depth in each population and excluded any window with a total sequencing depth higher than 1.5 times the mode of all sequencing depth values in the population. We joined the filtered Iberian lynx VCF with each of the Eurasian lynx filtered VCFs, creating two distinct population-pair datasets for downstream analysis: the Iberian-Western Eurasian (LPA-WEL) pair, and the Iberian-Southern Eurasian (LPA-SEL) pair. 129 Demographic inference In order to reconstruct the joint demographic history of the Iberian and Eurasian lynx population pairs we ran the software GADMA2 (Noskova et al., 2023). GADMA2 has the advantage of providing an option for dynamic demographic reconstruction through its structure model specification. This provides a more data-driven and unsupervised approach to demographic modeling by using genetic algorithms for global and local optimization (Momigliano et al., 2021), to dynamically explore the likelihood space under different numbers of epochs before and after population splits, as well as different population size change functions. We further filtered our dataset for running GADMA2 in optimal conditions. Specifically, to reduce the bias introduced by selection on demographic inference, we used bedtools (Quinlan & Hall, 2010) to remove any SNP found within annotated genes. Additionally, we pruned the resulting SNP dataset to reduce the amount of linkage disequilibrium between SNPs and to remove missing data. The latter step allowed direct likelihood comparison between model runs using the full Site Frequency Spectrum (SFS). We did so using plink v1.9 (Purcell et al., 2007), with --thin 50 000 to remove SNPs at physical distances of less than 50 000 bp between each other and --max-missing 1 to remove any SNP with missing data. We ran GADMA2 in structure model mode, specifying a final model structure with two epochs before and three after the population split, using the length of the different time spans, the effective population sizes, the growth rates of the active populations, and the migration rates among them as model parameters (Figure S3.1). We set a minimum of 20 000 generations for the divergence time between lynx species, allowing the population size changes to either be exponential or instantaneous. The default parameter bounds were changed to be [0.001 - 50] and [0 - 10] times the ancestral population size for each effective population size and event time respectively. We selected moments as the engine for frequency spectrum (FS) simulation and model evaluation, and set the per base per generation mutation rate to 6e-9 (Abascal et al., 2016). Following the suggestions included in Noskova et al. (2023), we selected a hyperparameter configuration for the genetic algorithm that best performed with the moments engine: Size of generation = 10, N elitism = 3, P mutation = 0.2, P crossover = 0.3, P random = 0.2, Mean mutation strength = 0.776, Const for mutation strength = 1.302, Mean mutation rate = 0.273, Const for mutation rate = 1.475. With the intention of thoroughly exploring the likelihood landscape, a total of 500 distinct runs of the genetic algorithm were performed. Finally, we compared the likelihoods of the best model of each of the 500 runs, selecting the highest scoring models. Local optimization was then performed on the best scoring models, using 100 block bootstrap replicates, with a chunk size of 200 kbp, of the dataset including all of the SNPs before pruning. The local optimization algorithm was run using GADMA’s run_ls_on_boot_data script, which implements dadi’s optimize_log function, and consists in using the BFGS method to optimize log transformed 130 parameter values. Median and 95% confidence intervals of each parameter were calculated from the optimized bootstrap replicates. Phasing variants In order to perform genome scans for introgression we needed phased haplotype data, which we obtained from genotype data of each lynx population. We first split our data into three distinct VCFs, one for each of the three populations, and then applied a two-step approach. We first used WhatsHap v1.1 (Martin et al., 2016) to obtain phase sets (--tag=PS) from individual reads and population data. Next, we used SHAPEIT v.4.2.1 (Delaneau et al., 2019) to infer each sample’s haplotype from the phase set data. An expected error rate of 0.0001 and "10b,1p,1b,1p,1b,1p,1b,1p,10m" MCMC iterations were used as suggested by the SHAPEIT4 manual. Once phased, the variants from each Eurasian lynx population were joined with the LPA variants, creating two population-pair phased datasets for the introgression scans (phased LPA-WEL and phased LPA-SEL). Introgression scans Methodological approach To identify genomic windows with signals of introgression in either direction between population pairs we used the introNets package (https://github.com/SchriderLab/introNets; Ray et al., 2024) to implement a discriminator model trained to look at haplotype SNP data for a particular genomic window and distinguish between different introgression scenarios. The discriminator model is a Residual Network (ResNet34), a type of deep neural network architecture that was developed for general image classification tasks (Koonce, 2021), and was recently implemented to detect introgressed regions between closely related species with high accuracy (Ray et al., 2024). In order to format the genomic windows for the model to classify, we used genotypes observed in phased biallelic SNPs to generate three dimensional tensors (l × n × m), which would be handled by the model the same way it handles multilayered images. The first dimension of the tensor (l) is the number of populations in the alignment, the second dimension (n) is the number of phased haplotypes observed in each population (two for each sequenced diploid individual), and the third dimension (m) is the number of SNPs in the window. As we are analyzing introgression between pairs of populations, the first dimension of the tensor is always equal to 2. Although the number of observed phased haplotypes changed in each population (LPA: 44, WEL: 40, SEL: 24), the tensor needs to have the same number of haplotypes in both populations for the image to be correctly formatted. To solve this issue and to not lose information on any of the sequenced haplotypes, we selected the highest number of observed haplotypes in any population (44) as the second dimension and upsampled the populations with missing haplotypes, selecting observed haplotypes at random to appear multiple times in the input tensor (following Ray et al., 2024). The third dimension was set to 128, which is the number of SNPs included in each 131 genomic window we analyzed. This resulted in tensors of shape (2, 44, 128) for all pairs of populations analyzed. The value of a particular entry in the tensor is the genotype observed at a particular SNP, in a particular haplotype, in a particular population. Since the SNPs are biallelic, the resulting tensors have binary values in their entries, 0 for the reference allele and 1 for the alternative. The discriminator model was trained to look at a particular tensor representing genotypes of a genomic window and assign a probability to each of 4 distinct scenarios that might have generated the resulting alignment: ‘ab’ where introgression from population 1 is observed in population 2; ‘ba’ where introgression from population 2 is observed in population 1; ‘bi’ where we observe introgression from the other population in both population 1 and population 2; and ‘none’ when no introgression is observed. Simulation and training To train the model to perform the discrimination task we ran coalescent simulations under each of the three selected demographic models of each population pair for each introgression scenario, using a modified version of the software ms (Hudson, 2002) called msmodified, developed by Durvasula & Sankararaman (2019) and available from the ArchIE (https://github.com/sriramlab/ArchIE) and introNets packages. When generating data for the cases without introgression (‘none’ scenario) the demographic history was left unmodified. When simulating cases with introgression we added to the demographic model a mass migration event in the appropriate direction, from population 1 to population 2 for ‘ab’ cases, from population 2 to population 1 for ‘ba’ cases and in both directions for ‘bi’ cases. The mass migration event was simulated by moving a proportion of lineages from one population into the other at a given time. The proportion was randomly selected to be between 5% and 50% of the source population lineages, while the timing of the event was randomly selected between 1 and 5000 generations before present. Tensors of the appropriate shape (see above) were generated from the simulated data, making sure, based on the required introgression case, that introgressed haplotypes were present in the receiving population; this information is recorded by msmodified. The final set of simulations for training consisted of 12000 simulations of each introgression scenario, simulated under each of the demographic models reconstructed for the population pair, for a total of 144000 training tensors. The total simulated dataset was split into training and validation sets, with the latter comprising 5% of the simulations for each introgression scenario. We performed training using the categorical cross entropy as loss function and the Adam optimizer with default settings (Kingma & Ba, 2014), as described in Ray et al. (2024). Maximum training time was set to 100 epochs, with 1500 training steps in each epoch, and a maximum of 10 consecutive epochs without validation loss decreasing before early termination of training. Throughout training, the set of weights that produced the lowest validation loss up to that point was recorded, and upon termination these weights were used as the final classification model. Two distinct models were trained, one for each population pair. 132 Model performance The performance of the model was evaluated using 100 additional simulations under each introgression scenario with each demographic model. The model’s output was transformed into a probability score for each introgression scenario using the softmax function. During model evaluation, we assigned a particular window to a particular introgression category using the following criteria: if the sum of the probabilities of ‘ab’ and ‘bi’ was higher than a probability threshold p, we assigned the window to the ‘ab’ category; if the sum of the probabilities for ‘ba’ and ‘bi’ was higher than p then the window would be assigned to ‘ba’; if a window is assigned to both ‘ab’ and ‘ba’, then it’s prediction is changed to ‘bi’; if the probability threshold was not exceeded for any of these cases, then the window was not assigned to any introgression category and was instead categorized as ‘none’. We evaluated model performance using five increasing values of p (0.75, 0.8, 0.85, 0.9, 0.95), comparing model precision and recall for each of the four introgression categories. Model precision is defined as the proportion of cases where a category was called when the data was really simulated under that introgression scenario, without considering the number of simulations under that scenario that were not correctly classified. Model recall on the other hand is defined as the proportion of cases where a category was correctly called compared to the total number of simulations in that category, without considering the number of simulations incorrectly assigned to that category. Introgressed windows and downstream analyses Finally, the model trained for each of the two population-pair datasets was run on phased data from the corresponding dataset. Overlapping sliding windows along the genome, with a window length of 128 SNPs and a step size of 64 SNPs were transformed into input tensors of the desired shape, generating a total of 79184 LPA-WEL genomic windows and 77565 LPA-SEL genomic windows for the introgression scan. We assigned genomic windows to different introgression categories using a criteria similar to the one used during model evaluation, by choosing an initial value of p of 0.9. We also assigned an introgression category to those windows with a probability of introgression (‘ab’ or ‘ba’ plus ‘bi’) greater than 0.7, if they were adjacent to another window with at least 0.9 probability for the same category. Consecutive windows assigned to the same introgression class were then merged, creating four sets of introgressed regions: (1) windows of introgression from LPA to WEL; (2) windows of introgression from LPA to SEL; (3) windows of introgression from WEL to LPA; (4) windows of introgression from SEL to LPA. Additionally, the final dataset of genomic windows of introgression from Eurasian lynx into LPA was generated by merging introgressed windows from SEL and WEL into LPA. We explored the effects of introgression on genetic and functional diversity , as well as the distribution of chromosomal locations of introgressed windows in each of the recipient populations (LPA, WEL and SEL). Specifically, we assessed 1) the impact of introgression 133 on nucleotide diversity (π), 2) whether introgressed regions were enriched or depleted of genes, and 3) whether introgressed regions where biased towards or away from the telomeres, by comparing each of these features of introgressed windows to a null model. The null model for a given population was composed of 100 datasets of windows of the same length as the introgressed ones, but placed at random locations along the genome using bedtools (Quinlan & Hall, 2010). We then calculated π, distance to the closest telomere and amount of overlap with annotated genes of each window in the random dataset and compared the differences in the distributions of these values in introgressed versus random windows. Effect sizes of the observed differences were obtained from log transformed data, to obtain normally distributed values, using Cohen’s D (Cd) in the R package lsr (Navarro, 2015). RESULTS We constructed a dataset of 22 Iberian lynx (LPA), 20 Western Eurasian lynx (WEL), and 12 Southern Eurasian lynx genomes (Figure 3.1), including 10 Iberian lynx genomes newly sequenced for this study (Methods). Our sequencing, reads alignment, and variant calling efforts resulted in an initial set of 33 660 575 raw autosomal variants. We then removed INDELs, any non-biallelic SNP, any SNP not meeting the quality thresholds and any SNP falling within low-complexity or repetitive regions of the reference genome. This left a set of 6 495 299 high-quality biallelic SNPs in our VCF file. The 15% missing data filter applied to each population resulted in a loss of around 3.9% of SNPs in Iberian lynx, 4.2% in Western and 1.8% in the Southern Eurasian lynx. The lower missing rate in Southern Eurasian lynx is probably given by the relatively high sequencing depth and quality obtained for all of the samples we have from this population, while other populations contain a mix of individuals with both relatively high and relatively low sequencing depth. After removing SNPs from possible collapsed paralogs in each population we generated datasets for each pair of Eurasian and Iberian lynx populations, obtaining 4 777 474 SNPs that were segregating in the Iberian-Western Eurasian (LPA-WEL) pair and 4 753 415 SNPs in the Iberian-Southern Eurasian (LPA-SEL) pair. Demographic inference To train a classifier to detect introgression, we first had to reconstruct the demographic history of each population pair to be used during training simulations. We sought to minimize any bias caused by natural selection by removing SNPs potentially affected by direct or linked selection (Johri et al., 2021; Schrider et al., 2016). After removing SNPs within genes, removing SNPs with missing data and pruning SNPs at less than 50k bp distance between each other we obtained a total of 31 628 SNPs in the LPA-WEL dataset, and 31 882 SNPs in the LPA-SEL dataset. We then used the joint site frequency spectra from these SNPs to estimate demographic histories with GAMDA2. We modeled the demographic history of LPA-WEL and LPA-EEL five hundred times for each population pair. From these, we selected each pair’s three highest likelihood reconstructions and ran local optimizations using 134 100 bootstrap datasets. The following parameters were optimized in each model: (1) the length in generations of each of the five epochs in the model, two before and three after divergence between populations; (2) effective population sizes of each population at the start and at the end of each epoch (these were equal if size was constant in the model, but differed if size changed exponentially during that epoch); (3) migration rates, calculated as the probability of choosing a parent in one population (destination) from the other population (source), in each direction between the populations during each post-divergence epoch. Divergence times between Iberian and Eurasian lynx are consistent in the reconstruction made with both population pairs, with both having one model placing the divergence around 150 kya, one model around 600 kya and one model around 800 kya (Figure 3.2). The recent effective population size of the LPA population follows a similar trend in all models, both when reconstructed in conjunction with WEL and with SEL (Figure 3.2). In all cases LPA experiences a decline, although this decline is inferred to be sharper (down to approximately one third of the size before the contraction) and more recent (between 3 kya and 20 kya) in models with WEL, and less intense (around a 50% reduction), and more prolonged (occurring over a period of between 30 kya and 50 kya), in models with SEL. Regardless of intensity and timing, all models recover a decline of the LPA population from around 5k to 7k individuals in ancient times to less than 3k individuals in more recent times (Figure 3.2). Differences in the population trajectories of WEL and SEL are evident in all model reconstructions (Figure 3.2). Although the effective population size of WEL in the most recent epochs is consistently estimated to be around 4k to 6k individuals and to have experienced a significant recent decline in all three LPA-WEL models, although the decline’s timing and intensity varies among models, decreasing from a starting size of 10k to 20k individuals around 3 kya to 20 kya. The trajectory through time of the SEL population is much more variable across models, slightly increasing in the last ten thousand years in one model, and staying constant or slightly decreasing during the last fifty thousand years in the other two models. In all models, the most recent population size is double that of WEL, with around 10k individuals. Migration rates between populations are more symmetrical and higher for all LPA-WEL models than between the LPA-SEL models (Figure 3.2). In two of the LPA-WEL models ancient migration rates were higher than most recent, with one of the models having a period of no migration between the two populations from 20 kya to 6 kya. The LPA-WEL model 1, which presents the lowest recovered divergence time among the three, shows a different pattern. In this model, migration from LPA to WEL around 2 kya has a dramatic increase, rising to ten times the levels recovered for the other models in recent times, followed by a great decline to around 4% of the other model’s migration rates. Migration from WEL to LPA 135 on the other hand is increasing slightly in that same time period. Reconstructions of ancient migration rates between LPA and SEL are more consistent among models. Ancient rates from LPA to SEL are always lower than what is estimated for recent times. The opposite is inferred for migration rates from SEL to LPA, which are always declining from ancient to most recent times in all models. Recent migration rates between LPA and WEL are generally higher than between LPA and SEL. 136 Figure 3.2 Summary of demographic modeling performed by running GADMA2’s structure model with three post-divergence epochs. Out of five hundred runs, the three highest scoring models are represented for each population pair, LPA-WEL on the left and LPA-SEL on the right. Each dot or line represents the optimized value of the parameter for one of one hundred bootstrap replicate datasets. Divergence times are represented by dots, color coded by model. Lines follow the same color code and represent changes in effective population size and directional migration rates in the last million years, across the three post-divergence epoch of all the models. Time is transformed to years assuming a generation time of 5 years (Lucena-Perez et al., 2018). Migration rates are defined in forward time as the probability of choosing a parent in the sink population from the source population. All values, except for the divergence times on top, are represented on a logarithmic scale. 137 regions of different recombination rate (Schumer et al., 2018), but is also true at the genome level in species with higher or lower average recombination rates (Veller et al., 2023), and in small versus large chromosomes (Edelman et al., 2019), with the former generally having higher per base pair recombination rates. We observe contrasting patterns of the interplay between selection and recombination on introgression in the genomes of Iberian and Eurasian lynx. On the one hand, in the introgressed windows of the three lynx populations we found similar amounts of gene sequences as random subsamples of the genomes. In principle this might indicate that, at least theoretically, selection in the receiving populations didn’t strongly favor the elimination of gene sequences coming from the other species. On the other hand, the effects of recombination on the probability of introgression can be observed in all three lynx populations: introgressed windows are more likely to be found closer to the telomeres, where recombination is higher (G. Li, Hillier, et al., 2016), and more introgression is observed in shorter chromosomes than longer ones in all populations. More precise and detailed inferences regarding variation in the strength of selection and the recombination rate along the genome might help clarify how the interaction between these two factors has shaped the genomic introgression landscape of lynxes. Demographic history and migration are also expected to have important consequences for the levels of genetic diversity detected in populations. In our results, genetic diversity in the regions with introgression of all lynx populations examined is higher than diversity measured at random locations along the genome. This positive correlation between introgression and diversity is expected, as divergent introgressing alleles increase the average distance between haplotypes in the region. This injection of diversity through migration, already observed by Lucena-Perez et al. (2024), might have represented an important offset to the loss of genetic diversity caused by population declines that happened throughout the history of the LPA and WEL populations. Although assisted gene flow between populations of different species is usually discouraged in conservation programs (Aitken & Whitlock, 2013; Frankham et al., 2011), natural migration probably played an important role in maintaining the genetic diversity of lynx populations through time. Conservation measures might allow higher flexibility with regards to source populations for assisted gene flow in circumstances similar to those faced recently by the Iberian lynx. When no additional conspecific populations exist, a recently divergent and not reproductively isolated species might represent a natural reservoir of novel genetic variation. Still, caution and constant genetic monitoring are warranted to avoid the genomic extinction risks associated with outbreeding depression and genetic swamping (Todesco et al., 2016). In conclusion, in this study we applied a novel deep learning approach to quantify and, for the first time, localize regions with signals of interspecific introgression in the genomes of three lynx populations. By comparing the amount of genome with introgression in two distinct populations of Eurasian lynx, we found evidence that the Western Eurasian lynx clade could 144 be the most closely related to the extinct Eurasian population that occupied the hybrid zone with the Iberian lynx. Although fewer, some signals of introgression are also observable in the Southern Eurasian lynx clade. The lower introgression in this population might be the result of biogeographical barriers isolating it more from the hybrid zone, or a stronger purging of introgressed material driven by its higher population size. We were also able to observe how recombination and selection influence the chromosomes and regions where introgression is more likely to happen, with higher recombination regions and smaller chromosomes being the most introgressed, and gene sequences just as likely as intergenic regions. Introgression is also boosting genetic diversity in the affected regions of all species, suggesting that assisted gene flow might be a viable solution for extreme conservation cases such as the one faced by the Iberian lynx, previous to its most recent conservation-driven recovery. REFERENCES Abascal, F., Corvelo, A., Cruz, F., Villanueva-Cañas, J. L., Vlasova, A., Marcet-Houben, M., Martínez-Cruz, B., Cheng, J. Y., Prieto, P., Quesada, V., Quilez, J., Li, G., García, F., Rubio-Camarillo, M., Frias, L., Ribeca, P., Capella-Gutiérrez, S., Rodríguez, J. M., Câmara, F., … Godoy, J. A. (2016). Extreme genomic erosion after recurrent demographic bottlenecks in the highly endangered Iberian lynx. Genome Biology, 17(1), 251. https://doi.org/10.1186/s13059-016-1090-1 Aitken, S. N., & Whitlock, M. C. (2013). Assisted Gene Flow to Facilitate Local Adaptation to Climate Change. Annual Review of Ecology, Evolution, and Systematics, 44(1), 367–388. https://doi.org/10.1146/annurev-ecolsys-110512-135747 Barton, N. H., & Hewitt, G. M. (1989). Adaptation, speciation and hybrid zones. Nature, 341(6242), 497–503. https://doi.org/10.1038/341497a0 Bazzicalupo, E., Lucena‐Perez, M., Kleinman‐Ruiz, D., Pavlov, A., Trajçe, A., Hoxha, B., Sanaja, B., Gurielidze, Z., Kerdikoshvili, N., Mamuchadze, J., Yarovenko, Y. A., Akkiev, M. I., Ratkiewicz, M., Saveljev, A. P., Melovski, D., Gavashelishvili, A., Schmidt, K., & Godoy, J. A. (2022). History, demography and genetic status of Balkan and Caucasian Lynx lynx (Linnaeus, 1758) populations revealed by genome‐wide variation. Diversity and Distributions, 28(1), 65–82. https://doi.org/10.1111/ddi.13439 Bierne, N., Lenormand, T., Bonhomme, F., & David, P. (2002). Deleterious mutations in a hybrid zone: Can mutational load decrease the barrier to gene flow? Genetical Research, 80(3), 197–204. https://doi.org/10.1017/S001667230200592X Breitenmoser, U., Breitenmoser-Würsten, C., Lanz, T., von Arx, M., & Antonevich, A. (2015). Lynx lynx (errata version published in 2017). IUCN Red List of Threatened Species. https://www.iucnredlist.org/en Brisbin, A., Bryc, K., Byrnes, J., Zakharia, F., Omberg, L., Degenhardt, J., Reynolds, A., Ostrer, H., Mezey, J. G., & Bustamante, C. D. (2012). PCAdmix: Principal Components-Based Assignment of Ancestry Along Each Chromosome in Individuals with Admixed Ancestry from Two or More Populations. Human Biology, 84(4), 145 343–364. https://doi.org/10.3378/027.084.0401 Casas-Marce, M., Marmesat, E., Soriano, L., Martínez-Cruz, B., Lucena-Perez, M., Nocete, F., Rodríguez-Hidalgo, A., Canals, A., Nadal, J., Detry, C., Bernáldez-Sánchez, E., Fernández-Rodríguez, C., Pérez-Ripoll, M., Stiller, M., Hofreiter, M., Rodríguez, A., Revilla, E., Delibes, M., & Godoy, J. A. (2017). Spatiotemporal Dynamics of Genetic Variation in the Iberian Lynx along Its Path to Extinction Reconstructed with Ancient DNA. Molecular Biology and Evolution, 34(11), 2893–2907. https://doi.org/10.1093/molbev/msx222 Chan, J., Perrone, V., Spence, J. P., Jenkins, P. A., Mathieson, S., & Song, Y. S. (2018). A Likelihood-Free Inference Framework for Population Genetic Data using Exchangeable Neural Networks. Advances in Neural Information Processing Systems, 31, 8594–8605. Chen, S., Zhou, Y., Chen, Y., & Gu, J. (2018). fastp: An ultra-fast all-in-one FASTQ preprocessor. Bioinformatics, 34(17), i884–i890. https://doi.org/10.1093/bioinformatics/bty560 Delaneau, O., Zagury, J.-F., Robinson, M. R., Marchini, J. L., & Dermitzakis, E. T. (2019). Accurate, scalable and integrative haplotype estimation. Nature Communications, 10(1), 5436. https://doi.org/10.1038/s41467-019-13225-y DeWoody, J. A., Harder, A. M., Mathur, S., & Willoughby, J. R. (2021). The long‐standing significance of genetic diversity in conservation. Molecular Ecology, 30(17), 4147–4154. https://doi.org/10.1111/mec.16051 Durvasula, A., & Sankararaman, S. (2019). A statistical model for reference-free inference of archaic local ancestry. PLOS Genetics, 15(5), e1008175. https://doi.org/10.1371/journal.pgen.1008175 Edelman, N. B., Frandsen, P. B., Miyagi, M., Clavijo, B., Davey, J., Dikow, R. B., García-Accinelli, G., Van Belleghem, S. M., Patterson, N., Neafsey, D. E., Challis, R., Kumar, S., Moreira, G. R. P., Salazar, C., Chouteau, M., Counterman, B. A., Papa, R., Blaxter, M., Reed, R. D., … Mallet, J. (2019). Genomic architecture and introgression shape a butterfly radiation. Science, 366(6465), 594–599. https://doi.org/10.1126/science.aaw2090 Edelman, N. B., & Mallet, J. (2021). Prevalence and Adaptive Impact of Introgression. Annual Review of Genetics, 55(1), 265–283. https://doi.org/10.1146/annurev-genet-021821-020805 Flagel, L., Brandvain, Y., & Schrider, D. R. (2019). The Unreasonable Effectiveness of Convolutional Neural Networks in Population Genetic Inference. Molecular Biology and Evolution, 36(2), 220–238. https://doi.org/10.1093/molbev/msy224 Frankham, R., Ballou, J. D., Eldridge, M. D. B., Lacy, R. C., Ralls, K., Dudash, M. R., & Fenster, C. B. (2011). Predicting the Probability of Outbreeding Depression: Predicting Outbreeding Depression. Conservation Biology, 25(3), 465–475. https://doi.org/10.1111/j.1523-1739.2011.01662.x Geneva, A. J., Muirhead, C. A., Kingan, S. B., & Garrigan, D. (2015). A New Method to Scan Genomes for Introgression in a Secondary Contact Model. PLOS ONE, 10(4), 146 e0118621. https://doi.org/10.1371/journal.pone.0118621 Gladieux, P., Ropars, J., Badouin, H., Branca, A., Aguileta, G., De Vienne, D. M., Rodríguez De La Vega, R. C., Branco, S., & Giraud, T. (2014). Fungal evolutionary genomics provides insight into the mechanisms of adaptive divergence in eukaryotes. Molecular Ecology, 23(4), 753–773. https://doi.org/10.1111/mec.12631 Goulet, B. E., Roda, F., & Hopkins, R. (2017). Hybridization in Plants: Old Ideas, New Techniques. Plant Physiology, 173(1), 65–78. https://doi.org/10.1104/pp.16.01340 Harris, K., & Nielsen, R. (2016). The Genetic Cost of Neanderthal Introgression. Genetics, 203(2), 881–891. https://doi.org/10.1534/genetics.116.186890 Hibbins, M. S., & Hahn, M. W. (2022). Phylogenomic approaches to detecting and characterizing introgression. Genetics, 220(2), iyab173. https://doi.org/10.1093/genetics/iyab173 Hudson, R. R. (2002). Generating samples under a Wright–Fisher neutral model of genetic variation. Bioinformatics, 18(2), 337–338. https://doi.org/10.1093/bioinformatics/18.2.337 Ingvarsson, P. K., & Whitlock, M. C. (2000). Heterosis increases the effective migration rate. Proceedings of the Royal Society of London. Series B: Biological Sciences, 267(1450), 1321–1326. https://doi.org/10.1098/rspb.2000.1145 Johnson, W. E. (2004). Phylogenetic and Phylogeographic Analysis of Iberian Lynx Populations. Journal of Heredity, 95(1), 19–28. https://doi.org/10.1093/jhered/esh006 Johnson, W. E., Eizirik, E., Pecon-Slattery, J., Murphy, W. J., Antunes, A., Teeling, E., & O’Brien, S. J. (2006). The Late Miocene Radiation of Modern Felidae: A Genetic Assessment. Science, 311(5757), 73–77. https://doi.org/10.1126/science.1122277 Johri, P., Riall, K., Becher, H., Excoffier, L., Charlesworth, B., & Jensen, J. D. (2021). The Impact of Purifying and Background Selection on the Inference of Population History: Problems and Prospects. Molecular Biology and Evolution, 38(7), 2986–3003. https://doi.org/10.1093/molbev/msab050 Juric, I., Aeschbacher, S., & Coop, G. (2016). The Strength of Selection against Neanderthal Introgression. PLOS Genetics, 12(11), e1006340. https://doi.org/10.1371/journal.pgen.1006340 Kim, B. Y., Huber, C. D., & Lohmueller, K. E. (2018). Deleterious variation shapes the genomic landscape of introgression. PLOS Genetics, 14(10), e1007741. https://doi.org/10.1371/journal.pgen.1007741 Kingma, D. P., & Ba, J. (2014). Adam: A Method for Stochastic Optimization (Version 9). arXiv. https://doi.org/10.48550/ARXIV.1412.6980 Kleinman-Ruiz, D., Lucena-Perez, M., Villanueva, B., Fernández, J., Saveljev, A. P., Ratkiewicz, M., Schmidt, K., Galtier, N., García-Dorado, A., & Godoy, J. A. (2022). Purging of deleterious burden in the endangered Iberian lynx. Proceedings of the National Academy of Sciences, 119(11), e2110614119. https://doi.org/10.1073/pnas.2110614119 147 Koonce, B. (2021). ResNet 34. In B. Koonce, Convolutional Neural Networks with Swift for Tensorflow (pp. 51–61). Apress. https://doi.org/10.1007/978-1-4842-6168-2_5 Kurtén, B., & Granqvist, E. (1987). Fossil pardel lynx (Lynx pardina spelaea Boule) from a cave in southern France. Annales Zoologici Fennici, 24(1), 39–43. Kyriazis, C. C., Wayne, R. K., & Lohmueller, K. E. (2021). Strongly deleterious mutations are a primary determinant of extinction risk due to inbreeding depression. Evolution Letters, 5(1), 33–47. https://doi.org/10.1002/evl3.209 Lawson, D. J., Hellenthal, G., Myers, S., & Falush, D. (2012). Inference of Population Structure using Dense Haplotype Data. PLoS Genetics, 8(1), e1002453. https://doi.org/10.1371/journal.pgen.1002453 Li, G., Davis, B. W., Eizirik, E., & Murphy, W. J. (2016). Phylogenomic evidence for ancient hybridization in the genomes of living cats (Felidae). Genome Research, 26(1), 1–11. https://doi.org/10.1101/gr.186668.114 Li, G., Figueiró, H. V., Eizirik, E., & Murphy, W. J. (2019). Recombination-Aware Phylogenomics Reveals the Structured Genomic Landscape of Hybridizing Cat Species. Molecular Biology and Evolution, 36(10), 2111–2126. https://doi.org/10.1093/molbev/msz139 Li, G., Hillier, L. W., Grahn, R. A., Zimin, A. V., David, V. A., Menotti-Raymond, M., Middleton, R., Hannah, S., Hendrickson, S., Makunin, A., O’Brien, S. J., Minx, P., Wilson, R. K., Lyons, L. A., Warren, W. C., & Murphy, W. J. (2016). A High-Resolution SNP Array-Based Linkage Map Anchors a New Domestic Cat Draft Genome Assembly and Provides Detailed Patterns of Recombination. G3&#58; Genes|Genomes|Genetics, 6(6), 1607–1616. https://doi.org/10.1534/g3.116.028746 Li, H. (2013). Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv:1303.3997 [q-Bio]. http://arxiv.org/abs/1303.3997 Lucena-Perez, M., Bazzicalupo, E., Paijmans, J., Kleinman-Ruiz, D., Dalén, L., Hofreiter, M., Delibes, M., Clavero, M., & Godoy, J. A. (2022). Ancient genome provides insights into the history of Eurasian lynx in Iberia and Western Europe. Quaternary Science Reviews, 285, 107518. https://doi.org/10.1016/j.quascirev.2022.107518 Lucena‐Perez, M., Kleinman‐Ruiz, D., Marmesat, E., Saveljev, A. P., Schmidt, K., & Godoy, J. A. (2021). Bottleneck‐associated changes in the genomic landscape of genetic diversity in wild lynx populations. Evolutionary Applications, 14(11), 2664–2679. https://doi.org/10.1111/eva.13302 Lucena‐Perez, M., Marmesat, E., Kleinman‐Ruiz, D., Martínez‐Cruz, B., Węcek, K., Saveljev, A. P., Seryodkin, I. V., Okhlopkov, I., Dvornikov, M. G., Ozolins, J., Galsandorj, N., Paunovic, M., Ratkiewicz, M., Schmidt, K., & Godoy, J. A. (2020). Genomic patterns in the widespread Eurasian lynx shaped by Late Quaternary climatic fluctuations and anthropogenic impacts. Molecular Ecology, 29(4), 812–828. https://doi.org/10.1111/mec.15366 Lucena-Perez, M., Paijmans, J. L. A., Nocete, F., Nadal, J., Detry, C., Dalén, L., Hofreiter, M., Barlow, A., & Godoy, J. A. (2024). Recent increase in species-wide diversity after interspecies introgression in the highly endangered Iberian lynx. Nature Ecology & 148 Evolution, 8(2), 282–292. https://doi.org/10.1038/s41559-023-02267-7 Mallet, J. (2005). Hybridization as an invasion of the genome. Trends in Ecology & Evolution, 20(5), 229–237. https://doi.org/10.1016/j.tree.2005.02.010 Martin, M., Patterson, M., Garg, S., O Fischer, S., Pisanti, N., Klau, G. W., Schöenhuth, A., & Marschall, T. (2016). WhatsHap: Fast and accurate read-based phasing [Preprint]. Bioinformatics. https://doi.org/10.1101/085050 Matute, D. R., Butler, I. A., Turissini, D. A., & Coyne, J. A. (2010). A Test of the Snowball Theory for the Rate of Evolution of Hybrid Incompatibilities. Science, 329(5998), 1518–1521. https://doi.org/10.1126/science.1193440 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 Research, 20(9), 1297–1303. https://doi.org/10.1101/gr.107524.110 Mecozzi, B., Sardella, R., Boscaini, A., Cherin, M., Costeur, L., Madurell-Malapeira, J., Pavia, M., Profico, A., & Iurino, D. A. (2021). The tale of a short-tailed cat: New outstanding Late Pleistocene fossils of Lynx pardinus from southern Italy. Quaternary Science Reviews, 262, 106840. https://doi.org/10.1016/j.quascirev.2021.106840 Mengüllüoğlu, D., Ambarlı, H., Barlow, A., Paijmans, J. L. A., Sayar, A. O., Emir, H., Kandemir, İ., Hofer, H., Fickel, J., & Förster, D. W. (2021). Mitogenome Phylogeny Including Data from Additional Subspecies Provides New Insights into the Historical Biogeography of the Eurasian lynx Lynx lynx. Genes, 12(8), 1216. https://doi.org/10.3390/genes12081216 Momigliano, P., Florin, A.-B., & Merilä, J. (2021). Biases in demographic modelling affect our understanding of recent divergence. Molecular Biology and Evolution, msab047. https://doi.org/10.1093/molbev/msab047 Moran, B. M., Payne, C., Langdon, Q., Powell, D. L., Brandvain, Y., & Schumer, M. (2021). The genomic consequences of hybridization. eLife, 10, e69016. https://doi.org/10.7554/eLife.69016 Nadeau, N. J., Ruiz, M., Salazar, P., Counterman, B., Medina, J. A., Ortiz-Zuazaga, H., Morrison, A., McMillan, W. O., Jiggins, C. D., & Papa, R. (2014). Population genomics of parallel hybrid zones in the mimetic butterflies, H. melpomene and H. erato. Genome Research, 24(8), 1316–1333. https://doi.org/10.1101/gr.169292.113 Navarro, D. J. (2015). Learning statistics with R: A tutorial for psychology students and other beginners. (Version 0.6) [Computer software]. https://learningstatisticswithr.com Noskova, E., Abramov, N., Iliutkin, S., Sidorin, A., Dobrynin, P., & Ulyantsev, V. I. (2023). GADMA2: More efficient and flexible demographic inference from genetic data. GigaScience, 12, giad059. https://doi.org/10.1093/gigascience/giad059 Okonechnikov, K., Conesa, A., & García-Alcalde, F. (2016). Qualimap 2: Advanced multi-sample quality control for high-throughput sequencing data. Bioinformatics, 32(2), 292–294. https://doi.org/10.1093/bioinformatics/btv566 149 Pedersen, B. S., & Quinlan, A. R. (2018). Mosdepth: Quick coverage calculation for genomes and exomes. Bioinformatics, 34(5), 867–868. https://doi.org/10.1093/bioinformatics/btx699 Purcell, S., Neale, B., Todd-Brown, K., Thomas, L., Ferreira, M. A. R., Bender, D., Maller, J., Sklar, P., de Bakker, P. I. W., Daly, M. J., & Sham, P. C. (2007). PLINK: A Tool Set for Whole-Genome Association and Population-Based Linkage Analyses. The American Journal of Human Genetics, 81(3), 559–575. https://doi.org/10.1086/519795 Quinlan, A. R., & Hall, I. M. (2010). BEDTools: A flexible suite of utilities for comparing genomic features. Bioinformatics, 26(6), 841–842. https://doi.org/10.1093/bioinformatics/btq033 Ray, D. D., Flagel, L., & Schrider, D. R. (2024). IntroUNET: Identifying introgressed alleles via semantic segmentation. PLOS Genetics, 20(2), e1010657. https://doi.org/10.1371/journal.pgen.1010657 Rodriguez, A. (2024). Lynx pardinus. IUCN Red List of Threatened Species. https://www.iucnredlist.org/species/12520/218695618 Rodríguez, A., & Delibes, M. (2002). Internal structure and patterns of contraction in the geographic range of the Iberian lynx. Ecography, 25(3), 314–328. https://doi.org/10.1034/j.1600-0587.2002.250308.x Rodríguez‐Varela, R., García, N., Nores, C., Álvarez‐Lao, D., Barnett, R., Arsuaga, J. L., & Valdiosera, C. (2016). Ancient DNA reveals past existence of Eurasian lynx in Spain. Journal of Zoology, 298(2), 94–102. https://doi.org/10.1111/jzo.12289 Rodríguez-Varela, R., Tagliacozzo, A., Ureña, I., García, N., Crégut-Bonnoure, E., Mannino, M. A., Arsuaga, J. L., & Valdiosera, C. (2015). Ancient DNA evidence of Iberian lynx palaeoendemism. Quaternary Science Reviews, 112, 172–180. https://doi.org/10.1016/j.quascirev.2015.01.009 Rosenzweig, B. K., Pease, J. B., Besansky, N. J., & Hahn, M. W. (2016). Powerful methods for detecting introgressed regions from population genomic data. Molecular Ecology, 25(11), 2387–2397. https://doi.org/10.1111/mec.13610 Schardl, C. L., & Craven, K. D. (2003). Interspecific hybridization in plant‐associated fungi and oomycetes: A review. Molecular Ecology, 12(11), 2861–2873. https://doi.org/10.1046/j.1365-294X.2003.01965.x Schrider, D. R., Ayroles, J., Matute, D. R., & Kern, A. D. (2018). Supervised machine learning reveals introgressed loci in the genomes of Drosophila simulans and D. sechellia. PLOS Genetics, 14(4), e1007341. https://doi.org/10.1371/journal.pgen.1007341 Schrider, D. R., Shanku, A. G., & Kern, A. D. (2016). Effects of Linked Selective Sweeps on Demographic Inference and Model Selection. Genetics, 204(3), 1207–1223. https://doi.org/10.1534/genetics.116.190223 Schumer, M., Xu, C., Powell, D. L., Durvasula, A., Skov, L., Holland, C., Blazier, J. C., Sankararaman, S., Andolfatto, P., Rosenthal, G. G., & Przeworski, M. (2018). Natural 150 selection interacts with recombination to shape the evolution of hybrid genomes. Science, 360(6389), 656–660. https://doi.org/10.1126/science.aar3684 Smith, C. C. R., Tittes, S., Ralph, P. L., & Kern, A. D. (2023). Dispersal inference from population genetic variation using a convolutional neural network. GENETICS, 224(2), iyad068. https://doi.org/10.1093/genetics/iyad068 Sommer, R. S., & Benecke, N. (2006). Late Pleistocene and Holocene development of the felid fauna (Felidae) of Europe: A review. Journal of Zoology, 269(1), 7–19. https://doi.org/10.1111/j.1469-7998.2005.00040.x Suvorov, A., Kim, B. Y., Wang, J., Armstrong, E. E., Peede, D., D’Agostino, E. R. R., Price, D. K., Waddell, P. J., Lang, M., Courtier-Orgogozo, V., David, J. R., Petrov, D., Matute, D. R., Schrider, D. R., & Comeault, A. A. (2022). Widespread introgression across a phylogeny of 155 Drosophila genomes. Current Biology, 32(1), 111-123.e5. https://doi.org/10.1016/j.cub.2021.10.052 Teixeira, J. C., & Huber, C. D. (2021). The inflated significance of neutral genetic diversity in conservation genetics. Proceedings of the National Academy of Sciences, 118(10), e2015096118. https://doi.org/10.1073/pnas.2015096118 Todesco, M., Pascual, M. A., Owens, G. L., Ostevik, K. L., Moyers, B. T., Hübner, S., Heredia, S. M., Hahn, M. A., Caseys, C., Bock, D. G., & Rieseberg, L. H. (2016). Hybridization and extinction. Evolutionary Applications, 9(7), 892–908. https://doi.org/10.1111/eva.12367 Vanderpool, D., Minh, B. Q., Lanfear, R., Hughes, D., Murali, S., Harris, R. A., Raveendran, M., Muzny, D. M., Hibbins, M. S., Williamson, R. J., Gibbs, R. A., Worley, K. C., Rogers, J., & Hahn, M. W. (2020). Primate phylogenomics uncovers multiple rapid radiations and ancient interspecific introgression. PLOS Biology, 18(12), e3000954. https://doi.org/10.1371/journal.pbio.3000954 Veller, C., Edelman, N. B., Muralidhar, P., & Nowak, M. A. (2023). Recombination and selection against introgressed DNA. Evolution, 77(4), 1131–1144. https://doi.org/10.1093/evolut/qpad021 Wang, L., Beissinger, T. M., Lorant, A., Ross-Ibarra, C., Ross-Ibarra, J., & Hufford, M. B. (2017). The interplay of demography and selection during maize domestication and expansion. Genome Biology, 18(1), 215. https://doi.org/10.1186/s13059-017-1346-4 Whitlock, M. C., Ingvarsson, P. K., & Hatfield, T. (2000). Local drift load and the heterosis of interconnected populations. Heredity, 84(4), 452–457. https://doi.org/10.1046/j.1365-2540.2000.00693.x 151 APPENDIX Table S3.1. Sampled individuals information. Sample’s name, species and population, publication in which it was first sequenced, Illumina platform used and average whole genome sequencing depth are provided. SAMPLE SPECIES POPULATION PUBLICATION PLATFORM DEPTH c_lp_sm_013 4 Lynx pardinus Iberian lynx (LPA) Kleinman-Ruiz et al. 2022 HiSeq2000,v3 6X c_lp_sm_015 5 Lynx pardinus Iberian lynx (LPA) Kleinman-Ruiz et al. 2022 HiSeq2000,v3 5X c_lp_sm_015 6 Lynx pardinus Iberian lynx (LPA) Kleinman-Ruiz et al. 2022 HiSeq2000,v3 6X c_lp_sm_016 1 Lynx pardinus Iberian lynx (LPA) Kleinman-Ruiz et al. 2022 HiSeq2000,v3 5X c_lp_sm_020 6 Lynx pardinus Iberian lynx (LPA) Kleinman-Ruiz et al. 2022 HiSeq2000,v3 5X c_lp_sm_020 8 Lynx pardinus Iberian lynx (LPA) Kleinman-Ruiz et al. 2022 HiSeq2000,v3 5X c_lp_sm_021 3 Lynx pardinus Iberian lynx (LPA) Kleinman-Ruiz et al. 2022 HiSeq2000,v3 5X c_lp_sm_022 6 Lynx pardinus Iberian lynx (LPA) Kleinman-Ruiz et al. 2022 HiSeq2000,v3 5X c_lp_sm_027 6 Lynx pardinus Iberian lynx (LPA) Kleinman-Ruiz et al. 2022 HiSeq2000,v3 5X c_lp_sm_032 0 Lynx pardinus Iberian lynx (LPA) Kleinman-Ruiz et al. 2022 HiSeq2000,v3 5X c_lp_sm_032 5 Lynx pardinus Iberian lynx (LPA) Kleinman-Ruiz et al. 2022 HiSeq2000,v3 5X c_lp_sm_045 0 Lynx pardinus Iberian lynx (LPA) Kleinman-Ruiz et al. 2022 HiSeq2000,v3 5X c_lp_sm_021 1 Lynx pardinus Iberian lynx (LPA) Chapter 3 NovaSeqXPlu s 33X c_lp_sm_022 0 Lynx pardinus Iberian lynx (LPA) Chapter 3 NovaSeqXPlu s 43X c_lp_sm_022 5 Lynx pardinus Iberian lynx (LPA) Chapter 3 NovaSeqXPlu s 33X c_lp_sm_023 9 Lynx pardinus Iberian lynx (LPA) Chapter 3 NovaSeqXPlu s 40X c_lp_sm_027 8 Lynx pardinus Iberian lynx (LPA) Chapter 3 NovaSeqXPlu s 38X c_lp_sm_037 7 Lynx pardinus Iberian lynx (LPA) Chapter 3 NovaSeqXPlu s 33X c_lp_sm_039 0 Lynx pardinus Iberian lynx (LPA) Chapter 3 NovaSeqXPlu s 32X c_lp_sm_045 2 Lynx pardinus Iberian lynx (LPA) Chapter 3 NovaSeqXPlu s 34X c_lp_sm_047 4 Lynx pardinus Iberian lynx (LPA) Chapter 3 NovaSeqXPlu s 40X 152 c_lp_sm_061 4 Lynx pardinus Iberian lynx (LPA) Chapter 3 NovaSeqXPlu s 53X c_ll_ca_0240 Lynx lynx Southern Eurasian lynx (SEL) Chapter 1 HiSeq2000,v4 39X c_ll_ca_0241 Lynx lynx Southern Eurasian lynx (SEL) Chapter 1 HiSeq2000,v4 24X c_ll_ca_0242 Lynx lynx Southern Eurasian lynx (SEL) Chapter 1 HiSeq2000,v4 24X c_ll_ca_0243 Lynx lynx Southern Eurasian lynx (SEL) Chapter 1 HiSeq2000,v4 29X c_ll_ca_0244 Lynx lynx Southern Eurasian lynx (SEL) Chapter 1 HiSeq2000,v4 20X c_ll_ca_0245 Lynx lynx Southern Eurasian lynx (SEL) Chapter 1 HiSeq2000,v4 10X c_ll_ca_0247 Lynx lynx Southern Eurasian lynx (SEL) Chapter 1 HiSeq2000,v4 14X c_ll_ca_0248 Lynx lynx Southern Eurasian lynx (SEL) Chapter 1 HiSeq2000,v4 12X c_ll_ca_0252 Lynx lynx Southern Eurasian lynx (SEL) Chapter 1 HiSeq2000,v4 18X c_ll_ca_0254 Lynx lynx Southern Eurasian lynx (SEL) Chapter 1 HiSeq2000,v4 15X c_ll_ca_0259 Lynx lynx Southern Eurasian lynx (SEL) Chapter 1 HiSeq2000,v4 13X c_ll_ca_0260 Lynx lynx Southern Eurasian lynx (SEL) Chapter 1 HiSeq2000,v4 17X c_ll_ki_0090 Lynx lynx Western Eurasian lynx (SEL) Lucena-Perez et al. 2020 HiSeq X-10 23X c_ll_ki_0091 Lynx lynx Western Eurasian lynx (SEL) Lucena-Perez et al. 2020 HiSeq2000,v3 6X c_ll_ki_0092 Lynx lynx Western Eurasian lynx (SEL) Lucena-Perez et al. 2020 HiSeq2000,v3 5X c_ll_ki_0093 Lynx lynx Western Eurasian lynx (SEL) Lucena-Perez et al. 2020 HiSeq2000,v3 5X c_ll_ki_0094 Lynx lynx Western Eurasian lynx (SEL) Lucena-Perez et al. 2020 HiSeq2000,v3 6X c_ll_ki_0095 Lynx lynx Western Eurasian lynx (SEL) Lucena-Perez et al. 2020 HiSeq2000,v3 5X c_ll_ki_0096 Lynx lynx Western Eurasian lynx (SEL) Lucena-Perez et al. 2020 HiSeq2000,v3 6X c_ll_ki_0097 Lynx lynx Western Eurasian lynx (SEL) Lucena-Perez et al. 2020 HiSeq2000,v3 6X c_ll_ki_0098 Lynx lynx Western Eurasian lynx (SEL) Lucena-Perez et al. 2020 HiSeq2000,v3 6X c_ll_ki_0099 Lynx lynx Western Eurasian lynx (SEL) Lucena-Perez et al. 2020 HiSeq2000,v3 6X c_ll_ki_0100 Lynx lynx Western Eurasian lynx (SEL) Lucena-Perez et al. 2020 HiSeq2000,v3 6X c_ll_ki_0101 Lynx lynx Western Eurasian lynx (SEL) Lucena-Perez et al. 2020 HiSeq2000,v3 6X 153