Identification of nonsense variants in the genomes of 15 Murciano-Granadina bucks and analysis of their segregation in parent-offspring trios

11224 ABSTRACT Nonsense variants can inactivate gene function by causing the synthesis of truncated proteins or by inducing nonsense mediated decay of messenger RNAs. The occurrence of such variants in the genomes of livestock species is modulated by multiple demographic and selective factors. Even though nonsense variants can have causal effects on embryo lethality, abortions, and disease, their genomic distribution and segregation in domestic goats have not been characterized in depth yet. In this work, we have sequenced the genomes of 15 Murciano-Granadina bucks with an average coverage of 32.92× ± 1.45×. Bioinformatic analysis revealed 947 nonsense variants consistently detected with SnpEff and Ensembl-VEP. These variants were especially abundant in the 3′ end of the protein-coding regions. Genes related to olfactory perception, ATPase activity coupled to transmembrane movement of substances, defense to virus, hormonal response, and sensory perception of taste were particularly enriched in nonsense variants. Seventeen nonsense variants expected to have harmful effects on fitness were genotyped in parent-offspring trios. We observed that several nonsense variants predicted to be lethal based on mouse knockout data did not have such effect, a finding that could be explained by the existence of multiple mechanisms counteracting lethality. These findings demonstrate that predicting the effects of putative nonsense variants on fitness is extremely challenging. As a matter of fact, such a goal could only be achieved by generating a high-quality telomere-to-telomere goat reference genome combined with carefully curated annotation and functional testing of promising candidate variants. Key words: knockout, stop gain variant, protein truncation, cost of domestication hypothesis INTRODUCTION There is evidence that deleterious variants accumulated in the genomes of domestic animals and plants at a higher rate than in their wild ancestors, a circumstance that imposed a cost on their fitness (Moyers et al., 2018). In populations experiencing bottlenecks, such as the domestic ones, purifying selection is less able to remove variants with adverse effects on fitness, genetic drift being the key evolutionary force deciding their fate (Moyers et al., 2018). Inbreeding, linkage with variants under artificial selection, and fast population growth combined with long-distance migration also contributed to increase the frequencies of deleterious variants in domestic populations (Moyers et al., 2018). Notably, the linkage between deleterious variants and those with causal effects on phenotypes of interest might impose limits on the response of such traits to artificial selection (Moyers et al., 2018). Among deleterious variation, nonsense variants are of particular importance because they can inactivate protein function by inducing either the premature cessation of translation or nonsense mediated decay of the corresponding transcript (MacArthur and Tyler-Smith 2010). Variants at canonical splice sites and missense substitutions can also have adverse effects on gene function, but in general these are less drastic than those of nonsense variants (Abramowicz and Gos 2018; González-Prendes et al., 2023). Indeed, nonsense variants are among the most frequent causes of genetic disease, representing ~11.5% of human inherited diseases and 20% of diseaseassociated single-nucleotide substitutions residing in gene coding regions (Mort et al., 2008). In addition, premature stop codons might also have negative effects on Identification of nonsense variants in the genomes of 15 MurcianoGranadina bucks and analysis of their segregation in parent-offspring trios Ke Wang,1,2,3 María Gracia Luigi-Sierra,1 Anna Castelló,1,3 Taina Figueiredo-Cardoso,1,3 Anna Mercadé,3 Amparo Martínez,4 Juan Vicente Delgado,4 Javier Fernández Álvarez,4 Antonia Noce,1 Mingjing Wang,1 Jordi Jordana,3 and Marcel Amills1,3* 1Centre de Recerca Agrigenòmica (CRAG), CSIC-IRTA-UAB-UB, Campus Universitat Autònoma de Barcelona, Bellaterra 08193, Spain 2Chinese Academy of Tropical Agricultural Sciences, Zhanjiang Experimental Station, Zhanjiang, Guangdong 524000, China 3Departament de Ciència Animal i dels Aliments, Facultat de Veterinària, Universitat Autònoma de Barcelona, Bellaterra 08193, Spain 4Departamento de Genética, Universidad de Córdoba, Córdoba 14071, Spain J. Dairy Sci. 107:11224–11238 https://doi.org/10.3168/jds.2024-24952 © 2024, The Authors. Published by Elsevier Inc. on behalf of the American Dairy Science Association®. This is an open access article under the CC BY license (https://creativecommons.org/licenses/by/4.0/). The list of standard abbreviations for JDS is available at adsa.org/jds-abbreviations-24. Nonstandard abbreviations are available in the Notes. Received March 25, 2024. Accepted July 31, 2024. *Corresponding author: marcel.amills@ uab .cat 11225 Journal of Dairy Science Vol. 107 No. 12, 2024 reproduction by impairing embryonic development and inducing abortions (Derks et al., 2019). The potential consequences of nonsense variants have often been predicted based on phenotypic data gathered in knockout mice (Majzoub and Muglia 1996). However, there is increasing evidence that phenotypes associated with gene inactivation in mice are not always recapitulated in humans homozygous for alleles inducing the truncation of the corresponding protein (Narasimhan et al., 2016). This could be due to multiple mechanisms, such as gene redundancy (El-Brolosy and Stainier 2017), alternative splicing (Narasimhan et al., 2016), ribosome readthrough (Li and Zhang 2019), and full linkage with suppressor variants (Matsui et al., 2017), to mention a few. In this work, we wanted to characterize the genomic landscape of nonsense variants in goats by sequencing the genomes of 15 Murciano-Granadina bucks. Moreover, we aimed to select 20 nonsense variant candidates to have harmful effects on viability, according to information collected in the Mouse Genome Informatics (MGI; Baldarelli et al., 2021) and Online Mendelian Inheritance in Man (OMIM; Hamosh et al., 2005) database annotations, and to genotype them in buck-dam-offspring trios to investigate if any of them displays an altered pattern of transmission. MATERIALS AND METHODS Ethics Statement The collection of blood is a routine procedure carried out by trained veterinarians working for the National Association of Murciano-Granadina Goat Breeders (CAPRIGRAN), so it does not require a permission from the Committee on Ethics in Animal and Human Experimentation of the Universitat Autònoma de Barcelona. Permission to extract ear samples was granted by the Review Committee for the Use of Animal Subjects of Northwest A&F University (contract number NWAFAC1008). Animal Material We have extracted blood samples from 6 MurcianoGranadina bucks, 284 dams mated to the 6 bucks, and 302 offspring produced by such matings. These 6 bucks were selected because they had ~40 to 50 offspring available. Moreover, blood was also retrieved from 9 additional Murciano-Granadina bucks as part of an experiment unrelated to the current study. Blood was kept in EDTA K3 coated vacuum tubes and stored at −20°C before processing. Murciano-Granadina is an autochthonous Spanish breed officially created in 1975 by the crossbreeding of 2 different Murciano and Granadina goat populations. Currently, the breed has a census of 115,105 heads (2020), and its adaptability to harsh environments as well as its good milking records (mean of 586 kg/ lactation; 5.1% of fat and 3.6% of protein in milk) have made it a highly appreciated breed in Spain and other countries (https: / / www .mapa .gob .es/ es/ ganaderia/ temas/ zootecnia/ razas -ganaderas/ razas/ catalogo -razas/ caprino/ murciano -granadina/ ). Part of the work developed in this project was carried out in the Zhanjiang Experimental Station of the Chinese Academy of Tropical Agricultural Sciences, located in Zhanjiang, Guangdong Province. Because of that, we had access to a collection of genomic DNA preparations from 250 ear tissues from 4 native Chinese breeds (Shaanbei White Cashmere goats, n = 50; Guanzhong goats, n = 50; Leizhou black goats, n = 50; Hezhang black goats, n = 50) and Nubian goats (n = 50) raised under the same feeding and management conditions (Gu et al., 2022), in the Demonstration Base of Tropical Herbage & Livestock Circular Agriculture (Zhanjiang, China). These goats were used to carry out the Sanger sequencing (Sanger et al., 1977) of the region containing a putative nonsense variant in the DAXX gene. Sequencing of Genomic DNA from 15 MurcianoGranadina Bucks and Bioinformatic Analysis Genomic DNA from the 15 bucks were sequenced at the Centre Nacional d’Anàlisi Genòmica (Barcelona, Spain; https: / / www .cnag .eu). The synthesis of pairedend multiplex libraries was performed using the KAPA HyperPrep kit (Roche, Sant Cugat, Spain), following the instructions of the manufacturer. Libraries were sequenced on a NovaSeq 6000 sequencing platform (Illumina) with a read length of 2 × 150 base pairs. Base calling and quality control analyses were carried out with the RTA v3 analysis pipeline following the guidelines of the manufacturer. Quality-checked filtered reads were mapped to the Capra hircus genome version ARS1 (Bickhart et al., 2017) using the BWA aligner v0.7.17 (Li and Durbin 2009). Aligned reads in Binary Alignment Map (BAM) files were sorted with Samtools v0.1.19 (Li et al., 2009) and duplicated reads were marked with Picard MarkDuplicates (https: / / broadinstitute .github .io/ picard/ ). Joint variant calling was performed with the Genome Analysis Toolkit (McKenna et al., 2010). The resulting variants were annotated with the database ARS1.99 using the SnpEff v5.0 (Cingolani et al., 2012) and EnsemblVEP (McLaren et al., 2016) packages. To extract the intersection of the nonsense variants detected with both methods, we used VCF-compare (https: / / vcftools .github .io/ perl _module .html #VCF -compare). The allelic frequencies of nonsense variants were calculated with VCFtools v4.2 (Danecek et al., 2011). Loci containing multiple potentially deleterious variants or missense Wang et al.: NONSENSE VARIANTS IN GOATS Journal of Dairy Science Vol. 107 No. 12, 2024 11226 variants in the same individual underwent in-depth examination through visual inspection of BAM alignments using the Integrative Genomics Viewer (https: / / igv .org/ ; Thorvaldsdóttir et al., 2013) to check if they are pseudogenes. Variants detected within reads mapping to putative pseudogenes were systematically removed. In addition, we used the multibase codon-associated variant re-annotation (MACARON) software (Khan et al., 2018) to identify, re-annotate, and predict the AA change resulting from multiple single-nucleotide variants within the same genetic codon. We did so to exclude nonsense variants mapping to codons containing other variants disrupting their expected protein truncation effect. Functional Annotation of Genes Containing Nonsense Variants The functional constraint of genes containing nonsense variants was assessed based on a loss-of-function observed/expected upper bound fraction (LOEUF) metric derived from the Genome Aggregation Database (gnomAD v3.1) sequencing data sets. The LOEUF metric (Karczewski et al., 2020) represents a conservative estimate of the ratio of observed to expected loss-of-function variants, and it can range from 0 (gene depleted from loss-of-function variants due to a strong evolutionary constraint) to 9 (gene not depleted from loss-of-function variants because of the absence of evolutionary constraint). Moreover, Gene Ontology (GO) enrichment analysis and the Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were performed with Annotation Visualization and Integration Discovery (https: / / david .ncifcrf .gov/ summary .jsp), gprofiler (Raudvere et al., 2019; https: / / biit .cs .ut .ee/ gprofiler/ gost), and KOBAS (Bu et al., 2021; http: / / kobas .cbi .pku .edu .cn/ anno _iden .php). Data from the MGI database (Baldarelli et al., 2021; http: / / www .informatics .jax .org/ ) and the OMIM (Hamosh et al., 2005; https: / / www .omim .org/ ) were used to infer the potential phenotypic effects of nonsense variants. TaqMan Open Array Genotyping of 20 Nonsense Variants with Predicted Harmful Effects We aimed to select a set of 20 nonsense variants with likely harmful consequences to be genotyped in parentoffspring trios. The potential deleteriousness of variants was predicted according to information from the MGI database (Baldarelli et al., 2021; http: / / www .informatics .jax .org/ ) and OMIM (Hamosh et al., 2005; https: / / www .omim .org/ ) database. Although nonsense variants within 20% of the 3′ end coding region are expected to be less harmful than those outside, we did not exclude such variants from our analysis because this criterium would have been too restrictive, leading to the exclusion of many variants of potential interest. The existence of multiple mRNA isoforms might also limit the effect of a nonsense variant (due to the existence of isoforms lacking the exon carrying the nonsense variant). However, we did not consider this criterium because deleteriousness can be strongly influenced by the expression levels of each isoform, a parameter that is currently unknown. Finally, each candidate variant was manually annotated and those that do not introduce stop codons were eliminated from the data set. The upstream and downstream flanking sequences (100 bp) of the 20 selected variants were submitted to the Custom TaqMan Assay Design Tool website (https: / / www5 .appliedbiosystems .com/ tools/ cadt, Life Technologies) to check if they could be typed in a TaqMan Open Array multiplex assay platform. Selected variants were genotyped in dams (n = 284) mated with each one of the 6 bucks and their offspring (n = 302). Genotyping was performed at the Servei Veterinari de Genética Molecular of the Universitat Autònoma de Barcelona (http: / / sct .uab .cat/ svgm/ en) by using a QuantStudio 12K Flex Real-Time PCR System (Thermo Fisher Scientific). Analysis and visualization of the results were carried out in a dedicated Thermo Fisher online platform (http: / / apps .thermofisher .com/ ). The Hardy-Weinberg equilibrium test was performed for each variant by using the Science Primer platform (https: / / scienceprimer .com/ hardy -weinberg -equilibrium -calculator). Sanger Sequencing of the Region Containing a Nonsense Variant in the DAXX Gene The nonsense XM_018038577.1:c.1134C > G variant associated with the goat DAXX gene was the only one showing a significant (P < 0.00001) depletion of homozygous genotypes for the G variant. When designing primers to amplify the region containing this variant and blasting them against the goat genome, we became aware that the region containing this variant is duplicated. So, we designed 2 pairs of primers (Supplemental Figure S1, see Notes) amplifying the chromosome-23 region containing the XM_018038577.1:c.1134C > G variant (represented by GenBank sequence XM_018038577.1) and the duplicated region on chromosome 6 (represented by GenBank sequence XM_018049164.1), which harbors a death domain-associated protein 6-like pseudogene (LOC108636199). Sequences of the primers were F_5′ CCC TCG TTC GTC AGG TGT GT-3′ and R_5′-AGG GAT CTG TTG CCC CGT TC-3′ for the partial amplification of the DAXX gene on chromosome 23, and F_5′-AGA TGA CTA TAG GCC AGG CAT TG-3′ and R_5′-TCC AGC GAG GAT TTA GGG GT-3′ for the partial amplification of the LOC108636199 locus on chromosome 6 (Supplemental Figure S1). Wang et al.: NONSENSE VARIANTS IN GOATS 11227 Journal of Dairy Science Vol. 107 No. 12, 2024 Regions flanked by these primers were amplified in a set of individuals comprising Shaanbei White Cashmere goats (n = 50), Guanzhong goats (n = 50), Leizhou black goats (n = 50), Hezhang Black goats (n = 50), and Nubian goats (n = 50). Genomic DNA was isolated from ear-tissue samples using the Genomic DNA Extraction Kit for animal tissues or cells (Solarbio, Beijing, China). The volume of the amplification reaction was 25 μL, and it contained 1.0 μL of genomic DNA (20 ng/μL), 1.0 μL of each primer (10 μM), and 12.5 μL of 2× Taq Master mix, which included 1.5 mM MgCl2, 200 μM dNTP, and 0.625 units of Taq DNA polymerase (BioLinker, Shanghai, China) The thermal profile was 94°C for 5 min; followed by 35 cycles of 30 s at 94°C, 30 s at 55°C, and 1 min at 72°C, plus a final extension step at 72°C for 5 min. Sanger sequencing of amplicons corresponding to these 2 regions in 250 individuals was performed by the AuGCT Company (Beijing, China) by using an Applied Biosystems 3730XL DNA sequencer (AuGCT, Beijing, China). RESULTS Identification of Nonsense Variants and Their Genomic Distribution The statistics obtained in the whole-genome sequencing experiment are shown in Supplemental Table S1 (see Notes). The average coverage of the 15 goat genomes was 32.92×, and the average percentage of uniquely mapped reads was 90.13%. With SnpEff (Cingolani et al., 2012) and Ensembl-VEP (McLaren et al., 2016), we detected 991 and 1,463 nonsense variants segregating in the genomes of the 15 Murciano-Granadina bucks, respectively. By using VCF-compare, we identified 947 nonsense variants consistently detected with both methods. After quality control and filtering (−max-missing 0.90 −minQ 30 −minDP 10), we retained 866 nonsense variants. From these, we filtered out 26 multinucleotide variants disrupting premature stop codons with MACARON (Khan et al., 2018). The allelic frequencies of nonsense variants were in general quite low (Figure 1A), because ~28% of them had frequencies below 0.05. Regarding the genomic distribution of nonsense variants, we observed that they are quite scattered along the protein-coding region, being particularly enriched in its 3′ end (Figure 1B). Functions and Biological Essentiality of Genes Containing Nonsense Variants Genes containing nonsense variants were enriched in GO terms related to the integral component of the membrane (P = 1.4E-15), the G-protein coupled receptor activity (P = 6.4E-11), and the olfactory receptor activity (Table 1; P = 1.5E-10). We also found a significant enrichment for the ATPase activity coupled to transmembrane movement of substances pathway (P = 0.0029), as the remaining pathways, including defense response to the virus (P = 0.0272), hormonal activity (P = 0.0294), and sensory perception of taste (P = 0.0417), were less significant. The KEGG pathway analysis also showed that, by far, olfactory transduction is the most significant pathway (P = 3.51E-29), followed by herpes simplex virus 1 infection (P = 1.61E-10) and antigen processing and presentation (Table 1; P = 3.34E-5). We also calculated the LOEUF scores of genes containing nonsense variants (Figure 1C). As expected, most of these genes had LOEUF scores of 5 or larger indicating that they are not under strong functional constraint, a result that is concordant with the functional enrichment in genes related to olfactory transduction. Genotyping of Variants Predicted to Have Harmful Effects in Parent-Offspring Trios By using the OMIM (Hamosh et al., 2005) and MGI (Baldarelli et al., 2021) databases, we identified 23 and 51 genes, respectively, containing nonsense variants and with LOEUF scores <6. The inactivation of several of these loci is expected to have harmful consequences including, but not limited to, lethality (Table 2). Based on such information plus bibliographic references, we initially selected 20 nonsense variants mapping to 19 genes (Table 2). Manual annotations of these variants (Supplemental Figure S2, see Notes) revealed that 3 of them do not introduce stop codons (Supplemental Figure S2), so they were removed from our data set. The remaining 17 variants were genotyped with an Open Array system in buck-dam-offspring trios formed by the 284 dams mated to the 6 bucks and 302 offspring (Table 2). After completing genotyping tasks, 4 nonsense variants in the ACVR2B, RNLS, NPC1, and OTUD6B genes did not yield valid results. As shown in Supplemental Figure S2, we observed that 12 out of the 13 successfully genotyped variants segregated in the genomes of a panel of 250 European goats (10 breeds), 249 African goats (4 breeds), and 109 Asian goats (4 breeds) that were sequenced in the VarGoats project (Denoyelle et al., 2021). Full details regarding how this whole-genome sequence data set was built can be found in Denoyelle et al. (2021). The only variant that did not segregate in any of these 608 individuals was XP_017894826.1:p.Arg204* in the BCL2 gene (Supplemental Figure S2), and the same variant had an extremely low frequency among the genotyped Murciano-Granadina individuals (Table 2). The remaining variants segregated in the investigated 608 goats from Europe, Africa, and Asia although, in general, at a minimum allele frequency (MAF) lower than 10%. ExWang et al.: NONSENSE VARIANTS IN GOATS Journal of Dairy Science Vol. 107 No. 12, 2024 11228 Wang et al.: NONSENSE VARIANTS IN GOATS Figure 1. (A) Allele frequencies of nonsense variants detected in the genomes of 15 Murciano-Granadina bucks. (B) Chromosome distribution (heatmap and absolute counts per chromosome) and location along the protein-coding region of nonsense variants segregating in the genomes of 15 Murciano-Granadina bucks. (C) Absolute counts of genes containing nonsense variants that have been classified according to their LOEUF score. A LOEUF score of 0 indicates that a gene is under strong selective constraint and has an essential function, whereas a LOEUF score of 9 corresponds to a gene for which inactivation is well tolerated. CDS = coding sequence. 11229 Journal of Dairy Science Vol. 107 No. 12, 2024 ceptions to this rule were NANOG NP_001272505.1:p. Gln273* in African (MAF = 0.205) and Asian (MAF = 0.1225) goats, MAP9 XP_017917128.1:p.Arg611* in Asian goats (MAF = 0.1875), PDE4D XP_005694732.3:p. Cys50* in European goats (MAF = 0.357), and DAXX XP_017894066.1:p.Tyr378* in European (MAF = 0.16) and African (MAF = 0.115) goats. We did not observe a significant absence of homozygous genotypes for several nonsense variants located in genes whose inactivation causes lethality in knockout mice (Table 3), such as MECR (XM_018057207.1:c.1152C > G, with the G-allele introducing a premature stop codon), BRCA2 (XM_018056629.1:c.8754T > A, with the A-allele introducing a premature stop codon), and NANOG (NM_001285576.1:c.817C > T, with the Tallele introducing a premature stop codon). For the XM_018063707.1: c.5419C > T (with the T-allele introducing a premature stop codon) variant in the TNRC6C gene, which is also a lethal gene in mice, only one individual displayed a TT genotype (Table 3). When checking the annotation of these 4 variants, we verified that all of them are located within the 10% 3′ end of the protein-coding region (Supplemental Figure S2). Nonsense variants in the PPA2, MAP9, and SLC25A26 genes also mapped to such region (Supplemental Figure S2). In Supplemental Figure S3 (see Notes), we show the repertoire of mRNA isoforms for genes with at least 2 distinct mRNA isoforms. In the case of the MECR and TNRC6C genes, the nonsense variants only affected one (MECR: mRNA X2, TNRC6C: mRNA X3) of the 6 (MECR) or 7 (TNRC6C) reported isoforms (Supplemental Figure S3). Moreover, the 2 affected isoforms showed a highly divergent sequence when compared with their counterparts. This isoform-specificity of nonsense variants was also observed for the BCL2, MAP9, PPA2, and SLC25A26 genes (Supplemental Figure S3). In addition, 6 variants in the NBAS, MAP9, SLC25A26, DAXX, and BCL2 genes showed a complete absence of individuals homozygous for the nonsense variant (Table 3), but such depletion was only significant (P < 0.00001) for the XM_018038577.1: c.1134C > G variant in the DAXX gene. Visual inspection of DAXX genotype clustering with the Thermo Fisher online platform (http: / / apps .thermofisher .com/ ) showed that the cluster of heterozygotes was quite scattered (Supplemental Figure S4, see Notes), suggesting that such cluster might be the result of the overlap of CG and GG clusters. Because the DAXX XM_018038577.1:c.1134C > G variant was a promising candidate variant to have lethal effects and the Open Array (Thermo Fisher) genotyping results were doubtful due to the scattered distribution of the heterozygous genotype, we decided to re-genotype this variant by Sanger sequencing. When designing primers, we became aware that the chromosome 23: 41192961–41193179 bp region containing the XM_018038577.1:c.1134C > G variant is duplicated on Wang et al.: NONSENSE VARIANTS IN GOATS Table 1. Gene ontology terms and KEGG pathways significantly (corrected P-value <0.05) overrepresented among genes containing predicted nonsense variants segregating in the genomes of 15 Murciano-Granadina bucks Term No. of genes P-value Corrected P-value1 Gene ontology terms GO: 0016021 ~integral component of membrane 144 2.67E-16 1.4E-15 GO: 0004930 ~G -protein coupled receptor activity 47 6.05E-12 6.4E-11 GO: 0004984 ~olfactory receptor activity 42 2.20E-11 1.5E-10 GO: 0042626 ~ATPase activity, coupled to transmembrane movement of substances 7 0.0011707 0.0029473 GO: 0015914 ~phospholipid transport 3 0.0120228 0.0148164 GO: 0006355 ~regulation of transcription, DNA-templated 12 0.0125267 0.0176249 GO: 0051607 ~defense response to virus 7 0.0127781 0.0272945 GO: 0005179 ~hormone activity 6 0.0167250 0.0294567 GO: 0030515 ~snoRNA binding 3 0.0168168 0.0294613 GO: 0022857 ~transmembrane transporter activity 7 0.0188329 0.0322914 GO: 0005886 ~plasma membrane 48 0.0238185 0.0341911 GO: 0050909 ~sensory perception of taste 3 0.0262560 0.0417436 GO: 0005044 ~scavenger receptor activity 50.0321463 0.0486462 KEGG pathways Olfactory transduction 105 1.56E-31 3.51E-29 Herpes simplex virus 1 infection 42 1.43E-12 1.61E-10 Antigen processing and presentation 13 4.46E-7 3.34E-5 Autoimmune thyroid disease 11 2.75E-6 1.55E-4 Graft-versus-host disease 9 1.29E-5 5.79E-4 Allograft rejection 9 2.22E-5 8.33E-4 Type I diabetes mellitus 9 5.24E-5 1.68E-3 Viral myocarditis 8 7.75E-4 2.00E-2 Complement and coagulation cascades 9 8.01E-4 2.00E-2 ABC transporters 8 1.16E-3 2.60E-2 1P-value corrected for multiple testing. Journal of Dairy Science Vol. 107 No. 12, 2024 11230 Wang et al.: NONSENSE VARIANTS IN GOATS Table 2. Location and main features of 20 predicted nonsense variants with potential harmful effects on biological viability CHIR1Position Gene symbol Phenotype AA and nucleotide replacement 1 79618036 RTP4 Homozygous knockout reduces susceptibility to parasitic and viral infection2 XP_005675197.2: p.Arg257* XM_005675140.3: c.769C > T 2 11614286 MECR Homozygous variants in this gene lead to embryonic lethality due to pleiotropic effects including abnormalities in fetal placental and cardiovascular development and embryonic growth retardation2 XP_017912696.1: p.Tyr384* XM_018057207.1: c.1152C > G 3 62127712 MCOLN3 Homozygotes are white with small patches of color and show severe behavioral abnormalities, poor postnatal viability and are almost infertile2 XP_017901176.1: p.Trp186* XM_018045687.1: c.557G > A 5 27520597 NR4A1 NR4A1 and NR4A3 abrogation lead to rapidly lethal acute myeloid leukemia in mice3 QHW05400.1: p.Arg588*6 MN197544: c.1762C > T6 6 20144367 PPA2 PP2A mutations may cause a spectrum of mitochondrial disease phenotypes; severe symptoms include seizures, lactic acidosis, cardiac arrhythmia, and death within days of birth4 XP_017904735.1: p.Arg327* XM_018049246.1: c.979C > T 11 82481680 NBAS The NBAS gene is recurrently mutated in pediatric hemophagocytic lymphohistiocytosis5 XP_017911109.1: p.Gln105* XM_018055620.1: c.313C > T 12 58071934 BRCA2 Homozygous null mutants are embryonic lethal with abnormalities including growth retardation, neural tube defects, and mesoderm abnormalities2 XP_017912118.1: p.Tyr2918* XM_018056629.1: c.8754T > A 14 81057228 PLEC Targeted variants of this gene result in neonatal death, skin blistering, impaired myofibril integrity, reduced hemidesmosome number, and disintegration of intercalated disks in the heart2 XP_017914185.1: p.Trp3754*6 XM_018058696.1: c.11261G > A6 14 81058774 PLEC XP_017914185.1: p.Arg3239*6 XM_018058696.1: c.9715A > T6 14 8887926 OTUD6B Mice homozygous for a knockout allele exhibit complete perinatal lethality, decreased fetal size, and ventricular septal defects2 XP_005689305.1: p.Trp10* XM_005689248.3: c.30G > A 15 81546465 NANOG Mice homozygous for a disruption in this gene die between E3.5 and E5.5 with abnormal embryonic and extraembryonic tissue development2 NP_001272505.1: p.Gln273* NM_001285576.1: c.817C > T 17 70604449 MAP9 Mice homozygous or heterozygous for a null allele exhibit increased aneuploidy, chromosomal instability, detached centrosomes, and pericentriolar material fragmentation in mouse embryonic fibroblasts2 XP_017917129.1: p.Arg611* XM_018061639.1: c.1831C > T 19 53253007 TNRC6C Mice homozygous for a gene trap allele exhibit complete neonatal lethality with cyanosis, respiratory distress, and thickened mesenchyme in air sacs2 XP_017919196.1: p.Arg1807* XM_018063707.1: c.5419C > T 20 20004713 PDE4D Homozygotes for targeted null variants exhibit delayed growth, female infertility associated with impaired ovulation, and reduced postnatal viability2 XP_005694732.3: p.Cys50* XM_005694675.3: c.150C > A 22 11721225 ACVR2B Mice homozygous for targeted variants that inactivate the gene show abnormal lateral asymmetry and homeotic transformation of the axial skeleton, and die shortly after birth with extensive cardiac defects2 XP_017922111.1: p.Lys213* XM_018066622.1: c.637A > T Continued 11231 Journal of Dairy Science Vol. 107 No. 12, 2024 chromosome 6: 3909498–3909716 bp (Figure 2). In consequence, we designed primers (Supplemental Figure S1) to specifically amplify each one of the 2 copies of the duplicated region. Two amplicons of 433 bp (chromosome 23) and 222 bp (chromosome 6) were sequenced with the Sanger method. Comparison of the sequences (Figure 2A) showed that the putative XM_018038577.1:c.1134 C > G variant in the DAXX gene (chromosome 23) is, in reality, a dinucleotide variant mapping to the LOC108636199 locus on chromosome 6 (XM_018049164.1:c.1270–1271 CA > GG variant). Indeed, the 2 orthologous sites in the DAXX gene are monomorphic. In consequence, for the DAXX gene only one genotype (CC) is possible for the putative XM_018038577.1:c.1134C > G variant on chromosome 23, whereas 3 genotypes (CC, CG, and GG) could eventually be found in the LOC108636199 locus for position 1270 of the XM_018049164.1:c.1270–1271 CA > GG variant on chromosome 6 targeted in the TaqMan assay (Figure 2B). These 2 regions have very similar nucleotide sequences (Figure 2A), so when they coamplify the resulting genotypes for the site targeted by the TaqMan assay can be exclusively CG (chromosome 23: CC + chromosome 6: CG or GG) or CC (chromosome 23: CC + chromosome 6: CC). This results in an apparent depletion of GG individuals. Indeed, no GG individuals were identified in the 250 Shaanbei White Cashmere, Guanzhong, Leizhou Black, Hezhang Black, and Nubian individuals for which the chromosome 23 region containing the nonsense variant was amplified with specific primers and sequenced with the Sanger method. DISCUSSION The Genomic Landscape of Nonsense Variants in 15 Capra hircus Genomes The sequencing of 15 Murciano-Granadina bucks with good coverage (32.9×) allowed us to identify 947 nonsense variants. Sequencing of 600 bovine exomes revealed the segregation of 1,377 nonsense variants (Charlier et al., 2016). In humans, Pelak et al. (2010) sequenced 20 genomes with a coverage similar to ours and detected 554 nonsense variants. Estimates obtained by the 1000 Genomes Project Consortium (2010, 2015) Wang et al.: NONSENSE VARIANTS IN GOATS CHIR1Position Gene symbol Phenotype AA and nucleotide replacement 22 35072050 SLC25A26 Embryos homozygous for a transposon insertion appear growth retarded and underdeveloped and die after E8.5 but before birth2 XP_005695813.1: p.Gln284* XM_005695756.3: c.850C > T 23 41193060 DAXX Mice homozygous for a targeted variant of this gene display extensive apoptosis and embryonic lethality2 XP_017894066.1: p.Tyr378* XM_018038577.1: c.1134C > G 24 33418088 NPC1 Homozygotes for spontaneous and chemically induced variants may exhibit lysosomal storage of nonesterified cholesterol, neurodegeneration, ataxia, presence of foam cells, sterility, and shortened lifespan2 XP_017894966.1: p.Arg1264* XM_018039477.1: c.3790C > T 24 61676529 BCL2 Homozygous null mutants show pleiotropic abnormalities including small size, increased postnatal mortality, polycystic kidneys, apoptotic involution of thymus and spleen, graying in the second hair follicle cycle, and reduced numbers of motor, sympathetic, and sensory neurons2 XP_017894826.1: p.Arg204* XM_018039337.1: c.610C > T 26 41247997 RNLS Mice homozygous for a knockout allele exhibit aggravated ischemic myocardial damage, increased heart rate, increased blood pressure, and increased serum levels of dopamine, adrenaline, and noradrenaline2 XP_017897022.1: p.Gln85* XM_018041533.1: c.253C > T 1CHIR = Capra hircus chromosome. 2MGI Mouse database. 3Mullican et al., 2007. 4Kennedy et al., 2016. 5Bi et al., 2022. 6Variants marked in bold are not correctly annotated and do not introduce any premature stop codon. Table 2 (Continued). Location and main features of 20 predicted nonsense variants with potential harmful effects on biological viability Journal of Dairy Science Vol. 107 No. 12, 2024 11232 indicate that each individual typically differs from the reference human genome sequence by 80 to 200 nonsense variants, and in pigs such number is ~30 (Groenen et al., 2012). In general, the ability and accuracy to detect nonsense variants depend strongly on parameters, such as population size and composition, sequencing technology and coverage, variant-calling algorithms, filtering criteria, and the availability of a high-quality and well-annotated genome assembly. Usually, these parameters differ across studies making it difficult to establish meaningful comparisons between them (MacArthur et al., 2012). Another important factor influencing the abundance of deleterious variants is effective population size (Charlier et al., 2016). Indeed, the increase of effective population size leads to an increment of the average number of embryonic lethal variants harbored per individual (Charlier et al., 2016). Oliveira et al. (2016) reported that the effective number of founders of the Murciano-Granadina breed encompassed 967 individuals, whereas the most recent estimate (from 2006, same authors) of effective population size is ~100 individuals. According to simulations reported by Charlier et al. (2016), this would correspond to 0.53 ± 0.21 recessive lethals carried on average per individual. We have observed that approximately one-fourth of the 947 nonsense variants segregating in the 15 buck genomes have frequencies below 5%. Extensive sequencing of human genomes revealed that 76% of the total variation has frequencies below 0.5%, and only 9.5% of the variants corresponded to common alleles (i.e., those with frequencies above 5%; 1000 Genomes Project Consortium, 2015). In addition, nonsense variants with harmful effects are expected to have lower frequencies than neutral alleles because they are purged by purifying selection. Our findings might be explained by the fact that we have sequenced a reduced set of individuals, making it very difficult to detect alleles with low frequencies. Indeed, the 1000 Genomes Project Consortium (2015) reported that the majority of variants identified in a single genome are common (only 1%–4% of the variants have a frequency below 0.5%), thus implying that common alleles (frequency ≥5%) are preferentially detected, over those with frequencies below 5%, when sequencing a small number of individuals. Wang et al.: NONSENSE VARIANTS IN GOATS Table 3. Absolute frequencies1 and departure from the Hardy-Weinberg equilibrium of 13 predicted nonsense variants successfully genotyped in dams mated with each one of the bucks and their offspring (in bold: nonsense variants never found in homozygosity)2 CHIR Gene Variant Buck3 Dam Offspring Frequency HW P-value Frequency HW P-value Frequency HW P-value 1 RTP4 XP_005675197.2: p.Arg257* 3/2/1 0.25 202/51/9 0.06 171/85/21 0.09 (8/6/1) (0.28) 2 MECR XP_017912696.1: p.Tyr384* 3/3/0 0.25 232/35/2 0.87 169/70/4 0.56 (10/5/0) (0.44) 3 MCOLN3 XP_017901176.1: p.Trp186* 4/2/0 0.44 165/88/23 0.09 167/108/11 0.45 (8/6/1) (0.28) 6 PPA2 XP_017904735.1: p.Arg327* 5/0/1 0.69 153/61/14 0.08 131/77/9 0.72 (10/3/2) (0.44) 11 NBAS XP_017911109.1 p.Gln105* 5/1/0 0.69 274/1/0 1.00 264/21/0 0.81 (14/1/0) (0.87) 12 BRCA2 XP_017912118.1: p.Tyr2918* 4/2/0 0.44 194/70/9 0.69 199/72/13 0.17 (9/5/1) (0.36) 15 NANOG NP_001272505.1: p.Gln273* 5/1/0 0.69 130/101/25 0.71 148/99/12 0.67 (13/2/0) (0.75) 17 MAP9 XP_0179171289.1: p.Arg611* 4/2/0 0.44 249/13/0 0.92 219/39/0 0.42 (11/4/0) (0.54) 19 TNRC6C XP_017919196.1: p.Arg1807* 5/1/0 0.69 254/20/1 0.68 259/26/0 0.72 (14/1/0) (0.87) 20 PDE4D XP_005694732.3: p.Cys50* 3/3/0 0.25 198/71/4 0.70 200/78/8 0.99 (11/4/0) (0.54) 22 SLC25A26 XP_005695813.1: p.Gln284* 5/1/0 0.69 262/10/0 0.95 252/28/0 0.68 (13/2/0) (0.75) 23 DAXX XP_017894066.1: p.Tyr378* 5/1/0 0.69 103/171/0 0.00 146/138/0 0.00 (11/4/0) (0.54) 24 BCL2 XP_017894826.1: p.Arg204* 5/1/0 0.69 246/1/0 1.00 237/0/0 1.00 (13/2/0) (0.75) 1Absolute frequencies are ordered as follows: homozygous individuals for the wild-type allele or heterozygous individuals for the nonsense variant or homozygous individuals for the nonsense variant. 2CHIR = Capra hircus chromosome; HW P-value = statistical significance of the departure from the Hardy-Weinberg equilibrium. 3Data from the 6 sequenced bucks mated to dams and, between parentheses, from the 15 sequenced bucks (please see the “Animal Material” section for more details).