Skip to main content

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

Abstract

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 [412]. 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 [2834]. Low genetic diversity at the Aime DRB and Aime DQA 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 DQBexon 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
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
figure 2

Alignment of segregating sites among the DQB 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 DQB variants detected from each conservation and management unit

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
figure 3

Amino acid sequences and the variability in the exon 2 of bear DQB 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
figure 4

Neighbor-joining tree and network based on amino acids in the PBR of DQB 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
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 DQBexon 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) [4850], 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 ( d N ) and synonymous ( d S ) substitution rates at the DQB exon 2 region
Table 3 Inference of balancing selection for 43 black bear DQB variants

Estimation of the divergence time of DQBvariants 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

We calculated four statistics of genetic differentiation (KST, FST, Hedrick’s GST 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 GST and Jost’s D, these values trended to be high between the WC and other units. The high values of GST 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 GST, 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, GST = 0.05, D = 0.03) were higher than those for eastern Japan (KST = 0.02, FST = 0.05, GST = 0.02, D = 0.01).

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

Discussion

Estimations of trans-species polymorphism in the DQBvariants 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 DQBgene 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” [5961]. 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 [6264], 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 GST[82] and Jost’s D[83] using the SMOGD 1.2.5 [84].

Measure of genetic diversity and distance among DQBvariants

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 DQBvariants

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

H e :

Expected heterozygosity

H O :

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.

