MHC class II DQB diversity in the Japanese black bear, Ursus thibetanus japonicus.

BACKGROUND: The major histocompatibility complex (MHC) genes are one of the most important genetic systems in the vertebrate immune response. The diversity of MHC genes may directly influence the survival of individuals against infectious disease. However, there has been no investigation of MHC diversity in the Asiatic black bear (Ursus thibetanus). Here, we analyzed 270-bp nucleotide sequences of the entire exon 2 region of the MHC DQB gene by using 188 samples from the Japanese black bear (Ursus thibetanus japonicus) from 12 local populations.
RESULTS: Among 185 of 188 samples, we identified 44 MHC variants that encoded 31 different amino acid sequences (allotypes) and one putative pseudogene. The phylogenetic analysis suggests that MHC variants detected from the Japanese black bear are derived from the DQB locus. One of the 31 DQB allotypes, Urth-DQB*01, was found to be common to all local populations. Moreover, this allotype was shared between the black bear on the Asian continent and the Japanese black bear, suggesting that Urth-DQB*01 might have been maintained in the ancestral black bear population for at least 300,000 years. Our findings, from calculating the ratio of non-synonymous to synonymous substitutions, indicate that balancing selection has maintained genetic variation of peptide-binding residues at the DQB locus of the Japanese black bear. From examination of genotype frequencies among local populations, we observed a considerably lower level of observed heterozygosity than expected.
CONCLUSIONS: The low level of observed heterozygosity suggests that genetic drift reduced DQB diversity in the Japanese black bear due to a bottleneck event at the population or species level. The decline of DQB diversity might have been accelerated by the loss of rare variants that have been maintained by negative frequency-dependent selection. Nevertheless, DQB diversity of the black bear appears to be relatively high compared with some other endangered mammalian species. This result suggests that the Japanese black bears may also retain more potential resistance against pathogens than other endangered mammalian species. To prevent further decline of potential resistance against pathogens, a conservation policy for the Japanese black bear should be designed to maintain MHC rare variants in each local population.

Background

In Japan, the Asiatic black bear (Ursus thibetanus) is classified as the subspecies (Ursus thibetanus japonicus). The Japanese black bear currently inhabits two main islands in Japan, Honshu and Shikoku. Historically, these bears were also distributed in Kyushu, but no native bears have been captured there for several decades [1,2]. However, there have been several reports of sightings of an animal believed to be a black bear in Kyushu [3]. Research on the genetic diversity in the black bear has begun in recent years and has been developed with the use of genetic markers such as mitochondrial DNA (mtDNA) sequences and microsatellites [4-12]. Using mtDNA markers we have found that the Japanese black bear is genetically quite distinct from other Asiatic black bears [8]. Sequences of the mtDNA control region indicated that the Japanese black bear could be divided genetically into three groups: western, eastern, and Shikoku/Kii-hanto populations [7,8]. Analyses of mtDNA and microsatellite markers [6,8] revealed lower genetic diversity and higher genetic differentiation in the western population than in the other populations. In fact, in most of the management units in western Japan the species has been designated as “endangered” in the Red Data Book of Japan [13] because of the fragmentation and isolation of their habitats.

The major histocompatibility complex (MHC) genes play an important role in the resistance to various pathogens, especially in wild animals. If MHC polymorphism declines, resistance to infectious diseases will decrease [14]. The MHC genes are thought to influence the survival of species or populations against pathogens. Therefore, MHC genetic variability in local populations is regarded as a potentially important index in conservation genetics.

In the human genome, the MHC region is located on the short arm of chromosome 6 at 6p21.3 and covers roughly 4 Mbp comprising 224 genes [15]. The human MHC region encodes at least three functionally and structurally different types of molecules, class I, class II, and class III. The MHC class II region is composed of five sub-regions (DP, DQ, DR, DM, and DO). Each sub-region contains at least one set of A and B genes, and these encode two polypeptide chains, α and β chains, respectively, which make up a heterodimeric molecule. For example, in a DQ molecule, the DQA gene encodes the α-chain, while DQB encodes the β-chain. The DP, DQ, and DR molecules are primarily important for the presentation of peptides and they are classified as classical class II molecules.

For the last three decades, studies of MHC allele variations have been performed not only in humans and laboratory animals but also in populations of wild mammals. Many such studies of wild mammals have demonstrated a relationship between low levels of MHC diversity and decreases in species population fitness (see below). MHC diversity was examined in pinniped species, which are believed to be closely related phylogenetically to ursine species [16,17]. Hoelzel and colleagues [18] reported that genetic variation at the DQB locus in the southern elephant seal (Mirounga leonina) did not appear to be low, whereas the northern elephant seal (Mirounga angustirostris), which was thought to have experienced a severe population bottleneck, exhibited much less variation at this locus. Lento et al. [19] and Weber et al. [20] demonstrated low genetic diversities at the MHC loci in the New Zealand sea lion (Phocarctos hookeri), and the northern elephant seal. The census sizes of these pinniped species have also diminished in the recent past probably due to human activity. In the case of non-carnivore mammals, heterozygosity of the DRB gene in a population of the striped mouse (Rhabdomys pumilio), has been shown to influence infection status [21]. A low level of diversity at the DAB locus (an MHC class II gene in marsupials) in the western barred bandicoot (Perameles bougainville) may influence its vulnerability against novel pathogens [22].

Many studies have revealed a correlation between MHC diversity and resistance against infectious disease; however, others have demonstrated that this is not the case universally. For example, genetic variation in DQ and DR sequences in California sea lions (Zalophus californianus), which is a thriving species, is lower than that in several other mammalian species [23,24]. The Tasmanian devil (Sarcophilus harrisii), which is an endangered species, has low MHC diversity, but individuals with restricted MHC repertoires (These individuals possessing only one of two phylogenetic groups) may be resistant to the devil facial tumor disease [25]. However, their MHC variants were not assigned to loci because of a lack of genomic information. Thereafter, the research group has shown that there is no direct evidence that MHC genetic difference between tumor-free and infected animals is associated with the resistance against the infection although certain genotypes at two MHC-linked microsatellite markers seem to be more frequently detected from healthy animals [26]. The DRB locus in the mountain goat (Oreamnos americanus) exhibited low genetic diversity, but this does not appear to have affected the resistance to infectious diseases in this species [27].

For Ursidae, studies of MHC diversity are somewhat limited. In the giant panda (Ailuropoda melanoleuca), MHC studies are relatively advanced and have been used to propose some possible implications for captive management [28-34]. Low genetic diversity at the AimeDRB and AimeDQA loci has been observed in the giant panda [28,30]. Wan et al. [31] attempted to understand the genomic structure and evolutionary history of Aime-MHC genes. The study has showed that the giant panda possesses a single DYB-like pseudogene (the DY subregion was thought to be ruminant-specific MHC genes, but may be shared by members of Laurasiatheria superorder.). Chen et al. [34] examined genetic diversity at seven Aime-MHC class II genes from 121 individuals and showed the presence of two monomorphic and five polymorphic loci. In brown bears (Ursus arctos) three DQA variants and 19 DRB variants were identified in 32 and 38 bears, respectively, although loci have not been assigned [35,36]. Kuduk and colleagues reported the expression of three MHC class I, two DRB, two DQB and one DQA loci in the brown bear [37]. Other than these studies on the giant panda and brown bear, there have been no published reports on the MHC genetic diversity in the Ursidae.

