Use of genomic tools for the study of the genetic basis of resistance to gastrointestinal nematode infections in adult sheep

Resumen

La actual tesis doctoral, que se ha desarrollado en el grupo de investigación de mejora genética animal de la Universidad de León (grupo MEGA-ULE), tiene como objetivo global el estudio de la base genética subyacente a la resistencia a los nematodos gastrointestinales (GIN) en ovejas adultas. Para ello se han utilizado metodologías empleadas en la genómica, tales como el chip de genotipado de SNPs de alta densidad (HD) y las tecnologías se cuenciación conocidas como next generation sequencing (NGS). Para el cumplimiento el objetivo global se han abordado dos objetivos específicos. El primero de ellos tuvo por objeto realizar un mapeo de alta densidad de regiones con influencia sobre la resistencia a los GIN en ovejas mediante el análisis de genotipos imputados para el chip HD-chip en una población de raza Churra previamente analizada con el chip de media densidad (50K-chip) y para la cual estaban disponibles datos de dos fenotipos relacionados con la carga parasitaria. El segundo objetivo específico se centró en la caracterización funcional de regiones genómicas asociadas con la resistencia a los GIN en ovejas adultas, utilizando el análsis del transcriptioma completo de diferentes tejidos (ganglio linfático abomasal y mucosa del abomaso) después de la infección experimental con T. circumcincta, mediante la metodología RNAseq. Dentro de las actividades necesarias para cumplir el primer objetivo concreto, el primer paso fue la imputación de genotipos del chip de 50K a la densidad del chip de HD (600K), realizándo previamente una estimación de la precisión de la estrategia de imputación utilizada. La imputación de genotipos se realizó en una población comercial de ovejas de raza Churra con genotipos disponibles para la densidad media de 50K-chip. Esta población estaba formada por 16 familias de medio hermanas, con un total de 1.670 hijas. Los 1.686 animales de la población se genotiparon con un chip de 50K SNP y un subconjunto de 240 animales de esa población, más otros 96 machos de inseminación artificial se genotiparon, además, con un chip HD (600K). Utilizando los dos conjuntos de datos genotípicos disponibles, los genotipos del 50K-chip para la población global fueron imputados al chip HD utilizando el software Beagle. La concordancia media obtenida en diferentes estudios y escenarios para los 26 cromosomas autosómicos, fue de alrededor el 90-92%, que está dentro del rango de precisión de imputación reportado por otros estudios genómicos de ovejas utilizando Beagle (83-93%).Los resultados sobre la evaluación de diferentes estrategias de imputación se han presentado en el Estudio 1 donde los genotipos imputados para el chip de alta densidad sirvieron para evaluar el nivel de desequilibrio ligamiento (LD) en el genoma de la raza Churra, identificar regiones conocidas como runs of homozygosity (ROH) y actualizar las estimaciones del tamaño efectivo (Ne) de la población con una mayor precisión aportada por la alta densidad de marcadores. Las estimaciones actualizadas de LD ponen de manifiesto que el grado de LD en la raza ovina Churra es aún más corto que el descrito en base al 50K-chip, siendo una de las extensiones más cortas descritas en razas ovinas. El coeficiente de consanguinidad, teniendo en cuenta la longitud total del ROH, mostró una estimación media (0,0404) inferior al nivel crítico. En general, la mayor precisión de las estimaciones actualizadas de LD sugiere que el HD-chip, combinado con una adecuada estrategia de imputación, ofrece una poderosa herramienta en estudios de mapeo genético o la aplicación de la selección genómica en ganado ovino. Aprovechando la estrategia de imputación optimizada, y considerando la limitación del uso de la selección genómica en ovejas lecheras debido a los costos de genotipado, el Estudio 2 describe el desarrollo in silico de un chip de baja densidad (LD-chip) y el potencial de esta herramienta genómica para imputar genotipos a media (50K-chip) y alta (HD-chip) densidad de marcadores. Para ello se utilizaron los genotipos del 50K-chip y HD-chip disponibles para la población global de ovejas Churra estudiada en el Estudio 1. Así, el LD-chip virtual diseñado incluyó un total de 2.935 marcadores comunes a los dos paneles y distribuidos uniformemente a través del genoma ovino. La precisión media para la imputación del LD-chip a 50K-chip osciló entre 93,53% y 93,58%, y para la imputación del LD-chip al HD-chip osciló entre 88,3% a 86,52% (dependiendo del nivel de enmascaramiento considerado, 10% y 30%, respectivamente). A nivel global, y aunque será necesario resolver ciertas limitaciones identificadas, el LD-chip descrito en el Estudio 2, junto con la sugerida estrategia de imputación, parecen aportar una aproximación apropiada para obtener información para pruebas genéticas básicas (por ejemplo, paternidad), y también para estimaciones de los valores de cría genómicos que podrían ser utilizados para guiar la toma de decisiones en relación a la selección de candidatos a la selección. La última actividad del primer objetivo fue la realización de estudios de mapeo genético de alta resolución. En primero lugar se realizó un estudio GWA con los genotipos imputados del HD-chip obtenidos para las 532 hijas de la población comercial y sus fenotipos disponibles para caracteres indicadores de resistencia a GIN. En el Estudio 3 se presentan los resultados del análisis GWAS para el caracter “niveles séricos de Inmunoglobulina A” (IgA), que identificó unicamente un SNP significativo al nivel de signficación 5% chromosome-wise, localizado en OAR15. Posteriormente, con el fin de explotar la estructura de la población de medio hermanas considerada en este trabajo realizamos un escaneo del genoma basado en el análisis de ligamiento (LA), que se describe en el Estudio 4 para los dos fenotipos indicadores de resistencia a la infección por GIN medidos en esta población, IgA y el recuento de huevos en heces (FEC). Este análisis ha permitido refinar los intervalos de confianza (CI) de los regiones QTL (del inglés Quantitative Trait Loci) previamente identificadas en oveja Churra para caracteres de resistencia a GIN en base al 50K-chip. Además, en el Estudio 4, la variabilidad genética del nuevo CI refinado para el QTL más prometedor, identificado en el cromosoma 6 (OAR6) para FEC, se estudió a alta resolución utilizando la información derivada de la WG-seq para un trío segregante para ese QTL (que incluyó el padre Qq y 2 hijas con fenotipos divergentes extremos para el rasgo FEC en concordancia con sus genotipos para el QTL, hija QQ e hija qq). En este estudio se han identificado un total de 433 variantes de alta calidad dentro de la región considerada en OAR6 (88,2-88,3 Mb), de las cuales 357 eran variantes intragénicas concordantes con el genotipo del QTL para los tres animales analizados. Basado en un análisis de anotación funcional de variantes, el Estudio 4 proporciona una lista de variantes que podrían considerarse como posibles mutaciones asociadas con la resistencia ovina a las infecciones por GIN. El segundo objetivo relacionado con la caracterización funcional de las regiones genómicas responsables de la resistencia a la GIN en ovejas adultas se basó en el estudio transcriptómico global a través de la secuenciación de RNA de muestras de tejido abomasal obtenidas para 12 ovejas adultas de raza Churra. Sobre la base de las medidas individuales de FEC en una infección experimental (EI1) con T. circumcincta , seis ovejas fueron clasificadas como susceptibles, y seis ovejas fueron clasificadas como resistentes a GIN. Posteriormente, estos animales fueron expuestos a una segunda infección experimental (EI2) con T. circumcincta. En el día 7 tras EI2, los animales fueron sacrificados y la necropsia permitió la recogida, para cada animal, de muestras de mucosa abomal y ganglio linfático abomasal, a partir de las cuales se obtuvieron muestras de mRNA que fueron secuenciadas mediante metodología RNA-seq. Utilizando las muestras de RNA-seq de la mucosa abomasal se realizó un análisis preliminar de expresión diferencial comparando las muestras de mucosa abomasal de ovejas clasificadas como resistentes o susceptibles a la infección experimental con T.circumcincta (Estudio 5). El análisis de expresión diferencial mediante los softwares DESeq2 y edgeR no mostró resultados significativos comunes, sugiriendo que a los 7 días tras la infección con T. circumcinta no hay una respuesta diferencial clara en relación al perfil transcripcional en la mucosa abomasal de animales resistentes y susceptibles. Debido a eso, el Estudio 6 se centró principalmente en describir los resultados del análisis DE realizado para las muestras de ganglio linfático abomasal. Para este tejido, el análisis de expresión realizado con los dos programas mencionados anteriormente identificó un total de 106 DEGs comunes, que fueron considerados como genes activados por GIN. Los análisis de enriquecimiento realizados para estos genes identificó como significativas algunas rutas relacionadas con la respuesta inmune mediada por citoquinas y la ruta de señalización de PPARG (Peroxisome proliferator-activated receptor gamma), además de términos de enfermedades relacionados con la inflamación y las enfermedades gastrointestinales. La comparación sistemática de los genes identificados como activados por GIN con resultados de estudios funcionales anteriores confirmó la implicación de genes como ITLN2, CLAC1 y galectinas, en los mecanismos inmunes activados en ovejas resistentes frente a T. circumcincta. Por otra parte, la información generada mediante la tecnología de RNA-seq para las muestras de ganglios linfáticos fue utilizada para identificar mutaciones con un posible impacto funcional en regiones previamente identificadas como asociadas a la resistencia parasitaria en el ganado ovino. Para ello, se realizó un análisis de identificación de variantes considerando conjuntamente todos los datos de transcriptoma de las muestras de ganglios linfáticos. El estudio inicial de la variabilidad genética identificada a lo largo de todo el transcriptoma y en los 106 genes previamente identificados como activados por GIN, descrito en Estudio 7, fue seguido por un análisis de anotación funcional más extenso descrito en el Estudio 8. En ambos estudios se realizó análisis de identificación de variantes con dos tipos diferentes de software, GATK_v3.7 y Samtools_v1.4 y las variantes comúnmente identificadas por los dos paquetes se consideraron "variantes de alta calidad". En el Estudio 7 se identificaron 1.326.960 variantes de alta calidad en todo el transcriptoma. De estas variantes de alta calidad, se extrajeron un total de 6.168 variantes dentro de los 106 DEG identificados en el Estudio 6. En el Estudio 8, además, se identificaron las variantes de alta calidad con un posible impacto funcional en dos tipos de regiones de interés, i) regiones QTL previamente descritas en la oveja para caracteres de resistencia a parásitos; ii) genes candidatos funcionales seleccionados a partir de estudios de expresión génica relacionados con la resistencia a las infecciones por GIN en ovejas. En general, este estudio proporciona una valiosa lista de posibles mutaciones causales que podrían considerarse como mutaciones causales candidatas relativas a la resistencia a las infecciones por GIN en ovejas. En el Estudio 9, los datos de RNA-seq obtenidos para ganglios linfáticos abomasales fueron sometidos a un nuevo análisis con el objetivo de identificar RNA largos no codificantes (lncRNAs) que pudieran asociarse a la resistencia a los GIN en ovejas. Utilizando un flujo de análisis bioinformático específico para la detección de este tipo de transcrito basado en el paquete FEELnc, se identificaron 44.203 genes anotados y un total de 77.039 transcritos, de los cuales 9.105 fueron clasificados como lncRNAs (2.092 de los cuales fueron considerados como novel). Un análisis adicional de expresión difirencial permitió identificar un total de 3.148 DEGs, 457 de los cuales estaban asociados a 683 de los lncRNA previamente identificados. Considerando los 367 lncRNAs para los que se identificaron 301 genes solapantes, este trabajo proporciona una nueva lista de genes activados por GIN que podrían estar regulados por lncRNA. Para algunos de esos genes (IRF1, CCL14, KRT8 y LGALS3) se ha descrito previamente una posible relación con la respuesta del huésped contra la infección por nematodos en ovejas o en ganado bovino. En conjunto, esta Tesis doctoral representa un avance hacia el conocimiento de la base genética de la resistencia a las infecciones por GIN en ovejas en dos sentidos: i) proporcionando un conocimiento más profundo sobre la organización estructural del genoma de las ovejas de raza Churra (mediante el análisis de la extensión de la LD, etc) y aumentando la densidad de marcadores de barridos genómicos previamnete realizaods por el grupo MEGA-ULE para caracteres de resistencia a GIN, y también (ii) al abordar, por primera vez en ovejas adultas, un enfoque funcional a través del análisis global del transcriptoma de tejidos afectados por la infección por GIN, para el estudio de los mecanismos genéticos subyacentes al fenotipo de resistencia. Como valor añadido, la estrategia optimizada para la imputación de genotipos podría suponer un primer paso hacia la implementación de la selección genómica en las poblaciones comerciales de razas lecheras españolas. Summary:The present PhD Thesis, which has been developed in the Animal Breeding research group of the University of León (MEGA-ULE group), was set out with the global objective to identify the genetic basis underlying resistance to gastrointestinal nematode (GIN) infection in adult sheep, by taking advantage of the latest available genomic technology such as the Illumina ovine High-Density Bead chip (HD-chip) and the costeffective next-generation sequencing (NGS) technologies. This global objective was achievedby addressing two parallel specific objectives. The first specific objective was focused on a high-density mapping of genomic regions underlying resistance to GIN in sheep using imputed HD SNP genotypes in a population previously analysed with the medium density SNP-chip (50K-chip) and for which two indicator traits of parasite infection level were available. The second objective was focused on the functional characterization of genomic regions that are responsible for resistance to GIN in adult sheep using RNA sequencing (RNA-seq) data extracted from abomasal tissues (abomasal mucosa and abomasal lymph node tissues) of adult Churra sheep after experimental infection with GIN T. circumcincta. Within the activities required to fulfil the first specific objective, the first step was the imputation of genotypes from the 50K-chip genotypes to the HD-chip density (600K), with a previous estimation of the accuracy of the imputation strategy utilized. The analysis was performed in a commercial population of Churra sheep with available genotypes for the medium density of 50K-chip. This population consisted of 16 half-sib families, with a total of 1,670 daughters. The 1,686 animals in the population were genotyped with a 50K-chip and a subset of 240 animals of that population, together with 96 artificial insemination sires, were also genotyped with a HD-chip (600K). Using the two genotype datasets, the 50K-chip genotypes of the global population were imputed to the HD-chip using the Beagle software. The average concordance obtained across different studies and scenarios for the 26 autosomal chromosomes, was around 90-92%, which is within the imputation accuracy range reported by other sheep genomic studies using Beagle (83-93%). The results on the evaluation of different imputation strategies have been presented in Study 1, where the imputed HD-chip genotypes were used to update the estimates of LD, identify regions involving runs of homozygosity (ROH) and update the estimation of effective population size (Ne) for the studied population. The updated LD estimations provided evidence that the extent of LD in Churra sheep is even shorter than that reported based on the 50K-chip and is one of the shortest extents compared with other sheep breeds. Through different comparisons, we have also assessed the impact of imputation on LD and Ne estimates. The inbreeding coefficient, considering the total length of the ROH, showed an average estimate (0.0404) lower than the critical level. In general, the higher accuracy of updated LD estimates suggests that the HD-chip, combined with an appropriate imputation strategy, offers a powerful tool for gene mapping studies or the implementation of genomic selection in sheep. By taking advantage of the imputation strategy optimized, and considering the limitation of using genomic selection in dairy sheep due to the genotyping costs, Study 2 reports the development of an in silico low-density chip (LD-chip) and the potential of this genomic tool to impute genotypes to medium (50K-chip) and high (HD-chip) marker densities. For that, we used the 50K and HD-chip genotypes available for the global Churra sheep population studied in Study 1. A total of 2,935 markers common to the two panels and evenly distributed across the sheep genome were included in the virtual LDchip. The average accuracy for the LD-chip to 50K-chip imputation ranged from 93.53% to 93.58%, and for the LD-chip to HD-chip imputation ranged from 88.3% to 86.52% (depending on the masking level considered, 10% and 30%, respectively). Globally, and although certain identified limitations need to be solved, the in silico LD-chip described in Study 2, together with the implemented imputation strategy, appears to provide an appropriate approach to obtain information for basic genetic tests (e.g. paternity), and also estimations of genomic breeding values that could guide breeding decisions. The last actitivity performed within the first objective was the performance of gene mapping studied based on the imputed HD-chip genotypes. First, we performed a GWA analysis based on the the imputed HD-chip genotypes obtained for the 532 daughters of the commercial population with available phenotypes for GIN resistance indicator traits. In Study 3 we present the results of the GWA study performed for the indicator trait Immunoglobulin A serum levels (IgA), which identified a single significant SNP at the 5% chromosome-wise level, located on OAR15. To exploit the half-sib structure of the studied resource population, we later performed a genome scan based on Linkage analysis (LA), which is described in Study 4 for two indicator traits of GIN resistance, IgA and faecal egg count (FEC). Also, in Study 4, the genetic variability of the new refined CIs for the most promising QTL, identified for FEC on chromosome 6 (OAR6), was studied at high resolution by using the information derived from the analysis of WG-seq datasets of a segregating trio for that QTL (which included the Qq sire, and 2 daughters with extremely divergent phenotypes for the FEC trait in concordance with their QTL inferred QTL genotypes, QQ and qq). In this analysis, a total of 433 high-quality variants were identified within the OAR6 CI QTL region (88.2-88.3 Mb), 357 of which were intragenic concordant variants for the three animals analysed. Based on a variant functional annotation analysis, Study 4 provides a list of variants that could be considered as potentially underlying GIN resistance in sheep. The second objective related with the functional characterization of genomic regions responsible for GIN resistance in adult sheep was based on the global transcriptomic study through RNA sequencing of abomasal tissue samples obtained for 12 adult ewes of Churra sheep after performing a first experimental infection (EI1) with GIN T.circumcincta. Based on the individual accumulated FEC measures after EI1, six sheep were classified as susceptible, and six sheep were classified as resistant to GIN. Later the selected animals were exposed to a second experimental infection (EI2) with T. circumcincta. At day 7 after EI2, the animals were sacrificed and at necropsy, abomasal mucosa and lymph node samples for all animals were immediately collected and the mRNA extracted for two tissues, abomasal mucosa and abomasal lymph node. Using the abomasal mucosal RNA-seq samples an initial preliminary differential transcriptomic analysis was performed between the abomasal mucosa samples obtained from two groups of sheep classified as resistant or susceptible to the experimental infection with T. circumcincta (Study 5). The differential expression analysis performed with two different R-based packages, DESeq2 and edgeR, did not show concordant results, suggesting that at 7 days post-infection, there is not a clear differential gene response in the abomasal mucosa between resistant and susceptible ewes. Because of that, Study 6 was focused on reporting the results of the DE analyses performed for the abomasal lymph samples. For this tissue, the expression analysis performed with the two previously mentioned software identified a total of 106 common differentially expressed genes (DEGs), which were considered as GIN-activated. The enrichment analyses performed for these GIN-activated genes identified some immune related and the PPARG (Peroxisome proliferator-activated receptor gamma) signaling pathway as well as disease terms related to inflammation and gastro-intestinal diseases. A systematic comparison with the results of previous studies confirmed the involvement of genes such as ITLN2, CLAC1 and galectins, in the immune mechanism activated against T. circumcincta in resistant sheep. In addition, the information generated through the RNA-seq technology for the lymph node samples was used to identify potential causal mutations related to GIN infection in sheep, by performing a variant calling analysis on this transcriptome dataset. The initial study of the genetic variability identified throughout the transcriptome and in the 106 genes previously identified as GIN-activated, described in Study 7, was followed by a more extensive functional annotation analysis described in Study 8. In both studies variant calling analyses were performed with two different types of software, GATK_v3.7 and Samtools_v1.4 and the variants commonly identified by the two different software packages were considered “high-quality variants”. In Study 7 we identified 1,326,960 high-quality variants across the whole transcriptome. From these high-quality variants, we extracted 6,168 variants within the 106 DEGs identified in Study 6. In Study 8, in addition, the high-quality variants with a potential functional impact were identified for two types of target regions, (i) QTL regions previously reported in sheep for parasite resistance and (ii) functional candidate genes selected from gene expression studies related to GIN resistance in sheep. Overall, the study provides a valuable list of potential causal mutations that could be considered as candidate causal mutations concerning GIN resistance in sheep. In Study 9, we continued exploiting the transcriptomic dataset obtained by RNA-seq for the abomasal lymph node tissue previously analysed to identify long noncoding RNAs (lncRNAs) that could be associated to GIN resistance in sheep. Using a specific bioinformatic pipeline for the identification of lncRNA based on the FEELnc package, the analysis identified 44,203 genes with 77,039 transcripts, 9,105 of which were identified as lncRNAs (2,092 of them were classified as novel). An additional differential gene expression analysis, identified a total of 3,148 DEGs, 457 of which were associated with 683 of the identified lncRNA transcripts. Focusing on the 367 lncRNA transcripts that were later found to include 301 overlapping genes, this work has highlighted a new list of GIN-activated genes that could be regulated by lncRNA, some of which (IRF1, CCL14, KRT8, and LGALS3) have been previously reported to influence on the host response against the nematode infection in sheep or cattle. Based on all the above mentioned, this PhD Thesis has provided a step forward towards the identification of genes underlying the resistance to GIN infections in sheep in two senses: (i) by providing a deeper knowledge on the structural organization of the Churra sheep genome (through the analysis of the LD extent, and related analyses) and increasing the marker density of previous genome scans performed by the MEGA-ULE group for GIN resistance, and (ii) by addressing, for the first time in adult sheep, a functional approach for the study of the genetic mechanisms underlying the GIN resistance phenotype in sheep. As an added value, the optimized imputation strategy of genotypes reported here may define a first step towards the implementation of genomic selection in the commercial populations of Spanish dairy sheep breeds.