References

  1. Japan Wildlife Research Center: Report on population status of Japanese bears in 1998. 1999, Japan Wildlife Research Center, in Japanese

    Google Scholar 

  2. Ohnishi N, Yasukochi Y: The origin of the last Asian black bear in Kyushu Island. Mamm Sci. 2010, 50: 177-180. in Japanese, with English summary

    Google Scholar 

  3. Kurihara T: Records of recent bear witnesses in Kyushu Island, Japan. Mamm Sci. 2010, 50: 187-193. in Japanese, with English summary

    Google Scholar 

  4. Saitoh T, Ishibashi Y, Kanamori H, Kitahara E: Genetic status of fragmented populations of the Asian black bear Ursus thibetanus in western Japan. Popul Ecol. 2001, 43: 221-227. 10.1007/s10144-001-8186-4.

    Article  Google Scholar 

  5. Ishibashi Y, Saitoh T: Phylogenetic relationships among fragmented Asian black bear (Ursus thibetanus) populations in western Japan. Conservat Genet. 2004, 5: 311-323.

    Article  CAS  Google Scholar 

  6. Ohnishi N, Saitoh T, Ishibashi Y, Oi T: Low genetic diversities in isolated populations of the Asian black bear (Ursus thibetanus) in Japan, in comparison with large stable populations. Conservat Genet. 2007, 8: 1331-1337. 10.1007/s10592-006-9281-z.

    Article  Google Scholar 

  7. Ohnishi N, Uno R, Ishibashi Y, Tamate HB, Oi T: The influence of climatic oscillations during the Quaternary Era on the genetic structure of Asian black bears in Japan. Heredity. 2009, 102: 579-589. 10.1038/hdy.2009.28.

    Article  PubMed  CAS  Google Scholar 

  8. Yasukochi Y, Nishida S, Han S-H, Kurosaki T, Yoneda M, Koike H: Genetic structure of the Asiatic black bear in Japan using mitochondrial DNA analysis. J Hered. 2009, 100: 297-308. 10.1093/jhered/esn097.

    Article  PubMed  CAS  Google Scholar 

  9. Ohnishi N, Yuasa T, Morimitsu Y, Oi T: Mass-intrusion-induced temporary shift in the genetic structure of an Asian black bear population. Mamm Stud. 2011, 36: 67-71. 10.3106/041.036.0204.

    Article  Google Scholar 

  10. Yoneda M, Mano T: The present situation and issues for estimating population size and monitoring population trends of bears in Japan. Mamm Sci. 2011, 51: 79-95. in Japanese, with English summary

    Google Scholar 

  11. Yamamoto T, Oka T, Ohnishi N, Tanaka H, Takatsuto N, Okumura Y: Genetic Characterization of Northernmost Isolated Population of Asian Black Bear (Ursus thibetanus) in Japan. Mamm Stud. 2012, 37: 85-91. 10.3106/041.037.0209.

    Article  Google Scholar 

  12. Uno R, Kondo M, Yuasa T, Yamauchi K, Tsuruga H, Tamate HB, Yoneda M: Assessment of genotyping accuracy in a non-invasive DNA-based population survey of Asiatic black bears (Ursus thibetanus): lessons from a large-scale pilot study in Iwate prefecture, northern Japan. Popul Ecol. 2012, 54: 509-519. 10.1007/s10144-012-0328-3.

    Article  Google Scholar 

  13. Ministry of the Environment: Mammalia. Threatened Wildlife of Japan—Red Data Book. 2002, Japan Wildlife Research Center, 1-in Japanese, 2

  14. Klein J, Takahata N: The major histocompatibility complex and the quest for origins. Immunol Rev. 1990, 113: 5-25. 10.1111/j.1600-065X.1990.tb00034.x.

    Article  PubMed  CAS  Google Scholar 

  15. Man B: Complete sequence and gene map of a human major histocompatibility complex. The MHC sequencing consortium. Nature. 1999, 401: 921-923. 10.1038/44853.

    Article  Google Scholar 

  16. Delisle I, Strobeck C: A phylogeny of the Caniformia (order Carnivora) based on 12 complete protein-coding mitochondrial genes. Mol Phylogenet Evol. 2005, 37: 192-201. 10.1016/j.ympev.2005.04.025.

    Article  PubMed  CAS  Google Scholar 

  17. Yu L, Luan PT, Jin W, Ryder OA, Chemnick LG, Davis HA, Zhang YP: Phylogenetic utility of nuclear introns in interfamilial relationships of Caniformia (order Carnivora). Syst Biol. 2011, 60: 175-187. 10.1093/sysbio/syq090.

    Article  PubMed  CAS  Google Scholar 

  18. Hoelzel AR, Stephens JC, O’Brien SJ: Molecular genetic diversity and evolution at the MHC DQB locus in four species of pinnipeds. Mol Biol Evol. 1999, 16: 611-618. 10.1093/oxfordjournals.molbev.a026143.

    Article  PubMed  CAS  Google Scholar 

  19. Lento GM, Baker CS, David V, Yuhki N, Gales NJ, O’Brien SJ: Automated single-strand conformation polymorphism reveals low diversity of a Major Histocompatibility Complex Class II gene in the threatened New Zealand sea lion. Mol Ecol Notes. 2003, 3: 346-349. 10.1046/j.1471-8286.2003.00445.x.

    Article  CAS  Google Scholar 

  20. Weber DS, Stewart BS, Schienman J, Lehman N: Major histocompatibility complex variation at three class II loci in the northern elephant seal. Mol Ecol. 2004, 13: 711-718. 10.1111/j.1365-294X.2004.02095.x.

    Article  PubMed  CAS  Google Scholar 

  21. Froeschke G, Sommer S: MHC class II DRB variability and parasite load in the striped mouse (Rhabdomys pumilio) in the Southern Kalahari. Mol Biol Evol. 2005, 22: 1254-1259. 10.1093/molbev/msi112.

    Article  PubMed  CAS  Google Scholar 

  22. Smith S, Belov K, Hughes J: MHC screening for marsupial conservation: extremely low levels of class II diversity indicate population vulnerability for an endangered Australian marsupial. Conservat Genet. 2009, 11: 269-278.

    Article  Google Scholar 

  23. Bowen L, Aldridge BM, Gulland F, Woo J, Van Bonn W, DeLong R, Stott JL, Johnson ML: Molecular characterization of expressed DQA and DQB genes in the California sea lion (Zalophus californianus). Immunogenetics. 2002, 54: 332-347. 10.1007/s00251-002-0472-6.

    Article  PubMed  CAS  Google Scholar 

  24. Bowen L, Aldridge BM, Gulland F, Van Bonn W, DeLong R, Melin S, Lowenstine LJ, Stott JL, Johnson ML: Class II multiformity generated by variable MHC-DRB region configurations in the California sea lion (Zalophus californianus). Immunogenetics. 2004, 56: 12-27. 10.1007/s00251-004-0655-4.

    Article  PubMed  CAS  Google Scholar 

  25. Siddle HV, Marzec J, Cheng Y, Jones M, Belov K: MHC gene copy number variation in Tasmanian devils: implications for the spread of a contagious cancer. Proc Biol Sci. 2010, 277: 2001-2006. 10.1098/rspb.2009.2362.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  26. Lane A, Cheng Y, Wright B, Hamede R, Levan L, Jones M, Ujvari B, Belov K: New insights into the role of MHC diversity in devil facial tumour disease. PLoS One. 2012, 7: e36955-10.1371/journal.pone.0036955.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  27. Mainguy J, Worley K, Côté SD, Coltman DW: Low MHC DRB class II diversity in the mountain goat: past bottlenecks and possible role of pathogens and parasites. Conservat Genet. 2007, 8: 885-891. 10.1007/s10592-006-9243-5.

    Article  CAS  Google Scholar 

  28. Zeng CJ, Yu JQ, Pan HJ, Wan QH, Fang SG: Assignment of the giant panda MHC class II gene cluster to chromosome 9q by fluorescence in situ hybridization. Cytogene Genome Res. 2005, 109: 533-

    Article  Google Scholar 

  29. Zeng CJ, Pan HJ, Gong SB, Yu JQ, Wan QH, Fang SG: Giant panda BAC library construction and assembly of a 650-kb contig spanning major histocompatibility complex class II region. BMC Genom. 2007, 8: 315-10.1186/1471-2164-8-315.

    Article  Google Scholar 

  30. Wan QH, Zhu L, Wu H, Fang SG: Major histocompatibility complex class II variation in the giant panda (Ailuropoda melanoleuca). Mol Ecol. 2006, 15: 2441-2450. 10.1111/j.1365-294X.2006.02966.x.

    Article  PubMed  CAS  Google Scholar 

  31. Wan QH, Zeng CJ, Ni XW, Pan HJ, Fang SG: Giant panda genomic data provide insight into the birth-and-death process of mammalian major histocompatibility complex class II genes. PLoS One. 2009, 4: 11-

    Article  Google Scholar 

  32. Zhu L, Ruan XD, Ge YF, Wan QH, Fang SG: Low major histocompatibility complex class II DQA diversity in the Giant Panda (Ailuropoda melanoleuca). BMC Genet. 2007, 8: 29-

    Article  PubMed  PubMed Central  Google Scholar 

  33. Pan HJ, Wan QH, Fang SG: Molecular characterization of major histocompatibility complex class I genes from the giant panda (Ailuropoda melanoleuca). Immunogenetics. 2008, 60: 185-193. 10.1007/s00251-008-0281-7.

    Article  PubMed  CAS  Google Scholar 

  34. Chen YY, Zhang YY, Zhang HM, Ge YF, Wan QH, Fang SG: Natural selection coupled with intragenic recombination shapes diversity patterns in the major histocompatibility complex class II genes of the giant panda. J Exp Zool B Mol Dev Evol. 2010, 314: 208-223.

    PubMed  Google Scholar 

  35. Goda N, Mano T, Masuda R: Genetic diversity of the MHC class-II DQA gene in brown bears (Ursus arctos) on Hokkaido, Northern Japan. Zoolog Sci. 2009, 26: 530-535. 10.2108/zsj.26.530.

    Article  PubMed  CAS  Google Scholar 

  36. Goda N, Mano T, Kosintsev P, Vorobiev A, Masuda R: Allelic diversity of the MHC class II DRB genes in brown bears (Ursus arctos) and a comparison of DRB sequences within the family Ursidae. Tissue Antigens. 2010, 76: 404-410. 10.1111/j.1399-0039.2010.01528.x.

    Article  PubMed  CAS  Google Scholar 

  37. Kuduk K, Babik W, Bojarska K, Sliwinska EB, Kindberg J, Taberlet P, Swenson JE, Radwan J: Evolution of major histocompatibility complex class I and class II genes in the brown bear. BMC Evol Biol. 2012, 12: 197-10.1186/1471-2148-12-197.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  38. Aguilar A, Roemer G, Debenham S, Binns M, Garcelon D, Wayne RK: High MHC diversity maintained by balancing selection in an otherwise genetically monomorphic mammal. Proc Nat Acad Sci USA. 2004, 101: 3490-3494. 10.1073/pnas.0306582101.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  39. Yasukochi Y, Kurosaki T, Yoneda M, Koike H: Identification of the expressed MHC class II DQB gene of the Asiatic black bear, Ursus thibetanus, in Japan. Gene Genet Syst. 2010, 85: 147-155. 10.1266/ggs.85.147.

    Article  CAS  Google Scholar 

  40. Brown JH, Jardetzky TS, Gorga JC, Stern LJ, Urban RG, Strominger JL, Wiley DC: Three-dimensional structure of the human class II histocompatibility antigen HLA-DR1. Nature. 1993, 364: 33-39. 10.1038/364033a0.

    Article  PubMed  CAS  Google Scholar 

  41. Stern LJ, Wiley DC: Antigenic peptide binding by class I and class II histocompatibility proteins. Structure. 1994, 2: 245-251. 10.1016/S0969-2126(00)00026-5.

    Article  PubMed  CAS  Google Scholar 

  42. Yoneda M: The local population division, conservation and management for the Asiatic black bear. Landscape Stud. 2001, 64: 314-317. in Japanese

    Google Scholar 

  43. Sawyer S: Statistical tests for detecting gene conversion. Mol Biol Evol. 1989, 6: 526-538.

    PubMed  CAS  Google Scholar 

  44. Satta Y: Balancing selection at HLA loci. The Proceedings of the 17th Taniguchi Symposium. Edited by: Takahata N. 1992, Science Society Press, Tokyo: Japan, 111-131.

    Google Scholar 

  45. Martin DP, Lemey P, Lott M, Moulton V, Posada D, Lefeuvre P: RDP3: a flexible and fast computer program for analyzing recombination. Bioinformatics (Oxford, England). 2010, 26: 2462-2463. 10.1093/bioinformatics/btq467.

    Article  CAS  Google Scholar 

  46. Saitou N, Nei M: The neighbor-joining method: a new method for reconstructing phylogenetic trees. Mol Biol Evol. 1987, 4: 406-425.

    PubMed  CAS  Google Scholar 

  47. Zhang J, Rosenberg HF, Nei M: Positive Darwinian selection after gene duplication in primate ribonuclease genes. Proc Nat Acad Sci USA. 1998, 95: 3708-3713. 10.1073/pnas.95.7.3708.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  48. Yang Z, Nielsen R, Goldman N, Pedersen AM: Codon-substitution models for heterogeneous selection pressure at amino acid sites. Genetics. 2000, 155: 431-449.

    PubMed  CAS  PubMed Central  Google Scholar 

  49. Yang Z, Wong WSW, Nielsen R: Bayes empirical Bayes inference of amino acid sites under positive selection. Mol Biol Evol. 2005, 22: 1107-1118. 10.1093/molbev/msi097.

    Article  PubMed  CAS  Google Scholar 

  50. Yang Z: PAML 4: phylogenetic analysis by maximum likelihood. Mol Biol Evol. 2007, 24: 1586-1591. 10.1093/molbev/msm088.

    Article  PubMed  CAS  Google Scholar 

  51. Wilson DJ, McVean G: Estimating diversifying selection and functional constraint in the presence of recombination. Genetics. 2006, 172: 1411-1425.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  52. Li W: Molecular Evolution. 1997, Sunderland, Massachusetts: Sinauer Associates

    Google Scholar 

  53. Krause J, Unger T, Noçon A, Malaspinas AS, Kolokotronis SO, Stiller M, Soibelzon L, Spriggs H, Dear PH, Briggs AW, Bray SC, O’Brien SJ, Rabeder G, Matheus P, Cooper A, Slatkin M, Pääbo S, Hofreiter M: Mitochondrial genomes reveal an explosive radiation of extinct and extant bears near the Miocene-Pliocene boundary. BMC Evol Biol. 2008, 8: 220-10.1186/1471-2148-8-220.

    Article  PubMed  PubMed Central  Google Scholar 

  54. Bon C, Caudy N, De Dieuleveult M, Fosse P, Philippe M, Maksud F, Beraud-Colomb É, Bouzaid E, Kefi R, Laugier C, Rousseau B, Casane D, Van Der Plicht J, Elalouf JM: Deciphering the complete mitochondrial genome and phylogeny of the extinct cave bear in the Paleolithic painted cave of Chauvet. Proc Nat Acad Sci USA. 2008, 105: 17447-17452. 10.1073/pnas.0806143105.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  55. Dobson M, Kawamura Y: Origin of the Japanese land mammal fauna: allocation of extant species to historically-based categories. Quat Res. 1998, 37: 385-395. 10.4116/jaqua.37.385.

    Article  Google Scholar 

  56. IUCN Red List of Threatened Species; Version. 2012, [www.iucnredlist.org], .8

  57. Frankham R, Ballou JD, Briscoe DA: Introduction to Conservation Genetics. 2002, Cambridge: Cambridge University Press

    Book  Google Scholar 

  58. Takahata N, Satta Y, Klein J: Polymorphism and balancing selection at major histocompatibility complex loci. Genetics. 1992, 130: 925-938.

    PubMed  CAS  PubMed Central  Google Scholar 

  59. Del Guercio MF, Sidney J, Hermanson G, Perez C, Grey HM, Kubo RT, Sette A: Binding of a peptide antigen to multiple HLA alleles allows definition of an A2-like supertype. J Immunol. 1995, 154: 685-693.

    PubMed  CAS  Google Scholar 

  60. Sidney J, Del Guercio MF, Southwood S, Engelhard VH, Appella E, Rammensee HG, Falk K, Rötzschke O, Takiguchi M, Kubo RT: Several HLA alleles share overlapping peptide specificities. J Immunol. 1995, 154: 247-259.

    PubMed  CAS  Google Scholar 

  61. Greenbaum J, Sidney J, Chung J, Brander C, Peters B, Sette A: Functional classification of class II human leukocyte antigen (HLA) molecules reveals seven different supertypes and a surprising degree of repertoire sharing across supertypes. Immunogenetics. 2011, 63: 325-335. 10.1007/s00251-011-0513-0.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  62. Hughes AL, Nei M: Pattern of nucleotide substitution at major histocompatibility complex class I loci reveals overdominant selection. Nature. 1988, 335: 167-170. 10.1038/335167a0.

    Article  PubMed  CAS  Google Scholar 

  63. Hughes AL, Nei M: Nucleotide substitution at major histocompatibility complex class II loci: evidence for overdominant selection. Proc Nat Acad Sci USA. 1989, 86: 958-962. 10.1073/pnas.86.3.958.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  64. Takahata N, Nei M: Allelic genealogy under overdominant and frequency-dependent selection and polymorphism of major histocompatibility complex loci. Genetics. 1990, 124: 967-978.

    PubMed  CAS  PubMed Central  Google Scholar 

  65. Takahata N: A simple genealogical structure of strongly balanced allelic lines and trans-species evolution of polymorphism. Proc Nat Acad Sci USA. 1990, 87: 2419-2423. 10.1073/pnas.87.7.2419.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  66. Satta Y: Balancing selection at HLA loci. Mechanisms of molecular evolution. Edited by: Takahata N, Clark A. 1993, Sunderland, Massachusetts: Sinauer Associates, 129-149.

    Google Scholar 

  67. Sutton JT, Nakagawa S, Robertson BC, Jamieson IG: Disentangling the roles of natural selection and genetic drift in shaping variation at MHC immunity genes. Mol Ecol. 2011, 20: 4408-20. 10.1111/j.1365-294X.2011.05292.x.

    Article  PubMed  Google Scholar 

  68. Ejsmond MJ, Babik W, Radwan J: MHC allele frequency distributions under parasite-driven selection: A simulation model. BMC Evol Biol. 2010, 10: 332-10.1186/1471-2148-10-332.

    Article  PubMed  PubMed Central  Google Scholar 

  69. Edwards SV, Grahn M, Potts WK: Dynamics of Mhc evolution in birds and crocodilians: amplification of class II genes with degenerate primers. Mol Ecol. 1995, 4: 719-729. 10.1111/j.1365-294X.1995.tb00272.x.

    Article  PubMed  CAS  Google Scholar 

  70. Alcaide M, Edwards SV, Negro JJ, Serrano D, Tella JL: Extensive polymorphism and geographical variation at a positively selected MHC class II B gene of the lesser kestrel (Falco naumanni). Mol Ecol. 2008, 17: 2652-2665. 10.1111/j.1365-294X.2008.03791.x.

    Article  PubMed  CAS  Google Scholar 

  71. Murray BW, Malik S, White BN: Sequence variation at the major histocompatibility complex locus DQ beta in beluga whales (Delphinapterus leucas). Mol Biol Evol. 1995, 12: 582-593.

    PubMed  CAS  Google Scholar 

  72. Tamura K, Peterson D, Peterson N, Stecher G, Nei M, Kumar S: MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol. 2011, 28: 2731-2739. 10.1093/molbev/msr121.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  73. Jones D, Taylor W, Thornton J: The rapid generation of mutation data matrices from protein sequences. Comput Appl Biosci. 1992, 8: 275-282.

    PubMed  CAS  Google Scholar 

  74. Guindon S, Dufayard J-F, Lefort V, Anisimova M, Hordijk W, Gascuel O: New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Syst Biol. 2010, 59: 307-321. 10.1093/sysbio/syq010.

    Article  PubMed  CAS  Google Scholar 

  75. Bryant D, Moulton V: Neighbor-net: an agglomerative method for the construction of phylogenetic networks. Mol Biol Evol. 2004, 21: 255-265.

    Article  PubMed  CAS  Google Scholar 

  76. Huson DH, Bryant D: Application of phylogenetic networks in evolutionary studies. Mol Biol Evol. 2006, 23: 254-267.

    Article  PubMed  CAS  Google Scholar 

  77. Ronquist F, Huelsenbeck JP: MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics. 2003, 19: 1572-1574. 10.1093/bioinformatics/btg180.

    Article  PubMed  CAS  Google Scholar 

  78. Rambaut A: Tree Figure Drawing Tool Version 1.3.1. 2009

    Google Scholar 

  79. Excoffier L, Laval G, Schneider S: Arlequin ver. 3.0: An integrated software package for population genetics data analysis. Evol Bioinformatics Online. 2005, 1: 47-50.

    CAS  Google Scholar 

  80. Hudson RR, Boos DD, Kaplan NL: A statistical test for detecting geographic subdivision. Mol Biol Evol. 1992, 9: 138-151.

    PubMed  CAS  Google Scholar 

  81. Librado P, Rozas J: DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009, 25: 1451-1452. 10.1093/bioinformatics/btp187.

    Article  PubMed  CAS  Google Scholar 

  82. Hedrick PW: A standardized genetic differentiation measure. Evolution. 2005, 59: 1633-1638.

    Article  PubMed  CAS  Google Scholar 

  83. Jost L: GST and its relatives do not measure differentiation. Mol Ecol. 2008, 17: 4015-4026. 10.1111/j.1365-294X.2008.03887.x.

    Article  PubMed  Google Scholar 

  84. Crawford NG: SMOGD: Software for the Measurement of Genetic Diversity. Mol Ecol Resour. 2010, 10: 556-557. 10.1111/j.1755-0998.2009.02801.x.

    Article  PubMed  Google Scholar 

  85. Wu TT, Kabat EA: An analysis of the sequences of the variable regions of Bence Jones proteins and myeloma light chains and their implications for antibody complementarity. J Exp Med. 1970, 132: 211-250. 10.1084/jem.132.2.211.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  86. Jukes T, Cantor C: Evolution of protein molecules. Mammalian protein metabolism. Edited by: Munro H. 1969, New York: Academic Press, 21-132.

    Chapter  Google Scholar 

  87. Nielsen R, Yang Z: Likelihood models for detecting positively selected amino acid sites and applications to the HIV-1 envelope gene. Genetics. 1998, 148: 929-936.

    PubMed  CAS  PubMed Central  Google Scholar 

  88. Sommer S: The importance of immune gene variability (MHC) in evolutionary ecology and conservation. Front Zool. 2005, 2: 16-10.1186/1742-9994-2-16.

    Article  PubMed  PubMed Central  Google Scholar 

  89. Satta Y, O’hUigin C, Takahata N, Klein J: Intensity of natural selection at the major histocompatibility complex loci. Proc Nat Acad Sci USA. 1994, 91: 7184-7188. 10.1073/pnas.91.15.7184.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  90. Nishida S, Tachi T, Baba Y, Hayashi K, Kawanami N, Koike H: DNA analysis of the Asiatic black bear in the Tohoku region using hair-trap samples. Report on model project of networking maintenance for habitat of wildlife. 2002, Wildlife Research Center, Tokyo: Japan, (in Japanese)

    Google Scholar 

  91. Marsden CD, Woodroffe R, Mills MGL, McNutt JW, Creel S, Groom R, Emmanuel M, Cleaveland S, Kat P, Rasmussen GS, Ginsberg J, Lines R, André JM, Begg C, Wayne RK, Mable BK: Spatial and temporal patterns of neutral and adaptive genetic variation in the endangered African wild dog (Lycaon pictus). Mol Ecol. 2012, 21: 1379-1393. 10.1111/j.1365-294X.2012.05477.x.

    Article  PubMed  Google Scholar 

  92. Castro-Prieto A, Wachter B, Sommer S: Cheetah paradigm revisited: MHC diversity in the world’s largest free-ranging population. Mol Biol Evol. 2011, 28: 1455-1468. 10.1093/molbev/msq330.

    Article  PubMed  CAS  Google Scholar 

  93. Becker L, Nieberg C, Jahreis K, Peters E: MHC class II variation in the endangered European mink Mustela lutreola (L. 1761)–consequences for species conservation. Immunogenetics. 2009, 61: 281-288. 10.1007/s00251-009-0362-2.

    Article  PubMed  CAS  Google Scholar 

  94. Niranjan SK, Deb SM, Kumar S, Mitra A, Sharma A, Sakaram D, Naskar S, Sharma D, Sharma SR: Allelic diversity at MHC class II DQ loci in buffalo (Bubalus bubalis): evidence for duplication. Vet Immunol Immunopathol. 2010, 138: 206-212. 10.1016/j.vetimm.2010.07.014.

    Article  PubMed  CAS  Google Scholar 

  95. Sena L, Schneider MPC, Brenig BB, Honeycutt RL, Honeycutt D, Womack JE, Skow LC: Polymorphism and gene organization of water buffalo MHC-DQB genes show homology to the BoLA DQB region. Anim Genet. 2011, 42: 378-385. 10.1111/j.1365-2052.2010.02157.x.

    Article  PubMed  CAS  Google Scholar 

  96. Diaz D, Naegeli M, Rodriguez R, Nino-Vasquez JJ, Moreno A, Patarroyo ME, Pluschke G, Daubenberger CA: Sequence and diversity of MHC DQA and DQB genes of the owl monkey Aotus nancymaae. Immunogenetics. 2000, 51: 528-537. 10.1007/s002510000189.

    Article  PubMed  CAS  Google Scholar 

  97. Meyer-Lucht Y, Sommer S: MHC diversity and the association to nematode parasitism in the yellow-necked mouse (Apodemus flavicollis). Mol Ecol. 2005, 14: 2233-2243. 10.1111/j.1365-294X.2005.02557.x.

    Article  PubMed  CAS  Google Scholar 

