Differential expression of Exaiptasia pallida GIMAP genes upon induction of apoptosis and autophagy suggests a potential role in cnidarian symbiosis and disease

Coral reefs provide a variety of important services including protection of shorelines from erosion, habitat for numerous organisms, and promoting fishing and tourism industries (Kirk and Weis, 2016; Moberg and Folke, 1999). Underlying the success of this productive and diverse ecosystem is the mutualistic relationship between corals and photosynthetic dinoflagellates of the family Symbiodiniaceae (Kirk and Weis, 2016; Roth, 2014). The cnidarian host relies on the photosynthetic capability of the dinoflagellate as a carbon source in nutrient-poor water, and, in return, the dinoflagellate receives a stable intracellular environment and inorganic nutrients such as nitrogen and phosphorus (Davy et al., 2012; Yellowlees et al., 2008). Currently, there are several threats to the vitality of coral reefs, including increased ocean temperatures, ocean acidification and coral disease, which ultimately lead to the destruction of this important ecosystem (Hoegh-Guldberg et al., 2007; Hughes et al., 2018; Randall and van Woesik, 2015). Many of these stressors can also cause massive coral bleaching, or the breakdown of symbiosis, which can not only lead to coral death but also negatively impacts organisms and industries that rely on the reef (Hoegh-Guldberg et al., 2017; Hughes et al., 2010). Coral disease is also a major threat to the reef ecosystem, with novel disease outbreaks occurring as a result of increased water temperatures (Randall and van Woesik, 2015). The molecular mechanisms involved in bleaching and disease are still not well understood. While it is known that some diseases cause bleaching or specifically target the dinoflagellate symbionts (Cervino et al., 2004; Rosenberg et al., 2009), less work has been done on how the symbiotic state of cnidarians impacts the cellular and molecular response to an invading microbe. As both bleaching and disease can simultaneously impact a reef, it is critical that more research be done to understand the common pathways involved in these responses.

A potential regulator of both cnidarian symbiosis and disease is GTPases of immunity associated proteins, or GIMAPs. GIMAPs represent a family of small G proteins that are most closely related to septin and dynamin GTPases and are characterized by the presence of the GTP binding AIG1 domain (Nitta and Takahama, 2007; Schwefel et al., 2010). These proteins are best characterized in vertebrates, where they show the highest expression in immune cells (Nitta et al., 2006). Specifically, GIMAPs regulate apoptosis, or programmed cell death, to determine lymphocyte fate during maturation and homeostasis, which has been proposed to occur through their interaction with Bcl-2 family proteins (Nitta et al., 2006; Wang and Li, 2009). During the initiation of apoptosis, the altered activity of Bcl-2 family proteins increases the permeability of the outer mitochondrial membrane, leading to the release of cytochrome c and ultimately activation of proteases called caspases that cleave key proteins, cytoskeletal elements and cell adhesion molecules (Czabotar et al., 2014; Moya et al., 2016). This in turn creates orderly death of the cell. In addition to regulating apoptosis, vertebrate GIMAPs have also been implicated in autophagy. Specifically, it has been shown that human GIMAP6 interacts with the autophagic protein GABARAPL2 (a mammalian Atg8 homolog), localizes to autophagosomes under starvation conditions, and is required for efficient autophagy (Pascall et al., 2013, 2018).

The mechanisms by which vertebrate GIMAPs operate are still poorly understood. Interestingly, it has been shown that not all GIMAPs have the catalytic ability to hydrolyze GTP, and that GTPase activity is regulated by the formation of homodimers and heterodimers (Schwefel et al., 2013). Specifically, it has been proposed that GIMAPs with GTPase activity form heterodimers with non-catalytic GIMAPs to stimulate GTP hydrolysis (Schwefel et al., 2013). In regards to apoptosis, it has been hypothesized that GIMAPs lacking GTPase activity are anti-apoptotic in function through the formation of GIMAP scaffolds, while those that possess catalytic ability promote apoptosis by disrupting these scaffolds (Ciucci and Bosselut, 2014; Schwefel et al., 2013). However, more work is needed to fully explore these hypotheses to understand the mechanisms by which GIMAPs regulate cellular processes.

While previously thought to be limited to vertebrates and plants, GIMAPs have recently been characterized in invertebrates including cnidarians and mollusks (Lu et al., 2020; McDowell et al., 2016; Weiss et al., 2013). In both phyla, studies indicated a conserved role for invertebrate GIMAPs in innate immunity. Specifically, when the coral Acropora millepora was exposed to muramyl dipeptide (MDP), a compound found in the bacterial cell wall, three GIMAP transcripts showed a massive increase in expression, indicating that these genes may play a role in cnidarian immunity (Weiss et al., 2013). Studies in two mollusk species, the snail Biomphalaria glabrata and the bivalve Crassostrea virginica, demonstrate upregulation of GIMAP genes in response to individual immune stimulants or exposure to live bacteria, respectively, also demonstrating a role in molluscan immunity (McDowell et al., 2014; Zhang et al., 2016). Therefore, investigating the role of invertebrate GIMAPs in immunity is an area that warrants further investigation.

