Genomic organization and adaptive evolution of IGHC genes in marine mammals.

Immunoglobulins are important elements of the adaptive immune system that bind to an immense variety of microbial antigens to neutralize infectivity and specify effector functions. In the present study, the immunoglobulin heavy chain constant region (IGHC) genes from marine mammals were identified and compared with those of their terrestrial relatives to explore their genomic organization and evolutionary characteristics. The genomic organization of marine mammal IGHC genes was shown to be conservative with other eutherian mammals. Stronger signals of positive selection on IGHC were revealed in terrestrial mammals than that in marine mammals with the branch-site model, displaying different selective pressure, which might suggest their divergent adaptations to contrasted environments.

Introduction

Immunoglobulins (Igs) are glycoprotein molecules that function as antibodies and are only expressed in jawed vertebrates (Flajnik and Kasahara, 2010). Typically, Igs compose two identical heavy (IgH) and two identical light chains (IgL). Both heavy and light chains can be divided into two regions, i.e., the constant region and the variable region, based on the variability of amino acid sequences. In general, the variable region recognizes antigens, whereas the constant region specifies critical effector or signaling functions, such as phagocytosis, binding Fc (fragment crystallizable) receptors and activation of complement (Schroeder and Cavacini et al., 2010). The classification of Igs is based on their heavy chain constant (IGHC) genes (Bengtén et al., 2000). To date, several IgH isotypes, including IgM, IgD, IgG, IgE, and IgA, have been identified in multiple vertebrates, and they are encoded by the IGHM, IGHD, IGHG, IGHE and IGHA genes respectively (Schroeder and Cavacini et al., 2010). Both IgM and IgD have been reported to be the most ancient IgH isotypes among jawed vertebrates, with IgM very conserved but IgD showing a high degree of structural and (presumably) functional diversity and even being lost in some species over evolutionary time (Sun et al., 2011; Ohta and Flajnik, 2006). IgG and IgE are both exclusively present in mammals (Bengtén et al., 2000; Flajnik and Kasahara, 2010), the former being the predominant serum antibody and binding receptors on phagocytic cells and the latter predominantly functioning in immunity against parasites (Maizels, 2005). Additionally, IgA was found in reptiles first, with similar functions in all animals (Flajnik and Kasahara, 2010). Participating in mucosal immunity, IgA coats pathogens to prevent proliferation and defend against local infection (Macpherson et al., 2008).

Extant marine mammals descended from terrestrial ancestors that re-entered in the aquatic environment at different times (Uhen, 2010). There are five main groups of extant marine mammals: Cetacea (dolphins, porpoises and whales), Sirenia (manatees, dugongs), Pinnipedia (walruses, sea lions and seals), polar bear (Ursus maritimus) and sea otter (Enhydra lutris) (Uhen, 2007). Cetacea and Sirenia are the only known completely aquatic mammals, and they appeared to return to the aquatic environment at approximately the same time, i.e. ∼50 Mya (million years ago) (Uhen, 2010; Thewissen et al., 2007). In contrast, semi-aquatic pinnipeds, with a close relationship to Musteloidea, later transferred to water environment at ∼28 Mya (Uhen, 2010; Flynn et al., 2005). Despite the independent evolutionary origins of these marine mammals groups, they have developed a series of specializations for aquatic lifestyles, such as loss of hair, thickened blubber, derivation of echolocation, and so on (Uhen, 2007).

During their transition from terrestrial environment to the sea, marine mammals must have faced dramatic changes in pathogens. Recent studies have shown that cetacean TLR4 and MHC were subjected to adaptive evolution, particularly in the lineages with dramatic habitat transition (e.g. the lineage leading to hippo + whale representing the early stage of aquatic adaptation) or rapid radiation lineages such as extant dolphins (Shen et al., 2012; Zhang et al., 2016). However, the evolution of Ig in the aquatic adaptation of marine mammals remains poorly explored, except for a few reports about the cloning and sequencing of IGHC genes in bottlenose dolphin (Tursiops truncatus) (Mancia et al., 2006, 2007; Lundqvist et al., 2002) and the Florida manatee (Trichechus manatus latirostris) (Breaux et al., 2017) and the identification of IgM, IgG and IgA and lineage-specific IGHG gene duplication in some cetaceans (Mancia et al., 2006), and so on. Therefore, in the present study, we first retrieved IGHC genes from the genomes of ten representative marine mammals. The genomic organization of IGHC cluster was identified and positive selection on the marine mammals was detected to address whether they experienced adaptive evolution in response to aquatic environments.