Although the direct contribution of MHC genes to the immunological fitness of animals is still unclear, the relatively high genetic diversity of the MHC reflects the abundance in repertoire of binding peptides, derived from pathogens. We thus believe that MHC diversity potentially contributes to resistance against environmental pathogens and species survival, indicating that study of these loci is essential for designing a plan for conservation management of animals. Because of large variability of MHC, the examination of MHC diversity is also useful for the management of endangered wild-animal populations that show no variation in genetic markers such as microsatellite loci [38]. In our previous study, we reported a functional DQB sequence (1,150 bp) in the Japanese black bear, identified using the rapid amplification of cDNA ends method [39]. Here, we investigated polymorphism in exon 2 of DQB that contains the hypervariable region involved in binding antigens for the Japanese black bear to consider the potential implications for resistance against pathogens and estimate the status of local populations on the basis of their MHC variation.

Results

Features of the DQB exon 2 variants in the Japanese black bear

The entire exon 2 region of the DQB locus (270 bp) was determined for the Japanese black bear. The amplified region corresponded to amino acids 6 to 94 of the β chain of the human DR molecule [40,41]. The frequency of variant typing success was lowest for tooth samples (Ca. 54%), and was highest for blood (Ca. 92%). This result probably reflects the difference in the amount of gDNA in each sample type. From 188 bear samples in 12 of 19 conservation and management units designed by Yoneda [42] (Figure 1), we detected 47 DQB variants. However, the application of criteria for genotyping used in this study (see Methods) reduced the number to 44 variants [DDBJ: AB634514–AB634536] and to 185 bear samples (Figures 1, 2 and Table 1). In the present study, no more than two different variants per individual were detected from the clones of the PCR products.

Figure 1

Geographical distribution of DQB allotypes and haplotypes of mtDNA control region. The inner pie chart indicates the observed frequencies of DQB allotypes in the Japanese black bear. The outer pie chart indicates and the observed haplotype frequencies of the mtDNA control region [8]. Twelve local populations based on the conservation and management units for the Japanese black bear [42] are shown as bold circles. The number of samples for DQB allotypes from each location is depicted in the center of each corresponding chart, and the number in parenthesis indicates the sample size used for mtDNA analysis of our previous study [8]. Superscript a represents the number of samples analyzed by Ishibashi and Saitoh [5] in addition to our previous study [8]. Superscript b indicates the number was obtained by identification of individuals using microsatellite analysis [90]. Superscript c indicates the putative pseudogene. Superscript d indicates haplotypes of mtDNA identified by Ishibashi and Saitoh [5].

Figure 2

Alignment of segregating sites among the exon 2 variants in the Japanese black bear. Dots indicate identity with the nucleotides of Urth-DQB*0101. The first two numeric characters after an asterisk in a variant name represent an allotype

Table 1

The observed number of nucleotide variants detected from each conservation and management unit

Conservation and management unit	CS	GAI	SO	EM	SA	NCA	HO	NK	EC	WC	KH	SK	Total
No. of sample	5	21	1	13	1	12	36	23	8	53	3	9	185
No. of variant	2	4	1	7	1	10	18	8	7	10	3	6
Urth-DQB*0101	8	29	2	10		8	34	23	6	86	2	9	217
Urth-DQB*0102										2			2
Urth-DQB*0103					2								2
Urth-DQB*0104									2				2
Urth-DQB*0105								2					2
Urth-DQB*0106										2			2
Urth-DQB*0107										1			1
Urth-DQB*0108													1
Urth-DQB*0109										1			1
Urth-DQB*0201		2		5		5	4	2			2	2	22
Urth-DQB*0202						1	8	4			2		15
Urth-DQB*0301				1		2	1			8		2	14
Urth-DQB*0302	2	2		2			5						11
Urth-DQB*0401				2			2	12	3				19
Urth-DQB*0402								1					1
Urth-DQB*0403							1						1
Urth-DQB*0501							2						2
Urth-DQB*0601							2						2
Urth-DQB*0701		9		5		2			1			2	19
Urth-DQB*0801						2							2
Urth-DQB*0901										2			2
Urth-DQB*1001									2				2
Urth-DQB*1101												2	2
Urth-DQB*1201							2						2
Urth-DQB*1301							2						2
Urth-DQB*1401							2						2
Urth-DQB*1501										2			2
Urth-DQB*1601							2						2
Urth-DQB*1701				1									1
Urth-DQB*1801							1						1
Urth-DQB*1901								1					1
Urth-DQB*2001								1					1
Urth-DQB*2101						1							1
Urth-DQB*2201						1							1
Urth-DQB*2301						1							1
Urth-DQB*2401						1							1
Urth-DQB*2501							1						1
Urth-DQB*2601									1				1
Urth-DQB*2701							1						1
Urth-DQB*2801							1						1
Urth-DQB*2901										1			1
Urth-DQB*3001												1	1
Urth-DQB*3101									1				1
Urth-DQB*3201										1			1

CS: Chokai Sanchi; GAI: Gassan-Asahi Iide; SO: Southern Ouu; EM: Echigo-Mikuni; SA: Southern Alps; NCA: North-Central Alps; HO: Hakusan-Okumino; NK: Northern Kinki; EC: Eastern Chugoku; WC: Western Chugoku; KH: Kii Hanto; SK: Shikoku.

Thirty-one distinct amino acid sequences (Urth-DQB*01 to Urth-DQB*31) were translated from nucleotide sequences of the 43 variants (Urth-DQB*0101 to Urth-DQB*3101) (Figure 3). Because one variant examined (Urth-DQB*3201) included a premature stop codon in the sequence, it was recognized as a pseudogene (Figure 3A). This Urth-DQB*3201 was detected from one individual which was heterozygous with this variant and Urth-DQB*0101. An NCBI BLAST search showed that sequence similarities between bear Urth-DQB and panda Aime-DQB variants were 84% to 91% and 93% to 96% similarity at the amino acid and nucleotide levels, respectively [34]. Of all the variants, Urth-DQB*0101, *0104, *0401, *1001, and *2601 correspond to the Urth-DQB*a to Urth-DQB*e variants detected in our previous study [39], at the DNA level. Thus, in the present study we additionally identified 39 novel variants at the DNA level in 178 samples. The nucleotide sequence of Urth-DQB*0101 was identical to that found in exon 2 of a cDNA sequence derived from mRNA of a bear sample in our previous study [39]. This suggests the amplification of a functional DQB locus in the present study, although there is no evidence of protein translation. Deletions and insertions were not detected among the sequences. There were 58 segregating sites among the 44 different nucleotide sequences (Figure 2). Amino acid substitutions were observed at 30 sites among the 31 allotypes (Figure 3). While 11 of these substituted sites were located in the putative peptide-binding residues (PBRs) [40], 19 corresponded to non-PBR positions. The presence of more substitution sites in the non-PBRs appears to suggest a loss of the open reading frame in the bear DQB gene. However, there were no frameshift or nonsense mutations in Urth-DQB sequences with the exception of the Urth-DQB*3201 although our study was limited to exon 2 (Figure 3). In addition, amino acids at 8 of 12 PBR codons that showed a monomorphic pattern were also conserved in the human DQB1 locus, and most of the segregating sites at the non-PBRs among bear DQB variants were singletons. Although we cannot completely exclude the possibility that pseudogenes are present in our data set, the fact that our results show more substitution sites in the non-PBRs than PBRs indicates that it is unlikely that we amplified pseudogenes. Therefore, we defined all bear DQB variants as the putative functional DQB molecule with the exception of Urth-DQB*3201. Thirty-one DQB allotypes could be classified into 15 types of PBR (i.e., 15 different amino acid sequences of the PBRs; Figure 4). In 184 bears (excluding the individual with the putative pseudogene Urth-DQB*3201), 57 of 64 heterozygotes possessed a combination of different types of PBR sequences.