In addition to looking at the role of cnidarian GIMAPs in immunity, we were also interested in assessing their conserved role in the processes of apoptosis and autophagy. These cellular pathways were of interest because of the known role that they play in both cnidarian symbiosis and disease. Apoptosis has previously been shown to be involved in both the onset and breakdown of symbiosis as it is activated upon colonization with heterologous symbionts and also promotes bleaching (Dunn and Weis, 2009; Dunn et al., 2007; Pernice et al., 2011; Tchernov et al., 2011). Therefore, the inhibition of apoptosis appears to be important for maintaining a stable symbiosis. Additionally, genes involved in apoptosis including pro-apoptotic Bcl-2 family members and caspases have been shown to be differentially expressed during several coral diseases or upon immune stimulation (Ainsworth et al., 2015; Fuess et al., 2017; Libro et al., 2013). There is also evidence to suggest that autophagy plays an active role in bleaching. Specifically, chemical induction of autophagy results in bleaching and has been shown to promote the formation of autophagic bodies within host tissue (Dunn et al., 2007; Hanes and Kempf, 2013). Autophagy has also been shown to be activated in more disease-resistant corals, therefore representing an important immune defence (Fuess et al., 2017). Together, these studies indicate that apoptosis and autophagy are involved in regulating cnidarian interactions with beneficial and pathogenic microbes in their environment. Thus, if GIMAPs regulate apoptosis and autophagy, they are likely important regulators of the bleaching and disease response, and therefore it is important to determine their function in cnidarians.

To investigate the role of cnidarian GIMAPs in apoptosis, autophagy and immunity, the emerging model anemone Exaiptasia pallida (Agassiz in Verrill 1864) (Grajales and Rodríguez, 2016) was used in this study. E. pallida is commonly used as a cnidarian model to study corals because it is relatively easy to rear and manipulate in the lab setting (Weis et al., 2008). In addition, E. pallida engages in a symbiotic relationship with Symbiodiniaceae, and can be maintained in both symbiotic and aposymbiotic states, allowing for studies on the dynamics of symbiosis (Weis et al., 2008). To better understand the function of cnidarian GIMAPs, the E. pallida genome was searched for GIMAP-like sequences, and subsequently expression of the four recovered genes, termed Ep_GIMAPs, was measured in response to chemical induction of apoptosis, autophagy and the immune response. This study represents the first targeted investigation of cnidarian GIMAPs and therefore provides important information regarding the function of these ancestral proteins.

For all experiments in this study, E. pallida of the lab genotype CC7 were used. All E. pallida were maintained at room temperature in Instant Ocean^® with a salinity of 35 ppt and were fed Artemia (brine shrimp) 3 times a week. Symbiotic anemones were kept under a 12 h light:dark cycle with a light intensity of approximately 25 μmol m⁻² s⁻¹. Aposymbiotic CC7 E. pallida were obtained using a menthol bleaching protocol modified slightly from previous studies (Bailey et al., 2019; Matthews et al., 2016). Symbiotic anemones were treated with a final concentration of 0.0913 mmol l⁻¹ menthol (dissolved in 100% ethanol) in 0.45 μm filtered sea water (FSW) for periods of 24 h until they were shown to be symbiont free based on imaging with a fluorescence microscope. Upon starting the menthol treatments, anemones were kept in the dark in FSW to prevent recolonization with symbionts. Before use in experiments, anemones were allowed to recover for a minimum of 3 months. Three to five days prior to the start of an experiment, anemones were moved into 6-well plates (1 animal per well) and were not fed to prevent Artemia contamination.