Methods

Genomic identification of IGHC genes

We first downloaded the IGHC genes of 13 terrestrial mammals, including cow (Bos taurus), pig (Sus scrofa), horse (Equus caballus), dog (Canis lupus familiaris), giant panda (Ailuropoda melanoleuca), black flying fox (Pteropus alecto), little brown bat (Myotis lucifugus), human (Homo sapiens), rhesus monkey (Macaca mulatta), rat (Rattus norvegicus), rabbit (Oryctolagus cuniculus), opossum (Monodelphis domestica) and platypus (Ornithorhynchus anatinus) from the National Center for Biotechnology Information (NCBI), Ensemble databases and the International ImMunoGeneTics information system (IMGT). The accession numbers of these sequences are compiled in Supplemental Table 1. Then we used IGHC sequences of these terrestrial relatives as a query to search for the orthologous genes from 10 marine mammal genomes with TBLASTN and BLASTN approaches, including bottlenose dolphin (Tursiops truncatus), killer whale (Orcinus orca), Yangtze finless porpoise (Neophocuena asiaeorientalis), baiji (Lipotes vexillifer), sperm whale (Physeter catodon), minke whale (Balaena acutorostrata), bowhead whale (Balaena mysticetus), Weddell seal (Leptonychotes weddellii), Pacific walrus (Odobenus rosmarus), and Florida manatee (Trichechus manatus latirostris). In addition, IGHC genes were also scanned from the genomes of three additional terrestrial mammals, i.e., hedgehog (Erinaceus europaeus), European shrew (Sorex araneus), and elephant (Loxodonta africana). The detailed information about the genome assemblies used in this study is given in Table 1 . Finally, we checked whether the newly identified putative IGHC genes were their best hit by a blast search against the non-redundant database from GenBank. The newly identified IGHC genes were further separated into intact genes, partial genes and pseudogenes based on the amino acid alignment and blast results. The genomic coordinates of the newly identified IGHC genes are listed in Supplemental Table 2.

Table 1

List of mammals genome included in this study.

Classification	Species name	Assembly	Date	Coverage and sequence technology	Contig N50 (×103)9))
Cetartiodactyla	Tursiops truncatus	turTru2	Jan.2012	2.5x Sanger; 3.5 × 454 FLX; 30x Illumina HiSeq	11.8 (2.5)
	Orcinus orca	Oorc_1.1	Jan.2013	200x Illumina HiSeq	70.3 (2.4)
	Neophocuena asiaeorientalis	Unreleased
	Lipotes vexillifer	Lipotes_vexillifer_v1	Jul.2013	115x Illumina Hiseq 2000	31.9 (2.4)
	Physeter catodon	Physeter_macrocephalus-2.0.2	Sep.2013	75x Illumina	35.3 (2.3)
	Balaena acutorostrata	BalAcu1.0	Oct.2013	92x Illumina HiSeq 2000	22.7 (2.4)
	Balaena mysticetus	http://www.bowheadwhale.org/	Jan.2015	150x Illumina HiSeq	34.8 (2.3)
Carnivora	Leptonychotes weddellii	LepWed1.0	Mar.2013	82x Illumina HiSeq	23.7 (3.2)
	Odobenus rosmarus	Oros_1.0	Jan.2013	200x Illumina	90.0 (2.4)
Chiroptera	Pteropus alecto	ASM32557v1	Dec.2013	110x Illumina HighSeq 2000	31.8 (2.0)
	Myotis lucifugus	Myoluc2.0	Sep.2010	7x Sanger	64.3 (2.0)
Eulipotyphla	Erinaceus europaeus	EriEur2.0	Sep.2012	79x Illumina Hi-Seq	21.4 (2.7)
	Sorex araneus	SorAra2.0	Aug.2012	120x Illumina Hi-Seq	22.6 (2.4)
Superorder Afrotheria	Trichechus manatus latirostris	TriManLat1.0	Jan.2012	150x Illumina HiSeq	37.8 (3.1)
	Loxodonta africana	Loxafr3.0	Jul.2009	7x Sanger; ABI	69.0 (3.2)