Figure 3

Amino acid sequences and the variability in the exon 2 of bear s. (A - B) Gray boxes indicate positions of the putative peptide binding residues (PBRs) deduced from the amino acid position in the human DR molecule, and the numbers across the sequence correspond to amino acid positions based on the chain of the human DR structure [40]. A filled asterisk represents the putative site under balancing selection identified by both of PAML and OmegaMap programs, while an open asterisk is that by the OmegaMap only. (A) Alignment of the predicted amino acid sequence translated from nucleotides of MHC DQB variants in the Japanese black bear. Dots indicate identity with amino acids of Urth-DQB*01. Urth-DQB*32 includes the putative premature stop codon within the sequence. (B) The variability level at amino acid residues among DQB variants of the Japanese black bear. The variability level was estimated by a Wu-Kabat plot [85]

Figure 4

Neighbor-joining tree and network based on amino acids in the PBR of allotypes. (A) A neighbor-joining (NJ) tree was constructed by using the number of differences method. DQB sequences of the giant panda are used as outgroups. Only bootstrap values over 50% are shown in this figure. (B) Neighbor-Net network is also based on the sequences of NJ tree

On the basis of the nucleotide sequences of the black bear DQB variants, we performed a search for intragenic recombination by using the GENECONV program [43]. Permutation testing indicated the possibility of recombination between some pairs of variants (Urth-DQB*0101, Urth-DQB*0104, Urth-DQB*0109, Urth-DQB*0402, Urth-DQB*0601, Urth-DQB*0701, Urth-DQB*0901, Urth-DQB*1101, Urth-DQB*1301, Urth-DQB*2001, Urth-DQB*2501, Urth-DQB*2701, Urth-DQB*2901, Urth-DQB*3001, Urth-DQB*3101, and Urth-DQB*3201). However, the Bonferroni-corrected BLAST-like test in GENECONV indicated no significant values (P = 1.0) for recombination events. In addition, not only the method described by Satta [44] but also the RDP4 program [45] indicated no significant signals for recombination or gene conversion. Thus, we used all the DQB variants of the Japanese black bear for further analyses.

Gene tree based on exon 2 sequences of ursine MHC class II genes

A maximum likelihood (ML) tree based on amino acid sequences of partial DQB and DRB exon 2 from the Ursidae species was constructed (Figure 5). In this ML tree, the DQB sequences formed a clearly separated cluster from the DRB sequences, in agreement with Yasukochi et al. [39]. These trees showed that the 31 sequences from the Japanese black bear formed a monophyletic group with the DQBs of the giant panda. The Bayesian and Neighbor-joining (NJ) [46] trees also showed a monophyletic clade of the black bear DQBs and panda DQBs (see Additional File 1: Figure S1).

Figure 5

Maximum likelihood tree for the exon 2 region of MHC class II genes. This tree is constructed based on 79 amino acids in the MHC class II DRB and DQB exon 2 region. The genetic distances are computed by the JTT model with gamma distribution. Only bootstrap values over 50% are shown in this figure. Numbers in parentheses are Genbank accession numbers

In the ML tree, DRBs of some ursine species were mixed with each other although the orthology of their sequences was unclear, indicating a possible trans-specific mode in the DRBs of ursine species. In the black bear/panda DQB clade, twenty black bear DQB allotypes (Urth-DQB*02 to *04, *06 to *09, *11 to *14, *17, *19 to *22, *26 to *28, and *30) appear to cluster more closely with panda rather than other black bear DQBs, suggesting the trans-species polymorphism. However, the bootstrap value is low (41%), and the pairwise distance among the twenty bears and other bear DQB variants showed closer relationships than those between black bear and panda DQBs. For example, the pairwise distances between the Urth-DQB*04 and other black bear DQBs ranged from 0.01 to 0.08 (mean value, 0.06), while those between the Urth-DQB*04 and panda DQBs ranged from 0.09 to 0.16 (mean value, 0.13). In addition, the Bayesian tree showed that DQB sequences of Japanese black bear formed a distinct clade from those of giant panda (see Additional File 1: Figure S1). Therefore, the allelic lineage of twenty black bear DQBs and panda DQBs may be different even though the ML tree appears to show the pattern of trans-species evolution.

Selection of the DQB exon 2 region

We estimated the mean numbers of synonymous and nonsynonymous substitutions within the 267-bp region of DQB exon 2 variants in the Japanese black bear by using a modified Nei–Gojobori method [47] (Table 2). For the black bear, the ratio of nonsynonymous to synonymous substitutions (dN/dS) was approximately 3.83 in the PBRs, 0.63 in the non-PBRs, and 1.48 in the entire exon 2 region. Although three statistical tests did not indicate recombination among all the DQB variants, according to the BLAST-like test without the Bonferroni-correction in GENECONV program (see above), 16 variants were shown to have possibly undergone recombination or gene conversion. Recombination would change the dN/dS value, however, even if these possible recombinants were removed, ratios of dN/dS showed slightly higher values: 4.20 in the PBRs, 0.58 in the non-PBRs, and 1.52 in the entire exon 2 region. This shows that inclusion of the possible recombinants does not affect the result. To evaluate the dN/dS of the black bear, we compared them with those of other mammalian DQB or DRB loci (Table 2). In the family Ursidae, the dN/dS ratio within the PBRs of DQB in the black bear was higher than that of the giant panda DQBs, while it was lower than that of brown bear DRBs (but these DRB variants are derived from at least two loci). Compared with other mammals, the dN/dS value of the black bear did not appear to be low. In addition, by using the PAML program with models that describe various selection pressures (M0 to M8) [48-50], ratios of ω were estimated for black bears (Table 3). Values of ω in the black bear were significantly higher than unity (log likelihood ratio test) under M2a and M8 models that allow for balancing selection (ω2 = 7.0 in M2a and ω = 5.7 in M8). The likelihood ratio test of hypotheses (M2a vs M1a and M8 vs M7) showed statistically significant values (P < 0.001), suggesting a balancing selection model is a better fit to the data compared to other models. Thus, the present results indicate that exon 2 of DQB in the black bear has been a target of balancing selection. The Bayesian analysis also showed that five codons were potentially under balancing selection. The mean value of ω, estimated with Bayesian inference by using the OmegaMap program [51], was also larger than unity (ω = 1.8). This analysis also suggested that 11 codons were probably under balancing selection. Of these 11 codons, five were estimated by both the PAML and the OmegaMap programs. A majority of the selected sites (β9, β28, β30, and β37) were located in the predicted PBRs (Figure 3).