The E. pallida genome was searched for GIMAP-like sequences using the Reef Genomics website (http://reefgenomics.org; Baumgarten et al., 2015; Liew et al., 2016). Previously characterized A. millepora GIMAP sequences (Weiss et al., 2013), as well as human and mouse GIMAP sequences obtained from NCBI (https://ncbi.nlm.nih.gov) (Table S2), were used as queries for a BLASTp search against the E. pallida genome using an e-value cutoff of 1×10⁻². A high e-value was used to account for the high sequence divergence between cnidarian and vertebrate GIMAP sequences. All obtained sequences were then entered into the NCBI conserved domain database (CDD) search tool (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi; Marchler-Bauer et al., 2017) to confirm the presence of the AIG1 domain (pfam04548 or cd01852) that is characteristic of GIMAPs, with an e-value cutoff of 1×10⁻⁴. Subsequently, all sequences containing the AIG1 domain were used as a query for a reciprocal BLASTp search in NCBI against Homo sapiens (taxid:9606), to ensure the top hit was a human GIMAP, further confirming the sequence identity. Lastly, as vertebrate GIMAPs also contain a coiled coil domain, each sequence was run through the program COILS (https://embnet.vital-it.ch/software/COILS_form.html; Lupas et al., 1991) using the unweighted MTIKD matrix to check for this structure. For all recovered proteins, the corresponding nucleotide sequence was also recovered from Reef Genomics. The BLAST searches revealed four GIMAP-like sequences in E. pallida: AIPGENE4644, AIPGENE4645, AIPGENE4712 and AIPGENE4714, which have been renamed Ep_GIMAP1, Ep_GIMAP2, Ep_GIMAP3 and Ep_GIMAP4, respectively.

Primers for each Ep_GIMAP were designed using primer3 either through the website version (https://www.primer3plus.com; Untergasser et al., 2012) or through the Geneious version 11.1.3 plugin (https://www.geneious.com; Kearse et al., 2012) to confirm the genome sequence (Table 1). PCRs for each primer set were performed using the GoTaq Flexi kit (Promega, Madison, WI, USA) with the following run protocol: 94°C for 2 min, 35 cycles of 94°C for 45 s, annealing temperature for 45 s (Table 1), and 72°C for 1 min, followed by a final extension at 72°C for 10 min. The resulting PCR products were verified on a 1% agarose gel and cleaned using the QIAquick PCR purification kit (Qiagen, Germantown, MD, USA), following the manufacturer's instructions, with a final elution step using 30 µl of water. Purified PCR products were sent to GenScript for sequencing (GenScript, Piscataway, NJ, USA). Any sequences that differed significantly from the genome sequences present in either Reef Genomics or NCBI were submitted to GenBank.

Wortmannin is a phosphoinositide 3-kinase (PI3K) inhibitor that has previously been used in E. pallida to induce apoptosis, while simultaneously inhibiting autophagy (Dunn et al., 2007). Therefore, this chemical was used for the purpose of inducing apoptosis and examining corresponding changes in Ep_GIMAP expression to understand their potential role in this process. Four treatment groups were used (n=3 per treatment group) including a 0.5 μmol l⁻¹ wortmannin treatment (product W3144, Sigma-Aldrich, St Louis MO, USA) and a vehicle control, which consisted of an equivalent amount of DMSO, for both symbiotic and aposymbiotic animals. Both symbiotic states were used to determine how the response to apoptosis induction was influenced by the presence of symbionts. For all treatments, the appropriate volume of either wortmannin or DMSO was added to FSW for a final treatment volume of 10 ml. Animals were sampled at 12, 24 and 48 h post-treatment for a caspase assay and RNA extractions. Anemones designated for the caspase assay were immediately processed according to the kit protocol (see details below) and those that were designated for RNA extractions were placed in 500 µl of RNAlater (Invitrogen, Carlsbad CA, USA) and stored at 4°C until further processing.

Immediately after sampling from the wortmannin treatment experiment, anemones were processed for a caspase activity assay to determine the level of apoptosis occurring at each time point. The caspase activity assay was completed with the ApoAlert Caspase-3 Colorimetric Kit (Takara, Mountain View, CA, USA). Briefly, anemones were homogenized in 300 of  µl cell lysis buffer, incubated on ice for 10 min, and centrifuged at 18,000 g for 10 min. The supernatant was then transferred to a new tube and stored at −20°C until further use. Subsequently, the manufacturer's protocol was followed to determine caspase 3 activity for each treatment–time point combination using a Synergy LX Multi-Mode Reader (BioTek, Winooski, VT, USA). Additionally, protein quantification using the Pierce BCA Protein assay kit (ThermoFisher, Waltham, MA, USA) was performed to account for the size of the animal. For both assays, each sample was run in duplicate. Caspase activity was normalized to the size of the animal by dividing the caspase activity by the amount of protein in the same volume of sample. Subsequently, the fold-change in caspase activity was calculated by dividing the size normalized caspase activity by the average control for the respective time point and symbiotic state. The fold-changes were analyzed using a 3-way ANOVA (factors: treatment, time and symbiotic state, with the model Caspase activity∼treatment×time×symbiotic state) followed by a Tukey's post hoc test for pairwise comparisons to determine whether there was a significant difference in caspase activity between treatments, across time and between symbiotic states. All statistical analyses and graph construction were completed through R version 3.6.1 (http://www.R-project.org/).

Rapamycin is an mTOR inhibitor that is commonly used to induce autophagy. This chemical has previously been used for this purpose in E. pallida and was shown to induce the formation of autophagic bodies at 12 h post-treatment (Hanes and Kempf, 2013). Therefore, we chose to use the same concentration and time point in this study. The four treatment groups included a 25 µmol l⁻¹ rapamycin treatment (product 13346, Cayman Chemicals, Ann Arbor, MI, USA) and a vehicle control with an equivalent amount of DMSO in both symbiotic and aposymbiotic animals (n=3 per treatment). These treatments were designated AC for the aposymbiotic vehicle control, AR for the aposymbiotic rapamycin, SC for the symbiotic vehicle control and SR for the symbiotic rapamycin. For all treatments, the appropriate amount of either rapamycin or DMSO was added to FSW for a final treatment volume of 10 ml. At 12 h post-treatment, the anemones were placed in 500 µl RNAlater (Invitrogen) for RNA extractions and stored at 4°C until processing.

Lipopolysaccharide (LPS), a component of the bacterial cell wall of gram-negative bacteria, was used to investigate the role of Ep_GIMAPs in the immune response. To decide on the appropriate time frame for sampling, a pilot study was conducted in symbiotic E. pallida in which anemones were exposed to LPS (L2880, Sigma-Aldrich) at a concentration of 10 µg ml⁻¹, dissolved in PBS for a total treatment volume of 10 ml. Vehicle controls were exposed to an equivalent volume of PBS in FSW. At 3 h post-treatment, anemones were placed in 500 µl of RNAlater (Invitrogen) for RNA extractions and stored at 4°C until processing. Based on the fact that no significant changes were observed at this time point (Fig. S1), a larger experiment with sampling at later time points was performed. For this experiment, both aposymbiotic and symbiotic E. pallida were exposed to either the LPS or vehicle control treatments as previously described for the pilot experiment. Treatments were refreshed every 24 h, and animals were sampled at 24, 48 and 72 h post-treatment (n=3 for each time point/treatment/symbiotic state combination).

For all experiments, RNA extractions were performed using a hybrid TRIzol reagent (Invitrogen) and RNeasy mini kit (Qiagen, Valencia, CA, USA) protocol as previously described (Bailey et al., 2019; Poole et al., 2016). The quality and quantity of the RNA obtained were determined using a Nanodrop-2000 (Thermo Fisher, Waltham, MA, USA) and by running samples on a 1% agarose gel. For samples with a Nanodrop A_260/230 value below 1.5, an RNA cleanup protocol was performed as previously described (Poole et al., 2016). RNA samples were stored at −80°C until future use. Following RNA extraction, cDNA synthesis was performed using the iScript gDNA Clear cDNA Synthesis Kit (Bio-Rad, Hercules, CA, USA) starting with 500 ng of total RNA. Success of cDNA synthesis was checked using L10 qPCR primers (Poole et al., 2016), and then samples were diluted with 13.3 μl of RNase-free water and stored at −20°C until future use in qPCR.

qPCR primers for all for Ep_GIMAPs were designed using the built-in primer tool in Geneious version 11.1.3 (https://www.geneious.com; Kearse et al., 2012) to amplify 100–200 bp regions of each transcript (Table 1). All products were cleaned and subsequently sequenced as previously described to ensure the correct product was amplified. A StepOne machine (Applied Biosystems, Foster City, CA, USA) was used to test primer efficiency, using cDNA dilutions of 1:4, 1:16, 1:64 and 1:256 (run protocol and reaction components described below). To calculate primer efficiency, dilutions were plotted against the Ct values, and the equation efficiency =10^−1/slope was used based on the best-fit line obtained for the data. Several primer sets were tested for each Ep_GIMAP, and the one with an efficiency closest to 100% (within a range of 90–110%) was selected for use in qPCR.

For all qPCR experiments performed, L10, L12 and PABP were used as reference genes. These reference genes were previously demonstrated to exhibit stable expression in studies related to symbiosis, which was relevant to all experiments performed in this study (Poole et al., 2016). Furthermore, the stability of these reference genes was verified for the wortmannin treatment by performing qPCR on the SC treatment at 12 h and AW treatment at 24 h. L10, L12, PABP and EF_1α were tested as candidate reference genes and geNorm (Vandesompele et al., 2002) analysis revealed that L10, L12 and PABP were the most stable.

All experimental qPCR plates were run on a QuantStudio7 Flex Real-Time PCR machine (Applied Biosystems, Foster City, CA, USA) under the comparative Ct (ΔΔCt) setting using the default run protocol with the addition of a melt curve. On each qPCR plate, the genes assessed included the three reference genes and four Ep_GIMAPs. Additionally, for the rapamycin experiment, LC3, a homolog of yeast Atg8, was included as an indicator of levels of autophagy. LC3 is recruited to autophagosomes during autophagy and initiates elongation (Weidberg et al., 2010). Although traditionally LC3 degradation rather than mRNA levels is used to assess levels of autophagy (Yoshii and Mizushima, 2017), this was not feasible without a custom E. pallida antibody. Therefore, although gene expression levels do not necessarily correlate with LC3 activity or turnover, this was the best available method for quantifying autophagy. Samples were run in triplicate and each plate included no-template and no-primer controls. A total reaction volume of 20 μl was used, which included 7.5 μl of water, 1 μl of each forward and reverse primer (5 μmol concentration), 0.5 μl of cDNA template and 10 μl of Power SYBR Green PCR Master Mix (Applied Biosystems). Resulting Ct values were exported from the software and relative gene expression levels were calculated as previously described (Bailey et al., 2019). Briefly, expression of the genes of interest was normalized to the geometric mean of the reference genes (L10, L12 and PABP) to produce a ΔCt value. Subsequently, the ΔΔCt was calculated by subtracting the ΔCt for each sample from the average of the reference sample ΔCt values. The reference sample was always the vehicle control for a particular time point, treatment and symbiotic state combination. The ΔΔCt, or relative quantities, were analyzed using an ANOVA for each gene (three-way ANOVA with factors of time, treatment and symbiotic state for wortmannin and LPS experiments with the model Relative quantity∼treatment*time*symbiotic state and two-way ANOVA with factors of treatment and symbiotic state for rapamycin experiment with the model Relative quantity∼treatment*symbiotic state) with Tukey's post hoc test for pairwise comparisons to determine whether there was a significant difference in gene expression between the vehicle controls and treated samples for a particular time point. Additionally, Pearson's correlation tests were performed to check for significant relationships between caspase activity and Ep_GIMAP expression for the wortmannin treatment and LC3 and Ep_GIMAP expression for the rapamycin treatment. Lastly, for the 3 h LPS treatment of symbiotic animals, a Wilcoxon rank sum test was performed. All statistical analyses and graph construction were completed through R version 3.6.1 (http://www.R-project.org/).

To identify potential interaction partners for the Ep_GIMAPs, the E. pallida genome was searched for Bcl-2 family proteins using the Reef Genomics website (http://reefgenomics.org; Baumgarten et al., 2015; Liew et al., 2016). To do so, protein sequences for Bcl-2 family members that were previously characterized in A. millepora (Moya et al., 2016) were used as queries in a BLASTp search against the E. pallida genome using an e-value cutoff of 1×10⁻²⁰, as it was at this point that a notable separation between BLAST results was observed, with the next closest hits having largely insignificant e-values. A much lower e-value was used than for the GIMAPs as the searches were performed exclusively with cnidarian sequences and therefore lower divergence between the query sequences and the hits was expected. Each E. pallida ortholog obtained underwent a domain search using the NCBI conserved domain database CDD search tool to ensure domain structure matched that of the corresponding A. millepora protein sequence (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi; Marchler-Bauer et al., 2017). Sequences were aligned in Geneious version 11.1.3 to remove any identical sequences (https://www.geneious.com; Kearse et al., 2012). All obtained sequences were then used as a query for a reciprocal BLASTp search in NCBI against Homo sapiens (taxid:9606), to ensure identity as a Bcl-2 family member.

Four GIMAP-like genes were identified in the E. pallida genome: AIPGENE4644, AIPGENE4645, AIPGENE4712 and AIPGENE4714, which will be referred to hereafter as Ep_GIMAPs. The Ep_GIMAP genes were found in clusters in the genome where Ep_GIMAP1 and Ep_GIMAP2 are located tandem to one another on scaffold 1042, whereas Ep_GIMAP3 and Ep_GIMAP4 are both located on scaffold 56. Based on the reciprocal BLAST results, the top hit for Ep_GIMAP1 was human GIMAP6 (e-value=0.001) and the top hit for Ep_GIMAP3 was human GIMAP7 (e-value=1×10⁻⁷). Ep_GIMAP2 and Ep_GIMAP4 had no significant hits to the human proteins in the reciprocal BLAST, but were still included because of their similarity to the previously characterized A. millepora GIMAPs (e-values around 1×10⁻³⁰). All Ep_GIMAPs contained the AIG1 domain that is characteristic of vertebrate GIMAP proteins. When run through the program COILS, all Ep_GIMAPs with the exception of Ep_GIMAP4 returned a probability of having a region with a coiled coil above 0.5 (Table S1).

Through verifying each genome sequence acquired from Reef Genomics through PCR and sequencing, discrepancies were discovered for Ep_GIMAP1, Ep_GIMAP2 and Ep_GIMAP4. The original sequence for AIPGENE4644 (Ep_GIMAP1) matched two protein sequences in the version of the genome on NCBI, thus indicating this predicted sequence may contain two genes. Our sequencing data confirmed a start codon in the middle of this sequence, thereby providing evidence for the NCBI version for this gene (accession number XM_021039897.2).

There were also difficulties encountered in amplifying the Ep_GIMAP2 sequence when using primers designed to the available AIPGENE4645 on Reef Genomics. As such, a reassembled version of the raw Illumina sequence reads for accession number SRR696721 that were downloaded from the sequence read archive (SRA) entry for the aposymbiotic CC7 transcriptome (http://www.ncbi.nlm.nih.gov/sra/SRX231866) using Trinity (Grabherr et al., 2011) was searched for a corresponding transcript, and a discrepancy in the 5′ region was present between the genome and reassembled transcriptome sequence. Primers were designed to the transcriptome sequence and the product amplified through PCR confirmed the transcriptome version of the 5′ region of the transcript and has been submitted to GenBank (accession number MT495601). Lastly, the AIPGENE4714 (Ep_GIMAP4) sequence on Reef Genomics corresponded to two different proteins and contained a mannosyltransferase domain at the N-terminus and AIG1 domain at the C-terminus, confirmed using the NCBI CDD search tool (Marchler-Bauer et al., 2017). When compared with the reassembled transcriptome sequence, two distinct transcripts were identified, matching the two halves of AIPGENE4714. The sequencing data confirmed the transcriptome sequence for Ep_GIMAP4, and has been submitted to GenBank (accession number MT495600).

The results from the caspase-3 activity assay indicated that wortmannin does induce apoptosis in E. pallida, but that this response does not occur immediately. Specifically, at 12 h post-treatment, caspase activity between the vehicle control and wortmannin treatment was not significantly different in either symbiotic or aposymbiotic E. pallida (Fig. 1). In contrast, at 24 h post-treatment, the caspase activity of the wortmannin-treated animals was significantly higher than that of the vehicle control in aposymbiotic and symbiotic anemones (P<0.001, P<0.017; ANOVA, Tukey's HSD), indicating elevated levels of apoptosis (Fig. 1). Similarly, at 48 h post-treatment, caspase activity was significantly higher in the aposymbiotic wortmannin-treated animals compared to the vehicle control (P<0.001, ANOVA, Tukey's HSD). The symbiotic anemones did not have a significant increase in caspase activity at 48 h as a result of wortmannin treatment (Fig. 1).

Furthermore, these data demonstrated a clear difference in the response between symbiotic and aposymbiotic animals. Over the course of the 48 h experiment, wortmannin-treated aposymbiotic animals showed a trend of increased caspase-3 activity over time, with the 24 and 48 h time points having significantly higher levels of caspase-3 activity than at 12 h post-treatment (P<0.001, ANOVA, Tukey's HSD). For symbiotic anemones, wortmannin-treated animals at 24 h had significantly increased caspase activity compared with the 12 h time point (P=0.021, ANOVA, Tukey's HSD) and there was a significant decrease in caspase activity between 24 and 48 h (P=0.006, ANOVA, Tukey's HSD). Lastly, there was a significant difference between caspase-3 activity of wortmannin-treated symbiotic and aposymbiotic animals at the 48 h time point (P<0.001, ANOVA, Tukey's HSD). Specifically, the aposymbiotic animals had approximately 5 times greater caspase activity than the symbiotic animals (Fig. 1). Overall, these data provide information on the time points at which gene expression differences could be expected if Ep_GIMAPs are involved in apoptosis.

The qPCR results for the wortmannin treatment indicated that when apoptosis is elevated, there is either a trend of or significant downregulation of most Ep_GIMAPs compared with the vehicle control. Specifically, at 12 h post-treatment, when there were not elevated levels of apoptosis, there were no significant differences in gene expression in either symbiotic or aposymbiotic E. pallida compared with the vehicle control (Fig. 2). At 24 h post-treatment, when elevated levels of apoptosis were detected, there were several significant changes in gene expression. Specifically, Ep_GIMAP1 was significantly downregulated in aposymbiotic and symbiotic E. pallida, as was Ep_GIMAP4 in symbiotic E. pallida compared with the vehicle control (Fig. 2A,D) (P=0.006, P=0.004, P=0.016, ANOVA, Tukey's HSD). While all other changes were not significant, there was a downward trend for all Ep_GIMAPs in both symbiotic states. At 48 h post-treatment, there were also several significant differences in gene expression. Specifically, Ep_GIMAP1, Ep_GIMAP2 and Ep_GIMAP4 were significantly downregulated in the aposymbiotic E. pallida compared to the vehicle control (Fig. 2A,B,D) (P<0.001, P<0.001, P<0.003, ANOVA, Tukey's HSD). However, at this time point, no significant differences in expression were detected in symbiotic animals.

While gene expression levels between symbiotic states were very similar at 12 and 24 h post-treatment, there were significant differences between wortmannin-treated aposymbiotic and symbiotic animals at 48 h post-treatment. Specifically, the aposymbiotic animals showed a much stronger response and had significantly lower levels of gene expression for both Ep_GIMAP1 and Ep_GIMAP2 (Fig. 2A,B) (P<0.001 for both, ANOVA, Tukey's HSD).

As there is no commercially available kit to detect cnidarian autophagy, knowledge from previous experiments (Hanes and Kempf, 2013) in combination with our LC3 gene expression data provide evidence for this process occurring at 12 h post-treatment with rapamycin. LC3 is an autophagic protein and previous work demonstrates that E. pallida LC3 (called Ep_LC3) localizes to autophagosomes in a mouse cell line (Flesher, 2013), consistent with patterns observed for vertebrate LC3. Given this functional evidence, it was decided LC3 was the best marker for autophagy in E. pallida and previously developed qPCR primers were used for this gene (Flesher, 2013). The qPCR results for rapamycin treatment indicated that Ep_LC3 was upregulated between the SC and SR groups (Fig. 3) (P=0.005, ANOVA, Tukey's HSD). There was still a trend of increased Ep_LC3 expression in the AR animals, but this was not significantly different from the AC group. The qPCR results also indicated that rapamycin treatment resulted in a downregulation of some Ep_GIMAPs. Specifically, Ep_GIMAP1 and Ep_GIMAP2 showed a significant downregulation between the AC and AR groups as well as the SC and SR groups (Fig. 3) (P=0.003, P<0.001 for aposymbiotic, P<0.001 for symbiotic, ANOVA, Tukey's HSD). Furthermore, there was a stronger downregulation in the symbiotic than in the aposymbiotic animals, which was significant for both Ep_GIMAP1 and Ep_GIMAP2 (P<0.001, P<0.020, ANOVA, Tukey's HSD). In contrast, there were no significant differences in expression for Ep_GIMAP3 and Ep_GIMAP4 in response to rapamycin treatment.

There were no significant changes in Ep_GIMAP gene expression observed in either the pilot study (Fig. S1) or symbiotic and aposymbiotic LPS-treated E. pallida compared with the vehicle control at 24, 48 and 72 h post-treatment (Fig. 4). However, a difference in the overall trend can be seen between the two symbiotic states, in which aposymbiotic anemones generally showed a trend of decreased Ep_GIMAP expression at 48 h post-treatment, while symbiotic E. pallida showed a trend of increased Ep_GIMAP expression at 48 h post-treatment (Fig. 4). However, the only significant difference in expression between the LPS-treated symbiotic and aposymbiotic E. pallida was for Ep_GIMAP2 at the 48 h time point (P=0.031, ANOVA, Tukey's HSD).

A total of 11 Bcl-2 family member orthologs were identified in the E. pallida genome, comprising six anti-apoptotic members and five pro-apoptotic members based on the functions determined in A. millepora (Moya et al., 2016) (Table 2). Specifically, there were E. pallida orthologs for the confirmed pro-apoptotic members AmBax and AmBClRAMBO, as well as the predicted pro-apoptotic members AmBokB, AmBokC and AmBak. In addition, there were E. pallida orthologs for the confirmed anti-apoptotic members AmBclWA, AmBclWB, AmBClWD, AmBokA and AmMcl1-like, as well as the predicted anti-apoptotic member AmBclWC (Table 2).

An increase in caspase activity at 24 and 48 h post-wortmannin treatment indicates that wortmannin treatment was not only successful at inducing apoptosis in E. pallida but also provides time points at which expression differences of Ep_GIMAPs may be observed if these proteins are involved in regulating apoptosis. Interestingly, symbiotic state also impacted the level of apoptosis observed in E. pallida. Overall, the wortmannin-treated aposymbiotic animals exhibited higher levels of caspase-3 activity than the symbiotic animals, a difference that was significant at the 48 h time point. Furthermore, although wortmannin treatment resulted in a significant increase in caspase-3 activity in the symbiotic animals at 24 h post-treatment, this activity decreased to levels similar to those of the vehicle control by 48 h. Together, these results suggest that the presence of symbionts decreases the duration and strength of the apoptotic response initiated by wortmannin. This result is not surprising as previous work indicates that the presence of symbionts suppresses cellular pathways including LPS-induced nitric oxide production, NF-κB expression and autophagy (Detournay et al., 2012; Mansfield et al., 2017; Voss et al., 2019 preprint). Therefore, this result adds to the growing body of knowledge on how symbiosis impacts the cellular response of the host cnidarian to various external stimuli. Future work that explores the molecular regulation of apoptosis in cnidarians could provide further information on how its initiation is impacted by symbiotic state.

Although not measured, no evidence of bleaching was observed during this study with wortmannin treatment, unlike previous findings (Dunn et al., 2007). Therefore, the slight elevation in apoptosis levels observed in the symbiotic animals at 24 h may not have been enough to initiate the bleaching response, or alternatively this suggests that under the experimental conditions used, apoptosis does not induce bleaching. Future studies that utilize different concentrations of wortmannin and measure both bleaching and apoptosis levels could provide more information on this topic.

Overall, the qPCR results indicate there are no expression changes in the Ep_GIMAPs when apoptosis is not elevated, but a general downregulation in gene expression was observed upon elevated levels of apoptosis at 24 and 48 h post-treatment. Specifically, this downregulation was significant for Ep_GIMAP1, Ep_GIMAP2 and Ep_GIMAP4 (Fig. 2). Additionally, this downregulation was generally stronger in aposymbiotic anemones particularly at the 48 h time point when levels of apoptosis were the highest. Therefore, there is a clear negative correlation between apoptosis induction and Ep_GIMAP expression (Pearson's correlation test, P<0.001 for Ep_GIMAP1 and Ep_GIMAP2, P=0.029 for Ep_GIMAP3, P=0.002 for Ep_GIMAP4). The decrease in Ep_GIMAP expression suggests the involvement of these proteins in apoptosis where they could be either enhancing activity of anti-apoptotic proteins or inhibiting activity of pro-apoptotic proteins. From the vertebrate literature, it is known that GIMAPs interact with several members of the Bcl-2 protein family to regulate apoptosis (Nitta and Takahama, 2007). Specifically, some GIMAPs interact with anti-apoptotic family members, such as Bcl-2, while others interact with pro-apoptotic members such as Bax or Bak (Nitta and Takahama, 2007). Therefore, it is possible that Ep_GIMAPs may also interact with Bcl-2 family members to regulate apoptosis. From the characterization of Bcl-2 family proteins in the E. pallida genome, it was shown that orthologs to all proteins that were previously characterized in A. millepora are present. Specifically, E. pallida has several proteins that in A. millepora have been shown to be anti-apoptotic including Ep_BclWA, Ep_BclWB, Ep_BclWD, Ep_BokA and Ep_Mcl1-like (Table 2) (Moya et al., 2016). Additionally, E. pallida has several proteins that in A. millepora have been shown to be pro-apoptotic including Ep_Bax and Ep_BclRAMBO. Therefore, future studies that focus on determining the interaction partners of Ep_GIMAPs will shed light on their potential regulation of apoptosis through Bcl-2 family proteins.

Although this study represents the first in invertebrates to functionally test the relationship between apoptosis and GIMAPs, this is not the first time it has been proposed. Specifically, in a transcriptomic study on the responses of the oyster Crassostrea virginica to exposure to the bacterium Roseovarius crassostreae, GIMAP genes and inhibitor of apoptosis (IAP) genes were both differentially expressed in resistant and susceptible oysters, providing a link between GIMAPs and apoptosis (McDowell et al., 2014). Additionally, in a study characterizing GIMAPs (also referred to as AIG1) in invertebrates, it was suggested that AIG1 proteins may inhibit endoplasmic reticulum stress induced apoptosis (Guerin et al., 2019). Therefore, although much of the current literature is speculative, there is growing interest in the link between invertebrate GIMAPs and apoptosis. As this study represents the first to directly explore this idea, it provides foundational information that can be expanded upon in future studies.

The qPCR data in response to wortmannin treatment also illustrate the potential functional diversification in cnidarian GIMAPs. In contrast to the other Ep_GIMAPs, Ep_GIMAP3 does not show any statistically significant differences in expression upon induction of apoptosis. Therefore, it may play a role in other cellular processes. In vertebrates, each GIMAP protein serves a unique role in the cell (Nitta and Takahama, 2007) and as such it would not be surprising if the same were true of cnidarian GIMAPs. However, as there was a trend for downregulation in the gene expression data and a small sample size (n=3) was used, it could be worth increasing the sample size in future studies to further confirm these results.

It is known that at least one vertebrate GIMAP interacts with an autophagic protein and therefore the same could be true of cnidarian GIMAPs (Pascall et al., 2013). In this study, rapamycin treatment of symbiotic anemones resulted in a significant increase in the expression of Ep_LC3. This provides evidence that the rapamycin treatment effectively induced autophagy in E. pallida, but this response was only significant in symbiotic animals, suggesting that the presence of symbionts amplified the effect of the rapamycin treatment, and resulted in higher levels of autophagy. This is an interesting result given a recent study demonstrating that autophagy is downregulated by symbiosis and that the presence of symbionts activates mTOR activity, which inhibits autophagy (Voss et al., 2019 preprint). Therefore, as rapamycin inhibits mTOR activity to promote autophagy, it is possible that the symbiotic animals saw a larger increase in autophagy upon rapamycin treatment than aposymbiotic animals, where mTOR activity was inherently lower and autophagy levels higher. Overall, these results demonstrate that symbiotic state impacts the response to autophagy induction by rapamycin.

The Ep_LC3 expression results inversely correlated with the expression changes in Ep_GIMAP1 and Ep_GIMAP2 (Pearson's correlation test, P<0.001 for Ep_GIMAP1 and P=0.031 for Ep_GIMAP2). Specifically, there was a significant downregulation in both the AR and SR treatments compared with the respective vehicle controls. This suggests that Ep_GIMAP1 and Ep_GIMAP2 may regulate proteins that are actively involved in autophagy. Interestingly, this response was stronger in the symbiotic animals, which corresponded to the higher Ep_LC3 expression levels, and therefore higher levels of autophagy. This provides further evidence for the involvement of Ep_GIMAP1 and Ep_GIMAP2 in autophagy. Future work that investigates the localization of Ep_GIMAPs during autophagy induction or protein interaction partners could shed further light on their role in this cellular process.

In contrast to Ep_GIMAP1 and Ep_GIMAP2, there was no significant effect of rapamycin treatment on Ep_GIMAP3 and Ep_GIMAP4 expression, indicating that they may not be involved in autophagy or may play a different role from the other Ep_GIMAPs. Interestingly, the response of the Ep_GIMAPs to rapamycin was directly correlated to genomic location, suggesting there could be co-regulation of Ep_GIMAP1 and Ep_GIMAP2 as well as Ep_GIMAP3 and Ep_GIMAP4.

The lack of significant differences in Ep_GIMAP expression in response to LPS treatment implies their function may not be directly involved in the intracellular signaling induced by this immune stimulant. These results are in stark contrast to those previously obtained in A. millepora, where challenge with the bacterial cell wall component MDP caused massive upregulation of three GIMAP transcripts (Weiss et al., 2013). Our results are also different from those obtained in an RNA-Seq experiment on immune stimulation in the snail Biomphalaria glabrata, where LPS treatment resulted in increased expression of four GIMAP genes, suggesting that GIMAPs play a role in molluscan immunity (Zhang et al., 2016). There could be several reasons for these differences. First, in this study, the animals were immersed with the immune stimulant rather than injected and therefore this difference may result in a weaker response. Additionally, in regards to the A. millepora study, although both are immune stimulants, MDP and LPS activate unique cellular responses (Girardin et al., 2003; Pålsson-McDermott and O'Neill, 2004) and therefore GIMAPs may not be a part of LPS-induced signaling in cnidarians. Lastly, other studies utilizing LPS exposure in cnidarians have shown a lack of differential expression for immune genes (Kvennefors et al., 2010; Mali et al., 2006). It is also important to consider the concentration of LPS used and the time points; the lack of Ep_GIMAP expression differences may imply that a greater concentration of LPS must be used, or that the time frame for exposure should be altered. Lastly, a small sample size (n=3) was used in this study and therefore increasing the number of anemones could help to clarify the trends observed in the data. Future studies that utilize different concentrations of LPS, a greater number of time points and a larger sample size will be necessary to definitively confirm that GIMAPs do not play a role in LPS-induced signaling. Additionally, studies that utilize live bacteria or a greater diversity of immune stimulants will help provide a better understanding of the potential role of GIMAPs in immunity.

Together, the results from treatment with wortmannin, rapamycin and LPS indicated functional divergence in Ep_GIMAPs and suggest that they do not play a direct role in LPS-induced signaling, but instead may be involved in the regulation of the related processes of apoptosis and autophagy. Specifically, the common pattern of downregulation of Ep_GIMAP1 and Ep_GIMAP2 under both induction of apoptosis and autophagy suggests that these proteins may be interacting with some common regulator of the two processes. The interrelated nature of cnidarian apoptosis and autophagy is one that has previously been proposed (Dunn et al., 2007) and is also becoming increasingly prevalent in the vertebrate literature. Specifically, in vertebrates, it has been shown that proteins that are traditionally associated with apoptosis, such as Bcl-2, caspases and beclin-1, also influence the expression and activity of autophagic proteins (Esteve and Knecht, 2011; Kang et al., 2011). Thus, given the similarities in gene expression profiles, it is possible that Ep_GIMAP1 and Ep_GIMAP2 may interact with one of these common regulators. As both apoptosis and autophagy contribute to bleaching (Dunn et al., 2007; Hanes and Kempf, 2013) and disease resistance (Fuess et al., 2017), GIMAPs may function in the breakdown of symbiosis and disease, despite the lack of response to an individual immune stimulant.

Overall, this work represents the first targeted functional study of cnidarian GIMAPs and reveals their potential role in regulating processes that are involved in cnidarian symbiosis and the disease response. Subsequent studies that investigate protein–protein interactions or localization of Ep_GIMAPs will help to elucidate their role in these processes and will increase our overall understanding of these proteins in basal metazoans.

We would like to extend a special thanks to Eli Meyer for the reassembly of the E. pallida transcriptome.