A robust sequencing assay of a thousand amplicons for the high-throughput population monitoring of Alpine ibex immunogenetics.

Polymorphism for immune functions can explain significant variation in health and reproductive success within species. Drastic loss in genetic diversity at such loci constitutes an extinction risk and should be monitored in species of conservation concern. However, effective implementations of genome-wide immune polymorphism sets into high-throughput genotyping assays are scarce. Here, we report the design and validation of a microfluidics-based amplicon sequencing assay to comprehensively capture genetic variation in Alpine ibex (Capra ibex). This species represents one of the most successful large mammal restorations recovering from a severely depressed census size and a massive loss in diversity at the major histocompatibility complex (MHC). We analysed 65 whole-genome sequencing sets of the Alpine ibex and related species to select the most representative markers and to prevent primer binding failures. In total, we designed ~1,000 amplicons densely covering the MHC, further immunity-related genes as well as randomly selected genome-wide markers for the assessment of neutral population structure. Our analysis of 158 individuals shows that the genome-wide markers perform equally well at resolving population structure as RAD-sequencing or low-coverage genome sequencing data sets. Immunity-related loci show unexpectedly high degrees of genetic differentiation within the species. Such information can now be used to define highly targeted individual translocations. Our design strategy can be realistically implemented into genetic surveys of a large range of species. In conclusion, leveraging whole-genome sequencing data sets to design targeted amplicon assays allows the simultaneous monitoring of multiple genetic risk factors and can be translated into species conservation recommendations.

INTRODUCTION

Biodiversity is currently undergoing a dramatic decline caused by ecological and anthropological pressures (Barnosky et al., 2011; Ceballos et al., 2015, 2017; WWF, 2018). Endangered species are particularly prone to genetic risks due to population bottlenecks and habitat fragmentation that lead to: genetic diversity loss, inbreeding, introgression from related species and accumulation of deleterious mutations (Allendorf et al., 2010; Frankham, 2005). Genetic diversity is generally thought to be declining in wild populations (Leigh et al., 2019). Low genetic diversity can affect fitness and survival (Reed & Frankham, 2003), as shown for instance in cheetah (Acinonyx jubatus, O’Brien et al., 1983; O’Brien et al., 1985), Florida panthers (Puma concolor coryi, Pimm et al., 2006; Roelke et al., 1993) and Alpine ibex (Capra ibex, Brambilla et al., 2015). Loss of genetic diversity at adaptive immune loci is particularly problematic because populations will lack genetic variants conferring disease resistance, and this has been shown to be problematic in amphibians (Kosch et al., 2019; Savage et al., 2011), Tasmanian devils (Sarcophilus harrisii, Siddle et al., 2007), giant panda (Ailuropoda melanoleuca, Zhu et al., 2020) and Alpine ibex (Brambilla et al., 2018). However, despite the importance for conservation management, the genetic underpinnings of disease susceptibility are largely unknown in most nonmodel species (Schoville et al., 2012). The major obstacles in investigating wildlife immunogenetics are the challenge of collecting quantitatively and qualitatively adequate phenotypic data on diseases in wild populations, as well as the lack of genetic tools suitable in the conservation framework (Holderegger et al., 2019). Hence, there is a pressing need to establish effective genetic monitoring tools for many wild species (Acevedo‐Whitehouse & Cunningham, 2006; Allendorf et al., 2010).

The major histocompatibility complex (MHC), a highly polymorphic region involved in foreign antigen recognition, is an important player for disease susceptibility in vertebrates. Co‐evolution with pathogens often causes balancing selection on MHC polymorphism, maintaining high genetic diversity including deeply divergent alleles. Genetic variants under selection often encode the ability to present a wider range of antigens to T‐cells and, thus, to recognize a greater variety of parasites. Consequently, heterozygotes often have elevated resistance and are favoured by selection (Bernatchez & Landry, 2003). Species that have undergone a strong bottleneck such as Alpine and Iberian ibex (Capra pyrenaica), Tasmanian devils, cheetahs or Galapagos penguins (Spheniscus mendiculus) show strongly reduced genetic diversity at the MHC compared to related species (Angelone et al., 2018; Bollmer et al., 2007; Brambilla et al., 2018). Low genetic diversity in cheetahs is thought to have contributed to high mortality rates (O’Brien et al., 1983, 2017). High prevalence of a fatal facial cancer in Tasmanian devil is thought to, at least, partially stem from low MHC diversity (Siddle et al., 2007). Low levels of MHC diversity are consequently a major threat to endangered species. Alongside with the adaptive immune system, the innate immune system plays an important role in the defence against a wide range of pathogens. Genetic variation within some genes, such as Toll‐like receptors (TLRs), has been shown to be linked to variation in immune competence (Ammerdorffer et al., 2014; Tschirren et al., 2013). Yet the field of wildlife immunogenetics is only just emerging and immune‐related genes outside of the MHC remain heavily understudied in wild species. The major challenge for monitoring nonmodel wild species is that inferring causal disease susceptibility from other species is only rarely possible and associations need to be established for each species–pathogen interaction (Acevedo‐Whitehouse & Cunningham, 2006).

To make informed conservation management decisions, the genetic health of populations should be taken in account (Allendorf et al., 2010; Biebach et al., 2016 ). Yet integrating information about immune‐related polymorphisms and genetic diversity across the genome remains a major challenge. Conservation genetics studies have traditionally relied on microsatellite markers (Allendorf et al., 2010; Brambilla et al., 2015; Ouborg et al., 2010; Witzenberger & Hochkirch, 2014), which are inexpensive and relatively easy to use. However, only in exceptional cases do microsatellites reveal adaptive genetic variation and are thus not applicable to most relevant immune‐related polymorphism (Allendorf et al., 2010; Steiner et al., 2013). Similarly, conservation genomics methods based on restriction site‐associated DNA sequencing (RAD; e.g., ddRAD, GBS, Davey et al., 2011) intrinsically lack the ability to target particular polymorphisms of interest (Garvin et al., 2010). Untargeted single nucleotide polymorphism (SNP) genotyping can also be achieved through the capture of RAD tags as in rapture (Ali et al., 2016; Komoroske et al., 2019; Sard et al., 2020) and radcap (Hoffberg et al., 2016). While the development and sample preparation time to capture RAD tags is short, the possibility to target specific loci (e.g., putative adaptive variation) is highly limited (reviewed in Meek & Larson, 2019). Whole genome sequencing resolves the issue of marker selection but is expensive, requires high‐quality input material to work reliably and imposes significant computational challenges due to the very large amounts of data (Davey et al., 2011). This is because most genome sequencing techniques are still relatively costly, and require bioinformatic expertise and powerful computing systems. These are major impediments for the successful implementation in conservation efforts (Shafer et al., 2015; Taylor et al., 2017).

