Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

Transcriptome analyses of immune tissues from three Japanese frogs (genus Rana ) reveals their utility in characterizing major histocompatibility complex class II

BMC Genomics volume 18, Article number: 994 (2017) Cite this article

Abstract

Background

In Japan and East Asia, endemic frogs appear to be tolerant or not susceptible to chytridiomycosis, a deadly amphibian disease caused by the chytrid fungus Batrachochytridium dendrobatidis (Bd). Japanese frogs may have evolved mechanisms of immune resistance to pathogens such as Bd. This study characterizes immune genes expressed in various tissues of healthy Japanese Rana frogs.

Results

We generated transcriptome data sets of skin, spleen and blood from three adult Japanese Ranidae frogs (Japanese brown frog Rana japonica, the montane brown frog Rana ornativentris, and Tago’s brown frog Rana tagoi tagoi) as well as whole body of R. japonica and R. ornativentris tadpoles. From this, we identified tissue- and stage-specific differentially expressed genes; in particular, the spleen was most enriched for immune-related genes. A specific immune gene, major histocompatibility complex class IIB (MHC-IIB), was further characterized due to its role in pathogen recognition. We identified a total of 33 MHC-IIB variants from the three focal species (n = 7 individuals each), which displayed evolutionary signatures related to increased MHC variation, including balancing selection. Our supertyping analyses of MHC-IIB variants from Japanese frogs and previously studied frog species identified potential physiochemical properties of MHC-II that may be important for recognizing and binding chytrid-related antigens.

Conclusions

This is one of the first studies to generate transcriptomic resources for Japanese frogs, and contributes to further understanding the immunogenetic factors associated with resistance to infectious diseases in amphibians such as chytridiomycosis. Notably, MHC-IIB supertyping analyses identified unique functional properties of specific MHC-IIB alleles that may partially contribute to Bd resistance, and such properties provide a springboard for future experimental validation.

Background

Changes in the host-pathogen-environment interactions can lead to differential susceptibility to infectious disease. Anuran amphibians experience marked changes in physiology and environment due to their biphasic life cycle which includes pre-metamorphic tadpoles that inhabit fully aqueous environments, and post-metamorphic adults that inhabit aqueous and/or terrestrial environments. Anurans undergo complete physiological reorganization at metamorphosis, including the immune system. Some ontological differences in immunity between tadpole and adult frogs include turnover of lymphocytes at metamorphosis, differences in antibody responses and expression of major histocompatibility complex (MHC) antigens (reviewed in [1]). Influence of life stage on immunity is supported by evidence that pathogens (e.g. protists, ranaviruses, fungi) differentially affect either tadpole or adult frog stages [2,3,4]. In addition, variation in habitat can impact the interaction of pathogens and host immune response [5], whereby pathogen diversity could differ between the aqueous environment of tadpoles and the terrestrial /aquatic lifestyle of adult frogs.

Chytridiomycosis is a devastating disease in amphibians caused by the chytrid fungus Batrachochytridium dendrobatidis (Bd). Bd infection results in epidermal disruption and pathophysiological changes sometimes leading to mortality [6], and has been attributed to the decline of amphibian populations worldwide [7,8,9,10]. Although variable Bd prevalence in Korea [11] and Japan [12] has been found, there is a lack of Bd related declines reported in East Asia as well as no evidence of Bd susceptibility in endemic East Asian frogs following experimental infection. There is an ongoing theory that Bd is endemic to Asia [11, 13], so East Asian and Japanese frogs may have had a lengthy co-evolution with the pathogenic fungi and acquired immune tolerance. Therefore, study of immune genes from Japanese amphibian species is important for further understanding amphibian-chytridiomycosis dynamics.

Major histocompatibility complex (MHC) are one of the most polymorphic gene families in vertebrates [14]. They code for membrane-bound glycoproteins that recognize, bind and present antigenic peptides to T lymphocytes, and thus are essential for adaptive immunity in jawed vertebrates. There are two major classes of MHC molecules: MHC class I (MHC-I) predominantly recognize and present endogenous antigenic peptides such as from viruses, while MHC class II (MHC-II) detect and present exogenous-derived peptides such as from bacteria and fungi [15]; because of this, studies of MHC genetics in East Asian frogs in the context of Bd fungus have largely focused on characterization of MHC-II genes [16,17,18]. MHC-II proteins are comprised of non-covalently associated alpha (α or MHC-IIA) and beta (β or MHC-IIB) chains, each with two extracellular domains (α1 and α2, and β1 and β2, respectively). The peptide binding region of the β1 chain has the highest variation, and this diversity governs the repertoire of antigenic determinants to which the host individual can respond [19].

MHC diversity is maintained predominantly by pathogen-mediated balancing selection in an evolutionary time scale [20], and several vertebrate studies have found associations between MHC genetic variation and infectious diseases (reviewed in [21]). In the case of amphibian-chytridiomycosis dynamics, MHC-IIB conformation was suggested to be associated with resistance to Bd [18]; this was based chiefly on amino acid properties at P9 binding pockets that were conserved between resistant Korean frog species and individuals of Australian alpine tree frogs (Litoria verreauxii alpina) that survived Bd infection. Other studies in non-Asian frog species have also identified associations between chytridiomycosis with MHC-IIB alleles or supertypes in wild and captive frogs [22, 23]; supertyping is a method for grouping alleles based on amino acid properties. While such studies have identified MHC class II-chytridiomycosis associations, resistance to disease is often elicited by polygenetic responses and this should be considered. Within Japanese frogs, immunogenetic studies have been limited to skin antimicrobial peptides (AMPs) and MHC class I from Ranidae species [24, 25]. With the increasing availability of next-generation sequencing technologies, the application of transcriptome analyses and genome-wide SNPs is now accessible to non-model organisms including frogs [26,27,28].

In Japan, brown frogs of true frog genus Rana are a common frog species distributed throughout the Japanese archipelago. Of these, three species, Japanese brown frog Rana japonica, the montane brown frog Rana ornativentris, and Tago’s brown frog Rana tagoi tagoi, are found on Honshu (Japanese mainland). Although they are distributed in neighboring areas, each species inhabits different habitats: grasslands from lowland to hillsides (R. japonica), lowland to montane area (R. ornativentris), and montane area close to torrents (R. t. tagoi). This can provide an intrinsic opportunity to compare the expression patterns of immune-related genes and validate selective traits that may be conserved across a variety of ecological and evolutionary backgrounds. The aim of this study is to perform transcriptome analyses of immune-related tissues from these three Japanese Rana frog species and compare expression of genes associated with immune function such as MHC and AMPs. We then focused on MHC-IIB, which may be important for chytridiomycosis resistance, for further analyses of expression and genetic characterization.

Results

Transcriptome data set and differential expression between tissues and life stages

We used our Illumina sequence results to assemble clean reads from each of the three Japanese Rana species into 303,238–646,586 transcripts with an average contig size of 561–650 bp (Table 1). Our BLAST search of all assembled transcripts against sequences of the Swissprot, human Ensembl, Protein family (Pfam), Kyoto Encyclopedia of Genes and Genomes (KEGG) Orthology (KO), and Gene Ontology (GO) databases indicated that among all assembled contigs, 16.61% to 19.45% were annotated using BLAST search in at least one database (Table 1).

