De la genética a la genómica en la conservación del lince ibérico

DE LA GENÉTICA A LA GENÓMICA EN LA CONSERVACIÓN DEL LINCE IBÉRICO DANIEL KLEINMAN RUIZ TESIS DOCTORAL 2020 MEMORIA PRESENTADA POR DANIEL KLEINMAN RUIZ Para optar al título de Doctor por la Universidad de Sevilla Manual de identidad corporativa “Estación Biológica de Doñana” Estación Biológica de Doñana, Consejo Superior de Investigaciones Científicas (EBD-CSIC). Departamento de Ecología Integrativa. Grupo de Biología Evolutiva. Programa doctoral de Biología Integrada de la Universidad de Sevilla. Línea de Biología Animal: Fisiología, Biotecnología, Biodiversidad, Evolución y Conservación. MEMORIA PRESENTADA POR DANIEL KLEINMAN RUIZ Para optar al título de Doctor por la Universidad de Sevilla Manual de identidad corporativa “Estación Biológica de Doñana” Estación Biológica de Doñana, Consejo Superior de Investigaciones Científicas (EBD-CSIC). Departamento de Ecología Integrativa. Grupo de Biología Evolutiva. Programa doctoral de Biología Integrada de la Universidad de Sevilla. Línea de Biología Animal: Fisiología, Biotecnología, Biodiversidad, Evolución y Conservación. CON EL VISTO BUENO DE Dr. JOSÉ ANTONIO GODOY LÓPEZ Director. Investigador en la Estación Biológica de Doñana del Consejo Superior de Investigaciones Científicas. Dr. JESÚS FERNÁNDEZ MARTÍN Codirector. Investigador en el Instituto Nacional de Investigación y Tecnología Agraria y Alimentaria. Dra. AURORA GARCÍA-DORADO GARCÍA Codirectora. Catedrática del Departamento de Genética, Fisiología y Microbiología de la Universidad Complutense de Madrid. Dr. JUAN FRANCISCO BELTRÁN GALA Tutor. Profesor Titular del Departamento de Zoología de la Universidad de Sevilla. CON EL VISTO BUENO DE Dr. JOSÉ ANTONIO GODOY LÓPEZ Director. Investigador en la Estación Biológica de Doñana del Consejo Superior de Investigaciones Científicas. Dr. JESÚS FERNÁNDEZ MARTÍN Codirector. Investigador en el Instituto Nacional de Investigación y Tecnología Agraria y Alimentaria. Dra. AURORA GARCÍA-DORADO GARCÍA Codirectora. Catedrática del Departamento de Genética, Fisiología y Microbiología de la Universidad Complutense de Madrid. Dr. JUAN FRANCISCO BELTRÁN GALA Tutor. Profesor del Departamento de Zoología de la Universidad de Sevilla. +01'00' Los trabajos que contiene esta tesis fueron financiados a través del Ministerio español de Ciencia e Innovacion (MICINN), anteriormente Ministerio de Economía y Competitividad (proyectos CGL2006-10853/BOS, CGL2010-21540/BOS, CGL2013-47755-P y CGL2017-84641-P a José A. Godoy López, y proyecto PGC2018-095810-B-I00 a Aurora García-Dorado García), mediante una beca “Proyectos Cero” con fondos del Banco Santander, CSIC y Fundación General CSIC, y a través de un contrato entre el CSIC y la Consejería de Medio Ambiente de la Junta de Andalucía para la caracterización genética de individuos y poblaciones andaluzas de lince ibérico. Soporte adicional lo proporcionaron los proyectos LIFE06NAT/E/000209 y LIFE10NAT/ES/570, el Centro Nacional de Ciencia de Polonia (proyecto 2014/15/B/NZ8/00212), la Fundación de Ciencia Rusa (proyecto 18-14-00093), y la Comisión Europea a través del proyecto BIOGEAST (PIRSESGA-2009-247652) enmarcado en las Acciones Marie Skłodowska-Curie del Séptimo Programa Marco de Investigación y Desarrollo Tecnológico. El autor, Daniel Kleinman Ruiz, recibió financiación en el periodo 2015-2019 del programa de ayudas de Fundación La Caixa (ID 100010434) para contratos de formación predoctoral en Centros de Excelencia Severo Ochoa (SEV-2012-0262) mediante el proyecto LCF/BQ/SO15/52260006, y en el año 2020 a través de un contrato de personal investigador vinculado al precitado proyecto PGC2018-095810-B-I00 del MICINN. Es la unidad formada por el símbolo y el logotipo en una de las dos únicas combinaciones válidas. El conjunto corporativo es el que normalmente constituye la unidad básica de identificación. La utilización de cada versión normalizada vendrá siempre determinada por las pautas reflejadas en este manual. Estos logotipos en color sólo pueden ir sobre fondo blanco o negro. Sobre otros fondos deberan colocarse estos logos en versión monocromo blanco o negro. Cuando la publicación no sea en cualquier idioma del estado español, se utilizará el logo en inglés. El conjunto corporativo de Obra Social ”la Caixa” Versión horizontal Versión vertical Versión horizontal monocromo Versión vertical monocromo Versión horizontal internacional Versión vertical internacional Manual de aplicación y convivencia de identidad corporativa para entidades beneficiarias del programa ”la Caixa” Empleo Joven Todos los materiales deben pasar por Normativa de la Fundación Bancaria ”la Caixa” para validar que son correctos según las normas. Para ello, se deben enviar a la dirección de correo electrónica: [email protected] CITACIÓN RECOMENDADA Kleinman-Ruiz, Daniel (2020). De la genética a la genómica en la conservación del lince ibérico. Tesis doctoral. Universidad de Sevilla, Sevilla, España. A mamá A papá Y a mis dos extraordinarias abuelas 16 RESUMEN En conjunto, los trabajos recogidos en esta tesis permiten validar el principal vehículo de gestión genética del lince ibérico en la actualidad, brindan nuevas herramientas genómicas con el potencial de vertebrar la conservación de la especie a largo plazo, y amplían el conocimiento de las dinámicas del lastre deletéreo en una especie al límite de la extinción, lo que puede resultar fundamental para predecir el destino de sus poblaciones y mejorar la eficacia de las acciones de conservación. 17 DE LA GENÉTICA A LA GENÓMICA EN LA CONSERVACIÓN DEL LINCE IBÉRICO ABSTRACT Our planet is in the midst of its sixth extinction crisis as a consequence of human activity. In the face of this global threat, conservation biology emerges with the aim of protecting and managing biodiversity, and studying the underlying causes of its loss. Given that the genetic erosion of populations —originally a consequence of their decline— may end up becoming one of these causes, its study and management is a pressing need. In recent decades, conservation genetics has devoted itself to this task driven by the refinement of molecular markers. Today, the development of massive sequencing technologies and genomics opens up new opportunities to advance the understanding of the genomic consequences of decline, and apply this knowledge and techniques to the preservation of species. In this era of integration of genomics into conservation genetics, those species that have undergone a more intense genetic erosion, and which, at the same time, are the target of active conservation and monitoring measures, and have more genomic resources available, are poised to become model species and drive forward both the theoretical and applied aspects of the discipline. The Iberian lynx (Lynx pardinus) meets all the requirements: the dramatic population decline that the species experienced during the 20th century earned it the distinction of being the world’s most endangered feline; there is an ever-growing extensive literature on the genetic deterioration of the species; in the past twenty years it has been subject to numerous preservation programmes that have led to an outstanding demographic recovery; and a reference genome for the species has recently been developed together with other genomic resources with enormous potential for impacting conservation. On this background, this thesis aims to integrate genomics into the conservation genetics of the model Iberian lynx. Chapter 1 evaluates the genetic status of the ex situ population, a key element in the extraordinary demographic recovery of the species, in order to confirm whether this demographic success is endorsed by appropriate and effective genetic management. The analyses, based mainly on the genotypes of the 36 microsatellites currently used to manage the population genetically, reveal a high representativeness of the genetic diversity found in the two remnant populations, as well as decreasing levels of inbreeding over time, and an expected heterozygosity close to the theoretical maximum. Altogether, these results support both the selection of founders and the genetic management of the population. At the same time, the verification of the 18 RESUMEN existence of a moderate degree of kinship among founders shows the convenience of integrating this information into management. Finally, this study provides some practical recommendations for the managers of the programme. Chapter 2 seeks to identify highly efficient (i.e. informative and economically profitable) subsets of Single-Nucleotide Polymorphisms (SNPs) for conservation purposes, by performing a screening of ~ 1,500 SNPs previously selected based on the analysis of whole genome sequences of 11 Iberian lynx individuals, and genotyped in an extended sample. With a modest increase in the number of markers (24-96), these panels offer greater power for individual identification, paternity assignment, kinship estimation or ancestry classification than the 12-36 currently used microsatellites. The combination of the nested nature of most of these subsets, the flexibility and automation afforded by microfluidic platforms, and the suitability of SNPs for the genetic analysis of non-invasive samples, place these panels as the optimal tool to establish a routine monitoring of individuals and comprehensive genetic management of all populations. Chapter 3 unlocks the potential of whole genome sequence data to characterise the deleterious variation of the Iberian lynx and the Eurasian lynx (Lynx lynx), two sister species with contrasting demographic histories and distribution ranges, and which together exhibit a wide array of recent population sizes. The comparative analysis of individual derived allele counts reveals a relative genomic load deficit in the Iberian lynx in the putatively most deleterious categories, which is supported by a statistical framework that accounts for both the sample and the evolutionary variance. This pattern is consistent with the demographic history of both species and suggests a higher prevalence of genetic purging in the Iberian lynx, which could have contributed towards the persistence of the species despite its extreme genetic erosion. Furthermore, this result is supported by the observation of less differential purging in genomic regions in which selection is expected to be less efficient (centromeres), or in which selection against recessive alleles is more efficient regardless of demography (X chromosome). Within each species, the genomic load does not vary among populations, probably because the census differences are too small and recent, but the smaller sized populations show more uniform allele frequency distributions, increased homozygosity, and a greater proportion of locally fixed alleles even within deleterious categories, which points to an accumulation of inbreeding depression due to recessive alleles. Finally, this study provides a genome-wide catalogue of thousands of potentially deleterious mutations that may become a valuable resource for future conservation actions. 19 DE LA GENÉTICA A LA GENÓMICA EN LA CONSERVACIÓN DEL LINCE IBÉRICO Altogether, the research work developed in this thesis validates the main actual resource for the genetic management of the Iberian lynx, provides new genomic tools with the potential to become the backbone of the long-term conservation of the species, and expands the knowledge of the deleterious load dynamics in a species at the verge of extinction, which can be essential to predict the fate of its populations and to improve the effectiveness of current conservation efforts. 20 21 INTRODUCCIÓN 22 INTRODUCCIÓN 23 DE LA GENÉTICA A LA GENÓMICA EN LA CONSERVACIÓN DEL LINCE IBÉRICO EL VALOR DE LA BIODIVERSIDAD Y SU CONSERVACIÓN Año 2020. Nos encontramos inmersos en la sexta crisis de extinción masiva, cuyas tasas en la actualidad superan en cientos o incluso miles de veces las basales (Barnosky et al. 2011; Pimm et al. 2014; Ceballos et al. 2015). Se trata de un problema global de primer orden, puesto que todas las especies juegan un papel en el entramado de sistemas que sustenta la vida, y su pérdida es un proceso irreversible. De forma directa o indirecta, la causa subyacente es la actividad humana, cuyo impacto a gran escala se inició con la revolución agrícola, hace unos 11.000 años (Zeder 2011), y no ha dejado de acelerarse desde entonces (Naciones Unidas, Departmento de Asuntos Económicos y Sociales, División de Población 2019). Por lo general, las especies no son capaces de adaptarse a cambios tan drásticos, por lo que no es de extrañar que muchas se hayan extinguido recientemente, y hasta un 25-50% de las remanentes afrontan un riesgo elevado de extinción en las próximas décadas (Koh et al. 2004; Diaz et al. 2019). Uno de los factores imprescindibles para hacer frente a esta crisis consiste en cultivar una amplia comprensión y concienciación pública del valor de la biodiversidad, y promover la participación ciudadana en acciones que puedan prevenir su declive. Para ello, el Convenio sobre la Diversidad Biológica (CDB) fijó en el año 2010 como la primera de sus Metas de Aichi para 2020 que “las personas tengan conciencia del valor de la diversidad biológica y de los pasos que pueden seguir para su conservación y utilización sostenible” (Convenio sobre la Diversidad Biológica 2010). El valor que tienen las especies, y por tanto aquel asociado a su pérdida, se ha tratado de ilustrar habitualmente mediante sus múltiples funciones o beneficios. Tal es el caso de la conocida clasificación de Pearce & Moran (1994), que distingue cuatro categorías de valores económicos o tangibles aplicables a la biodiversidad: el de uso directo (recursos que producen beneficio económico o social cuantificable, como por ejemplo aquellos relacionados con la alimentación o el turismo), el de uso indirecto (servicios ecosistémicos como la polinización o los ciclos de nutrientes), el potencial (que recoge los posibles usos futuros directos e indirectos), y el de legado (que garantiza la disponibilidad de la diversidad biológica y el funcionamiento de los ecosistemas para las generaciones futuras). Todos estos valores son sin duda importantes y útiles para concienciar a la opinión general, y por ello se están fomentando en la actualidad desde numerosos organismos gubernamentales y no gubernamentales. Pero existe otro valor relegado al último lugar en la clasificación de Pearce y Moran que es difícilmente mensurable en términos económicos, pero no así en términos éticos: el valor de existencia o intrínseco de la biodiversidad, que es el resultado de miles o 24 INTRODUCCIÓN millones de años de contingencias y confluencias irrepetibles que caracterizan la evolución. La extinción de especies es un proceso permanente que conlleva la pérdida de un inmenso abanico de caracteres únicos, incluyendo genes, fenotipos, funciones, comportamientos, roles ecológicos, y las complejas interacciones entre todos ellos. Entender el entramado de la vida únicamente desde un punto de vista antropocentrista, en términos de usos y beneficios directos o indirectos, potenciales o realizados es, como discuten Kopnina et al. (2018), un fracaso ético que debe ser resuelto. Afortunadamente, la conservación de la biodiversidad, que además de difundir sus valores engloba también su manejo y protección y el estudio de las causas subyacentes a su pérdida, ha recibido un enorme impulso académico desde que hacia 1980 comenzase a desarrollarse el cuerpo teórico de la biología de la conservación (Soulé & Wilcox 1980). Se trata de una disciplina científica de crisis en la que confluyen conocimientos de las ciencias naturales y sociales, así como aspectos prácticos del manejo de recursos naturales, y que busca fundamentar las políticas de conservación en la teoría biológica. Sin conservación, más especies aún estarían amenazadas en el presente o se habrían extinguido (Butchart et al. 2006; Hoffmann et al. 2015). Bajo su paraguas se han identificado y confirmado numerosos factores de amenaza de la biodiversidad relacionados con la actividad antrópica. Entre ellos destacan la pérdida y degradación de hábitats, la sobreexplotación, la contaminación, el cambio climático, el comercio ilegal de especies, y la introducción de especies foráneas (ya sea domésticas o silvestres) que terminan por convertirse en invasoras (Dirzo et al. 2014; Ceballos et al. 2015; Ceballos et al. 2017; Lister & Garcia 2018; Diaz et al. 2019; Unión Internacional para la Conservación de la Naturaleza 2020). Por añadidura, muchos de estos factores pueden interactuar y retroalimentarse entre sí. TEORÍA DE LA GENÉTICA DE LA CONSERVACIÓN El sendero hacia la extinción de una especie cursa a menudo con el aislamiento y la reducción de sus poblaciones. Las poblaciones pequeñas son extremadamente susceptibles a eventos estocásticos de índole demográfica o ambiental (Lande 1993), pero también a numerosos procesos de deterioro genético que pueden precipitar su extinción. En muchos de ellos juega un papel clave la deriva genética, que es el cambio intergeneracional en las frecuencias alélicas de una población debido al efecto de muestreo que se introduce en cada evento de reproducción (Wright 1939; Masel 2011) y, como fuerza estocástica, es tanto más intensa cuanto menor sea el tamaño poblacional. Por un lado, la deriva conduce a una pérdida generalizada de diversidad genética y, por ende, de potencial evolutivo de las poblaciones, que es el salvoconducto necesario para adaptarse a los cambios medioambientales (Frankham 2005); 25 DE LA GENÉTICA A LA GENÓMICA EN LA CONSERVACIÓN DEL LINCE IBÉRICO esta situación cobra particular relevancia en el escenario actual de cambio global acelerado. Por otro lado, el número esperado de mutaciones deletéreas por individuo que se han acumulado durante un periodo evolutivo cualquiera está gobernado por el equilibrio entre la magnitud de sus efectos deletéreos y el tamaño efectivo poblacional, o dicho de otro modo, entre las fuerzas de la selección natural y la deriva genética (Kimura 1962; Ohta 1973; Lande 1995), como se explica un poco más adelante. Otro factor de deterioro genético que cobra relevancia en poblaciones pequeñas de especies con reproducción sexual es la endogamia, es decir, la producción de descendencia a partir de individuos relacionados genéticamente. En esta descendencia, consanguínea, la homocigosis aumentada conlleva una mayor manifestación de los efectos deletéreos recesivos. La consiguiente pérdida de eficacia biológica de la población se conoce como depresión consanguínea (Darwin 1876; Ralls et al. 1979; Charlesworth & Charlesworth 1987). En conjunto, todos estos procesos que afectan negativamente a los componentes de la eficacia biológica conforman un funesto bucle de retroalimentación positiva que conduce a tamaños poblacionales progresivamente menores, y que se conoce como vórtice de extinción (Gilpin & Soulé 1986). Una vez que una población entra en esta espiral, los riesgos de naturaleza genética pueden persistir incluso si desaparecen las causas primeras que propiciaron el declive (Frankham 2003). Esto significa, en definitiva, que los procesos de deterioro genético de poblaciones pequeñas, originalmente consecuencia del declive, pueden acabar convirtiéndose en causa de intensificación del mismo. Numerosas evidencias empíricas avalan este paradigma de la genética de la conservación (Fagan & Holmes 2006; Blomqvist et al. 2010; Palomares et al. 2012; Norén et al. 2016). Así pues, incluso en aquellos casos en que se pudiesen implementar medidas de conservación a gran escala como la preservación o restauración de hábitats, o la lucha contra el cambio climático, esenciales para revertir el deterioro global y prevenir que miles de especies entren en vórtice de extinción, estas medidas podrían resultar insuficientes por sí mismas para evitar la desaparición de muchas especies que ya están en una situación especialmente vulnerable. Por tanto, el estudio, comprensión y abordaje de los procesos de deterioro genético resulta una necesidad creciente y apremiante en el actual contexto de crisis de la biodiversidad. Si bien los mecanismos del deterioro genético de la eficacia biológica están bien establecidos, aún quedan cuestiones por esclarecer en relación con las dinámicas del lastre genético, que incluye: i) el lastre atribuible a deletéreos segregantes, tanto expresado (esto es, la reducción en la eficacia promedio de una población por causas genéticas respecto de un genotipo sin ninguno de estos deletéreos) como oculto en heterocigosis (lastre de consanguinidad) (Morton et al. 1956; Kimura et al. 1963), y 32 INTRODUCCIÓN Numerosos estudios publicados en los últimos años documentan exhaustivamente la erosión genética del lince ibérico. En primer lugar, su diversidad genética actual, que ya se sospechaba muy reducida en base a datos de ADN mitocondrial y microsatélites (Johnson et al. 2004; Casas-Marce et al. 2013), es extremadamente baja a nivel de genoma completo (π = 0.026%), de las menores reportadas hasta la fecha en especies animales, y aproximadamente la mitad en DON que en AND (Abascal et al. 2016). Además se ha inferido que los niveles de diversidad mitogenómica de la especie milenios atrás, si bien más altos que en la actualidad, ya eran bajos de por sí (Rodríguez-Varela et al. 2015; Casas-Marce et al. 2017). El hecho de que ambas poblaciones exhiban una diferenciación genética elevada (Casas-Marce et al. 2013), moldeada por la deriva genética reciente (Casas-Marce et al. 2017), sugiere que sus acervos representan hasta cierto punto subconjuntos relativamente complementarios de la diversidad histórica. Por otro lado, se han observado señales de deterioro genético de su eficacia biológica, que se concretan en: i) un lastre de consanguinidad elevado, como sugiere la correlación positiva entre la proporción de esperma anormal y el nivel de consanguinidad individual (Ruiz-López et al. 2012); ii) una elevada incidencia de depleción linfoide (Peña et al. 2006), glomerulonefritis membranosa (Jiménez et al. 2008), y otros rasgos deletéreos con probable base genética como una epilepsia idiopática (Martínez et al. 2013); y iii) una abundancia relativa de variantes no sinónimas frente a las sinónimas que podría sugerir una relajación de la selección purificadora (Abascal et al. 2016), si bien dicha señal podría estar modulada por otras fuerzas (Simons & Sella 2016). Como factor de riesgo adicional, en la población de DON se han descrito apareamientos entre parientes próximos (Lucena-Perez et al. 2018), lo que por otro lado es esperable en poblaciones pequeñas incluso en condiciones de panmixia. En paralelo, se ha observado una disminución a lo largo del tiempo de la supervivencia y el tamaño de la camada concomitante con un aumento de la homocigosidad, lo que se puede considerar un indicio de depresión consanguínea (Palomares et al. 2012). Por ello, se ha sugerido que la población de DON ha entrado en vórtice de extinción (Palomares et al. 2012). Este escenario tan dramático propició a principios del siglo XXI la implementación de un amplio abanico de estrategias activas de conservación, financiadas fundamentalmente por cuatro proyectos sucesivos del programa europeo LIFE en el periodo 1994-2018 (Simón et al. 2012). El repertorio incluye numerosas medidas no relacionadas con la genética, tales como: i) la restauración del hábitat; ii) la recuperación del conejo de monte (Oryctolagus cuniculus), que constituye la presa principal del lince ibérico; iii) la concienciación ciudadana y el establecimiento de acuerdos con propietarios de fincas particulares, cazadores y gestores; o iv) la asistencia sanitaria ante enfermedades infecciosas, entre otras muchas. Asimismo, 33 DE LA GENÉTICA A LA GENÓMICA EN LA CONSERVACIÓN DEL LINCE IBÉRICO incluye una monitorización intensiva de las poblaciones de lince (Simón et al. 2012), consistente principalmente en una combinación de radio-seguimiento y fototrampeo (Garrote et al. 2011; Gil-Sánchez et al. 2011; Garrote et al. 2014) (Garrote et al. 2011; Gil-Sánchez et al. 2011; Garrote et al. 2014). Esta monitorización, si bien resulta muy útil para estimar censos y analizar las dinámicas demográficas de las poblaciones, hasta la fecha no incorpora las valiosas capas de información genética a las que permiten acceder los marcadores moleculares. De cualquier manera, la conservación del lince ibérico también incluye algunas importantes medidas cimentadas en la genética, a saber: el rescate genético de DON mediante la translocación de individuos de AND (Simón et al. 2012); un programa de reproducción ex situ que gestiona genéticamente (en base a 36 microsatélites y la genealogía) una población cautiva híbrida (Vargas et al. 2008; Godoy et al. 2009; Vargas, Rivas, et al. 2009; Vargas, Breitenmoser, et al. 2009); y un programa de reintroducción de individuos procedentes de cautividad destinado a recuperar la distribución histórica de la especie en la península ibérica (Vargas et al. 2008), que hasta la fecha ha fundado seis poblaciones. Como denominador común de estas medidas se sitúa la mezcla de ambos acervos remanentes, AND y DON, a fin de minimizar la endogamia y recuperar parcialmente la diversidad perdida localmente durante los más recientes cuellos de botella, como avalan las inferencias extraídas de una reconstrucción de la historia evolutiva reciente de la especie (Casas-Marce et al. 2017). El éxito indiscutible de todas estas acciones de conservación ha sido reflejado en la reclasificación de la especie de la categoría “En peligro crítico (CR)” a la categoría “En peligro (EN)” en la Lista Roja de 2015 (Rodríguez & Calzada 2015) que elabora la Unión Internacional para la Conservación de la Naturaleza (UICN). También por el número de individuos en libertad, que en los últimos años está creciendo a un ritmo notable: 686 en el censo de 2018 (http://www.iberlince.eu/images/docs/3_ InformesLIFE/Informe_Censo_2018.pdf), 894 en 2019 (según las notas de prensa distribuidas por el Ministerio para la Transición Ecológica y el Reto Demográfico), y se estima que por encima del hito de 1.000 individuos en 2020, es decir, unas diez veces más que al inicio de las medidas. Con todo, los datos demográficos y genéticos del lince ibérico distan mucho de los de especies fuera de peligro, y su supervivencia a largo plazo sigue dependiendo de acciones de conservación continuadas en el tiempo. Estas medidas deben adaptarse a un escenario de creciente complejidad que involucra un número cada vez mayor de áreas de reintroducción, individuos migrantes, y grados de admixture o mezcla genética. Semejante entramado de conservación plantea nuevos desafíos técnicos, y exige herramientas moleculares altamente informativas y eficientes para alcanzar una monitorización y una gestión integral de la especie. 34 INTRODUCCIÓN Más allá de esta batería de medidas activas de conservación y de las oportunidades de investigación que comportan, más allá de la exhaustiva literatura sobre los patrones de deterioro genético de la especie, edificada en la estela de un profundo conocimiento acerca de su ecología (Delibes 1980; Fernández et al. 2006; Revilla & Wiegand 2008) y su declive demográfico (Rodríguez & Delibes 1992; Rodríguez & Delibes 2002; Calzada et al. 2007), el lince ibérico se erige como modelo para la genética de la conservación por la convergencia de factores adicionales. Cuenta con un genoma de referencia (Abascal et al. 2016), el primer genoma complejo cuya secuenciación se ha liderado y financiado principalmente desde España junto con los del olivo (Olea europaea) y el rodaballo (Scophthalmus maximus) (Cruz et al. 2016; Figueras et al. 2016). Además, como recursos genómicos adicionales la especie dispone de una anotación profusa, la sintenia con el genoma de referencia del gato doméstico, un chip Illumina GoldenGate de 1.536 SNPs seleccionados por sus frecuencias alélicas intermedias, y varios individuos resecuenciados procedentes de ambas poblaciones (Abascal et al. 2016). A estos se suma alrededor de un centenar de individuos resecuenciados (Lucena-Perez et al. 2020) de la especie hermana, el lince euroasiático o boreal (Lynx lynx). Se trata de uno de los felinos más ampliamente distribuidos del mundo, cuyos censos históricos son muy superiores a los del lince ibérico (Abascal et al. 2016), y cuyas poblaciones cubren toda una gama de zonas climáticas, tamaños y estados de conservación, proporcionando un contrapunto así como un rango extendido de escenarios demográficos respecto al lince ibérico. Así pues, el estudio comparativo de ambas especies puede arrojar luz sobre numerosos aspectos evolutivos relevantes para la conservación de sus poblaciones. Por último, cabe destacar como factor adicional el carácter emblemático del lince ibérico, que promueve el apoyo ciudadano, político e incluso científico a las inversiones necesarias para su conservación. En conjunto, estos factores le otorgan el rol de especie paraguas, y lo sitúan en una posición privilegiada para avanzar el campo de la conservación desde los nuevos enfoques, tanto teóricos como aplicados, a los que proporciona acceso la genómica. 35 DE LA GENÉTICA A LA GENÓMICA EN LA CONSERVACIÓN DEL LINCE IBÉRICO REFERENCIAS Abascal F, Corvelo A, Cruz F, Villanueva-Cañas JL, Vlasova A, Marcet-Houben M, Martínez-Cruz B, Cheng JY, Prieto P, Quesada V, et al. 2016. Extreme genomic erosion after recurrent demographic bottlenecks in the highly endangered Iberian lynx. Genome Biol. 17:251. Agrawal AF, Whitlock MC. 2011. Inferences about the distribution of dominance drawn from yeast gene knockout data. Genetics 187(2):553–566. Allendorf FW, Hohenlohe PA, Luikart G. 2010. Genomics and the future of conservation genetics. Nat. Rev. Genet. 11(10):697–709. Avila V, Amador C, García-Dorado A. 2010. The purge of genetic load through restricted panmixia in a Drosophila experiment. J. Evol. Biol. 23(9):1937–1946. De Barba M, Waits LP, Garton EO, Genovesi P, Randi E, Mustoni A, Groff C. 2010. The power of genetic monitoring for studying demography, ecology and genetics of a reintroduced brown bear population. Mol. Ecol. 19(18):3938–3951. Barnosky AD, Matzke N, Tomiya S, Wogan GOU, Swartz B, Quental TB, Marshall C, McGuire JL, Lindsey EL, Maguire KC, et al. 2011. Has the Earth’s sixth mass extinction already arrived? Nature 471(7336):51–57. Beckstein C, Böcker S, Bogdan M, Bruelheide H, Bücker HM, Denzler J, Dittrich P, Grosse I, Hinneburg A, König-Ries B, et al. 2014. Explorative analysis of heterogeneous, unstructured, and uncertain data - A computer science perspective on biodiversity research. En: Proceedings of 3rd International Conference on Data Management Technologies and Applications - Volume: 1 DATA. p. 251–257. Beja-Pereira A, Oliveira R, Alves PC, Schwartz MK, Luikart G. 2009. Advancing ecological understandings through technological transformations in noninvasive genetics. Mol. Ecol. Resour. 9(5):1279–1301. Bersabé D, García-Dorado A. 2013. On the genetic parameter determining the efficiency of purging: an estimate for Drosophila egg-to-pupae viability. J. Evol. Biol. 26(2):375–385. Blomqvist D, Pauliny A, Larsson M, Flodin L-Å. 2010. Trapped in the extinction vortex? Strong genetic effects in a declining vertebrate population. BMC Evol. Biol. 10:33. 36 INTRODUCCIÓN Butchart SHM, Stattersfield AJ, Collar NJ. 2006. How many bird extinctions have we prevented? Oryx 40(3):266–278. Calzada J, Guzmán JN, Rodríguez A. 2007. Lynx pardinus (Temminck, 1827). Ficha libro rojo. En: Palomo L, Gisbert J, Blanco J, editores. Atlas y libro rojo de los mamíferos terrestres de España. Madrid: Dirección General para la Biodiversidad- SECEM-SECEMU. p. 345–347. De Cara MAR, Villanueva B, Toro M., Fernández J. 2013. Purging deleterious mutations in conservation programmes: combining optimal contributions with inbred matings. Heredity (Edinb). 110(6):530–537. Carroll EL, Bruford MW, DeWoody JA, Leroy G, Strand A, Waits L, Wang J. 2018. Genetic and genomic monitoring with minimally invasive sampling methods. Evol. Appl. 11(7):1094–1119. 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, et al. 2017. Spatiotemporal dynamics of genetic variation in the Iberian lynx along its path to extinction reconstructed with ancient DNA. Mol. Biol. Evol. 34(11):2893–2907. Casas-Marce M, Soriano L, López-Bao J, Godoy J. 2013. Genetics at the verge of extinction: insights from the Iberian lynx. Mol. Ecol. 22(22):5503–5515. Ceballos G, Ehrlich PR, Barnosky AD, García A, Pringle RM, Palmer TM. 2015. Accelerated modern human-induced species losses: entering the sixth mass extinction. Sci. Adv. 1(5):e1400253. Ceballos G, Ehrlich PR, Dirzo R. 2017. Biological annihilation via the ongoing sixth mass extinction signaled by vertebrate population losses and declines. Proc. Natl. Acad. Sci. U. S. A. 114(30):E6089–E6096. Charlesworth D, Charlesworth B. 1987. Inbreeding depression and its evolutionary consequences. Annu. Rev. Ecol. Syst. 18:237–268. Collins FS, Brooks LD, Chakravarti A. 1998. A DNA polymorphism discovery resource for research on human genetic variation. Genome Res. 8(12):1229–1231. Convenio sobre la Diversidad Biológica. 2010. Decisión adoptada por la conferencia de las partes en el convenio sobre la diversidad biológica en su décima reunión. X/31. The Strategic Plan for Biodiversity 2011-2020 and the Aichi Biodiversity Targets. En: Tenth meeting of the Conference of the Parties to the Convention on Biological Diversity. Nagoya, Japón. p. 249–266. 37 DE LA GENÉTICA A LA GENÓMICA EN LA CONSERVACIÓN DEL LINCE IBÉRICO Cruz F, Julca I, Gómez-Garrido J, Loska D, Marcet-Houben M, Cano E, Galán B, Frias L, Ribeca P, Derdak S, et al. 2016. Genome sequence of the olive tree, Olea europaea. Gigascience 5:29. Darwin C. 1876. The effects of cross and self fertilisation in the vegetable kingdom. Londres: Murray, John. Dawson TP, Jackson ST, House JI, Prentice IC, Mace GM. 2011. Beyond predictions: biodiversity conservation in a changing climate. Science 332(6025):53–58. Delibes M. 1980. Feeding ecology of the Spanish lynx in the Coto Doñana. Acta Theriol. (Warsz). 25(24):309–324. Diaz S, Settele J, Brondizio E, Ngo HT, Gueze M, Agard J, Arneth A, Balvanera P, Brauman K, Butchart S, et al. 2019. Summary for policymakers of the global assessment report on biodiversity and ecosystem services of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services. DiBattista JD. 2008. Patterns of genetic variation in anthropogenically impacted populations. Conserv. Genet. 9(1):141–156. Dirzo R, Young HS, Galetti M, Ceballos G, Isaac NJB, Collen B. 2014. Defaunation in the Anthropocene. Science 345(6195):401–406. Doyle JM, Katzner TE, Roemer GW, Cain JW, Millsap BA, McIntyre CL, Sonsthagen SA, Fernandez NB, Wheeler M, Bulut Z, et al. 2016. Genetic structure and viability selection in the golden eagle (Aquila chrysaetos), a vagile raptor with a Holarctic distribution. Conserv. Genet. 17:1307–1322. Van Dyke F, Lamb RL. 2020. Conservation Biology: Foundations, Concepts, Applications. 3ª ed. Springer International Publishing. Eyre-Walker A, Keightley PD. 2007. The distribution of fitness effects of new mutations. Nat. Rev. Genet. 8(8):610–618. Fagan WF, Holmes EE. 2006. Quantifying the extinction vortex. Ecol. Lett. 9(1):51– 60. Fernández N, Delibes M, Palomares F. 2006. Landscape evaluation in conservation: molecular sampling and habitat modeling for the Iberian lynx. Ecol. Appl. 16(3):1037– 1049. 38 INTRODUCCIÓN Figueras A, Robledo D, Corvelo A, Hermida M, Pereiro P, Rubiolo JA, Gómez- Garrido J, Carreté L, Bello X, Gut M, et al. 2016. Whole genome sequencing of turbot (Scophthalmus maximus; Pleuronectiformes): a fish adapted to demersal life. DNA Res. 23(3):181–192. Frankham R. 2003. Genetics and conservation biology. C. R. Biol. 326(Suppl 1):S22– S29. Frankham R. 2005. Genetics and extinction. Biol. Conserv. 126(2):131–140. Frankham R. 2010. Where are we in conservation genetics and where do we need to go? Conserv. Genet. 11(2):661–663. Frankham R, Ballou JD, Briscoe DA. 2010. Introduction to conservation genetics. Cambridge: Cambridge University Press. García-Dorado A. 2012. Understanding and predicting the fitness decline of shrunk populations: inbreeding, purging, mutation, and standard selection. Genetics 190(4):1461–1476. Garrote G, de Ayala RP, Pereira P, Robles F, Guzman N, García FJ, Iglesias MC, Hervás J, Fajardo I, Simón M, et al. 2011. Estimation of the Iberian lynx (Lynx pardinus) population in the Doñana area, SW Spain, using capture-recapture analysis of camera-trapping data. Eur. J. Wildl. Res. 57(2):355–362. Garrote G, de Ayala RP, Tellería JL. 2014. A comparison of scat counts and cameratrapping as means of assessing Iberian lynx abundance. Eur. J. Wildl. Res. 60(6):885– 889. Genome 10K Community of Scientists, Haussler D, O’Brien SJ, Ryder OA, Keith Barker F, Clamp M, Crawford AJ, Hanner R, Hanotte O, Johnson WE, et al. 2009. Genome 10K: a proposal to obtain whole-genome sequence for 10,000 vertebrate species. J. Hered. 100(6):659–674. Gil-Sánchez JM, Moral M, Bueno J, Rodríguez-Siles J, Lillo S, Pérez J, Martín JM, Valenzuela G, Garrote G, Torralba B, et al. 2011. The use of camera trapping for estimating Iberian lynx (Lynx pardinus) home ranges. Eur. J. Wildl. Res. 57(6):1203– 1211. Gilpin ME, Soulé ME. 1986. Minimum viable populations: processes of species extinction. En: Soulé ME, editor. Conservation biology: the science of scarcity and diversity. Sunderland (Massachusetts): Sinauer. p. 19–34. 39 DE LA GENÉTICA A LA GENÓMICA EN LA CONSERVACIÓN DEL LINCE IBÉRICO Godoy JA, Casas-Marcé M, Fernández J. 2009. Genetic issues in the implementation of the Iberian lynx ex situ conservation programme. En: Vargas A, Breitenmoser C, Breitenmoser U, editores. Iberian lynx ex situ conservation: an interdisciplinary approach. Madrid: Fundación Biodiversidad, en colaboración con: Unión Internacional para la Conservación de la Naturaleza, Comisión de Supervivencia de Especies (IUCN/SSC) Cat Specialist Group. p. 86–99. Grossen C, Guillaume F, Keller LF, Croll D. 2020. Purging of highly deleterious mutations through severe bottlenecks in Alpine ibex. Nat. Commun. 11(1):1001. Guichoux E, Lagache L, Wagner S, Chaumeil P, Léger P, Lepais O, Lepoittevin C, Malausa T, Revardel E, Salin F, et al. 2011. Current trends in microsatellite genotyping. Mol. Ecol. Resour. 11(4):591–611. Guzmán JN, García FJ, Garrote G, Pérez de Ayala R, Iglesias C. 2004. El lince ibérico (Lynx pardinus) en España y Portugal. Censo diagnóstico de sus poblaciones. Madrid: Dirección General para la Biodiversidad. Hayden EC. 2009. 10,000 genomes to come. Nature 462(7269):21. Hayden EC. 2014. Technology: the $1,000 genome. Nature 507(7492):294–295. Hedrick PW. 1994. Purging inbreeding depression and the probability of extinction: full-sib mating. Heredity (Edinb). 73(4):363–372. Hedrick PW, García-Dorado A. 2016. Understanding inbreeding depression, purging, and genetic rescue. Trends Ecol. Evol. 31(12):940–952. Hoban S, Arntzen JA, Bruford MW, Godoy JA, Rus Hoelzel A, Segelbacher G, Vilà C, Bertorelle G. 2014. Comparative evaluation of potential indicators and temporal sampling protocols for monitoring genetic erosion. Evol. Appl. 7(9):984–998. Hoffman EA, Blouin MS. 2004. Historical data refute recent range contraction as cause of low genetic diversity in isolated frog populations. Mol. Ecol. 13(2):271–276. Hoffmann M, Duckworth JW, Holmes K, Mallon DP, Rodrigues ASL, Stuart SN. 2015. The difference conservation makes to extinction risk of the world’s ungulates. Conserv. Biol. 29(5):1303–1313. Hunter ME, Hoban SM, Bruford MW, Segelbacher G, Bernatchez L. 2018. Next-generation conservation genetics and biodiversity monitoring. Evol. Appl. 11(7):1029–1034. 40 INTRODUCCIÓN Jiménez MA, Sánchez B, Alenza MDP, García P, López J V, Rodríguez A, Muñoz A, Martínez F, Vargas A, Peña L. 2008. Membranous glomerulonephritis in the Iberian lynx (Lynx pardinus). Vet. Immunol. Immunopathol. 121(1–2):34–43. Johnson WE, Godoy JA, Palomares F, Delibes M, Fernandes M, Revilla E, O’Brien SJ. 2004. Phylogenetic and phylogeographic analysis of Iberian lynx populations. J. Hered. 95(1):19–28. Jump AS, Hunt JM, Peñuelas J. 2006. Rapid climate change-related growth decline at the southern range edge of Fagus sylvatica. Glob. Chang. Biol. 12(11):2163–2174. Kimura M. 1962. On the probability of fixation of mutant genes in a population. Genetics 47(6):713–719. Kimura M, Maruyama T, Crow JF. 1963. The mutation load in small populations. Genetics 48(10):1303–1312. Koepfli KP, Paten B, O’Brien SJ, Antunes A, Belov K, Bustamante C, Castoe TA, Clawson H, Crawford AJ, Diekhans M, et al. 2015. The genome 10K project: a way forward. Annu. Rev. Anim. Biosci. 3:57–111. Koh LP, Dunn RR, Sodhi NS, Colwell RK, Proctor HC, Smith VS. 2004. Species coextinctions and the biodiversity crisis. Science 305(5690):1632–1634. Kopnina H, Washington H, Taylor B, J Piccolo J. 2018. Anthropocentrism: more than just a misunderstood problem. J. Agric. Environ. Ethics 31(1):109–127. Kraus RHS, VonHoldt B, Cocchiararo B, Harms V, Bayerl H, Kühn R, Förster DW, Fickel J, Roos C, Nowak C. 2015. A single-nucleotide polymorphism-based approach for rapid and cost-effective genetic wolf monitoring in Europe based on noninvasively collected samples. Mol. Ecol. Resour. 15(2):295–305. Lande R. 1993. Risks of population extinction from demographic and environmental stochasticity and random catastrophes. Am. Nat. 142(6):911–927. Lande R. 1995. Mutation and conservation. Conserv. Biol. 9:782–791. Latter BDH, Mulley JC, Reid D, Pascoe L. 1995. Reduced genetic load revealed by slow inbreeding in Drosophila melanogaster. Genetics 139(1):287–297. Lewin HA, Robinson GE, Kress WJ, Baker WJ, Coddington J, Crandall KA, Durbin R, Edwards S V., Forest F, Gilbert MTP, et al. 2018. Earth BioGenome Project: sequencing life for the future of life. Proc. Natl. Acad. Sci. U. S. A. 115(17):4325–4333. 41 DE LA GENÉTICA A LA GENÓMICA EN LA CONSERVACIÓN DEL LINCE IBÉRICO Lister BC, Garcia A. 2018. Climate-driven declines in arthropod abundance restructure a rainforest food web. Proc. Natl. Acad. Sci. U. S. A. 115(44):E10397– E10406. López-Cortegano E, Vilas A, Caballero A, García-Dorado A. 2016. Estimation of genetic purging under competitive conditions. Evolution 70(8):1856–1870. Lucena-Perez M, Marmesat E, Kleinman-Ruiz D, Martínez-Cruz B, Węcek K, Saveljev AP, Seryodkin IV, Okhlopkov I, Dvornikov MG, Ozolins J, et al. 2020. Genomic patterns in the widespread Eurasian lynx shaped by Late Quaternary climatic fluctuations and anthropogenic impacts. Mol. Ecol. 29(4):812–828. Lucena-Perez M, Soriano L, López-Bao JV, Marmesat E, Fernández L, Palomares F, Godoy JA. 2018. Reproductive biology and genealogy in the endangered Iberian lynx: implications for conservation. Mamm. Biol. 89:7–13. Marschall T, Marz M, Abeel T, Dijkstra L, Dutilh BE, Ghaffaari A, Kersey P, Kloosterman WP, Mäkinen V, Novak AM, et al. 2018. Computational pan-genomics: status, promises and challenges. Brief. Bioinform. 19(1):118–135. 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). J. Zoo Wildl. Med. 44(4):845–852. Masel J. 2011. Genetic drift. Curr. Biol. 21(20):R837–R838. McGowan PJK, Traylor-Holzer K, Leus K. 2017. IUCN guidelines for determining when and how ex situ management should be used in species conservation. Conserv. Lett. 10(3):361–366. McMahon BJ, Teeling EC, Höglund J. 2014. How and why should we implement genomics into conservation? Evol. Appl. 7(9):999–1007. Morin PA, Luikart G, Wayne RK, The SNP Workshop Group. 2004. SNPs in ecology, evolution and conservation. Trends Ecol. Evol. 19(4):208–216. Morton NE, Crow JF, Muller HJ. 1956. An estimate of the mutational damage in man from data on consanguineous marriages. Proc. Natl. Acad. Sci. U. S. A. 42(11):855–863. Murchison EP, Tovar C, Hsu A, Bender HS, Kheradpour P, Rebbeck CA, Obendorf D, Conlan C, Bahlo M, Blizzard CA, et al. 2010. The Tasmanian devil transcriptome reveals Schwann cell origins of a clonally transmissible cancer. Science 327(5961):84–87. 48 OBJETIVOS 49 DE LA GENÉTICA A LA GENÓMICA EN LA CONSERVACIÓN DEL LINCE IBÉRICO OBJETIVOS La presente tesis ilustra sobre el caso modélico del lince ibérico la ampliación de la genética a la genómica de la conservación. El capítulo 1 presenta la evaluación genética de la población ex situ de lince ibérico desde su establecimiento en el año 2003 hasta el año 2015, como pieza clave en la extraordinaria recuperación de la especie. Este programa ha tenido un éxito rotundo desde el punto de vista demográfico, puesto que su producción de individuos para reintroducción ha superado con creces incluso las previsiones más optimistas. No obstante, es importante constatar en qué grado se han cumplido los objetivos genéticos que se establecieron al principio del programa y, por tanto, cuál es el estado genético de los individuos que se están reintroduciendo en la naturaleza. El objetivo de este capítulo es evaluar el estado genético de la población cautiva de lince ibérico y, concretamente, i) si sus fundadores son representativos de los acervos genéticos remanentes (AND y DON), ii) el grado de relaciones genéticas entre ellos, iii) cómo ha evolucionado la diversidad genética y la consanguinidad con el transcurso del tiempo, y iv) el grado de estructura entre centros de cría, todo ello en base al pedigrí de la población y los genotipos de 36 microsatélites, las herramientas que hasta el día de hoy se han utilizado para su gestión genética. El capítulo 2 presenta baterías de marcadores genómicos altamente eficientes, es decir, informativos y económicamente rentables, que permitan establecer un seguimiento de individuos y una gestión genética integral de todas las poblaciones con vistas a la monitorización genética y conservación de la especie a largo plazo. Partiendo de un total de ~1.500 SNPs previamente seleccionados sobre secuencias de genoma completo de 11 individuos de lince ibérico, y genotipados en una muestra ampliada, los objetivos concretos de este capítulo consisten en: i) extraer una serie de subconjuntos anidados de SNPs altamente informativos para la identificación de individuos, la asignación de paternidades, la estimación de parentescos, y la determinación de la ascendencia, y ii) comparar su poder con el de los 36 microsatélites en los que se basa el manejo en la actualidad. 50 OBJETIVOS En el capítulo 3 se aprovecha el potencial de los datos de secuencias de genoma completo para caracterizar la variación deletérea del lince ibérico y el lince boreal, dos especies hermanas con historias demográficas y áreas de distribución contrastadas y que en conjunto exhiben un amplio rango de tamaños poblacionales recientes. La metodología se basa fundamentalmente en el cálculo del promedio de alelos derivados deletéreos por individuo (la carga genómica), y otros estadísticos basados en recuentos individuales de alelos derivados y anotados funcionalmente, y se respalda en un marco estadístico que da cuenta tanto de la varianza muestral como de la evolutiva. El objetivo de este capítulo consiste en realizar un análisis comparativo de los patrones de la variación deletérea entre el lince ibérico y el lince boreal, a fin de desvelar las siguientes incógnitas: i) si existen diferencias en carga genómica entre ambas especies, ii) si dichas diferencias ocurren entre poblaciones dentro de cada especie, iii) si dichas diferencias están originadas por diferencias en la eficacia de la selección purificadora o de la purga genética, y iv) cómo varían estas diferencias entre regiones genómicas y categorías de anotación y de efecto mutacional sometidas a distintas presiones selectivas. 51 DE LA GENÉTICA A LA GENÓMICA EN LA CONSERVACIÓN DEL LINCE IBÉRICO CAPÍTULO 1 GENETIC EVALUATION OF THE IBERIAN LYNX EX SITU CONSERVATION PROGRAMME 54 CAPÍTULO 1 55 GENETIC EVALUATION OF THE EX SITU POPULATION ABSTRACT Ex situ programmes have become critical for improving the conservation of many threatened species, as they establish backup populations and provide individuals for reintroduction and reinforcement of wild populations. The Iberian lynx was considered the most threatened felid species in the world in the wake of a dramatic decline during the second half of the 20th century that reduced its numbers to around only 100 individuals. An ex situ conservation programme was established in 2003 with individuals from the two well-differentiated, remnant populations, with great success from a demographic point of view. Here we evaluate the genetic status of the Iberian lynx captive population based on molecular data from 36 microsatellites, including patterns of relatedness and representativeness of the two remnant genetic backgrounds among founders, the evolution of diversity and inbreeding over the years, and genetic differentiation among breeding facilities. In general terms, the ex situ population harbours most of the genetic variability found in the two wild populations and has been able to maintain reasonably low levels of inbreeding and high diversity, thus validating the applied management measures, and potentially representing a model for other species in need of conservation. 56 CAPÍTULO 1 INTRODUCTION Human activity is driving species to extinction at exceptionally high rates (Pimm et al. 2014). Ever-growing threats to biodiversity, such as overexploitation, habitat destruction, and pervasiveness of invasive species (Wilson 1992), pose a challenging ordeal for species persistence (Butchart et al. 2010). At the same time, it is becoming evident that conservation actions can have a positive impact (Hoffmann et al. 2010). In situ conservation is the main focus of biological conservation based on the multiple benefits derived from the integrated protection of habitats and ecosystems. However, in the current scenario of rapid climate change and increased habitat destruction, reliance solely on in situ approaches may not ensure the long-term –or, under critical scenarios, even the short-term– persistence of many species. In contrast, the potential of ex situ conservation, once relegated to educational, research and fund-raising roles (Pritchard et al. 2012), to contribute more directly to biological conservation is rapidly increasing, thanks to improved global coordination, infrastructures, technical knowledge and practices (Pritchard et al. 2012; Redford et al. 2012). Captive populations act as an insurance reserve in case disaster strikes wild populations. Individuals in captive populations are protected from threats such as predation and poaching, they are subject to continuous surveillance, and they can receive timely medical assistance. As collateral benefits, ex situ programmes provide valuable insights into species’ biology by supporting and enhancing research on behaviour, reproduction, and physiology, among others (McGowan et al. 2017). Indeed, ex situ conservation has played a major role in the improved conservation of many vertebrate species, including the California condor (Gymnogyps californianus), Arabian oryx (Oryx leucoryx), whooping crane (Grus americana), and blackfooted ferret (Mustela nigripes) (Hoffmann et al. 2010). Still, ex situ conservation should be regarded as an additional means of preserving species, so whenever possible it should be the source for individual reintroductions into the wild and be part of an integrative approach to species-wide management in conjunction with in situ actions (Russello and Amato 2007; Pritchard et al. 2012; Redford et al. 2012; McGowan et al. 2017). More importantly, captive populations can be used to boost genetic diversity and minimize inbreeding through managed and monitored admixture of differentiated remnant populations when available. This is particularly relevant since species under extinction risk are often genetically eroded (i.e. show loss of genetic diversity and accumulation of inbreeding), as the consequence of species decline and fragmentation. However, captive populations –which are generally founded by a small number of individuals, are usually small themselves, and have an environment different from 57 GENETIC EVALUATION OF THE EX SITU POPULATION wild– carry their own additional genetic risks, namely: i) loss of genetic diversity, ii) accumulation of inbreeding, iii) incidence of genetic diseases and other deleterious traits, and iv) adaptation to captivity (Frankham 2008; Frankham et al. 2010). Genetic risks in captivity, though, can be somewhat ameliorated through adequate genetic management. Towards this goal, the commonly recommended and optimal strategy (as proven both theoretically and by simulations) is based on the minimisation of average kinship (Ballou and Lacy 1995; Meuwissen 1997; Fernández and Toro 1999; Caballero and Toro 2000; Fernández et al. 2004). The implementation of the minimum kinship strategy has traditionally been based on studbook records and pedigree analyses (Ballou et al. 2010). However, molecular data can assist and improve genetic management by assessing the representativeness of founders, detecting studbook errors and, most importantly, unravelling the typically unknown relatedness among founders. By doing so, they allow for higher accuracy in routine tasks such as identifying optimal breeding pairs, exchanging individuals between facilities, or selecting suitable candidates for reintroduction (Frankham et al. 2010). Of the four extant species in the lynx genus, the Iberian lynx (Lynx pardinus) has the smallest distribution area and population size. During the second half of the 20th century the species underwent a drastic decline that saw the majority of its populations become extinct in the span of five decades (Rodríguez and Delibes 1992; Rodríguez and Delibes 2002; Calzada et al. 2007). By 2002, less than 100 individuals remained, distributed among two isolated populations in Andalusia (southern Spain): a peripheral population in the protected area of Doñana (DON), and a central one in Sierra de Andújar in the Sierra Morena range (AND). Separated by ca. 240 km, the two populations have been effectively isolated for almost two centuries (ca. 40 generations), until the start of translocations in 2007 (Guzmán et al. 2004; Casas-Marce et al. 2017). This scenario led to the recognition of the Iberian lynx as the most endangered felid in the world (Nowell and Jackson 1996) and to its classification as “critically endangered” in the 2002 and 2008 IUCN red lists (Rodríguez and Calzada 2015). Genetic erosion in the Iberian lynx has been well documented. Studies based on mitochondrial, microsatellite and genome-wide data have revealed a high inbreeding rate, high differentiation between the two populations, and extremely low whole-ge- nome level diversity (π = 0.026%), which is among the lowest ever reported – in the same ballpark as that of the cheetah (Acinonyx jubatus) and notably lower than that of other felids and humans (Johnson et al. 2004; Palomares et al. 2012; Casas-Marce et al. 2013; Dobrynin et al. 2015; Abascal et al. 2016). Moreover, several observations indicative of inbreeding depression have been reported, including an increase 64 CAPÍTULO 1 PEDIGREE ANALYSES Pedigree analyses for a version of the studbook dated the 21st of April, 2017 (Supplementary Table S1) were carried out with the software Population Management x (PMx; Lacy et al. 2012). We estimated values for gene diversity (GD; measured as the reduction in HE with respect to the founding population) and mean inbreeding (F) over time taking into account or not the kinships among founders (empirical kinships ON vs. empirical kinships OFF, respectively), which were estimated from molecular marker data as detailed in the previous section. OPTIMAL AND REALISED ANCESTRY PROPORTIONS AMONG FOUNDERS AND THEORETICAL MAXIMUM HE We estimated the optimal proportion of founders (i.e. the contribution of each genetic stock to the pool that would maximise its genetic diversity) following the equation on page 1,369 of Toro & Caballero (2005), and on the basis of the genetic diversity within, and Nei’s minimum distance (Nei 1973) between, the two remnant populations empirically estimated from microsatellite genotypes. We often refer to the HE of this ideal mix as the theoretical maximum HE, and use it as a reference point to gauge the actual estimated HE. For comparison purposes, we also inferred from the census the actual proportion of CAP-WB individuals from each pool each year, as well as the proportion of contributing founders (defined as those CAP-WB individuals that had offspring before the end of 2017). For all practical purposes, CAP-WB individuals coming from a first backcross of an F1 individual with an individual of DON ancestry contributed to calculations as 0.75 DON and 0.25 AND founders. 65 GENETIC EVALUATION OF THE EX SITU POPULATION RESULTS GENOTYPING All samples were genotyped at a minimum of 25 (out of the 36 total) loci, with 96% of them genotyped at 30 or more. Thus, no sample was discarded because of a high proportion of missing data. The proportion of genotypes with Mendelian inconsistencies across all CAP sire-dam-offspring trios in the dataset was 1.35%, indicating a low overall error rate. COMPARING WILD AND CAPTIVE GENETIC VARIATION Genetic diversity and differentiation. Average number of alleles observed per locus for CAP (3.81; range: 2-11) was similar to previously reported values for wild populations (3.75; Casas-Marce et al. 2013), although the average HO was higher (0.509; range: 0.103-0.792) than in any of the remnant populations (0.313 and 0.467 in FRL-WB-DON and FRL-WB-AND, respectively), as expected for an admixed population (Table 1). As previously reported, the genetic diversity of CAPWB-AND was higher than that of CAP-WB-DON (Casas-Marce et al. 2013). We found no significant differences in any of the diversity statistics in comparisons between CAP-WB-AND and FRL-WB-AND populations (Wilcoxon test: P > 0.05), but NA and AR were significantly lower in CAP-WB-DON than in FRLWB-DON (P < 0.012 and P < 1.226 * 10-06, respectively). The sets of alleles in each pool of CAP-WB individuals and in their respective FRLWB population were largely overlapping, more so for AND (where only two and three alleles were private to the FRL-WB and CAP-WB populations, respectively) than for DON (where 12 alleles sampled in FRL-WB were not represented in CAP-WB, and one was exclusive to the latter). Differences between FRL-WB and CAP-WB sets were noticeably reduced when alleles with only one occurrence (singletons, most likely arising from genotyping errors or mutations) were ignored: only one and three alleles were private to FRL-WB-AND and FRL-WB-DON populations, respectively, and none was private to the CAP-WB population (Figure 1). Finally, genetic differentiation was extremely low between CAP-WB-DON and FRL-WB-DON populations (FST = -0.001; 99% CI [-0.022, 0.029]), and even lower for AND (FST = -0.009; 99% CI [-0.011, -0.006]). These slightly negative values suggest that variance is higher within than between groups. 66 CAPÍTULO 1 PPooppuullaattiioonnFFRRLL--WWBB--DDOONNCCAAPP--WWBB--DDOONNFFRRLL--WWBB--AANNDDCCAAPP--WWBB--AANNDDCCAAPP NN 50 9 54 36 239 HHEE 0.294 ± 0.037 0.300 ± 0.038 0.454 ± 0.030 0.460 ± 0.031 0.510 ± 0.026 HHOO 0.313 ± 0.041 0.358 ± 0.047 0.467 ± 0.034 0.466 ± 0.035 0.509 ± 0.027 NNAA 2.250 ± 0.134 1.944 ± 0.119 * 3.306 ± 0.202 3.333 ± 0.203 3.806 ± 0.308 AARR 1.890 ± 0.101 1.874 ± 0.108 ** 2.700 ± 0.131 2.739 ± 0.135 2.907 ± 0.144 FFIISS -0.054 [-0.093, 0.015] -0.135 [-0.216, -0.039] -0.02 [-0.058, 0.016] 0.002 [-0.057, 0.053] 0.003 [-0.019, 0.022] FFSSTT NA DD NA -0.001 [-0.018, 0.021] -0.009 [-0.010, -0.007] -0.003 ± 0.006 -0.009 ± 0.001 Table 1. Genetic diversity and differentiation statistics. N: number of genotyped individuals; HE: expected heterozygosity; HO: observed heterozygosity; NA: mean number of alleles per locus; AR: allelic richness; FIS: population inbreeding coefficient; FST: Weir & Cockerham’s fixation index; D: Jost’s genetic distance; CAP: captive population; CAPWB: captive wild-born individuals; FRLWB: free-living wild-born individuals; DON: Doñana; AND: Andújar. AR was rarefied to the smallest number of alleles seen in any of the five populations (16 alleles). Asterisks accompanying CAPWB diversity statistics denote significant differences when compared to the respective FRLWB population using a two-tailed paired Wilcoxon signed-rank test (* P < 0.012, ** P < 1.226 * 10-06). FST was estimated for the comparison between each CAPWB group and its respective FRLWB population. Standard Error is given for the diversity statistics, and 95% Confidence Interval for the F statistics. 67 GENETIC EVALUATION OF THE EX SITU POPULATION In summary, lack of structure between founders and their source population indicate that both genetic stocks are reasonably well represented in captivity, whereas diversity statistics suggest a slightly incomplete representation of wild allelic variation in founders from DON. Optimal and realised ancestry proportion among founders and theoretical maximum HE. According to theoretical predictions, and given the unbiased expected homozygosity of each genetic stock (DON: 0.703; AND: 0.541) and observed Nei’s minimum distance between the two (0.3), a proportion of 0.36 founders from DON and 0.64 from AND (~1:2 DON:AND) would yield the theoretical maximum HE of 0.54 (Supplementary Figure S3). The actual proportion of founders from each remnant population has changed over the years and is approaching –but has not yet reached– this optimum, showing a relative deficit of CAP-WB-DON individuals for most of the study period (Supplementary Figure S4). This may have contributed to a realized HE a little below the theoretical maximum at the end of the study period (see below). Molecular estimation of relatedness among CAP-WB individuals. Relatedness values among CAP-WB individuals are reported in Supplementary Table S3. Average R was notable (R = 0.098 after excluding self-relatedness). The actual distribution exhibited a marked bimodality in which the mode closer to zero amassed relatedness between DON and AND individuals and the higher mode encompassed relatedness within each genetic pool (Figure 2). The subsequent shoulder and a long tail suggest the existence of certain pairs of moderately related founders within both stocks, but mainly within the DON subset. Figure 1. Comparison of the number of alleles (in parentheses: number of alleles after removing singletons) between each captive wild-born set (CAP-WB) and its respective source population (FRL-WB). Left: Doñana population (DON); right: Andújar population (AND). Four admixed CAP-WB-DON individuals were excluded. 68 CAPÍTULO 1 0 50 100 150 200 0.0 0.1 0.2 0.3 0.4 relatedness pairs of individuals within CAP−WB−AND within CAP−WB−DON between CAP−WB−AND & CAP−WB−DON Figure 2. Relatedness between pairs of CAP-WB individuals. Different shades of grey depict comparison within the AND subset, within the DON subset, and between DON and AND individuals. CAP-WB: captive wild-born individuals; DON: Doñana; AND: Andújar. To characterise relationships between founders better, we obtained maximum likelihood inferences of relatedness among pairs within each genetic CAP-WB stock using the observed allelic frequencies of each population (Supplementary Tables S4 and S5). We obtained moderate levels of ambiguous kinship; in other words, many individuals could not be assigned unambiguously as unrelated or related, especially in the DON stock, where diversity was lower and relationships are known to be more complex (Lucena-Perez et al. 2018). Nevertheless, both stocks were confirmed to include a number of related pairs (resolved with a confidence of 95%), including some instances of FS and PO relationships (Supplementary Figure S5). EVOLUTION OF DIVERSITY OVER THE YEARS Molecular statistics. Average HO of the CAP-CB cohorts fluctuated over the years, mostly driven by the degree of admixture of the two differentiated genetic stocks of DON and AND (Figure 3a). HO peaked in years when most of the offspring were F1, decreased later when successive generations of admixture (backcrosses) became more frequent, and then stabilized in more recent years around the theoretical maximum value for HE. HO consistently stayed above HE, resulting in negative FIS throughout the entire history of the captive population, an expected outcome of the inbreeding avoidance strategy implemented by the genetic management programme, which has favoured the mating of least related couples (Figures 3a and 3b). The whole population (census approach) followed a pattern similar to that of cohorts: while HO exceeded HE in CAP-CB individuals, reaching higher values in years with a higher proportion of F1 individuals, CAP-WB individuals exhibited the opposite pattern, i.e. HE –which seemed to flatten above 0.5 following the stabilisation of allelic frequencies– was higher than average HO (Figures 4a and 4b). Since 69 GENETIC EVALUATION OF THE EX SITU POPULATION the CAP-WB group resulted from the pooling of two differentiated populations, such a Wahlund effect (Wahlund 1928) was expected. Consistent with this, global FIS for the whole CAP population was positive in the early phase of the ex situ programme (when CAP-WB was the major, or even only, component), and subsequently steadily decreased as the proportion of admixed individuals continued to increase, finally reaching a negative value in the last year of our study (Figure 4c). Regarding the structure between facilities in 2015, FST values were low between all pairs of breeding facilities, although still significant for the Acebuche-Granadilla, Granadilla-Olivilla, and Acebuche-Olivilla pairs. The highest value (0.018; 99% CI [0.009, 0.028]) was found between Acebuche and Granadilla, while the lowest (-0.005; 99% CI [-0.012, 0.001]) was observed for the Silves and Olivilla comparison (Supplementary Figure S6). Figure 3. Evolution of diversity and inbreeding statistics over time for each year’s CAP-CB births (cohorts approach). a: observed heterozygosity (HO) and expected heterozygosity (HE); dashed lines represent the theoretical maximum HE (OPT), and each FRL-WB population’s HE. b: fixation index (FIS); the dashed line indicates FIS = 0. DON: Doñana; AND: Andújar. ●● ●● ●● ●● ●● ●● ●● ●●●● ●● ●● −0.6 −0.4 −0.2 0.0 0.2 2001 2003 2005 2007 2009 2011 2013 2015 year FIS H E OPT H E AND H E DON ● ● ● ●● ● ● ●● ● ● 0.2 0.3 0.4 0.5 0.6 0.7 2001 2003 2005 2007 2009 2011 2013 2015 year H ● H O H E a b ●● ●● ●● ●● ●● ●● ●● ●●●● ●● ●● −0.6 −0.4 −0.2 0.0 0.2 2001 2003 2005 2007 2009 2011 2013 2015 year FIS H E OPT H E AND H E DON ● ● ● ●● ● ● ●● ● ● 0.2 0.3 0.4 0.5 0.6 0.7 2001 2003 2005 2007 2009 2011 2013 2015 year H ● H O H E a b 70 CAPÍTULO 1 HE OPT HE AND HE DON ● ● ● ●●●● ●●●● 0.2 0.3 0.4 0.5 0.6 0.7 2001 2003 2005 2007 2009 2011 2013 2015 year H E ● CAP−CB CAP−WB CAP HE OPT HE AND HE DON ●● ● ●●●● ●● ●● 0.2 0.3 0.4 0.5 0.6 0.7 2001 2003 2005 2007 2009 2011 2013 2015 year H O ● CAP−CB CAP−WB CAP ●● ●● ●●●●●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● −0.6 −0.4 −0.2 0.0 0.2 2001 2003 2005 2007 2009 2011 2013 2015 year F IS ●● ●● CAP−CB CAP b c a Figure 4. Evolution of diversity statistics over time for each year’s total population (census approach). a: observed heterozygosity (HO); dashed lines represent the theoretical maximum HE (OPT) and each FRL-WB population’s HE. b: expected heterozygosity (HE), with the same dashed lines as in a. c: fixation index (F); the dashed line indicates FIS = 0. CAP: captive population; CAP-CB: captive-born individuals; CAP-WB: captive wild-born individuals; FRL-WB: free-living wild-born individuals; DON: Doñana; AND: Andújar. HE OPT HE AND HE DON ● ● ● ●●●● ●●●● 0.2 0.3 0.4 0.5 0.6 0.7 2001 2003 2005 2007 2009 2011 2013 2015 year H E ● CAP−CB CAP−WB CAP HE OPT HE AND HE DON ●● ● ●●●● ●● ●● 0.2 0.3 0.4 0.5 0.6 0.7 2001 2003 2005 2007 2009 2011 2013 2015 year H O ● CAP−CB CAP−WB CAP ●● ●● ●●●●●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● −0.6 −0.4 −0.2 0.0 0.2 2001 2003 2005 2007 2009 2011 2013 2015 year F IS ●● ●● CAP−CB CAP b c a HE OPT HE AND HE DON ● ● ● ●●●● ●●●● 0.2 0.3 0.4 0.5 0.6 0.7 2001 2003 2005 2007 2009 2011 2013 2015 year H E ● CAP−CB CAP−WB CAP HE OPT HE AND HE DON ●● ● ●●●● ●● ●● 0.2 0.3 0.4 0.5 0.6 0.7 2001 2003 2005 2007 2009 2011 2013 2015 year H O ● CAP−CB CAP−WB CAP ●● ●● ●●●●●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● −0.6 −0.4 −0.2 0.0 0.2 2001 2003 2005 2007 2009 2011 2013 2015 year F IS ●● ●● CAP−CB CAP b c a 71 GENETIC EVALUATION OF THE EX SITU POPULATION 0.4 0.6 0.8 1.0 2001 2003 2005 2007 2009 2011 2013 2015 year GD Empirical kinships OFF ON 0.00 0.05 0.10 0.15 0.20 2001 2003 2005 2007 2009 2011 2013 2015 year F Empirical kinships OFF ON a b 0.4 0.6 0.8 1.0 2001 2003 2005 2007 2009 2011 2013 2015 year GD Empirical kinships OFF ON 0.00 0.05 0.10 0.15 0.20 2001 2003 2005 2007 2009 2011 2013 2015 year F Empirical kinships OFF ON a b Figure 5. Evolution of pedigree-based diversity and inbreeding statistics over time when enabling (ON) or not enabling (OFF) the empirical kinship matrix. a: genetic diversity (GD); b: inbreeding coefficient (F). Pedigree statistics. GD increased swiftly in the first few years of the breeding programme as average kinship diminished with new recruitments, but stabilised from 2007 onwards at a value higher than 0.97. When accounting for the empirical relatedness between founders, however, GD tends towards a value around 0.87 (Figure 5a). Under the assumption of zero relatedness between founders, the breeding programme has been able to keep inbreeding at bay (F < 10-4 up to 2015). When considering molecular kinships, however, the picture is starkly different (Figure 5b). In earlier years values up to nearly 0.2 were reached, likely reflecting recruitment of wild-born individuals with varying degrees of inbreeding. In later years, inbreeding has fluctuated around 0.1. 72 CAPÍTULO 1 DISCUSSION Since its inception in 2003, the Iberian lynx ex situ programme has become a major pillar on which conservation efforts to save the species from extinction rest. Although its course of action has been clearly validated in demographic terms, evaluating its success from a genetic point of view is not as straightforward. A thorough study was necessary to either back the adopted set of measures and decisions or shed light on potential concerns or amendable issues. In this study, we discovered that the species’ remnant genetic diversity is fairly well represented in the CAP population, more so for the AND genetic background than for the DON one. Considering that the heterozygosity of the CAP population has been approaching the theoretical maximum level, and its mean inbreeding has been constantly shrinking, the programme’s genetic management has been satisfactory. Moreover, differentiation between breeding facilities is also negligible. In short, the programme has thus far been able to achieve its goals from a genetic perspective. TO MIX OR NOT TO MIX A first important decision to be made when genetically differentiated populations persist is whether to mix them and manage them as a single unit or to keep them separated in captivity. In those cases where remnant populations have different local adaptations or a long history of isolation, mixing stocks could potentially lead to a reduction in the biological fitness of the hybrid offspring (outbreeding depression) due to extrinsic or intrinsic factors (Templeton 1986; Lynch 1991; The SSC Re-introduction Specialist Group 1995). In contrast, when differentiation is due to stochastic processes, and/or when the risk of inbreeding in the source populations far outweighs the risk of outbreeding, mixing affords the potential advantages of maximising genetic diversity and boosting fitness of hybrid individuals (Chesser 1983; Lacy 1987). For the Iberian lynx ex situ programme, it was decided that the two remnant populations of DON and AND should be mixed in order to restore the highest possible amount of historical diversity and to limit inbreeding depression. This decision, originally informed by accumulated signs of inbreeding depression (Peña et al. 2006; Jiménez et al. 2008; Palomares et al. 2012; Ruiz-López et al. 2012) and historical range reconstructions showing likely connections in the recent past, has recently been further validated by an extensive analysis of the historical variation in the species, which showed a progressive loss of genetic diversity within populations and the accumulation of genetic differentiation among them, resulting from the contraction and fragmentation that took place in the latter half of the 20th 73 GENETIC EVALUATION OF THE EX SITU POPULATION century (Casas-Marce et al. 2017). Even though the estimated date of divergence of the two remnant populations (ca. 200 ybp) was earlier than previously thought based on range reconstructions (around 50 ybp), the absence of major environmental differences and the evidence of gene flow in the recent past do conclusively support the management of the two remnant populations as a single unit in captivity as well as in the wild (Frankham et al. 2011). These points also validate the genetic reinforcement of DON through the translocation of animals from AND which was initiated in 2007 (Simón et al. 2012). The fitness consequences of admixture in the Iberian lynx have not been formally analysed yet, but as often observed (e.g., Frankham, 2016), the fitness of outbred individuals coming from admixture of the two inbred lines appears to be higher, and at the very least it is not lower than that of their pure, inbred counterparts, both in captivity and in the wild (Simón et al. 2012). REPRESENTATIVENESS One of the main goals of ex situ conservation is to safeguard genetic diversity of threatened species or populations, and so capturing as much of that variation as possible should be the first objective. Both genetic stocks of DON and AND seem to have been represented well in captivity, as revealed by our analyses on private alleles and particularly by the extremely low FST values between each founder subset and their source population. As expected, our results suggest that of the two CAP-WB stocks, the smaller CAP-WB-DON (N=9 pure and 4 backcross individuals) provides a slightly worse representation of its source than its counterpart from AND (N=36), as indicated by its higher proportion of unsampled alleles (Figure 1) and significant differences in NA and AR (Table 1). Thus, recruiting pure DON individuals should be a priority in the near future. On the other hand, once a source population is well represented genetically, as is now the case for AND, randomly selected new founders will likely make a relatively small contribution to the global genetic diversity and may rank low in breeding priority. Given this situation, together with current space limitations, carefully selecting potential founders based on their average kinship to the current captive population becomes a priority. A related but different matter concerns the proportions in which the two genetic stocks should be mixed in captivity, considering the difference in their genetic diversity and the level of genetic differentiation observed between them. Our estimated optimal contribution of about two AND lynx for every DON individual (0.64 vs. 0.36) has not quite yet been realised, as the AND stock has been overrepresented since 2003 to the detriment of the DON one. Empirical proportions did approach 80 CAPÍTULO 1 Caballero A, Toro MA. 2000. Interrelations between effective population size and other pedigree tools for the management of conserved populations. Genet. Res. 75(3):331–343. Calzada J, Guzmán JN, Rodríguez A. 2007. Lynx pardinus (Temminck, 1827). Ficha libro rojo. In: Palomo L, Gisbert J, Blanco J, editors. Atlas y libro rojo de los mamíferos terrestres de España. Madrid: Dirección General para la Biodiversidad- SECEM-SECEMU. p. 345–347. 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, et al. 2017. Spatiotemporal dynamics of genetic variation in the Iberian lynx along its path to extinction reconstructed with ancient DNA. Mol. Biol. Evol. 34(11):2893–2907. Casas-Marce M, Soriano L, López-Bao J, Godoy J. 2013. Genetics at the verge of extinction: insights from the Iberian lynx. Mol. Ecol. 22(22):5503–5515. Caughley G. 1994. Directions in conservation biology. J. Anim. Ecol. 63(2):215–244. Chen H, Boutros PC. 2011. VennDiagram: a package for the generation of highlycustomizable Venn and Euler diagrams in R. BMC Bioinformatics 12:35. Chesser RK. 1983. Isolation by distance: relationship to management of genetic resources. In: Schonewald-Cox C, Chambers S, MacBryde B, Thomas W, editors. Genetics and conservation: a reference for managing wild animal and plant populations. London & Amsterdam: Benjamin Cummings Publishing Co. p. 66–77. Conway WG. 1986. The practical difficulties and financial implications of endangered species breeding programmes. Int. Zoo Yearb. 24(1):210–219. Csilléry K, Johnson T, Beraldi D, Clutton-Brock T, Coltman D, Hansson B, Spong G, Pemberton JM. 2006. Performance of marker-based relatedness estimators in natural populations of outbred vertebrates. Genetics 173(4):2091–101. Dobrynin P, Liu S, Tamazian G, Xiong Z, Yurchenko AA, Krasheninnikova K, Kliver S, Schmidt-Küntzel A, Koepfli KP, Johnson W, et al. 2015. Genomic legacy of the African cheetah, Acinonyx jubatus. Genome Biol. 16:277. Fernández J, Toro MA. 1999. The use of mathematical programming to control inbreeding in selection schemes. J. Anim. Breed. Genet. 116(6):447–466. Fernández J, Toro MA. 2006. A new method to estimate relatedness from molecular markers. Mol. Ecol. 15(6):1657–1667. 81 GENETIC EVALUATION OF THE EX SITU POPULATION Fernández J, Toro MA, Caballero A. 2004. Managing individuals’ contributions to maximize the allelic diversity maintained in small, conserved populations. Conserv. Biol. 18(5):1358–1367. Frankham R. 2008. Genetic adaptation to captivity in species conservation programs. Mol. Ecol. 17(1):325–333. Frankham R. 2016. Genetic rescue benefits persist to at least the F3 generation, based on a meta-analysis. Biol. Conserv. 195:33–36. Frankham R, Ballou JD, Briscoe DA. 2010. Introduction to conservation genetics. Cambridge: Cambridge University Press. Frankham R, Ballou JD, Eldridge MDB, Lacy RC, Ralls K, Dudash MR, Fenster CB. 2011. Predicting the probability of outbreeding depression. Conserv. Biol. 25(3):465–475. Godoy JA, Casas-Marcé M, Fernández J. 2009. Genetic issues in the implementation of the Iberian lynx ex situ conservation programme. In: Vargas A, Breitenmoser C, Breitenmoser U, editors. Iberian lynx ex situ conservation: an interdisciplinary approach. Madrid: Fundación Biodiversidad in collaboration with: International Union for Conservation of Nature, Species Survival Commission (IUCN /SSC) Cat Specialist Group. p. 86–99. Gómez-Romano F, Villanueva B, de Cara MA, Fernández J. 2013. Maintaining genetic diversity using molecular coancestry: the effect of marker density and effective population size. Genet. Sel. Evol. 45(1):38. Gonçalves Da Silva A, Lalonde DR, Quse V, Shoemaker A, Russello MA. 2010. Genetic approaches refine ex situ lowland tapir (Tapirus terrestris) conservation. J. Hered. 101(5):581–90. Goudet J. 1995. FSTAT (version 1.2): a computer program to calculate F-statistics. J. Hered. 86(6):485–486. Goudet J. 2002. FSTAT, a program to estimate and test gene diversities and fixation indices (version 2.9.3). https://www2.unil.ch/popgen/softwares/fstat.htm. Guzmán JN, García FJ, Garrote G, Pérez de Ayala R, Iglesias C. 2004. El lince ibérico (Lynx pardinus) en España y Portugal. Censo diagnóstico de sus poblaciones. Madrid: Dirección General para la Biodiversidad. 82 CAPÍTULO 1 Henkel JR, Jones KL, Hereford SG, Savoie ML, Leibo SP, Howard JJ. 2012. Integrating microsatellite and pedigree analyses to facilitate the captive management of the endangered Mississippi sandhill crane (Grus canadensis pulla). Zoo Biol. 31(3):322–35. Hoffmann M, Hilton-Taylor C, Angulo A, Böhm M, Brooks TM, Butchart SHM, Carpenter KE, Chanson J, Collen B, Cox NA, et al. 2010. The impact of conservation on the status of the world’s vertebrates. Science 330(6010):1503–1509. Hogg CJ, Grueber CE, Pemberton D, Fox S, Lee AV, Ivy JA, Belov K. 2017. “Devil Tools & Tech”: a synergy of conservation research and management practice. Conserv. Lett. 10(1):133–138. Jiménez MA, Sánchez B, Alenza MDP, García P, López J V, Rodríguez A, Muñoz A, Martínez F, Vargas A, Peña L. 2008. Membranous glomerulonephritis in the Iberian lynx (Lynx pardinus). Vet. Immunol. Immunopathol. 121(1–2):34–43. Johnson WE, Godoy JA, Palomares F, Delibes M, Fernandes M, Revilla E, O’Brien SJ. 2004. Phylogenetic and phylogeographic analysis of Iberian lynx populations. J. Hered. 95(1):19–28. Jones KL, Glenn TC, Lacy RC, Pierce JR, Unruh N, Mirande CM, Chavez-Ramirez F. 2002. Refining the whooping crane studbook by incorporating microsatellite DNA and leg-banding analyses. Conserv. Biol. 16(3):789–799. Jost L. 2008. GST and its relatives do not measure differentiation. Mol. Ecol. 17(18):4015–4026. Kalinowski ST, Taper ML, Marshall TC. 2007. Revising how the computer program CERVUS accommodates genotyping error increases success in paternity assignment. Mol. Ecol. 16(5):1099–1106. Kalinowski ST, Wagner AP, Taper ML. 2006. ML-RELATE: a computer program for maximum likelihood estimation of relatedness and relationship. Mol. Ecol. Notes 6(2):576–579. Kardos M, Luikart G, Allendorf FW. 2015. Measuring individual inbreeding in the age of genomics: marker-based measures are better than pedigrees. Heredity (Edinb). 115(1):63–72. Kleinman-Ruiz D, Martínez-Cruz B, Soriano L, Lucena-Perez M, Cruz F, Villanueva B, Fernández J, Godoy JA. 2017. Novel efficient genome-wide SNP panels for the conservation of the highly endangered Iberian lynx. BMC Genomics 18:556. 83 GENETIC EVALUATION OF THE EX SITU POPULATION Lacy RC. 1987. Loss of genetic diversity from managed populations: interacting effects of drift, mutation, immigration, selection, and population subdivision. Conserv. Biol. 1(2):143–158. Lacy RC, Ballou JD, Pollak JP. 2012. PMx: software package for demographic and genetic analysis and management of pedigreed populations. Methods Ecol. Evol. 3(2):433–437. Laikre L. 1999. Hereditary defects and conservation genetic management of captive populations. Zoo Biol. 18(2):81–99. Lucena-Perez M, Soriano L, López-Bao JV, Marmesat E, Fernández L, Palomares F, Godoy JA. 2018. Reproductive biology and genealogy in the endangered Iberian lynx: implications for conservation. Mamm. Biol. 89:7–13. Lynch M. 1991. The genetic interpretation of inbreeding depression and outbreeding depression. Evolution 45(3):622–629. Marshall TC, Slate J, Kruuk LEB, Pemberton JM. 1998. Statistical confidence for likelihood-based paternity inference in natural populations. Mol. Ecol. 7(5):639–655. 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). J. Zoo Wildl. Med. 44(4):845–852. McGowan PJK, Traylor-Holzer K, Leus K. 2017. IUCN guidelines for determining when and how ex situ management should be used in species conservation. Conserv. Lett. 10(3):361–366. Meuwissen TH. 1997. Maximizing the response of selection with a predefined rate of inbreeding. J. Anim. Sci 75(4):934–940. Nei M. 1973. Analysis of gene diversity in subdivided populations. Proc. Natl. Acad. Sci. 70(12–I):3321–3323. Norman AJ, Putnam AS, Ivy JA. 2019. Use of molecular data in zoo and aquarium collection management: benefits, challenges, and best practices. Zoo Biol. 38(1):106– 118. Nowell K, Jackson P. 1996. Wild cats: status survey and conservation action plan. Gland: IUCN/SSC Cat Specialist Group. IUCN. 84 CAPÍTULO 1 Oliehoek PA, Windig JJ, Van Arendonk JAM, Bijma P. 2006. Estimating relatedness between individuals in general populations with a focus on their use in conservation programs. Genetics 173(1):483–496. Palomares F, Godoy JA, López-Bao JV, Rodríguez A, Roques S, Casas-Marce M, Revilla E, Delibes M. 2012. Possible extinction vortex for a population of Iberian lynx on the verge of extirpation. Conserv. Biol. 26(4):689–697. Peakall R, Smouse P. 2012. GenAlEx 6.5: Genetic analysis in Excel. Population genetic software for teaching and research – an update. Bioinformatics 28(19):2537– 2539. Peña L, Garcia P, Jiménez MA, Benito A, Pérez-Alenza MD, Sánchez B. 2006. Histopathological and immunohistochemical findings in lymphoid tissues of the endangered Iberian lynx (Lynx pardinus). Comp. Immunol. Microbiol. Infect. Dis. 29(2–3):114–126. Pimm SL, Jenkins CN, Abell R, Brooks TM, Gittleman JL, Joppa LN, Raven PH, Roberts CM, Sexton JO. 2014. The biodiversity of species and their rates of extinction, distribution, and protection. Science 344(6187):1246752. Pritchard DJ, Fa JE, Oldfield S, Harrop SR. 2012. Bring the captive closer to the wild: redefining the role of ex situ conservation. Oryx 46(1):18–23. R Core Team. 2020. R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. Rahbek C. 1993. Captive breeding – a useful tool in the preservation of biodiversity? Biodivers. Conserv. 2:426–437. Redford KH, Jensen DB, Breheny JJ. 2012. Integrating the captive and the wild. Science 338(6111):1157–1158. Rodríguez A, Calzada J. 2015. Lynx pardinus. The IUCN Red List of Threatened Species. 2015:e.T12520A174111773. Rodríguez A, Delibes M. 1992. Current range and status of the Iberian lynx Felis pardina Temminck, 1824 in Spain. Biol. Conserv. 61(3):189–196. Rodríguez A, Delibes M. 2002. Internal structure and patterns of contraction in the geographic range of the Iberian lynx. Ecography (Cop.). 25(3):314–328. 85 GENETIC EVALUATION OF THE EX SITU POPULATION Ruiz-López MJ, Gañán N, Godoy JA, Del Olmo A, Garde J, Espeso G, Vargas A, Martínez F, Roldán ERS, Gomendio M. 2012. Heterozygosity-fitness correlations and inbreeding depression in two critically endangered mammals. Conserv. Biol. 26(6):1121–1129. Russello MA, Amato G. 2007. On the horns of a dilemma: molecular approaches refine ex situ conservation in crisis. Mol. Ecol. 16(12):2405–2406. Ryder O, Miller W, Ralls K, Ballou JD, Steiner CC, Mitelberg A, Romanov MN, Chemnick LG, Mace M, Schuster S. 2016. Whole genome sequencing of California condors is now utilized for guiding genetic management. In: International Plant and Animal Genome XXIV Conference, 8-13 January 2016. San Diego, CA. Sambrook J, Russell DW. 2006. Purification of nucleic acids by extraction with phenol:chloroform. Cold Spring Harb. Protoc. 2006(1):pdb.prot4455. Schäfer F, Reiners TE. 2017. Long term vs short term impact of founder relatedness on gene diversity and inbreeding within the European Endangered Species Programme (EEP) for the Nepalese red panda (Ailurus f. fulgens). J. Zoo Aquarium Res. 5(2):86–91. Seal US, Foose TJ, Ellis S. 1994. Conservation Assessment and Management Plans (CAMPs) and Global Captive Action Plans (GCAPs). In: Olney PJ, Mace G, Feistner A, editors. Creative conservation: interactive management of wild and captive animals. Dordrecht: Springer Netherlands. p. 312–325. Simón MA, Gil-Sánchez JM, Ruiz G, Garrote G, Mccain EB, Fernández L, López- Parra M, Rojas E, Arenas-Rojas R, Rey T Del, et al. 2012. Reverse of the decline of the endangered Iberian Lynx. Conserv. Biol. 26(4):731–736. Snyder NFR, Derrickson SR, Beissinger SR, Wiley JW, Smith TB, Toone WD, Miller B. 1996. Limitations of captive breeding in endangered species recovery. Conserv. Biol. 10(2):338–348. Templeton A. 1986. Coadaptation and outbreeding depression. In: Soule M, editor. Conservation biology. The science of scarcity and diversity. Sunderland, MA: Sinauer Assoc. p. 105–121. The SSC Re-introduction Specialist Group. 1995. IUCN/SSC guidelines for reintroductions. In: Forty-first meeting of the IUCN Council. Gland. 86 CAPÍTULO 1 Tokarska M, Marshall T, Kowalczyk R, Wójcik JM, Pertoldi C, Kristensen TN, Loeschcke V, Gregersen VR, Bendixen C. 2009. Effectiveness of microsatellite and SNP markers for parentage and identity analysis in species with low genetic diversity: the case of European bison. Heredity (Edinb). 103(4):326–332. Toro MA, Caballero A. 2005. Characterization and conservation of genetic diversity in subdivided populations. Philos. Trans. R. Soc. B Biol. Sci. 360(1459):1367–1378. Vargas A, Breitenmoser C, Breitenmoser U. 2009. Iberian lynx ex situ conservation: an interdisciplinary approach. Madrid: Fundación Biodiversidad in collaboration with: International Union for Conservation of Nature, Species Survival Commission (IUCN /SSC) Cat Specialist Group. Vargas A, Rivas A, Sánchez I, Martínez F, Godoy JA, Roldán E, Serra RC, Pérez MJ. 2009. Iberian lynx conservation breeding program - update. Cat News 50:17. Vargas A, Sánchez I, Martínez F, Rivas A, Godoy JA, Roldán E, Simón MA, Serra R, Pérez MJ, Enseñat C, et al. 2008. The Iberian lynx Lynx pardinus conservation breeding program. Int. Zoo Yearb. 42(1):190–198. Wahlund S. 1928. Zusammensetzung von populationen und korrelationserscheinungen vom standpunkt der vererbungslehre aus betrachtet. Hereditas 11(1):65–106. Weir BS, Cockerham CC. 1984. Estimating F-statistics for the analysis of population structure. Evolution 38(6):1358–1370. Wilson EO. 1992. The diversity of life. Cambridge: Belknap Press. Witzenberger KA, Hochkirch A. 2011. Ex situ conservation genetics: a review of molecular studies on the genetic consequences of captive breeding programmes for endangered animal species. Biodivers. Conserv. 20(9):1843–1861. 87 GENETIC EVALUATION OF THE EX SITU POPULATION APÉNDICE CAPÍTULO 1 GENETIC EVALUATION OF THE IBERIAN LYNX EX SITU CONSERVATION PROGRAMME 96 APÉNDICE CAPÍTULO 1 0.018% [0.009,%0.028]% 0.013% [0.004,%0.022]% 0.010% [0.002,%0.019]% 0.004% [-0.004,%0.013]% -0.005% [-0.012,%0.001]% 0.006% [-0.003,%0.017]% Granadilla( Olivilla(Silves( Acebuche( Supplementary Figure S6. Diagram of the genetic structure among the main breeding facilities that comprise the Iberian lynx ex situ programme, measured as FST values [99% CI]. 97 GENETIC EVALUATION OF THE EX SITU POPULATION CAPÍTULO 2 NOVEL EFFICIENT GENOMEWIDE SNP PANELS FOR THE CONSERVATION OF THE HIGHLY ENDANGERED IBERIAN LYNX 100 CAPÍTULO 2 101 EFFICIENT GENOME-WIDE SNP PANELS FOR CONSERVATION ABSTRACT Introduction: The Iberian lynx (Lynx pardinus) has been acknowledged as the most endangered felid species in the world. An intense contraction and fragmentation during the 20th century left less than 100 individuals split in two isolated and genetically eroded populations by 2002. Genetic monitoring and management so far have been based on 36 STRs, but their limited variability and the more complex situation of current populations demand more efficient molecular markers. The recent characterization of the Iberian lynx genome identified more than 1.6 million SNPs, of which 1,536 were selected and genotyped in an extended Iberian lynx sample. Methods: We validated 1,492 SNPs and analysed their heterozygosity, Hardy-Wein- berg equilibrium, and linkage disequilibrium. We then selected a panel of 343 minimally linked autosomal SNPs from which we extracted subsets optimized for four different typical tasks in conservation applications: individual identification, parentage assignment, relatedness estimation, and admixture classification, and compared their power to currently used STR panels. Results: We ascribed 21 SNPs to chromosome X based on their segregation patterns, and identified one additional marker that showed significant differentiation between sexes. For all applications considered, panels of autosomal SNPs showed higher power than the currently used STR set with only a very modest increase in the number of markers. Conclusions: These novel panels of highly informative genome-wide SNPs provide more powerful, efficient, and flexible tools for the genetic management and non-in- vasive monitoring of Iberian lynx populations. This example highlights an important outcome of whole-genome studies in genetically threatened species. 102 CAPÍTULO 2 INTRODUCTION Genetic markers provide powerful tools for the study and conservation of biodiversity. They provide valuable insights into population sizes, mating systems, relatedness, population structure, and dispersal rates (Allendorf et al. 2010; Frankham et al. 2010), which can be effectively integrated into management and monitoring programmes of wildlife populations. Such programmes typically involve one or several of the following tasks: (i) individual identification, which allows for the estimation of census sizes, the identification of individual ranges and movements, and some forensic applications (Morin et al. 2004); (ii) parentage assignment, which is important in the study of mating systems and reproductive success and contributes to pedigree reconstruction (De Barba et al. 2010); (iii) relatedness estimation, which in addition to informing on important aspects of species biology, is of highest relevance for the genetic management of populations, as it forms the basis of the commonly used kinship minimization strategy (Caballero and Toro 2000); and (iv) detection of hybridization or admixture, as classification of individuals into specific ancestry categories might prove useful for comparing the fitness consequences of admixture (e.g. Miller et al. 2012). In addition, the genetic information gained can be used for the estimation of effective population size, genetic diversity, gene flow, and to infer demographic changes. Most importantly, all these tasks can be performed on genotypes obtained from non-invasively collected samples, avoiding capture and minimizing risks to animals and people (Goossens and Bruford 2009). Microsatellites or short-tandem repeats (STR) have been the markers of choice for the genetic monitoring and management of populations for the last two decades (Schlötterer 2004; Selkoe and Toonen 2006), but Single-Nucleotide Polymorphisms (SNPs) have been advocated as superior markers (Bensch et al. 2002; Belfiore et al. 2003). SNPs are single base nucleotide variants which make up a stable polymorphism in a species (Brookes 1999), and represent the most abundant and widespread source of sequence variation within genomes (Brumfield et al. 2003). Their appeal with respect to STRs also builds on: (i) their lower mutation rates (Estoup et al. 2002), which in turn imply less homoplasy (Morin et al. 2004), (ii) their lower expected error rates both in genotyping and allele calling (e.g. Quillery et al. 2014), and (iii) their automation potential (reviewed in Olsen et al. 2011). Finally, and contrary to what happens with STRs, SNP based assays are highly repeatable and easy to standardize across collaborators (Kraus et al. 2015). At the same time, their biallelic nature often implies lower resolution and statistical power per marker compared to the multiallelic STRs, but this can be easily counterbalanced by increasing the number of SNPs (Gärke et al. 2012). 103 EFFICIENT GENOME-WIDE SNP PANELS FOR CONSERVATION All in all, SNPs allow for a notorious increase in the power and accuracy of most genetic analyses when available in sufficient numbers. SNP discovery in non-model species, although still a major hurdle, has been greatly facilitated by recent advancements in Next Generation Sequencing (NGS) (Davey et al. 2011). Nanoscale genetic analyses on microfluidic platforms (Senapati et al. 2009) have streamlined, reduced costs, and added flexibility to the genotyping assays, minimizing at the same time the DNA quality and quantity requirements of chip-based technologies (Bhat et al. 2012; Kraus et al. 2015). Furthermore, the genotyping costs and the effort necessary to achieve the required power can be drastically reduced by selecting the most informative SNP for each particular application. Indeed, high-throughput SNP genotyping has been successfully applied to the non-invasive genetic monitoring of some wild mammal populations based on scat samples (Kraus et al. 2015; Norman and Spong 2015), and to genetic studies based on museum samples (Spitzer et al. 2016). The Iberian lynx (Lynx pardinus) has been highlighted as the most endangered felid species in the world (Nowell and Jackson 1996). Seven out of the nine populations extant in the 1980s became extinct by the end of the 20th century (Rodríguez and Delibes 2004), and the remaining two, Doñana (DON) and Andújar (AND), separated by ca. 240 km (Guzmán et al. 2004), have been effectively isolated at least since the 1950s (ca. 14 generations) until the start of translocations in 2007 (Simón et al. 2012). Several studies have revealed extremely low mitochondrial and STR diversity, record low genome-wide diversity, and excess of potentially deleterious variants, especially in the smallest population (DON), as well as high differentiation between the two populations (Johnson et al. 2004; Palomares et al. 2012; Casas-Marce et al. 2013; Abascal et al. 2016). Most importantly, evidence of inbreeding depression is accumulating (Peña et al. 2006; Jiménez et al. 2008; Meli et al. 2010; Ruiz-López et al. 2012) and the occurrence of an extinction vortex has been suggested for the DON population (Palomares et al. 2012). Early in the 21st century the Iberian lynx dramatic decline prompted the implementation of active conservation strategies with the main funding of three successive European LIFE projects (Simón et al. 2012). Conservation actions include (i) the close-up monitoring of individuals, (ii) an ex situ breeding programme that manages a captive population (CAP) resulting from the admixture of the two remnant differentiated genetic pools (DON and AND), (iii) genetic rescue through translocations between these two populations, and (iv) a reintroduction programme aimed at recovering the species’ historical distribution in Spain and Portugal. Population monitoring so far has mostly relied on radio-tracking and camera trapping (Garrote et al. 2011; Gil-Sanchez et al. 2011; Garrote et al. 2014), 104 CAPÍTULO 2 while marker-assisted genetic evaluation and management have been based on 36 STR markers (Godoy et al. 2009; Casas-Marce et al. 2013). The remarkable success of these actions has been endorsed by the recent downlisting of the species from “Critically Endangered” to “Endangered” in the 2015 International Union for Conservation of Nature (IUCN) Red List (Rodríguez and Calzada 2015), and also by the latest census (2015), which estimates a total of 406 free-living individuals in the Iberian Peninsula across all remnant and reintroduced populations (http://www. iberlince.eu/images/docs/3_InformesLIFE/Informe_Censo_2015.pdf). In spite of this, the species is still in high risk of extinction, and its survival heavily dependent on continued conservation measures. Given that in recent years many reintroduced populations have been founded by individuals from the admixed CAP population and several translocation events have been carried out between DON and AND, currently all Iberian lynx populations harbour varying degrees of genetic admixture. Such an increasingly complex and interwoven conservation landscape poses new challenges and demands tools as informative and efficient as possible for an integrative management and monitoring of the species. In this study, we analysed data for ca. 1,500 SNPs previously selected from whole genome sequences and previously genotyped in 380 Iberian lynx samples from five populations (Abascal et al. 2016). The resulting genotypes are analysed here for heterozygosity, Hardy-Weinberg equilibrium (HWE), and linkage disequilibrium (LD). We also identify sex chromosome markers and estimate genotyping error rates. Finally, we rank reliable autosomal SNPs by their informativeness for individual identification, parentage assignment, relatedness, and admixture estimation, and we identify optimal subsets of SNPs for each of these tasks. 105 EFFICIENT GENOME-WIDE SNP PANELS FOR CONSERVATION MATERIAL & METHODS SNP SELECTION, SAMPLES AND GENOTYPING Out of > 1.6 million SNPs identified through whole genome sequencing of 11 Iberian lynx males (Abascal et al. 2016), a total of 1,536 SNPs were selected based on their global Minor Allele Frequency (MAF ≥ 0.4), inter-SNP distance (d > 0.6 Mb), and genotyping quality score predicted by the Illumina’s Assay Design Tool algorithm (ADT > 0.9). A total of 380 Iberian lynx blood and muscle samples, including 28 replicates as internal controls, were genotyped for these SNPs using Illumina GoldenGate technology. Genotypes were successfully obtained for 349 samples corresponding to 329 individuals (and 20 replicates). Markers and samples with more than 20% of missing genotypes were identified with PLINK v1.90b3.36 (Purcell et al. 2007) and discarded from further analyses. To estimate an error rate for the final dataset we first verified replicate genotypes (Supplementary Table S1) with PLINK’s Identity By State (IBS) test, and then we tallied discordant genotypes within individual and divided them by the total number of actual pairwise comparisons using a custom R script. Replicate genotypes were discarded from further analyses. PEDIGREE REFINEMENT AND MENDELIAN INCONSISTENCIES With the purpose of refining the pedigree and, ultimately, identifying markers departing from autosomal Mendelian inheritance (including X-linked markers), we assessed Mendelian inconsistencies based on known parent-offspring (PO) relationships using the --mendel function in PLINK. Some of these discrepancies allowed us to identify a few misassignments and amend the current pedigree. Additionally, some previously unknown relationships were uncovered by maximum likelihood analyses using ML-RELATE (Kalinowski et al. 2006) and COLONY (Jones and Wang 2010). Markers in the refined version of the pedigree that systematically yielded Mendelian incompatibilities between sires and their male offspring were flagged as X-linked markers and treated separately from the main autosomal database unless otherwise specified. Next, we performed the Fisher’s exact test in R (R Core Team 2020) to assess differences in allele and genotype frequencies between sexes, which could be indicative of sex-linked inheritance or just genotyping artefacts. Finally, Mendelian inconsistencies in PO dyads and triads in the filtered autosomal database were quantified as a proxy for error rate. 112 CAPÍTULO 2 S S e e t t N N u u m m b b e e r r o o f f S S N N P P s s A A n n c c e e s s t t r r y y N N u u m m b b e e r r o o f f m m o o n n o o m m o o r r p p h h i i c c S S N N P P s s H H E E H H O O F F S S T T ( ( D D O O N N - - A A N N D D ) ) F F I I S S PDON 471 0.235 0.240 -0.013 PAND 5 0.398 0.408 -0.022 PDON 7 0.216 0.236 - PAND 0 0.368 0.352 - A A u u t t o o s s o o m m a a l l 1,471 0.28 X X - - l l i i n n k k e e d d 21 0.309 Table 1. Number of autosomal and X-linked SNPs per ancestry, and average variability and differentiation statistics. HE: expected heterozygosity; HO: observed heterozygosity; FIS: population inbreeding coefficient; FST: Weir & Cockerham’s fixation index; DON: Doñana; AND: Andújar; PDON: pure DON; PAND: pure AND. 113 EFFICIENT GENOME-WIDE SNP PANELS FOR CONSERVATION ● ● ● ● ● ● ● ● ● ● ● ● 6×10−11 3.7 ×10−6 ● ●● ● 5.1 ×10−10 9.7 ×10−5 10−120 10−80 10−40 100 Probability of identity PID PIDs a ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●4.7 ×10−5 1.7 ×10−3 1.3 ×10−7 ● ●● ● ● ● 5.6 ×10−6 2.2 ×10−3 1.6 ×10−9 10−40 10−30 10−20 10−10 100 Probability of non−exclusion PnE1 PnE2 PnE3 b ● ● ● ● ● ● ● ● ● ● ● ● ● ● 6.8 ×10−5 5.5 ×10−2 1.8 ×10−1 ● ● ● ●● ● 7.5 ×10−3 2.5 ×10−1 3.3 ×10−1 12 24 48 96 192 343 12 36 10−6 10−4 10−2 100 Subset size 1 − PWR 1−PWRFS_HS 1−PWRFS_UR 1−PWRHS_UR c SNPs STRs Figure 1. Power of different-sized panels of SNPs and STRs for individual identification (a), parentage exclusion (b), and relationship discrimination (c). White circles indicate the smallest panel informative enough for the respective application, and their power is indicated. 1 – PWR value for the FS vs. UR comparison with 343 SNPs is not shown due to KININFOR’s lack of sensitivity beyond PWR = 0.9999991 (DON: Doñana population; AND: Andújar population; PDON: pure DON; PAND: pure AND). 114 CAPÍTULO 2 Admixture and ancestry estimation. Genetic markers can be used to classify individuals in each of the possible ancestry or admixture categories generated following the mixing of the two lynx genetic pools (Supplementary Figure S3). For this application, we ranked markers according to their Ir and their FST between PDON and PAND. Then, using NEWHYBRIDS, we estimated the power of differently-ranked and -sized panels to classify individuals in ancestry categories as the rate of correct classifications of simulated genotypes generated with HYBRIDLAB (Supplementary Table S17). Generally speaking, the 343 SNP set was the best performing set. It showed the same power as the set of 36 STRs for the categories of easier discrimination (PDON: 98% vs. 100%; PAND: 94% vs. 94%; F1: 100% vs. 99%) and starkly outperformed it for the remaining categories (F2: 98% vs. 86%; BcDON: 94% vs. 69%; BcAND: 88% vs. 77%; 2BcDON: 87% vs. 78%; 2BcAND: 84% vs. 50%). Notwithstanding, the performance of reduced SNP panels of highest Ir was not satisfactory (Supplementary Table S17). In contrast, panels based on the rank of FST between PDON and PAND performed comparatively much better. Results for the 96 SNP set show a perfect classification of F1 individuals and a highly accurate (>95%) classification of pure and F2 individuals, whereas correct assignments were slightly lower for first backcross (BcDON: 86%; BcAND: 84%) and second backcross categories (2BcDON: 85%; 2BcAND: 78%). Only 33 out of the 800 simulated individuals couldn’t be unambiguously assigned to any ancestry level. For smaller panels the classification accuracy fell as expected with size, and particularly so for the less distinguishable backcross and second backcross categories. Overall, the performance of the panel of 48 SNPs was similar to that of the 36 STRs (Figure 2). The results of the analyses of empirical data only partially confirmed the patterns observed for simulated data. Larger (SNP) panels showed higher rates of misclassification of pure individuals as second backcrosses than expected from simulations, whereas a better correlation was observed for smaller panels. The sets of 48 SNPs and 36 STRs behaved similarly to each other and to expectations from simulations, with a slightly lower classification of pure DON individuals with SNPs (Supplementary Table S18). All F1 individuals were perfectly classified across all panels. Unfortunately, admixed categories beyond F1 are poorly represented or completely absent from the empirical data, preventing a more thorough comparison of simulated and empirical results. 115 EFFICIENT GENOME-WIDE SNP PANELS FOR CONSERVATION ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 20 40 60 80 100 PDON PAND F1 F2 BcDON BcAND 2BcDON 2BcAND Ancestry levels Correct classification rate (%) Panels ● ● ● ● ● ● ● ● ● ● ● ● 343 SNPs 192 SNPs 96 SNPs 48 SNPs 24 SNPs 12 SNPs 36 STRs 12 STRs Figure 2. Rate of correct classification of simulated populations with eight distinct degrees of admixture (ancestry levels) for different-sized SNP and STR panels of highest FST markers (PDON: pure DON; PAND: pure AND; F1: filial generation of PDON x PAND cross; F2: offspring of F1 x F1 cross; BcDON: first backcross DON; BcAND: first backcross AND; 2BcAND: second backcross AND). 116 CAPÍTULO 2 DISCUSSION Our extensive curation and evaluation of a dataset of genotypes for 1,536 SNPs allowed us to select maximally informative autosomal SNP panels for individual identification, parentage assignment, relatedness estimation, and admixture estimation. Proposed panels showed higher power than the currently used 36 STR set with only a very modest increase in the number of markers, while adding at the same time the many benefits inherent to SNPs (e.g., their applicability to low quality samples). In combination with flexible genotyping technologies they do provide novel, efficient and cost-effective tools for the demographic and genetic monitoring and management of this highly endangered and actively managed species. Furthermore, we identified 21 X-linked markers that may prove useful for complementing parentage analyses, determining sex, or genetic studies focused on this chromosome. Genetic monitoring and most forensic applications involve a first step in which genotypes are obtained from new samples and, by checking them against the reference database, are either assigned to already recorded individuals, or identified as previously unsampled individuals. The proposed panel of 24 SNPs has sufficient power to discriminate even between closely related individuals. In a subsequent step in genetic monitoring, novel individuals may be assigned to parents so that a population genealogy is progressively built (e.g. De Barba et al. 2010). We found that a panel of only 48 SNPs provides sufficient power to assign parentage with high confidence, even in the most demanding situation where none of the parents are known. It must be noted, though, that more markers (e.g. 96 SNPs) could be needed in specific scenarios with multiple, closely related putative sires (Fernández et al. 2013), as observed in a previous study on reproduction patterns in the DON population (Lucena-Perez et al. submitted). One important feature of genetic monitoring is that it can be based on non-invasively collected samples, such as faeces, hairs or feathers (Schwartz et al. 2007), and here SNPs offer important advantages over STRs (Morin et al. 2004; Seddon et al. 2005). Degraded DNA in low quality samples is prone to low rates of amplification and high rates of genotyping errors, including false alleles and allelic dropout (Bonin et al. 2004). Reported error rates for STR markers vary greatly among studies, ranging from 0% to above 30% when considering just allelic dropout (reviewed in Broquet et al. 2007), thus imposing the need for multiple genotyping replicates of the same extract (Navidi et al. 1992). However, degraded DNA can be typed more efficiently 117 EFFICIENT GENOME-WIDE SNP PANELS FOR CONSERVATION and accurately for SNPs due to their considerably shorter target DNA sequence (~50-70 bp) compared to STRs’ (80-300 bp) (Morin et al. 2004). Indeed, the few existing studies on SNP genotyping of non-invasively collected wildlife samples report low rates of genotyping failure and errors, e.g. <10% and ~1%, respectively, in wolves (Kraus et al. 2015), and 0.36% and 0.038% in bears (Norman and Spong 2015), eliminating the need for costly and laborious systematic replicates. The other very important use of genetic markers in conservation is genetic management, which is typically based on a kinship minimization strategy (Caballero and Toro 2000). This strategy is presently being applied to the management of the captive Iberian lynx population using STR-based estimations of relatedness among founders in combination with the completely known genealogy (Godoy et al. 2009). However, the currently used 36 STR set in fact provides limited information on relatedness, since it correctly discriminates <80% –a threshold used by Wang (2006)– of FS from HS dyads at a significance level of p < 0.05. On the contrary, 96 SNPs provide enough power to correctly discriminate >80% of HS from UR dyads, which is the most demanding of the three tested comparisons. Still, even our whole 343 SNP panel does not provide enough density to completely replace genealogy in kinship minimization strategies, which has been estimated in 3Ne SNPs/Morgan (Gómez- Romano et al. 2013). On the other hand, even when quite complete genealogies are available, genetic markers can be used to estimate the unknown relatedness among the founders. Furthermore, genetic markers do estimate the realized –rather than the expected– kinship or inbreeding (which can be more useful measures for most applications) for the whole genome or for specific regions (Gómez-Romano et al. 2013; Kardos et al. 2015). The increasingly complex situation of Iberian lynx conservation demands the extension of the genetic management to the remnant and reintroduced populations, ideally within a single integrative genetic management programme for the species. Relatedness estimation based on the novel SNP panels presented here will become most critical for this objective, since with the exception of DON (Lucena-Pérez et al. submitted), other populations lack reliable and reasonably complete genealogies. The fact that PID, PIDs, PnE, and Ir share a nearly identical ranking of markers means that panels designed for a less demanding purpose (e.g. individual identification) are nested within those for more demanding applications (e.g. parentage assignment). This enables a module workflow that grants great flexibility to genetic monitoring and management of Iberian lynx populations. For example, non-invasively collected samples would be assigned to specific individuals based on the 12-24 SNP panels. The 118 CAPÍTULO 2 novel unique genotypes could then be assigned to parents and incorporated to the pedigree after genotyping one of its replicate samples for the next most informative 24-72 SNPs. Finally, individual lynx that are candidates for being translocated or incorporated into captivity could be additionally typed for the remaining SNPs for relatedness estimation. We also considered admixture estimation as a potential application of the novel SNP markers, since the two highly differentiated genetic pools have recently been admixed in captivity and in the wild through translocations. Differentiating individuals by ancestry might be most relevant for the assessment of the fitness consequences of admixture, i.e. for the evaluation of genetic rescue or outbreeding depression (e.g. Miller et al. 2012). Since most of the markers of highest HE are poorly differentiated between the two pure ancestries, admixture tests would require their own particular sets of ancestry informative markers, with little overlap with those selected for other applications. For this purpose, using loci with different alleles fixed in each population would be optimal, but these were excluded by the MAF requirement in the initial selection. Fortunately, a sufficient number of highly differentiated SNPs were available to provide enough power for ancestry classification. Simulated individuals of PDON, PAND, F1, and F2 categories are reasonably well classified with the 48 highest FST SNP set and the 36 STR set, whereas less differentiated categories require at least 96 SNPs. The observed accuracy of ancestry classification of empirical genotypes with SNP panels became lower than that inferred from simulated genotypes as panel size increased. This is probably due to the additional LD derived from the relatively high relatedness structure and, especially, the recent admixture in current Iberian lynx populations (admixture linkage disequilibrium) (Loh et al. 2013). This effect gains relevance as the density of markers increases, and is expected to dilute out with successive generations of random mating following admixture (Pfaff et al. 2001). This should also result in the overestimation of power of individual discrimination, parentage assignments, and relationship inference, which assume markers combined in the same panel are completely unlinked. However, individual identification and parentage assignment shouldn’t be much affected by this, since i) the proposed panels are small, and ii) the threshold of estimated power that we chose for these applications is well above the minimum required, an advice that could be extended to similar studies in other species to account for any possible deviation of theoretical or simulation assumptions. Finally, even if not completely realistic, simulations and theoretical power estimations are still useful for the comparison among marker types and panel sizes. 119 EFFICIENT GENOME-WIDE SNP PANELS FOR CONSERVATION Our observed SNPs to STRs ratio (e.g. around 1.3 for individual identification) is at the lower end of or even below the range of 1.7 to 5.56 SNPs per STR reported in humans (see Gärke et al. 2012 and the references therein). The high performance of SNPs when compared to STRs found here can be explained by the restrictive selection of the most informative among the more than 1.6 million SNPs found in the analyses of 11 Iberian lynx genomes (Abascal et al. 2016). In contrast, the 36 STRs were selected among far fewer markers, which were also originally developed for other felid species: 25 out of ~250 from domestic cat, seven out of eight from bobcat (Lynx rufus), and four out of six from Canada lynx (Lynx canadensis) (Casas-Marce et al. 2013). This more limited search and the ascertainment bias associated to the heterologous nature of these markers have limited the variability of the selected STRs, to the point where some of the worst performing STRs have lower HE than any of the selected SNPs. These results highlight the importance of the variant discovery and screening phases. Once considered a major limitation of SNP adoption, NGS approaches have greatly facilitated these tasks in non-model organisms, especially when applied to reduced genome subsets or transcriptomes (Kumar et al. 2012). Still, such screening effort will be especially challenging for genetically eroded species, requiring a more extensive sampling of the genome (as exemplified here) or yielding larger panels for any given power than in genetically diverse species. The decision on how many markers to use for the most demanding applications (e.g. relatedness estimation) will mostly rely on cost-budget considerations, and these heavily depend on the genotyping technology used. The excessive cost of chip-based technologies –but also their limited flexibility and high input sample requirement– have deterred the application of non-invasive SNP-based approaches in wildlife. On the other hand, technologies based on fluorescent detection provide the advantage of a single PCR reaction, thus avoiding multiple manipulation steps (Fabbri et al. 2012), and when implemented with nanofluidic systems, the sample, reagent, and labour requirements are substantially reduced. For example, the Fluidigm system enables flexible and high-throughput SNP genotyping through the use of dynamic arrays with Integrated Fluidic Circuits (IFCs) of different sizes, such as 48.48 (48 samples against 48 markers per run) or 192.24 (192 samples against 24 markers). The philosophy behind this technology fits well the needs of a conservation programme, where samples may not accrue at a uniform pace, and questions requiring different SNP sets might need to be addressed as soon as they arise. Moreover, hands-on laboratory costs with the Fluidigm system are substantially cheaper than for STR genotyping, and initial investments in synthesis of assay probes are now similar – but 120 CAPÍTULO 2 SNP probes can be applied to a much larger number of assays (Kraus et al. 2015). Thus, in the long term, the Fluidigm system stands as an efficient and cost-effective option for continued monitoring of wildlife. The Iberian lynx example highlights one valuable practical outcome of whole-genome studies in heavily eroded species, i.e. the identification of scarce highly informative markers that will most likely be missed in smaller scale screenings. Indeed, highly efficient and cost-effective marker panels for genetic monitoring and management are likely the most immediate and feasible contribution of genomics to endangered species conservation. The Iberian lynx is a quite extreme example of genetically eroded and intensively managed species, with ongoing conservation actions including captive breeding, translocations and reintroductions. The long-term success of such actions, and ultimately the viability of the species, will largely depend on the implementation of sound, efficient, and science-driven monitoring and management programmes. As the novel SNP markers presented here provide higher power, efficiency, and flexibility than currently used STRs, they can make a substantial contribution toward this goal by guiding the comprehensive genetic management of captive, remnant and reintroduced populations and, in combination with non-invasive sampling, by complementing –and eventually substituting– monitoring programmes currently based almost exclusively on radiotracking and camera-trapping. 121 EFFICIENT GENOME-WIDE SNP PANELS FOR CONSERVATION ACKNOWLEDGMENTS Funding for this project was provided by the Spanish Dirección General de Investigación Científica y Técnica (CGL2013-47755-P), and Banco de Santander, CSIC and Fundación General CSIC through a “Proyectos Cero” grant. DKR and MLP were supported by PhD contracts from Programa Internacional de Becas “La Caixa-Severo Ochoa”. FC and BMC were hired under the EcoGenes project funded by the European FP7 programme (264125-FP7-REGPOT-2010-1). EBD-CSIC received support from the Spanish Ministry of Economy and Competitiveness under the ‘Centro de Excelencia Severo Ochoa 2013-2017’ program, grants SEV- 2012-0262. Iberian lynx samples were contributed by the Iberian lynx Ex Situ Breeding Programme, the Centro de Análisis y Diagnóstico de la Fauna Silvestre and the LIFE projects LIFE06NAT/E/000209 and LIFE10NAT/ES/570, with authorization of Consejería de Medio Ambiente de la Junta de Andalucía. Logistical support was provided by Laboratorio de Ecología Molecular, Estación Biológica de Doñana, CSIC (LEM-EBD). We appreciate useful comments and suggestions from Ana Píriz, Elena Marmesat and two anonymous reviewers. DATA ACCESSIBILITY STATEMENT The datasets analysed in this study are available in our institutional repository, DIGITAL.CSIC: http://digital.csic.es/handle/10261/153267. 128 CAPÍTULO 2 Waits LP, Luikart G, Taberlet P. 2001. Estimating the probability of identity among genotypes in natural populations: cautions and guidelines. Mol. Ecol. 10(1):249–256. Wang J. 2006. Informativeness of genetic markers for pairwise relationship and relatedness inference. Theor. Popul. Biol. 70(3):300–321. 129 EFFICIENT GENOME-WIDE SNP PANELS FOR CONSERVATION APÉNDICE CAPÍTULO 2 NOVEL EFFICIENT GENOMEWIDE SNP PANELS FOR THE CONSERVATION OF THE HIGHLY ENDANGERED IBERIAN LYNX 132 APÉNDICE CAPÍTULO 2 133 EFFICIENT GENOME-WIDE SNP PANELS FOR CONSERVATION SUPPLEMENTARY FIGURES Figure S1. Location of syntenic SNPs in the domestic cat’s chromosomes. Figure S2. Family trees with Mendelian errors. Figure S3. Ancestry of Iberian lynx individuals inferred under the assumption of two genetic clusters. SUPPLEMENTARY TABLES * Table S1. Genotypes of replicate samples. Table S2. Results of the Fisher’s exact test on differences in allele and genotype frequencies between sexes for the pure DON pool. Table S3. Results of the Fisher’s exact test on differences in allele and genotype frequencies between sexes for the pure AND pool. Table S4. Basic population genetics statistics per autosomal marker. Table S5. Basic population genetics statistics per X-linked marker. Table S6. Results for HWE analysis per autosomal marker for the pure DON pool. Table S7. Results for HWE analysis per autosomal marker for the pure AND pool. Table S8. Results for HWE analysis per X-linked marker for the pure DON pool. Table S9. Results for HWE analysis per X-linked marker for the pure AND pool. Table S10. Hemimatrix of pairwise r2 values between autosomal markers for the pure DON pool. Table S11. Hemimatrix of pairwise r2 values between autosomal markers for the pure AND pool. Table S12. Hemimatrix of pairwise r2 values between X-linked markers for the pure DON pool. Table S13. Hemimatrix of pairwise r2 values between X-linked markers for the pure AND pool. 134 APÉNDICE CAPÍTULO 2 Table S14. Hemimatrix of Kendall rank correlation τ values between HE and statistics selected to rank markers. Table S15. Informativeness and rankings for individual identification, parentage assignment, and relatedness estimation of SNPs in the 343 SNP panel. Table S16. Informativeness and rankings for individual identification, parentage assignment, and relatedness estimation of STRs in the 36 STR panel. Table S17. Summary matrices with NEWHYBRIDS’ correct, ambiguous, and cross-classification rates of simulated individuals in eight ancestry levels. Table S18. Summary matrices with NEWHYBRIDS’ correct, ambiguous, and cross-classification values (and percentages) of empirical individuals in eight ancestry levels. * Due to the their huge volume, supplementary tables are provided in digital format and can be accessed from https://bmcgenomics.biomedcentral.com/articles/10.1186/ s12864-017-3946-5#Sec18. 135 EFFICIENT GENOME-WIDE SNP PANELS FOR CONSERVATION Supplementary Figure S1. Location of syntenic SNPs in the domestic cat’s chromosomes. 136 APÉNDICE CAPÍTULO 2 Figure S2. Family trees accruing all eight Mendelian errors concerning SNP 1317372 (A); SNP 334649 (B); SNP 619378 (C); SNP 2119932 (D); and SNP 2057563 (E). Dark background: allele A; light background: allele B; white background: unknown allele.          Supplementary Figure S2. Family trees accruing all eight Mendelian errors concerning SNP 1317372 (A); SNP 334649 (B); SNP 619378 (C); SNP 2119932 (D); and SNP 2057563 (E). Dark background: allele A; light background: allele B; white background: unknown allele. 137 EFFICIENT GENOME-WIDE SNP PANELS FOR CONSERVATION 2BcAND F2 BcDON Supplementary Figure S3. Ancestry of Iberian lynx individuals in each of the two differentiated genetic pools, as inferred from the analysis of 1,471 autosomal SNPs in STRUCTURE under the assumption of two genetic clusters (K = 2). Each individual is represented by a vertical line with the height of each colour representing the estimated fraction of their genome belonging to each cluster. The black vertical lines delimit groups of individuals of different admixture categories, know a priori from their provenance or recent genealogy (PDON: pure DON; PAND: pure AND; F1: filial generation of PDON x PAND cross; F2: offspring of F1 x F1 cross; BcDON: first backcross DON; BcAND: first backcross AND; 2BcAND: second backcross AND). Ten runs of 100,000 MCMC iterations (of which the first 10,000 were discarded as burn-in) were performed in STRUCTURE v2.3.4 with highly concordant results. A representative run is depicted here. CLUMPAK main pipeline - Job 1606237876 summary Major modes for the uploaded data: K=2 Minor modes for the uploaded data: Division of runs by mode: K=2 1/1 PDON BcAND PAND F1 144 CAPÍTULO 3 free-living lynxes in 2018 (http://www.iberlince.eu/images/docs/3_InformesLIFE/ Informe_Censo_2018.pdf), and over the (estimated) milestone of 1,000 in 2020. Genetic erosion in the species is well documented, with studies on mitochondrial, microsatellite and genome-wide data providing evidence of a high inbreeding rate, high differentiation between the two populations, and a whole-genome species-wide nucleotide diversity (π = 0.026%) that is among the lowest reported so far (Palomares et al. 2012; Casas-Marce et al. 2013; Abascal et al. 2016). Furthermore, there are documented signs of high inbreeding load in the species (Ruiz-López et al. 2012) and inbreeding depression in the DON population (Palomares et al. 2012). On the other hand, the Eurasian lynx is one of the most widespread felids in the world, with a distribution area that extends from Eastern Asia to Central Europe, and covers a wide range of habitats. It is currently classified as “least concern” by the IUCN at the global level. However, recently inferred demographic trajectories have revealed a generalised decline across its range during the last millennia (Lucena- Perez et al. 2020). This decline has been particularly intense in the European side as a consequence of higher anthropogenic pressure, to the point where many of the westernmost populations became extirpated during the 19th and 20th centuries, and some of the remaining ones, such as those in Northeast Poland (POL) and Norway (NOR), underwent severe bottlenecks and complete genetic isolation for decades (Jedrzejewski et al. 1996; Linnell et al. 2010). Meanwhile the Kirov (European Russia) population (KIR), which is part of a large continuous patch of the species in Eastern Europe, has maintained larger size and good connectivity with neighbouring populations (Ratkiewicz et al. 2014; Lucena-Perez et al. 2020). Even so, genetic diversity for KIR (π = 0.046%) is barely twice as high as that of the Iberian lynx, and still below that of severely bottlenecked populations of other species, such as the Apennine brown bear, the Scandinavian wolverine, or the Amur tiger (Lucena-Perez et al. 2020). The contrasting demographic histories of the Iberian and the Eurasian lynx, coupled with the availability of a reference genome for the former species, showcases this dyad as a suitable model for carrying out genome-wide studies of the dynamics of genomic load in bottlenecked populations. We specifically intend to address four questions: i) whether there are differences in the amount or the patterns of genomic load between the Iberian and the Eurasian lynx; ii) whether such differences extend to populations within each species; iii) if differences in genomic load are driven by the relaxation of purifying selection or by genetic purging; and iv) how do these differences vary across genomic features, mutation effect categories, and genomic regions under contrasting efficiency of natural selection. Beyond broadening our understanding of genetic processes at the verge of extinction, such insights could be critical to predict the fate of small lynx populations, and to improve the efficacy of on-going conservation programmes. 145 A COMPARATIVE ANALYSIS OF GENOMIC LOAD MATERIALS & METHODS SAMPLING Our sample included 31 Iberian lynx individuals from the last two remnant populations of Andújar-Cardeña (AND; N=19) and Doñana-Aljarafe (DON; N=12), in Southern Spain, before they were reconnected through translocations and reintroductions, and 29 Eurasian lynx from the Russian Kirov region (KIR; N=13), Norway (NOR; N=8) and Białowieża and Knyszyn Primeval Forests in North-Eastern Poland (POL; N=8). Additionally, we included one individual of bobcat (Lynx rufus), from the Jerez Zoo (Spain), to be used as outgroup for the identification of the ancestral and derived states of the detected variants. A summary of samples is provided in Supplementary Table S1. For details on the demography of the Iberian and Eurasian lynx populations beyond those provided in the Introduction, please refer to Casas-Marce et al. (2017) and Lucena-Perez et al. (2020), respectively, and to Abascal et al. (2016) for additional information on the divergence and demographic trajectories of the two species. GENOMIC DATA We included both newly sequenced samples (as detailed in the next section) and previously reported whole-genome sequences, and arranged them in four datasets (Supplementary Table S1) for depth filtering purposes. The Iberian lynx highcoverage (ILH) dataset, previously reported in Abascal et al. (2016), comprised a total of six AND and four DON individuals originally re-sequenced at a depth of ~26X (range: 24.3X to 28.4X) using the Illumina GAIIx platform. We included an additional individual from AND, the one providing the reference genome, which was subsampled to a depth in the range of the re-sequenced individuals (Abascal et al. 2016). The Eurasian lynx high-coverage (ELH) dataset consisted of one sample from KIR sequenced at a depth of 22X using Illumina HiSeq X-10, as detailed in Lucena-Perez et al. (2020). All remaining individuals were re-sequenced at a depth in the range of 5-7X using Illumina HiSeq2000. Of these, 20 newly sequenced samples (12 from AND and eight from DON) formed the Iberian lynx low-coverage (ILL) dataset, while the remaining 28 samples (12 from KIR, eight from NOR and eight from POL) composed the Eurasian lynx low-coverage (ELL) dataset. Wholegenome sequences of all Eurasian lynx samples, whose intergenic sequences were 146 CAPÍTULO 3 previously published in Lucena-Perez et al. (2020), are reported here for the first time. To avoid any possible batch effect, including depth and sequencing technology biases, for most analyses we only considered individuals from the ILL and ELL datasets. Finally, we also included the one bobcat individual re-sequenced at the whole-genome level in Lucena-Perez et al. (2020). DNA EXTRACTION, LIBRARY PREPARATION & SEQUENCING DNA was extracted from good quality tissue or blood samples (which were digested overnight using proteinase K) using NucleoMag B-beads (NucleoMag DNA from tissue kit) in the LEM-EBD (Seville, Spain). For each sample, 2.0 micrograms of genomic DNA were sheared on a Covaris™ E210, and then end-repaired, adenylated and ligated to Illumina specific indexed paired-end adaptors using the NO-PCR protocol, TruSeq™ DNA Sample Preparation Kit v2 (Illumina Inc.), and the KAPA Library Preparation kit (Kapa Biosystems). Then, AMPure XP beads (Agencourt, Beckman Coulter) were used to select fragments with insert size of 220-550bp. Library Quantification Kit (Kapa Biosystems) was used to quantify the final libraries, which were sequenced at Centro Nacional de Análisis Genómico (CNAG-CRG; Barcelona, Spain) in a fraction of a sequencing lane of HiSeq2000 flowcell v3 (Illumina Inc.), in paired end mode 2x101bp, and using TruSeq SBS Kit v3-HS (Illumina Inc.). Sequencing, primary data and image analysis, base calling, and quality scoring of the runs all followed the standard Illumina procedures. READ PROCESSING AND MAPPING Adapter trimming was performed using SeqPrep (https://github.com/jstjohn/SeqPrep) before reads were mapped to the Iberian lynx reference genome using BWAMEM (Li 2013) with default parameters. Duplicates were marked with MarkDuplicates from Picard (https://broadinstitute.github.io/picard/ v2.4.1). GATK v3.4 (McKenna et al. 2010) was then used to perform a local realignment and a Base Quality Score Recalibration (BQSR). To minimise the risk of identifying potential contaminant variants as endogenous, we used bamtools filter (Barnett et al. 2011) to discard reads with an edit distance (NM) higher than 2 (NM ≥ 3) and 3 (NM ≥ 4) in Iberian lynx and Eurasian lynx samples, respectively. See the Supplementary Material for more detailed information on read processing and mapping. 147 A COMPARATIVE ANALYSIS OF GENOMIC LOAD VARIANT CALLING AND POLARISATION GATK v3.7 HaplotypeCaller (McKenna et al. 2010) was used to generate variation data for the pool of 60 Iberian and Eurasian lynx samples. Repetitive, low complexity and low mappability regions, which together accounted for 46.24% of the genome (Abascal et al. 2016), were excluded from the calling. We imported the ancestral allele annotation from the bobcat genome into the vcf file using fill-aa from VCFtools v0.1.13 (Danecek et al. 2011) as a stepping-stone to the polarisation of variants into the ancestral and derived states (Supplementary Material), for which we employed VcfFilterJdk (Lindenbaum and Redon 2018). VARIANT FILTERING A stepwise filtering framework was designed to keep only high confidence biallelic SNPs. First, INDELs and multiallelic SNPs (with three or more alleles) were removed from the dataset using GATK v3.7 SelectVariants (McKenna et al. 2010). Second, derived positions that were fixed in both species, as well as variants that could not be polarized (i.e. those lacking the ancestral allele information), were filtered out using BCFtools (Li 2011). Third, GATK v3.7 SelectVariants was again used to apply hard-filters following the GATK best practices manual (Van der Auwera et al. 2013). Standard thresholds were accepted for all filtering statistics after inspection of their genome-wide distributions; in detail, variants were discarded if they met any of the following criteria: QUAL < 30.0, QD < 2.0, FS > 60.0, MQ < 40.0, MQRankSum < -12.5, and ReadPosRankSum < -8.0. Fourth, to exclude sites that were covered unequally across the different sequencing technologies or the two species, and those that were poorly covered in general, the command filter –e from BCFtools was applied with the flags F_MISSING > 0.15 & DP < 200. By doing so, only variants that were genotyped in at least 85% of all individuals and showed high overall depth were kept. Finally, all polymorphic sites with depth exceeding the average plus three times the standard deviation in any of the four datasets of samples (i.e. ILH, ILL, ELH, ELL) were filtered out from the global dataset to discard false variants caused by collapsed paralogs. For this purpose, depth distribution data was computed using ANGSD (Korneliussen et al. 2014). Additionally, we used GATK v3.7 SelectVariants to extract for each species the set of locally fixed variants, which were eventually used for fixation rate analyses, and the set of segregating variants, used to build population allele count histograms. 148 CAPÍTULO 3 VARIANT ANNOTATION We used SnpEff v4.3i (Cingolani et al. 2012) to annotate variants based on an annotation file that included only the principal isoforms, as identified in the Iberian lynx genome project (Abascal et al. 2016). This file was subsequently modified to incorporate intergenic and intronic regions (which were delimited by exclusion criteria), as well as ultra conserved non-coding elements (UCNEs), whose annotation was carried from the human genome assembly hg19 (Dimitrieva and Bucher 2013) to cat coordinates (assembly felCat5) using liftOver (http://rohsdb.cmb.usc.edu/ GBshape/cgi-bin/hgLiftOver), and then to the Iberian lynx reference genome positions by using the available lynx-to-cat synteny (Abascal et al. 2016). In the end, this custom annotation resource enabled the classification of variants into intergenic, intronic, coding (CDS), and UCNE categories. The ability of SnpEff to predict the functional effects of variants allowed us to define three further mutation categories within CDS: synonymous, missense (i.e. non-synonymous), and loss-of-function (LoF). The latter group comprises variants whose predicted impact on function is classified by SnpEff as high effect, including stop-gained, start-gained, start-lost, splice-acceptor and splice-donor variants. The realisation that some of these categories constitute a wide amalgam of gene models and selection coefficients encouraged us to further break them down into subcategories. Missense tolerated vs. deleterious. We split missense variants based on their degree of conservation, the biochemical distance between the encoded aminoacids, and the sequence context, as inferred by Provean v1.1.5 (Choi and Chan 2015), and further detailed in the Supplementary Material. Mutations with a Provean score higher than -2.5 were classified as “missense tolerated”, whereas those with a score equal to or lower than -2.5 were classified as “missense deleterious”. UCNE moderate vs. UCNE high. Variants in regions annotated as UCNE were categorised into two groups according to their presumed evolutionary constraint, estimated through Genomic Evolutionary Rate Profiling (GERP) scores (Davydov et al. 2010) that were remapped from human to lynx (Supplementary Material). Mutations with 2 < GERP ≤ 5 were grouped into “UCNE moderate”, and those with GERP > 5 were classified as “UCNE high”. Variants with GERP ≤ 2 were discarded from further analyses on account of their low number and presumed low constraint. 149 A COMPARATIVE ANALYSIS OF GENOMIC LOAD GENOMIC VARIABLES BASED ON INDIVIDUAL COUNTS GATK v3.7 SelectVariants (McKenna et al. 2010) with the option –env (exclude non-variants) was used to subsample all three vcf files (containing all, only fixed, and only segregating variants, respectively) for each population and subsequently for each individual in the dataset. Then, the alternative allele count was recorded across sites and annotation categories to obtain the individual counts of derived alleles (i.e. the genomic load), heterozygous derived alleles, and homozygous derived alleles. Similarly, individual counts of fixed derived alleles were tallied to compare the fixation rates between the two species. For each of these four variables, we then calculated the population arithmetic mean of the individual counts, and the species average of the corresponding population means weighted by their corresponding sampling sizes. Meanwhile, the alternative allele count column from the population vcf files containing segregating variants was used to construct the allele count histograms. Only samples from the ILL and ELL datasets were considered for these analyses. We also obtained the number of derived alleles per individual separately for the autosomal subtelomeric and pericentromeric regions. These were defined at the contig level by including those with at least 1,000 bp annotated previously as subtelomeric or pericentromeric, by taking 10 Mb from the chromosome end or around the centromere, respectively, based on the synteny with cat (Abascal et al. 2016). Similarly, we counted the number of derived alleles per individual for the X chromosome, defined at the scaffold level based on the synteny with cat (Abascal et al. 2016), after discarding all bases assigned as subtelomeric in order to exclude the pseudoautosomal region. PER SITE ESTIMATES Average per site estimates for each category require an estimation of the total number of sites in the category following the filtering steps. Under the assumption that the effect of our filters on variant sites is proportional to that on monomorphic sites, we calculated the proportion of variants that passed the filters for the intergenic, intronic, coding, and UCNE categories, and then applied it to the total number of sites in the respective category. This number was calculated as the sum of the lengths of all corresponding tracks in the annotation file after subtracting repetitive and low mappability regions, as these had previously been excluded from the variant calling. In the case of synonymous, missense, and LoF mutations, which are potential states of coding variants rather than classes of sites, we used their relative densities as estimated by Yamaguchi-Kabata et al. (2008) for the human genome under a model of 150 CAPÍTULO 3 nucleotide change, according to which these mutations account for 28.5%, 68.1%, and 3.4% of the coding variants, respectively. Finally, to compute the per-site estimates of heterozygosity, the number of sites in the category was used as the denominator for the number of heterozygous derived alleles in the category. ERROR ESTIMATES AND SIGNIFICANCE TESTING In order to estimate the errors of the population and species averages for the different genomic variables, which allow testing whether the observed differences between them are significant, we need to consider two sources of variation. One is the sampling variance within each population or species, which produces a sampling error due to the limited number of individuals that have been sequenced. However, even if this sampling error was negligible (e.g. if a very high number of sequences from each population were available), populations that are evolving under identical selective pressures and demographic scenarios are expected to be different due to the stochastic nature of the evolutionary process (i.e. the random nature of breeding events, recombination, segregation, mutation, etc.), and thus, observed differences between the means of the populations or the two species may not be ascribed to differences in the forces governing evolution. In other words, evolutionary stochasticity must be considered as the second source of variation. For each variable analysed, the variance within (Vw) was estimated as the sample variance (for populations) or as the weighted average of the within population variances (for species). The sampling error of the average was then estimated as SE = s(Vw/T), where T is the total number of individuals available for that population or species. The stochastic evolutionary variance of each mean (Vs) was obtained by randomly resampling a set of genomic blocks 1,000 times to simulate additional sampling data (as detailed in the Supplementary Material), a procedure referred to as bootstrapping in Simons and Sella (2016). We then obtained the overall errors of the average of each population or species (M) as the square root of the corresponding overall variance, computed as s(M) = s(SE2 + Vs). Finally, at each of the following five categories: synonymous, missense tolerated, missense deleterious, LoF, and UCNE, we used this s(M) estimate for the average number of derived alleles per individual to assess the significance of the differences between species using a Z-test. We controlled the overall type-I error of the five categories using the sequential Bonferroni test, and we also identified the contrasts that are significant after controlling for the False Discovery Rate (FDR) (Benjamini and Hochberg 1995). 151 A COMPARATIVE ANALYSIS OF GENOMIC LOAD Error estimates and significance for tests focused on subtelomeric and pericentromeric regions were estimated analogously to those of the whole genome, although only blocks that had at least 85% of their sites annotated as subtelomeric or pericentromeric (respectively) were considered for the resampling simulation procedure (Supplementary Material). In addition, we used the Z-test to evaluate the significance of the differences between the two species for the variable Dtel-cen, which results from subtracting the pericentromeric from the subtelomeric species average for the number of derived alleles per individual. We plotted all aforementioned population and species averages, together with the estimated error bars, using the R (R Core Team 2020) package ggplot2 (Wickham 2016). OTHER POPULATION STATISTICS We computed and report as Supplementary Material other population statistics that have been used typically to address genetic load. These include the RX/Y statistic (Xue et al. 2015; Narasimhan et al. 2016), the πN/πS ratio, and the PN/PS ratio (McDonald and Kreitman 1991). 152 CAPÍTULO 3 RESULTS GENOMIC DIVERSITY We analysed 60 genomes (48 for most calculations) sampled from the two remnant populations of the highly endangered Iberian lynx and from three demographically distinct populations of the widespread Eurasian lynx in order to unravel the dynamics of their genomic load. The mean genome-wide heterozygosity was 2.34 x 10-4 (± 3.92 x 10-6) per base pair in the Iberian lynx (N = 20), which is roughly half that in the Eurasian lynx (4.24 x 10-4 ± 6.18 x 10-6; N = 28). The genetic diversity of the populations (Figure 1) ranked parallel to the available estimates of recent effective population sizes (Casas-Marce et al. 2013; Casas-Marce et al. 2017; Lucena-Perez et al. 2020): KIR, 4.90 x 10-4 (± 4.55 x 10-6); POL, 4.09 x 10-4 (± 1.02 x 10-5); NOR, 3.40 x 10-4 (± 1.78 x 10-5); AND, 2.84 x 10-4 (± 5.99 x 10-6); and DON, 1.60 x 10-4 (± 3.91 x 10-6). VARIANT ANNOTATION As explained in the introduction, bottlenecks and persistent small population sizes are predicted to shape the efficacy of selection, leaving an uneven footprint across regions under different pressures from natural selection. To analyse the effects of diverse demographic scenarios across populations, we considered several mutation categories that are expected to differ on their effects on fitness. Out of the total 4,388,391 high-confidence (i.e. post-filtering) biallelic SNPs that we identified in the pool of 60 samples, we annotated 2,838,376 (64.6%) in intergenic regions, 1,449,287 (33.0%) in introns, 54,736 (1.2%) in CDS sequences, and 893 (<<0.1%) in UCNEs. Within coding regions, SnpEff identified 31,023 (56.7%) synonymous; 23,452 (42.8%) missense; and 576 (0.01%) LoF mutations. Among missense variants, 14,378 mutations with low predicted impact (Provean score > -2.5) were subclassified as tolerated, and 9,264 with high predicted impact (Provean score ≤ -2.5) as deleterious. Similarly, UCNE variants were grouped according to sequence conservation-based GERP scores (Davydov et al. 2010) into two subcategories: moderate (moderately constrained; 2 < GERP ≤ 5), with 341 mutations, and high (highly constrained; GERP > 5), with 397 mutations. GENETIC DIVERSITY WITHIN ANNOTATION CATEGORIES Genetic diversity matched expectations of reduced heterozygosity in functional regions relative to neutral regions, particularly for the putatively more constrained 153 A COMPARATIVE ANALYSIS OF GENOMIC LOAD categories (missense deleterious, LoF and UCNE high). Although between population differences were smaller in less diverse categories, the population ranking for genome-wide nucleotide diversity held across all annotation categories (Figure 1; Supplementary Figure S1). Patterns of derived homozygosity were the inverse of diversity patterns in the case of neutral categories, but for the putatively deleterious ones, differences between populations disappeared or even changed sign (Supplementary Figure S2). POPULATION ALLELE COUNT HISTOGRAMS The frequency or count distribution of derived alleles provides a window into several aspects of the genetic state of a population. Focusing only on sites that are segregating within each species, we found a general trend of downwardly shifted frequency distributions for mutations with a higher predicted impact on function, such as those in the missense deleterious, LoF and UCNE categories, compared to intergenic and intronic mutations (Figure 2; Supplementary Table S2). In contrast, the histograms are roughly uniform across all categories in the highly bottlenecked DON population. Additionally, in NOR, POL, and DON we found a higher proportion of derived alleles that are locally fixed compared to the larger populations of each species, KIR and AND, particularly in the more neutral categories, as expected from classical theory (Kimura 1962). GENOMIC LOAD ESTIMATES The overall number of deleterious derived alleles (segregating and fixed) is the natural measurement of the genomic load. We counted the number of derived alleles per individual for each annotation category and represented the population and species averages. Using this statistic, we found similar average counts in both species for intergenic and introns, as expected for putatively neutral categories. In contrast, we found a deficit of LoF (~25%) and missense deleterious (~11%) –but not missense tolerated– mutations in Iberian lynx with respect to Eurasian lynx (Figure 3; Table 1), both of which remained significant even after applying Bonferroni and FDR adjustments (Benjamini and Hochberg 1995). This observation suggests that purging has played a major role in shaping the current genetic load of Iberian lynx. We also observed less UCNE high mutations and an unexpected slight deficit of synonymous derived alleles in Iberian with respect to Eurasian individuals, although neither reached statistical significance.