Table 2

Comparison of non-synonymous () and synonymous () substitution rates at the exon 2 region

Species	Gene	No. of sample	No. of variant	Length	PBRs				Non-PBRs				Entire exon 2 region
					dN/dxS	dN			dN/dS	dN			dN/dS	dN			dS
Ursidae
Japanese black beara	DQB	184k	43k	267 bp	3.83	0.11	±	0.04	0.63	0.02	±	0.01	1.48	0.04	±	0.01	0.03	±	0.01
Giant pandab	DQB1	121	6	267 bp	1.33	0.02	±	0.06	0.47	0.01	±	0.02	0.71	0.02	±	0.05	0.01	±	0.03
Brown bear*c	DRB	38	19	267 bp	4.09	0.06	±	0.29	0.72	0.02	±	0.05	1.53	0.02	±	0.07	0.02	±	0.11
Other carnivores
African wild dogd	DRB	320	19	267 bp	5.28	0.06	±	0.21	0.56	0.01	±	0.02	1.72	0.02	±	0.04	0.02	±	0.07
Namibian cheetah*e	DRB	139	4	237 bp	2.01	0.06	±	0.25	0.37	0.01	±	0.05	0.76	0.04	±	0.13	0.02	±	0.10
European minkf	DRB	20	9	228 bp	1.78	0.06	±	0.16	0.26	0.01	±	0.02	0.65	0.03	±	0.09	0.02	±	0.06
Non-carnivores
Water buffalo*g, h	DQB	210	22	243 bp	2.65	0.07	±	0.37	0.92	0.03	±	0.13	1.34	0.03	±	0.14	0.03	±	0.19
Ma's night monkeyi	DQB1	19	14	246 bp	2.71	0.05	±	0.19	0.95	0.02	±	0.07	1.39	0.02	±	0.07	0.02	±	0.10
Yellow necked mousej	DRB	146	27	216 bp	7.75	0.09	±	0.38	1.13	0.02	±	0.06	2.46	0.02	±	0.05	0.02	±	0.12

The dN and dS are estimated by the modified Nei-Gojobori method.

a: present study; b: Chen et al. [34]; c: Goda et al. [36]; d: Marsden et al. [91];

e: Castro-Prieto et al. [92]; f: Becker et al. [93]; g: Niranjan et al. [94]; h: Sena et al. [95];

i: Diaz et al. [96]; j: Meyer-Lucht and Sommer [97]; k: An individual that has the putative pseudogene, Urth-DQB*3201, is excluded.

*: The number of MHC variants may be derived from multiple loci.

Table 3

Inference of balancing selection for 43 black bear variants

Method	Model code	Parameter estimates	Likelihood or confidential interval	Positively selected sites
Maximum likelihood (PAML)	M0 (one ratio)	ω = 0.41	−1570
	M1a (nearly neutral)	p0 = 0.85, p1 = 0.15	−1386
	M2a (balancing selection)	ω2 = 7.00	−1326	β9**, β26**, β28**, β30**, β37**
	M3 (discrete)	ω2 = 7.98	−1325
	M7 (beta)	p = 0.05, q = 0.07	−1394
	M8 (beta and ω)	ω = 5.71	−1337	β9**, β26**, β28**, β30**, β37**
Bayesian inference (OmegaMap)	−	ω = 1.80a	ω: 0.30 - 12.49a(95% HPD)	β9**, β22**, β26**, β28*, β29**, β30**, β37*, β50*, β57*, β68*, β84*

The putative pseudogene (Urth-DQB*3201) is excluded from this analysis. ω = dN / dS. p0: the proportion of conserved sites with 0 <ω < 1. p1 (=1 - p0): the proportion of neutral sites. Values of p and q are parameters of the beta distribution. *: posterior probability, P < 0.95; **: posterior probability, P < 0.99. a: a mean value. Bold characters indicate positions of the putative peptide binding residues (PBRs) deduced from the amino acid position in the human DR molecule [40].

Estimation of the divergence time of DQB variants in the Japanese black bear

To estimate the divergence time of DQB variants in the Japanese black bear, we calculated the maximum genetic distance at synonymous sites that corresponded to the time to most recent common ancestor (TMRCA). The resulting value was 0.08 ± 0.03, and, by assuming 3.5 × 10-9 per site per year as the average synonymous substitution rate for mammalian nuclear DNA [52], the divergence time of the bear DQB variants was estimated to be about 11 million years ago (MYA). Even if the possible recombinants described above were removed, the result did not change.

Geographical distribution and genetic variation in the MHC exon 2 variants in the conservation and management units

Five DQB allotypes (Urth-DQB*01, Urth-DQB*02, Urth-DQB*03, Urth-DQB*04, and Urth-DQB*07) were shared among the management units (Figure 1 and Table 1). In particular, Urth-DQB*01 was observed throughout all the management units. This allotype, Urth-DQB*01, corresponds to six nucleotide DQB variants (Urth-DQB*0101 to Urth-DQB*0109). The variant Urth-DQB*0101 was detected in a particularly large number (133) of bears. This suggests that Urth-DQB*01 may have been dominantly maintained in the black bear population in recent times. This is probably because of its important immune role against pathogens or for demographic reasons.

On the basis of the variant frequency in each of the six management units from which more than 10 individuals were sampled, the expected heterozygosity (He) ranged from 0.33 to 0.85 (mean He, 0.65) (Table 4). Values of observed heterozygosity (Ho) ranged from 0.19 to 0.69 (mean Ho, 0.42), which were lower than the He for all six units. This suggests non-random mating and/or population subdivision in each unit or the existence of null alleles. Among the six units, He and Ho were the highest in the North-Central Alps (NCA) unit (He = 0.85) and the Echigo-Mikuni (EM) unit (Ho = 0.69), respectively, while those of Western Chugoku (WC) were the lowest (He = 0.33, Ho = 0.19). In the NCA unit, the Ho value (Ho = 0.33) was significantly low (P < 0.001) compared with the highest He value (He = 0.85). The observed and expected heterozygosities of bears in western Japan (He = 0.54, Ho = 0.31) were lower than that of bears in eastern Japan (He = 0.73, Ho = 0.41). Nucleotide diversity (πS) based on synonymous substitutions within a unit ranged from 0.005 to 0.024 among the six units, and the overall value for black bears was 0.016. The πS value of the WC unit (0.005) was the lowest among the six units.

Table 4

Gene diversities and nucleotide diversities for each conservation and management unit

		MHC DQB exon 2 region (267 bp)					mtDNA control region (Yasukochi et al. 2009)
	Conservation and management unit	No. of samples	No. of alleles	He	Ho	πS	No. of samples	No. of haplotypes	h	π
Eastern Japan	Gassan-Asahi Iide (GAI)	21	4	0.48 ± 0.08	0.33	0.018 ± 0.009	6	5	0.93 ± 0.12	0.011 ± 0.008
	Echigo-Mikuni (EM)	13	7	0.79 ± 0.05	0.69	0.024 ± 0.011	16	7	0.82 ± 0.07	0.009 ± 0.006
	North-Central Alps (NCA)	12	10	0.85 ± 0.05*	0.33*	0.022 ± 0.010	23	8	0.85 ± 0.05	0.011 ± 0.007d
	Hakusan-Okumino (HO)	36	18	0.76± 0.05*	0.39*	0.019 ± 0.008	74	8	0.72 ± 0.03	0.010 ± 0.006