Table 1 Summary statistics of Illumina sequencing, assembly and annotation of normalized transcriptomes from three Japanese Ranidae species, including AMP transcripts identified
Full size table

Using the transcriptome data sets derived from RNA samples, we extracted 917–32,363 differentially expressed (DE) transcripts from pairwise comparisons within each of the three focal species. This corresponded to 531–7160 DE genes after encoding with Ensembl gene IDs and 8–74 GO-slim (GO level two) categories, with highest number of DE genes and GO categories in adult-tadpole pairwise comparisons (Additional file 1: Table S1). Several level two GO categories were enriched in specific adult tissues across all three frog species studied (Fig. 1a, Additional file 1: Table S2): for example, in all species spleen was enriched in ‘immune system process’, ‘cell death’ and ‘cell proliferation’, while blood was enriched for several terms within the ‘cellular component’ parent category (‘nucleus’, ‘cell’, ‘intracellular’), and skin was enriched for ‘small molecule metabolic process’ and ‘oxidoreductase activity’. In addition, several GO categories were enriched in tadpoles compared to adult tissues in both R. japonica and R. ornativentris (Fig. 1a, Additional file 1: Table S2).These included ‘cytoskeleton organization’, ‘cell-cell signaling’, and ‘cytoskeletal protein binding’.

Fig. 1
figure 1

Transcriptome data was utilized to identify differentially expressed gene ontology (GO) categories between tissue types and life stages. a Summary of all enriched transcripts (proportion > 0.01) distributed into level two GO categories, run independently in each of the three focal frog species. Only GO categories enriched in a specific adult tissue (blood, skin, spleen) across all three species, or enriched in tadpoles compared to adult tissues in both R. japonica and R. ornativentris, are presented (full details in Additional file 1: Table S2). b Venn diagrams of full GO term enriched groups in blood, skin, blood (and tadpole) from the three species. Number of immune-related enriched GO terms are in black (summarized in Additional file 1: Tables S3 - S5), while total number of enriched GO terms belonging to the ‘biological process’ parent category are in grey. Position of term GO:0019886 (antigen processing and presentation of exogenous peptide antigen via MHC class II) is indicated by boxed number. Photo source: Q. Lau

Full size image

Differential expression of immune-related genes, including MHC and antimicrobial peptides (AMPs)

Further examination centering on the biological process parent term ‘immune system process’ confirmed that the spleen contained the most enriched immune-related GO terms in all three species, with 19 to 46 GO terms (Fig. 1b, Additional file 1: Tables S3 to S5). Nine of the GO terms were common in spleen of all three species: inflammatory response (GO:0006954) and its regulation and negative regulation (GO:0050727 and GO:0050728), leukocyte cell-cell adhesion (GO:0007159) and migration (GO:0050900), T cell differentiation (GO:0030217) and homeostasis (GO:0043029), T cell receptor signaling pathway (GO:0050852), and positive regulation of B cell differentiation (GO:0045579). In blood and skin samples of the focal species, two to seven immune-related GO terms were enriched; of these, only regulation of inflammatory response (GO:0050727) in the skin was enriched across all species. Among adult frogs, the MHC class II related GO term (GO:0019886) was enriched in various samples, including spleen and skin of R. ornativentris and R. japonica, and blood of R. t. tagoi (Fig. 1B, Additional file 1: Tables S3 - S5). Similarly, expression of MHC class II transcripts/variants was the highest in the spleen, with the exception of the Rata01 variant (highest expression in blood, Table 2). In addition, MHC class I normalized expression was higher in spleen and blood and lower in tadpoles (Table 2).

Table 2 Normalized expression of MHC class IIB variants as well as select MHC class I and AMP sequences identified from transcriptome data of three Japanese Rana species, represented by TMM (trimmed mean of log expression ratio) values
Full size table

When comparing adult with tadpole samples, there were no overlapping immune-related GO terms enriched in either R. japonica or R. ornativentris tadpoles. In the R. japonica tadpole, two terms related to complement activation were enriched (GO:0030449 and GO:0006957). In the R. ornativentris tadpole, terms related to intracellular transport of virus (GO:0075733), and interestingly, antigen processing and presentation of exogenous peptide antigen via MHC class I and II (GO:0002479 and GO:0019886, respectively) were enriched. An in-depth look into the GO:0019886 term revealed that genes coding for the MHC class II protein were not upregulated in R. ornativentris tadpoles, but rather genes related to intracellular assembly and transport of the MHC class II protein (Additional file 1: Table S6). This was further validated by relatively low expression levels of MHC class II transcripts in tadpoles (Table 2).

While no specific GO terms related to antimicrobial peptides (AMPs) were enriched, we isolated several AMP transcripts from the three Rana species using tblastn search (17–24 AMP genes per species, Table 1). Seven AMP families were identified in all three species- bradykinin, histone, japonicin, kininogen, preprobrevinin, and preprochensinin. From normalized expression values of AMP genes (Table 2, Additional file 1: Table S7), highest expression was predominant in adult skin relative to all other samples including all 17 AMP genes of R. t. tagoi. Meanwhile, some AMP genes had variable expression: for example, histone was not highly expressed in the skin of all species, and some AMP genes had higher expression in tadpoles, like brevinin-1Toa and andersonin–R in R. japonica and R. ornativentris, respectively (Table 2).

MHC class II sequence analyses

From molecular cloning we identified nine, 11, and 13 MHC-IIB sequences from R. japonica, R. ornativentris and R. t. tagoi, respectively (n = 7 individuals per species, accession numbers MF555153-MF555185). Amplification of multiple loci was only evident in R. japonica, where two to three sequences were identified per individual (Table 3). Seven of the nine R. japonica MHC-IIB variants were shared within the population studied. For the other two focal species, only Raor01, Raor02, and Rata01 variants were shared between individuals. However, identical nucleotide sequences were shared between a R. ornativentris and a R. tagoi tagoi individual (Raor09 and Rata04, Fig. 2), and these two variants were identical at the amino acid level with Raor01 and Raor02 found in other individuals. Interestingly, we also found that the Rata01 variant had high exon 3 sequence similarity to that of the American bullfrog (L. catesbeianus, Fig. 2); this divergent variant clusters phylogenetically with L. catesbeianus at β2 domain (exon 3) but not β1 domain (exon 2) sequences (Fig. 3, Additional file 1: figure S1). There was limited species-specific phylogenetic clustering in either β1 or β2 domain; some variants from different species clustered with strong bootstrap support, for example Raor07 with Rata05/11, and Raor08 with Rata02 (Fig. 3). The number of segregating sites as well as sequence divergence within individuals was highest in R. t. tagoi, but following exclusion of the divergent Rata01 variant that was found in a single individual, the parameters were comparable between the three focal species (Table 3).

Table 3 Summary of MHC-IIB variants, genetic divergence and codon-based Z tests for selection in R. japonica, R. ornativentris, and R. t. tagoi
Full size table
Fig. 2
figure 2

