Population genomics as a tool for management and conservation of brown trout (Salmo trutta) in the Iberian Peninsula

DECLARACIÓN DO AUTOR/A DA TESE D./Dna. Adrián Casanova Chiclana Título da tese: Population genomics as a tool for management and conservation of brown trout (Salmo trutta) in the Iberian Peninsula Presento a miña tese, seguindo o procedemento axeitado ao Regulamento, e declaro que: 1) A tese abarca os resultados da elaboración do meu traballo. 2) De ser o caso, na tese faise referencia ás colaboracións que tivo este traballo. 3) Confirmo que a tese non incorre en ningún tipo de plaxio doutros autores nin de traballos presentados por min para a obtención doutros títulos. 4) A tese é a versión definitiva presentada para a súa defensa e coincide a versión impresa coa presentada en formato electrónico E comprométome a presentar o Compromiso Documental de Supervisión no caso de que o orixinal non estea na Escola. En Lugo, 07 de Xuño de 2021. Sinatura electrónica Este traballo foi financiado grazas as axudas de apoio á etapa predoutoral da Consellería de Cultura, Educación e Ordenación Universitaria e a Consellería de Economía, Emprego e Industria a través da Axencia Galega de Innovación (GAIN), cofinanciadas parcialmente polo programa operativo FSE Galicia 2014-2020, na convocatoria do ano 2017 (Nº expediente: ED481A-2017/091). ABSTRACT Brown trout (Salmo trutta) is a cold-water salmonid with ecological, commercial and recreational importance. Previous genetic studies highlighted a notable genetic structure in natural populations, which is a fundamental biological feature for the conservation of the species. Nevertheless, many factors threaten this genetic richness. One of the most important is genetic introgression, that originated from the introduction of aquaculture individuals in the wild environment in the last decades. This genetic erosion can have adverse consequences in brown trout adaptation capacity, especially in the context of a reduction of available habitats driven, among other causes, by global warming. This doctoral thesis had two main strands: (1) A bioinformatic benchmark to evaluate the biological conclusion robustness drawn from diverse SNP panels derived from different building-loci pipelines. For this strand, a broad of five aquatic species representing different genomic and/or population structure scenarios was used. This strand was performed in a context in which, due to the irruption of new techniques as RAD-seq (Restricion site-associated DNA sequencing) for the library preparation and Next Generation DNA Sequencing (NGS), the volume of data generated has grown exponentially in the last decade, promoting the development of new bioinformatics tools to process it. On the other side, (2) a genomic approach was used for the first time on brown trout populations from the Iberian Peninsula, to evaluate genetic diversity levels, population structure, natural hybridization patterns, evaluation of population introgression with aquaculture individuals of the same species and the detection of candidate genomic regions under selection pressure. The information obtained will allow the improvement of current management and conservation strategies of wild brown trout genetic resources. KEYWORDS Brown trout, bioinformatic benchmark, RAD-seq, Stacks 2, conservation genomics. RESUMO A troita común (Salmo trutta) é un salmónido de auga fría con importancia ecolóxica, comercial e recreativa. Estudos xenéticos anteriores puxeron de manifesto unha notable estrutura xenética nas poboacións naturais, o que constitúe unha característica biolóxica fundamental para a conservación da especie. Con todo, moitos factores ameazan esta riqueza xenética. Un dos máis importantes é a introgresión xenética, orixinada pola introdución de individuos de acuicultura no medio natural nas últimas décadas. Esta erosión xenética pode ter consecuencias adversas na capacidade de adaptación da troita común, especialmente no contexto dunha redución dos hábitats dispoñibles impulsada, entre outras causas, polo quecemento global. Esta tese de doutoramento presentou dúas liñas principais: (1) Unha referencia bioinformática para avaliar a robustez das conclusións biolóxicas derivadas da análise de diversos paneis de SNPs procedentes de diferente software bioinformático para a construción de loci. Para esta liña, empregáronse cinco especies acuáticas que representaban diferentes escenarios xenómicos e/ou de estrutura poboacional. Esta investigación realizouse nun contexto no que, debido á irrupción de novas técnicas como RAD-seq (Restricion siteassociated DNA sequencing) para a preparación de librarías e Next Generation DNA Sequencing (NGS), o volume de datos xerados creceu exponencialmente na última década, promovendo o desenvolvemento de novas ferramentas bioinformáticas para o seu procesamento. Doutra banda, (2) utilizouse por primeira vez unha aproximación xenómica en poboacións de troita común da Península Ibérica, para avaliar os niveis de diversidade xenética, estrutura poboacional, patróns de hibridación natural, avaliación da introgresión poboacional con individuos de acuicultura da mesma especie e a detección de rexións xenómicas candidatas a estar baixo presión selectiva. A información obtida permitirá mellorar as actuais estratexias de xestión e conservación dos recursos xenéticos da poboacións naturais de troita común. PALABRAS CHAVE Troita común, referencia bioinformática, RAD-seq, Stacks 2, xenómica da conservación. RESUMEN La trucha común (Salmo trutta) es un salmónido de agua fría con importancia ecológica, comercial y recreativa. Estudios genéticos anteriores pusieron de manifiesto una notable estructura genética en las poblaciones naturales, lo que constituye una característica biológica fundamental para la conservación de la especie. Sin embargo, muchos factores amenazan esta riqueza genética. Uno de los más importantes es la introgresión genética, originada por la introducción de individuos de acuicultura en el medio natural en las últimas décadas. Esta erosión genética puede tener consecuencias adversas en la capacidad de adaptación de la trucha común, especialmente en el contexto de una reducción de los hábitats disponibles impulsada, entre otras causas, por el calentamiento global. Esta tesis doctoral presentó dos líneas principales: (1) Una referencia bioinformática para evaluar la robustez de las conclusiones biológicas derivadas del análisis de diversos paneles de SNPs procedentes de diferente software bioinformático para la construcción de loci. Para esta línea, se emplearon cinco especies acuáticas que representaban diferentes escenarios genómicos y/o de estructura poblacional. Esta investigación se realizó en un contexto en el que, debido a la irrupción de nuevas técnicas como RAD-seq (Restricion site-associated DNA sequencing) para la preparación de librerías y Next Generation DNA Sequencing (NGS), el volumen de datos generados ha crecido exponencialmente en la última década, promoviendo el desarrollo de nuevas herramientas bioinformáticas para su procesamiento. Por otro lado, (2) se ha utilizado por primera vez una aproximación genómica en poblaciones de trucha común de la Península Ibérica, para evaluar los niveles de diversidad genética, estructura poblacional, patrones de hibridación natural, evaluación de la introgresión poblacional con individuos de acuicultura de la misma especie y la detección de regiones genómicas candidatas bajo presión selectiva. La información obtenida permitirá mejorar las actuales estrategias de gestión y conservación de los recursos genéticos de las poblaciones naturales de trucha común. PALABRAS CLAVE Trucha común, referencia bioinformática, RAD-seq, Stacks 2, genómica de la conservación. RESUMO AMPLO 1 INTRODUCCIÓN Troita común (Salmo trutta) A troita común (Salmo trutta) é unha especie de peixe de auga doce, pertencente ao orde Salmoniformes. Esta orde está formada por unha única familia, Salmonidae dentro da cal existen dúas especies autóctonas en España, pertencentes ao xénero Salmo: a troita común (Salmo trutta) e o salmón atlántico (Salmo salar) (Doadrio 2001). Varias especies do xénero Salmo son filoxenéticamente próximas polo que adoita considerarse a S. trutta coma un complexo multiespecie (Gratton et al. 2014; Splendiani et al. 2019; Vera et al. 2011). Esta especie vive en regatos rápidos, ríos e lagos fríos e ben osixenados (7 mg O2/L e unha saturación de osíxeno do 80%; Cobo et al. 2015; Eklöv et al. 1999). O seu rango de distribución natural abrangue principalmente Eurasia (Freyhof 2011), desde o norte de Noruega e o nordeste de Rusia oriental (Bernatchez et al. 1992) ata o Cordal do Atlas no norte de África (Sanz 2018). Debido ás introducións mediadas polo ser humano, a troita común ten un área de distribución foránea que inclúe o extremo oriental de Asia, Oceanía, África e América (Casal 2006). A troita común presenta unha ampla variedade fenotípica relacionada coa diversidade do medio no cal vive e o seu comportamento migratorio. Seguindo criterios ecolóxicos e de ciclo vital identificáronse tradicionalmente tres formas de troita común: (1) residente, (2) anádroma (reo) e (3) lacustre. Os reos pódense atopar na vertente atlántica, non así na mediterránea, sendo o río Limia (42o N) o seu límite meridional de ocorrencia (Bouza et al. 1999). Estes exemplares aliméntanse principalmente preto da costa, non moi lonxe da desembocadura dos ríos de procedencia (Kottelat e Freyhof 2007). A troita é un peixe de tamaño medio que en España pode acadar os 60 cm e os 10 kg, acadando os individuos anádromos e lacustres maiores tamaños (45-60 cm de lonxitude estándar, SL), que os residentes (20-30 cm SL, Kottelat y Freyhof 2007). A troita común pode presentar manchas escuras no dorso e nos laterais do corpo, as veces rodeadas por un halo pálido. Os reos tralo esguinado aforrar recursos de soporte físico e tempo de computación. A mellor configuración de parámetros sería aquela que conducise aos máis consistentes a traveso das diferentes réplicas e empregando os parámetros poboacionais que correspondan aos obxectivos de investigación. A pesares de implicar un consumo de tempo, este paso preliminar melloraría a solidez dos resultados e as conclusións biolóxicas derivadas ao mellorar o uso das ferramentas bioinformáticas. 2. Análise da introgresión xenómica das poboacións silvestres por parte das poboacións de piscifactoría nos ecosistemas inestables da vertente mediterránea. 3. Estudo preliminar da hibridación natural nun contacto secundario na vertente atlántica. 4. Mellora das estratexias de conservación e xestión para manter a integridade dos liñaxes autóctonos e recuperar as poboacións naturais en zonas con alta taxa de introgresión. Para acadar estes tres obxectivos diferentes, mostras de troita común procedentes da Península Ibérica foron tomadas tanto da vertente atlántica (concas do ríos Miño e Duero) coma da vertente mediterránea (rexión Pirenaica), ademais de troitas procedentes dun centro piscícola (Bagà). Nesta mostraxe atopábanse representadas as diferentes liñaxes mitocondriais da Península Ibérica (AD, AT, DU, ME), poboacións naturais afectadas con diferente intensidade por soltas desde piscifactorías e zonas de hibridación natural. Empregouse Bowtie 1 (Langmead 2009) e Stacks 2 xunto co xenoma de referencia de troita común para a obtención de loci e a detección e xenotipado dos SNPs. Trala obtención dun panel robusto con miles de SNPs realizouse unha serie de análises bioinformáticas para determinar o grao de afección por soltas de piscifactoría desde un enfoque xenómico, introgresión ao longo do xenoma, diversidade xenética, estrutura poboacional, estimación do censo efectivo, deteccións de pegadas de selección e anotación funcional. Deseñouse unha ferramenta molecular para avaliar os efectos da repoboación nas poboacións naturais de troita común, permitindo detectar a introgresión a nivel individual e poboacional dun xeito máis barato e cunha maior resolución, grazas ao elevado incremento de marcadores moleculares empregados, respecto á metodoloxía empregada ata o momento coas poboacións da Península Ibérica (emprego do locus diagnóstico LDH-C*). Coa resolución xenómica acadada puidéronse detectar tamén pequenas rexións de introgresión ao longo do xenoma, inferindo a ancestralidade local a nivel cromosómico. Todas as poboacións afectadas por repoboación estiveron en equilibrio Hardy-Weinberg, o que suxire que as variantes alélicas dos stocks de piscifactoría introducidas durante as últimas décadas integráronse e aleatorizáronse no patrimonio xenómico das poboacións silvestres. Detectáronse grandes diferencias nos valores de diversidade xenética entre as diferentes poboacións de troita común na Península Ibérica. Confirmouse o patrón de aumento da diversidade xenética seguindo un gradiente leste-oeste e sur-norte, cos valores máis altos na conca hidrográfica do Miño. Este patrón estaría en consonancia coas rexións xeográficas cunhas condicións ambientais máis adecuadas para esta especie. Detectouse unha fonda estruturación entre as poboacións pertencentes a diferentes concas, incluso entre poboacións próximas. A maior diferenciación xenética detectouse entre as vertentes atlántica e mediterránea, consecuencia da ausencia de fluxo xénico debida ao illamento xeográfico existente. Dentro da vertente atlántica, as concas do Duero e do Miño apenas terían fluxo xénico debido á escaseza de reos que poidan comunicar ambas concas, de habelos. Empregando diferentes subconxuntos de SNPs nucleares, non se detectaron as zonas de hibridación natural previamente reportadas en estudos anteriores. Isto suxeriría que só quedan rastros dos eventos de hibridación por contacto secundario debido á forte deriva xenética asociada ao pequeno censo efectivo nas poboacións de troita común. Detectáronse SNPs que poden estar baixo presión selectiva (i.e. outliers) empregando diferentes niveis xerárquicos, na súa maioría de selección estabilizadora entre as vertentes mediterránea e atlántica, namentres que baixo selección diverxente na vertente atlántica tanto dentro como entre concas (i.e. Miño-Duero). De todos xeitos, a detección bioinformática destes outliers sería o paso previo a análises máis exhaustivas para a validación de presión selectiva nesas rexións xenómicas, mediante análises de expresión xénica en experimentos de ambiente común (common-garden). A modo de peche, a información obtida da análise xenómica poboacional realizada nesta tese será de gran valor para deseñar as directrices de conservación das poboacións de troita común da Península Ibérica. Por exemplo, poderíanse aplicar diferentes políticas de conservación en función do grao de ancestralidade cos stocks de piscifactoría nas poboacións naturais de troita común: (1) erradicación das poboacións naturalizadas das soltas de piscifactoría que aínda poidan existir; e (2) establecemento de diferentes cotas de pesca (individuos por persoa e día) en función dos diferentes graos de introgresión detectados, reducindo as cotas nos ríos con poboacións naturais non afectadas xenéticamente polas soltas de piscifactoría. A detección de zonas de hibridación natural, non inducidas directamente polas actividades humanas, sería interesante desde o punto de vista da conservación, para protexer os procesos naturais que poidan aumentar a diversidade xenética e a viabilidade das poboacións. Algunhas das poboacións analizadas mostraron un baixo censo efectivo, o que podería levar á definición das poboacións de troita común como vulnerables segundo este criterio e posteriormente ao establecemento de medidas de conservación como moratorias temporais de pesca ou redución de cotas en unidades hidrográficas ben establecidas (e.g. arroios, tramos fluviais). A alta diferenciación xenética detectada nas poboacións de troita común pertencentes a unha mesma conca hidrográfica, implicaría que a repoboación baixo principios de conservación (utilizar individuos xenéticamente similares ás poboacións repoboadas) contemplada nos plans e leis de conservación pode ser difícil de implementar dunha forma estrita. En caso de aplicarse, deberían utilizarse poboacións non afectadas por eventos de introgresión, segundo os enfoques de ascendencia global e local. A definición de refuxios de troita común, nos que debería prohibirse a repoboación ou a pesca, en tramos de diferentes ríos con boas condicións ecolóxicas podería ser unha acción complementaria para a conservación. De todos os xeitos, para adoptar medidas de xestión e conservación, sería necesario traballar con información actualizada sobre unha serie de poboacións representativas de diferentes unidades da Península Ibérica (e.g. ríos, concas hidrográficas, Unidades de Xestión). O ideal sería que as estratexias fosen consensuadas e compartidas entre as diferentes Administracións para un plan de conservación integrado máis aló dos límites administrativos. A información xenómica obtida, combinada cos resultados doutras disciplinas (e.g. ecoloxía), debería axudar a tomar as mellores decisións de xestión. En calquera caso, ningunha medida de conservación debería abordarse sen os contextos ecolóxicos (e.g. o estado de conservación dos hábitats) e socioeconómicos asociados (e.g. a lexislación vixente e os usos económicos e culturais existentes). Relación de artigos que aportan contido á presente Tese de doutoramento. Os contidos corresponden ao capítulo 3 e capítulo 5.1 da presente Tese de doutoramento. Título: Low impact of different SNP panels from two building-loci pipelines on RAD-Seq population genomic metrics: case study on five diverse aquatic species Ano: 2021 Revista: BMC Genomics Volumen, Páx.: 22, Article number 150 DOI: https://doi.org/10.1186/s12864-021-07465-w Contribución do doutorando: Adrián Casanova Chiclana codeseñou o estudo, realizou a meirande parte das análises bioinformáticas e escribiu a versión inicial do manuscrito ademais das diferentes versións coa contribución do resto de coautores Índices/s de calidade: 4.093 5-year Citation Impact Factor Este é un artigo de acceso aberto distribuído baixo os termos da licencia Creative Commons CC BY, que permite o uso, distribución e reprodución sen restricións en calquera medio, sempre que o traballo orixinal estea debidamente citado. Non é necesario a obtención de permiso para reempregar este artigo. Aplícase CC0 para o material complementario relacionado con este artigo e non se require atribución. AGRADECEMENTOS En primeiro lugar gustaríame dar as grazas aos meus directores de tese Manuel Vera e Paulino Martínez. Por toda a súa axuda, dedicación e formación nesta etapa. Todas as palabras de agradecemento para vós serían escasas. Para todas as persoas coas que tiven o privilexio de formar equipo neste período no Departamento de Xenética. Profesionais das que puiden aprender e persoas coas que puiden compartir. A Hurzos, por aquela chamada. A todas as miñas amizades, a miña benquerida familia, A Asociación Cultural Náufragos do Paradiso e a todas as súas integrantes; por compartir a vida. A Tamara, por ser unha estrela fugaz feita realidade. E por suposto, e con todo o agarimo que pode conter un fillo. Mil grazas a miña mai, María. Por todo o bo que teño. This page is intentionally left blank. No tracks in the road, nothing living anywhere. The fire blackened boulders like the shapes of bears on the starkly wooded slopes. He stood on a stone bridge where the waters slurried into a pool and turned slowly in a grey foam. Where once he had watched trout swaying in the current, tracking their perfect shadows on the stones beneath. The Road (2006). Cormac McCarthy. A vida resúmese en saltar presas. Troita común. INDEX 1 General introduction ................................................................... 1 1.1 Brown trout............................................................................. 1 1.1.1 Taxonomy and distribution ............................................. 1 1.1.2 Morphology ..................................................................... 2 1.1.3 Anadromy of Iberian brown trout ................................... 4 1.1.4 Biology and ecology ....................................................... 6 1.1.5 Threats to brown trout ..................................................... 8 1.1.6 Human perspective ........................................................ 10 1.2 Brown trout genomics .......................................................... 11 1.2.1 Whole genome duplication and diploidization of salmonids ..................................................................................... 11 1.2.2 Genetic structure and phylogeography ......................... 13 1.2.3 Hybridization and introgression scenarios .................... 15 1.3 Conservation genetics/genomics .......................................... 17 1.4 Genetic and genomic tools ................................................... 21 1.4.1 DNA sequencing technologies ...................................... 21 1.4.2 Molecular markers ........................................................ 24 1.5 Brown trout population genomics ........................................ 25 2 Objectives ................................................................................... 27 3 Low impact of different SNP panels from two building-loci pipelines on RAD-seq population genomic metrics: case study on five diverse aquatic species .............................................................. 28 3.1 Introduction .......................................................................... 28 3.2 Material and Methods .......................................................... 35 3.2.1 Sampling ....................................................................... 35 3.2.2 Library preparation ....................................................... 36 3.2.3 Bioinformatic analysis .................................................. 36 3.3 Results .................................................................................. 43 4 Population genomics for conservation of brown trout resources in the Mediterranean and Atlantic slopes of Iberian Peninsula: structure, selection and patterns of hybridization and introgression in Atlantic and Mediterranean river drainages ............................. 53 4.1 Introduction .......................................................................... 53 4.2 Material and Methods .......................................................... 61 4.2.1 Sampling ....................................................................... 61 4.2.2 Identification and genotyping of SNPs: defining a reference panel for all populations studied ................................. 65 4.2.3 Identification of farmed and wild samples ................... 66 4.2.4 Genetic diversity ........................................................... 69 4.2.5 Genetic differentiation and structure ............................ 70 4.2.6 Effective population size estimation ............................ 72 4.2.7 Adaptative variation ..................................................... 73 4.3 Results .................................................................................. 75 4.3.1 Identification of individuals of farmed ancestry........... 78 4.3.2 Genetic diversity ........................................................... 87 CASANOVA CHICLANA, ADRIÁN 6 1.1.4 Biology and ecology Brown trout is a sexual dimorphic species with differences in secondary sexual characters (Reyes-Gavilán et al. 1997), but even in habitat use (Klemetsen et al. 2003). Brown trout normally reach sexual maturity between two and three years of life, later in females than in males (Alp et al. 2003). The reproduction is in autumn or winter and is phenotypically detectable by an increase in the gonadosomatic index (GSI; Jamalzadeh et al. 2013). The reproduction period is sooner at higher altitudes or latitudes due to the lower temperature of the water involving longer incubation periods (Cobo et al. 2015). The development time is usually measured in degree-days that accumulates the number of oC through days of development (Ojanguren and Braña 2003) and that would be around 450 degree-days on average (Vøllestad and Lillehammer 2000). The incubation of eggs takes more than 40 days with normal freshwater temperatures in Iberian Peninsula (Doadrio 2001). Brown trout is highly selective in their choice of the spawning area and the gravel granulometry is an important factor (Barlaup et al. 2008). Females select the most appropriate spawning sites and dig the redd where the eggs are deposited (Kottelat and Freyhof 2007), generally in well-oxygenated waters because hatching of eggs requires a constant flow of water through the gravels to oxygenate them and remove waste products. Fertilization is external and occurs immediately after the eggs and sperm are expelled into the aquatic environment. In the Northern Hemisphere, the main period of growth is from April to November, when water temperatures are mostly between 5 and 20 oC. In this range, salmonids feed enough, not only to maintain themselves but also for fattening and maturation (Cobo et al. 2015). In fish, standard metabolic rate (SMR) is primarily limited by the capacity of the gills to extract oxygen from water. An increasing in oxygen uptake is necessary to do activities such as feeding, swimming, General introduction 7 etc (Eddy and Handy 2012). Like other fish, brown trout present skin breathing whose extent depends on environmental conditions (e.g. water temperature), but usually lower than 10% (Hill et al. 2016). Brown trout is susceptible to human activities that may affect water quality (e.g. pollution) and requires well oxygenated freshwater with higher concentrations than 7 mg O2/L and oxygen saturation around 80% (Cobo et al. 2015; Eklöv et al. 1999). Lower concentrations can lead to delayed fry development and growth (Dumas et al. 2007) up to massive fish deaths. This species was evaluated as possible biomonitor in different studies (Lamas et al. 2007; Linde et al. 1998; Schmidt et al. 1999). Brown trout has been considered a flagship (Denic and Geist 2010) or umbrella species (Lobón-Cerviá 2018). Although these terms have been used interchangeably in different studies (Caro 2010), they have different meanings. A flagship species would be an emblematic species selected to act as icon or symbol for a defined habitat, environmental campaign, etc; while an umbrella species could be species whose conservation confers protection to many naturally co-occurring species (Roberge and Angelstam 2004). The use of these terms to refer to the brown trout shows its importance for the functioning of the ecosystem and the value for humans. One of the main ecological relationship would be predator-prey relationship. Brown trout is a top fish predator wherever it lives (Sánchez-Hernández 2016). Furthermore, it is an organism with a very broad diet (i.e euriphagic species; Alonso et al. 2017). The composition of brown trout diet can vary among individuals, populations, age classes, habitats, and seasons. Habitually, the size of the prey increases with the size of the individual due to the increase of fish length and mouth gape. This implies that dietary changes during the ontogenetic development of top fish predators involve changes in structural properties of food webs (Sánchez-Hernández 2016). CASANOVA CHICLANA, ADRIÁN 8 Beyond its trophic relationship as an upper-level predator in ecosystem it should be noted the condition of brown trout as a host (i.e parasitism relationship) with the freshwater pearl mussel (Margaritifera margaritifera), which is an endangered species (Moorkens et al. 2017). In their life cycle, freshwater mussels have a specialized larva (i.e. glochidium), which must usually parasitize a host fish upon which it encysts and metamorphoses. Salmo trutta and S. salar are the specific host fish on which the M. margaritifera glochidia develop. This species is an integral component of the river ecosystem, filtering particles, serving as a substratum or habitat for algae and benthic invertebrates, respectively. Therefore, these ecological functions performed by freshwater mussels can be compromised in regions with lower mussel specific richness as the Iberian Peninsula (see Fig. 7 in Lopes-Lima et al. 2017), highlighting the conservation interest of its host species. Brown trout, outside its native distribution range is involved in different drawbacks with native fish species such as (1) reducing population size by different ways (e.g. food competition, predation), (2) hybridization (e.g. with Salvelinus fontinalis in United States of America, USA; see Blanc and Chevassus 1986; Sorensen et al. 1995). Being this species considered by the Invasive Species Specialist Group (ISSG) as one of the 100 World's Worst Invasive Alien Species (http://www.issg.org/worst100_species.html). 1.1.5 Threats to brown trout Brown trout is as Least Concern in the IUCN Red List category (Freyhof 2011); however, it is affected by habitat fragmentation, water pollution, global warming, overfishing and hatchery introgression. One characteristic in any ecosystem is the connectivity. Referred to river ecosystems there are many definitions (see Table 1 in Wohl 2017). One of them would be the circulation of organisms, matter, and General introduction 9 energy through the riverine landscape via the aqueous medium (Amoros and Roux 1988). As a multidimensional space, rivers have different connectivity components (i.e. vertical, lateral, and longitudinal), nevertheless for salmonids the most studied was longitudinal connectivity. The most prominent constructions in rivers are dams, built for different human purposes (e.g. hydroelectric power, water level regulation and irrigation). Nevertheless, dams affect fish migration, isolating populations, and can affect water properties and habitat conditions, upstream and downstream (e.g. mean summer temperature in downstream water; see Lessard and Hayes 2003). Heggenes and Røed (2006) detected a significant correlation between years since dam construction and FST values, with an increasing local differentiation presumably caused by genetic drift. Birnie-Gauvin et al. (2018) documented the high increase of smolt number in Villestrup River (Denmark) after the removal of six dams with fishways. To improve river connectivity, nature-like bypass channels have been built in different rivers with good results, although they require an important investment (Dodd et al. 2017). For study and compare river connectivity, the Catalonian Water Agency developed in 2006 the Index de Connectivitat Fluvial (ICF; Solà et al. 2011). Water pollution is the introduction of contaminants, related with human-mediated activities, in water bodies. Brown trout is considered a sensitive species to water pollution (see Luckenbach et al. 2001; Pickering 1989), although there have been reported some cases of tolerance in highly polluted ecosystems (e.g. metal pollution; Durrant et al. 2011), related to molecular mechanism of tolerance (see Webster et al. 2013). There is a wide variety of contaminants studied with the brown trout, from pesticides to microplastics (see Nusbaumer et al. 2020; O’Connor et al. 2020) and the bioaccumulation in different tissues due to the importance for human consumption (see Linde et al. 2004). CASANOVA CHICLANA, ADRIÁN 10 The ecological impacts of human societies have increased over the centuries, especially since the industrial revolution. These impacts are carried out in all scales (i.e. from local to global scale). A new geological epoch-term has been coined to reflect the magnitude of these impacts: “The Anthropocene” (Crutzen and Stoermer 2000). Global warming is one of the phenomena best known by society, with a range of consequences for the climate, economy, biodiversity, etc. The increase of temperatures would lead to the loss of suitable thermal brown trout habitats in lower latitudes (Almodóvar et al. 2012) or mountain regions (Hari et al. 2006). Although, adaption can play a role to changing temperature regimes (Jensen et al. 2008). Restocking can endanger native populations by disrupting local adaptations with the introduction of domestication selection (Pinter et al. 2019) and influencing over different traits as dispersal (Saint-Pé et al. 2018). The current trends are conservation aquaculture approaches (e.g. enhance wild populations with individuals of the highest genetic similarity; see Froehlich et al. 2017) or the release of triploid individuals for fishing (i.e. sterile individuals). Nevertheless, triploidy is not always 100% successful and effective methods to detect diploid individuals would be necessary to avoid the release of thousands of fertile hatchery individuals (Sanz et al. 2020). 1.1.6 Human perspective Recreational fishing is the main use of wild fish stocks in freshwaters belonging to developed countries (Arlinghaus et al. 2002), although it also exists in developing countries. According to the FAO, recreational fishing would be “fishing of aquatic animals (mainly fish) that do not constitute the individual’s primary resource to meet basic nutritional needs and are not generally sold or otherwise traded on export, domestic or black markets" (FAO 2012). When this activity implies a displacement outside the region of origin of anglers, it could be spoken about angling tourism, linked with active tourism, General introduction 11 ecotourism, and village tourism (FAO 2012). There are some studies about the economic impact of recreational fishing in some regions (e.g. Scotland, Butler et al. 2009; Gotland, Blicharska and Rönnbäck 2018), that can reach estimated incomes of the order of hundreds millions euros (see Inland Fisheries Ireland 2013). The number of freshwater fishing licences issued in Spain in 2018 was around 490,000 with an economic value close to seven million euros (Ministerio para la Transición Ecológica y el Reto Demográfico, MITERD 2020). To face up the angler’s demand, every year releases of individuals or eggs of different freshwater species are performed. Among them, the brown trout represents the highest number of individuals or eggs released, exceeding 60% of the total (MITERD 2020). These data are not strange if we bear in mind that brown trout is one of the most valued species by anglers (FAO 2012), being one of the fish species providing an ecosystem service of utilitarian value (Millennium Ecosystem Assessment 2005). 1.2 BROWN TROUT GENOMICS 1.2.1 Whole genome duplication and diploidization of salmonids Brown trout is a species of tetraploid origin in process of diploidization (Leitwein et al. 2017). Polyploidy can represent advantages or disadvantages for organisms. For instance, in sexual organism advantages would be related to heterosis and gene redundancy. The former could render more vigorous fish than their diploid counterparts, while the last would shield polyploids from the deleterious effect of mutations. On the other hand, polyploidy can originate problems in meiosis due to homologous mispairing and bad segregation, thus giving rise to gamete unbalance and lower viability (Comai 2005). CASANOVA CHICLANA, ADRIÁN 12 There have been four whole genome duplication (WGD) events from the common ancestor of vertebrates to salmonid divergence respect to the rest of teleost fish. Two whole genome duplications (referred to as the 1R and 2R WGDs) occurred in the common ancestor of all vertebrates, determining an increase in chromosome number from 10-13 in the protovertebrate karyotype to 40-52 in most vertebrates (Nakatani et al. 2011; Ohno 1970). A third (3R) WGD occurred 400 million years ago (Ma) in the ancestor of teleost fish, while a fourth salmonid-specific WGD (Ss4R) took place ~90 Ma in the ancestor of salmonid fish (Allendorf and Thorgaard 1984; Macqueen and Johnston 2014). According to Glover et al. (2016), chromosomes or genes originated by whole-genome duplication events within species are called ohnologs, whilst those originated by allopolyploidization of a hybrid between different species would be homoeologs. Nowadays, the Ss4R is considered an autotetraploid genome duplication event (Campbell et al. 2020; Christensen and Davidson 2017). Salmonids would be in a diploidization process after the duplication event (Pouzadoux et al. 2017), which would imply disomic inheritance. This diploidization would not have been completed since multivalent pairing is observed at salmonid meiosis (Allendorf and Thorgaard 1984) supporting the tetrasomic inheritance observed (Allendorf et al. 2015), specially toward the telomeres (Nugent et al. 2017). The Atlantic salmon is the closest relative of brown trout and the number of chromosomes of their karyotypes largely differ (S. salar 2n = 58 vs S. trutta 2n = 80). There have been reported signals of intense chromosomal rearrangements through the evolution of these species, consistent in 13 fusions and two fissions which occurred in the Atlantic salmon branch after S.salar/S.trutta divergence (Pouzadoux et al. 2017). General introduction 13 Salmonids are gonochoristic (separated sexes), and sex would be genetically determined following a XX/XY system (Davidson et al. 2009; Yano et al. 2013). Nevertheless, sex determination can be influenced by environmental factors (e.g. water temperature), including the presence of pollutants or hormone treatments (Johnstone et al. 1978). X and Y chromosomes in salmonids display a large pseudoautosomal region and a small sex determining region (Davidson et al. 2009), so chromosomes are essentially homomorphic when inspected with optical microscope. In the last years, a novel master sex-determining gene has been characterized in rainbow trout (Yano et al. 2012), named sdY (sexually dimorphic on the Y chromosome). The presence of sdY only in phenotypic males of most salmonid species strongly suggests its conservation in the group (Araneda et al. 2019; Yano et al. 2013). 1.2.2 Genetic structure and phylogeography Brown trout is a species with high genetic differentiation among populations, even those belonging to nearby sections of the same river (Bouza et al. 1999; Fernández-Cebrián et al. 2014), being one of the most structured vertebrate species (FST > 0.60; Ferguson 1989). However, this general rule is not always fulfilled (Heggenes and Røed 2006). In the Iberian Peninsula, this high differentiation was observed between river basins of the Mediterranean and Atlantic slopes using different molecular markers, such as allozymes (Bouza et al. 1999, 2001), Restriction Fragment Length Polymorphisms (RFLPs; Machordom et al. 2000) and microsatellites (Campos et al. 2006; Fernández-Cebrián et al. 2014; Martínez et al. 2007; Vera et al. 2010). Brown trout displays a complex phylogeography, tangled by human-induced processes of introgression between lineages (Sanz 2018). Five major mitochondrial lineages were initially identified in brown trout using the control region of mitochondrial DNA (CR mtDNA; Bernatchez et al. 1992): Adriatic (AD), Atlantic (AT), CASANOVA CHICLANA, ADRIÁN 14 Danubian (DA), Marmoratus (MA) and Mediterranean (ME). Since then, several studies reported additional lineages restricted to the Iberian Peninsula, the Duero lineage (DU; Vera et al. 2010); the Tigris lineage in Turkey (TI; Sušnik et al. 2005); and an eighth North African lineage (Tougard et al. 2018), previously considered a different species (i.e. Salmo macrostigma). The mtDNA control region has been the main molecular marker used to identify geographically structured phylogenetic assemblages (Fig. 5) and mtDNA haplogroups have been used as reliable proxies of evolutionary significant units (ESUs) combined with nuclear molecular markers such as restriction fragment length polymorphisms (RFLPs; Castro et al. 1999; Schöffmann et al. 2007), internal transcribed spacers (ITS) of rDNA (Presa et al. 2002), microsatellites (Vera et al. 2013) or single nucleotide polymorphisms (SNPs; Pustovrh et al. 2011). Figure 5. Geographical distribution of the S. trutta species complex lineages based on mitochondrial DNA data. Four lineages are represented in the Iberian Peninsula: Atlantic (black squares), Duero (Iberian Peninsula endemism; white square), Adriatic (black circle) and Mediterranean (white circle). Figure from Sanz (2018). General introduction 15 1.2.3 Hybridization and introgression scenarios Within fish, hybridization is facilitated by reproductive features such as external fertilization and common spawning grounds (Hubbs 1955). When talking about hybridization and introgression in brown trout three different scenarios should be highlighted: (1) genetic introgression of wild populations with hatchery stocks, (2) natural hybridization between evolutionary lineages (secondary contacts), and (3) hybridization with closely related species (e.g. S. salar). Releasing farmed brown trout into the wild has been a management practice applied since the early 1900's. Stocking in European river basins have been mainly performed with a stock of Central European origin (Martínez et al. 1993; Morán et al. 1991; Vera et al. 2018). This practice has been performed with other salmonid species as well (e.g. Salvelinus fontinalis, Lehnert et al. 2020; S. salar, Saltveit 2006). The second aforementioned scenario is related to natural hybridization between lineages and/or subspecies, even though the validity of the Salmo trutta complex taxonomy is controversial. Experiments performed to test reproductive isolation between S.trutta subspecies belonging to the same mitochondrial lineage finally demonstrated a single biological entity (four subspecies belonging to Danubian lineage; Kalayci et al. 2018). In the Iberian Peninsula two main natural hybridization scenarios have been reported: (1) between AT and DU lineages in Duero River Basin and Miño-Sil River Basin with a parapatric distribution pattern (Martínez et al. 2007; Vera et al. 2010; Vilas et al. 2010); and (2) between AD and MED with a mosaic distribution pattern (Sanz et al. 2002). The final scenario outlined before would include: (1) natural interspecific hybridization, and (2) human-mediated interspecific hybridization. In the first one, probably the most relevant case would be the hybridization between S. trutta and S. salar in the north of Spain CASANOVA CHICLANA, ADRIÁN 22 Dideoxi-nucleotides are labelled with different fluorescent labels, activated by the laser into automated sequencers (Freeland 2011). Second-generation sequencing (better known as Next Generation Sequencing, NGS) is a high-throughput DNA sequencing protocol where billions of DNA strands are sequenced in parallel after an amplification step of each individual DNA molecule. This technology has reduced prices drastically and generated an enormous amount of data (reads: sequences of nucleotides from the sequenced DNA molecules). This volume of information requires new bioinformatic tools and more hardware requirements (e.g. multicore processors, higher volume of memory), that can be achieved easily in supercomputing centres (e.g. BSC, Barcelona Supercomputing Center; CESGA, Centro de Supercomputación de Galicia). Third generation sequencing (TGS) does not require amplification and single DNA molecules are sequenced in real time (SMRT; Wong et al. 2019). Recently, fourth-generation sequencing has been described as the compendium of techniques to conduct genomic analysis directly in the cell and this is expected to be useful in more specific applications than other sequencing generations with broader applicability (Ari and Arikan 2016). Since the beginning of NGS technology, different platforms have emerged with specific sequencing features such as read length, number of DNA molecules sequenced in parallel, etc (von Bubnoff 2008). Related to read length, there are two broad categories: Short-read NGS (usually < 300 pb) and long-read NGS (> 10 kb on average) (Mantere et al. 2019). Short-read NGS can be classified in two general types: sequencing by synthesis (SBS) and sequencing by ligation (SBL) (Slatko et al. 2018). Briefly, SBS uses DNA polymerase to incorporate complementary nucleotides to the elongating strand (Illumina platform) and SBL uses the mismatch sensitivity of DNA ligase to determine the underlying sequence of nucleotides in each DNA sequence (SOLiD General introduction 23 platform; Ho et al. 2011; Huang et al. 2012). According to Goodwin et al. (2016) the sequencing by synthesis can be classified as: cyclic reversible termination (CRT; Illumina, Qiagen platforms) and singlenucleotide addition (SNA; 454 pyrosequencing devices, Ion Torrent). In this thesis a sequencing by synthesis approach was used, provided by Illumina technology, which will be explained in more detail. For more information, it is encouraged the reading of these reviews: Goodwin et al. (2016) and Wong et al. (2019). Illumina sequencing workflow can be summarized in three steps: (1) DNA library preparation by fragmentation of DNA followed by 5’ and 3′ adapter ligation with barcodes to identify individuals that are sequenced together (multiplexed); (2) amplification of each fragment using bridge-PCR in a flow cell, which generates clonal clusters from bound fragments; and (3) read sequencing with fluorophore-labelled nucleotides, terminally blocked with a 3’-O-azydomethyl group, which hybridize to complementary bases linked to the emission of a fluorescent colour. Finally, fluorophores are washed, and nucleotides 3’-OH group are regenerated, by this way a new cycle begin with the addition of new nucleotides. This cycle is repeated “n” times to create a read length with “n” nucleotides. Sequencing data is exported to an output file, usually FASTQ format with the nucleotide sequences of each read with a measure of the quality of the nucleotide identification. Once, the DNA of samples is extracted and prepared, different genomic approaches can be followed to prepare the DNA library. On one hand, the whole genome sequencing (WGS) procedure determines the nucleotide order in the target genome in one run (genome resequencing). Nevertheless, this option was until recently prohibitive on a population scale (Luca et al. 2011). On the other hand, using reduced representation libraries (RRLs) to work with a fraction of the whole genome in a more cost-effective way. With this approach reduced subsets of loci can be obtained with exome capture, transcriptome CASANOVA CHICLANA, ADRIÁN 24 sequencing (i.e. RNA-seq), etc (Hirsch et al. 2014). A subset of loci can be obtained through restriction enzyme digestion, followed by highthroughput sequencing of genomic sequences adjacent to the enzyme cut-site (restriction site-associated DNA sequencing, RAD-seq). A variety of RAD-seq techniques has been developed (2bRAD, ddRAD, etc; see Andrews et al. 2016) from the original publication (Baird et al. 2008). Some reasons that fuelled the popularity of RAD-seq were its feasibility for different genomes and the non-need of a reference genome, enabling studies of non-model organisms (Davey et al. 2011). 1.4.2 Molecular markers Molecular population genetics was born 50 years ago with the first measures of genetic variation in allozyme loci (i.e. protein variants), the first molecular markers (Casillas and Barbadilla 2017). Molecular markers in sexually reproducing organisms can have biparental inheritance (nuclear DNA) or uniparental inheritance (mitochondrial DNA and plastids). Molecular markers are usually co-dominant meaning higher accuracy and informativeness, enabling to obtain genotypic and allelic frequencies in populations, as the basic information for their later application using the statistical methodologies of population genetics. Microsatellite loci are tandem repetitive DNA of two to six base pairs motifs occurring at different genomic locations in a high frequency (on average 1 microsatellite/30 kb). SNPs could theoretically have a maximum of four possible variants, although nearly all SNPs are biallelic due to the low mutation rate per nucleotide site (Brown 2018). Microsatellites show much higher variation due to their much higher mutation rate (10-3 - 10-5; Castro et al. 2006) thus having a much higher number of alleles per locus. Nevertheless, SNPs are much more frequent than microsatellites throughout the genome (1 SNP/300 bp in the human genome), thus accounting for approximately 90% of the total genetic variation (Collins et al. 1999). General introduction 25 To date, the markers used to estimate different population metrics in Spanish brown trout populations (allozymes, microsatellite loci or mitochondrial DNA control region; see Garcia-Marín et al. 2018) are assumed to be selectively neutral, reflecting only demographic processes from the balance between gene flow and genetic drift. 1.5 BROWN TROUT POPULATION GENOMICS Different studies using genomic strategies have been reported to date in European brown trout, some of them using SNP-arrays (see Bekkevold et al. 2019; Linløkken et al. 2017). The diverse topics addressed are related to assessing genetic structure and identification of conservation units, restocking evaluation and monitoring, detection of adaptive variation (i.e. signatures of selection) in a wide range of features (e.g. individual growth) and linkage mapping, among others. Linløkken et al. (2017) using a 3,781 SNP panel detected outlier SNP loci under selection linked to genes involved in growth. Bekkevold et al. (2019) with a similar number of SNPs refined population structure and identified signatures of selection and their relationships to environmental factors. Leitwein et al. (2018) with 75,684 SNPs evaluated the impact of restocking from domestic Atlantic and Mediterranean lineages into supplemented wild populations across the genome using a local ancestry approach. Lemopoulos et al. (2018) with ~20,000 SNPs identified eight candidate genes potentially associated with brown trout migratory behaviour. Leitwein et al. (2017) with about 20,000 SNPs constructed a linkage map for brown trout to study chromosomal rearrangements, recombination rates and the effect of selection on neutral diversity. Bernaś et al. (2020) with 3,843 SNPs studied population structure and detected outlier loci in Southern Baltic region. Finally, Jacobs et al. (2018) analysed the evolutionary history between two reproductively isolated life history morphs in Scotland obtaining genes potentially related to observed phenotypic differences. CASANOVA CHICLANA, ADRIÁN 26 To our knowledge, no research following a population genomics approach has been addressed to date with brown trout in the Iberian Peninsula. Until 2015, the main way to build loci and call genotypes with bioinformatic tools in brown trout using RAD-seq was a de novo approach (i.e. without reference genome). The genome of Atlantic salmon (Salmo salar), a congeneric species, was available since 2015 from the International Cooperation Project to Sequence the Atlantic Salmon Genome initiative (ICSASG). Reference genome approaches were carried out in brown trout RAD-seq data using the Atlantic salmon as reference genome (see Lemopoulos et al. 2018, 2019b; Paris et al. 2017; Valette et al. 2020). While reference genomes of closely related species can be used to obtain SNPs to estimate different population genetics parameters (e.g. nucleotide diversity; see Galla et al. 2019), accurate requires the species' own reference genome. For instance, conspecific genomes are recommendable to obtain the highest possible resolution (i.e. number and position of SNPs) to study patterns of hybridization along genome. In June 2019 was available at National Center for Biotechnology Information (NCBI) the first brown trout genome. This achievement was made within the framework of the “25 genomes for 25 years” initiative and promoted by the Wellcome Sanger Institute (https://www.sanger.ac.uk/collaboration/25-genomes-25years/; Hansen et al. 2021). Because it has been available for just over a year, there are very few published research using brown trout reference genome (see De-Kayne et al. 2020; Grimholt and Lukacs 2021; Sudhagar et al. 2020). Objectives 27 2 OBJECTIVES The present thesis aimed at studying the genetic structure of brown trout (Salmo trutta) in the Atlantic and Mediterranean slopes of the Iberian Peninsula using a population genomics approach. Natural and artificial hybridization and introgression patterns were studied and all the information was considered for the management and conservation of native genetic resources of this species in the Iberian Peninsula. The specific objectives of the study were: 1. Evaluation of different bioinformatic pipelines to identify and genotype SNPs to achieve consistent biological interpretations through a population genomics approach. 2. Analysis of genomic introgression of wild populations by hatchery stocks in the unstable ecosystems of the Mediterranean drainage. 3. Preliminary study of natural hybridization in a secondary contact in the Atlantic slope. 4. Improvement of conservation and management strategies to maintain the integrity of native lineages and recover natural populations in areas with high introgression rate due to restocking. CASANOVA CHICLANA, ADRIÁN 28 3 LOW IMPACT OF DIFFERENT SNP PANELS FROM TWO BUILDING-LOCI PIPELINES ON RAD-SEQ POPULATION GENOMIC METRICS: CASE STUDY ON FIVE DIVERSE AQUATIC SPECIES 3.1 INTRODUCTION Next-generation sequencing technologies have represented a breakthrough for genomic studies (Quail et al. 2012) due to the huge reduction of sequencing cost (less than 0.02$ per Mb; Wetterstrand 2020) and the development of a broad and versatile range of techniques for different genomic approaches (Goodwin et al. 2016). By harnessing the possibilities of NGS, diverse reduced-representation genome sequencing approaches, useful to identify and genotype thousands of markers for genomic screening, were suggested and quickly became popular (Baird et al. 2008; Davey et al. 2011). One of these approaches is the restriction site-associated DNA sequencing (RAD-seq), currently in a more mature phase, which includes different methods (e.g. ddRADseq, ezRAD-seq, 2bRAD-seq) whose performances have been compared using simulations and real data (Andrews et al. 2016). RADseq methods require specific library preparation protocols, which exploit the ability of restriction enzymes (REs) to cut at specific genomic targets rendering a collection of fragments representative of a Low impact of different SNP panels from two building-loci pipelines on RAD-seq population genomic metrics: case study on five diverse aquatic species 29 genome fraction to be compared among samples. These collections can be screened to identify and genotype a variable number of SNPs depending on the goals of the study for population genomics, linkage mapping or genome wide association studies, among others. The 2bRAD method here used exploits the properties of IIB REs which produce a collection of short DNA fragments (between 33 and 36 bp) by cutting at both sides of the recognition site (Wang et al. 2012). This method has the advantages of simple library preparation, short-reads to be sequenced (single-end 50 bp) and, as other methods, the number of loci can be adjusted both using REs with different recognition site frequency or by fixing nucleotides in the adaptors during library construction (i.e. selective-base ligation) (Barbanti et al. 2020; Wang et al. 2012). Genomic laboratory protocols have been set up and optimized through years by introducing modifications on the original RAD-seq methodology to get better results using different laboratory protocols for different scenarios (e.g. samples with low DNA quality, genome size, etc; see Fig. 5 in Barbanti et al. 2020). Similarly, the bioinformatic pipelines starting from raw data, a critical issue in RAD-seq methodologies, have undergone an important refinement and diversification. Nevertheless, there is not a consensus about what is the best strategy for each scenario, despite the increasing number of studies addressed to evaluate the impact of technical and/or bioinformatic protocols (Díaz-Arce et al. 2019; O’Leary et al. 2018). In a typical 2bRAD library, hundreds of millions of reads are generated, and they need to be allocated to each multiplexed individual (dozens to hundreds in the same lane) and to each genomic position or locus in the reference genome (or RAD-tag catalogue). The rationale behind this is stacking raw reads belonging to the same locus, while discerning and separating at the same time the reads belonging to different loci. Results could be improved if a reference genome, belonging to the species itself or to CASANOVA CHICLANA, ADRIÁN 30 other congeneric species, is available. This would enable to avoid mixing of reads pertaining to paralogous loci. In November 2020, there were reference genomes for 25 bivalve species and subspecies (22 genera) and 583 fish species (338 genera) with different assembly confidence at the NCBI database (https://www.ncbi.nlm.nih.gov/datasets/). Nevertheless, there are about 9,200 species within the 1,260 bivalve genera (Huber 2010) and 35,672 recognized species within the 5,212 documented fish genera (Fricke et al. 2020). All in all, less than 0.2% of the genomes of the known eukaryotic species have been sequenced to date (Lewin et al. 2018). Although full genome sequencing assembly is becoming progressively more robust thanks to the long-read sequencing methods and assembling strategies, most of the species will have to wait for long before their genomes are assembled. Therefore, de novo approaches (i.e. stacking reads without a reference genome) will be the only option for many studies, although some initiatives are trying to change this perspective (e.g. Earth Biogenome Project; https://www.earthbiogenome.org/). For this reason, one of the strengths of a RAD-based method is its applicability without a reference genome (Rochette and Catchen 2017). There are different bioinformatic pipelines to identify a high number of SNPs and achieve confident genotypes using a RAD-seq approach. The most popular one is Stacks (Catchen et al. 2011, 2013) with around 3,000 citations, at Google scholar in Nov. 2020, but several other alternatives, including Meyer’s 2b-RAD pipeline (Wang et al. 2012), which was the original building-loci pipeline for 2b-RAD data, have been recently published. Some of these pipelines are able to perform a de novo approach (dDocent; Puritz et al. 2014), whereas others need a reference genome for alignment (Fast-GBS, Torkamaneh et al. 2017; TASSEL-GBS v2, Glaubitz et al. 2014) or can address both approaches (Stacks, Meyer’s 2b-RAD v2.1 pipeline, ipyrad; Eaton and Overcast 2020). Several of these alternative pipelines merge and Low impact of different SNP panels from two building-loci pipelines on RAD-seq population genomic metrics: case study on five diverse aquatic species 31 concatenate pre-existing applications, making their design flexible and customized according to the data managed and the goals of the study, but also providing upgrading and reliable bug-fix (e.g. dDocent, FastGBS, Meyer’s 2b-RAD v2.1 pipeline). Several factors should be considered for the selection of the bioinformatic pipeline to be used, among which sampling variance (number of samples and reads across them), population structure, genome architecture of the species studied and budget (e.g. read coverage) are the most relevant. The genome of each species has its particular size, history (e.g. duplication events), polymorphism, complexity and interindividual variability, which can hinder the identification of stacks of reads (putative RAD-loci) and their variants, circumstances that should be considered when choosing the appropriate building-loci pipeline and its parameters. Studies comparing bioinformatic pipelines and strategies already exist. Some comparisons between de novo and reference-based approaches are available (Shafer et al. 2017; Torkamaneh et al. 2016), and one of them tested the performance of the different strategies used to obtain accurate population genetics inferences (Shafer et al. 2017). Noticeable differences were observed among bioinformatic pipelines in the number of detected SNPs, sometimes resulting in distinct values for population descriptors and inference (Shafer et al. 2017). Other studies have evaluated the same software with different species to optimise the selection of bioinformatic parameters (Stacks 1.42; Paris et al. 2017); (Stacks 1.44; Díaz-Arce et al. 2019), making a common advice of doing preliminary trials to optimize the building-loci pipeline selected parameters. Published step-by-step protocols with a single species also exist (Rochette and Catchen 2017). A number of SNP calling comparison between Stacks 1.08 and dDocent 1.0 has been carried out using three fish species (Puritz et al. 2014), while Sovic et al. (2015) tested a novel pipeline (i.e. AftrRAD) vs Stacks and PYRAD using simulated and species datasets to assess computational efficiency and CASANOVA CHICLANA, ADRIÁN 38 A reference genome-based pipeline was used for Manila clam and brown trout, the species with chromosome assembly level reference genomes. In this case, the differences in the pipelines of Stacks 2 and Meyer’s 2b-RAD v2.1 are lower. For instance, Stacks 2 reduces the number of modules necessary from six to two when using reference genome and the building of putative loci is conditioned by the step using a short-read aligner. Our goal was not to compare this option between pipelines, but to take as reference the genome-based approach to be compared with the de novo approach as a gold standard within each pipeline. Building-loci pipelines: analysis After demultiplexing raw data, several filtering criteria were applied: (1) all reads were trimmed and filtered by the RE recognition site to retain only those sequences of 36 bp (or 32 bp from CII RE catshark) centred on the RE recognition site using our own Perl scripts and Trimmomatic 0.38 (Bolger et al. 2014); and (2) process_radtags (module belonging to Stacks) was used to remove reads with uncalled bases (-c option). Other parameters were species-specific (e.g. window sliding size, -w; score limit, see Table S2A). To check raw and filtered reads quality, FastQC 0.11.7 (Andrews 2010) was used. Stacks input sequences were oriented in the same orientation using our own Perl script to avoid oversplitting. To say, the set of input reads for both building-loci pipelines in each species was always the same. At the building-loci pipeline step, the parameters considered were: (1) minimum number of identical reads to create a stack (default values were used); (2) maximum number of mismatches between RAD-loci within and between individuals (-M 2/-n 2 for fish or -M 3/-n 3 for bivalves and their analogous parameters with ALT pipeline); (3) indels were discarded (i.e. --disable-gapped at different modules); (4) SNP calling model and alpha cut-off: their default values were used in the STA pipeline (Stacks 2.0), while in the ALT pipeline (Meyer’s 2b-RAD Low impact of different SNP panels from two building-loci pipelines on RAD-seq population genomic metrics: case study on five diverse aquatic species 39 v2.1) we considered a range of 0.1–0.2 to determine the genotype at each position (default values are 0.01–0.25); to say, when the frequency of the less frequent allele was lower than 0.1 the genotype was called as homozygous while frequencies higher than 0.2 were called as heterozygous; intermediate allele frequencies for the less frequent allele were called as uncertain (see Tables S2B and S2C). When a reference genome was available, Bowtie 1.1.2 was used as short read aligner. The number of mismatches allowed between reads and the reference genome using the -v alignment mode was 2 mismatches for brown trout and 3 mismatches for Manila clam, the same as mentioned above. Only reads which aligned to a single site in the reference genome were considered (-m 1). The same parameter values were used in Stacks modules shared between reference-based genome and de novo approaches (i.e. gstacks and populations modules). SNP filtering steps and creation of datasets The raw SNP panels of the STA and ALT pipelines were filtered using the same parameters for consistency, retaining only a set of markers and alleles represented across the individuals genotyped. Filters were applied in the same order for each dataset (Fig. 8), as recommended (O’Leary et al. 2018): (1) Uniquely biallelic SNPs were included (BIAL filter); (2) a minimum locus coverage of eight reads was chosen (Min coverage filter); (3) RAD-loci with ≤ 3 SNPs were retained; (4) SNPs were retained only if the less frequent allele was represented at least three times at the whole species sample (MAC, minimum allele count, filter); (5) genotyped in at least in 60% of the individuals in each population for a SNP to be retained (POP filter); (6) SNPs were included when they adjust to Hardy-Weinberg expectations (p-value > 0.01) in more than half of the populations analysed (HW filter); and (7) only the first SNP per RAD-locus was retained when several SNPs were called in the same RAD-locus to avoid redundant information. CASANOVA CHICLANA, ADRIÁN 40 Figure 8. Scheme of filtering steps to obtain the different SNP panels. STA, ALT, COM and MER, representing Stacks, Alternative, Shared and Merged panels, respectively. According to the aforementioned criteria, four SNP panels were tested and compared: (1) STA and (2) ALT SNP panels were further used to obtain (3) common (COM) and (4) merged (MER) panels. Low impact of different SNP panels from two building-loci pipelines on RAD-seq population genomic metrics: case study on five diverse aquatic species 41 When reference genome was available three additional SNP panels were obtained (i.e. RG, RG-STA, RG-ALT). RAD-loci of these panels from both pipelines were compared to identify shared and private RADloci. In order to do this, cd-hit-est was used to cluster similar RAD-loci taken from the two pipeline catalogues with the same threshold of similarity (-c) used in the clustering steps of building-loci pipelines (i.e. two mismatches maximum for fish species and three mismatches maximum for bivalve species). Furthermore, we used a -g value of 1 when clustering sequences to meet the established similarity threshold and a band alignment width (-b) of 1 to avoid previously separated sequences due to indels at specific clusters. This procedure rendered COM and MER SNP panels created using a customized Perl script. SNPs at shared RAD-loci between STA and ALT panels were selected after the sixth filtering step, since the first SNP at shared RAD-loci could be different after the full filtering pipeline (Fig. 8). MER panel was finally created by taking the SNPs from “private” ALT and STA pipelines RAD-loci (i.e. those from cdhit-est clusters with only STA or ALT pipelines RAD-tags) plus the COM SNP panel previously obtained. Again, only the first SNP per RAD-locus was retained to avoid redundant information. The Genepop files of all shared SNP panels were compared to quantify their genotyping differences using own Perl script (see https://github.com/abhortas/USC-RAD-seq-scripts). Comparison of outputs and population genetics analyses The results of the aforementioned pipelines were compared using both different quantitative (i.e. number of SNPs) and population genomics metrics. Filtered Genepop files were obtained using customized Perl scripts. These files were transformed for subsequent analyses using the PGDSpider 2.1.1.5 software (Lischer and Excoffier 2012). Firstly, the consequences of filtering over the number of RADloci and SNPs were evaluated for each combination of species-pipeline; CASANOVA CHICLANA, ADRIÁN 42 secondly, common, and private RADloci/SNPs between the two pipelines were obtained for each species. Finally, biological interpretations from each pipeline/species were compared, including basic population genetics results (i.e. genetic diversity levels and population structure). Observed and expected heterozygosity (Ho and He, respectively), inbreeding coefficient (FIS, using 1,000 bootstrap iterations to estimate their 95% confidence intervals) and allelic richness were calculated per population using diveRsity R package 1.9 (Keenan et al. 2013). Global FST calculation and HW tests were performed with Genepop R package 1.1.7 (Rousset 2008). STRUCTURE 2.3.4 software (Pritchard et al. 2000), using R package ParallelStructure 1.0 (Besnier and Glover 2013), was used to define the most likely number of population units (K) present with LOCPRIOR model with correlated allele frequencies model, testing K values from 1 to the number of sampling localities in the species dataset + 1 with 10 replicates composed by 100,000 Markov chain Monte Carlo (MCMC) replicates and a burn-in period of 10,000 steps. STRUCTURE results were parsed with STRUCTURE HARVESTER (Earl and vonHoldt 2012), which implements the Evanno’s method (Evanno et al. 2005) to detect the most likely number of clusters according to the data. CLUMPAK (Kopelman et al. 2015) was used to merge runs with the same K that suggested similar patterns of structuring and to obtain cluster membership plots. As a second approach to detect population structure, a Discriminant Analysis of Principal Components (DAPC) analysis was performed based on genetic data, as implemented in R package Adegenet 2.1 (Jombart 2008; Jombart and Ahmed 2011). The optimal number of principal components to be used was estimated with the cross-validation method implemented in the package and from one to three discriminant components were retained according to the amount of population structure variation they explained. Finally, outlier loci potentially under Low impact of different SNP panels from two building-loci pipelines on RAD-seq population genomic metrics: case study on five diverse aquatic species 43 selection (OL), i.e. those showing higher or lower differentiation values (i.e. FST) across populations than the neutral background, were detected using the Bayesian approach implemented in BayeScan 2.1 (Foll and Gaggiotti 2008) with default parameters. Loci with a Log10 posterior odds (PO) higher than 1.5 were retained as potential outliers for later comparison among the four datasets resulting from the two pipelines. 3.3 RESULTS The number of filtered reads loaded into building-loci pipelines using the reference genome approach was lower than with the de novo approach. A percentage of 22.1 and 25.4% were finally used in Manila clam and in brown trout, respectively. This reduction mostly due to those reads aligning to more than one place (59.6 and 72.6%, respectively) that were filtered out (i.e. -m 1 in Bowtie 1.1.2), the remaining reads failed to align with the mismatch criteria applied (18.3 and 2.0%, respectively). The number of initial SNPs, after the building-loci step with the de novo approach, ranged from 56,074 in brown trout (STA) and 125,823 in silver catfish (ALT) to 356,389 (STA) and 426,317 (ALT) in common cockle (Tables S3-S7). These figures dropped throughout the successive filtering steps (Fig. 8) up to finally being retained from 0.2% in Manila clam STA panel to 20.5% in silver catfish STA panel (Fig. 9). CASANOVA CHICLANA, ADRIÁN 44 Figure 9. Number of SNPs from the initial building-loci pipelines (blue bars) to the final panels (green bars) through the different SNP filtering steps. There was a remarkable difference at the initial number of SNPs obtained between STA and ALT pipelines in brown trout and smallspotted catshark, although the outcome after filtering was rather similar (Table 1). No comparison was made at BIAL filter (i.e. SNPs with more than two alleles excluded; Fig. 8), since triallelic SNPs are removed with Stacks by default. The proportion of missing genotypes after applying the minimum coverage filtering step was higher with ALT panel than with STA in almost all species (Fig. 10). Filtering patterns varied among species due to the different weight of each filtering step. For instance, in bivalves, where genetic polymorphism is higher, the SNP retention after the third filtering step (i.e. RAD-loci with ≤ 3 SNPs per RAD-locus) was much lower than in fish species (Fig. 9), while in silver catfish and small-spotted catshark, was due to the minimum allele count (MAC). This was related to the smaller sampling size of those fish species (N = 21 and N = 28, respectively) and the higher frequency of missing data, especially in small-spotted catshark (83% after depth filter in STA and ALT panels). When comparing pipelines, no clear differences on the filtering pattern were observed through the different filtering steps (Fig. 9), except for MAC in brown trout, where more SNPs were pruned in the ALT panel. This could be related to the Low impact of different SNP panels from two building-loci pipelines on RAD-seq population genomic metrics: case study on five diverse aquatic species 45 increment on the missing data after the minimum coverage filter, with higher frequency in ALT (59.3% vs 43.7% for ALT and STA, respectively (Fig. 10). After the third filtering step, in brown trout (N = 52) there were significantly more missing genotypes per SNP (pvalue < 0.01) at ALT panel on average (28.74 ± 15.24) as compared to STA panel (21.07 ± 17.54). The two species with the lowest median coverage at this step were small-spotted catshark (median = 10x for ALT and STA panels), brown trout (median = 11x and 14x for ALT and STA panels, respectively) and Manila clam (21x for the STA panel). The final number of SNPs ranged from 479 (STA) and 956 (ALT) in Manila clam to 21,468 (STA) and 22,481 (ALT) in silver catfish. These figures were always higher with the ALT pipeline for all species. The number of SNPs in COM panels, those from STA panel shared with ALT panel, ranged from 206 in Manila clam to 17,459 SNPs in silver catfish, while the percentage of SNPs called in STA that were found in ALT ranged from 23.9% in small-spotted catshark to 81.3% in silver catfish. The main source of variation when considering the whole COM panels genotype dataset (for instance in silver catfish Ntotal COM genotypes (366,639) = Nsamples (21) x NCOM SNP (17,459)) was missing data, ranging from 2.6% in common cockle to 14.8% at small-spotted catshark (Figs. S1-S4). Roughly speaking, there would be three types of missing calling data with a de novo approach: (1) SNPs not genotyped by a building-loci pipeline due to too low coverage (< 3x with used configuration); (2) SNPs genotyped by a building-loci pipeline but with a coverage lower than 8x (Min coverage 8x filtering step; Fig. 8); and (3) SNPs with enough coverage but ambiguous alternative nucleotide depth (see http://elimeyer.github.io/2bRAD_utilities/#genotype). The main source of missing genotype differences between pipelines was related to COM SNPs from STA pipeline that passed the minimum coverage filter (Min coverage 8x), but not were genotyped by ALT pipeline due to having a coverage lower than 3x. Both pipelines never genotyped with a CASANOVA CHICLANA, ADRIÁN 46 coverage lower than 3x with the used configuration. This situation was found in most species, causing from 46.8% of the total missing genotype differences between ALT and STA genotypes in shared SNPs (COM panel) in brown trout to 63.2% in small-spotted catshark. Figure 10. Percentage of missing genotypes after the Min coverage filter (8x). Boxplots were obtained with the percentage of missing genotypes through the different samples. Up and down triangles represent the percentage of missing genotypes at different SNP panels after and before 8x coverage filter, respectively. Nevertheless, the main source of missing data differences in Manila clam was genotypes removed by coverage filter (Min coverage Low impact of different SNP panels from two building-loci pipelines on RAD-seq population genomic metrics: case study on five diverse aquatic species 47 8x) in ALT but not in STA homozygous-heterozygous differences between pipelines at the same SNP and individual ranged from 0.5% in brown trout to 2.6% in small-spotted catshark with respect to the whole COM genotype panel (Table 2). The frequency of genotyping differences caused by different homozygotes at the same SNP and individual (e.g. AA for one pipeline and GG for the other) was negligible in almost all cases (from 0 to 0.098%). The number of SNPs obtained with reference genome (RG) approach was always lower than with both de novo building-loci pipelines (see Tables S3 and S5), however, the number of SNPs shared between RG and each of the de novo pipelines was similar (i.e. RG-STA and RG-ALT). The percentage of SNPs obtained with RG approach detected as well in both de novo pipelines was 47.7% (RG-STA/ RG) and 37.7% (RG-ALT/RG) in Manila clam and 73.0% (RG-STA/RG) and 80.9% (RG-ALT/RG) in brown trout, unlike the reverse way where the percentages of shared SNPs were lower due to the higher number of SNPs obtained with ALT: 28.8% (RG-STA/STA) and 11.4% (RGALT/ALT) for Manila clam with (Table S3) and 46.5% (RG-STA/STA) and 32.9% (RG-ALT/ALT) for brown trout (Table S5). Both de novo building-loci pipelines performed similarly when compared with reference genome (RG) approach genotyping (RG-STA and RG-ALT genotype comparisons; Table 2). The population parameters evaluated (e.g. diversity levels, global FST) were roughly similar when using the different SNP panels in each species, showing higher differences in FIS values (Table 1). The most notable differences would be found in the FIS values and when comparing de novo and reference genome approaches in brown trout, especially regarding Hardy-Weinberg tests and related population parameters (i.e. FIS, Ho vs He; Table S5). Here, unlike Manila clam, the proportion of SNPs with extremely low FIS (≤ - 0.5) greatly differed between both approaches. The structure patterns obtained using CASANOVA CHICLANA, ADRIÁN 54 exploitation. According to ichthyoarchaeological analyses, brown trout consumption was practiced since the Neolithic in human communities of the Sierra de Atapuerca, where brown trout is common nowadays (Blanco-Lapaz and Vergès 2016). The existence of brown trout as a species occurred around 10-14 Ma (Lecaudey et al. 2018) and the divergence in the S.trutta complex took place during the Pliocene, around 2.5-5 Ma (Crête-Lafrenière et al. 2012). Since then, brown trout populations occupied the Iberian rivers through the recurrent glaciation events associated to the Pleistocene, when the Iberian Peninsula was an important refugia for many species due to its warmer climatic conditions, which later expanded northwards from this southern European refugium (Weiss 2010). The analysis of mtDNA sequence in the last three decades has allowed identifying four Iberian native lineages: Duero (DU), Atlantic (AT), Adriatic (AD) and Mediterranean (ME). The first one would be an endemism of the Iberian Peninsula (Vera et al. 2010), located in the inner sections of the Duero (Bouza et al. 2008; Hermida et al. 2009; Martínez et al. 2007; Vera et al. 2010) and Miño-Sil basins (Vera et al. 2015; Vilas et al. 2010). These mitochondrial lineages are characterized by different haplotypes (Cortey et al. 2009; Vera et al. 2010). Moreover, these four matriarchal lineages are geographical distributed. DU and AT lineages are distributed in the western slope of Iberian Peninsula, while AD and ME lineages in the Mediterranean one. DU and AT brown trout are distributed in the Miño and Duero basins following a spatial segregation in parapatry with a not well-defined hybrid zone between them (Bouza et al. 2001; Martínez et al. 2007). This hybrid zone would be the consequence of secondary contacts (Bouza et al. 2001) between the endemic DU lineage, ancestral in the Iberian Peninsula, and the nowadays more abundant AT arrived in later expansions from the north (García-Marín et al. 2018). On the other hand, AD and ME lineages are sympatrically distributed, although following a mosaic patchiness in the Population genomics for conservation of brown trout resources in the Mediterranean and Atlantic slopes of Iberian Peninsula 55 Mediterranean slope (Cortey et al. 2004; Sanz et al. 2002; Vera et al. 2013). Resident and migratory (sea trout) forms of brown trout are present in Galicia, since the Miño outlets at parallel 42º N, the distribution limit of this migratory form (Bouza et al. 1999). In addition to the phylogeographic variation of brown trout in the Iberian Peninsula previously described, there is a foreign genetic lineage due to releases of hatchery brown trout using a stock from Central European origin commonly used for restocking practices throughout the European continent. During decades, millions of fertile brown trout individuals have been released as eggs or fry into Spanish rivers to satisfy angling demand and counterbalance the decline of wild populations. The hatchery stock used was founded with the aforementioned strain belonging to an AT clade different than the wild AT lineage living in Iberian Peninsula (Cortey et al. 2009; Machordom et al. 2000). The genetic consequences of this restocking activity have been studied with allozyme markers by different authors, who detected different degrees of introgression throughout the Iberian Peninsula (Almodóvar et al. 2006; Martínez et al. 1993; Morán et al. 1991). The genetic introgression detected in Galician rivers using the diagnostic marker LDH-C*90 was low, suggesting very low viability of hatchery individuals in Galician rivers (Arias et al. 1995; Martínez et al. 1993). Nevertheless, a much higher impact of restocking was detected in nonflowing waters (lagoons and reservoires; Martínez et al. 1993). The lower impact of restocking in the Atlantic Drainage has been related to the stable hydrological conditions of this region regarding Mediterranean Drainage, where the impact of restocking is higher (Almodóvar et al. 2006; Vera et al. 2013). For conservation purposes the first is to define what to conserve. There are two main concepts: Management Units (MUs; i.e. populations within species that are genetically distinct enough to require separate management) and Evolutionary Significant Units CASANOVA CHICLANA, ADRIÁN 56 (ESUs), “whose divergence can be measured or evaluated by putting differential emphasis on the role of evolutionary forces at varied temporal scales” (Casacci et al. 2014). For the latter there are different definitions from the original by Ryder (1986) that can be inspected in Casacci et al. (2014; see Table 1), but all have in common the sentence in quotation marks outlined above. MUs and ESUs are not synonymous concepts (see Moritz 1994). To define specific management units, population genetic structure approaches are necessary. Furthermore, in species such as brown trout characterized by high population structure due to different factors (e.g. habitat patchiness; Ferguson 1989), defining MUs is a matter of controversy, since populations separated by a few kilometres show significant differentiation (Bouza et al. 1999; Carlsson and Nilsson 2001). This issue was the goal of previous studies with allozymes (Bouza et al. 1999, 2001; Martínez et al. 1993), microsatellites (Vilas et al. 2010) and mtDNA markers (Bouza et al. 2008). Two different mitochondrial lineages of brown trout co-occur in Galicia parapatrically segregated, DU and AT, such as in Duero Basin (Bouza et al. 2001; Vera et al. 2010). This distribution has been hypothesized as the consequence of a secondary contact between those two divergent lineages, being particularly remarkable in both basins, where putative hybrid zones would exist (Bouza et al. 2001; Martínez et al. 2007; Vilas et al. 2010). Within DU lineage, an endemism of the Iberian Peninsula (Machordom et al. 2000), there are different haplotypes with significant divergence, advising against treating this lineage as a single ESU (Vera et al. 2015). Within Mediterranean Slope, in Catalonia there is the figure of conservation of ‘genetic refuges’, where releases from hatchery stock are banned (Araguas et al. 2017). In other basins, it has been proposed to define additional refuges to protect remnant native brown trout lineages (AD and ME; Vera et al. 2013). Population genomics for conservation of brown trout resources in the Mediterranean and Atlantic slopes of Iberian Peninsula 57 Over the last century, specific legislation has been developed in Spain concerning the regulation of angling, environmental protection, and the protection of biodiversity. Different laws were developed from the most specific, as Salmon Protection Law (1912), to those with larger spectrum such as Inland Fisheries Law (1907, 1929, 1942). These laws reflected the need for the regularization of the activity with a conservationist perspective. Also, within the European Union, the Council Directive 78/659/EEC of 18 July 1978 reflected the need for good quality fresh water to make life possible for the different fishes, with mention to salmonid and cyprinid water bodies. This Directive was replaced by Directive 2000/60/EC (Directiva Marco del Agua; DMA), designed to protect and improve freshwater quality, setting standards for defining water quality based on its "ecological status”, using among other biological quality indicators. Over the years, the preservation of native genetic biodiversity has been incorporated into increasingly advanced laws as research has progressed (Box 1). By the same way, regulatory and conservationist legislation has been implemented in different autonomous communities (see Ley 22/2009 de 23 de diciembre, de ordenación sostenible de la pesca en aguas continentales de Cataluña; Box 1). Regarding freshwater angling in Galicia, during almost thirty years, the Ley 7/1992 del 24 de julio, de pesca fluvial de Galicia have been active. Recently, a new law was approved at the Galician Parliament (i.e. Ley 2/2021, de 8 de enero, de pesca continental de Galicia). There are some highlights respect to previous one (Box 2). CASANOVA CHICLANA, ADRIÁN 58 Box 1. Selection of legislative highlights at different Public Administration levels. Artículo 45 de la Constitución Española (1978). “1. Todos tienen el derecho a disfrutar de un medio ambiente adecuado para el desarrollo de la persona, así como el deber de conservarlo. 2. Los poderes públicos velarán por la utilización racional de todos los recursos naturales, con el fin de proteger y mejorar la calidad de la vida y defender y restaurar el medio ambiente, apoyándose en la indispensable solidaridad colectiva. 3. Para quienes violen lo dispuesto en el apartado anterior, en los términos que la ley fije se establecerán sanciones penales o, en su caso, administrativas, así como la obligación de reparar el daño causado.” Directiva del Consejo 78/659/EEC de 18 de julio de 1978, relativa a la calidad de las aguas dulces que deben protegerse o mejorarse para mantener la vida de los peces. Artículo 1.1: “La presente Directiva trata de la calidad de las aguas continentales y se aplicará a las aguas que requieren protección o mejora para ser aptas para la vida de los peces, declaradas como tales por los Estados miembros.” Ley 42/2007, de 13 de diciembre, del Patrimonio Natural y de la Biodiversidad. Artículo 67: “El Inventario Español de Caza y Pesca, dependiente del Ministerio de Medio Ambiente, mantendrá la información más completa de las poblaciones, capturas y evolución genética de las especies cuya caza o pesca estén autorizadas, con especial atención a las especies migradoras.” Ley 22/2009 de 23 de diciembre, de ordenación sostenible de la pesca en aguas continentales de Cataluña. Artículo 29.2: “Las zonas de pesca controlada intensiva deben ubicarse en cursos, tramos de cursos o masas de agua transformados artificialmente, en especial los embalses, y fuera de las aguas de reserva genética, para evitar la degradación biológica y, en especial, genética de las poblaciones de especies autóctonas.” Artículo 47.3: “No pueden emplazarse nuevos centros industriales de producción de fauna en aguas continentales en derivaciones de tramos de cursos de agua que hayan sido declarados refugios de pesca ni en aguas de reserva genética.” Population genomics for conservation of brown trout resources in the Mediterranean and Atlantic slopes of Iberian Peninsula 59 Box 2. Selection of highlights from Ley 2/2021, de 8 de enero, de pesca continental de Galicia. Artículo 10 (Investigación en materia de pesca continental): “La consejería competente en materia de pesca continental impulsará la mejora del conocimiento sobre la etología, biología y dinámica poblacional de las especies de la fauna acuática, en especial de las pescables. Asimismo, impulsará el conocimiento genético de las poblaciones ictícolas y del impacto de las especies exóticas invasoras sobre el ecosistema acuático, y la mejora de los métodos de gestión de la pesca continental.” Artículo 61.4 (Sueltas): “Las sueltas se realizarán con especies autóctonas y con ejemplares nacidos en libertad o procedentes de centros ictiogénicos dependientes de la consejería competente en materia de pesca continental y obtenidos de reproductores capturados en la misma cuenca hidrográfica en la que se va a realizar la suelta o, en su defecto, con ecotipos de la mayor similitud genética posible.” Artículo 62.2 (Repoblaciones piscícolas): “Las repoblaciones piscícolas se realizarán con especies autóctonas y con ejemplares nacidos en libertad o procedentes de centros ictiogénicos dependientes de la consejería competente en materia de pesca continental y obtenidos de reproductores capturados en la misma cuenca hidrográfica en la que se va a realizar la repoblación o, en su defecto, con ecotipos de la mayor similitud genética posible.” Artículo 62.4 (Repoblaciones piscícolas): “No podrán repoblarse aquellos tramos de agua en los que habiten poblaciones piscícolas de interés por sus peculiaridades biológicas o genéticas, así como aquellos tramos de agua en los que exista algún régimen de protección especial, salvo por razones de defensa de las poblaciones, debidamente justificadas.” Artículo 63.1 (Centros ictiogénicos): “Se declaran de interés general los centros ictiogénicos para el fomento de la recuperación y conservación de las poblaciones piscícolas salvajes y del medio en el que se desarrollan.” CASANOVA CHICLANA, ADRIÁN 60 The specific aims of this chapter were: (1) Identification of farmed and wild samples and development of diagnostic SNP panels; (2) to analyse patterns of human-mediated introgression across the genome; (3) to analyse the pattern of genetic diversity and structure through populations, and (4) to find traces of natural selection. Population genomics for conservation of brown trout resources in the Mediterranean and Atlantic slopes of Iberian Peninsula 61 4.2 MATERIAL AND METHODS 4.2.1 Sampling Thirteen populations of brown trout (Salmo trutta) from Miño and Duero River basins draining into the Atlantic Ocean, and Ter Basin into the Mediterranean Sea from the Iberian Peninsula were sampled (Table 3). All populations from Miño and Duero basins were collected in tributaries (Fig. 11). Populations belonging to the Atlantic Slope included representatives of the endemic DU lineage (LE, P2 and P3), AT lineage (VI, AG1, and BL), and presumably hybrid populations (FE, CH, CE, and OM) according to previous isozyme, mtDNA and microsatellite data (Bouza et al. 1999, 2001, 2008; Martínez et al. 2007; Vilas et al. 2010). Mediterranean populations belonged to some of tributaries of Ter Basin (Núria and Freser), including two temporal replicates from Núria and Ter rivers (NU04-NU14 and TE04-TE14, respectively; Table 3, Fig. 12). Different restocking incidence with the Central European hatchery stock had been reported previously on the Ter Watershed (Araguas et al. 2017). Two representative samples from the Bagà hatchery stock used for restocking in Catalonia were included in our analysis. Despite there are different hatchery stocks across the Iberian Peninsula they are quite homogeneous due to their recent common origin (García-Marín et al. 1991), and furthermore, even across Europe (Bohling et al. 2016). So, this hatchery was used as reference to detect individuals of hatchery ancestry in wild populations in the samples analysed in this study. Duero Basin is interrupted by big hydroelectric dams built mainly in the 50'-60's of XX century, while the Mediterranean populations mostly by dams lower than three metres from similar construction dates. A total number of 299 individuals captured by electrofishing were analysed, mainly 0+, 1+ and 2+ classes, and included the four native mitochondrial lineages previously identified in the Iberian Peninsula (i.e. AD, AT, DU and ME; Table 3). Furthermore, the samples studied included the two distribution patterns CASANOVA CHICLANA, ADRIÁN 62 known in the studied area: (1) parapatry or spatial segregation with a hybrid contact zone between AT and DU lineages in Miño and Duero basins; and (2) mosaic-sympatry between AD and ME lineages in Mediterranean Drainage. Population genomics for conservation of brown trout resources in the Mediterranean and Atlantic slopes of Iberian Peninsula 63 Table 3. Characteristics of the brown trout (S. trutta) samples from Iberian Peninsula analysed. Origin Code No. individuals mtDNA lineage Miño-Sil Basin 59 Viñao River (2003) VI 16 (AT) Ferreira River (2003) FE 13 (AT/DU) Chamoso River (2003) CH 14 (AT/DU) Lea River (2003) LE 16 (DU) Duero Basin 119 Águeda River (2002) AG1 20 (AT) Porto do Rei Búbal River (2002) BL 19 (AT) Cega River (2002) CE 20 (AT/DU) Omaña River (2002) OM 20 (DU) Pisuerga River 2 (2002) P2 20 (DU) Pisuerga River 3 (2002) P3 20 (DU) Hatchery 39 Hatchery release individuals (2014) BA14 19 (AT) Hatchery spawners (2014) S 20 (AT) Catalonia river basins 82 Núria River (2004) NU04 16 (AD/ME) Núria River (2014) NU14 16 (AD/ME) Queralbs, in Freser River (2014) QB14 18 (AD/ME) Ter River (2004) TE04 14 (AD/ME) Ter River (2014) TE14 18 (AD/ME) mtDNA lineage: Atlantic (AT), Duero (DU), Adriatic (AD) and Mediterranean (ME) CASANOVA CHICLANA, ADRIÁN 70 4.2.5 Genetic differentiation and structure Global and pairwise coefficients of population differentiation (FST values; see Weir and Cockerman 1984) were estimated using different hierarchical criteria. Pairwise FST between populations was calculated using StaMPP R package 1.6 (Pembleton et al. 2013) with the ‘stamppFst’ function. For this, 10,000 bootstrap replicates across loci to generate 95% confidence intervals and p-values regarding the null hypothesis (FST = 0) were used. Global FST for the whole dataset and for each region considered was calculated using R package Genepop 1.1.7 with the ‘Fst’ function. Four different statistical methods were applied to investigate genetic structure in wild populations using all SNPs, including nonparametric approaches (e.g. Discriminant Analysis of Principal Components, DAPCs) approaches, as recommended by Linck and Battey (2019): (1) Analysis of the MOlecular VAriance (AMOVA; Excoffier et al. 1992), (2) bayesian clustering method with STRUCTURE, (3) DAPCs with cross validation and (4) DAPCs with the number of Principal Components (PCs) comprehending ≥ 90% of the variance (e.g. Ríos et al. 2020; Vera et al. 2018). Further STRUCTURE analyses were performed for wild hybrid zones previously identified (Bouza et al. 2001; Martínez et al. 2007; Vilas et al. 2010) with a selection of AIMs (Ancestry Informative Markers) with FST > 0.50 between reference populations. AMOVAs were performed with Arlequin 3.5.2.2 (Excoffier and Lischer 2010), computing F-statistics derived from different hierarchical partitions. Statistical significance of F-statistics for each scenario was tested with 10,000 permutations. With this approach, the best grouping scenario maximizes FCT value (relative component of diversity among groups) by reducing FSC one (idem among populations within groups). Nine a priori grouping scenarios were independently analysed: (1) Atlantic vs Mediterranean slope populations; (2) Population genomics for conservation of brown trout resources in the Mediterranean and Atlantic slopes of Iberian Peninsula 71 populations grouped according to their basin (i.e. Miño, Duero and Ter basins); (3) Miño and Duero basin populations; (4) across temporal replicates and localities in Mediterranean Slope; (5a) Duero Basin populations grouped according to STRUCTURE cluster composition (see Results section); (5b) populations grouped according to Bouza et al. (2001); (5c) populations grouped according to Martínez et al. (2007); (6a) Miño basin populations grouped according to STRUCTURE determined cluster composition and (6b) populations grouped according to Bouza et al. (2008) and Vilas et al. (2010). For STRUCTURE and DAPC analyses, samples within the temporal range 2002-2004 were used. STRUCTURE analyses were run with the same admixture model and the same length in burn-in and MCMC described above. Non a priori population information was used (POPINFO = 0). When only populations from one basin were included, correlated allele frequency models were used, nevertheless, when populations from different basins were included independent allele frequency models were applied due to the high differentiation reported between populations belonging to different basins. StructureSelector web based software (Li and Liu 2018) was used to obtain K estimators and CLUMPAK outputs (Kopelman et al. 2015). Three K estimators were used to identify the most likely number of clusters: the deltaK ad hoc estimator (Evanno 2005), Mean LnP(K) (Pritchard et al. 2000) and MedMeaK (Puechmaille 2016). CLUMPAK output files rendered STRUCTURE bar plots illustrating membership of individuals to inferred genomic clusters. DAPC 2-step approach: A Principal Component Analysis (PCA) from the matrix of the genotypes is performed and then, a selected number of Principal Components (PCs) instead of the original SNP genotypes is used as input for the linear discriminant analysis (LDA). Initially, the ‘find.cluster’ function, implemented in the software Adegenet 2.1 (Jombart and Ahmed 2011), working with all PCs was CASANOVA CHICLANA, ADRIÁN 72 applied to determine the best supported number of genetic clusters using the bayesian information criterion (BIC). The ‘find.cluster’ function runs successive K-means clustering with increasing number of clusters (K) and provides a BIC value for each evaluated K-value, where the lowest BIC is the “optimal” number of clusters. The maximum number of clusters assayed was the double the number of populations and 100,000 iteration per run and 100 starting centroids were used. The selection of the optimal number of PCs to be further used in the LDA was done via cross-validation method where the data are split into: a training set (90% of the data) and a validation set (10% of the data). Cross validation was carried out in two steps (‘xvalDapc’ function): (1) a maximum number of 300 PCs were tested with 100 replicates; (2) with the results obtained, a second cross-validation was run by specifying a narrow range of PCs with 1,000 replicates. The best number of PCs retained was associated with the lowest Root Mean Square Error (RMSE). Additionally, DAPCs retaining at least 90% of the cumulative variation of the data were performed. The resultant clusters were represented in a 2D-scatterplot using the best linear components of DAPC. After this, the software GENECLASS2 (Piry et al. 2004) was used to ascertain the presence of first-generation immigrants (F0) at populations. Since plausibly not all immigrant recipient populations are sampled, the L-home option (Lh test statistic; Paetkau et al. 2004) of GENECLASS2 was used (Lh; the likelihood of drawing that individual’s genotype from the population in which it was samples). The simulation algorithm of Paetkau et al. (2004) with 1,000 simulated individuals was applied. 4.2.6 Effective population size estimation Contemporary effective population size (Ne) was estimated for each population using NeEstimator 2.1 software (Do et al. 2014) with the Linkage Disequilibrium (LDNe) method under a random mating model (Waples 2006). In the two Mediterranean locations where, Population genomics for conservation of brown trout resources in the Mediterranean and Atlantic slopes of Iberian Peninsula 73 temporal replicates were available (Núria and Ter), three different formulations of the temporal method (Jorde and Ryman 2007; Nei and Tajima 1981; Pollak 1983) and generation sets were used. The 95% confidence intervals were determined using the non-parametric jackknife method, more recommendable when the number of loci is large (> 100; see NeEstimator Help file). To prevent potential biases introduced by low frequency alleles, singleton alleles were removed in all analyses to ensure a minimum allele frequency (MAF) > 1/2N. Due to linkage disequilibrium can be generated by the sampling process itself (England et al. 2006) and the small sample size (N~15), when LDNe method was applied, seven minimum allele frequency thresholds were used (between 0.10 and 0.40). Accordingly, the effect of different MAF thresholds over Ne estimates was evaluated as suggested by Marandel et al. (2020). Ne estimates were considered confident when this subsampling SNP method according to MAF reached a plateau of Ne estimates. Finally, the impact of physical linkage among markers was evaluated with Ne estimates from comparisons between SNPs placed in different chromosomes to calculate r2. 4.2.7 Adaptative variation SNPs in genomic regions under selection (outlier dataset) were identified by both Arlequin 3.5.2.2 and BayeScan 2.1. Different subsets of the total dataset were set up by removing SNPs that were monomorphic or singletons in sample subsets (e.g. Miño Basin) in the same way as in the structure analyses. Different statistical approaches were applied to identify a confident set of outlier loci as recommended Narum and Hess (2011). The BayeScan procedure removes the effects shared by all loci (beta), influenced by genetic-drift, from locus-specific effects (alpha), potentially driven by selection. Departure from neutrality at a given locus is assumed when alpha is significantly different from 0. The values of alpha can be informative about the type of selection (i.e. positive values suggest divergent selection and CASANOVA CHICLANA, ADRIÁN 74 negative values suggest balancing selection). For these analyses, very small sample size can be used assuming the lower power of the tests performed, but with no particular risk on bias. The BayeScan analyses were carried out for 20 pilot runs, 5,000 iterations, 100,000 burn-in steps and -pr_odds flag 10 (i.e. the odds for a neutral evolution model were 10 times higher than a model including selection). Loci with qvalues (False Discovery Rate, FDR) < 0.05 were considered significant outliers. In Arlequin, two models for detection of loci under selection were implemented: (1) the finite island model and (2) the hierarchical island model (as defined by Slatkin and Voelm 1991) to avoid a large fraction of false positives when populations share a recent history or belong to a hierarchically subdivided population (Excoffier et al. 2009). For the finite island model, Arlequin was set up for testing 100,000 simulations with 1,000 demes, and when using a hierarchical finite island model 100,000 simulations with 1,000 simulated demes and 100 groups were set up. For detecting outlier loci, FST was used in both models, as recommended by the Arlequin manual. Loci with p-value < 0.01 were considered as significant outliers considering the tendency of this program to false positives (Narum and Hess 2011). For the hierarchical finite island model different grouping with significant FCT in AMOVA analyses were used. Outliers were analysed (1) for all samples belonging to the temporal range 2002-2004; (2) among populations belonging to Atlantic Slope (i.e. Miño and Duero basins); (3) among populations from Duero Basin; (4) among populations from Miño Basin. The obtained outliers were classified into two categories: (1) suggestive outliers, those detected in any of the methods applied; (2) consistent outliers, those detected with all methods. Non-synonymous substitutions due to allelic variants at SNPs located in exons were evaluated. Open reading frames (ORFs) were checked with ORFfinder (NCBI) and BioEdit 7.2.5 (Hall 1999) using coding DNA sequences Population genomics for conservation of brown trout resources in the Mediterranean and Atlantic slopes of Iberian Peninsula 75 (CDSs) from the brown trout database in Ensembl (https://www.ensembl.org/Salmo_trutta/). Gene mining to identify candidate genes under selection was performed on those genomic regions where two or more consistent outlier loci were placed in the same chromosome region within a range < 500 kb. Gene annotation was performed with the Ensembl database 103. GO terms were obtained with Blast2GO (Götz et al. 2008), module of the software OmicsBox 1.4.12 (https://www.biobam.com/omicsbox-update-1-4/). 4.3 RESULTS The number of potential RAD-loci detected in silico in the brown trout genome was 557,761, representing ~0.85% of the total genome assembly size. Among them, five RAD-loci were in the mitochondrial genome. The number of 2bRAD-loci per chromosome was highly correlated with chromosome length (R2 = 0.98, p-value = 0.000). On average one 2bRAD-locus occurred every 4,000 nucleotides. A total of 2,054,115,335 raw reads from 299 individuals were produced on the NextSeq500 sequencing platform, averaging 6,869,950 reads per individual. After filtering steps 1,410,823,259 reads were retained (68.7%). The most astringent filtering step was RE recognition site presence (~20% removed reads). Out of those reads, 610,183,196 aligned against the brown trout reference genome with -v 3 (43.2%). Sixteen individuals with less than 0.9 M aligned reads were discarded (Table 4). CASANOVA CHICLANA, ADRIÁN 76 Table 4. Number of individuals per sample used for subsequent analyses (N = 283), classified into natural basins or hatchery. Origin Code No. individuals No. final individuals Miño-Sil Basin 59 56 Viñao River VI 16 15 Ferreira River FE 13 13 Chamoso River CH 14 13 Lea River LE 16 15 Duero Basin 119 106 Águeda River AG1 20 16 Porto do Rei Búbal River BL 19 16 Cega River CE 20 19 Omaña River OM 20 20 Pisuerga River 2 P2 20 18 Pisuerga River 3 P3 20 17 Hatchery 39 39 Hatchery release individual BA14 19 19 Hatchery spawners R 20 20 Catalonia river basins 82 82 Núria River NU04 16 16 Núria River NU14 16 16 Queralbs, in Freser River QB14 18 18 Ter River TE04 14 14 Ter River TE14 18 18 Population genomics for conservation of brown trout resources in the Mediterranean and Atlantic slopes of Iberian Peninsula 77 A total of 361,129 RAD-loci were built by Stacks 2.41, comprising 606,037,805 aligned reads and representing an average 9.0x coverage per locus and individual. After Stacks 2.41 parsing, 191,406 SNPs were retained. Among these, 1,150 SNPs (0.6%) were within overlapping RAD-loci, being discarded for further analyses. Among the next filtering steps, the main SNP dropping was due to population representation ≥ 60% (66.7%; Table 5) and then the MAC ≥ 3 filtering step (41.1%). Table 5. Filtering steps for the final SNP panel. Filtering step Number of SNPs Stacks 2.41 output 191,406 (6 mitochondrial SNPs) MinDP ≥ 6x 191,406 MAC ≥ 3 112,652 SNPs represented 60%/pop 37,552 HW (p-value > 0.05) 34,583 No overlapping SNPs 34,548 RAD-loci with ≤ 3 biallelic SNPs 32,981 1 SNP/RAD-locus selected (more pol.) 25,247 (3 mitochondrial SNPs) Well-located SNPs 24,830 Three mitochondrial SNPs were conserved after filtering steps, all of them located within coding regions. None of the final RAD-loci, either nuclearor mtDNA-linked, had uncalled nucleotides in the reference genome to which they were aligned (i.e. hard-masked DNA sequences). Thus, the final dataset was composed by 24,830 nuclear SNPs (Table 5). CASANOVA CHICLANA, ADRIÁN 78 4.3.1 Identification of individuals of farmed ancestry A total of 34 individuals of hatchery ancestry were identified with the outlined criteria in the wild populations studied using the whole genomic information. All individuals from VI, CH, AG1, OM, P2 and P3, TE04 and TE14 were considered wild, while, at the other end, most individuals from QB14 were of hatchery ancestry (only four wild individuals). Therefore, this population was removed for further characterization of Iberian populations. One individual of hatchery ancestry from FE, three from LE, five from CE, five from NU04, four from NU14 and two from BL (as reported by Martínez et al. 2007), were identified according to the qW criterion. However, the populations with hatchery ancestry were in Hardy-Weinberg Equilibrium (HWE; see next section). All the ~ 25,000 SNPs were ranked by FST between hatchery (BA14 and S, samples for release and spawners, respectively) vs all wild to identify the SNPs with highest diagnostic power to develop a costeffective and high-resolution tool to elucidate hatchery ancestry. The number of SNPs with FST > 0.95, FST > 0.99 and FST = 1 were 214, 38 and nine, respectively (Table 6). Then, the performance of these different SNP subsets obtained were evaluated as the percentage of correct classification regarding the whole SNP dataset (Fig. 13). In ten populations the three subsets showed the same classification success as the whole SNP dataset. Despite the subset with the highest number of SNPs (214 SNPs panel) showed the best performance, the 38 SNPs panel performed very similarly and only one individual from NU14, excluding CE, was not correctly classified. Furthermore, with only nine diagnostic SNPs the classification success was encouraging. In all evaluations, CE was remarkable because the classification success dropped to ~80% regarding the whole SNP dataset (Fig. 13). Population genomics for conservation of brown trout resources in the Mediterranean and Atlantic slopes of Iberian Peninsula 79 Table 6. Ancestry informative markers (AIMs) selected. W: Wild populations, QB14 population was not employed; H: Hatchery samples, BA14 and S. FST W-H Number of SNPs Number of chromosomes represented >0.95 214 37 >0.99 38 18 1 9 5 Figure 13. Correct individual hatchery ancestry classification in all populations studied using the whole SNP dataset. The 80% of correct classification is highlighted with a red line. The lower performance in CE of most SNP subsets has likely to do with the qW values of those misclassified individuals, very close to the threshold established, and the same occurred with NU14-01 individual (Table 7). The qH obtained among the different SNP subsets showed a highly significant correlation between each SNP subset and the whole SNP dataset, especially when using the panel of 38 and 214 SNPs (r = 0.88 and 0.92, respectively; p-value < 0.001). The incidence of restocking measured as the mean hatchery ancestry (qH) was between low and moderate (range: FE 0.01 - NU04 0.11; Table 8) excluding QB14, the most affected population (qH = 0.20). The Atlantic Slope CASANOVA CHICLANA, ADRIÁN 86 Figure 14. (Cont.) Population genomics for conservation of brown trout resources in the Mediterranean and Atlantic slopes of Iberian Peninsula 87 4.3.2 Genetic diversity Only wild individuals were used for estimation of genetic diversity and population structure. Genetic diversity estimators were rather heterogeneous across populations and river basins (Table 9). Allelic richness (Ar) ranged from 1.067 (P2) to 1.38 (FE). Expected heterozygosity ranged between Pisuerga river populations (P2: 0.023 and P3: 0.031) and two Miño basin populations (CH: 0.121 and FE: 0.123). Within the Atlantic Slope expected heterozygosity was more than double in Miño than in Duero Basin (average range: 0.106 vs 0.045), while the Ter Basin in the Mediterranean Slope showed the lowest value (0.040), being slightly lower than Duero Basin. Between temporal replicates the highest difference was found in Núria River (He 0.035 vs 0.044 for NU04 and NU14, respectively). Observed and expected heterozygosities were very similar in each population as reflected by the intrapopulation fixation index (FIS), which was low but negative for almost all populations. These values were significant in most populations using the 95% confidence interval approach (FIS ≠ 0; value not included in the confidence interval), despite all of them did not globally deviated from random mating according to exact tests (see below Hardy-Weinberg equilibrium). Both hatchery stock samples showed high genetic diversity figures, very close to the upper range detected in Miño Basin. A notable amount of private alleles were detected, in all cases at very low frequencies, reflecting the structure of brown trout populations. CASANOVA CHICLANA, ADRIÁN 88 Table 9. Genetic diversity in brown trout populations from Iberian Peninsula and a hatchery stock. Ar, mean allelic richness; Ho, mean observed heterozygosity; He, mean expected heterozygosity; FIS, global intrapopulation fixation index; BC 95% CI, bias corrected 95% confidence intervals of FIS values; PA, number of private alleles. No hatchery ancestry individuals were used, with the exception of hatchery individuals from Bagà fish farm (i.e. BA14 and S). A total of 245 individuals were used. The values were rounded to three decimal numbers. Population codes Ar Ho He FIS lower and upper BC 95% CI PA VI 1.239 0.088 0.080 -0.100 -0.155 -0.065 340 FE 1.380 0.134 0.123 -0.085 -0.148 -0.047 119 CH 1.358 0.129 0.121 -0.065 -0.124 -0.031 170 LE 1.321 0.111 0.101 -0.102 -0.190 -0.050 68 AG1 1.135 0.045 0.045 0.001 -0.086 0.104 457 BL 1.276 0.097 0.091 -0.069 -0.114 -0.043 471 CE 1.203 0.047 0.043 -0.093 -0.158 -0.050 179 OM 1.112 0.038 0.035 -0.075 -0.102 -0.053 316 P2 1.067 0.024 0.023 -0.050 -0.085 -0.023 86 P3 1.090 0.035 0.031 -0.143 -0.232 -0.089 156 BA14 1.275 0.09 0.085 -0.059 -0.089 -0.038 197 S 1.301 0.098 0.090 -0.084 -0.144 -0.040 312 NU04 1.105 0.037 0.035 -0.067 -0.185 0.022 35 NU14 1.152 0.047 0.044 -0.067 -0.134 -0.021 37 TE4 1.120 0.043 0.041 -0.041 -0.091 -0.004 10 TE14 1.114 0.041 0.039 -0.061 -0.107 -0.018 17 Population genomics for conservation of brown trout resources in the Mediterranean and Atlantic slopes of Iberian Peninsula 89 No deviations from Hardy-Weinberg equilibrium (HWE) were detected at population level (global Fisher's test; Table 10). The proportion of loci that showed HWE deviation at p-value < 0.05 ranged from 0.9% in NU04 to 10.7% in AG1 with an average of 3.0%, below to expected 5% expected by chance. The proportion of deviations increased when the individuals with hatchery ancestry were included in the analysis, as expected due to a Wahlund effect. CASANOVA CHICLANA, ADRIÁN 90 Table 10. Deviation from Hardy–Weinberg (HW) expectations for each population over all loci using Fisher's exact tests. In grey when individuals with hatchery ancestry were included in analysis. Population codes No. individuals N SNPs per analysis HW (global p-value) HW (% p-values <0.05) VI 15 5,260 1.00 2.1 FE 12 8,006 1.00 1.5 FEH 13 8,314 1.00 1.8 CH 13 8,033 1.00 1.6 LE 12 6,166 1.00 1.1 LEH 15 7,193 1.00 1.9 AG1 16 3,145 1.00 10.7 BL 14 6,192 1.00 1.8 BLH 15 6,693 1.00 4.3 CE 14 3,448 1.00 1.1 CEH 19 6,510 1.00 2.7 OM 20 2,479 1.00 3.0 P2 18 1,609 1.00 3.5 P3 17 2,119 1.00 5.5 BA14 19 6,923 1.00 3.0 S 20 7,686 1.00 1.9 NU04 11 2,247 1.00 1.3 NU04 H 16 7,679 1.00 0.9 NU14 12 2,751 1.00 2.3 NU14 H 16 6,374 1.00 10.7 QB14H 18 9,690 1.00 1.8 TE04 14 2,790 1.00 2.0 TE14 18 2,754 1.00 2.9 Population genomics for conservation of brown trout resources in the Mediterranean and Atlantic slopes of Iberian Peninsula 91 4.3.3 Genetic differentiation and structure All SNPs were included in these analyses. All pairwise FST comparisons were significant (Table 11), no confidence intervals included zero. Using wild and hatchery populations, the lowest pairwise FST was found between BA14 and S (hatchery) and between temporal replicates from Ter (FST = 0.015 and FST = 0.017, for BA14-S and TE04TE14 respectively). Higher differentiation values were obtained in Núria temporal replicates. Between Núria replicates, FST value was far higher when exclusively wild samples were used (FST from 0.030 to 0.102 when using all and wild individuals, respectively). High pairwise differentiation between populations belonging to the same river basin was observed in Duero Basin (FST = 0.607 between AG1-P2), including populations that belonged to the same tributary (i.e. Pisuerga River, FST = 0.294 between P2-P3). Putative hybrid populations from Miño Basin showed lower FST values with DU lineage population (i.e. LE) than with the AT representative population (VI). The global FST was 0.680 (pvalue = 0.000). Pairwise basin comparisons between Mediterranean (Ter) and the Atlantic Slope (Miño and Duero) were higher (FST values > 0.660; Table 12) than within Atlantic Slope (FST = 0.300). Global FST was higher in Duero Basin than in Miño Basin (0.460 vs 0.192, respectively). CASANOVA CHICLANA, ADRIÁN 92 Table 11. Pairwise FST values between populations (below diagonal) obtained with N = 245 individuals from wild and hatchery populations and the whole SNP panel with 24,830 SNPs. Miño Basin Duero Basin Hatchery Mediterranean Basin VI FE CH LE AG1 BL CE OM P2 P3 BA14 S NU04 NU14 TE04 TE14 VI FE 0.265 CH 0.267 0.050 LE 0.326 0.071 0.106 AG1 0.584 0.419 0.445 0.425 BL 0.403 0.218 0.243 0.222 0.401 CE 0.575 0.402 0.429 0.407 0.536 0.431 OM 0.627 0.466 0.494 0.458 0.538 0.471 0.372 P2 0.671 0.515 0.539 0.522 0.607 0.523 0.428 0.390 P3 0.636 0.474 0.501 0.470 0.556 0.481 0.340 0.288 0.294 BA14 0.567 0.473 0.475 0.552 0.703 0.584 0.695 0.739 0.762 0.743 S 0.557 0.466 0.468 0.543 0.692 0.576 0.684 0.728 0.750 0.731 0.015 NU04 0.784 0.708 0.709 0.761 0.856 0.764 0.860 0.876 0.903 0.886 0.735 0.721 NU14 0.768 0.693 0.694 0.744 0.840 0.750 0.842 0.862 0.887 0.871 0.719 0.706 0.102 TE04 0.780 0.711 0.711 0.759 0.846 0.762 0.848 0.866 0.889 0.874 0.732 0.719 0.194 0.156 TE14 0.793 0.732 0.732 0.777 0.852 0.778 0.854 0.869 0.890 0.876 0.747 0.734 0.223 0.188 0.017 Population genomics for conservation of brown trout resources in the Mediterranean and Atlantic slopes of Iberian Peninsula 93 Table 12. Pairwise FST values between river drainages (below diagonal). FST values were obtained with N = 176 individuals from wild contemporary populations and the 20,293 SNPs panel. Atlantic Slope: Miño Basin populations: VI, FE, CH, LE; Duero Basin: AG1, BL, CE, OM, P2, P3; Mediterranean Slope: Ter Basin: NU04, TE04. The global FST values for each basin are shown on the diagonal. MIÑO DUERO MED MIÑO 0.192 DUERO 0.300 0.460 MED 0.661 0.771 0.151 AMOVA results revealed significant genetic structuring among river basins (Table 13) belonging to Atlantic vs Mediterranean slopes or between different river basins. However, only marginal significant genetic structuring was found in Duero Basin according to two of the criteria established and no intergroup (FCT) significant values were detected in Miño Basin in the two scenarios tested. In the Mediterranean Slope, the geographical variance (locations) within Ter Basin was higher than temporal variance (temporal samples). In all hypotheses tested the highest variance was found within populations. CASANOVA CHICLANA, ADRIÁN 94 Table 13. Analyses of AMOVA with Arlequin 3.5.2.2. In light green, grouping hypotheses with populations coming from different river basins. Atlantic Slope with populations from Miño and Duero basins and Mediterranean Slope with two samples (NU04 and TE04) from Ter Basin. In dark yellow, population structure between localities and temporal replicates in Mediterranean Slope. In light yellow, hypotheses tested according to STRUCTURE results and previous reports using populations from Duero Basin (Bouza et al. 2001; Martínez et al. 2007). Atlantic mtDNA lineage (AT): AG1 and BL populations. Putative hybrid zone (PHZ: CE and OM populations). Duero mtDNA lineage (DU): P2 and P3 populations. In light blue, hypothesis tested using populations from Miño Basin according to STRUCTURE results and previous references (Bouza et al. 2008 and Vilas et al. 2010). Hypotheses df Variance % Variation 1. Two groups (Atlantic vs Mediterranean Slope) Among groups 1 1994.74 67.35 Among populations within groups 10 420.40 14.19 Within populations 340 546.79 18.46 F statistics: FCT = 0.67**, FSC = 0.43***, FST = 0.82*** 2. Three groups (Miño, Duero and Ter basins) Among groups 2 1064.83 56.26 Among populations within groups 9 281.162 14.85 Within populations 340 546.79 28.89 F statistics: FCT = 0.56***, FSC = 0.34***, FST = 0.71*** 3. Two groups (Miño vs Duero basins) Among groups 1 298.05 25.69 Among populations within groups 8 298.91 25.77 Within populations 292 563.11 48.54 F statistics: FCT = 0.26**, FSC = 0.35***, FST = 0.51*** +p-value < 0.1, **p-value < 0.05, ***p-value < 0.01. Population genomics for conservation of brown trout resources in the Mediterranean and Atlantic slopes of Iberian Peninsula 95 Table 13. (Cont.) Hypotheses df Variance % Variation 4. Two groups within MED Slope (NU04+NU14 vs TE04+TE14) Among groups (rivers) 1 80.46 14.65 Among temporal replicates (2004-14) 2 24.89 4.53 Within populations 106 443.99 80.82 F statistics: FCT = 0.15, FSC = 0.05***, FST = 0.19*** 5a. Two groups within Duero Basin (AT vs PHZ+DU) Among groups 1 248.54 27.52 Among populations within groups 4 233.32 25.84 Within populations 192 421.15 46.64 F statistics: FCT = 0.28+, FSC = 0.36***, FST = 0.53*** 5b. Three groups within Duero Basin (AT, PHZ, DU) Among groups 2 142.72 17.64 Among populations within groups 3 245.40 30.32 Within populations 192 421.15 52.04 F statistics: FCT = 0.176+, FSC = 0.37***, FST = 0.48*** 5c. Four groups within Duero Basin (AT, PHZ, P2, P3) Among groups 3 49.75 6.32 Among populations within groups 2 316.13 40.17 Within populations 192 421.152 53.51 F statistics: FCT = 0.06, FSC = 0.43***, FST = 0.46*** +p-value < 0.1, **p-value < 0.05, ***p-value < 0.01. CASANOVA CHICLANA, ADRIÁN 102 Figure 20. Scatterplots of individuals on the two principal DA eigenvalues of DAPC. All populations belong to Duero Basin. The graph represents the individuals as dots and the groups as inertia ellipses. PCAs and DAs eigenvalues are displayed inset. A: Scatterplot with 60.6% of the variance. B: Scatterplot with 91.0% of the variance. Wild hybrid zones (Duero and Miño basins) Within the Duero Basin, the most likely number of populations units varied among K estimators: ΔK, MedMean K, and Mean LnP(K), two, five, and six, respectively. With K = 2, BL appeared as a distinct mixed population, not so CE (Fig. 21). Once more, the results obtained do not show OM as a hybrid population. With the highest K each population was represented as a separate cluster. Population genomics for conservation of brown trout resources in the Mediterranean and Atlantic slopes of Iberian Peninsula 103 Figure 21. CLUMPAK plot of STRUCTURE assignment from K = 2 to 6. N = 99 wild individuals from six populations from Duero Basin and 1,778 AIMs were used. mc: minor cluster. Each individual is represented as a vertical bar partitioned into segments according to the proportion of the genome belonging to each of the clusters identified (K) by STRUCTURE. Within the Miño Basin, the most likely number of populations units was similar among K estimators: ΔK with two and MedMean K, and Mean LnP(K) with three. With K = 2, two clusters were delineated corresponding to VI in the outlet and the inner Miño Basin populations, with some degree of hybridization observed in FE and CH (qVI ~ 0.10; Fig. 22). With K = 3, FE and CH appeared to be hybrid populations but one of the components would not belong to any of the reference populations used. CASANOVA CHICLANA, ADRIÁN 104 Figure 22. CLUMPAK plot of STRUCTURE assignment from K = 2 to 3. N = 52 wild individuals from four populations from Miño Basin and 1,429 AIMs were used. Each individual is represented as a vertical bar partitioned into segments according to the proportion of the genome belonging to each of the clusters identified (K) by STRUCTURE. 4.3.4 Effective population size The estimates of effective population size (Ne) using NeEstimator software yielded quite often finite values (Tables 14 and 15). A plateau in Ne estimates was obtained with some populations (e.g P2 and P3). Nevertheless, the upper boundaries of the 95% confidence intervals were infinite with many populations and the different MAF thresholds applied, indicating that the estimates may not be very robust in these cases. Ne estimates of the same order of magnitude were obtained between LDNe and temporal methods when temporal replicates were available, showing that Ter population would have a higher effective population size than Núria population (Tables 14 and 16), the last one with higher hatchery ancestry proportion. Ne estimates using pairs of SNPs placed in distinct chromosomes showed higher values (Table 15) than obtained considering all the markers available, which makes sense since LD is reduced when using markers on different chromosomes. In any case, the values obtained between both methodologies were similar. Population genomics for conservation of brown trout resources in the Mediterranean and Atlantic slopes of Iberian Peninsula 105 Table 14. Effective population size calculated with the Linkage Disequilibrium method (Waples 2006) using all SNP pairs and considering different MAF thresholds. In parentheses, 95% CI. Inf.: Infinite. ≥0.10 ≥0.15 ≥0.20 ≥0.25 ≥0.30 ≥0.35 ≥0.40 VI 82.1 (4.7-Inf.) 83.0 (4.1-Inf.) 96.2 (3.8-Inf.) 116.1 (4.0-Inf.) 197.7 (4.5-Inf.) 559.6 (4.6-Inf.) Inf. (5.5-Inf) FE 35.1 (9.6-Inf.) 37.7 (9.8-Inf.) 40.2 (9.9-Inf.) 44.0 (10.0-Inf.) 48.7 (9.8-Inf.) 57.9 (10.3-Inf.) 73.8 (10.7-Inf.) CH Inf. (6.7-Inf.) 7032.0 (6.0-Inf.) 5607.3 (5.4-Inf.) 160297.9 (5.3-Inf.) Inf. (5.2-Inf.) Inf. (4.9-Inf.) Inf. (5.5-Inf.) LE Inf. (1.9-Inf.) Inf. (1.9-Inf.) Inf. (1.8-Inf.) Inf. (1.7-Inf.) Inf. (1.7-Inf.) Inf. (1.7-Inf.) Inf. (1.9-Inf.) AG1 19.0 (6.7-Inf.) 21.2 (7.6-Inf.) 21.7 (7.7-Inf.) 22.1 (7.7-Inf.) 23.7 (8.1-Inf.) 24.0 (8.1-Inf.) 28.1 (8.4-Inf.) BL 496.2 (24.0-Inf.) 562.2 (22.2-Inf.) 484.4 (21.9-Inf.) 443.4 (21.6-Inf.) 582.9 (21.2-Inf.) 1301.1 (22.7-Inf.) Inf. (21.3-Inf.) CE 3.9 (1.8-20.1) 9.1 (2.6-105.8) 18.8 (4.9-Inf.) 33.5 (6.7-Inf.) 117.4 (9.8-Inf.) Inf. (12.1-Inf.) Inf. (15.0-Inf.) OM 58.1 (23.1-Inf.) 56.0 (22.0-Inf.) 57.5 (22.4-Inf.) 57.2 (22.3-Inf.) 58.7 (22.6-Inf.) 60.0 (23.0-Inf.) 62.9 (23.4-Inf.) P2 33.0 (18.0-105.9) 33.5 (17.8-119.1) 32.5 (17.8-101.3) 34.9 (18.8-120.3) 37.9 (19.5-167.2) 41.1 (18.8-743.2) 46.6 (19.3-Inf.) P3 10.4 (6.4-17.8) 11.3 (7.3-18.6) 11.9 (7.7-19.8) 11.7 (7.3-20.2) 11.9 (7.3-21.2) 11.6 (6.9-21.9) 10.1 (6.0-18.2) NU04 4.2 (2.1-13.2) 3.3 (2.0-11.9) 3.1 (1.9-10.5) 3.0 (1.8-9.8) 2.8 (1.7-10.4) 2.9 (1.7-10.8) 4.6 (2.1-15.0) NU14 27.0 (8.3-Inf.) 35.0 (10.1-Inf.) 36.9 (10.3-Inf.) 35.7 (10.1-Inf.) 34.2 (9.2-Inf.) 34.7 (9.1-Inf.) 43.4 (10.5-Inf.) TE04 59.1 (26.2-Inf.) 55.1 (23.1-Inf.) 52.9 (23.6-Inf.) 51.3 (23.4-Inf.) 45.0 (21.6-507.8) 46.1 (22.7-383.9) 39.1 (19.7-226.3) TE14 24.1 (14.0-56.2) 22.9 (13.5-51.5) 21.4 (12.6-46.8) 22.0 (13.2-46.6) 20.5 (12.2-43.7) 19.4 (11.5-41.1) 18.0 (10.1-42.3) CASANOVA CHICLANA, ADRIÁN 106 Table 15. Effective population size calculated with the Linkage Disequilibrium method (Waples 2006) using pairwise comparisons between SNPs located in different chromosomes and considering different MAF thresholds. In parentheses, 95% CI. Inf.: Infinite. ≥0.10 ≥0.15 ≥0.20 ≥0.25 ≥0.30 ≥0.35 ≥0.40 VI 99.7 (4.8-Inf.) 101.4 (4.3-Inf.) 121.6 (4.0-Inf.) 158.2 (4.2-Inf.) 380.9 (4.6-Inf.) Inf. (4.8-Inf.) Inf. (5.6-Inf.) FE 37.0 (9.8-Inf.) 39.9 (10.0-Inf.) 42.9 (10.1-Inf.) 47.6 (10.2-Inf.) 53.7 (10.0-Inf.) 66.3 (10.5-Inf.) 88.9 (11.0-Inf.) CH Inf. (6.8-Inf.) Inf. (6.1-Inf.) Inf. (5.5-Inf.) Inf. (5.4-Inf.) Inf. (5.2-Inf.) Inf. (4.9-Inf.) Inf. (5.5-Inf.) LE Inf. (1.9-Inf.) Inf. (1.9-Inf.) Inf. (1.8-Inf.) Inf. (1.7-Inf.) Inf. (1.7-Inf.) Inf. (1.7-Inf.) Inf. (1.9-Inf.) AG1 20.3 (6.9-Inf.) 23.1 (7.9-Inf.) 23.8 (8.1-Inf.) 24.4 (8.0-Inf.) 26.8 (8.5-Inf.) 28.0 (8.7-Inf.) 34.5 (9.1-Inf.) BL 4592.4 (25.3-Inf.) Inf. (23.3-Inf.) 3874.7 (22.9-Inf.) 3665.2 (22.7-Inf.) Inf. (22.3-Inf.) Inf. (24.0-Inf.) Inf. (22.5-Inf.) CE 6.6 (2.157.6) 12.6 (2.9-Inf.) 27.2 (6.0-Inf.) 59.4 (7.8-Inf.) Inf. (11.5-Inf.) Inf. (14.4-Inf.) Inf. (18.6-Inf.) OM 64.5 (24.3-Inf.) 63.1 (23.2-Inf.) 66.3 (23.9-Inf.) 66.6 (23.8-Inf.) 68.9 (24.2-Inf.) 72.3 (24.8-Inf.) 77.8 (25.5-Inf.) P2 38.3 (19.5182.9) 39.6 (19.5243.6) 39.3 (19.8205.9) 43.3 (21.1329.5) 48.4 (22.31522.4) 54.4 (21.6-Inf.) 68.4 (23.1-Inf.) P3 11.2 (7.0-19.1) 12.3 (8.0-20.4) 13.0 (8.4-21.9) 12.7 (8.0-22.4) 12.8 (7.9-23.5) 12.6 (7.3-24.6) 10.9 (6.420.3) NU04 4.4 (2.1-13.9) 3.3 (2.0-12.6) 3.2 (2.011.1) 3.1 (1.9-10.5) 2.9 (1.7-11.2) 3.0 (1.7-11.6) 5.2 (2.2-16.7) NU14 29.9 (8.7-Inf.) 39.7 (10.6-Inf.) 42.9 (11.0-Inf.) 41.9 (10.8-Inf.) 41.3 (9.8-Inf.) 43.4 (9.8-Inf.) 56.9 (11.0-Inf.) TE04 65.2 (27.5-Inf.) 61.3 (24.3-Inf.) 59.7 (25.0-Inf.) 58.4 (24.9-Inf.) 50.4 (23.0-Inf.) 52.3 (24.1-13273.5) 43.1 (20.6490.4) TE14 26.0 (14.865.3) 25.0 (14.460.4) 23.5 (13.655.5) 24.5 (14.356.7) 22.9 (13.353.1) 21.8 (12.650.3) 20.6 (11.254.9) Population genomics for conservation of brown trout resources in the Mediterranean and Atlantic slopes of Iberian Peninsula 107 Table 16. Effective population size calculated with the temporal method. In parentheses, 95% CI (loci jackknife method). With ten years between samples the most plausible number of generations in brown trout would be three (in bold). Núria (2004-2014) Pollak Nei/Tajima Jorde/Ryman 1 generation 2.5 (2.3-2.6) 3.1 (2.9-3.3) 2.3 (2.1-2.4) 2 generations 4.9 (4.6-5.2) 6.2 (5.8-6.5) 4.5 (4.3-4.8) 3 generations 7.4 (6.9-7.8) 9.2 (8.7-9.8) 6.8 (6.4-7.2) Ter (2004-2014) 1 generation 15.1 (13.1-17.6) 18.8 (16.1-22.4) 14.8 (12.7-17.7) 2 generations 30.3 (26.2-35.3) 37.7 (32.2-44.8) 29.6 (25.5-35.3) 3 generations 45.4 (39.352.9) 56.5 (48.367.2) 44.4 (38.2-53.0) CASANOVA CHICLANA, ADRIÁN 108 4.3.5 Adaptative variation In all cases studied, a higher number of outlier loci were detected with Arlequin than with BayeScan (Table 17). Consistent outliers were found with all analyses (i.e. outliers detected with all statistical approaches). When using samples from both slopes, a big difference was obtained in the number of outliers identified between BayeScan and Arlequin approaches. Further, the proportion of outliers under divergent and balancing selection showed opposite proportions between both approaches, being balancing outliers more frequent with Bayescan (86.4%) while divergent with Arlequin (92.9%), this pattern being even more accentuated in the hierarchical scenario. Some SNPs were shared between approaches for within each type of selection, but never between different types of selection. Balancing selection would be a priori more easily detected using FST outlier tests in a highly divergent genetic scenario, such as that occuring between Atlantic and Mediterranean slopes. Most consistent outliers detected here showed balancing selection (22 of 26 consistent outliers). It is likely that a large part of the divergent outliers detected with Arlequin may be false positives. Three SNP pairs related to divergent selection were found within a 500 kb window in the brown trout genome, and gene mining was performed to identify candidate genes and their associated GO terms (Table 18). Population genomics for conservation of brown trout resources in the Mediterranean and Atlantic slopes of Iberian Peninsula 109 Table 17. Outlier loci detected using different subsets of samples and models. In bold type, the number of outlier loci under divergent selection, in italic type number of outlier loci under balancing selection. MI: Miño Basin; DU: Duero Basin; MED (Ter Basin): Mediterranean Basin. With the first subset of samples two Arlequin hierarchical analyses were performed: (A) hierarchical groups according to slopes (Atlantic and Mediterranean slopes) and (B) hierarchical groups according to river basins (i.e. Miño, Duero and Ter). ID Subset of samples SNPs BayeScan (q-value < 0.05) Arlequin (p-value < 0.01) Arlequin hierarchical (p-value < 0.01) Consistent outliers (detected in all analyses) 1 MI+DU+MED (2004) 20,293 184 (25;159) 3,114 (2,267;847) (A) 2,466 (2,433;33) (B) 2,256 (2,185;71) 26 (4;22) 2 Atlantic Slope 15,158 217 (96;121) 461 (379;82) 311 (256;55) 84 (44;40) 3 Duero Basin 10,777 95 (47;48) 249 (202;47) — 67 (38;29) 4 Miño Basin 11,332 13 (13;0) 608 (263;345) — 4 (4;0) CASANOVA CHICLANA, ADRIÁN 110 Table 18. Gene mining in the brown trout genome on the three windows (500 kb), where pairs of outlier loci were detected when comparing all samples of Mediterranean and Atlantic slopes: (1) chromosome 4 between 32,491,083 and 32,991,083 bp (light yellow); (2) chromosome 21 between 1,125,412 and 1,625,412 bp (light blue); and (3) chromosome 22 between 2,163,580 and 2,663,580 bp (light green). ID GENE DESCRIPTION GO TERMS FOR BIOLOGICAL PROCESS (LEVEL 3) pde4ba cAMP-specific 3',5'-cyclic phosphodiesterase 4B-like GO:0006807 nitrogen compound metabolic process GO:0071704 organic substance metabolic process GO:0050794 regulation of cellular process sgip1a SH3-containing GRB2-like protein 3interacting protein 1 GO:0071840 cellular component organization or biogenesis GO:0051234 establishment of localization ENSSTUG00000015164 relaxin-3-like GO:0050794 regulation of cellular process GO:0051716 cellular response to stimulus GO:0007154 cell communication ENSSTUG00000015172 Possibly ORF2p gene in Danio rerio ttc4 tetratricopeptide repeat protein 4-like GO:0007275 multicellular organism development GO:0048856 anatomical structure development GO:0009653 anatomical structure morphogenesis lurap1 leucine rich adaptor protein 1-like GO:0007275 multicellular organism development GO:0007165 signal transduction GO:0051716 cellular response to stimulus ENSSTUG00000015254 Possibly col24a1 gene in Danio rerio GO terms level 3 for Biological Process category were obtained with Blast2GO. Population genomics for conservation of brown trout resources in the Mediterranean and Atlantic slopes of Iberian Peninsula 111 Table 18. (Cont.) ID GENE DESCRIPTION GO TERMS FOR BIOLOGICAL PROCESS (LEVEL 3) ankhb ANKH inorganic pyrophosphate transport regulator GO:0055085 transmembrane transport GO:0051234 establishment of localization RNF170 ring finger protein 170 GO:0044238 primary metabolic process GO:0006807 nitrogen compound metabolic process GO:0071704 organic substance metabolic process lrrc6 leucine rich repeat containing 6 GO:0006928 movement of cell or subcellular component GO:0007275 multicellular organism development GO:0009653 anatomical structure morphogenesis kcnq3 potassium voltage-gated channel subfamily Q member 3 GO:0055085 transmembrane transport GO:0050794 regulation of cellular process GO:0032879 regulation of localization efr3a EFR3 homolog A GO:0051641 cellular localization GO:0033036 macromolecule localization adgrl2a adhesion G protein-coupled receptor L2a GO:0007165 signal transduction GO:0007154 cell communication GO:0007275 multicellular organism development GO terms level 3 for Biological Process category were obtained with Blast2GO. CASANOVA CHICLANA, ADRIÁN 118 due to the short length of 2b-RAD reads before going to the buildingloci pipeline. Some population parameters between both approaches showed relevant differences (see FIS values in brown trout; Table S5). A practical approach to decide between pipelines with different building-loci strategies for handling genotyping by sequencing (GBS) data is to assay trials with a small subset of data and check for their results using a meaningful set of population parameters, previously selected according to the objectives of the study. Indeed, using the number of SNPs obtained as the main criterion (Puritz et al. 2014) to decide the best building-loci pipeline to be used is not advisable, since a higher number of SNPs does not necessarily indicate a better stacking and confident RAD-seq data (Díaz-Arce et al. 2019), and consequently, it might have a negative impact on the confidence of results and biological inferences. The initial number of SNPs obtained with STA and ALT pipelines across the different species tested was rather similar except for the brown trout and the small-spotted catshark. These species showed the lowest median coverage, near to the selected threshold coverage filter (8x) hence the differences on the number of putative loci from input data. While Stacks 2, with a de novo approach, starts with individual data demanding a number of identical reads to build a locus (see Material and Methods), Meyer’s 2b-RAD v2.1 works with a combined subset of confident reads from all samples to build a global reference panel to which align every read. In cases with low coverage, a less demanding criterion to build loci can produce large differences in the initial number of SNPs. Nevertheless, through the filtering steps, the SNP number from building-loci pipelines converged in both species and importantly, the population parameters between SNPs panels were similar. Once chosen the pipeline, it would be recommendable to run several trials with different parameters to properly adjust them to the dataset. For instance, -M in Stacks (which defines the maximum nucleotide differences allowed between intraindividual putative loci) Discussion 119 depends on the levels of polymorphism of the species and -m in Stacks on the existing coverage (Paris et al. 2017). In the same way as for the choice of the building-loci pipeline, it would be advisable to choose parameters taking into account the results from population outcomes, since there is not a unique pipeline suited to every situation, as already indicated Torkamaneh et al. (2016). After the building-loci pipeline, it is important to adjust filtering (criteria and order; O’Leary et al. 2018) according to the particular scenario of each species (e.g. sequencing and genotyping errors, duplicated loci; Benestan et al. 2016). Since the filtering parameters are dataset dependent (Hendricks et al. 2018), the filtering criteria should be adjusted accordingly (e.g. the stringency of MAC filtering step is sample size dependent). For instance, the number of SNPs was markedly reduced through filtering steps and the highest difference in the percentage of retained SNPs was found among species. In the study by O’Leary et al. (2018) the percentage of retained SNPs ranged from 0 to 63% using the same filtering pipeline with four marine fish species. In our study, the three SNP/RAD-locus filter used to avoid inconsistent RAD-loci could not work well for highly polymorphic species or taxa (e.g. bivalves). Furthermore, the POP filter (i.e. 60% call rate per population) could be applied not so stringently since in previous studies qualitative interpretations of population parameters were maintained in most cases (Shafer et al. 2017; Wright et al. 2019), and sometimes even improved (Hodel et al. 2017). Notwithstanding, the drawback could be using larger SNP panels for similar information. If well, for some studies it is fundamental to achieve the highest density of SNPs possible (e.g. linkage disequilibrium, outlier detection and gene mining, Genome-wide association study; GWAS). The biggest difference between pipelines was found with the MAC filtering step in brown trout which could be explained by the higher average of missing genotypes per SNP (MAC is sample size dependant) and the lower coverage per CASANOVA CHICLANA, ADRIÁN 120 RAD-locus (misclassification of heterozygotes) in the ALT pipeline. Finally, more filtering steps might be necessary, especially when working without a reference genome (e.g. FIS SNP filtering step when paralogs or null alleles can be a problem to avoid misinterpretations); this is the case of HWE deviations in brown trout caused by potential paralogs whose impact can be reduced using a reference genome. Attention should be paid to the order of the different filtering steps because this can alter the final SNP panel. When adjusting the filtering parameters, it would be advisable to consider not exclusively the number of removed SNPs at each step separately, since they could result from the interaction among filtering steps. For instance, the coverage filter determines the increase of missing data which influences the percentage of SNPs eliminated by MAC and population representation filters, according to the stringency of the coverage threshold used. Furthermore, missing data may be due to a lower coverage than the selected threshold or for not being genotyped by the building-loci software with the genotyping options selected (e.g. previously selected nucleotide frequencies range to genotype in ALT pipeline). We found that the last could be the main source of COM SNPs genotyping differences between both building-loci pipelines excluding Manila clam. This means that the ALT pipeline genotyping parameter should be improved by choosing appropriate ranges for each species. The objective of any filtering strategy is removing SNPs that are not reliable without losing informative SNPs. Different factors can influence the filtering criteria, e.g. to achieve the number of SNPs required to meet the research goals. In this sense, a panel made up with markers found by two different pipelines should ensure reliability. It was found that 67% of SNPs from Stacks panel were common with UNEAK panel, using a de novo approach in soybean (Glycine max L.) data (Torkamaneh et al. 2016). With reference genome the overlap percentages among Stacks and other building-loci pipelines ranged Discussion 121 from 76 to 96% (Torkamaneh et al. 2016). Using a reference genome approach the percentages of shared SNPs between Stacks with SAMtools and GATK ranged from 7.3 to 71.4% (Wright et al. 2019). The lowest values could be partially explained because Stacks panel recruited many more SNPs than the other building-loci pipeline. The lowest percentage of COM SNPs taking STA panel as genotyping reference in our study (i.e. 23.9% in small-spotted catshark and 43.0% in Manila clam) were in panels with less than 1,000 SNPs. These low values may be explained by a strong filtering effect, on shared SNPs between pipelines. The highest number of COM SNPs were detected when STA panels included the highest number of SNPs, around 74% in brown trout and 81% in silver catfish. Despite including a lower number of SNPs, the COM panels provided roughly similar results to the larger ones. This suggests that most informative markers are retained downstream, with the advantage of working with a reduced panel that can simplify and speed-up analyses. In the study by Díaz-Arce et al. (2019) the possible effect of SNP number on FST estimation using reduced SNPs subsets was tested and similar values regarding the full panel were obtained. Moreover, estimated genotyping accuracy may be higher with SNPs shared by more than one building-loci pipeline according to Torkamaneh et al. 2016. The impact of genotypic differences between shared SNP panels was low, such as those obtained by Wright et al. (2019). Summarizing, the results obtained suggest that both building-loci pipelines are adequate and provide more confident results adjusting parameters and SNP filtering steps to the research context. Despite the differences observed in the number of SNPs among de novo approach panels, this seems not to affect dramatically the conclusions, at least in the biological scenarios managed in this study. When there is no reference genome, a COM panel could be interesting in terms of SNP panel consistency with species with high genomic complexity. In a CASANOVA CHICLANA, ADRIÁN 122 general way for some population parameters, to have less SNPs do not imply loss of biological information and a COM panel could increase data reliability in these cases. In the case of choosing this option differences in genotyping between pipelines should be checked, although in this study genotyping differences between pipelines were infrequent. The main source of genotyping differences in COM SNP panels was missing data and had different sources. On one hand, different genotypes can be obtained due to the different building-loci pipeline parameters to call genotypes (e.g. --alpha at Stacks) and the different alignment strategies (e.g. reads with alternative allele can be stacked into another putative loci). For missing genotype differences, it should be considered that even RAD-loci showing high coverage, might have missing data if building-loci pipeline parameters involved in genotyping are not properly set up. Furthermore, small differences in the building-loci pipeline could have more influence in the number of missing genotypes when working with low coverage loci. Anyway, it would be advisable to use a few intralibrary and inter-library sample-replicates to estimate genotyping errors (Mastretta-Yanes et al. 2015) to increase the confidence in our data, especially if the RAD-seq libraries are designed with low estimated coverage per locus (e.g. around 10x), due to the impact of coverage in genotyping error rates (Fountain et al. 2016). 5.2 BROWN TROUT POPULATION GENOMICS Impact of restocking in natural populations River drainages of Spain have been massively restocked to counterbalance population depletion since the 70's using a hatchery stock of Central-European origin. This original stock was distributed among all hatcheries across Spanish geography since then, and consequently a slight genetic differentiation can be detected nowadays within an essentially homogeneous genetic background (FST < 0.05; Discussion 123 Almodóvar et al. 2006; García-Marín et al. 2018). Anyway, the genetic divergence between the hatchery stock and Spanish brown trout wild populations is much larger than among hatchery stocks used for restocking (García-Marín et al. 1991; Martínez et al. 1993; Vera et al. 2013), which supports the use of a single hatchery (BA14) as reference in our study to check the incidence of restocking in Atlantic and Mediterranean drainages. The impact of restocking has been evaluated until recently using the nuclear diagnostic marker LDH-C*, fixed for the *90 allele in hatchery stocks and for the *100 allele in wild populations (Morán et al. 1991). This method showed important limitations in scenarios where a moderate or high ancient introgression occurred because of the viability of hybrids and their offspring. In this case, the use of LDH-C* could be inaccurate at population level, and not useful at individual level since LDH-C*100/100 individuals could have an important hatchery genomic background and viceversa. Estimates of introgression based on mtDNA markers have been used as well, nevertheless hatchery mtDNA haplotypes and hatchery nuclear markers do not necessarily match either at individual (Sanz et al. 2006) or at population level (see Table 6 in Plan de gestión de la trucha común en Castilla la Mancha 2019). The use of a high number of SNPs distributed across the whole brown trout genome and the availability of a reference hatchery sample from Bagà, enabled us a more accurate classification of individuals by hatchery ancestry using the probability of membership assignment with STRUCTURE, as previously reported (Hansen et al. 2001; Prado et al. 2018). We applied a conservative q threshold to assign individuals as wild (qW) or hatchery (qH) ancestry by applying a self-assignment test in a population of known genomic background (Bagà). A slightly higher conservative threshold than that used in other investigations on brown trout (Sanz et al. 2009) and turbot CASANOVA CHICLANA, ADRIÁN 124 (Scophthalmus maximus; Prado et al. 2018) was used in accordance with our population scenario. The incidence of restocking detected in our study was variable across the populations studied in the Mediterranean and Atlantic slopes. Although a higher impact was observed in the Mediterranean drainage than in the Atlantic one, as previously reported (Almodóvar et al. 2006), the differences among Mediterranean populations were remarkable, sometimes between populations separated by a few kilometres. Temporal replicates from Ter River (TE) did not show any hatchery incidence, while samples from Queralbs, pertaining to the same river drainage (Ter Basin), were the most affected. A similar observation was reported by Araguas et al. (2017) using five microsatellites and the diagnostic locus LDH-C*. Núria and Queralbs locations are separated by five kilometres in the same basin, but currently could be partially isolated by different barriers (e.g. dam in Daió hydroelectric power plant). Such barriers might favour the genetic integrity of the brown trout from Ter River. Dams can act as a barrier for alien species invasion and in some cases, they have been built to protect native species (Dana et al. 2011). In the Atlantic drainage, the same individuals from BL were removed by in Martínez et al. (2007) being considered as pure hatchery (BL-15) and F1 (BL-24) using the LDH-C* locus, while here a more refined genomic constitution of those individuals was achieved. Usually, to classify individuals as wild or hatchery introgressed to take management decisions, a global ancestry approach using tools as STRUCTURE was employed. Nevertheless, local ancestry inference would be more powerful to detect small chromosomic regions affected by hatchery introgression. As in the study of Leitwein et al. (2018), some individuals previously identified as “pure wild” showed “introgression signals” with local ancestry approach. Chromosomic regions that may be prone or resistant to introgression were detected in our study. However, to avoid interpreting sampling effects as Discussion 125 introgression patterns, more admixed samples should be analysed and compared to establish consistent resistant or prone regions to introgression. We intended to develop a more powerful, informative and cheaper molecular tool to evaluate restocking in Spanish drainages by using different panels of progressively higher power to detect hatchery ancestry individuals. We assumed that the most consistent qH values were those obtained with the whole SNP dataset, thus being used as reference for evaluating the performance of the different subsets. However, low informative SNPs might be filtered according to FST values defining subsets of high-resolution SNPs. So, we identified the most informative markers using FST between wild and hatchery samples and detected nine diagnostic SNPs capable to identify most hatchery ancestry individuals in our sample. A second and a third subset of 38 (FST > 0.99) and 214 (FST > 0.95) SNPs were also evaluated for the higher resolution and the lower cost as possible. The performance of the three SNP subsets for the classification of individuals according to hatchery ancestry using, either the STRUCTURE approach or the proportion of hatchery alleles at individual level using diagnostic loci, was remarkable and the correlations obtained with the whole genomic data were highly significant for the three subsets (p-value < 0.001). Some discrepancies at population and individual level were related to individuals with a slight hatchery ancestry, especially in CE population. Further, the hybridization and introgression across the genome is not homogeneous (see for instance Wang et al. 2019 for soybean) and a small number of AIMs (Ancestry Informative Markers) may not be enough. These SNPs, identified in silico from RAD-seq, should be validated in the future with techniques well fitted to handle small SNP panels (e.g. SNaPshot, Sequenom, TaqMan) to devise the best molecular tool combining statistical power and low price. The current price of RE-digestion for LDH-C* genotyping is close to 10€ / CASANOVA CHICLANA, ADRIÁN 126 individual, while a single multiplex for Sequenom could even be cheaper. Similar molecular tools have been successfully applied for the identification of hybrids in aquatic organisms (Maroso et al. 2018, 2019). Nevertheless, these SNPs should be tested with new individuals and more samples covering a wider distribution of the Iberian Peninsula to evaluate their performance at individual and population level regarding the previous results obtained with the LDH-C* locus. Genetic diversity Preservation of genetic diversity is a key point to maintain the potential for adaptation of natural populations to environmental changes (Frankham et al. 2010), which are rapidly affected by human activities and the ongoing climate change. Since the 90's until recently, the most commonly genetic markers used to estimate genetic diversity for conservation and management of bioresources were microsatellites (Saint-Pé et al. 2019). A lot of studies have been performed with microsatellite markers even during the last 10 years (Araguas et al. 2017; Berrebi et al. 2019; Vera et al. 2013, 2018; Vilas et al. 2010). The arrival of NGS techniques and related techniques have allowed a quick identification and cheap genotyping of thousands of SNPs, even without reference genomes, the so-called genotyping by sequencing (GBS) techniques (Davey et al. 2011; Robledo et al. 2018). The arrival of SNPs has let population and individual screening of genomes for more accurate estimation of genetic diversity and structure, and especially, the identification of footprints of selection (Bernatchez 2016). The fact that SNPs are usually biallelic markers while microsatellites are hypervariable with tens of alleles per locus, makes that genetic diversity estimators cannot be directly comparable within populations (Bouza et al. 2001; Saint-Pé et al. 2019; Vilas et al. 2010), but still the relative genetic diversity among populations can be compared with both types of markers. Discussion 127 The average figures of genetic diversity found in the present study for the Iberian Peninsula were always much lower than those reported with microsatellites (Bouza et al. 2001; Martínez et al. 2007) for all populations (He = 0.064 ± 0.034), with Galician populations (0.106 ± 0.020) showing higher diversity than Duero (0.045 ± 0.024) and Mediterranean populations (0.040 ± 0.004). Low genetic diversity in Mediterranean Slope is expected due the important impact of genetic drift due to unstable hydrology and isolation of river basins (Araguas et al. 2017; Vera et al. 2013). Conversely, populations above the parallel 42º N, including the Galician region, are expected to show higher genetic diversity due to stable hydrology and inter-basin connection through the anadromous trout, the so-called sea trout (Antunes et al. 2006; Bouza et al. 1999; García-Marín et al. 2018; Östergren and Nilsson 2012). The migration of sea trout would explain a higher effective population size (Ne), and accordingly, higher genetic diversity. Our results, therefore, meet to previous observations. However, comparison with other studies in northern regions have shown much higher genetic diversity than that found here. For instance, a wide survey across Northern European populations of S. trutta (72 locations placed in Great Britain, Germany, and Scandinavia) with 3,872 SNPs showed an average He > 0.30 (Bekkevold et al. 2020), three times higher than that found in Galicia in our study. These differences found among populations of the same species with SNPs has usually to do with the filtering performed in each study. Thus, the SNPs used by Bekkevold et al. (2020) were genotyped with a SNP-chip where highly polymorphic SNPs had been selected, so rendering an average MAF of 0.28 and only 1.75% showing MAF < 0.05; this unavoidably upwards He estimations. Similar genetic diversity to that found by Bekkevold et al. (2020) was obtained in Southern Baltic area by Bernaś et al. (2020) with 3,843 SNPs in a region where sea trout also occurs. Average MAF in our study was 0.13, with a much higher proportion of loci with MAF < 0.05 (44.0%). Many of these alleles were private of specific CASANOVA CHICLANA, ADRIÁN 230 Figure S9. Comparison between CLUMPAK and DAPC outputs for small-spotted catshark (S. canicula) samples (N = 28). Two locations are included: IS (Irish Sea), NS (North Sea).