Improved haplotype resolution of highly duplicated MHC genes in a long-read genome assembly using MiSeq amplicons

Genomic reference of the great reed warbler MHC region

The great reed warbler genome was characterized through a de novo genome assembly reconstructed from long-read sequencing (PacBio, Menlo Park, CA, USA), short-read sequencing (Illumina, San Diego, CA, USA), and optical mapping data (BioNano, San Diego, CA, USA) from a female individual, our focal individual (Sigeman et al., 2021). The final genome assembly comprised 3,013 scaffolds, with a total size of 1.2 Gb of which half of the assembled genome was included in scaffolds that were at least 21.4 Mbp long (N50). As a first step, the genome was assembled using the Falcon assembler (v.0.4.2) (Chin et al., 2013), which separates sequences into primary contigs and associated contigs. Associated contigs represent alternative haplotypes and the primary contigs were used for subsequent scaffolding steps. As some primary contigs may represent incorrectly assigned associated contigs, which is more likely to happen in species with high heterozygosity, the Purge Haplotigs pipeline (Roach, Schmidt & Borneman, 2018) was used to remove sequences that were putative alternate haplotypes of other sequences based on read coverage and pairwise sequence similarities. Full-length MHC variants were identified in both versions of long-read assemblies, GRW Falcon-2017 and Purge Haplotigs, by blasting complementary DNA Sanger sequences, exons 2–4 as described in Westerdahl et al. (2022), from eight MHC-I and seven MHC-IIB loci previously characterized from expressed messenger RNA in the focal individual (Westerdahl, Wittzell & von Schantz, 1999; Westerdahl, Wittzell & von Schantz, 2000). Note, MHC ‘variants’ are called alleles when we reference to the Falcon assembly and genes when we reference to the Purge Haplotigs.

Amplicon sequencing

Amplicon sequencing was performed using high-throughput sequencing technology (Illumina MiSeq, San Diego, CA, USA) in order to genotype MHC-I and MHC-IIB genes in a great reed warbler family including the focal individual, two siblings as well as both parents. The birds are part of a long-term monitored wild population at Lake Kvismaren in Sweden (Hansson et al., 2007). The focal individual of this study was sacrificed in 1996 and all individuals in this study were blood sampled with the permission from the Swedish Environmental Protection Agency (Permit number M98-96). Genomic DNA (gDNA) was extracted using a modified phenol-chloroform DNA extraction method of blood samples (Sambrook, Fritsch & Maniatis, 1989). MHC-I loci were genotyped by PCR amplification of either a 246 bp or a 262 bp fragment within exon 3 with previously designed primers HNalla/HN46 (Westerdahl et al., 2004a; Westerdahl et al., 2004b) as described in Roved et al. (2018) and newly designed primers HNalla-1/R3Ex3b, respectively (Table S1). For MHC-IIB genotyping, four newly designed primer pair (PP) combinations (PP1, PP2, PP3, and PP5) were used to amplify MHC-IIB exon 2 and the fragments length obtained varied between 256 bp and 351 bp (Table S1).

Filtering method and allele selection

For both MHC-I and MHC-IIB, only the amplicons sequenced from the focal individual and its parents were used in the present study. All three samples were run in replicates for all primer pair combinations except for HNalla/HN46. Obtained sequences were trimmed of the adapters and primer sequences using Cutadapt (Martin, 2011). Trimmed sequences were imported in RStudio v.1.2.1335 (R Studio Core Team, 2020) and filtered with the R packages DADA2 (Callahan et al., 2016) as described in Stervander et al. (2020) and Drews & Westerdahl (2019) using the function filterAndTrim and adjusting settings for expected error rate and quality cut-off for each primer pair combination.