Western Japan	Northern Kinki (NK)	23	8	0.68 ± 0.06	0.56	0.016 ± 0.008	96 (64)	4	0.75 ± 0.03	0.011 ± 0.007
	Eastern Chugoku (EC)	8	7	0.83±0.07*	0.38*	0.013 ± 0.007	31 (22)	2	0.25 ± 0.10	0.003 ± 0.002
	Western Chugoku (WC)	53	11	0.33 ± 0.06*	0.19*	0.005± 0.003	154 (33)	6	0.56 ± 0.04	0.005 ± 0.004
	Shikoku (SK)	9	6	0.73 ± 0.10*	0.33*	0.021 ± 0.010	11c	2	0.55 ± 0.07	0.009 ± 0.006
	Overall	175	44	0.65 ± 0.03*	0.36*	0.016 ± 0.007	455	24	0.90 ± 0.01	0.013 ± 0.008

Bold character indicates conservation and management units which have more than 10 samples. Chokai Sanchi, Southern Ouu, Southern Alps and Kii Hanto units are excluded due to small number of samples. He: the expected heterozygosity; Ho: the observed heterozygosity. πS or π: the nucleotide diversity. h: the haplotype diversity. *: a significant departure from Hardy-Weinberg equilibrium (P < 0.01). ( ): the number of haplotype data analyzed by [5]. a: the number includes subhaplotypes. b: the number was taken from the result of an identification of individuals using microsatellite analysis [90] because samples were obtained by hair trap. c: the number is samples that are newly added to samples used in our previous study [8]. d: the π value described in [8] (π = 0.001) was incorrect, and 0.01 was a correct value.

We calculated four statistics of genetic differentiation (KST, FST, Hedrick’s G’ST and Jost’s D) for pairs of management units by using the nucleotide sequences of DQB variants (Table 5). Among the six management units, most KST and FST values indicated statistically significant levels of genetic differentiation, including those between the WC and other units (mean KST, 0.12; mean FST, 0.17), and the Gassan-Asahi Iide (GAI) and other units (mean KST, 0.05; mean FST, 0.08). With respect to Hedrick’s G’ST and Jost’s D, these values trended to be high between the WC and other units. The high values of G’ST and D (more than 0.05) and significant values of KST and FST were observed only between the WC and EM/NCA/Hakusan-Okumino (HO) units (mean KST, 0.13; mean FST, 0.20; mean G’ST, 0.12; mean D, 0.07). We calculated average values for each of four statistics between management unit pairs in western and eastern Japan. The results indicated that the values for western Japan (KST = 0.05, FST = 0.10, G’ST = 0.05, D = 0.03) were higher than those for eastern Japan (KST = 0.02, FST = 0.05, G’ST = 0.02, D = 0.01).

Table 5

Four statistics of genetic differentiation for variants between the pairs of management units

	Eastern Japan					Western Japan
Conservation and management unit		Gassan-Asahi Iide (GAI)	Echigo-Mikuni (EM)	North-Central Alps (NCA)	Hakusan-Okumino (HO)	Northern Kinki (NK)	Eastern Chugoku (EC)	Western Chugoku (WC)
Eastern Japan	GAI	KST / FST/ Hedrick's G'ST/ Jost's D
	EM	0.040*/ 0.068*/ 0.036**/ 0.022**
	NCA	0.057**/ 0.112**/ 0.049**/ 0.030**	-0.006/ -0.012/ -0.022**/ -0.015**
	HO	0.030**/ 0.063**/ 0.078*/ 0.049*	0.008/ 0.029*/ -0.002*/ -0.002**	-0.001/ 0.026*/ -0.013**/ -0.010**
Western Japan	NK	0.057**/ 0.104**/ 0.063*/ 0.037*	0.017/ 0.054*/ 0.016**/ 0.011**	-0.002/ 0.069**/ 0.021/ 0.015**	-0.001/ 0.029**/ 0.023**/ 0.016**
	EC	0.009/ 0.116**/ 0.056*/ 0.033**	0.037**/ 0.020/ 0.031**/ 0.021**	0.029/ 0.031/ 0.047*/ 0.033*	0.004/ 0.028*/ 0.062*/ 0.044*	0.006/ 0.011/ 0.002**/ 0.001**
	WC	0.051**/ 0.073*/ 0.026**/ 0.010**	0.151**/ 0.229**/ 0.114*/ 0.064*	0.141**/ 0.242**/ 0.114*/ 0.065*	0.098**/ 0.118**/ 0.132/ 0.076*	0.140**/ 0.163**/ 0.105*/ 0.055*	0.031**/ 0.236**/ 0.097*/ 0.052*
	SK	0.000/ 0.032 / -0.003**/ -0.002**	-0.010/ -0.008 / -0.001**/ 0.000**	0.008/ -0.005/ 0.003**/ 0.002**	0.005/ 0.008/ 0.026**/ 0.018**	0.021/ 0.047*/ 0.029*/ 0.018**	0.016/ 0.025/ 0.025**/ 0.016**	0.077**/ 0.127**/ 0.036**/ 0.017**

From top line, the values in Table indicate the KST and FST, Hedrick's G'ST and Jost's D statistics. Bold italic shows values which are statistically significant for both of KST and FST, with positive values of G'ST and D. Underline represents statistical significant KST and FST with more than 0.05 of G'ST and D. Bold characters indicate conservation and management units which have more than 10 samples. Chokai Sanchi, Southern Ouu, Southern Alps and Kii Hanto units are excluded due to small number of samples. *: P < 0.05; **: P < 0.01 (permutation test, 1,000 replicates).

Discussion

Estimations of trans-species polymorphism in the DQB variants of the Japanese black bear

In this study, we detected 44 DQB variants in 185 Japanese black bears. These 44 nucleotide variants translated to give 31 allotypes and one putative, and the 31 allotypes were used to construct an ML tree to examine relationships among MHC class II sequences of ursine species (Figure 5). The tree showed that the DQB sequences of the Japanese black bear formed a monophyletic group with those of the giant panda. In our previous study, we reported that black bear DQB sequences cluster with DQB sequences of various mammalian species [39]. These results suggest that all MHC sequences of black bears are likely to have been derived from the DQB gene.

The ML tree suggested a pattern of DRB trans-species evolution among ursine species, in agreement with the Goda et al. [36] (Figure 5). The estimated divergence time of black bear DQB variants was approximately 11 ± 4.6 MYA. Previous studies have estimated that a rapid radiation of the genus Ursus has occurred at around 5 to 3 MYA [53,54]. Therefore, even if our estimate overestimated divergence time because of the limited number of synonymous sites, it is possible that the DQB allelic lineage of the Japanese black bear was shared with other species in the genus Ursus. In fact, Urth-DQB*0101 has been detected in the Asiatic black bear in Russia (YY, unpublished observation), although the PCR amplification indicating this was performed with the primer pair used in our previous study [39]. Fossil records suggest that the ancestral Japanese black bear migrated to the Japanese archipelago during the Middle Pleistocene (0.5 MYA - 0.3 MYA) [55]. Therefore, in the Asiatic black bear, including the Japanese black bear, the Urth-DQB*01 lineage may have been maintained for at least 0.3 million years (MY).

