Adaptations to Endosymbiosis in a Cnidarian-Dinoflagellate Association: Differential Gene Expression and Specific Gene Duplications

Summary

Trophic endosymbiosis between anthozoans and photosynthetic dinoflagellates forms the key foundation of reef ecosystems. Dysfunction and collapse of symbiosis lead to bleaching (symbiont expulsion), which is responsible for the severe worldwide decline of coral reefs. Molecular signals are central to the stability of this partnership and are therefore closely related to coral health. To decipher inter-partner signaling, we developed genomic resources (cDNA library and microarrays) from the symbiotic sea anemone Anemonia viridis. Here we describe differential expression between symbiotic (also called zooxanthellate anemones) or aposymbiotic (also called bleached) A. viridis specimens, using microarray hybridizations and qPCR experiments. We mapped, for the first time, transcript abundance separately in the epidermal cell layer and the gastrodermal cells that host photosynthetic symbionts. Transcriptomic profiles showed large inter-individual variability, indicating that aposymbiosis could be induced by different pathways. We defined a restricted subset of 39 common genes that are characteristic of the symbiotic or aposymbiotic states. We demonstrated that transcription of many genes belonging to this set is specifically enhanced in the symbiotic cells (gastroderm). A model is proposed where the aposymbiotic and therefore heterotrophic state triggers vesicular trafficking, whereas the symbiotic and therefore autotrophic state favors metabolic exchanges between host and symbiont. Several genetic pathways were investigated in more detail: i) a key vitamin K–dependant process involved in the dinoflagellate-cnidarian recognition; ii) two cnidarian tissue-specific carbonic anhydrases involved in the carbon transfer from the environment to the intracellular symbionts; iii) host collagen synthesis, mostly supported by the symbiotic tissue. Further, we identified specific gene duplications and showed that the cnidarian-specific isoform was also up-regulated both in the symbiotic state and in the gastroderm. Our results thus offer new insight into the inter-partner signaling required for the physiological mechanisms of the symbiosis that is crucial for coral health.

Introduction

The mutualistic symbiosis of anthozoans (Cnidaria), such as corals and sea anemones, with their intracellular dinoflagellate symbionts, mostly of the genus Symbiodinium, forms both trophic and structural foundation of coral reef ecosystems. Anthozoans have a very simple body plan and are composed of two tissue layers, the epidermis and the gastroderm (also called ectoderm and endoderm, respectively [1], Figure S1A). They host their unicellular symbionts, also called zooxanthellae, inside vacuoles (symbiosomes) within the gastrodermal cell layer. Safely localized inside the animal host cells, the photosynthetic symbionts fix large quantities of carbon dioxide. Most of the reduced organic carbon produced is then translocated to the host as mobile compounds, such as glycerol, lipids and amino acids [2]. In return, the host provides suitable conditions for symbiont photosynthesis: inorganic nitrogen, phosphorus and inorganic carbon, as well as a favorable high light environment [2].

This partnership generates many constraints, however, resulting in physiological and cellular adaptations (for review, [3], [4]). For example, the presence of photosynthetic zooxanthellae within the gastrodermal cells requires the host to transport inorganic carbon from the surrounding seawater to the symbionts, a process in which carbonic anhydrases (CAs) are central [5], [6]. To prevent possible cellular damage resulting from photosynthesis-induced hyperoxia, both partners also express a high diversity of antioxidant enzymes, including catalases, peroxidases and superoxide dismutases (SODs), [7]–[9]. SOD and CA isoforms specific to symbiotic anthozoans have been demonstrated [8], [10].

Environmental perturbations such as an increase in seawater temperature may induce dysfunction and collapse of the symbiosis, leading to zooxanthellae loss or so-called “bleaching”, and this phenomenon has led to severe worldwide decline of coral reefs [11]. The establishment and maintenance of this partnership must therefore be dependent on intimate molecular communications between the partners, including recognition and tolerance of symbionts, as well as adaptations for mutual transport and exchange of nutritional resources. A small number of candidate genes for this molecular dialogue has been proposed, including CAs and the cell adhesion protein Sym32 [10], [12], [13]. These have so far been examined by targeted protein analyses (for review, see [14]). Large-scale gene expression studies have tended to highlight a modulation of the host transcriptome, in particular genes involved in cell adhesion, lipid metabolism, cell cycle regulation, or cell death [15]–[17]. While most transcriptomic approaches have been performed in symbiotic cnidarians under thermal stress, imposed in order to understand the molecular and cellular basis of the early events leading to bleaching [18], [19], we rather focused our present experimental approach on two stable states: symbiotic and aposymbiotic individuals.

How do anthozoans maintain a stable partnership with their photosynthetic dinoflagellate symbionts? To decipher the molecular dialogue inferred by the presence of the dinoflagellate symbionts (Symbiodinium clade temperate A) within the sea anemone Anemonia viridis, we compared transcriptomes of symbiotic and aposymbiotic specimens using a symbiosis-dedicated microarray. This oligonucleotide microarray (2,000 features) was developed from the A. viridis 40,000 EST collection [20] and is dedicated to genes potentially involved in symbiosis regulatory pathways.

Our two main goals were to identify gene expression patterns characteristic of: i) the symbiotic and aposymbiotic conditions, ii) the epidermis and gastroderm tissue layers. The main advantage of our model, A. viridis, is it allows efficient separation of the two tissue layers of the animal (epidermis and gastroderm) with minimal cross contamination [8]. Concerted DNA microarrays and quantitative RT-PCR (qPCR) analyses outlined characteristic gene expression signatures for the symbiotic and aposymbiotic states. Among newly identified genes, some appeared to result from anthozoan-specific gene duplications. Analyses of tissue-specific expression demonstrated that most of the host genes involved in symbiotic interactions are preferentially expressed in the gastroderm (i.e. the zooxanthellate tissue layer). We detailed several genetic regulatory pathways involved in dinoflagellate-cnidarian recognition, carbon transfer to the intracellular symbionts and mesoglea constitution. Finally, we propose a model where the aposymbiotic state triggers vesicular trafficking whereas the symbiotic state favors metabolic exchanges between host and symbionts.

Results

Molecular quantification of dinoflagellates

The “bleached” phenotype is shared by many stressed symbiotic cnidarians and is the consequence of a massive upstream loss of symbionts and hence their photosynthetic pigments. Most approaches to quantify symbiont loss use manual zooxanthellae counts. Here, we developed a fast and accurate approach using real-time quantitative PCR on total genomic DNA extracts to quantify the relative number of nuclei and hence cells.