Targeted SNP‐based genotyping overcomes many of the constraints of the above‐mentioned methodologies, because the approach typically tolerates a wide range of input material and is cost‐effective (Eriksson et al., 2020; Kleinman‐Ruiz et al., 2017; Natesh et al., 2019; Wright et al., 2015). Early applications of targeted SNP genotyping methods relied on hybridization arrays and fluorescent probes. Applications such as TaqMan, Sequenom and SNaPshot arrays have enabled the analyses of tens to hundreds of SNPs in hundreds of samples (reviewed in Syvänen, 2001). The targeted sequencing of amplicons evolved from array‐based systems and takes advantage of the large amounts of sequencing data produced by next‐generation sequencers. Nonetheless, two important limitations apply to existing targeted sequencing approaches. First, unknown polymorphisms in flanking sequences of the targeted SNPs can lead to amplification dropouts. This shortcoming can be remedied by integrating whole‐genome sequencing information from a limited number of individuals. Second, the parallel amplification of a large number of loci requires optimization and normalization steps to ensure balanced amounts of PCR product for each locus‐by‐individual combination (Meek & Larson, 2019).

Amplicon sequencing by genotyping‐in‐thousands (GT‐seq) offers a highly flexible and customizable approach to analyse hundreds of markers (typically 200–500; Campbell et al., 2015; reviewed in Meek & Larson, 2019). GT‐seq saw broad applications across systems including in large‐scale salmonid surveys (Matala et al., 2016; May et al., 2020). However, the development of a new GT‐seq assay requires significant efforts to optimize multiplex amplifications. A number of similar methods for amplicon sequencing have also been developed (e.g., MTA‐seq, three‐round multiplex, Onda et al. 2018; Chen et al. 2016). Some alternatives to GT‐seq rely partially on commercial solutions, including AmpliSeq, TruSeq Custom Amplicon and the microfluidics‐based Juno system. Microfluidics technology is notable because it allows for parallel multiplex amplifications and automated sample–primer pool pairing in miniaturized wells. This lowers overall reaction volumes, ensures accurate mixing and reaction conditions, and reduces primer incompatibility issues. These factors stabilize amplification yields across amplicons and, as a consequence, increase the total number of usable amplicons per assay (Chiu et al., 2017; Hess et al., 2020).

Here, we develop and validate a high‐throughput amplicon sequencing assay based on a microfluidics platform, which enables the simultaneous monitoring of population structure, immunogenetics and introgression from a domestic relative in Alpine ibex. The species is an excellent model to assess the broad applicability of novel genotyping assays for conservation. Because of a dramatic recent bottleneck caused by a near extinction event in the 19th century (Grodinsky & Stüwe, 1987), the census size was reduced to fewer than 100 individuals limited to one single population in northern Italy (Gran Paradiso). A captive breeding programme restored the species to a current census size of over 50,000 individuals spread across the European Alps (Brambilla et al., 2020; Toïgo et al., 2020). However, the historical bottleneck left substantial genome‐wide signatures of low heterozygosity, in particular at MHC loci, which may threaten long‐term species viability (Grossen et al., 2014, 2018, 2020). Hybridization events with domestic goat (Capra hircus) produced introgression tracts at the MHC re‐establishing some genetic variation, which may have been lost due to the species bottleneck (Grossen et al., 2014). Recently, concerns were raised over population declines as direct and indirect consequences of epizootic disease outbreaks (e.g., sarcoptic mange, respiratory diseases, infectious keratoconjunctivitis, brucellosis). A notable brucellosis outbreak was recently observed in the area of Bargy in France. Cases of brucellosis reported in cattle and humans in the region were allegedly traced back to the ibex population, so that French authorities undertook a massive eradication programme leading to the culling of more than 250 individuals (44% of the estimated population) within 2 years following the outbreak (Mick et al., 2014; Quéméré et al., 2020). In summary, Alpine ibex represent a model species where high‐throughput monitoring of inbreeding levels, potential hybridization events and immunogenetic diversity can substantially improve genetic analyses and conservation management efforts.

We established a high‐throughput assay of nearly 1,000 loci covering largely genome‐wide (putatively neutral) polymorphisms, variants found at the MHC and other immune‐related loci, as well as diagnostic variants useful to detect recent hybridization events with domestic goat. We used 51 whole‐genome sequences of Alpine ibex and domestic goat to identify relevant SNPs. Furthermore, we used whole genome sequences to mask polymorphism at primer sites to maximize amplification success. Based on highly parallel Illumina amplicon sequencing of 172 Alpine ibex, Iberian ibex and domestic goat samples, we assessed the accuracy and robustness of the assay across variable input sample quality. Finally, we compared the high‐throughput assay with RAD‐seq and low‐coverage whole‐genome sequencing data sets on the same populations.

MATERIALS AND METHODS

Collection of samples and DNA extraction

Genotyping was performed on 158 Alpine ibex samples, representing eight populations, each with 19–20 sampled individuals. The samples also included an individual suspected to be an albino. Sample material consisted either of tissue or blood. Samples were collected by biopsy darting or during captures carried out in the framework of the EU Interreg‐Alcotra 1664 LEMED‐ibex project. Further samples were previously collected during captures or legal hunting (Biebach & Keller, 2009; Brambilla et al., 2015; Willisch et al., 2012). In addition to the 158 samples described above, we evaluated the performance of four DNA samples collected from faeces of Alpine ibex in the Gran Paradiso National Park (Italy). The faecal samples originated from two mother–kid pairs sampled on the same day. Visual identification and microsatellite analyses helped confirm that the samples were from different individuals (Willisch et al., 2012). No replicates were performed for these samples in our study. We added six further tissue samples from suspected hybrids (based on field observations), five tissue samples from domestic goats and three from Iberian ibex. Detailed information about the origin, collection method and sampling year of each sample is provided in Table S1. Tissue samples were stored in 95%–100% ethanol at room temperature. Blood samples were collected in EDTA Vacutainer tubes and stored at −20℃. Faecal samples were collected opportunistically immediately after defecation, frozen within 3–4 hr and stored at −20℃ until analysis. DNA extractions were carried out using the DNeasy Blood & Tissue and the QIAmp DNA Stool kits (Qiagen) following the manufacturer's protocol.

DNA concentrations of all samples were assessed using a Qubit fluorometer assay (ThermoFisher). The manufacturer's protocol recommends DNA concentrations of 100–200 ng µl−1 for downstream applications. Because these concentrations could not be reached for most of the samples, we chose an ad hoc concentration of 50 ng µl−1 for normalization of all samples. Samples with initial concentrations below 50 ng µl−1 were purified with AMPure magnetic beads (Beckman Coulter) following the manufacturer's protocol and eluted in smaller volumes to achieve the required final concentration. Samples, for which the standard concentration could not be reached despite AMPure concentration, were also processed but excluded for assay validation steps (see Table S1). To explore effects of low DNA input, we diluted three of the Alpine ibex samples in a dilution series starting from 100 ng µl−1 down to 20 ng µl−1 (5‐fold), 4 ng µl−1 (25‐fold), 0.8 ng µl−1 (125‐fold) and 0.16 ng µl−1 (625‐fold). We selected samples with a concentration of at least 100 ng µl−1, among which we included one corresponding to an individual for which whole‐genome sequencing data are available and two were randomly chosen.

