Functional Genomic Characterization of Viral Nervous Necrosis Disease resistance in farmed European seabass

In this study, we aimed to functionally characterize the genomic background of viral nervous necrosis resistance in farmed European seabass. Offspring from a cross of 25 males and 25 females were challenged with nervous necrosis virus and their response to the infection was recorded as days to death or binary survival (dead or alive) for about 29 days. A total of 1390 European seabass fish including 50 parents (25 dams and 25 sires) and 1340 juvenile offspring were genotyped on a MedFish Axiom Affymetrix SNP array for 29,888 SNP markers. The SNP array genotype profiles were first subjected to initial genotype call quality control using Axiom Analysis Suite v4.0.3.3, and at a threshold of 97% genotype call rate a total of 1362 fish and 28,404 were retained. This genotype data was further subjected to quality control analysis using Plink1.9 software, where SNP SNPs with genotype call rate across samples (--geno) < 95%, minor allele frequency (--maf) < 5%, and significantly deviating from Hardy-Weinberg equilibrium (HWE < 1E-4) were removed from the analysis. Consequently, a total of 24,740 SNP remained for subsequent analyses and this data was utilized to reconstruct the pedigree of all the offspring using the APIS v1.0.1 package in R. A section of this fish, 50 parents, and 40 offspring were whole-genome sequenced (paired, 150bp) at a coverage of 15X. Firstly, sequence data was assessed for sequencing quality using FastQC v0.11.9 and cleaned of adaptors using Trimmomatic v0.39 software. Then cleaned sequence data was aligned to the most recent European seabass reference genome (GCA_905237075.1) using the BWA v0.7.16 short read aligner and subsequently SNP and Indels variants and their genotypes were called from the alignments using BCFtools v1.9. Using VCFtools v0.1.13 the genotype data was cleaned of variants that had: 1) more than two alleles (--max-alleles 2), 2) genotype calls in less than half of the samples (--max-missing 0.5), 3) genotype quality less than 30 (--minQ 30), 4) low read coverage (--minDP 8), 5) very high coverage (--maxDP 31.05), and 6) minor allele count greater or equal to two (--mac 3). Also, the data was cleaned of variants that significantly deviated from Hardy-Weinberg equilibrium (--hwe 1E-7) in the parents. As a result, a total of 8,368,178 variants (8.4M variants) were retained after filtering. The whole genome sequence variant genotypes of the 90 samples were utilized as the reference used to impute all the offspring (1 fish) from SNP array genotype dataset to whole genome variant genotypes of 8.4M variants (SNPs and Indels) using Fimpute v3 software with the pedigree relationships between the animals supplied to the software. The whole genome data generated from this imputation together with phenotype data of survival and days to death were then used for WGS genome-wide association study analysis using gcta v1.93.2beta software. These genotypes were further utilized in a targeted expression quantitative loci identification for the candidate genes located within the genomic region identified as associated with nervous necrosis resistance in this study via eGWAS analyses implemented in gcta v1.93.2beta software. This work was funded by the European Union's Horizon 2020 research and innovation program under grant agreement No 817923 (AQUAFAANG).