Symbiodinium (temperate A clade) nuclear Elongation Factor 2 (EF2), Sucrose Phosphate Synthase (SPS) and Ascorbate Peroxidase (APX) and the A. viridis nuclear Coatomer subunit gamma (COP-γ Regulator of Chromosome Condensation protein 2 (RCC2) and Niemann-Pick disease type C1 (NPC1) gene copy numbers were assessed for total genomic DNA extracts from 5 symbiotic anemones (Sy1–Sy5) and 6 aposymbiotic anemones (AS1–AS6) collected along the Mediterranean coasts around Nice, France (Figure S2) and kept in laboratory culture, as well as in vitro cultured Symbiodinium (CZ) and epidermal tissue fraction (Ep). Although gene locus number per nucleus is assumed to remain constant for a given species, dinoflagellates have been shown to undergo gene specific amplifications [21], [22]. In Amphidinium carterae, the EF2 gene is present in tandem repeat contrary to APX [21]. Whilst the 3 A. viridis gene loci were present in a nearly 1∶1 ratio in all specimens, in Symbiodinium the number of genes was relatively variable from one specimen to another despite the proximity of collection sites to one another. The SPS∶APX ratio was around 1∶1 in all specimens except for AS2 and AS6. EF2∶APX and EF2∶SPS ratios displayed strong polymorphism, likely due to EF2 variable gene amplification (Figure S3A). Gene copy numbers between A.viridis and its symbiont were also measured. Figure 1 shows the relative symbiont to host nucleus ratio for APX/COP-γ SPS/COP-γ and EF2/COP-γ (similar results were obtained with RCC2 and NPC1, Figure S3B). Whatever gene locus used, the pattern of symbiont to host cell ratio was reproducibly similar with ratios between aposymbiotic and symbiotic specimens ranging from 8.2×10⁻² (AS4 versus Sy1) to over 10⁻⁴ fold (AS6 versus Sy5).

Inter-individual variability in symbiotic and aposymbiotic expression profiles

In order to identify gene expression patterns characteristic of symbiosis, the gene expression profiles of symbiotic and aposymbiotic anemones were compared using a symbiosis-dedicated oligoarray. 60-mer oligonucleotides were designed from 2,000 sequences putatively involved in symbiosis selected from the large clustered and annotated A. viridis EST collection [20]. Figure S1 gives an overview of the Gene Ontology (GO) functional annotations of the selected sequences.

Microarray experiments were performed between symbiotic and aposymbiotic anemones, using a dual-dye protocol, one condition being labeled with Cy3, and the second by Cy5. The experimental design of hybridizations is shown in Figure S4. Microarray results were analyzed according to two methods. First, in order to have an overall estimate of the genes differentially expressed between the two states, we treated the 11 individual specimens as two batch categories: symbiotic anemones (Sy1–Sy5) and aposymbiotic anemones (AS1–AS6). Statistical analysis was performed with the limmaGUI package [23] that defined for each gene in each experiment an average fold change (M, corresponding to a log2 ratio between the 2 experimental conditions) and a statistical value (called B, with positive values for the more significant variations). 58 and 78 genes were significantly up-regulated (|M|>0.59 and B>0) in the symbiotic (hereafter called “SY gene set”) and aposymbiotic states (hereafter called “APO gene set”), respectively (Figure 2A, Tables S1 and S2). Functional annotation of these 136 genes using GO terms and statistical analysis by Gossip (Fisher's exact test, p<0.05) showed that the terms “plastid”, “calcium ion binding”, and “protein modification process” were over-represented in the SY gene set, in contrast to “cytosol”, “cytoplasmic membrane-bounded vesicle” and “transcription”, which were over-represented in the APO gene set (Figure S1C).

Second, in order to gain insights about inter-individual variation in gene expression, we compared each individual symbiotic anemone (Sy1 to Sy5) with the batch of 6 aposymbiotic anemones (Apo), and inversely, compared each aposymbiotic anemone (AS1 to AS6) with the batch of 5 symbiotic anemones (Sy). This analysis (Figure 2A) showed unexpectedly high inter-individual variability of differential gene expression patterns within these 136 genes. Figure 2B shows the number of individual anemones for which a given gene is significantly differentially expressed between the SY and APO gene sets. For instance, only 8 SY genes were differentially expressed in all 11 anemones. Conversely, one APO gene (M = −0.63; B = 1.15) encoding for Interferon Regulatory Factor 1 was significantly differentially expressed in only 3 out of 11 anemones. To obtain an indication of the source of variability, we clustered the different anemones according to their expression. Different clustering parameters were tested and showed similar cluster trees; the most representative is shown underneath the heatmap (Figure 2A). Whatever the parameters used, specimens Sy4 and Sy5, AS1 and AS2, and AS4 and AS5 always paired together. Remarkably, even if the anemones were maintained in the same culture conditions at least 3 months before sampling, such co-clustering seemed to correlate with their previous life history: different collection dates/areas for symbiotic specimens, and different bleaching causes for aposymbiotic sea anemones (see Material and Methods and Figure S2). Nonetheless, these variable transcriptomic profiles define the same stable symbiotic or aposymbiotic phenotype.

Definition of the Kern gene subset

The “Kern” gene subset (named after the German word for nucleus) was where the remainder of the analysis was focused, and is defined as those genes of cnidarians origin only (excludes genes of unicellular or prokaryotic origin) found in the SY or APO gene set, and differentially expressed in at least 8 out of the 11 specimens. Table 1 lists the representative 19 and 20 Kern genes found in the SY and APO gene sets respectively. A blast homology search against the Nematostella vectensis genome database and the Uniprot generalist database, and searches for specific protein signatures (signal peptide, trans-membrane domains or others domains) allowed the genes to be named and functions and cellular localizations assigned (Table 1). Among the different functional categories found, cell adhesion proteins (8/39) were the most represented, a process indeed expected to play a key role in signaling events between partners. A total of 7 other genes were involved in metabolism (4 SY and 3 APO). The four SY genes were specifically involved in fatty acid metabolism, indicating that fatty acid metabolism is likely a preponderant metabolic pathway of the symbiotic condition. Although all genes listed in Table 1 would merit further investigation, only several selectively targeted genes are detailed in this study (see below).

In order to obtain support for the functional implication of some of the Kern genes, we monitored the expression of genes of interest directly after thermal stress, which is known to disrupt symbiosis. Anemones Sy3–5 were subjected to an 8°C temperature increase, and the expression of CA2-c, CA2-m and NPC2-D dropped by around 3 fold after 24 h and 48 h (Figure 3). This result is interesting, considering their up-regulation in symbiotic individuals. This immediate response to environmental stress precedes symbiosis breakdown and thus strengthens support for the possible role of these 3 genes in endosymbiosis in A. viridis.

Gastroderm-specific gene expression

As the two tissue layers can be separated in A. viridis, it is a powerful biological model for studying tissue-specific gene expression. We used two experimental designs to compare gastrodermal versus epidermal gene expression profiles (Figure S4). Firstly we directly hybridized the cDNAs of the two tissues against each other on the same array, and secondly we compared each tissue sample to the aposymbiotic AS6 reference and thus defined tissue expression profiles by transitivity. Only the genes with at least 1.50 fold increase in transcript abundance (M>0.59) in one tissue in both experiments were assigned as differentially expressed (Table 2, Tables S1 and S2). As expected, the expression of most Symbiodinium genes was restricted to the gastroderm. The Venn diagram in Figure 4 summarizes distribution of the 1,715 cnidarian-specific genes (prokaryote and zooxanthellae genes were excluded from the analysis) according to their gastrodermal (Ga), epidermal (Ep), and preferential symbiotic (SY) versus aposymbiotic (APO) expression. We further confirmed preferential expression of several genes in the different tissues using real-time quantitative PCR (Figure 5). Taken together, microarray and qPCR results showed that many more A. viridis genes were preferentially expressed in the gastroderm (71%) than in the epidermis (29%). More interestingly, among the 17 genes which were both preferentially expressed in a given tissue and differentially regulated under symbiotic/aposymbiotic status, the large majority (12) were in fact up-regulated within the gastroderm of symbiotic anemones (i.e. the zooxanthellate tissue). This suggests that the presence of symbionts directly modulates the gene expression of their hosting gastrodermal cells. In addition, out of these 12 genes, 9 belong to the Kern subset of genes, supporting the importance of this gene set in inter-partner communication and regulation (Figure 5 and Table 2). Hence, we have identified different categories of genes potentially involved in sequential aspects of symbiosis, and showed for the first time that the animal tissue-specific molecular response to the presence/absence of zooxanthellae was restricted almost entirely to preferential expression of genes within the compartment hosting the symbionts. These results support the crucial role of the gastroderm in this symbiotic interaction.

The gastroderm-specific gene expression disclosed interesting hallmarks. Carbonic anhydrases are primordial enzymes necessary for the transport of inorganic carbon through biological membranes. Two isoforms were identified in our library and had increased transcript abundance in zooxanthellate anemones. Based on their respective protein signatures and blast homology, one isoform (Av_CA2-c, Kern # 3) is cytoplasmic whereas the other isoform (Av_CA2-m, Kern # 5) probably localizes to the outer plasma membrane, since it contains a signal peptide sequence as well as a trans-membrane and GPI anchor domains at its NH₂ and COOH termini, respectively. The cytoplasmic CA2 was equally distributed in the two compartments, but the membrane-anchored CA2 was principally expressed in the gastroderm (Figure 5). Another interesting finding on the tissue distribution of gene expression was that collagen biosynthesis was broadly under the control of the gastroderm. In cnidarians, fibrillar collagen genes were recently shown to be much more represented than first expected from bilaterian evolutionary comparisons, with 8 different genes present in the genome of N. vectensis [24]. In A. viridis, we found 13 different cDNAs corresponding to portions of the N. vectensis homologs. Additionally, we monitored the expression of 5 genes regulating the post-translational processing of collagen synthesis. Interestingly, 14 of these cDNAs showed increased transcript abundance in the gastroderm (Table 2). The mesoglea, the acellular layer between epidermis and gastroderm, which mostly consists of collagen, may partly originate from the gastroderm in A. viridis. However, none of these genes showed differential expression between aposymbiotic and symbiotic states, except for Prolyl-4-hydroxylase alpha. This enzyme is a key chaperone for the biosynthesis of collagen, catalyzing the hydroxylation of proline residues of procollagen chains necessary for their correct three dimensional folding [25]. Two different isoforms were identified in our sequence library and the expression of one of them (Kern # 14) was specifically enhanced in zooxanthellate anemones. Thus, in addition to showing that the gastroderm synthesized much of collagen fibers, our results infer that Symbiodinium may exert, directly or indirectly, post-translational control on collagen synthesis and modulate the formation of the mesoglea in A. viridis.

Cnidarian-specific gene duplication of symbiosis-related genes

During the course of our analysis, we noticed that some of the genes involved in symbiosis (Kern genes) had multiple paralogs in our dataset. Such was the case with the previously described MERP gene family [20], but the same held true for other Kern genes including Niemann Pick type C2 (NPC2), Calumenin, Sym32 and C3 Complement (C3) families (Table 2). We conducted phylogenetic analyses on the NPC2, Calumenin and Sym32 gene families in order to assign an evolutionary origin for each member. A. viridis isoforms present in our dataset were compared to their homologs in N.vectensis and other representative eukaryotes with a complete genome sequence. Both maximum likelihood (Figure 6 and Figure 7A) and Bayesian (Figures S5, S6, S7) methods gave very similar trees.

In vertebrates, Calumenin belongs to the CREC gene family, which encompasses 5 members (CAB45, ERC-55, Reticulocalbin 1 and 3, and Calumenin) [26]. From phylogenetic analysis across metazoans (Figure 6A), 3 groups of homologs can be defined; the CAB45 homologs, the ERC55 homologs and the CALU (Calumenin, Reticulocalbin 1 and 3) homologs. Cnidarian homologs are found within these 3 groups, however additional cnidarian homologs (NvCALUa,c,f–h, AvCALUa,c) are also present outside the CALU and ERC55 groups, representing cnidarian specific CALU gene duplications.

NPC2 (Figure 6B) is a single copy gene in human, nematode, sea urchin, Placozoa and yeast and is independently duplicated in ascidian and fish. Two copies of the NPC2 gene were found in sea anemones. Phylogenetic analysis supports a gene duplication in Anthozoans, with NvNPC2a and AvNPC2a representing the orthologs of the chordate NPC2, and NvNPC2b and AvNPC2-D defining an anthozoan specific duplication.

In the symbiotic sea anemones A. viridis and Anthopleura elegantissima, Sym32 is composed of two adjacent FasI domains (Figure 7A). In A. viridis, we identified two FasI-containing proteins: Sym 32 and the related periostin gene (PN), which is conserved across metazoans (only cnidarians and human are shown). AvSym32 and AvPN are composed of 2 and 4 FasI domains, respectively. In N. Nematostella, there are also two comparable FasI-containing proteins; both are composed of 4 FasI domains. One corresponds to the PN homolog (NvPN) and the second, which we named Nv2Sym, is similar to a tandem duplication of the AvSym32 sequence with 2×2 FasI domains. In human, two related FasI-containing proteins are characterized: the cognate Periostin (HsPN) and the “Transforming growth factor-beta-induced protein ig-h3” precursor (HsBGH3), both with 4 FasI domains. As proteins did not have the same number of FasI domains, we conducted a phylogenetic analysis on the alignment of all single FasI domains in order to gain insight into the domain evolution of this protein family (Figure 7A and Figure S7B). Human and cnidarian PN genes likely evolved from a common ancestor while Hs_PN and Hs_BGH3 would have duplicated after the Cnidaria - Bilateria separation. Consequently, AvSym32 and Nv2Sym are cnidarian-specific genes. Whether a Sym32 version containing only 2 FasI domains is specific to symbiotic anthozoans remains to be clarified.

Comparison of the expression of the A. viridis CALU, NPC2 and sym32 homologs showed that the isoform which was more highly expressed in the symbiotic state was also preferentially expressed in the gastroderm, whereas the other members were ubiquitously expressed (Table 2). For instance, AvCalu-a (Kern # 1) was the most up regulated gene in the symbiotic condition and localized in the gastroderm, whereas expression of both the other isoforms (AvCalu-b and AvCalu-c) was neither different between symbiotic versus aposymbiotic conditions nor between the tissue layers (Table 2). The same held true for the two isoforms of NPC2 (AvNPC2-D (Kern # 2) versus AvNPC2-a) and for Sym32 (Kern # 4) versus the Periostin homologs (Table 2). Most interestingly, the isoform that was differentially expressed was always member of the Kern gene set.

Thus, analysis of N. vectensis gene families identified cnidarian-specific gene duplications (NvCALUa,c,f–h, NvNPC2b, Nv2sym) and A. viridis ortholog analysis showed that among the cnidarian-specific isoform, the Kern genes (AvCALUa, AvNPC2-D and AvSym32) were up-regulated both in symbiotic state and in the gastroderm.

Discussion

Many anthozoans rely on photosynthetic endosymbionts (mostly Symbiodinium sp.) to grow in oligotrophic environments. The symbiotic relationship involves regulatory crosstalk between partners that allows the association to persist. This interpartner communication includes: i) the recognition of the partners, ii) the ability of symbionts to colonize host cells without being rejected by the host immune system, iii) the regulation of symbiont population, and iv) adaptations for mutual transport and exchange of nutritional resources [27]. In order to identify the genes potentially involved in the molecular dialog supporting this endosymbiosis, we used a microarray approach to compare the gene expression profiles of 11 A. viridis anemones representative of the symbiotic and the aposymbiotic states. We identified a subset of A. viridis genes, which we named “Kern”, characteristic of the symbiotic (SY) and aposymbiotic (APO) states in the sense that they were preferentially expressed in symbiotic or aposymbiotic anemones. Consequently, the Kern genes are potentially important candidate genes for the maintenance of the symbiosis process. Since in A. viridis, the two tissue layers can be separated, we were also able to assign tissue-specific gene expression. Many genes or gene functions that were expressed more highly in symbiotic anemones were also expressed more highly in the gastroderm, the symbiont hosting tissue. Functional annotation associated with sub-cellular localization of the Kern gene products allowed us to draw a map of the 39 genes differentially regulated in the symbiotic and aposymbiotic states (Figure 8).