Amino acid alignment of selected MHC-IIB variants from three Japanese Rana species (R. japonica- Raja, R. ornativentris- Raor, and R. t. tagoi- Rata) spanning from exons 1 to 4. All β1 domain (exon 2) sequences from these species were allocated to one of three supertypes (ST-D, ST-E and ST-F), indicated on the right. Positive selected codon sites detected using omegaMap are shaded in grey. Amphibian peptide binding residues [18] are indicated by coloured boxes that represent pocket residues (P4, P6,P7 or P9). For supertyping analyses, the first three pocket residues were excluded due to absent sequence information in some other species (position 9, 11, and 13). Sequences were also included for human-HosaDQB (Genbank accession M33907.1), and transcriptome sequences obtained from R.pirica- Rapi and L. catesbeianus-Lica. **The Raor-09 variant is identical at the amino acid level with Rata-04, Raor-01 and Raor-02 variants

Full size image
Fig. 3
figure 3

Phylogenetic relationships between MHC-IIB variants (amino acid sequences) identified in R. japonica (Raja), R. ornativentris (Raor), and R. t. tagoi (Rata) and other amphibians using neighbour-joining method and 1000 bootstrap replicates. Phylogenies were constructed based from (a) entire β1 domain encoded by exon 2, and (b) entire β2 domain encoded by exon 3. Nodes with bootstrap support >70% are indicated. Accession numbers from other species are indicated on labels, and ‘tr’ indicates sequences obtained from transcriptomic data of R. pirica (Mori, T., unpublished) and L. catesbeianus (DRA accession number SRP051787). The unique variant Rata-01, which clusters differently in each exon, is emphasised by a box. Supertypes allocated to β1 domain variants are indicated by coloured lines; other species used for supertype analyses are not included due to incomplete sequences

Full size image

Selection tests over the entire β1 or β2 domains (exons 2 or 3, respectively) indicated that only the β1 domain of R. japonica is under balancing selection (Z-value = 3.641, p < 0.0001), while the β2 domain of all three focal species is under purifying selection (Table 3). Nevertheless, we identified 22 to 29 codon sites with an excess of non-synonymous substitutions by omegaMap; most sites contain multiple amino acid replacements, which is supportive of balancing selection. The majority of the ‘positive selected sites’ (18–21 sites) were located within exon 2 which encodes the functionally important β1 domain (Fig. 2). Of the 16 important peptide-binding residues in amphibians defined by Bataille et al. [18], 11 sites were also identified as under positive selection in all three Japanese Rana species studied. Recombination breakpoints were identified using GARD with significant topological incongruence (p < 0.05) in R. ornativentris (breakpoint at nucleotide positions 272 and 352) and R. t. tagoi (breakpoint at position 344).

Supertyping analyses

MHC-IIB alleles and variants from the three Japanese Rana species along with alleles from resistant and susceptible frogs of other studies were assigned to one of six supertypes (ST), ST-A to ST-F, by high membership probability (Fig. 4, Additional file 1: figure S2). Within each of the Japanese Rana species, MHC-IIB variants were allocated to at least two supertypes: R. japonica in ST-D and -F, R. ornativentris to ST-E and -F, and both R. t. tagoi and R. pirica assigned to ST-D, −E and -F (Additional file 1: Table S8). Variants from resistant South Korean frogs were allocated to ST-C, −D, and -F, while those from resistant carrier Lithobates catesbeianus were confined to ST-F. Three supertypes (ST-A to ST–C) from this study were similar to four supertypes (ST1 to ST4) identified in Lithobates yavapaiensis by Savage and Zamudio [23]: ST-A equivalent to ST1, ST-B equivalent to ST2 and ST3 which could not be differentiated in our dataset, and ST-C equivalent to ST4. ST-A and -B were specific to L. yavapaiensis, while ST-C included two variants from B. orientalis (South Korea) and ten variants from Litoria verreauxii alpina (Australia). The remaining 19 variants from L. verreauxii alpina were allocated to ST-D, and include the eight alleles found in survivors following Bd infection [18].

Fig. 4
figure 4

a Supertype scatterplot of MHC-IIB alleles or variants (dots) from the three Japanese Rana species and other species studied allocated to one of six supertypes (ellipses, ST-A to ST-F); the bottom-left graph represents the cumulative variance retained by 20 principal components, and the top-left graph represents eigenvalues retained for the discriminant analysis. b The membership probability of each allele or variant from the nine species used to allocate supertypes (alleles from L. verreauxii alpina considered associated with Bd resistance by Bataille et al. [18] are indicated by ‘R’)

Full size image

Discussion

This study is one of the first transcriptome reports from Japanese frogs and includes three species considered resistant to the deadly worldwide fungal disease caused by Batrachochytridium dendrobatidis. We have produced a transcriptome data set focusing on immune genes of healthy frogs, which provides a beneficial framework for future studies. Here, we discuss some applications of this transcriptomic resource and then present one example of such application: an in-depth examination of MHC-IIB, an immune gene potentially important for chytrid fungus resistance.

A comparison of differentially expressed immune genes between adult tissues indicated that immune genes are most highly expressed in the spleen, the major lymphatic tissue in amphibians [29]. Spleen was also found to be more enriched for immune function relative to skin in a transcriptome study of Lithobates yavapaiensis [28]. The skin of all three focal species was enriched for the GO term for regulation of inflammatory response. This supports the idea that amphibian skin may be a critical first line of defense against pathogens in the environment. Amphibian skin is a rich source of anti-microbial peptides (or AMPs), which have been shown to have an important role in innate defense against pathogens such as Bd [30]. Each species excretes an unique array of AMPs, and previous studies in all the focal Japanese frogs have isolated a number of AMPs from skin including japonicin-1 and -2, brevinin-1 and -2, palustrin-2, ranatuerin-2, and temporin families [24, 31, 32]. The transcriptomes we have generated in this study have provided an initial overview of the full AMP repertoire of the three species, including several previously described AMP families, and confirmed high expression predominantly in the skin. Our transcriptome data set provides a robust framework for future in-depth analyses of AMP diversity as well as comparative genomics between amphibian species that differ in Bd susceptibility.

Another focus of this study is the characterization of the MHC-IIB gene, due to its potential role in recognizing and binding of Bd-related peptides. This is supported by recent chytridiomycosis-MHCII association studies [18, 22, 23]. We characterized MHC-IIB in the three focal Japanese Rana species to explore whether selection has occurred at this gene. In addition, we have compared MHC-IIB expression levels within different adult tissues and with tadpoles. A comparison within adults showed that tissues with enriched MHC class II-related GO terms, which include genes related to peptide binding and MHC class II assembly and transport, also had higher expression levels of MHC-IIB genes. Unlike MHC class I, which is expressed in all nucleated somatic cells, class II is expressed predominantly on antigen-presenting cells such as dendrocytes, macrophages and some lymphocytes. MHC-IIB appears to be most highly expressed in the spleen of the Japanese Rana. This aligns with Savage et al. [28] where enrichment of specific GO terms in L. yavapaiensis spleen transcriptome suggested a role for MHC-based immunity. However, the limited number of samples here prevents further conclusions (n = 3). Indeed, expression of specific MHC-IIB loci could vary across tissues; for example, the Rata01 variant has highest expression in blood of the single R. t. tagoi transcriptome (Table 2), and this variant also is of interest due to its seemingly divergent phylogenetic history (Fig. 3).