SNP discovery based on whole‐genome sequenced individuals across species

To identify SNPs segregating within the species, we used 29 Alpine ibex whole genome sequences (representing seven different populations). Twenty‐two additional genome sequences (16 domestic goats and six bezoar, Capra aegagrus) were used to identify species diagnostic markers for introgression analyses. Furthermore, 14 additional genome sequences (four Iberian ibex; two Siberian ibex, C. sibirica; two Nubian ibex, C. nubiana; one Markhor, C. falconeri and five sheep, Ovis sp., see Grossen et al., 2020 for details) were included in the analyses to mask polymorphic positions in the domestic goat reference genome (see below). Data from domestic goat, bezoar and sheep were produced by the NextGen Consortium focusing on cattle, sheep and goat genome sequencing (https://www.epfl.ch/labs/nextgen). Raw reads were trimmed using trimmomatic version 0.36 (Bolger et al., 2014). Trimmed reads were mapped with default settings using bwa‐mem (Li et al., 2009) to the domestic goat reference genome (version ARS1, Bickhart et al., 2017) and duplicates were marked using markduplicates from picard (http://broadinstitute.github.io/picard, version 1.130). After genotype calling using haplotypecaller and genotypegvcf (McKenna et al., 2010; gatk, version 4.0.8.0, Van der Auwera et al., 2013), SNPs were removed using VariantFiltration of gatk if: QD < 2.0, FS > 40.0, SOR > 5.0, MQ < 20.0, −3.0 > MQRankSum > 3.0, −3.0 > ReadPosRankSum > 3.0 and AN < 46 (80% of all Alpine ibex individuals). Genome‐wide polymorphism was used in three ways to design amplicons for the high‐throughput sequencing assay: (a) to discover SNPs variable in Alpine ibex; (b) to discover SNPs with fixed alleles between Alpine ibex and domestic goats (i.e., diagnostic markers); and (c) to mask polymorphic sites near targeted SNPs to prevent designing primers spanning polymorphic sites potentially causing amplification drop‐outs.

SNP effect prediction and selection

We annotated SNPs using snpeff version 4.3t 2017‐11‐24 (Cingolani et al., 2012) with gene annotations produced for the domestic goat genome ARS1 as reference (Bickhart et al., 2017). SNPs located in repeat‐masked regions of the reference genome were excluded to avoid designing amplicons in repetitive regions. We designed three main categories of marker sets:

Marker type 1—Genome‐wide, putatively neutral markers

For the design of genome‐wide markers, we focused on SNPs segregating only among Alpine ibex individuals from the Gran Paradiso population, the source of all current Alpine ibex populations (four million SNPs). Our aim was to most accurately reflect segregating neutral polymorphism in the species. Hence, including nonsource individuals could lead to ascertainment bias. Loci were selected with a genotyping rate of ≥70%, a minimal genotype quality of 20 and at a minimal distance of 2.3 Mb between adjacent SNPs (in vcftools, Danecek et al., 2011).

Marker type 2—Diagnostic markers to detect hybridization

To detect recent hybridization events, we identified loci with fixed alleles distinguishing domestic goats and Alpine ibex. Specifically, we chose loci with fixed allele frequency differences between all Alpine ibex and all goats (both domestic goat and bezoar), with a minimal genotyping rate of 80%, a minimal genotype quality of 20 and a minimal distance between adjacent SNPs of 13 Mb (in vcftools, Danecek et al., 2011).

Marker type 3—Immunogenetic markers

Markers in genes relevant to the immune system were selected for two subgroups (here called MHC region and Immune, outside MHC) with a minimal allele count of 1 among Alpine ibex.

3a Immunogenetic markers in the MHC region

Because of the importance of the MHC region for immune functions and evidence for introgression from the domestic goat, we covered the entire MHC region on chromosome 23 (positions 20,892,916–23,588,623 bp). The chromosomal location of the MHC was identified using blastn version 2.7.1+ (Altschul et al., 1990) based on gene sequences reported as belonging to the MHC in a previous goat genome assembly version (CHIR1, Dong et al., 2013; Table S17). Additionally, we used gene ontology (GO) and gene homology to search for all MHC‐related genes on chromosome 23. Matching gene sequences within 250 kb of the homologous MHC region of the CHIR1 assembly were considered as part of the MHC region for further analyses. We split the MHC region into 270 windows of 10 kb using bedtools version 2.27.1 (Quinlan & Hall, 2010) to identify SNPs for the amplicon assay. Evidence for introgression from domestic goat is particularly strong at the DRB exon 2 (Grossen et al., 2014). We hence designed the amplicon assay with more dense SNPs in the DRB region (positions 23,411,211–23,511,211). The DRB region was localized between positions 23,451,944–23,470,477 using blastn version 2.7.1+ (Altschul et al., 1990) and the DRB sequence provided by Dong et al. (2013). A buffer zone of 40,733 bp before and after the gene was added thereby also including the gene ENSCHIG00000008942, which encodes immune‐related functions. We divided the DRB region into 20 windows of 5 kb. Using R, we randomly selected a single SNP from every window in the defined MHC and DRB by prioritizing moderate impact mutations based on snpeff and SNPs in coding regions. We manually selected a marker in the DRB exon 2 at position 23,460,796 bp.

3b Other immunogenetic markers outside MHC region

We targeted immune‐related genes, outside of the MHC region on chromosome 23, based on a candidate gene approach using GO, gene homology and literature reports (Acevedo‐Whitehouse & Cunningham, 2006; Turner et al., 2011). We searched for all MHC‐related GO terms on Ensembl Biomart in the goat genome. Then, we performed interproscan version 5.31–70.0 analyses and searched for protein domain annotations matching the term “histocompatibility.” For the lists by Acevedo‐Whitehouse and Cunningham (2006) and Turner et al. (2011), we focused on 211 genes encoding important immune functions for mammals including interferons, interleukins, TLRs and MHC‐related genes outside of the MHC region on chromosome 23. We prioritized SNPs in coding sequences with a moderate impact prediction based on snpeff and retained only a single SNP per gene (Table S2).

Amplicon design

We extracted a 1,001‐bp sequence from the domestic goat reference genome centred around each target SNP using bcftools version 1.9 ( https://www.htslib.org/ ). The extracted sequences had masked positions for all repetitive regions (masked reference genome) and positions which were polymorphic in Alpine ibex and/or with a minor allele count of three among all other Capra species (samtools version 1.9 with option ‐m; Li et al., 2009). Sequences with more than 50% masked bases or sequences overlapping between different amplicon sets were excluded. We obtained a set of 1,589 sequences for assay primer design by Fluidigm. About 25% of the sequences were rejected by Fluidigm due to the lack of primer options matching the design criteria, leaving 1,265 sequences for oligonucleotide primer synthesis. The targeted amplicon length was 200 bp.

Targeted DNA sequencing library preparation and SNP calling

Libraries were prepared following the manufacturer's protocol for the Juno LP 192.24 integrated fluidic circuits plate (IFC). After loading all reagents on the IFC, target amplicons were generated for each sample through PCR amplification on a specialized thermocycler (Juno system; Fluidigm). A total of 1,265 primer pairs subdivided into 24 multiplex pools were combined in the IFC with each reagent mix consisting of a DNA sample and an individual barcode. After amplification, samples were pooled into a single tube and purified thrice using AMPure XP magnetic beads (Beckman Coulter) to ensure removal of excess primers before adapter ligation. Finally, sequencing adapters were added by PCR to the purified library followed by a final round of purification according to the manufacturer's protocol. The quantity and quality of the library were assessed using a Qubit fluorometer assay (ThermoFisher) and a 4200 TapeStation electrophoresis instrument (Agilent). The final library was sequenced on a single lane of a NextSeq 500 system (Illumina) in midoutput mode adding ~30% PhiX to avoid potential problems due to low sequence complexity. We demultiplexed raw read data using bcl2fastq version 2.19.0.316 and used trimmomatic version 0.38 (Bolger et al., 2014) for quality trimming. Forward and reverse reads were merged using flash version 1.2.11 (Magoč & Salzberg, 2011) and aligned to the goat reference genome ARS1 using bowtie2 version 2.3.5 (Langmead & Salzberg, 2012). Read depths for each step were estimated with multiqc version 1.7 (Ewels et al., 2016). We called SNPs using haplotypecaller, combinegvcfs and genotypegvcfs from gatk version 4.0.1 (McKenna et al., 2010; Van der Auwera et al., 2013). Variant sites were excluded if matching any of the following conditions: QD < 5, MQ < 20, −2 > ReadPosRankSum > 2, −2 > MQRankSum > 2, −2 > BaseQRankSum > 2.

Marker system performance in Alpine ibex populations

To evaluate key performance metrics of the assay, we compared the outcome against two major classes of current population genomics sequencing approaches: RAD‐seq and low‐coverage whole genome sequencing. We used data sets reporting analyses on Alpine ibex populations (Grossen et al., 2018, 2020; Leigh et al., 2018) including the founder population of Gran Paradiso and the three Swiss populations Pleureur, Brienzer Rothorn and Albris. Because high‐coverage whole genome sequencing is not feasible for large‐scale applications, we generated realistic, low‐coverage data sets at approximately 1× coverage. For this, we downsampled 15 whole genome Illumina sequencing data sets (Grossen et al., 2020). We used sambamba version 0.6.6 (Tarasov et al., 2015) to downsample individual bamfiles to a fraction of 0.05 (producing a final coverage of ~1×). We used the software angsd (Korneliussen et al., 2014) to calculate genotype likelihoods with the following options: ‐doGlf 2, ‐doMajorMinor 1, ‐doMaf 1, ‐minMaf 0.05, ‐SNP_pval 1e‐6, ‐minMapQ 20, ‐minQ 20, ‐skipTriallelic 1, ‐uniqueOnly 1, ‐remove_bads 1, ‐only_proper_pairs 1. The resulting likelihoods were used to run pcangsd (Meisner & Albrechtsen, 2018) and ngsadmix (Skotte et al., 2013). For the RAD‐seq data set, we first trimmed reads using trimmomatic version 3.6 (Bolger et al., 2014) and performed read mapping using hisat2 version 2.1 (Kim et al., 2019) on the ARS1 reference genome. We deduplicated bam files using markduplicates from picard version 2.5 (http://broadinstitute.github.io/picard/) and called SNPs on all autosomes (1–29) using the gatk version 4.1 pipeline with haplotypecaller, genomicsdbimport and genotypegvcfs (McKenna et al., 2010). SNPs were flagged using the gatk VariantFiltration tool and removed if any of the conditions were matched: QD < 2.0, FS > 60.0, SOR > 3.0, MQ < 30.0, −12.5 > MQRankSum > 12.5 and −8.0 > ReadPosRankSum. Next, we filtered for SNPs falling within 100 bp of an SbfI restriction cut site identified by in silico analyses with the ENSEMBL tool restrict (Yates et al., 2019). To compare genotyping performance between RAD‐seq and the amplicon sequencing data sets, we filtered for a minimal individual genotyping rate of 0.5 and kept polymorphic sites only. For all further analysis (population differentiation), we kept only bi‐allelic SNPs with a minimal genotyping rate of 0.9 and a minor allele frequency of 0.01. Furthermore, we removed SNPs with a heterozygosity above 0.8 from the RAD‐seq data set, and for the amplicon sequencing, we only retained SNPs from the genome‐wide (neutral) set and removed sites on the Xa/b sex chromosome.

Population genetic data analyses

Genetic data analyses were done using R 4.0.2 (R Core Team 2018). The R package {BioCircos} was used to generate the circular plot. Principal component analyses (PCAs) were performed using the glPCA function from the R package {adegenet}. The R package {hierfstat} was used for F ST calculations. The SNP intersection matrix was visualized with {UpSetR}. Genotype assignment plots were generated using sparse non‐negative matrix factorization algorithms as implemented in the R package {LEA}. For each marker set, we ran 100 repetitions per K (K = 1–10) with entropy = TRUE to find the most likely number of clusters (i.e., K with the lowest entropy). Tajima's D estimates in coding sequences (i.e., all immune‐related genes targeted by the amplicon) were calculated using the 29 Alpine ibex whole‐genome sequencing data sets and the R package {PopGenome}.

RESULTS

Assessment of locus quality across the targeted sequencing assay

We analysed whole genome sequencing data sets comprising a total of 138 million SNPs segregating among all samples (65 whole genome sequences, Capra and Ovis), 5.3 million SNPs segregating among the 29 sequenced Alpine ibex, and 4 million SNPs segregating within the Alpine ibex founder population Gran Paradiso (Italy). Based on this SNP set, we designed 1,256 amplicons and performed targeted amplicon sequencing covering genome‐wide polymorphisms, markers to detect recent hybrids and variants related to immune functions (Figure 1b,c). In contrast to the GT‐seq approach, amplifications were performed in 24 pools of ~52 primer pairs using a microfluidics platform (Figure 1a). We analysed a total of 187 samples (representing 172 individuals) in a single run including 158 Alpine ibex, five domestic goats, three Iberian ibex and six suspected hybrids between Alpine ibex and domestic goats (Table S1). Illumina NextSeq 550 sequencing generated a total of 108 million read pairs after removal of PhiX spike‐ins (32.4 Gb total sequence length). Of the 1,265 target SNPs, 1,189 SNPs passed hard quality filtering in gatk. To assess genotype quality across loci, we first focused only on samples that satisfied our threshold DNA concentration of ≥50 ng µl−1 (see Methods). Furthermore, we required that each analysed individual produced at least 100,000 mapped reads across all loci (Figure 2a). The 65 retained individuals represented seven different Alpine ibex populations, two domestic goat breeds as well as three Iberian ibex. The number of mapped reads ranged from 230,801 to 1,576,615 reads (Figure 2a; Table S3). Based on these high‐quality samples, we found that the median read depth per locus was high with nearly all loci having >20 reads (median across loci =389 reads, Figure 2b). We found seven loci with a median read depth of zero. The highest median read depth was 1,931 for a marker designed in the MHC. A total of 172 loci were invariant among the 65 high‐quality samples and were hence removed (retaining 1,017 SNPs). Next, we analysed the genotyping rate across loci and found that 940 loci were genotyped in all high‐quality samples and 989 loci were genotyped in more than 75% of samples (Figure 2c). We discarded a further 28 of the 1,017 remaining SNPs due to a genotyping rate (among individuals) below 75% (including the seven loci with median read depth of zero) to prioritize loci providing the highest information content across individuals. We retained a total of 989 high‐quality loci for further analysis (Figure 2d).

FIGURE 1

Design and distribution of targeted amplicon sequencing loci. (a) Design principles of the genotyping‐in‐thousands (GT‐seq) and the microfluidics‐based Juno system approaches. GT‐seq relies on the optimization of primer mixes for the simultaneous amplification of hundreds of loci (Campbell et al., 2015; reviewed in Meek & Larson, 2019). By design, microfluidics separates primer pairs into pools to reduce the number of amplifications per microfluidic well (Chiu et al., 2017). (b) Distribution along chromosomes for each marker category (genome‐wide, immune loci, MHC region, diagnostic). (c) Regularly spaced loci across the MHC region on chromosome 23 (20.89–23.59 Mb). Blue segments identify coding sequences in the MHC. The DRB gene region of the MHC (zoom view) was targeted by a denser array of loci

FIGURE 2

Assessment of amplification consistency and genotyping rates. (a) Number of mapped reads per individual sample, colored by population (N = 65). Seven Alpine ibex populations are shown in addition to domestic goats (N = 2) and Iberian ibex (N = 3). Low‐quality input DNA samples (faecal samples etc.) were excluded. (b) Median read depth per locus across 1,017 amplicon loci. Markers are coloured according to their category. The blue line represents the median (389) and the blue area represents the 95% confidence interval for each marker. (c) Variation in genotyping rates across the 1,017 loci. (d) Outcomes of the different filtering stages: <75% genotyping rate: loci removed due to low genotyping rate among the 65 high‐quality individuals (28 loci); Filtered out (gatk): removed based on SNP hard filtering in gatk (76 loci); Non‐variant: monomorphic among the 65 high‐quality samples (172 loci); Passed: a total of 989 loci were retained for further analyses. The marker category of the retained markers is shown below. Marker categories are detailed in Figure 1

Assessment of genotype quality, rate and accuracy

The mean genotype quality (GQ, i.e., the Phred‐scaled confidence that the genotype assignment is correct) across the 65 high‐quality samples was on average 96 (Figure 3a). Taking advantage of four Alpine ibex individuals, which were both whole genome sequenced and genotyped using our assay, we analysed the overall accuracy of genotypes. Among the four individuals, a total of 62 loci were not assigned a genotype based on whole‐genome sequencing. On average, 901 loci (range: 883–910, Figure 3b) showed perfectly matching genotypes and 20 loci (range: 17–27, Figure 3b) showed one mismatching allele (i.e., heterozygote vs. homozygote mismatch). The average genotype quality of the mismatched genotypes was 63 while it was 96 for matching genotypes. We found no complete allelic mismatch in any of the four individuals (i.e., homozygous calls for distinct alleles). We found that genotyping rates >90% across loci are retained by diluting samples 25‐fold from 100 ng µl−1 down to 4 ng µl−1 (Figure 3c). At 0.8 ng µl−1, the genotyping rate was >75% for two out of three samples. The genotyping rate was >23% for 0.16 ng µl−1 of DNA (625‐fold dilution). Furthermore, we analysed four faecal samples (without replication). Three of the four faecal samples had a 25% genotyping rate (i.e., 248–270 genotyped loci). The faecal sample with the lowest quality had a 10.3% genotyping rate corresponding to 102 genotyped loci. Samples with low genotyping rates typically show inconsistently genotyped loci across samples. We found that this was indeed the case with 41–133 of the loci being genotyped in only one out of four faecal samples. A total of 41 loci were genotyped in at least three fecal samples (Figure 3d).

FIGURE 3

Genotyping accuracy and performance. (a) Mean genotype quality (GQ, i.e., the Phred‐scaled confidence that the genotype assignment is correct) for each high‐quality locus across the 65 high‐quality DNA samples. (b) Genotype accuracy assessed by matching recovered alleles from the targeted amplicon sequencing and whole‐genome sequencing of the same four individuals. (c) Assessment of genotyping rates for serial dilutions of three individuals and four faecal samples. (d) Analysis of genotyped loci shared among the four faecal samples. The overlaps show recovered genotypes out of a total of 989 loci

High‐resolution population structure

A major reason for genotyping species of conservation concern is the identification of population subdivisions and admixture events. We expanded our genotyping assay to the full set of Alpine ibex samples (n = 158) spanning the extant distribution range across the Alps (Figure 4a). For this, we focused on the specifically designed genome‐wide markers excluding immune‐related and diagnostic markers as well as sites on the Xa/b sex chromosome (total 617 SNPs). Based on a PCA, we identified three major genotype clusters. The largest cluster was composed of the Gran Paradiso source population and populations reintroduced directly from Gran Paradiso to Italy or Switzerland (Figure 4b). A second cluster grouped the two French populations Haute Maurienne and Champsaur (the latter founded with individuals coming from the former). The third cluster was composed of the isolated Alpi Marittime population, which is thought to have only six effective founder individuals (Terrier & Rossi, 1994). The identified population structure was also supported by individual ancestry coefficients using a sparse non‐negative matrix factorization algorithm (K = 3, Figure S1; Frichot et al. 2014). At the K with the lowest entropy (K = 5, Figure S2), the analyses revealed a fine‐scale population structure: all populations except for the two French populations Haute Maurienne and Champsaur were clearly distinct (Figure 4d). The genotyping assay performed also well for population‐level assignments of the low‐input/quality samples. All four Gran Paradiso faecal samples had similar principal component values (Figure 4b) and population assignments (Figure 4d) as other, high‐quality samples from the same population. We also analysed population differentiation of Alpine ibex using pairwise F ST (Figure 4c). We found a consistent pattern separating Alpi Marittime from all other populations. Furthermore, the two French populations showed relatively low differentiation consistent with their foundation history: the Champsaur population was founded only 25 years ago (fewer than four ibex generations) with 31 individuals coming from Haute Maurienne. The three Swiss populations (Albris, Brienzer Rothorn and Pleureur) were only weakly differentiated from the Gran Paradiso source population as expected from former analysis based on microsatellites (Biebach & Keller, 2009).

FIGURE 4

Population genetic analyses of Alpine ibex. (a) Map showing the extant Alpine ibex distribution (grey) and sampled populations (red). (b) Principal component analysis of 158 Alpine ibex individuals based on 617 genome‐wide SNPs designed for capturing genome‐wide population differentiation. (c) Pairwise F ST matrix of all Alpine ibex populations based on genome‐wide SNPs. Gran Paradiso was the only population surviving the near extinction and hence is the source population of all existing Alpine ibex populations. The Swiss populations Albris, Pleureur and Brienzer Rothorn and the Italian population Alpi Marittime were founded in early 1900 from Gran Paradiso individuals. The populations Tarvisio, Champsaur and Haute Maurienne were later founded independently from Gran Paradiso individuals. (c) Structure‐like analysis (based on sparse non‐negative matrix factorization algorithms) of all Alpine ibex with K = 5. *Faecal samples

Performance contrasts among next‐generation sequencing methods

To objectively assess the performance of the newly developed assay, we analysed genotyping outcomes of four core Alpine ibex populations (Gran Paradiso, Pleureur, Albris and Brienzer Rothorn; Figure 5; Figure S2). After quality filtering of each data set to ensure objective comparisons (see Methods), we retained 892 loci from targeted amplicon sequencing (n = 75 individuals), 26,547 RAD‐seq loci (n = 82 individuals) and 3 million low‐coverage whole genome sequencing loci (n = 15 individuals). Overall, 97% of the individuals had a per‐individual genotyping rate of ≥90% for targeted amplicon sequencing loci contrasting with 23% of the individuals genotyped at ≥90% for RAD‐seq loci (Figure 5a). Enforcing a per‐locus genotyping rate of ≥90% over all individuals, 96% of targeted amplicon sequencing loci but only 39% of RAD‐seq loci were retained (Figure 5b). Locus and individual‐level genotyping rates cannot meaningfully be retrieved from genotype likelihood‐based analyses (low‐coverage whole genome sequencing data set). We performed comparative population differentiation analyses and found that the global F ST ranged between 0.071 (RAD‐seq, 8,316 SNPs) and 0.077 (targeted amplicon sequencing, 588 SNPs). Pairwise F ST estimates were also similar among marker systems (Figure 5c). PCAs constructed from targeted amplicon sequencing and RAD‐seq markers clearly resolved the four populations (Figure 5d,e). However, the low‐coverage whole genome sequencing did not resolve Albris and Brienzer Rothorn populations (Figure 5f). The first and second principal component axes explained 6.2% and 4.9% for the targeted amplicon sequencing (Figure 5d), 5.2% and 4.2% for the RAD‐seq (Figure 5e) and 23.8%–6.1% for the low‐coverage whole genome sequencing (Figure 5f), respectively. Genotype assignments to clusters showed clear population differentiation for the targeted amplicon sequencing markers (Figure 5g) and slightly weaker resolution for RAD‐seq markers (Figure 5h). Low‐coverage whole genome sequencing genotyping clearly separated Gran Paradiso and Pleureur populations but again failed to resolve Albris and Brienzer Rothorn populations (Figure 5i).

FIGURE 5

Performance comparisons between targeted amplicon sequencing, RAD‐seq and low‐coverage whole genome sequencing on four Alpine ibex populations. (a) Individual genotyping rate ranked by individual for amplicon sequencing and RAD‐seq. (b) Proportion of SNPs retained at a genotyping rate of 90% as a function of the percentage of individuals included. The comparison is only meaningful for amplicon sequencing and RAD‐seq. (c) Boxplot of all pairwise FST estimates among populations for the amplicon sequencing and RAD‐seq. (d–f) Principal component analysis of population differentiation for (d) the amplicon sequencing (plotted are −PC1 against −PC2), (e) RAD‐seq and (f) low‐coverage whole‐genome sequencing. Structure‐like analysis (based on sparse non‐negative matrix factorization algorithms) of (g) the amplicon sequencing, (h) RAD‐seq and (i) low‐coverage whole‐genome sequencing data sets; (i) is based on an ngsadmix analysis (Skotte et al., 2013)

Detection of recent hybrids and introgression tracts

We analysed the performance of the SNP markers specifically designed to detect introgression from domestic goats into Alpine ibex. To assess the power to discriminate genotypes suspected to be from hybrid individuals, we performed a PCA including all 158 Alpine ibex individuals, suspected hybrids based on field reports of unusual phenotypes (n = 6), domestic goats (n = 5) and Iberian ibex (n = 3; Figure 6a). As expected, the first PC clearly separated domestic goats, Iberian ibex as well as Alpine ibex. The second PC differentiated the Alpi Marittime population from all other Alpine ibex. The suspected hybrid individual from Alpi Marittime (AM‐H) clearly clustered with domestic goats. Two additional suspected hybrids (GR‐ib1, GR‐ib2), as well as a potential albino individual (FR‐blanc has largely whitish fur, but no red eyes) clustered with Alpine ibex. Three suspected hybrids (TI‐ib, GPHB1 and GP‐ib‐V02‐17) were located near the midpoint between domestic goats and Alpine ibex matching expectations for recent (F1 or backcross) hybrid genotypes.

FIGURE 6

Performance of targeted amplicon sequencing to detect recent hybridization events. (a) Principal component analysis of 168 individuals (excluding faecal samples) including domestic goat and Iberian ibex (N = 5 and 3) and suspected hybrids based on phenotypic observations (N = 6). The PCA was performed based on 617 genome‐wide SNPs designed for capturing putatively neutral population differentiation. (b) Analysis of diagnostic markers for the detection of hybrids for six suspected hybrids and one potential albino. The proportion of diagnostic markers assigned to domestic goat are shown, where the proportion is calculated as (1 × homozygous goat genotype + ½ × heterozygous genotype)/(total genotypes at diagnostic markers). (c) Chromosome painting using diagnostic markers along the 29 autosomes. Colours identify Alpine ibex (black), heterozygous (purple) and domestic goat (red) genotypes. Sample GPHB1 was of low DNA quality and quantity (drop of blood conserved on FTA filter paper; see Methods)

Using the goat–Alpine ibex diagnostic marker set, we analysed the recency of the hybridization event by identifying contributions from each parental species (Figure 6b). Of the six suspected hybrids and one potential albino, two individuals were confirmed to be Alpine ibex (FR‐blanc, GR‐ib1) and one a domestic goat (AM‐H). Our results confirm that the unusual phenotype of FR‐blanc reported from the field was not caused by recent domestic goat introgression. Individual GR‐ib1 was suspected to be a hybrid because it was behaving in very unusual ways, seeking proximity to buildings. Individual AM‐H was reported to resemble domestic goat but was living among Alpine ibex. Another suspected hybrid (GR‐ib2) with white hoofs showed a weak sign of domestic goat introgression (only one marker on chromosome 10 was heterozygous) suggesting a potential backcross. Note that GR_ib2 showed no clear differentiation from Alpine ibex based on the genome‐wide marker PCA underlining the usefulness of specifically designed diagnostic loci. Two individuals (GP_ib_V02_17 and GPHB) showed signs of ~25% domestic goat introgression (i.e., probable F 2 backcrosses). GP_ib_V02_17 had darker fur than commonly seen in Alpine ibex and an unusual horn shape (no nodes and a triangular transverse section). Chromosome painting showed one individual was being heterozygous for all diagnostic markers (TI_ib, Figure 6c), hence representing very probably an F 1 hybrid. This individual was observed going into a stable following goats and had an unusual horn shape. The DNA from GPHB was from a small amount of blood stored on FTA paper explaining the poor genotype quality. Field reports suspected an F 1 hybrid but both the contribution plot (Figure 6b) and the chromosomal painting (Figure 6c) suggest an F 2 hybrid. Our analyses show the power of highly discriminatory markers to detect recent hybridization and introgression events in a pool of individuals with field‐reported, suspected admixture.

Immunogenetics of Alpine ibex populations

The targeted amplicon sequencing specifically focused on 297 polymorphisms in immune‐related genes within and outside the MHC. Using whole‐genome sequencing data for Alpine ibex (n = 29), we found that Tajima's D in Alpine ibex ranged from −2.4 to 4.5 across all genes encoded in the MHC region and from −1.4 to 2.5 for all immune‐related genes outside of the MHC region targeted by the amplicon sequencing. The median Tajima's D was lower in the MHC (0.45) compared to immune loci outside of the MHC (1.3, Figure 7a,b). The high Tajima's D values in some immune‐related loci suggest long‐term maintenance of alleles through balancing selection (Figure 7a). Comparing different sets of amplicon targets in our assay, we found that the genome‐wide markers aimed at resolving population structure showed overall the highest average heterozygosity (Figure 7c). The lowest genome‐wide heterozygosity was found in Alpi Marittime consistent with the severe founding bottleneck. Immune‐related and MHC loci showed consistently lower levels of heterozygosity compared to genome‐wide loci. A notable exception was the Alpi Marittime population where the MHC showed higher levels of heterozygosity (0.25) than other immune‐related and genome‐wide markers (Figure 7c). The MHC of the Albris population showed a surprisingly low average heterozygosity (0.11) compared to genome‐wide markers. On a PCA, the MHC in domestic goat and Iberian ibex showed a low degree of differentiation (Figure 7d). This is probably explained by the focus on segregating polymorphism in Alpine ibex only. MHC genotypes showed tight clusters among Alpine ibex individuals but only weak population signatures, which is in marked contrast to genome‐wide markers (Figures 4b and 7d). Clusters of nearly identical Alpine ibex genotypes were generally composed of genotypes from multiple populations. The two French populations Haute Maurienne and Champsaur shared most MHC genotypes. Individuals from the Alpi Marittime were largely distinct from all other populations with the exception of a shared genotype with the Italian population Tarvisio. Some individuals from Alpi Marittime were translocated to Tarvisio in 1993, which may explain our finding.

FIGURE 7

Genome‐wide immunogenetics survey of Alpine ibex populations. (a) Boxplot of Tajima's D estimates based on whole‐genome sequencing of 29 Alpine ibex of immune‐related loci covered by the targeted amplicon sequencing assay. The MHC region is shown separately from other immune‐related loci. (b) Tajima's D estimates shown per each gene represented by immune‐related loci covered by the targeted amplicon sequencing assay. (c) Targeted amplicon sequencing analyses of heterozygosity across Alpine ibex populations (N = 154, excluding faecal samples). Heterozygosity is shown separately for genome‐wide SNPs, the MHC region and other immune‐related SNP loci. (d) Principal component analysis of Alpine ibex and sister species based on 138 SNPs genotyped in the MHC region using targeted amplicon sequencing

DISCUSSION

Genetic monitoring is central to many population surveys and conservation efforts. Here, we developed an accurate and versatile microfluidics‐based tool for the monitoring of Alpine ibex population health. The set of nearly 1,000 loci enables the concurrent assessment of population structure, the detection of recent hybridization events and immune function‐related genotypes. We show that the assay performs well in comparison to restriction‐site based DNA sequencing and low‐coverage whole genome sequencing with a probable wider range of acceptable input material. Assessments of immunity‐related genotypes can be integrated in individual translocation recommendations for conservation management of the species.

Effective high‐throughput genotyping relies on the consistent amplification of a large number of loci, robustness to variation in input DNA quality and affordable costs. We show that the microfluidics‐based targeted sequencing approach produces largely uniform coverage both across 989 loci and individuals. Integrating whole‐genome sequencing derived SNPs from different species into the design of amplicons allowed us to amplify across Alpine ibex, Iberian ibex and the domestic goat. This extends the usability of the assay with the caveat that the potential ascertainment bias in other species should be considered because SNPs were selected based on polymorphism in Alpine ibex only. Furthermore, the identified genotypes were well validated by cross‐referencing with the whole‐genome sequencing data sets. For Alpine ibex, uniform, high‐accuracy genotyping with low input DNA has previously only been achieved through microsatellite marker analyses (Biebach & Keller, 2009) or SNP chips designed on related model species (Grossen et al., 2014). An SNP chip genotyping study based on a 52 K Illumina Goat SNP Chip recovered 677 polymorphic markers out of a total of ~52,000 markers known to be polymorphic among domestic goats (Grossen et al., 2014; Tosser‐Klopp et al., 2015). This shows how even SNP chips designed for closely related species can be largely unsuitable and may suffer from substantial ascertainment bias. High‐throughput methods such as RAD‐seq and GBS sequencing are widely used next‐generation sequencing approaches for nonmodel species often producing at least tens of thousands of SNPs. However, such restriction‐based reduced representation approaches do not allow the targeting of specific loci (e.g., immune‐relevant loci) and often produce highly uneven read depth across loci (Andrews et al., 2016; Jiang et al., 2016). Furthermore, we show that replicating RAD‐seq and low‐coverage whole genome sequencing on the same set of populations analysed in our targeted amplicon sequencing approach produces a highly similar resolution of the genetic structure. The comparatively low number of SNPs assayed in the targeted amplicon sequencing is probably compensated for by the highly consistent genotyping rates across loci and individuals.

We show that the microfluidics approach for amplicon sequencing produces consistent amplification across nearly all loci if minimum DNA requirements are fulfilled. We have set a minimum threshold of 50 ng µl−1 (or a total of 100 ng), which is usually achievable for well‐preserved tissue samples. We have also investigated the potential to recover genotypes from degraded and low‐concentration samples such as faecal material. Faecal sampling is a widely used noninvasive sampling method and sometimes the only option for elusive species (Beja‐Pereira et al., 2009). Faecal samples produce low‐quality (i.e., degraded) and low‐quantity DNA. In addition, faecal DNA samples can be heavily contaminated with bacteria, plant or prey DNA. Hence, the actual endogenous DNA is often much below the measured total DNA concentration. We analysed faecal DNA quantities from 4–23 ng µl−1 but found that the proportion of successfully amplified loci was comparable to diluted nonfaecal DNA at around 0.1 ng µl−1. Our preliminary analyses using four faecal samples suggest that the new amplicon assay produces sufficient genotyping data to assign the origin of faecal samples to individual populations. To assess the resolution of the genome‐wide set of amplicon markers, we have genotyped a collection of Alpine ibex samples representing all major reintroduction events. Alpine ibex were limited to the Gran Paradiso National Park in Italy in the 19th century. Starting in the early 20th century, populations were introduced independently from Gran Paradiso to Switzerland, France and Italy and from there further populations were founded. The population structure of extant Alpine ibex populations is dominated by signals of population reintroductions and translocations (Biebach & Keller, 2009; Grossen et al., 2018). The Alpine ibex genotypes assessed by our targeted sequencing approach confirmed all major aspects of the reintroduction history including placing Gran Paradiso at the centre of the extant genetic diversity. Populations reintroduced directly from Gran Paradiso were at the closest genetic distances to the source population yet showed distinct clustering as expected from the strong bottlenecks imposed by the translocation of few individuals. A major concern for wild species co‐occurring with closely related domestic animals is the potential for hybridization and introgression. Alpine ibex populations across Europe are monitored for the presence of individuals with atypical phenotypes (e.g., Giacometti et al., 2004; Steyer et al., 2016; Todesco et al., 2016). We assessed genotypes of seven suspected hybrids (including one potential albino) and were able to show that only three individuals are clearly identifiable as recent hybrids. One individual showed possible signs of past introgression and three suspected hybrids clustered either with domestic goat or Alpine ibex genotypes, highlighting the importance of genetic monitoring of suspected hybridization events.

A major genetic factor for the long‐term survival of endangered species is genetic diversity at immune loci, in particular the MHC. We successfully amplified hundreds of loci involved in key immune functions focusing on nonsynonymous polymorphisms. In parallel, we amplified a dense array of loci spanning the MHC. This 2‐Mb locus on chromosome 23 recently received genetic material from domestic goats, regenerating heterozygosity at the DRB locus (Grossen et al., 2014). Genotypes across the MHC clustered tightly among populations, revealing that Alpine ibex were genetically impoverished across the entire MHC. Comparative analyses with domestic goat breeds will enable a refinement of our understanding of historical and potentially ongoing introgression into the Alpine ibex gene pool. Our findings of low diversity at the MHC are also in accordance with previous studies based on microsatellite and RAD genotyping (e.g., Alasaad et al., 2012; Brambilla et al., 2018; Grossen et al., 2014). Underlining the relevance of performing large‐scale targeted amplicon sequencing of immune loci are recent findings that high heterozygosity at the MHC was correlated with higher resistance to infectious keratoconjunctivitis, a major disease factor in some Alpine ibex populations (Brambilla et al., 2018). Surveillance of disease susceptibility through genetic analyses of immune loci will also significantly improve conservation strategies by informing the choice of founder individuals prior to reintroductions.

Our microfluidics‐based amplicon sequencing approach joins GT‐seq, the currently most widely used method for amplicon sequencing (Campbell et al., 2015; Meek & Larson, 2019), for the survey of wild populations. The major advantage of GT‐seq relative to our approach is the possibility to multiplex thousands of samples (relative to the hundreds currently possible with microfluidics). In particular for conservation genetics applications, the cost per sample is an important factor to consider. We refer here to Campbell et al. (2015) to compare costs with the GT‐seq approach. To provide useful comparisons, we assume identical costs for DNA extraction, library purification, quality checks and sequencing, because we expect these steps to be identical between the procedures. For 2,000 samples, GT‐seq costs were estimated to be 1.23 USD per sample for points 1 and 6–8 (Appendix S1; Campbell et al., 2015). Our protocol was based on 384 single‐index sequencing barcodes, but this could easily be modified to reach higher plexity. The remaining costs for oligo synthesis and amplification were estimated to be 2.75 USD per sample (Campbell et al., 2015; Table S4). In direct comparison, our microfluidics approach costs 14.40 USD per sample, including the cost of outsourced primer design, Juno IFC chips and kits (Table S4). Note that our cost estimates are for a Swiss research facility in 2020 compared to Campbell et al. (2015). The Juno microfluidic system from Fluidigm used in this study is also not standard molecular laboratory equipment. In summary, GT‐seq is more cost‐effective per sample and is independent from commercial instruments. Note though that GT‐seq may imply more development time (~4 months according to Meek & Larson, 2019) to ensure sufficiently balanced amplicons in a single pool (or multiple pools). To our knowledge, the sensitivity to input DNA quantity and quality of GT‐seq and the microfluidics approach has not been systematically compared. We expect both methods to perform similarly well, however. Recent targeted sequencing applications focused on amplifying fragments of the MHC, sex‐specific markers and marker sets to resolve parentage in large sample pools (Bootsma et al., 2020, Fuselli et al., 2018, Razali et al., 2017, Yu et al., 2020, reviewed in Meek & Larson, 2019). Our approach integrating whole‐genome resequencing data sets both to identify the most valuable polymorphisms to target and to mask polymorphisms interfering with amplification should be universally beneficial. Provided that at least a draft genome is available from a closely related species, even low‐coverage data sets will provide important benefits to the assay development and validation.

In conclusion, our microfluidics‐based targeted amplicon assay allows the simultaneous monitoring of multiple genetic risk factors and demonstrates how next‐generation sequencing techniques can be adapted for the needs of population genetic surveys and conservation management. The whole genome sequencing‐informed approach enabled us to select a highly specific set of loci to simultaneously address questions of population structure, recent hybridization events and how polymorphism is shaped across major components of the immune system. Efficient and precise characterization of individual genotypes can be translated into recommendations on how to prioritize translocation events and replenish genetic diversity at immune loci. The versatility to amplify the same loci across related species also enables powerful screens for recent introgression events. Our study shows the relevance of bridging population genomic investigations with assays that can be realistically implemented into population genetic surveys and decision‐making for conservation management.

CONFLICT OF INTEREST

The authors declare that there is no conflict of interest.

AUTHOR CONTRIBUTIONS

C.K., C.G. and D.C. conceived the study, C.K., D.W., G.C., I.B. and C.G. performed analyses, A.B. contributed samples, D.L. contributed data sets, C.K., C.G. and D.C. wrote the manuscript. All authors commented on the text.

Supplementary Material

Click here for additional data file.

Table S1‐S4

Click here for additional data file.