Sequence alignment and phylogenetic analysis

To confirm that the regions of alignment had sufficient confidence to test for positive selection and determine phylogenetic relationships, the unaligned multiple sequence alignment (MSA) of IGHC genes translated to amino-acid sequences and further back-translated into the corresponding codon were submitted to GUIDANCE (Penn et al., 2010a, b) to obtain alignment confidence scores. GUIDANCE used a PRANK alignment program (version 140,110) and returned a colored MSA that allowed to delimit aligned regions ambiguously. Then the columns with confidence scores lower than the threshold were excluded from further analyses.

Phylogenetic trees were reconstructed using Bayesian inference (BI) in MrBayes 3.2.2 (Huelsenbeck and Ronquist, 2001) and Maximum likelihood (ML) in RAxML 8.0.26 (Stamatakis, 2006). The best-fit evolutionary model for each dataset was chosen by MrModeltest 2.3 (Nylander, 2004) with the Akaike Information Criterion (AIC). Bayesian inference analysis was run for 106 generations with one cold and three heated Markov chains to ensure the analysis remaining in global rather than local optima and trees were sampled every 100 generations. The first 1000 trees were discarded as the "burn-in". Because the samples from the beginning of the chain are unreliable. For ML analysis, the maximum likelihood tree was reconstructed with 1000 bootstrap replications.

Detection of positive selection

The rate of nonsynonymous (d N) and synonymous (d S) substitution (d N/d S or omega, ω) measures the direction and intensity of natural selection, with ω = 1, ω > 1 and ω < 1 indicating neutral evolution, positive and purifying selection, respectively. To detect the impact of natural selection on IGHC genes, we performed a number of statistical tests based on calculating the ω values using the CODEML program implemented in PAML 4.7 package (Yang, 2007). The selection analyses using the gene tree reconstructed by BI and ML methods generated nearly the same results as those from the well-supported species phylogeny (Perelman et al., 2011; Zhou et al., 2011a); thus, only the latter analysis is shown here.

To detect positively selected sites on IGHC genes, we first performed the site specific models, i.e., M8a (nearly neutral, beta distribution: 0 < ω0 < 1 and ω1 = 1) versus M8 (positive selection, beta distribution: 0 < ω0 < 1 and ω1 > 1), as implemented in CODEML program. In addition, positively selected sites were identified using Datamonkey, which incorporated the rate of synonymous substitution, whereas the d S was fixed in PAML. Three ML methods, i.e. the single likelihood ancestor counting (SLAC), fixed-effect likelihood (FEL), and random-effect likelihood (REL) implemented in the Datamonkey website (Pond and Frost, 2005a), were used to infer the positively selected sites with default settings and significance levels of 0.1, 0.1, and 50, respectively. Although SLAC, a conservative method, is suitable for large data sets, it usually underestimates the substitution rate. FEL method seems to capture the pattern of rate variation better than SLAC and REL, but it tends to be less conservative than SLAC to detect selection in data sets of intermediate size. REL approach is more powerful than SLAC and FEL, but it has the higher rates of false positives for small data sets (Pond and Frost, 2005a,b). These three predictions were conducted to estimating ω value of each site.

To further detect whether the possibility of positive selection was limited on specific lineages, we compared the free-ratio model, which estimated a separate ω ratio for each lineage, with the one-ratio model that estimated a same ω ratio for all lineages in the tree (Yang, 1998). In addition, the branch-site model were conducted to detect evidence of positive selection on individual codon along a specific lineage (Zhang et al., 2005). In branch-site model, we compared Ma (positive selection, 0 < ω0 < 1, ω1 = 1 and ω2 ≥ 1) with Ma0 (neutral selection, 0 < ω0 < 1, ω1 = 1 and ω2 = 1) to detect positive selection. The significance of the differences between two nested models was tested using likelihood ratio tests (LRT) statistic (-2[LogLikelihood1-LogLikelihood2]) with a chi-square distribution. Bayes Empirical Bayes (BEB) approach (Yang et al., 2005) implemented in the CODEML was used to test the posterior probabilities (PPs) of positively selected sites, where sites with PPs > 0.90 were considered as candidates undergoing positive selection. Finally, to support the PAML results, we employed a complementary protein-level approach implemented in TreeSAAP (Woolley et al., 2003) that detected the significant changes in amino acid physicochemical properties.