If sequences of other ursine species are obtained in the future, it will be important to carefully distinguish orthologous genes from paralogous ones. It is possible that variants of paralogous genes can aberrantly reflect the trans-specific pattern of alleles in orthologous genes.

Selective pressure on the DQB gene of the Japanese black bear

To estimate whether balancing selection have acted on the DQB gene of black bears, we calculated dN/dS (or ω) ratios by using three methods: the modified Nei–Gojobori model, ML, and Bayesian inferences. All three methods indicated that the dN/dS values were significantly higher than unity, supporting the influence of balancing selection on the DQB locus of the Japanese black bear. According to the modified Nei–Gojobori method, the dN/dS value in the PBR was 3.8. This value is relatively higher than those of some other mammalian species (Table 2). The values of the giant panda DQB (1.3), the Namibian cheetah DRB (2.0), and the European mink DRB (1.8) are all lower than that of the black bear. In addition, the number of MHC variants in the black bear is considerably higher than that of the three species even though the variants of the Namibian cheetah may be derived from multiple loci (Table 2). Since these mammalian species are all designated as an endangered species in the International Union for Conservation of Nature (IUCN) Red List [56], the effective population sizes of these species are expected to be quite small. In theory, the effect of selective pressure decreases with diminishing the effective population size, due to a larger effect of random genetic drift [57]. Thus, since the census population size of the Japanese black bear is relatively large, (more than 10,000 individuals; the exact number of individuals is still unclear.), the black bear may show the higher dN/dS value, compared with other endangered species.

Peptide-binding repertoires of the DQB molecules in the Japanese black bear

The amino acids in the PBRs for the Japanese black bear DQBs can be divided into 15 different types (Figure 4). Theoretically, the average number of nonsynonymous changes within the PBRs should be equal to the number of functionally different alleles that segregate in a population [58]. These groups might correspond to so-called “super-types” [59-61]. A super-type is a group of alleles with a similar repertoire of peptide binding. The mean number of nonsynonymous substitution in the PBRs (KB) of the black bear DQB was 10.3, and this translated to an expected number of allelic lineages of approximately ten. This is slightly inconsistent with the number of different types in the PBRs in Figure 4. This suggests that some PBRs may be functionally similar (i.e., a same super-type). Alternatively, since the length of DQB sequences is limited to exon 2 only, the resolution of the KB estimate may be lower.

Effects of demographic events on the Japanese black bear population

We compared results obtained from analyses of the mtDNA [8] and DQB genetic markers to estimate demographic events that have occurred in Japanese black bear populations. MHC alleles are maintained for long periods by balancing selection [62-64], while the coalescence times of mtDNA are short because of its maternal inheritance mode and ploidy level. Due to this characteristic, mtDNA has a higher susceptibility to the effect of recent demographic factors. Therefore, the time scales in evolutionary history represented by these two types of genetic markers are different.

The extent of sharing DQB allotypes among Japanese black bear management units was higher than that of mtDNA haplotypes (Figure 1). In particular, Urth-DQB*01 (Urth-DQB*0101) was frequently found throughout Japan. This observation could reflect the difference in coalescence time between the two systems and a time-dependent difference in the degree of population fragmentation. Alternatively, if the same dominant pathogen affects different populations (i.e., a similar selective pressure), a DQB allotype that confers resistance to this pathogen is expected to be shared by different populations. On the other hand, the management unit-specific allotypes were found from most of the units. This suggests that rare allotypes also have been maintained in the Japanese black bear population, possibly due to negative frequency-dependent selection [64].

Values of observed heterozygosity (Ho) were lower than the He for all six units. Some management units showed a significant departure from Hardy-Weinberg equilibrium (HWE) (Table 4). The excess of He appeared not to be associated with the quality of DNA samples as there was no correlation between the number of homozygotes and the quality of DNA samples (P = 0.90). In addition, to examine the possibility of null alleles, we determined nucleotide sequences of the primer target region in 96 samples (YY, unpublished data). While all individuals showed only one degenerated site in the reverse primer-annealing region, 10 individuals possessed seven degenerated sites in the forward primer region. In the forward primer region of the 10 individuals, there were two sets of nucleotide sequences that differed by less than 7 bp. The seven degenerated sites are derived from these two nucleotide sequences. However, the forward primer used in the present study (URDQBex2-F4) is designed to specifically anneal to only one particular sequence (see Methods) in the primer region. Therefore, the seven degenerated sites are unlikely to be the cause of null allele. This suggests that the excess of He is not caused by technical problems. Since cross-locus amplification potentially decreases the possibility of the He excess (i.e., it increases the Ho values because of false heterozygotes), the He excess suggests that our designed primer amplifies a single DQB locus. Although we cannot completely exclude the possibility of null alleles, the He excess is more likely caused by demographic events such as non-random mating and/or population subdivision in a local population.

In our previous study [8], FST values for mtDNA haplotypes among the management unit pairs in western Japan were larger than those among the pairs in eastern Japan. Since mtDNA haplotype distributions reflect recent demographic events, the large FST in western Japan indicates a more infrequent gene flow within western Japanese populations. The level of genetic differentiation at DQB among the units in western Japan was also high. This increase may be due to non-random mating and genetic drift in western Japan in recent history. Non-random mating and genetic drift may be caused by fragmentation and isolation of habitat.

The genealogical structure of MHC alleles under strong symmetric balancing selection is similar to that of neutral genes, with the exception of a difference in time scales [65]. To estimate the effect of population subdivision in the Japanese black bear in the distant past, we calculated the ratio of the largest genetic distance (TMRCA, dTMRCA) to the mean distance (dS or d) in bear DQB variants. At neutral genes, dTMRCA is expected to be twice the value of d[66]. The ratio of dTMRCA to dS in the DQB variants was approximately 3 (dTMRCA = 0.078 ± 0.032, dS = 0.028 ± 0.013), suggesting an effect of population fragmentation. Thus, the Japanese black bear population might have been isolated and fragmented not only in the recent past but also in the relatively distant past.

Genetic degradation of black bears in the WC unit

For the WC unit in western Japan, values of genetic diversity, for not only the MHC DQBs (He = 0.33, π = 0.005) but also mtDNA (h = 0.56, π = 0.005; [8]), were significantly lower than those for the other management units (Table 4). In addition, both FST and KST values calculated between the WC unit and other units were significantly high, and also Hedrick’s G’ST and Jost’D tended to be higher than those among the other units. This suggests that genetic drift and/or non-random mating also influence MHC diversity in the WC unit, and it implies that black bears in the WC unit might be gradually diminishing in their potential to adapt to environmental adversity, such as the presence of various pathogens. Thus, even if bears in the WC unit currently have resistance against infectious disease, we need to design a conservation plan to maintain MHC diversity in the bears in this unit. The habitat area of the Eastern Chugoku (SK), Kii Hanto (KH), and Shikoku (SK) units in western Japan are smaller than that of the WC unit. Since previous studies that examined genetic diversities of black bears by using mtDNA and microsatellite markers have suggested a lower level of genetic diversities in the populations of western Japan [6,8], MHC variation in these management units might be lower than that in the WC unit. Therefore, it is essential to investigate MHC variation in these units further by using an increased number of samples.