Surprisingly, GO analyses suggested that MHCII-related terms were upregulated in tadpoles of R. ornativentris (n = 2); this contrasts with transcript expression levels in R. japonica and R. ornativentris where expression was lowest in tadpoles (Table 2). This was resolved following further examination of the enriched GO terms which indicated that differentially expressed genes were related to transport and assembly of MHC-II proteins rather than the actual molecule. One possibility is that during pre-metamorphosis, MHC class II is immature and the MHC class II-making machinery is up-regulated during the tadpole life stage. Studies in tadpoles of the unrelated Xenopus frogs suggest a different pattern whereby MHC class II is ubiquitously expressed during pre-metamorphosis and MHC class I is low until post-metamorphosis [33, 34]. Class I transcripts in the tadpole of R. ornativentris were also lower than adults samples, but limited sample size in this study prevents any conclusions for MHC expression during development. To elucidate whether there are true developmental differences in the Japanese Rana relative to Xenopus, in the future we can investigate MHC class I and II expression levels throughout the tadpole life cycle until metamorphosis. We can utilize the transcriptome resources established in this study to further examine genes important for innate immunity (e.g. complement system, antimicrobial peptides) and adaptive immunity (e.g. MHC class I and II, and genes related to MHC assembly and transport) for quantitative studies of expression level changes and investigation of developmental immunogenetics.

Using primers developed from our transcriptome data, we characterized MHC-IIB from multiple individuals using molecular cloning. Subsequently, we analyzed and identified evidence of the major evolutionary mechanisms that have driven MHC-IIB diversity in the three focal species: balancing selection (e.g. trans-species polymorphism, detection of codon sites with excess of non-synonymous substitutions), recombination and gene duplication. In the three species, signatures of balancing selection were most prominent, especially at the functionally important β1 domain, with identification of many positively selected codon sites in all three species, as well as global excess of non-synonymous substitutions in R. japonica. There was also evidence of trans-species polymorphism, with variants across species phylogenetically clustered with high bootstrap support. In particular there were identical MHC-IIB sequences (Rata04 and Raor09) shared between a R. t. tagoi and a R. ornativentris individual, despite these two species being more evolutionarily divergent relative to R. japonica [35]. Of the other potential mechanisms driving MHC diversity, gene duplication was detected in R. japonica only, while recombination was detected in the other two focal species. The limited evidence of gene duplication and recombination at MHC-IIB in these species is in accordance with the variable number of alleles amplified (one to four) and lack of recombination detection in other ranid frogs [36]. The initial characterization of MHC-IIB in Japanese Rana frogs here can drive future diversity studies of populations across Japan.

This is among the first studies to compare MHC-IIB supertypes across a range of Bd-resistant and Bd-susceptible frogs. While direct experimental evidence of the role of MHC in Bd resistance in Japanese ranids is lacking, we can still speculate that immune genes, including MHC-IIB, contribute to the lack of Bd-related demographic declines in Japan and East Asia. Through supertyping analyses, we obtained an overview of the potential physiochemical properties that may contribute to Bd resistance. Three supertypes (ST-D, ST-E and ST-F) were found in ‘resistant’ frogs from Japan, and ST-D also includes alleles with the advantageous P9 binding pocket conformation (i.e., aromatic β37, acidic Aspβ57, Proβ56 and hydrophobic β60 residues) from resistant frogs of Australia and Korea [18]. Although MHC sequences from Japanese ranids and other Bd-resistant frogs [18] share similar P9 pocket conformations, our observation of MHC sequences (ST4, referred as ST-C in this study) from Bd-resistant North Amrerican Lithobates yavapaiensis [23] clustering most closely with MHCII-B alleles from Australian L. verreauxii alpina frogs that did not survive Bd-infection [18] suggests that protective effects of MHC supertypes may vary across species. This may indicate that distinct Bd-resistant MHC conformations have independently evolved in different populations or continents (America versus Asia-Pacific). Such differences are plausible considering the long evolutionary distance between Japanese and American ranids [37]. In addition, the complex interaction between host-pathogen-environment could drive variation in the protective effect of MHC-IIB, and background differences in MHC structure from neighboring non-peptide binding residues could cause differences in peptide binding.

Future investigations should focus on the impact of amino acid sequence on MHC peptide binding pocket structure [38] and how this influences MHC-Bd peptide binding affinity, using computational [39] or in vitro [40] functional binding assays. Another approach is to investigate the impact of MHC-IIB supertypes on Bd resistance with laboratory challenge studies. In addition, further characterization of other non-MHC genes is important since resistance to diseases like Bd are likely to be polygenic.

Conclusions

Here we sequenced the transcriptome of three Japanese Rana frog species and investigated potential immunogenetic factors contributing to these species’ resistance to the deadly worldwide disease chytridiomycosis. We characterized MHC-IIB for the first time in three species of Japanese ranids, and also explored its expression and molecular evolution. In addition, we expanded on previous studies of functional interpretation of MHC-IIB diversity, by utilizing supertyping analyses to explore the physiochemical properties of MHC-IIB that may be important for recognizing and binding Bd-related antigens and initiating subsequent resistance. By characterizing MHC-IIB in Japanese frogs and exploring its expression and molecular evolution, we have contributed to further understanding complex host-chytridiomycosis dynamics.

Methods

Tissue collection and nucleic acid isolation, library construction and sequencing

We collected tissues from three common Ranidae frog species from Japan: the Japanese brown frog (Rana japonica), the montane brown frog (Rana ornativentris), and Tago’s brown frog (Rana tagoi tagoi). To ensure sufficient sequence coverage, species-specific transcriptome data for the three species were each represented by a single adult individual, with additional tadpole samples for two species. A total of 12 samples were used for transcriptome analysis (Table 1, Additional file 1: Table S7): blood, skin and spleen from a single male adult individual of each species (total nine samples), and three tadpole samples from R. japonica and R. ornativentris. Adult frogs originated from Hiroshima prefecture, Japan, and previously described in [25] (Additional file 1: Table S9); all animals were housed in laboratory conditions for minimum five weeks before euthanasia and exhibited no clinical signs of any disease, thus considered ‘healthy’. Animals were euthanized through immersion in tricaine methanesulfonate (MS222, 0.5–3 g/L water), and samples were collected from the three individuals for immediate RNA extraction. We also included three tadpole samples, euthanized in the same manner, at different larval stages [41]: (i) stage 29 (s29) R. japonica whole body, (ii) s29 R. ornativentris whole body (both from Hiroshima prefecture, Japan, and stored in RNAlater solution (Applied Biosystems, Carlsbad, CA, USA) at -20 °C prior to RNA extraction), and (iii) stage 24 (s24) R. ornativentris skin (three individuals combined) originating from Yokohama Nature Observation Forest park (Kanagawa prefecture, Japan). We used ISOGEN (Nippon Gene, Tokyo, Japan) to extract total RNA following manufacturer instructions from the 12 samples. Sequencing libraries were created from 1 to 3 μg RNA of each sample using NEBNext® Poly(A) mRNA Magnetic Isolation Module and NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (New England Bio Labs, MA, USA) and used to perform short DNA sequencing (paired-end 100-bp) with a cDNA Illumina Hiseq2000 sequencing system (12 samples in one lane).

For the MHC-IIB genetic characterization study, we euthanized and collected spleen samples from an additional six adult individuals per species and extracted RNA using the same methods, giving a total of n = 7 spleen samples per species. PrimeScript™ RT reagent kit (Takara Bio Inc., Otsu, Japan) was used to synthesize first-strand complementary DNA (cDNA). These cDNA samples were previously used to study MHC class I genes [25].