Results and discussion

Genomic organization of IGHC genes in marine mammals

In the present study, we retrieved IGHC genes from the genomes of ten marine mammals (including seven cetaceans, one manatee and two pinnipeds) and compared them with those of terrestrial relatives to provide insights into the evolution of IGHC genes in marine mammals. We found that the genomic regions containing an IGHC cluster are conservative with special flanking sequences, i.e., the TMEM121 and IGH diversity gene cluster (Eguchi-Ogawa et al., 2012), which helped us obtain the complete IGHC genes in mammalian genomes. A total of 182 IGHC genes were identified in the 26 mammals in our study (Fig. 1 ), including 79 sequences that were newly identified in 13 species. For the 10 marine mammals, a total of 59 sequences were identified in this study.

Fig. 1

Positive selection at IGHC genes across the mammal phylogeny and the presence/absence patterns of mammal IGHC genes.

A widely accepted phylogeny of mammals was used for selective pressure analysis in the free ratio model and the branch-site model. The branches indicated by vertical lines or circles show evidence of undergoing positive selection. The red font indicates marine mammals; black font indicates terrestrial mammals. The IGHC genes of corresponding species were listed in the right part of Fig. 1. Each circle represents a gene member; circles represent a complete gene sequence; semicircles indicate a partial gene sequence; red cross means no sequences found by TBLASTN and BLASTN approaches; black cross means the gene really lost; the species highlighted with a black line mean a pseudogene.

To define the orthologous relationships among IGHC genes, we reconstructed phylogenetic trees for each gene, as well as a phylogenetic tree incorporating all members of IGHC genes from 26 mammals using two different methods (BI and ML). Similar tree topologies were inferred (Fig. 2 , Supplemental Figs. 1–11). The phylogenetic relationship showed that the IGHC genes of marine mammals could be clustered well with those of their terrestrial relatives, suggesting that marine mammals possess the typical isotype repertoire as other eutherian mammals. Similar to other mammals, marine mammals were found to contain a single IGHM, IGHD, IGHE, and IGHA genes, and several IGHG genes based on the phylogenetic reconstruction. The order of IGHC genes in the IgH locus of marine mammals was determined to be 5′-IGHM-IGHD-IGHG-IGHE-IGHA-3′, which is in agreement with their terrestrial relatives (Supplemental Fig. 12). Furthermore, we found that the first domain (CH1) of IGHM and the CH1 of IGHD of cetaceans shared nearly identical sequences, as shown in artiodactylas (Zhao et al., 2002) (Supplemental Fig. 13). The genomic organization of marine mammal IGHC genes is similar to that of other species examined, which suggested a constraint acting on their genomic arrangement and the conserved genomic organization of the IGHC genes.

Fig. 2

Phylogenetic tree of mammal IGHG genes produced using MrBayes.

The CH domains of IGHG genes were used in the tree reconstruction. Tutr, Tursiops truncatus; Oror, Orcinus orca; Neas, Neophocuena asiaeorientalis; Live, Lipotes vexillifer; Phca, Physeter catodon; Baac, Balaena acutorostrata; Bamy, Balaena mysticetus; Lewe, Leptonychotes weddellii; Odro, Odobenus rosmarus; Trma, Trichechus manatus latirostris; Mylu, Myotis lucifugus; Ptal, Pteropus alecto; Ereu, Erinaceus europaeus; Soar, Sorex araneus; Loaf, Loxodonta africana; Bota, Bos taurus; Susc, Sus scrofa; Eqca, Equus caballus; Calu, Canis lupus familiaris; Aime, Ailuropoda melanoleuca; Hosa, Homo sapiens; Mamu, Macaca mulatta; Rano, Rattus norvegicus; Orcu, Oryctolagus cuniculus; Modo, Monodelphis domestica; Oran, Ornithorhynchus anatinus.