The effects of natural selection and genetic drift on patterns of MHC variation

We observed significant He excess and genetic differentiation among the management units due to population fragmentation. Ohnishi et al. [6] examined genetic diversity of Japanese black bears by using 10 microsatellite makers. In the previous study, the deviation of Ho from He was not statistically significant in all examined local populations by using the HWE test. However, microsatellite markers that were not used in the study of Ohnishi et al. [6] showed such a deviation in the Eastern Chugoku (EC) and WC units in western Japan [4]. In addition, the genetic diversities within these units were similar to that of the isolated population of the brown bear [6]. Sutton et al. [67] proposed that negative frequency-dependent selection may be predominant in pre-bottlenecked populations and that greater loss of MHC diversity compared to neutral genetic diversity is accelerated because of uneven allele distributions. Such skewed MHC allele frequencies can occur under negative frequency-dependent selection or joint action of this selection [67,68]. Indeed, numerous MHC rare variants of the Japanese black bear were observed in the present study. Therefore, Japanese black bears might have lost some DQB rare variants due to a bottleneck event in their evolutionary history.

Comparison of two criteria for MHC genotyping

Due to the possibility of PCR artifacts, criteria for MHC genotyping on the basis of partial sequences (e.g., exon 2 only) remain a matter of debate. For instance, previous studies have identified MHC variants according to more conservative criteria: If a rare sequence is detected in one individual alone (heterozygote) and the sequence differs by less than 3 bp from a known sequence of the same PCR product, the sequence is considered to be an artifact introduced by PCR error [69,70]. Although we attempted to carefully identify DQB variants, the methods might not be definitive. Therefore, we compared the results of analyses between the two criteria. According to our criteria, we detected 44 variants from 185 bears. Application of more conservative criteria reduced the number to 23 variants and to 179 bear samples. Nevertheless, results of analyses such as estimates of dN/dS ratio, genetic differentiation among units and deviation from the HWE did not largely change even though the number of variants was different between two criteria. This suggests that our conclusion is not affected by the differences between criteria.

Future studies

In this study, we detected a maximum of 44 variants including one putative (23 at a minimum). However, there is a possibility that some of these variants are pseudogenized (i.e., contain frameshift or nonsense mutations in other exons or promoter regions), although we confirmed that the nucleotide sequence of exon 2 in a near full-length cDNA sequence was identical to that of the Urth-DQB*01 [39]. Since a majority of our samples was collected from tooth, the size of amplified DNA fragments was limited to the 270 bp of exon 2. In addition, the sample size of our study was not sufficient to investigate MHC diversity in each management unit. To evaluate our results precisely, further analysis based on the sequences of the entire coding region of DQB should be performed by using more samples.

We designed two pairs of primers to amplify the bear DQB exon 2. When one of these pairs (URDQBex2-F1 and -R1) was used, three different sequences were detected from an individual by TA cloning, suggesting that at least two DQB loci were amplified. Wan et al. [31] have reported that the giant panda genome has three DQB loci, two functional and one non-functional. Kuduk et al. [37] have also confirmed expression of two DQB genes in the brown bear. The PCR results of our study suggest that the Japanese black bear may also possess at least two DQB loci in the genome (See Methods). It will be desirable to determine the genomic sequence across the wider MHC region in order to assure the specificity of each MHC locus in the bear genome.

Conclusions

In the present study, we detected 44 MHC variants from 185 Japanese black bears. The phylogenetic analysis suggests that their sequences are derived from the DQB loci. The value of dN/dS ratio indicates that the bear DQB variants have been maintained by balancing selection. The estimation of TMRCA among their variants suggests that the DQB lineages arose before the origin of genus Ursus. At least, the Urth-DQB*0101 variant has been shared between Japanese and continental black bears for 0.3 MY (before the migration of Asiatic black bear from the Asian continent into Japanese archipelago). However, DQB genotyping for the Japanese black bear showed considerably lower Ho values than He values in each local population. One may hypothesize that skewed allelic distribution due to negative frequency-dependent selection and the bottleneck or population fragmentation in each management unit might accelerate loss of DQB rare variants. Therefore, the results of our study imply that MHC variation in Japanese black bears might have declined due to demographic events such as habitat fragmentation and isolation at the population or species level. Nevertheless, DQB genetic diversity in Japanese black bears appears to be relatively high compared to some other endangered mammalian species. This result suggests that for the present the Japanese black bears appear to retain more potential resistance against pathogens than other endangered mammalian species. To prevent further decline of potential resistance against pathogens, a conservation policy for the Japanese black bear should be designed to maintain MHC rare variants in each of management units. To achieve this objective, it is important to stop further fragmentation and isolation of bear habitats.

Methods

Collection of samples from Japanese black bears

In this study, we collected 307 Japanese black bear samples throughout the Japanese archipelago. The samples were collected from individuals that were hunted or captured for pest control. Of these samples, 185 were used for MHC genotyping. These 185 samples were obtained from 12 conservation and management units designated by Yoneda [42] (Figure 1) between 1970 and 2006. Local populations of the black bears in Japan are divided into 19 conservation and management units, taking into account the prevention of bear movement, such as rivers, railways and so on. The division of local populations in our study followed the divisions of the bear management units. Five samples were obtained from the Chokai Sanchi (CS) unit; 21 from the Gassan-Asahi Iide (GAI) unit; one from the Southern Ouu (SO) unit; 13 from the Echigo-Mikuni (EM) unit; one from the Southern Alps (SA) unit; 12 from the North-Central Alps (NCA) unit; 36 from the Hakusan-Okumino (HO) unit; 23 from the Northern Kinki (NK) unit; 8 from the Eastern Chugoku (EC) unit; 53 from the Western Chugoku (WC) unit; 3 from the Kii Hanto (KH) unit; and 9 from the Shikoku (SK) unit (Figure 1). We obtained the following types of samples: muscle and liver tissue, 42 samples; blood, 10; hair, 12; oral mucous membrane, 4; and tooth or mandibula bone, 117. Based on the mtDNA analysis (see above), there were three genetically distinct groups of Japanese black bears. However, the sample sizes of the management units in the Shikoku/Kii-hanto populations were quite small in this study (three samples in the KH unit and nine in the SK unit). Therefore, for the purpose of discussion of population structure, we combined the Shikoku/Kii-hanto population with the western Japanese population, as is the case with our previous study [8].

Experimental procedures for MHC genotyping