For the parents and the focal individual, the numbers of filtered reads for both MHC-I primer pairs (PPs) ranged between 19,193 and 34,310 (median 26,907). The maximum percentage of reads for an MHC-I amplicon allele (number of reads for a particular allele divided by the total number of reads for all alleles in an individual) was 10.8–12.9% and the minimum percentage of reads was 0.12–0.23% for all individuals. In addition, to keep an allele we took into account that (i) it should be found in both replicates of the same individual within PP combination or (ii) it should be found with both PPs in the same individual or (iii) it should be found in at least two individuals in the family independently of the PP combination (because not all alleles are expected to be amplified by all PPs). Following these rules, we retained two additional MHC-I alleles with low number of reads (0.06%), but which were amplified in the focal individual and its parents.

For MHC-IIB, we selected data from the two best-performing primer pair combinations (PPs) i.e., for which we obtained the largest number of different MHC-IIB amplicon alleles post-trimming (PP1: 93 and PP5: 89). Only sequences of the expected length or that were shorter by a combination of a multiple of 3 bp (considering the length of a full codon) were kept. The two other MHC-IIB PPs gave a low number of unique amplicon alleles in the focal individual (PP2: 74 and PP3: 65) and data from both PPs were discarded (Table S1). The number of filtered reads for both PPs ranged between 29,960 and 51,190 (median 40,571). The maximum percentage of reads for an MHC-IIB amplicon allele was 2.9–5.8% and the minimum percentage of reads was 0.12–0.23% for all individuals. MHC-IIB genes are highly duplicated in the great reed warbler and many more paralogs are amplified compared to MHC-I (Westerdahl et al., 2022). Consequently, the read depth observed (i.e., percentage of reads) for all MHC-IIB amplicon alleles is lower as many more alleles are amplified simultaneously in a single individual. In the focal individual, the number of unique MHC-I and MHC-IIB amplicon alleles amplified was assessed for each primer pair combination (Table S1). Finally, we identified identical amplicon sequences amplified by both PPs for MHC-I and for MHC-IIB using the online program Seqeqseq (http://130.235.244.92/apps/seqeqseq.html).

Mapping of MHC-I and MHC-IIB HTS-amplicon alleles

In order to identify paternally and maternally inherited annotated MHC alleles, MHC-I and MHC-IIB amplicon alleles found in the focal individual were sorted into three categories. Amplicon alleles found in the focal individual and the father were identified as paternal (P) whereas amplicon alleles found in the focal individual and the mother were identified as maternal (M). Amplicon alleles found in the focal individual and both parents were identified as unresolved alleles (U).

For MHC-I, amplicon alleles amplified with both PPs represent a fragment of 246 or 262 bp within exon 3, which has a total length of 270 bp (Table S1). For MHC-IIB, the amplicon alleles amplified using PP5 represent a fragment of 256 bp within exon 2 whereas amplicon alleles amplified using PP1 were trimmed to correspond to the length of exon 2, which is 270 bp (Table S1). All MHC-I and MHC-IIB amplicon alleles found in the focal individual were mapped using Geneious Prime^® 2020.1.2 (v. 11.0.6) to the original great reed warbler long-read genome assembly (GRW Falcon-2017) and to post-processed haploid genome assembly (Purge Haplotigs) which has gone through the post-assembly filtering procedure (Fig. 1, Tables S2–S4). We first used the Geneious RNA mapper with the following settings: Custom Sensitivity and allowing 0% mismatches per read. We also performed a complementary mapping step with more relaxed settings using the Geneious mapper (custom sensitivity, allowing 1–4% mismatches per read). Allowing mismatches can help to identify putative alternate (second) alleles for a gene copy if the alleles only differ by one or two nucleotides. It is also a way to account for remaining sequencing errors or indels in genome assemblies during the mapping step. Finally, nucleotide sequence similarities (based on pairwise similarity matrices) were compared between amplicon alleles and the corresponding coding exon sequences of annotated MHC-I and MHC-IIB alleles in both long-read genome assemblies.

Amplicon alleles were assigned to annotated MHC alleles when they shared identical nucleotide sequences (Fig. 1, Tables S3–S4). Thus, an amplicon allele could be assigned to multiple annotated MHC alleles if the coding sequence was identical with each of them (based on nucleotide sequence of exon 3 for MHC-I alleles and exon 2 for MHC-IIB alleles). Additional verifications were performed to confirm the assignment of amplicon alleles that were not identical to any annotated MHC alleles based on nucleotide sequence similarity. Two different conditions were tested: (1) to be assigned, a single amplicon allele should have at least 99% nucleotide sequence similarity (corresponding to one or two nucleotide substitutions) with a unique annotated MHC allele; (2) to be considered as two allelic variants of a unique annotated MHC gene, two amplicon alleles should map and share high nucleotide sequence similarity (>99%) with each other and with a unique annotated MHC allele. Amplicon alleles mapping with lower nucleotide sequence similarity were not considered.

Coding sequence similarity of MHC-IIB haplotypes

In the great reed warbler, the MHC-IIB genes have expanded rather recently and most MHC-IIB gene copies are similar to each other (Westerdahl et al., 2022). However, MHC-IIB alleles originating from a single gene copy are expected to be more similar than alleles originating from different gene copies. Thus, we expect MHC-IIB alleles from two complementary haplotypes to be more similar to each other than MHC-IIB alleles from different genomic regions, as measured by pairwise sequence distances. We extracted all 56 full-length MHC-IIB genes in open reading frame (ORF, putatively functional) found in the haploid representation of the great reed warbler genome assembly Purge Haplotigs and aligned the sequences using the Geneious alignment in Geneious Prime^® 2021.1.1. We computed between-group mean distances by calculating nucleotide pairwise distance between all MHC-IIB gene copies in scaffolds which had >5 tandemly duplicated genes (scaffolds: Aaru-DAB*554, Aaru-DAB*357, Aaru-DAB*120, Aaru-DAB*301 and Aaru-DAB*45) in MEGA version 11 (Tamura, Stecher & Kumar, 2021), considering pairwise deletion and performing 1,000 bootstraps replicates. The same analyses were performed for full-length amino acids sequences. Westerdahl et al. (2022) suggested that the two scaffolds Aaru-DAB*554 (maternally inherited) and Aaru-DAB*357 (paternally inherited) may represent the two parental haplotypes of a single MHC-IIB genomic region, i.e., the two plausible complementary haplotypes. To investigate this further we computed the mean nucleotide pairwise distance between each annotated MHC-IIB gene belonging to scaffold Aaru-DAB*554 and each of the MHC-IIB genes at the four scaffolds with >5 tandemly duplicated genes. The same analysis was also performed for Aaru-DAB*357.

MHC diversity in the long-read genome assemblies

In the original long-read genome assembly of the focal individual, GRW Falcon-2017, 25 MHC-I and 96 MHC-IIB full-length alleles were annotated in primary contigs (four MHC-IIB alleles were found in associated contigs), hence the primary contigs of GRW Falcon-2017 represent a diploid version of the MHC region. Following the post-assembly procedure (Purge Haplotigs), which aimed to remove sequences that represent alternative haplotypes of other sequences and give a haploid representation of the genome, 18 MHC-I and 66 MHC-IIB full-length genes were included and of these, 15 MHC-I and 56 MHC-IIB genes contained a predicted full-length open reading frame (ORF, Westerdahl et al., 2022). In total, 16 scaffolds with 1–19 MHC genes were identified in the great reed warbler Purge Haplotigs assembly (three MHC-I scaffolds, 11 MHC-IIB scaffolds and two scaffolds with both MHC-I and MHC-IIB genes). Of the 16 MHC scaffolds, five scaffolds with MHC-IIB genes were larger (103, 178, 175, 510 and 755 kbp) and had 6–13 tandemly duplicated MHC-IIB genes, i.e., an MHC gene with at least one MHC gene as nearest neighbor within a short intergenic distance (on average 6,443 bp between tandemly duplicated MHC-IIB genes).

MHC diversity based on amplicon HTS

Twenty-nine MHC-I alleles and 95 MHC-IIB alleles were amplified in the focal individual using amplicon HTS (22 MHC-I exon 3 alleles and 85 MHC-IIB exon 2 alleles contained an ORF; Table 1). The primers for amplicon sequencing were carefully designed to amplify the majority of all available MHC-I and MHC-IIB alleles in the great reed warbler (Table S1). Twelve out of 29 MHC-I amplicon alleles amplified in the focal individual were identified as paternal alleles (P; blue), seven were identified as maternal alleles (M; yellow) and nine were unresolved (U; turquoise; Fig. 2A). For MHC-IIB, 36 out of 95 amplicon alleles in the focal individual were paternal, 34 alleles were maternal, and 25 alleles were unresolved (Fig. 2A).

MHC diversity in amplicons compared to MHC diversity in the GRW Falcon-2017 assembly

For MHC-I, 84% (21 alleles of the 25 full-length alleles) of the annotated alleles in GRW Falcon-2017 were successfully detected by amplicon alleles in the focal individual (Fig. 3, Fig. S1, Table S3). For MHC-IIB, the success rate was slightly lower, 78% (78 of the 100 full-length alleles (96 on primary and four on associated contigs)) of the annotated MHC alleles in GRW Falcon-2017 were detected using amplicon sequencing (Fig. 3, Fig. S1, Table S4). The total MHC diversity in amplicons, 29 MHC-I and 95 MHC-IIB alleles, was comparable to the diversity of the annotated MHC alleles in GRW Falcon-2017 assembly (25 MHC-I and 100 MHC-IIB, including associated contigs; Table 1, Fig. S1).

Sixteen MHC-I exon 3 amplicon alleles mapped to 21 annotated MHC-I genes in GRW Falcon-2017 (Fig. 3, Tables S2–S3). Four of these 16 MHC-I amplicon alleles mapped multiple times to a total of nine annotated MHC-I alleles (one amplicon allele mapped to three annotated alleles and three amplicon alleles mapped to two annotated alleles each; Fig. S1), and therefore most of the annotated alleles in GRW Falcon-2017 (84%) were detected. However, 13 of 29 MHC-I amplicon alleles (45%) did not map to any full-length annotated MHC-I allele, suggesting that the genome assembly lacks annotations for many MHC-I alleles and/or that several MHC-I exon 3 amplicon alleles are PCR products from gene fragments. The latter is supported by two of the 13 MHC-I amplicons, which map to locations without full-length annotated MHC-I alleles in the GRW Falcon-2017 assembly (Fig. S1).

Fifty-eight MHC-IIB exon 2 amplicon alleles mapped to 78 annotated full-length MHC-IIB alleles in GRW Falcon-2017 (Fig. 3, Table S4). Ten of these amplicon alleles mapped multiple times to a total of 32 annotated MHC-IIB alleles that shared identical nucleotide sequences for exon 2 but differed elsewhere in the gene sequence (Fig. S1). However, 37 of 95 MHC-IIB amplicon alleles (39%) did not map to any full-length annotated MHC-IIB allele, suggesting an incomplete genome assembly and/or that MHC-IIB exon 2 amplicons are PCR-products from gene fragments. Again, the latter scenario is supported by seven of the 37 amplicons, which map to locations without full-length annotated MHC-IIB alleles in the GRW Falcon-2017 assembly (Fig. S1).

MHC diversity in amplicons compared to MHC diversity in the Purge Haplotigs assembly

For MHC-I, 89% (16 of the 18 full-length genes) of the annotated MHC genes in Purge Haplotigs were successfully detected by amplicon alleles in the focal individual (Fig. 2B, Table S3). The success rate was again slightly lower for MHC-IIB, where 82% of the annotated MHC genes in Purge Haplotigs (54 of the 66 full-length genes) were detected using amplicon sequencing (Fig. 2B, Table S4). Thirteen MHC-I amplicon alleles mapped perfectly to 16 annotated full-length MHC-I genes (Acar-UA genes, Fig. 2B), of which two MHC-I amplicon alleles mapped multiple times (Table S3; note that one amplicon allele in the focal individual was not found in either parent, see Tables S2–S3). Fifty-three MHC-IIB exon 2 amplicon alleles mapped to 54 annotated MHC-IIB genes (Acar-DAB genes, Fig. 2B, Table S4) of which three MHC-IIB amplicon alleles mapped multiple times.

Haplotype sorting of tandemly duplicated MHC genes in the Purge Haplotigs assembly

The parental amplicon alleles that mapped to the annotated full-length tandemly duplicated MHC-I and MHC-IIB genes in the Purge Haplotigs assembly show which scaffolds that are dominated by paternal (P) and maternal (M) MHC genes. Thus, since seven MHC-I genes were assigned as paternal and only one was assigned as maternal on scaffold 508 (Acar-UA*508), this scaffold is likely to represent a paternal haplotype (Fig. 2B, Table S3).

There were 54 MHC-IIB genes in the Purge Haplotigs assembly with information about parental origin (13P, 19M, and 22U, Fig. 2B, Table S4), and among the five larger scaffolds (103–755 kbp) with 6–13 tandemly duplicated MHC-IIB genes, one scaffold was putatively paternal (Acar-DAB*357), two were putatively maternal (Acar-DAB*554 and *120) and two were mixed (Acar-DAB*45 and *301).

One way of assessing whether the number of MHC-IIB genes is overestimated in the Purge Haplotigs assembly is to investigate whether any of the paternal and maternal scaffolds with tandemly duplicated genes are likely to represent haplotypes of the same genomic region. The nucleotide sequence of the MHC-IIB alleles is expected to be more similar within than between paralogs, so the average nucleotide distance between alleles from complementary haplotypes should be lower than between non-complementary haplotypes. We thus compared the average nucleotide pairwise distances (p-distances) of the tandemly duplicated MHC-IIB genes containing an ORF between the five larger scaffolds to investigate whether any pair-wise comparison among the scaffolds represent putatively complementary haplotypes. The mean nucleotide p-distance between genes in the two scaffolds Aaru-DAB*554 and Aaru-DAB*357 was lower (0.057) than the mean nucleotide p-distances of other scaffold comparisons (p-distance: 0.077–0.084; Fig. 4A). The same between scaffold difference was observed for the mean amino-acid pairwise distance (Aaru-DAB*554 and Aaru-DAB*357 p-distance: 0.093; other scaffolds p-distance: 0.123–0.134; Table S3).

Then we set out to compare the allelic distances on the two scaffolds Aaru-DAB*554 and Aaru-DAB*357 in more detail. First, we compared the nucleotide p-distance between each annotated gene on scaffold Aaru-DAB*357 and each gene on the other four larger scaffolds. The nucleotide distance per gene was always smaller when compared to genes on scaffold Aaru-DAB*554 than compared with genes on the three other scaffolds (Fig. 4B). Second, we compared the nucleotide p-distance between each annotated gene on scaffold Aaru-DAB*554 and each gene on the four larger scaffolds. Likewise, the mean nucleotide p-distance between each annotated gene on scaffold Aaru-DAB*554 was smaller when compared to genes on scaffold Aaru-DAB*357 than compared with genes on the three other scaffolds (Fig. S2). Both the average p-distance of all genes per scaffold and the p-distance per gene are smaller between scaffolds Aaru-DAB*554 and Aaru-DAB*357 compared with the other larger scaffolds, suggesting that scaffolds Aaru-DAB*554 and Aaru-DAB*357 represent both haplotypes of the same genomic region.