The expansion of the IgG subclass is primarily a feature of mammals (Hsu et al., 2006). The IgG-like Igs of chicken and lizard showed no subclass diversification (Butler et al., 2011). The gene duplication of IGHG genes is believed to have occurred most recently, and it belongs to a lineage-specific event (Wagner et al., 2002). In the present study, we found that there are usually two or four IGHG genes in marine mammals (Fig. 1). Like other mammals, the IGHG gene in marine mammals comprises three CH domains and a hinge segment. Furthermore, they have close sequence similarities to the IGHG genes in their terrestrial relatives. Based on the phylogenetic tree, it was found that IGHG genes within one species clustered together rather than with those of other species (Fig. 2). However, the IGHG genes of five species grouped together with their true orthologues in closely related species, including bottlenose dolphin, killer whale, Yangtze finless porpoise, Weddell seal and Pacific walrus, likely due to the species of Delphinidae (e.g., dolphin, killer whale) rapidly diverging during a very narrow time frame of ∼3 million years (Zhou et al., 2011b) from the last common ancestral IGHG genes. These observations provided evidence to support the view that species-specific duplication of IGHG genes occurred during the evolution of mammalian species (Wagner et al., 2002). Overall, the presence of IgG subclasses and species-specific duplication of IGHG genes in marine mammals further confirmed the finding that IGHC genes of marine mammals fit the canonical pattern of their terrestrial relatives.

The number of IGHC genes has been noted to vary among different mammals. For example, 2 IGHA genes were present in humans and platypuses (Mage et al., 2006), whereas 13 IGHA genes were found in rabbits (Mccallum et al., 2004). For IGHG genes, seven copies were found in horses (Wagner et al., 2004), five in humans (Takahashi et al., 1982), and only one in rabbits (Mage et al., 2006) and opossum (Wang et al., 2009). No IGHG gene was detected in the bowhead whale and little brown bat (Fig. 1), which may due to the bad genomic assembly of both species considering that IgG subclasses have been identified in other whales (Mancia et al., 2006) and partial sequence encoding of IgG in little brown bat (Myotis lucifugus) was recovered by PCR (Butler et al., 2011). We also failed to retrieve IGHD genes from the genomes of sperm whale, hedgehog, European shrew, little brown bat, black flying fox and elephant despite considerable efforts (Fig. 1). One possible explanation is the poor genome quality, and the an alternative explanation is that the IGHD gene might have really been lost in evolution, similar to the situation in opossums (Wang et al., 2009) and rabbits (Mage et al., 2006). Further investigations to recover IGHD transcript by PCR will be needed in the future. Despite necessary further verification of IGHC genes in some species(e.g. with transcript evidence), changes in the copy number of IGHC genes among different mammalian lineages might reflect the rapid evolution of Ig genes (Hsu et al., 2006), corresponding to the marked variability found in the pathogens in different environments.

Adaptive evolution of IGHC genes in mammals

To detect the possible roles of positive selection during the evolution of IGHC genes in mammalian lineages, evolutionary analyses were conducted using a series of selection models implemented in the CODEML program and DataMonkey website. Because of the absence of intact IGHD genes in many mammals and the species-specific duplication of IGHG genes, only datasets of the IGHM, IGHE and IGHA genes were used to test the selection. Site model analysis of all mammals revealed that the models incorporating positive selection (M8) fit the data significantly better than the neutral model (M8a) for both the IGHM and IGHA genes, whereas no evidence of positive selection was identified at the IGHE genes (Table 2 ). Specially, M8 model detected nine positively selective sites at the IGHM and IGHA genes, whereas a total of 39 codons from the IGHM, IGHE and IGHA genes were also found to be under positive selection using the FEL, REL, and SLAC methods implemented in Datamonkey. Eleven codons from three genes (5 in IGHM, 2 in IGHE, 4 in IGHA) were predicted to be robust sites under positive selection identified by at least two ML methods. Furthermore, TreeSAAP analysis showed that these 11 positively selected sites have undergone radical changes in their physicochemical properties (Table 3 ), and this may be regarded as another evidence to support significant positive selection. Then, we marked the putatively selected sites in the alignment and found that 36.4% (4/11) of these sites were located in, or close to, the functionally important regions of IGHC genes (Supplemental Fig. 14). This result suggested that positive selection on mammals may be the major driving force for the evolution of IGHC genes to adapt to environmental pathogens. Besides, for expanded subclasses, such as the 13 IGHA genes of rabbit and multiple copies of IGHG gene in mammals, we tested whether this species-specific expansion has been influenced by positive selection. The site model analysis (M8 vs. M8a, M1a vs. M2a implemented in PAML) showed that the IGHA genes of rabbits were under positive selection with similar positively selected sites (Supplemental Table 3). The branch-site model was used to examine whether some branches leading to gene duplication of IGHG genes (branch a-o) have been positively selected. The results showed that signals of positive selection were detected in branches a (Delphinoidea), c (Physeter catodon), e (Bos taurus), f (Sus scrofa), g (Equus caballus), and h (Macaca mulatta) (Supplemental Table 4, Supplemental Fig. 15). Overall, pervasive positive selection was detected on IGHC genes of mammals.