Download references

Acknowledgements

The authors wish to acknowledge institutes that helped in sample collection, including the Ministry of the Environment of Japan; the Japan Wildlife Research Center; the Aso Cuddly Dominion and Culture; Shikoku Institute of Natural History; the International Affairs Division, Department of Environment and Civic Affairs, Shimane Prefectural Government; Forestry and Forest Products Research Institute; Forestry Preservation Division, Tokushima prefecture; Agriculture, Forestry and Fisheries Department in Tokushima prefecture; Yamagata University; Keio University; and Hakusan Nature Conservation Center for their kind contribution to this work. This work was supported in part by the Center for the Promotion of Integrated Sciences (CPIS) of Sokendai. The authors also thank Toshinao Okayama, Tomoyuki Uchiyama, Naoki Okabe, Bungo Kanazawa, Masahiro Fujita, Seigo Sawada, Hidetoshi Tamate and Reina Uno for their kind contribution in sample collection and Chris Wood and John A. Eimes for the critical checking of the English within the manuscript. We owe special thanks to Hidenori Tachida, Ryuichi Masuda, Makoto Komoda and John A. Eimes for providing valuable comments.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Yoshiki Yasukochi.

Additional information

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.

Electronic supplementary material

12862_2012_2253_MOESM1_ESM.pdf

Additional file 1: The NJ and Bayesian trees of ursine DQB and DRB 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. (PDF 606 KB)

12862_2012_2253_MOESM2_ESM.pdf

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. (PDF 256 KB)

Authors’ original submitted files for images

Rights and permissions

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Yasukochi, Y., Kurosaki, T., Yoneda, M. et al. MHC class II DQB diversity in the Japanese black bear, Ursus thibetanus japonicus. BMC Evol Biol 12, 230 (2012). https://doi.org/10.1186/1471-2148-12-230

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/1471-2148-12-230

Keywords