Autotrophy and heterotrophy

One of the first conclusions from our present model is that many genes that are present in higher abundance in the aposymbiotic state are associated with active vesicle trafficking, both in the endocytotic and secretory pathways, as opposed to the many membrane-bound proteins associated with cellular recognition and adhesion coded by genes in higher abundance in the symbiotic state (Table 1 and Figure 8). Under aposymbiotic conditions, transcription of several proteins controlling trafficking was markedly enhanced. These include TRAPPC2 and TRAPPC10 (Kern # 22 and 23), which are members of the TRAPP complexes that regulate endoplasmic reticulum (ER) to post-golgi vesicular trafficking [28], thus attesting for more active de novo biosynthesis in the aposymbiotic anemone than in the symbiotic host. As part of the same line of evidence, the highly differentially expressed gene AP-2 sigma-1 (Kern # 38) codes for a key protein in the formation of clathrin-coated vesicle [29], implying increased endocytic trafficking. Among other aposymbiotic marker genes, HVCN1 and SLC30a (Kern # 26 & 36) are two ion channels shown to be present in phagosomes and late endosomal/lysosomal vesicles, respectively [30], [31]. Of note, our dedicated array comprises several oligonucleotide probes for various rab mRNAs (including rab1a, 2a, 7a, 8a and 11a). Rab proteins have been shown to control vesicular trafficking in various species including symbiotic anemones [32]–[34] where the symbionts have been shown to detour host membrane trafficking resulting in failure of host lysosomes fuse with the symbiosomal membrane [34]. However, in our array results, none of the Av_rab genes were differentially expressed. This does not preclude a role of these Rab proteins in A. viridis symbiosis, as expression regulation may be at the translational level. Interestingly, the fact that vesicular trafficking appears to be much more active in aposymbiotic A. viridis corroborates with heterotrophy, which mostly rely on predation, digestion and de novo synthesis, as opposed to autotrophic symbiotic anemones, in which 60% of carbon flux is provided by symbionts [35]. Numerous studies have characterized glycerol, lipids and amino acids as the major mobile compounds transferred from the symbiont (reviewed in [2], [4]). At the level of metabolism pathways, Phosphoenolpyruvate carboxykinase (PEPCK-C; Kern # 28) had enhanced expression in aposymbiotic conditions. This enzyme catalyzes the rate-controlling step of gluconeogenesis when pyruvate is used as the substrate. Two metabolic pathways can initiate gluconeogenesis: the PEPCK-C dependent pathway, which is fuelled by pyruvate (product of glycolysis), and the PEPCK-C independent pathway, which is fuelled by glycerol [36]. As glycerol is acquired from the symbiont, PEPCK-C is indeed expected to be down-regulated in the symbiotic state.

In zooxanthellate anemones, lipids are the second most abundant transferred mobile compounds from the symbiont [2], [4] therefore requiring up-regulation by the host of several key regulatory enzymes for fatty acid metabolism. Although mobile compounds have been identified using radioactive tracers [2], compound transfer pathways from the symbiont remain to be characterized. Four Kern genes were involved in lipid processing pathways, LSD-2 (lipid storage droplet), SDR12 (terpenoid metabolism), 2-hydroxyacyl-CoA lyase 1 and NPC2. These 4 genes were up-regulated in the symbiotic state (Kern # 2, 13, 15, 17 in Table 1). 2-hydroxyacyl-CoA lyase 1 is responsible for fatty acid alpha-oxidation, a modification specific to fatty acids synthesized from chloroplastic organisms [37]. NPC2 (AvNPC2-D, Kern # 2) is among the most up-regulated genes, substantiating another study which noted that a NPC2 transcript was up-regulated in symbiotic Aiptasia pulchella anemones [38]. NPC2 acts in synergy with NPC1 in the transport of sterols (cholesterol) using the late endosomal/lysosomal system [39]–[41]. Noticeably, NPC1 and NPC2 are preferentially expressed in the gastroderm (Table 2). Based on their tissue expression profile and symbiosis-related regulation, NPC1 and NPC2 appear to be major candidate genes for the transport in A. viridis of sterol compounds produced by Symbiodinium.

Another notable difference between aposymbiotic and symbiotic states is that aposymbiotic-specific Kern genes showed no tissue expression preference. This is in marked opposition to the numerous symbiosis-related Kern genes expressed principally in the gastroderm. Hence, expression of the latter is both limited to the symbiont-containing cells and correlated with the presence of Symbiodinium. Although the mechanisms underlying such control are yet to be characterized, several molecular dialogs between host and symbiont can be envisaged.

Recognition and non-self tolerance by A. viridis

The Complement C3 (C3) is a precursor protein involved in adaptive immunity in vertebrates and also in cellular recognition, inflammatory processes and phagocytosis in invertebrates [42], [43]. C3 thus appears central to non-self response in metazoans and homologs have been identified throughout the metazoans, including cnidarians [42], [44]. In the coral Acropora millepora, one C3-like protein has been shown to localize in close association with the symbiosome in the gastroderm layer, supporting a role in recognition of Symbiodinium [44]. In sea anemones, there are three C3 isoforms in the non-symbiotic N. vectensis genome and a minimum of four different C3 isoforms were identified in the A. viridis EST dataset (PG and CS, personal communication). We monitored the expression of two of these: AvC3-1 and AvC3-2. Although AvC3-1 expression was mainly restricted to the gastroderm and may well represent the A. millepora functional homolog, no differential expression was observed between symbiotic and aposymbiotic anemones (validated by qPCR, Figure 5). On the other hand, AvC3-2 (Kern # 37) was expressed in both tissue layers, but was strongly repressed in the presence of symbionts. The functional divergence between the different C3 isoforms and their relative participation in Symbiodinium tolerance remains to be determined. However, based on their respective expression profiles in A. viridis, it is conceivable that recognition of and response to Symbiodinium may be carried out by different C3 paralogs.

Another recognition process that we uncovered in our experiment is related to the vitamin K-dependent (VKD) γ-carboxylation of Sym32, under the control of the Calumenin protein. In human, Calumenin proteins are composed of 6 to 7 Ca²⁺ binding EF hand domains. They are principally present in the Endoplasmic Reticulum (ER) due to specific ER retention motifs at their COOH termini [26], [45], [46]. Their differences in function are unclear, but they are associated with Ca²⁺ dependent processes, especially with post-translational VKD γ-carboxylation of several proteins [26]. The increase of intracellular vitamin K (a cofactor produced from plant) activates two proteins conserved throughout metazoans: Vitamin K1 2,3-Epoxide Reductase (VKOR) and γ-Carboxylase. The latter recognizes carboxylase recognition sites (CRS) in newly synthesized proteins and adds a CO₂ group to adjacent glutamic acid (Glu) residues resulting in the production of γ-carboxyglutamic acid (Gla)-containing proteins [47], [48]. VKD proteins have been mainly explored in mammals and include bone marrow proteins such as Osteocalcin, Matrix Gla Protein and Periostin (PN) which are involved in bone formation [49]. The vitamin K cycle is itself under the negative control of Calumenin [26], [46] (Figure 7B). In A. viridis, the Calumenin homolog AvCALU-a (Kern # 1) was the most up-regulated gene of the symbiotic state detected in the microarray experiments (Table 1). It is also preferentially expressed in gastrodermal cells. It contains an ER retention motif and is thus expected to down-regulate the VKD γ-carboxylation in the ER of zooxanthellate A. viridis cells.

Sym32 was first characterized by Weis and colleagues as a symbiosis-specific protein that is over-represented at the perisymbiotic membrane of the zooxanthellate anemone A. elegantissima (AeSym32), but also present at the surface of gastrodermal vesicles in the aposymbiotic state [13], [50]. In A. viridis, the Sym32 ortholog (AvSym32, Kern # 4) is up-regulated by more than 40 fold in the symbiotic state and is principally expressed in the gastroderm (Figure 5, Table 1 and Table 2). Both AeSym32 and AvSym32 are composed of two fasciclin I (FasI) domains (Figure 7), functionally associated with cell-cell recognition and thus potentially involved in anemone-zooxanthellae interaction [13], [51]. Importantly, in both sea anemones, the first FasI domain of Sym32 is highly similar in sequence to the first FasI domain of PN, and contains a sequence motif very similar to the CRS domain of PN characterized in vertebrates [49]. Sym32 is thus an implicit substrate for the γ-Carboxylase in sea anemones. As part of the same line of evidence, a 2D immunoblot of zooxanthellate A. elegantissima extract hybridized with an anti-Sym32 antibody showed the presence of two spots of 32 KDa, one at PI = 7.9 and the other at PI = 8.2 [51]. Such a PI difference for the same protein could very well correspond to difference between Glu and Gla containing Sym32 proteins since γ-carboxyglutamic acid decreases the PI of a protein (e.g. [52]). However, analysis of freshly isolated Symbiodinium (FIZ, where most of the peri-symbiosomal membrane - of host origin - remains attached around the isolated Symbiodinium cell) of A. elegantissima showed the disappearance of the PI = 7.9 spot in favor of the PI = 8.2 spot [13]. These experiments infer that the Glu-Sym32 protein would localize to the perisymbiotic membrane. Since γ-carboxyglutamic acid residues have been shown to modify the three dimensional structure and Ca²⁺ binding affinity, Glu-Sym32 and Gla-Sym32 would show different ligand properties [53]. Moreover, up-regulation of Calumenin shown in our work suggests an inhibition of γ-carboxylation which would favor the production of Glu-Sym32 in zooxanthellate A. viridis cells (Figure 7B). Thus, the vitamin K-dependent γ-carboxylation of Sym32, with its effect on interpartner recognition and the symbiotic process, is definitively a pathway that should be investigated.

Carbonic anhydrases and the transport of bicarbonate from seawater to the zooxanthellae

The photosynthetic symbionts are separated from the surrounding seawater by several host membranes: membranes of the epidermal cell layer, the collagenous basal membrane, gastrodermal cells, and perisymbiotic vesicles. The main source of inorganic carbon (Ci) for photosynthesis is seawater bicarbonate (HCO₃⁻), which implies transport of exogenous inorganic carbon through these layers of animal tissue [5], [54]. In seawater (pH 8.2), most inorganic carbon is in the form of HCO₃⁻, a form that needs carrier-mediation to cross membranes, and that is not readily converted to CO₂ in the absence of enzymatic action [55], [56]. The currently accepted model for external Ci uptake by the host involves an H⁺-ATPase acidifying the boundary layer where bicarbonate is converted to CO₂ by an external (likely membrane-bound) carbonic anhydrase isoform [3]. The uncharged CO₂ molecule then diffuses into the epidermal cell following the concentration gradient created by the extrusion of H⁺ in the external medium. Once in the animal cytoplasm, CO₂ is equilibrated with HCO₃⁻ according to the intracellular pH by another CA isoform, which prevents back-diffusion of CO₂ (for review see [3]). The mechanism of transport of Ci through the other membranes to the symbionts is currently debated (for reviews, see [5], [57], [58]). However, previous works have highlighted the role of a CA localized on the perisymbiotic [59] or algal membrane [60].

According to this model, carbonic anhydrases are crucial enzymes for carbon supply to symbiont photosynthesis. In the present study, the expression of two different carbonic anhydrases (Av_CA2-c and Av_CA2-m) was monitored. Both show highly enhanced expression in symbiotic specimens compared with aposymbiotic ones (4 and 2.9 fold, respectively), suggesting that both isoforms could be involved in the symbiosis/metabolic exchanges between partners. This result is consistent with the relevant work of Weis and her collaborators [6], [59], [61], showing that enzyme activity and transcript quantity are higher in symbiotic than in non-symbiotic specimens of sea anemones A. pulchella and A. elegantissima.

AvCA2-c isoform, the cytosolic isoform, is equally expressed in both tissue layers. We suggest that AvCA2-c catalyzes the intra-cellular reversible hydration/dehydration of CO₂ into HCO₃⁻ to facilitate the transport of CO₂ through the membranes and cells (Figure 8). AvCA2-m is a membrane-bound isoform specifically expressed in the gastrodermal layer of anemones where symbionts are located. This membrane-bound isoform can be located either on the plasma membrane of the gastrodermal cells, or on the perisymbiotic membrane surrounding the symbiont. In the first case, AvCA2-m would favor the transfer of CO₂ from one cell layer to another by preventing back-diffusion of CO₂ through membranes. In the second case, this isoform would catalyze the final conversion of HCO₃⁻ into CO₂ for photosynthetic needs. It should be noted that, since in N. vectensis at least six different isoforms of carbonic anhydrases have been identified (AM and D. Zoccola, personal communication), other A. viridis CA isoforms are expected to contribute in the Ci transport from seawater to Symbiodinium.

The previous model of Ci transport assumed that CO₂ crosses membranes by diffusion through concentration gradients between both sides of the cell plasma membrane. For most of the past century, gas molecules such as CO₂ were presumed to cross biological membranes merely by diffusing through the lipid phase. This view was challenged recently by studies demonstrating the permeability of water channel aquaporins and certain Rh-family members to CO₂ [62]–[64]. The physiological function of those channels for CO₂ transport seems to be particularly important since, for instance, one Rh protein (RhAG) accounts for up to 50% of the CO₂ transport of human red blood cells [65]. In A. viridis, we identified two RhAG isoforms, one of which (Av_RhAG1) is preferentially expressed in the gastroderm of zooxanthellate anemones (Table 2). We suggest that these proteins could have a role within the holobiont to facilitate CO₂ uptake, possibly in conjunction with the membrane-bound CA2. In human erythrocytes, RhAG and CA2 are part of the same Band 3 multiprotein complex involved in anion exchanges [66], [67]. It is worth mentioning that a role of Av_RhAG1 in NH₄ transport is equally valid, as in the cnidarian-dinoflagellate symbiosis, ammonium resulting from host metabolism is not excreted into the surrounding water but is immediately re-assimilated by the algae, either through diffusion or by a transporter, and then recycled [4], [68].

Is the combination of gene duplication and gastroderm-specific expression a symbiosis adaptation by anthozoans?

Within the set of Kern genes up-regulated in the symbiotic condition, we noticed that some of the genes with a proposed function in symbiosis, including Sym32, Calumenin, and NPC2, were expressed as 2 or more related copies. Phylogenetic analysis allowed the conclusion to be made that there were cnidarian-specific gene duplications. Remarkably, our expression results in A. viridis showed that these cnidarian-specific duplicates were both preferentially expressed in the gastroderm (hosting zooxanthellae) and in the symbiotic condition.

Such cnidarian-specific gene duplications could correlate with the amenability of various cnidarians to have accepted photosynthetic endosymbionts during evolution. However, these gene duplications are not restricted to symbiotic anemones (e.g. N. vectensis). Their selective advantage with regard to symbiosis therefore remains to be determined, together with the origin of endosymbiosis in cnidarians: was photosynthetic symbiosis acquired or lost in various branches of the phylum? Moreover, at the transcriptional level, we have shown for several cases that expression of one isoform among the duplicated gene copies was specifically up-regulated both in the gastroderm and in the presence of symbionts. Thus, we suggest that these neofunctionalizations would be associated with the physiological constraints of endosymbiosis and would be tuned to the presence of zooxanthellae by transcriptional control. It is not known how such control is exerted or whether such gene regulation is restricted to symbiotic cnidarians. Comparative tissue expression of the different orthologs in non-symbiotic cnidarians, such as N. vectensis, would provide insight into the adaptive origin of symbiosis.

Coral bleaching and transcription plasticity

One of our most unexpected results was that outside the subset of Kern genes, many differentially expressed genes were up- or down-regulated in only few individual anemones, despite the fact that anemones were in a fixed symbiotic or aposymbiotic state before sampling. Similar inter-individual variable expression profiles have already been quoted in two studies comparing the response of A. millepora corals to environmental changes [15], [69].

On the other hand, several groups studying symbiosis breakdown in diverse cnidarian species showed that different cellular mechanisms were observed during loss of zooxanthellae, including apoptosis, necrosis, exocytosis, or phagocytosis [11]. Such variability in the cellular processes involved may reflect distinct causes of bleaching, i.e. different responses to variable environmental changes. Indeed, bleaching can be caused by a multitude of environmental stressors including changes in seawater temperature, salinity, ultraviolet radiation, increased sedimentation, nutrients and pollutants [70]. In response to different stresses, the host-associated microbiota is greatly and specifically modified within the whole holobiont (i.e. the community composed of cnidarian host, dinoflagellates and associated microbes) [71].

In the case of the Mediterranean symbiotic cnidarian used here, Symbiodinium belongs to the same clade temperate A. However, our results on the gene copy numbers of only 3 genes (EF2, SPS and APX) showed unexpected high polymorphism in the Symbiodinium hosted by anemones collected from neighboring location. In addition, we showed that the abundance of zooxanthellae per host cell is variable within our set of symbiotic sea anemones. Thus, variation in the associated microbiota or in Symbiodinium sub-clades may have stable extended effects on the host gene expression profile. This underlines that, although the term “bleached” defines one visual phenotype resulting from various responses, it can result from diverse expression profiles, and probably from various conditions of aposymbiosis. Indeed, not all bleached cnidarians die.

A. viridis can be sampled at different time points without harm to the anemone, and is therefore an excellent model organism for kinetics experiments. Additionally, since sea anemones are non-calcifying anthozoans, they allow the study of symbiosis-associated processes outside the cross-regulatory pathways of mineralization found in corals. It would be of interest, using the same array technique, to monitor the kinetics of the specific gene expression response to one or combination of different environmental stressors. In addition to giving a transcriptional map of the gene-specific response to stress, it would highlight whether different bleaching expression profiles are fixed in time or whether they tend to stabilize at a unique Kern profile after an adaptation period. Moreover, the unique amenability of A. viridis to separate the ectoderm from the endoderm will permit us to partition the scope of the environmental response in the cellular layer hosting the symbionts (gastroderm) and in the tissue in direct contact to the milieu (epidermis).

Materials and Methods

Collection and maintenance of A. viridis specimens

Mediterranean sea anemone specimens, Anemonia viridis (Forskål, 1775), were collected in five locations on and around the French Riviera (Figure S2): Antibes (Salis and Croutons sites), Villefranche-sur-Mer, Monaco and Menton. A total of 5 symbiotic and 6 aposymbiotic specimens were used in this study.

Symbiotic specimens Sy1 and Sy3–5 were collected from Antibes Croutons (Figure S2), Sy2 was collected from Antibes Salis. These anemones were maintained for several months in seawater aquaria at 17.0±0.5°C with weekly water renewal. A metal halide lamp (HQI-TS 400 W, Philips) provided light at a constant saturating irradiance of 250 µmol m⁻² s⁻¹ on a 12/12 h light/dark cycle. Naturally-occurring aposymbiotic animals were sampled from the public aquaria of the Oceanographic Museum of Monaco and originated from different sites in Monaco (AS1 and AS2) or Menton (AS3). The three other aposymbiotic specimens were obtained after either a thermal stress-induced bleaching (AS4, AS5, collected in Antibes-Crouton) or a treatment with the catalase inhibitor aminotriazol (AS6, collected in Villefranche-sur-Mer, [7]). These latter 3 anemones were maintained in the dark (necessary to keep them bleached) for a minimum of 2 months. Both symbiotic and stress induced aposymbiotic anemones were fed twice a week with frozen Artemia salina.

Sy3, Sy4 and Sy5 specimens were later subjected to an 8°C-temperature increase (17 to 25°C) to assess imposed stress response. Tentacles were sampled over a 48 h kinetic period (t0, t24h and t48h) and then subjected to RT-qPCR experiments.

Specimen sampling was always at 10:00 (2 hours after light start) to avoid circadian effect of gene expression fluctuations.

Nucleic acid extraction

RNA extraction

Total RNA was extracted from whole tentacle samples (epidermis plus gastroderm and zooxanthellae) or separated tissues using Trizol Reagent (Invitrogen), as described previously [20]. Gastrodermal (gastroderm and zooxanthellae) and epidermis layers were separated as described previously [9]. A DNase treatment (RQ1-DNAse, Promega) was performed on the RNA sample to avoid genomic DNA contamination.

gDNA extraction

Genomic DNA was extracted as follows: samples (whole tentacles, separated tissue layers (epidermis or gastroderm) or cultured zooxanthellae) were first immersed in [100 mM Tris pH = 8.2; 100 mM EDTA; 500 mM NaCl; 2.5% triton] to block DNases, and then transferred into [4 M Thiocyanante Guanidium; 25 mM Sodium citrate; 1% Triton X100; 50 µg/ml Proteinase K; 5 mM EDTA; 50 mM DTT] before potter homogenization. Homogenates were incubated at 55°C for 1 hour. Genomic DNA was purified by phenol/chloroform/isoamyl alcohol extraction and isopropanol precipitation.

A. viridis oligoarray hybridizations

Over 2,000 genes were selected from the A. viridis clustered and annotated EST collection [20] for their putative participation in symbiotic processes. Genes were selected by keyword search in their functional annotation for matches with the following terms: intracellular transport, metabolic processes (lipids, proteins), signal transduction, organelle organisation and biogenesis, response to stress, trans-membrane, apoptosis and cell death. Two thousand 60-mer oligonucleotides were designed (Eurogentec, Belgium) and printed on slides in triplicate (Eurogentec, Belgium). A luciferase oligonucleotide was also spotted as an external control. This A. viridis Oligo2K version 1.0 oligoarray is fully described under platform record GPL10546, stored in the Gene Expression Omnibus (GEO) at NCBI (http://www.ncbi.nlm.nih.gov/geo).

cDNA synthesis and labelling was performed from total RNA using the ChipShot Direct labelling and Clean-up system (Promega), according to the manufacturer's instructions. Five ng of a Lux mRNA exogenous standard was added to each mRNA sample before labelling. RNA quality was evaluated using the Agilent Bioanalyzer 2100 and quantified on a ND-1000 Spectrophotometer (NanoDrop). Then, 500 ng of each labeled sample were mixed in Hi-RPM hybridization buffer (Agilent), and incubated on the array overnight at 47°C. Slides were scanned after post-hybridization washes, using a GenePix 4200A scanner (Axon Instruments). Data acquisition and quality control were performed using Genepix Pro software.

Experimental design and statistical analyses

Dye-swap experiments were performed for each hybridization condition. Two different sets of hybridizations were performed (Figure S3): i) symbiotic versus aposymbiotic specimens and ii) epidermis versus gastroderm tissues. The same AS6 extract was used as the sample reference in all experiments. Experimental data and associated microarray designs were deposited in the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo/) under SuperSeries record: GSE22375 and platform record GPL10546.

Background (positive offset of 50) was evaluated according as described in [72]. The data were normalized by the print-tip loess method (within-array normalization) and by quantile method (between-array normalization) using the LimmaGUI package from Bioconductor [23]. Means of ratios from all comparisons were calculated for each gene, and B test analysis was done using LimmaGUI. Differential gene expression was determined by the Bayesian statistical method (B value), with a cut-off of zero for the B value and a log ratio |M|>0.590 as significant. All normalized data sets were registered in the GEO database under the accession number GSE22375.

Cluster 3.0 software [73] was used to estimate the hierarchical clustering between individual anemone array results. The following parameters (complete, average, and centroid linkage) were tested and gave similar results.

Real-time quantitative PCR experiments

Specific primers amplifying around 100 bp were designed using the software Primer3 [74]. The primer sequences used in this study are listed in Table S3. Amplicon specificity for either A. viridis or Symbiodinium was tested against epidermal (extracts without zooxanthellae) or CZ (zooxanthellae culture) extracts, respectively. Expected length of the amplicons was checked by agarose gel electrophoresis after regular PCR amplification. Primer efficiencies were determined using standard curve analysis with a 10-fold dilution series of pooled cDNA from both control and treated samples (data not shown), and ranged from 1.8 to 2. The qPCR products were sequenced (Macrogen Inc, Korea) and all matched the expected product identities.

cDNAs were prepared using SuperScriptII reverse transcriptase (Invitrogen) and a mixture of oligodT and random primers, according to the manufacturer's instructions. Transcript level quantification was performed using the SYBR green fluorescence method and a Light Cycler 480 (Roche). The PCR conditions were as follows: 1× SYBR green mix (LC480 SYBR Green Master Mix, Roche), 100 nM primers and 2.5 ng of cDNA in a total volume of 15 µl. Each sample was run in triplicate using the following PCR parameters: 94°C for 10 min, followed by 40 cycles of 15 s at 94°C, 20 s at 60°C and 15 s at 72°C, then a dissociation curve step (60 to 95°C) to confirm the absence of non-specific products. The dissociation curves showed a single amplification product and no primer dimers.

Several control genes were chosen based on the microarray results as a whole (most stably expressed genes, in all tested conditions), and the expression stability of ten putative control genes was evaluated using the GeNorm software [75]. A reliable normalization factor was calculated based on the expression level of the most stable control genes. The control genes finally selected in this study are RPLP0, RCC2, and COP-γ. Expression levels of target genes were normalized using the normalization factor described above and the results given as expression relative to the aposymbiotic specimen (AS6) value as calibrator (reference). The significance of the results was tested using t-tests (software Jump 5.1, Cary, USA). Results were considered statistically significant when P<0.05.

Bleaching quantification

The relative abundance of zooxanthellae in A. viridis cells was quantified using a real-time quantitative qPCR method. For each organism, three nuclear genes were chosen from available datasets: COP-γ (Coatomer subunit gamma), RCC2 (Regulator of Chromosome Condensation protein 2) and NPC1 (Niemann Pick type C1) for A. viridis, and EF2 (Elongation factor 2), Sucrose phosphate synthase (SPS) and Ascorbate peroxidase (APX) for Symbiodinium clade temperate A. Primer design was performed using Primer3 and the same parameters described above; the primer sequences are given in the Table S3. Standard curves were generated with six logarithm dilutions of corresponding cloned sequences. Results are expressed as relative quantification of Symbiodinium nuclei (or nuclear genes) to A.viridis nuclei.

Sequence analyses

Sequencing of either A.viridis library clones or newly cloned cDNAs was performed by Macrogen Inc.

Signal peptides and Trans-membrane domains were predicted using SignalP and TMHMM, respectively, from the CBS prediction servers (www.cbs.dtu.dk/services). GPI anchors were predicted with PredGPI (http://gpcr.biocomp.unibo.it/predgpi). Other Domains were inferred from PFAM (pfam.sanger.ac.uk). Blast analysis tools used were from NCBI (blast.ncbi.nlm.nih.gov).

Functionnal annotation of the A. viridis EST dataset (GO annotations) has been performed using Blast2GO [20]. Statistical assessment of annotation differences between SY-genes (test group) and APO-genes (reference group) was performed using the Gossip package [76], that employs Fisher's Exact Test to estimate the significance of associations between two categorical variables.

Sequence alignments were performed using MultAlin [77] and ClustalW with Blosum62 default parameters. Alignments were optimized manually using the 2 computer generated alignments as a model. Using the segment of the sequence alignment conserved in all sequences (bordered by the bar above alignments in Figures S5, S6, S7), the best-fitted substitution model was evaluated using ProtTest [78]. Using parameters indicated in the Figure legend of each alignment tree, a Maximum Likelihood tree was determined using phyML and branches support were calculated using aLRT [79]. Alternatively, Bayesian analysis using MrBayes 3-1.2 (mrbayes.csit.fsu.edu) was conducted with the following settings: fixed rate amino acid model was set to mixed (prset aamodelpr = mixed) and proportion of invariable sites model was combined with the Gamma model (lset rates = invgamma).

Supporting Information

Attachment 1

Attachment 2

Attachment 3

Attachment 4

Attachment 5

Attachment 6

Attachment 7

Attachment 8

Attachment 9

Attachment 10

Attachment 11