Table 2

Selective pressure analyses of mammal IGHC genes by site model.

Gene	Models	-lnLa	P value	M8 ω value	Positively Select Sitesb(P > 0.90)
IGHA	M8	10704.474		2.045	7 T 0.997; 203 A 0.937; 213 A 1.000; 237 R 1.000; 254 T 0.961; 321 N 0.977
	M8a	10716.250	<0.001
IGHE	M8	13389.996	0.065	1.371
	M8a	13391.701
IGHM	M8	14298.945	0.006	1.888	3 P 0.946; 101 Q 0.974; 294 L 0.990
	M8a	14302.713

lnL is the log-likelihood score.

PPs of Bayes Empirical Bayes (BEB) analysis with P > 0.9 was regarded as candidates for selection.

Table 3

Positively selected sites of mammal IGHC genes detected by PAML, Datamonkey and TreeSAAP.

Genes	AA positions	PAML		Datamonkey			TreeSAAP
		Branch-site model (P > 0.90)	Site model (P > 0.90)	SLAC (P < 0.1)	FEL (P < 0.1)	REL (BF > 50)	Radical Changes in Amino Acid Propertiesa	Total
IGHA	200			0.0634811	0.0304991	77.8394	pHi, αc, Et,	3
	213		1.000		0.0406899		Pα, Ns, Bl, Br, RF, Pc, K°, h, F, Pr, αc, αm, αn, Ra, Hp, Ht, Et, P	18
	328			0.0351638	0.0935704		Pr, αc, Et	3
	364			0.0610423		52.2109	Pα, Ns, Pβ, Br, RF, Pc, K°, h, F, Pr, p, αc, αn, Ra, Hp, Ht, Et, P	18
	365			0.034845	0.0453226	52.0384	Pα, Pc, K°, Pr, p, αc, Et, P	8
IGHE	201			0.0950233	0.098536		Ns, Pβ, Br, RF, Pc, pK’, Ca, h, pHi, El, F, Hnc, V°, Pr, p, Esm, Ra, Hp, Ht, Et, P	21
	224			0.099105	0.00221453	312.686	Pα, K°, Ht	3
IGHM	48			0.0880392	0.0671209		Pα, Ns, Bl, RF, Pc, pHi, F, αc, αn, P	10
	254			0.0585277	0.0527633		Pα, K°, Ht	3
	294		0.990		0.0526114		Pα, Pβ, Pc, K°, pK’, pHi, Pr, αc, αn, Ra, P	11
	345			0.0287169	0.0361345	57.5523	Pα, Ns, Pr, p, αc, Ra, Ht, P	8

Radical changes in amino acid properties under categories 6–8 were detected in TreeSAAP. = α-helical tendencies, = Average number of surrounding residues, = β-Structure tendencies, = Bulkiness, = Buriedness, = Chromatographic index, = Coil tendencies, ° = Compressibility, = Equilibrium constant of ionization for COOH, = Helical contact energy, = Hydropathy, = Isoelectric point, = Long-range non-bonded energy, = Mean r.m.s. fluctuation displacement, = Normal consensus hydrophobicity, ° = Partial specific volume,  = Polar requirement, = Polarity, = Power to be C-terminal, α-helix, = Power to be middle, α-helix, = Power to be N-terminal of an α-helix, = Short- and medium-range non-bonded energy, = Solvent accessible reduction ratio, = Surrounding hydrophobicity, = Thermodynamic transfer hydrophobicity, = Total non-bonded energy, = Turn tendencies.

As is well known, many physicochemical properties, such as temperature, pressure, illumination and osmotic pressure, are different between water and land environments, possibly causing the difference in the diversity and abundance of pathogens. To date, there has been no deterministic conclusion on whether the pathogenic microbes in marine environment or in terrestrial environment are more diverse. Some early studies have supported weaker pathogenic pressure in the water environment (Zhang et al., 2016; Trowsdale et al., 1989; Slade, 1992), while several recent studies have examined the microbiome of marine mammals, suggesting that although the taxonomic composition of microbial communities in marine mammals is distinct from those of terrestrial mammals, its diversity is not lower than that of terrestrial mammals (Bik et al., 2016; Nelson et al., 2015).