DNA extraction was performed using the procedure described by our previous study [8]. To amplify the entire exon 2 of the DQB in the Japanese black bear, we initially used a single primer set URDQBex2-F1 and URDQBex2-R1 [39]. However, when this primer set was used, three variants per individual were detected from several individuals (i.e., at least two loci were amplified). When nucleotide sequences of the three variants were compared, two of three variants showed nearly identical sequences. Therefore, for amplification of a single DQB locus, a pair of internal primers for introns 1 and 2, URDQBex2-F4 (5′-GAGGCTGTGGGTTGGCCTGAGGCGCGGA-3′) and URDQBex2-R4 (5′-CCGCGGGGCGCGAACGGCCTGGCTCA-3′), was designed to specifically amplify the sequences in intron 1 region of their two similar variants. The PCR amplification, identification of heterozygous or homozygous individuals by TA cloning, and sequencing were performed as described by Yasukochi et al. [39]. Murray et al. [71] showed that the probability that both of two variant sequences in a heterozygote are not detected by sequencing five clones (assuming that every variant is equally amplified at an equal frequency) is very low (P = 0.06). In our study, all final PCR products (188) were cloned, and 8 to 22 clones per individual were sequenced to decrease the possibility of false homozygotes caused by amplification bias. In addition, to confirm sequencing reliability, we compared the result of direct-sequencing at the polymorphic sites in uncloned PCR product with cloned sequences (see Additional File 2: Figure S2). Since such a direct-sequencing are derived from the PCR product amplified from an original gDNA, it decreases the possibility of error incurred during cloning or PCR of the plasmid insert after cloning. If the frequencies of two variants detected from a heterozygote were extremely different and a double peak corresponding to the correct nucleotides was not observed after direct sequencing, that individual was excluded from our analyses (see Additional File 2: Figure S2). In addition, we used high fidelity polymerase during all PCRs (Ex Taq Hot Start Version TaKaRa Bio Inc.).

Phylogenetic relationships and population differentiation

The nucleotide sequences were aligned and translated into amino acids by using MEGA v5.0 [72]. The NJ and ML trees were constructed on the basis of Jones-Taylor-Thornton (JTT) model with gamma distribution [73]. The NJ tree construction and estimation of the best fit substitution model were performed by using the same software. The ML tree was estimated by using the PhyML 3.0 [74]. Bootstrap analysis was performed (1,000 replications in the NJ tree and 500 replications in the ML tree). A network of MHC variants was constructed with the Neighbor-Net method [75] by using Splits-tree4 ver. 4.11.3 [76]. A Bayesian phylogenetic analysis was conducted with MrBayes v3.2.1 [77]. The analysis was implemented considering 11,500,000 generations and tree sampling every 100 generations. The first 28,750 trees were discarded as burn-in. The genetic distance was estimated by using the JTT model. The Bayesian and ML trees were visualized with FigTree version 1.3.1 (http://tree.bio.ed.ac.uk/software/figtree/, [78]). The expected heterozygosity (He) for each management unit and Wright’s F-statistics (FST) between the pairs of management units were calculated using Arlequin v3.11 [79]. In addition, we calculated KST[80] that, like FST, is an index of genetic differentiation, using the DnaSP 5.10 [81]. FST is based on frequencies at polymorphic sites, while KST uses the average number of differences between sequences. The FST and KST are possibly misleading regarding the differentiation when heterozygosity is high. Thus, we also calculated Hedrick’s G’ST[82] and Jost’s D[83] using the SMOGD 1.2.5 [84].

Measure of genetic diversity and distance among DQB variants

The extent of amino acid variability was evaluated using a Wu–Kabat plot [85]. The mean numbers of nonsynonymous substitutions per nonsynonymous site (dN) and that of synonymous substitutions per synonymous site (dS) were calculated for excluding identical sequences of DQB variants. We estimated nucleotide diversity (πS) among all the sequences in each management unit or across the entire population. The divergence time of the Japanese black bear DQB variants was estimated using synonymous substitutions alone. We assumed that the maximum genetic distance at synonymous sites roughly corresponded to the TMRCA of DQBs. These genetic distance parameters (πS, dN, and dS) were calculated using the modified Nei–Gojobori method [47] with Jukes–Cantor correction [86] by using MEGA v5.0. The transition/transversion bias in the calculation was assumed as R = 1.58 by the ML method. Standard errors were calculated from 1,000 bootstrapping replicates.

Estimation of natural selection for DQB variants

Ratios of dN/dS were calculated for each of the putative PBRs, non-PBRs, and for the overall exon 2 region. The extent of selective pressure at the amino acid level was measured by the ratio of nonsynonymous to synonymous rates (ω = dN/dS). We estimated ω ratios with the maximum likelihood method by using CODEML in the PAML program version 4.2b [50]. Several models that differed in their hypotheses concerning the statistical distributions of the ω ratio [48,49] were considered in this study. Nested models can be compared in pairs by using a likelihood-ratio test [87]. A Bayesian approach (Bayes empirical Bayes) in CODEML is usually used to predict codons that were potentially under positive selection or balancing selection in models M2a and M8 [49]. Although it is difficult to distinguish between positive selection and balancing selection, we use this Bayesian approach to estimate sites under balancing selection because MHC genes are thought to be affected by balancing selection rather than positive selection [65,88]. The ML method has the potential to overestimate ω values because of the existence of recombination [70]. Therefore, we estimated the ω value by using a Bayesian method with OmegaMap [51], which enabled us to infer balancing selection in the presence of recombination. The average number of nonsynonymous substitutions in the PBRs (KB) was also estimated using method II, as described by Satta et al. [89].

Test for reciprocal recombination or gene conversion

The GENECONV program was used with default settings. The global test for reciprocal recombination or gene conversion events was performed with 10,000 permutations of the sequence alignment to assess significance. Moreover, we examined the possibility of recombination by using the method described by Satta [44]. In the method, it is assumed that the number of substitutions in a particular region to the entire region is binomially distributed. We divided exon 2 into two regions, β sheet and α helix, for the examination of recombination. The RDP4 software was also used for detecting a possible recombination [45].

Abbreviations

MHC: Major histocompatibility complex; PBR: Peptide-binding residue; NJ: Neighbor-joining; TMRCA: Time to most recent common ancestor; MY: Million years; MYA: Million years ago; He: Expected heterozygosity; HO: Observed heterozygosity; ML: Maximum likelihood; JTT: Jones-Taylor-Thornton; HWE: Hardy-Weinberg equilibrium; CS: Chokai Sanchi; GAI: Gassan-Asahi Iide; SO: Southern Ouu; EM: Echigo-Mikuni; SA: Southern Alps; NCA: North-Central Alps; HO: Hakusan-Okumino; NK: Northern Kinki; EC: Eastern Chugoku; WC: Western Chugoku; KH: Kii Hanto; SK: Shikoku.

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

YY carried out the experiments and data analysis, and wrote the draft of manuscript. TK and MY participated in sample collection. HK and YS were the primary supervisors. All authors have read and approved the final manuscript.

Additional file 1

The NJ and Bayesian trees of ursine and genes. The NJ and Bayesian trees were constructed based on amino acid sequences of partial exon 2 region of ursine DQB or DRB genes. The JTT model with gamma distribution is used for the NJ and Bayesian trees, respectively. Numbers in parentheses are Genbank accession numbers. Only bootstrap values over 50% and posterior probabilities over 80% are shown in the NJ and Bayesian trees, respectively.

Click here for file

Additional file 2

The comparison of nucleotides at the polymorphic sites between uncloned and cloned PCR products. (A) An example that an individual is identified as a heterozygote. The upper figure is the alignment of nucleotide sequences of uncloned and cloned PCR products. The lower figure is the result of direct-sequencing for uncloned PCR product. (B) An example that an individual is not identified as a heterozygote. Such an individual is excluded from analyses of the present study.

Click here for file