De novo transcriptome assembly and gene annotation

From the short-read data generated by Illumina HiSeq2000 sequencing, adaptor sequences and low-quality reads were removed, and we then performed de novo short-read transcriptome assembly using default settings in Trinity version 2.2.0 [42]. Three independent reference transcript data sets were created for each of the three focal species, combining all samples from the same species: R. japonica (blood, skin, spleen, s29 tadpole body), R. ornativentris (blood, skin, spleen, s24 tadpole skin, s29 tadpole body), and R. t. tagoi (blood, spleen, skin).

To obtain open reading frames and amino acid sequences, assembly results were further processed using the Transdecoder program within the Trinotate platform [43]. Assembled transcripts were annotated with NCBI-BLAST-2.2.30 with E-values ≤1e-5 against five databases: the Swissprot protein database (http://www.expasy.ch/sprot), human amino acid sequence dataset GRCh38.p5 from Ensembl (http://www.ensembl.org/), the Pfam database (http://pfam.xfam.org), the KO database (http://www.genome.jp/kegg/ko.html), and the GO database (http://www.geneontology.org/). Top hits from BLAST search were regarded as annotations and compiled into annotation reports. We also conducted preliminary identification of antimicrobial peptides (AMPs), using NCBI-BLAST-2.2.30 with E-values ≤1e-5 for tblastn search of 1875 publically available anuran defense peptides collected from the Uniprot database (http://www.uniprot.org/) against our transcriptome data sets.

Analysis of differential expression (DE)

Differential expression analysis was conducted between samples within each of the species-specific reference transcripts, with a goal to extract (a) all genes and (b) immune-related genes that have tissue- or stage-specific expression. Reads were aligned to their respective reference transcript using RSEM v 1.3.0 [44] to estimate transcript abundance. EdgeR (v3.4), which permits analysis in the absence of biological replicates, was used to extract differentially expressed (DE) transcripts with dispersion parameter of 0.1. Pairwise sample comparisons were then conducted within each of the three species using false discovery rate (FDR) cut-off of 0.001, and transcripts at least 4-fold differentially expressed were extracted. All DE transcripts collected in each of the 12 samples were encoded with functional information from the annotation reports, including Ensembl ID terms. We then performed gene ontology (GO) enrichment analyses using GOseq [45] on significantly DE genes to get a general overview of tissue- or stage-specific expression. From pairwise comparisons, enriched DE genes were collated with GO-slim categories which represent level two GO parent categories. For immune-related genes, collated DE genes with Ensembl gene IDs were submitted to GO analysis implemented in DAVID [46], and enriched non-parental GO terms under the parent GO term ‘immune system process’ or related to interaction with microorganisms were collected using threshold of gene count over five and p < 0.05.

Identification and characterization of MHC-IIB sequences

We scanned the NCBI-BLAST search results for all MHC-IIB sequences independently in the three species. All MHC-IIB transcript sequences that were full coding sequences, spanning from exons one to four, were retained and confirmed by additional search using NCBI-BLAST-2.4.0 against the nucleotide collection (nr/nt) database. To compare MHC-IIB expression between tissue types and between adult and tadpole frogs, we extracted trimmed mean log expression ratio, or TMM-normalized values [47], for MHC-IIB transcripts in each of the 12 samples generated by RSEM v 1.3.0. We also collected TMM-normalized values for MHC class I transcripts and candidate AMP genes identified from prior tblastn search.

From the alignment of MHC-IIB transcript sequences from R. japonica (2 sequences), R. ornativentris (1 sequence), and R. t. tagoi (2 sequences), we designed degenerate primers using Primer3 [48] to amplify partial coding sequence of MHC-IIB. The forward 5’-GMAKYAYHWGCAGMAASATG-3′ and reverse 5’-CACWCCRGCAAYRATAARYA-3′ primers are located in exons 1 and 4, respectively, and amplify all of exons 2 (β1 domain) and 3 (β2 domain). Polymerase chain reaction (PCR) amplification, cloning reactions and sequencing were conducted on spleen cDNA samples (n = 7 per species) following Lau et al. [25], except using a lower PCR annealing temperature of 48 °C.

MHC-IIB sequence analyses, phylogenetics, and selection tests

Criteria for inclusion of MHC-IIB sequences followed Lau et al. [25] to avoid PCR and cloning artifacts, briefly: (1) amplified from more than one clone, and (2) differed to other sequences by more than two nucleotides. Exceptions were allowed for sequences that differed by one or two non-synonymous substitutions, as they were confirmed by repeat PCR and cloning reactions. For phylogenetic analyses, partial MHC-IIB coding sequences from the three focal species were aligned in MEGA7 [49] with that of other amphibians: Rana pirica (Mori Tsukasa, unpublished transcriptome data), L. catesbeianus [50] (transcriptome data), Rhinella marina [51], Xenopus tropicalis [52], X. laevis [53], Ambystoma mexicanum [54], A. tigrinum [55], Quasipaa spinosa [16], and Andrias davidianus [56]. Neighbor joining trees (p-distance, 1000 bootstrap replicates) were constructed from amino acid alignments for the coding sequence of exon 2 (β1 domain) and exon 3 (β2 domain). Maximum likelihood trees (WAG + G model, 1000 bootstrap replicates) were also constructed in MEGA7 to confirm phylogenetic relationships.

Molecular evidence of balancing selection can be inferred by an excess of nonsynonymous (dN) to synonymous (dS) substitutions, and we tested dN – dS to assess if historical balancing selection has shaped MHC-IIB sequence diversity in Japanese Rana. Firstly, we used Z-tests for positive (dN > dS), neutral (dN = dS) and purifying (dN < dS) selection in MEGA7 to assess if selection is acting globally in β1 and β2 domain sequences, for each species independently. We also used a codon-based approach, omegaMap version 5.0 [57], to identify individual codons under selection (dN/dS or ω > 1 with a posterior probability of >0.95) even in the presence of recombination. This was conducted for sequence alignments for each of the species independently, and selection parameter (ω) and population recombination rate (ρ) was co-estimated and allowed to vary along the sequence. Two independent runs were conducted following Lau et al. [58] with 1 × 106 Markov chain Monte Carlo iterations, and results were visualized and interpreted using R v 3.3.2 [59]. To test for presence of recombination, we used GARD (Genetic Algorithm Recombination Detection) [60], implemented in the Datamonkey website (http://www.datamonkey.org).

Supertyping analyses

Supertyping analyses were used to characterize the functional diversity of the MHC-IIB variants identified from the focal Japanese Rana frog species compared to that of six other frog species: resistant Bombina orientalis and Bufo gargarizans from South Korea [18]; resistant carrier L. catesbeiana (transcriptome data, DRA accession number SRP051787); resistant R. pirica from Hokkaido Japan (Mori, T., unpublished transcriptome data); and Lithobates yavapaiensis from North America [23] and Litoria verreauxii alpina from Australia [18], both of which are susceptible species with resistant populations or individuals. The supertyping method groups alleles based on similar physiochemical properties at specified amino acid sites. Following Savage and Zamudio [23], we aligned all MHC-IIB amino acid sequences identified from the focal species with that of the other species, and extracted 13 codon positions from the β1 domain (exon 2) identified as peptide-binding region in frogs [18]. Each of the 13 sites were characterized for five physiochemical descriptor variables: z1 (hydrophobicity), z2 (steric bulk), z3 (polarity), z4 and z5 (electronic effects) [61, 62]. Using adegenet 2.0.1 in R v 3.3.2, discriminant analysis of principle components (DAPC) was used to define functional genetic clusters. Using K-means clustering algorithm with Bayesian information criterion (BIC), alleles or variants were classified into clusters that represented distinct MHC supertypes.

Abbreviations

AMPs:

antimicrobial peptides

Bd:

Batrachochytridium dendrobatidis

GO:

gene ontology

MHC:

major histocompatibility complex

MHC-IIB:

major histocompatibility complex class II beta chain

ST:

supertype

References

  1. Rollins-Smith LA. Metamorphosis and the amphibian immune system. Immunol Rev. 1998:221–30.

  2. Chambouvet A, Gower DJ, Jirků M, Yabsley MJ, Davis AK, Leonard G, et al. Cryptic infection of a broad taxonomic and geographic diversity of tadpoles by Perkinsea protists. Proc Natl Acad Sci U S A. 2015;112:E4743–51. doi:10.1073/pnas.1500163112.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Haislip NA, Gray MJ, Hoverman JT, Miller DL. Development and disease: how susceptibility to an emerging pathogen changes through anuran development. PLoS One. 2011;6:e22307. doi:10.1371/journal.pone.0022307.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Blaustein AR, Romansic JM, Scheessele EA, Han BA, Pessier AP, Longcore JE. Interspecific variation in susceptibility of frog tadpoles to the pathogenic fungus Batrachochytrium dendrobatidis. Conserv Biol. 2005;19:1460–8. doi:10.1111/j.1523-1739.2005.00195.x.

    Article  Google Scholar 

  5. James TY, Toledo LF, Rödder D, da Silva Leite D, Belasen AM, Betancourt-Román CM, et al. Disentangling host, pathogen, and environmental determinants of a recently emerged wildlife disease: lessons from the first 15 years of amphibian chytridiomycosis research. Ecol. Evol. 2015;5:4079–4097. doi:10.1002/ece3.1672.

  6. Voyles J, Young S, Berger L, Campbell C, Voyles WF, Dinudom A. Pathogenesis of chytridiomycosis, a cause of catastrophic amphibian declines. Science. 2009;326:582–5.

    Article  CAS  PubMed  Google Scholar 

  7. Longcore J, Longcore J, Pessier A, Halteman WA. Chytridiomycosis widespread in anurans of northeastern United States. J Wildl Manag. 2007;71:435–44. doi:10.2193/2006-345.

    Article  Google Scholar 

  8. Skerratt LF, Berger L, Speare R, Cashins S, McDonald KR, Phillott AD, et al. Spread of chytridiomycosis has caused the rapid global decline and extinction of frogs. EcoHealth. 2007;4:125–34. doi:10.1007/s10393-007-0093-5.

  9. Daszak P, Cunningham AA, Hyatt AD. Infectious disease and amphibian population declines. Divers Distrib. 2003:141–50. doi:10.1046/j.1472-4642.2003.00016.x.

  10. Wake DB, Vredenburg VT. Are we in the midst of the sixth mass extinction? A view from the world of amphibians. Proc Natl Acad Sci U S A. 2008;105:11466–73. doi:10.1073/pnas.0801921105.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Bataille A, Fong JJ, Cha M, Wogan GOU, Baek HJ, Lee H, et al. Genetic evidence for a high diversity and wide distribution of endemic strains of the pathogenic chytrid fungus Batrachochytrium dendrobatidis in wild Asian amphibians. Mol Ecol. 2013;22:4196–209. doi:10.1111/mec.12385.

    Article  CAS  PubMed  Google Scholar 

  12. Goka K, Yokoyama J, Une Y, Kuroki T, Suzuki K, Nakahara M, et al. Amphibian chytridiomycosis in Japan: distribution, haplotypes and possible route of entry into Japan. Mol Ecol. 2009;18:4757–74. doi:10.1111/j.1365-294X.2009.04384.x.

    Article  CAS  PubMed  Google Scholar 

  13. Fisher MC. Endemic and introduced haplotypes of Batrachochytrium dendrobatidis in Japanese amphibians: sink or source? Mol Ecol. 2009;18:4731–3. doi:10.1111/j.1365-294X.2009.04385.x.

  14. Hedrick PW. Evolutionary genetics of the major histocompatibility complex. Am Nat. 1994;143:945–64. doi:10.1086/285643.

    Article  Google Scholar 

  15. Braciale TJ, Morrison L a, Sweetser MT, Sambrook J, Gething MJ, Braciale VL. Antigen presentation pathways to class I and class II MHC-restricted T lymphocytes. Immunol Rev 1987;98:95–114. doi:10.1111/j.1600-065X.1987.tb00521.x.

  16. Yu X, Zheng R, Zhang J, Shen B, Dong B. Genetic polymorphism of major histocompatibility complex class IIB alleles and pathogen resistance in the giant spiny frog Quasipaa spinosa. Infect. Genet. Evol. 2014;28:175–82. doi:10.1016/j.meegid.2014.09.028.

    Article  CAS  PubMed  Google Scholar 

  17. Li F, Shu Y, Polymorphism WH. Of exon 2 of MHC class II B gene in the Chinese concave-eared torrent frog (Odorrana tormota). Biodivers Sci. 2013;20:184–92. doi:10.3724/SP.J.1003.2012.09211.

    Article  Google Scholar 

  18. Bataille A, Cashins SD, Grogan L, Skerratt LF, Hunter D, McFadden M, et al. Susceptibility of amphibians to chytridiomycosis is associated with MHC class II conformation. Proc R Soc B Biol Sci. 2015;282:–20143127. doi:10.1098/rspb.2014.3127.

  19. Brown JH, Jardetzky TS, Gorga JC, Stern LJ, Urban RG, Strominger JL, et al. Three-dimensional structure of the human class II histocompatibility antigen HLA-DR1. Nature. 1993;364:33–9. doi:10.1038/364033a0.

    Article  CAS  PubMed  Google Scholar 

  20. Spurgin LG, Richardson DS. How pathogens drive genetic diversity: MHC, mechanisms and misunderstandings. Proc R Soc B Biol Sci. 2010;277:979–88. doi:10.1098/rspb.2009.2084.

    Article  CAS  Google Scholar 

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

    Article  PubMed  PubMed Central  Google Scholar 

  22. Savage AE, Zamudio KR. MHC genotypes associate with resistance to a frog-killing fungus. Proc Natl Acad Sci U S A. 2011;108:16705–10. doi:10.1073/pnas.1106893108.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Savage AE, Zamudio KR. Adaptive tolerance to a pathogenic fungus drives major histocompatibility complex evolution in natural amphibian populations. Proc Biol Sci. 2016;283:20153115. doi:10.1098/rspb.2015.3115.

    Article  PubMed  PubMed Central  Google Scholar 

  24. Ohnuma A, Conlon JM, Iwamuro S. Differential expression of genes encoding preprobrevinin-2, prepropalustrin-2, and preproranatuerin-2 in developing larvae and adult tissues of the mountain brown frog Rana ornativentris. Comp Biochem Physiol - C Toxicol Pharmacol. 2010;151:122–30. doi:10.1016/j.cbpc.2009.09.004.

    Article  PubMed  Google Scholar 

  25. Lau Q, Igawa T, Komaki S, Satta Y. Characterisation of major histocompatibility complex class I genes in Japanese Ranidae frogs. Immunogenetics. 2016;68:797–806. doi:10.1007/s00251-016-0934-x.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Price SJ, Garner TWJ, Balloux F, Ruis C, Paszkiewicz KH, Moore K, et al. A de novo assembly of the common frog (Rana temporaria) transcriptome and comparison of transcription following exposure to Ranavirus and Batrachochytrium dendrobatidis. PLoS One. 2015;10:1–23. doi:10.1371/journal.pone.0130500.

    CAS  Google Scholar 

  27. Huang L, Li J, Anboukaria H, Luo Z, Zhao M, Wu H. Comparative transcriptome analyses of seven anurans reveal functions and adaptations of amphibian skin. Sci Rep. 2016;6:24069. doi:10.1038/srep24069.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Savage AE, Kiemnec-Tyburczy KM, Ellison AR, Fleischer RC, Zamudio KR. Conservation and divergence in the frog immunome: pyrosequencing and de novo assembly of immune tissue transcriptomes. Gene. 2014;542:98–108. doi:10.1016/j.gene.2014.03.051.

    Article  CAS  PubMed  Google Scholar 

  29. Du Pasquier L, Schwager J, Flajnik MF. The immune system of Xenopus. Annu Rev Immunol. 1989;7:251–75. doi:10.1146/annurev.iy.07.040189.001343.

    Article  CAS  PubMed  Google Scholar 

  30. Rollins-Smith LA. The role of amphibian antimicrobial peptides in protection of amphibians from pathogens linked to global amphibian declines. Biochim Biophys Acta Biomembr. 2009;1788:1593–9. doi:10.1016/j.bbamem.2009.03.008.

    Article  CAS  Google Scholar 

  31. Conlon JM, Coquet L, Jouenne T, Leprince J, Vaudry H, Iwamuro S. Evidence from the primary structures of dermal antimicrobial peptides that Rana tagoi okiensis and Rana tagoi tagoi (Ranidae) are not conspecific subspecies. Toxicon. 2010;55:430–5. doi:10.1016/j.toxicon.2009.09.010.

    Article  CAS  PubMed  Google Scholar 

  32. Isaacson T, Soto A, Iwamuro S, Knoop FC, Conlon JM. Antimicrobial peptides with atypical structural features from the skin of the Japanese brown frog Rana japonica. Peptides. 2002;23:419–25. doi:10.1016/S0196-9781(01)00634-9.

    Article  CAS  PubMed  Google Scholar 

  33. Du Pasquier L, Flajnik MF. Expression of MHC class II antigens during Xenopus development. Dev Immunol. 1990;1:85–95.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Flajnik MF, Kaufman JF, Hsu E, Manes M, Parisot R, Du Pasquier L. Major Histocompatibility complex-encoded class I molecules are absent in immunologically competent Xenopus before metamorphosis. J Immunol. 1986;137:3891–9.

    CAS  PubMed  Google Scholar 

  35. Yuan ZY, Zhou WW, Chen X, Poyarkov NA, Chen HM, Jang-Liaw NH, et al. Spatiotemporal diversification of the true frogs (genus Rana): a historical framework for a widely studied group of model organisms. Syst Biol. 2016;65:824–42. doi:10.1093/sysbio/syw055.

    Article  PubMed  Google Scholar 

  36. Kiemnec-Tyburczy KM, Richmond JQ, Savage AE, Zamudio KR. Selection, trans-species polymorphism, and locus identification of major histocompatibility complex class IIb alleles of new world ranid frogs. Immunogenetics. 2010;62:741–51. doi:10.1007/s00251-010-0476-6.

    Article  CAS  PubMed  Google Scholar 

  37. Bossuyt F, Brown RM, Hillis DM, Cannatella DC, Milinkovitch MC. Phylogeny and biogeography of a cosmopolitan frog radiation: late cretaceous diversification resulted in continent-scale endemism in the family ranidae. Syst Biol. 2006;55:579–94. doi:10.1080/10635150600812551.

    Article  PubMed  Google Scholar 

  38. Jones EY, Fugger L, Strominger JL, Siebold C. MHC class II proteins and disease: a structural perspective. Nat Rev Immunol. 2006;6:271–82. doi:10.1038/nri1805.

    Article  CAS  PubMed  Google Scholar 

  39. Hoof I, Peters B, Sidney J, Pedersen LE, Sette A, Lund O, et al. NetMHCpan, a method for MHC class i binding prediction beyond humans. Immunogenetics. 2009;61:1–13. doi:10.1007/s00251-008-0341-z.

    Article  CAS  PubMed  Google Scholar 

  40. Salvat R, Moise L, Bailey-Kellogg C, Griswold KE. A high throughput MHC II binding assay for quantitative analysis of peptide epitopes. J Vis Exp. 2014;(85):1–11. doi:10.3791/51308.

  41. Gosner KLA. Simplified table for staging anuran embryos larvae with notes on identification. Herpetologica. 1960;16:183–90. doi:10.2307/3890061.

    Google Scholar 

  42. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Trinity: reconstructing a full-length transcriptome without a genome from RNA-Seq data. Nat Biotechnol. 2013;29:644–52. doi:10.1038/nbt.1883.Trinity.

    Article  Google Scholar 

  43. Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, et al. De novo transcript sequence reconstruction from RNA-seq using the trinity platform for reference generation and analysis. Nat Protoc. 2013;8:1494–512. doi:10.1038/nprot.2013.084.

    Article  CAS  PubMed  Google Scholar 

  44. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323. doi:10.1186/1471-2105-12-323.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11:R14. doi:10.1186/gb-2010-11-2-r14.

    Article  PubMed  PubMed Central  Google Scholar 

  46. Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4:44–57. doi:10.1038/nprot.2008.211.

    Article  CAS  Google Scholar 

  47. Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010;11:R25. doi:10.1186/gb-2010-11-3-r25.

    Article  PubMed  PubMed Central  Google Scholar 

  48. Koressaar T, Remm M. Enhancements and modifications of primer design program Primer3. Bioinformatics. 2007;23:1289–91.

  49. Kumar S, Stecher G, Tamura K. MEGA7: molecular evolutionary genetics analysis version 7.0 for bigger datasets. Mol Biol Evol. 2016;33:1870–4. doi:10.1093/molbev/msw054.

    Article  CAS  PubMed  Google Scholar 

  50. Birol I, Behsaz B, Hammond SA, Kucuk E, Veldhoen N, Helbing CC. De novo transcriptome assemblies of Rana (Lithobates) catesbeiana and Xenopus laevis tadpole livers for comparative genomics without reference genomes. PLoS One. 2015;10:1–18. doi:10.1371/journal.pone.0130720.

    Article  Google Scholar 

  51. Lillie M, Cui J, Shine R, Belov K. Molecular characterization of MHC class II in the Australian invasive cane toad reveals multiple splice variants. Immunogenetics. 2016;68:449–60. doi:10.1007/s00251-016-0919-9.

    Article  CAS  PubMed  Google Scholar 

  52. Klein SL, Strausberg RL, Wagner L, Pontius J, Clifton SW, Richardson P. Genetic and genomic tools for Xenopus research: the NIH Xenopus initiative. Dev Dyn. 2002:384–91. doi:10.1002/dvdy.10174.

  53. Sato K, Flajnik MF, Pasquier D, Katagiri M, Kasahara M, Sato K, et al. Evolution of the MHC: isolation of class II beta-chain cDNA clones from the amphibian Xenopus laevis. J Immunol. 1993;150:2831–43.

    CAS  PubMed  Google Scholar 

  54. Laurens V, Chapusot C, del Rosario OM, Bentrari F, Padros MR, Tournefier A. Axolotl MHC class II beta chain: predominance of one allele and alternative splicing of the beta1 domain. Eur J Immunol. 2001;31:506–15. doi:10.1002/1521-4141(200102)31:2<506::AID-IMMU506>3.0.CO;2-P.

    Article  CAS  PubMed  Google Scholar 

  55. Bos DH, DeWoody JA. Molecular characterization of major histocompatibility complex class II alleles in wild tiger salamanders (Ambystoma tigrinum). Immunogenetics. 2005;57:775–81. doi:10.1007/s00251-005-0038-5.

    Article  CAS  PubMed  Google Scholar 

  56. Zhu R, Chen ZY, Wang J, Di YJ, Liao XY, Gui JF, et al. Extensive diversification of MHC in Chinese giant salamanders Andrias davidianus (Anda-MHC) reveals novel splice variants. Dev Comp Immunol. 2014;42:311–22. doi:10.1016/j.dci.2013.10.001.

    Article  CAS  PubMed  Google Scholar 

  57. Wilson DJ, McVean G. Estimating diversifying selection and functional constraint in the presence of recombination. Genetics. 2006;172:1411–25. doi:10.1534/genetics.105.044917.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Lau Q, Jaratlerdsiri W, Griffith JE, Gongora J, Higgins DP. MHC class II diversity of koala (Phascolarctos cinereus) populations across their range. Heredity (Edinb). 2014;113:1–10. doi:10.1038/hdy.2014.30.

    Article  Google Scholar 

  59. R Core Team. R: A language and environment for statistical computing. Vienna, Austria; 2016. Available from: https://www.r-project.org

  60. Kosakovsky Pond SL, Posada D, Gravenor MB, Woelk CH, Frost SDW. GARD: a genetic algorithm for recombination detection. Bioinformatics. 2006;22:3096–8. doi:10.1093/bioinformatics/btl474.

    Article  PubMed  Google Scholar 

  61. Sandberg M, Eriksson L, Jonsson J, Sjöström M, Wold S. New chemical descriptors relevant for the design of biologically active peptides. A multivariate characterization of 87 amino acids. J Med Chem. 1998;41:2481–91. doi:10.1021/jm9700575.

    Article  CAS  PubMed  Google Scholar 

  62. Didinger C, Eimes JA, Lillie M, Waldman B. Multiple major histocompatibility complex class I genes in asian anurans: ontogeny and phylogeny of expression. Dev Comp Immunol. 2016;70:69–79. doi:10.1016/j.dci.2016.12.003.

    Article  PubMed  Google Scholar 

Download references

Acknowledgements

We thank staff of Hiroshima University Institute for Amphibian Biology for animal husbandry, Yohei Terai for assistance with cDNA library preparation, Mori Tsukasa for contribution of transcriptome data from Rana pirica, and Anna Savage for providing detailed MHC-II supertype information from L. yavapaiensis.

Funding

This work was funded by the Japan Society for the Promotion of Science (JSPS). QL was financially supported by JSPS postdoctoral research fellowship program coordinated by the Australian Academy of Science.

Availability of data and materials

All MHC-IIB sequences characterized in this study are deposited in GenBank under accession numbers MF555153-MF555185 (https://0-www-ncbi-nlm-nih-gov.brum.beds.ac.uk/nuccore/MF555170 until https://0-www-ncbi-nlm-nih-gov.brum.beds.ac.uk/nuccore/MF555185). Transcriptome data are deposited in DDBJ under accession number DRA006004 (http://trace.ddbj.nig.ac.jp/DRASearch/submission?acc=DRA006004). MHC-IIB phylogenetic data including alignments are deposited in TreeBASE study accession no. S21988 (http://purl.org/phylo/treebase/phylows/study/TB2:S21988).

Author information

Authors and Affiliations

  1. Department of Evolutionary Studies of Biosystems, Sokendai, The Graduate University for Advanced Studies, Kamiyamaguchi 1560-35, Hayama, Kanagawa, 240-0193, Japan

    Quintin Lau & Yoko Satta

  2. Amphibian Research Center, Hiroshima University, 1-3-1, Higashi-Hiroshima, Hiroshima, 739-8526, Japan

    Takeshi Igawa

  3. Department of Bioscience, Nagahama Institute of Bio-Science and Technology, Tamura-cho 1266, Nagahama, Shiga, 526-0829, Japan

    Ryuhei Minei

  4. One Health Research Group, College of Public Health, Medical and Veterinary Sciences, James Cook University, Townsville, QLD, 4811, Australia

    Tiffany A. Kosch

Authors
  1. Quintin Lau
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Takeshi Igawa
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Ryuhei Minei
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Tiffany A. Kosch
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Yoko Satta
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

QL, TI and YS conceived the study. QL conducted the experiments, analyzed data and drafted the manuscript. TI arranged samples and assisted in assembly and annotation of transcriptome data. RM wrote custom python scripts and conducted differential expression analyses. TK assisted with MHC sequence and supertype analyses and contributed to drafting the manuscript. YS contributed to all evolutionary analyses and contributed to drafting the manuscript. All authors read and approved the manuscript.

Corresponding author

Correspondence to Quintin Lau.

Ethics declarations

Ethics approval and consent to participate

All samples from adult frogs are identical to those used in a prior study [25] and approved by Hiroshima University Animal Research Committee (approval number G14–2). Depending on the species, adult frogs were collected at different life stages at different locations within Hiroshima prefecture and then raised in captivity at Hiroshima University Institute for Amphibian Biology: in brief, R. japonica originated from multiple egg clusters collected from Etajima, R. ornativentris originated from tadpoles collected from Yoshiwa, and R. tagoi tagoi originated from adults collected from Shobara (full details in Additional file 1: Table S9). S29 tadpole samples were collected from the same location within Hiroshima. The S24 tadpole samples collected from Yokohama were approved by Sokendai Research committee (approval number 46) and with permission from Yokohama Nature Observation Forest Park. All handling of adult and tadpole frogs were conducted under strict animal handling guidelines set by both research committees.

Consent for publication

All authors declare consent to publish.

Competing interests

The authors declare that they have no potential conflict of interest.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional file

Additional file 1:

Supplementary information, including Tables S1–S9, and Figures S1–S2. (DOCX 267 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Lau, Q., Igawa, T., Minei, R. et al. Transcriptome analyses of immune tissues from three Japanese frogs (genus Rana ) reveals their utility in characterizing major histocompatibility complex class II. BMC Genomics 18, 994 (2017). https://0-doi-org.brum.beds.ac.uk/10.1186/s12864-017-4404-0

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://0-doi-org.brum.beds.ac.uk/10.1186/s12864-017-4404-0

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us