To further test whether marine and terrestrial mammals were subjected to different pathogenic pressures during their evolution, we used the branch model in our study. The LRTs showed that the free-ratio model fits better than one-ratio model for IGHM, IGHE and IGHA genes, indicating heterogeneous selective pressures on different lineages. The ω values greater than 1 are limited to some branches in marine and terrestrial mammals, including the terminal branch of minke whale, the last common ancestral (LCA) branch of baleen whales, the LCA branch of Glires, and the LCA branch of Euarchontoglires at IGHA gene; the terminal branch of bowhead whale, and the LCA branch of Eulipotyphla at IGHE genes; the LCA branch of Laurasiatheria at IGHM gene (Fig. 1), suggesting similar selective pressures acting on IGHC genes in marine and terrestrial mammals. This result could provide evidence to support that the diversity and abundance of pathogens in water was comparable to that in land.

We also used more stringent branch-site model to predict positive selection acting on specific sites in each lineage. Almost all signs of positive selection were identified in terrestrial lineages (e.g., artiodactyla, perissodactyla, chiroptera, eulipotyphla, primates, Glires, Fig. 1, Supplemental Table 4). By contrast, for marine mammal lineages, only the terminal branch of Florida manatee was subjected to positive selection (Fig. 1, Supplemental Table 4). This result supported that pathogenic pressure was weaker in aquatic environment than in terrestrial environment. Moreover, recent studies showed a very low number of variable genes that directly binds to antigens in marine mammals, which also supported this conclusion (Olivieri et al., 2014). But further study on the variable region will be needed, providing insights to understand whether the difference in selective pressure between marine and terrestrial mammals existed and whether such a difference could respond to the different pathogenic pressure between aquatic and land environments.

In the present study, IGHC genes were determined to have undergone positive selection in Chiroptera, which is now known as notorious reservoir hosts for a number of pathogenic viruses, including Nipah, Hendra, Ebola, and severe acute respiratory syndrome (SARS) responsible for severe human and livestock disease outbreaks. However, they rarely display clinical symptoms (Zhang et al., 2013; Wang, 2009). The ability of bats to coexist with viruses presents an interesting immunological problem. Previous studies have founded that some immunoglobulin gene superfamilies underwent expansion in Brandt’s bat (Seim et al., 2013) and that some positively selected genes were involved in the immune response (Zhang et al., 2013). Thus, IGHC genes subject to positive selection in chiropterans may be related to the long-term coexistence of bats and viruses.

Evidences for positive selection has also been detected in some domesticated mammals such as pig, horse and cattle. Approximately 10,000 years ago, the domestication of farm animals by selection for biological, agricultural and biomedical importance began, leading to animals that slowly adapted to local condition. Approximately 200 years ago, selection pressures on livestock increased with the rise of the concept of breeds (Taberlet et al., 2008). It has been previously proposed that artificial selection indirectly drove for rapid adaptation of domestic animals and pathogen genes (Smith et al., 2012). Therefore, we postulated that artificial selection resulted in the adaptive evolution observed in IGHC genes of domestic animals.

Conclusions

This study has provided the first characterization of IGHC genes in marine mammals and explored their potential molecular genetic basis in the secondary aquatic adaptation. Our analyses revealed that the IGHC genes in marine mammals fit the canonical pattern of their terrestrial relatives. In addition, different levels of selection were detected between marine and terrestrial mammalian lineages with the branch-site model, suggesting a divergent adaptation to land and aquatic environments probably due to contrasting pathogenic pressures. Further investigation on the variable region of IgH genes is necessary to improve our understanding of the immune adaptation of marine mammals.

Funding

This work was supported by the Key Project of National Natural Science of China (grant number 31630071 to GY); the National Natural Science Foundation of China (NSFC) (grant numbers 31570379; 31772448 to SX); and the National Key Research and Development Programme, Ministry of Science and Technology (Grant number 2016YFC0503200 to GY and SX; the Priority Academic Program Development of Jiangsu Higher Education Institutions to GY and SX).

Declarations of interests

None.