Sea anemones (Exaiptasia diaphana previously known as Aiptasia pallida) were purchased from a supplier (Wards Natural Science, Rochester, NY, United States) and acclimated to laboratory conditions (45 L recirculating artificial seawater at 22°C and 36 ppt salinity) for 1 week. Anemones (100 individuals) were randomly subdivided into a control group or 20 ppb treatments of 17β-estradiol (E2), testosterone [T], cholesterol, benzyl butyl phthalate (BBP), or oxybenzone (BP-3). Toxicant purity was pharmaceutical secondary standard grade for E2, T, cholesterol, and BP-3 while BBP was analytical standard grade (Sigma-Aldrich, St. Louis, MO, United States). Acetone (150 μg/L) was used to solubilized toxicants were prior to dilution in seawater. Acetone has previously been used as a wetting agent (Morgan et al., 2001; Morgan et al., 2012). The 20 ppb concentration was chosen to induce responses and does not necessarily reflect environmental conditions even though BP-3 has been detected at concentrations >20 ppb (Downs et al., 2016) and cholesterol at 11ppb (French et al., 2015). All treatments were nominal concentrations for 4 h in 1L artificial seawater under ambient laboratory lighting during the early spring season. The 4-h exposure was chosen to detect responses at an early stage of exposure.

Complementary DNA sequences previously generated by Representational Difference Analysis (RDA) (Edman et al., 1997; Hubank and Schatz, 1999; Pastorian et al., 2000) were screened for candidate ESTs to be used in this investigation. Candidate ESTs were selected based on sequence homology (BLASTX analysis) to genes of interest (GOI) with known functions involving sex hormone synthesis, sterol binding, lipid transport, and/or indirect sterol interactions. Subsets of these ESTs have previously been used to characterize cnidarian stress responses (Morgan et al., 2012; Morgan et al., 2015).

Anemones from each treatment were pooled and total RNA was isolated using Trizol (Invitrogen, United States). Two milliliter phase-lock gels (5′Prime, United States) aided in the recovery of the aqueous phase. Total RNA was DNase I digested and purified (New England BioLabs, United States). Messenger RNA (mRNA) was isolated (Oligotex, Qiagen, United States) from 100 μg of DNase I treated total RNA. First stand synthesis used SuperScript IV (Invitrogen, United States) along with random hexamers and oligo-dT primers to reverse transcribe 1 μg of poly-A enriched mRNA. The RT reaction conditions were 1 h at 37°C, followed by 1 min at each temperature between 42 and 50°C.

A QuantStudio 7 Flex Real-Time PCR system (Applied Biosystems Waltham, MA, United States) used a SYBR Green-based assay to perform qPCR. A 1/100 dilution of first-strand synthesis reactions were used as templates for all qPCR reactions. Primers for each GOI were created using Primer3 (https://primer3.org) (Table 1). Components for each 20 µl qPCR reaction included: 10 µl Luna^® Universal qPCR Mix (New England BioLabs, United States), 2.5 µl Forward primer (10 µM), 2.5 µl Reverse primer (10 µM), 2.5 µl dH₂O, 2.5 µl sample. Thermocycling conditions were 1 cycle at 95°C for 1 min; 40 cycles of 95°C for 15 s and then 60°C for 30 s; and concluding with 1 cycle of melt curve analysis. Four replicate reactions were used for analyzing the relative expression of each GOI. Ribosomal protein L11 (RPL11) was used as qPCR reference gene (Kenkel et al., 2011). Melt-curve analysis, primer efficiencies, and gel electrophoresis confirmed specificity of priming. Replicate Cq values were averaged to determine ∆Cq and ∆∆Cq for each treatment and GOI. All ∆Cq and ∆∆Cq values are based on the consistent expression of RPL11 across all treatments. The ∆∆Cq method was used to determine the differences between targeted GOIs and a single reference gene (Bustin et al., 2009). Univariate ANOVA was performed on ∆∆Cq data to identify significant expression of individual GOIs across treatments. Similarities in variance between treatments were determined by Levene’s Test of Equality of Error Variances. If the variance between treatments was similar, then the Student-Neuman-Keuls (SNK) posthoc test was performed to determine which treatment(s) were significantly different from the rest. If variance between treatments was different, then Tamhane’s T2 posthoc test was applied since Univariate ANOVA is generally insensitive to heteroscedasticity.

Human protein homologs have previously been used as substitutes for cnidarian proteins in molecular modelling (Khalturin et al., 2018). Some cnidarian proteins have such significant homology to human proteins that they are even capable of stimulating human responses (Dunn et al., 2006; Wood-Charlson and Weis, 2009; Duffy and Frank, 2011; Dani et al., 2014; Brennan et al., 2017; Mansfield et al., 2017). The most similar vertebrate structures (according to BLASTX E-value) were used for selecting which Protein Databank (PDB) files should be used for initial docking experiments. Sequence alignments of cnidarian and relevant homologs determined which relevant protein sequences (i.e., PDB files) were to be used as representatives of corresponding cnidarian proteins. Protein coordinates (PDB files) for homologous constructs were downloaded from the PDB, and any existing ligands or small molecules (e.g., water or salt ions) were removed. Cnidarian protein structures were created with the SWISS-MODEL web server (https://swissmodel.expasy.org) using the homologs as templates. Ligands were created in Avogadro (Hanwell et al., 2012) and energy-minimized prior to docking. Docking was performed using AutoDock Vina (Trott and Olson, 2010); the search space was restricted to the known ligand binding site. After docking, the results were analyzed in UCSF Chimera (Pettersen et al., 2004) with the Find Clashes/Contacts function which was set to detect any atoms within 0.4 Å of the ligand with the best docking score.

BLASTX analysis of ESTs identified eleven candidate GOIs with of functions associated with steroidogenesis, gametogenesis, cholesterol transport, immunity, phagocytosis, or the Hedgehog signaling pathway (Table 2). These GOIs were: 17β-hydroxysteroid dehydrogenase type 14 (17βHSD14), 17β-hydroxysteroid dehydrogenase type 12 (17βHSD12), Niemann-Pick C type 2 (NPC2), Equistatin (EI), Complement component C3 (C3), Cathepsin L (CTSL), Patched domain-containing protein 3 (PTCH3), Smoothened (SMO), Desert Hedgehog protein A (DHH), zinc finger protein Gli2 (GLI2), and Vitellogenin (VTG). Two GOIs (17βHSD14; 17βHSD12) are associated with steroidogenesis. 17βHSD14 converts E2 into E1 and T into androstenedione (A4), while 17βHSD12 performs the reverse reactions, E1 into E2, and A4 into T (Tarrant et al., 2003; Luu-The et al., 2006; Shikina et al., 2016; Salah et al., 2019). Three GOIs (CTSL, EI, NPC2) are associated with lysosome/endosomes and/or symbiosome/phagosome (Dani et al., 2014; Mohamed et al., 2016; Dani et al., 2017; Hambleton et al., 2019). One GOI (C3) is associated with innate immunity (Miller et al., 2007; Kvennefors et al., 2010; Palmer and Traylor-Knowles, 2012; Ocampo et al., 2015). Four GOIs (PTCH3, SMO, DHH, GLI2) are representative of Hedgehog (HH) signaling which is highly conserved within vertebrate mammals (Sinkovics, 2015). One GOI (VTG) is the precursor of the egg yolk protein vitellin (Matozzo et al., 2008; Shikina et al., 2013; Lotan et al., 2014).

For the GOIs associated with steroidogenesis, 17βHSD14 expression was significantly up-regulated (p < 0.01) in each exogenous chemical treatment compared to the control, whereas 17βHSD12 expression was significantly up-regulated (p < 0.01) only in the BP-3 and BBP treatments (Figure 1). VTG exhibited significant down-regulated expression (p < 0.01) in E2, T, and cholesterol treatments (see Figure 1).

GOIs associated with phagocytosis had significant differences in expression in various treatments when compared to the control condition (Figure 2). Cathepsin (CTSL) expression was significantly up-regulated in the BP-3 and BBP treatments (p < 0.01). Equistatin (EI) was significantly up-regulated in E2, BP-3, and BBP treatments (p < 0.01) and significantly down-regulated in the T and cholesterol treatments (p < 0.01) (see Figure 2). NPC2 was significantly down-regulated in cholesterol (p < 0.01) and significantly up-regulated (p < 0.01) in BP-3 and BBP treatments (see Figure 2). Lastly, C3 was significantly up-regulated in E2, T, and BP-3 treatments (see Figure 2).

For the GOIs associated with HH signaling, PTCH3 was significantly down-regulated (p < 0.01) in response to all three sterols, but significantly up-regulated (p < 0.01) in BP-3 and BBP treatments (Figure 3). SMO was significantly down-regulated (p < 0.01) in all treatments except E2 (see Figure 3). DHH was significantly down-regulated (p < 0.01) in T and cholesterol treatments, but significantly up-regulated (p < 0.01) in BP-3 and up-regulated in BBP treatments (p < 0.05) (see Figure 3). GLI2 was significantly up-regulated (p < 0.01) in all treatments except E2 (see Figure 3).

Univariate ANOVA reveals all sterol treatments were significantly different from each other as well as the xenobiotics (p < 0.01). While the BP-3 and BBP treatments were significantly different from all the sterols (p < 0.01), they were not significantly different from each other (Figure 4).

In order to better understand how sterols and xenobiotics are able to interact with targeted proteins, in silico modelling was used to predict the most energetically favorable ligand/protein binding interactions. The binding affinities of each ligand to the proteins of interest were estimated using docking simulations. Comparative genomic analyses indicate that cnidarian genomes show important similarities to vertebrates in gene content, genomic structure, and organization (Putnam et al., 2007). However, the cnidarian proteome is not as well characterized as the human proteome. Therefore, macromolecular structures needed for the in silico modelling were generated in one of two ways. The first model was the structure in the Protein Databank with the highest degree of homology to the GOI. These PDB files were: 6oeu (PTCHD3), 5kwy (NPC2), 6xbm (SMO), 5hs6 (17βHSD14), and 2et6 (17βHSD12). The best PDB structures corresponded to the human sequences in all cases except for 17βHSD12, where the protein with the closest homology was from Candida tropicalis. The second model was generated using SWISS-MODEL, with the Exaiptasia sequence as input and the closest PDB structure as a template. Each ligand was docked to both protein models. The docking algorithm used flexible fitting with an energy-minimized ligand to determine possible binding conformations and relative binding affinities. All ligands bind to each protein, although the strength of the interactions varied. The binding free energy of the best docking pose for each ligand/receptor combination is shown in Table 3.

Since a crystal structure of 17βHSD14 with a similar ligand (estrone) is available on the PDB (5HS6), it is possible to compare the binding site of an experimentally-determined protein-ligand structure with our docked models (Laskowski and Swindells, 2011). Figure 5 shows a comparison between human 17βHSD14 with bound estrone, human 17βHSD14 with docked estradiol, and Exaiptasia 17βHSD14 with docked estradiol. The search space was confined to the binding pocket for each protein-ligand interaction, and the orientation of the ligand within the pocket was fairly well-conserved (see Figure 5). Additionally, the residues that interact with each ligand were well-conserved between cnidarian and homologous models. For example, estrone and estradiol formed hydrophobic contacts with His93, Gln148, Trp192, and Leu195 in both crystal structure and the docked model of human 17βHSD14. In the cnidarian protein, these residues correspond to His102, Gln157, Trp201, and Leu204, and are predicted to interact the estradiol ligand in the docked model as well (see Figure 5).

Endocrine disruption involves chemicals that alter steroidogenesis which interfere with the hormonal signaling leading to down-stream physiological abnormalities in reproduction, development, increased risk of disease, and/or interference with immune and nervous system functions (Heindel et al., 2012). The 17β-hydroxysteroid dehydrogenases (17βHSDs) are involved with steroidogenesis of both estrogens and androgens. 17βHSDs have been detected in cnidarians (Tarrant et al., 2003; Blomquist et al., 2006; Shikina et al., 2016). 17βHSD12 and 17βHSD14 represent the last step in the steroidogenesis pathway prior to the aromatization of T into E2. 17βHSD12 reduces estrone (E1) into E2 as well as androstenedione (A4) into T, while 17βHSD14 performs the opposite reaction by oxidizing E2 into E1, and T into A4 (Luu-The et al., 2006; Caldwell et al., 2012; Sivik et al., 2012; Thomas and Potter, 2013; Shikina et al., 2016; Salah et al., 2019). Concentrations of competing sterols may drive which sterol gets converted by a specific 17βHSD (Lathe and Kotelevtsev, 2014). E2 concentration is known to regulate transcription of 17βHSDs (Rotinen et al., 2009). Transcriptional profiles herein demonstrate the organism (regardless of sex) responds to the environmental “overload” of E2 or T by modulating the transcription of the enzyme in order to transform both sex hormones into less potent sterols. These profiles also reveal that cholesterol and xenobiotics (BBP and BP-3) interfere with steroidogenesis by altering normal transcriptional demands for 17βHSDs (see Figure 1).

Cnidarians get cholesterol and other sterols by either dietary intake and/or from their endosymbiont (Revel et al., 2016; Hambleton et al., 2019). Cholesterol transport and metabolism is tightly regulated (Maxfield et al., 2016). Due to lipid-rich tissues, it has been proposed that lipophilic compounds will easily diffuse into cnidarian tissues (Peters et al., 1997). Sterols are most concentrated in cnidarian epidermis and some sterols are known to increase in concentration after exposure to environmental pollutants (Revel et al., 2016; Stien et al., 2020). The lack of significant expression of C3 in the cholesterol treatment is consistent with the fact that it is a natural constituent of animal membranes.

Five GOIs (C3, PTCH3, CTSL, EI, NPC2) have functional associations with the endocytic processes whereas a sixth GOI (DHH) is a secreted morphogen associated with exocytosis/secretion. C3 is associated with phagocytosis of foreign particles into cnidarians (Palmer and Traylor-Knowles, 2012), while PTCH regulates functions of: receptor-mediated endocytosis of sterol-protein ligands; down-stream morphogens; and the removal of oxysterols (Incardona et al., 2000; Strutt et al., 2001; Corcoran and Scott, 2006; Zhong et al., 2014). Three GOIs (CTSL, EI, NPC2) are associated with lysosomes and phagocytosis. Cathepsins are lysosomal proteases and EI is the cnidarian homolog of an inhibitor of cathepsins (Lenarčič et al., 1997; Štrukelj et al., 2000). Lysosomes contribute to the termination of hormone effects by degrading proteins interacting with the hormone (Qualmann et al., 2000; Totta et al., 2014). NPC2 functions as a key regulator of sterol transport between Symbiodinium and cnidarian hosts which is representative of stable symbioses in a structure known as the symbiosome (Infante et al., 2008; Vance, 2010; Ganot et al., 2011; Dani et al., 2014; Bucher et al., 2016; Revel et al., 2016; Hambleton et al., 2019; Lu et al., 2020). Seven GOIs (NPC2, 17βHSD14, C3, PTCH3, SMO, DHH, VTG) have known linkages to endoderm. (Kvennefors et al., 2010; Shikina et al., 2013; Sanders et al., 2014; Shikina et al., 2016).

Five GOIs (17βHSD14, 17βHSD12, NPC2, PTCH3, SMO) that have sterol binding pockets were used in silico modelling. In silico modelling has previously demonstrated phthalates and steroids without an aromatic A-ring are capable of binding to sterol-binding proteins (Sheikh and Beg, 2017; Khalturin et al., 2018). Patched genes have multiple functionally redundant homologs that cannot discriminate between ligands (Carpenter et al., 1998; Zhulyn et al., 2015; Hu and Song, 2019). Justification for model docking is based on the fact that cnidarian genomes have representative vertebrate homologs of relevant proteins (Miller et al., 2007; Putnam et al., 2007; Kvennefors et al., 2010). In silico modelling demonstrates the lowest energy docking poses for the crystal structures of the homologous proteins and the models of the Exaiptasia proteins correlate well with each other (see Table 3). The binding scores for the Exaiptasia protein models generally follow the same trends as the homologs from the PDB. In both cases, the steroid ligands (cholesterol, E2, and T) bound more strongly than BP-3 or BBP, although the differences in binding energy seem to be smaller for the 17βHSD12 and 17βHSD14 proteins. This result was expected, since BBP and BP-3 are considerably smaller and as such, would make less contact with binding pocket residues (see Table 3). Although there are a few minor differences between homologs, such as hydrogen bonding patterns, the orientation of the ligand in the 17βHSD14 binding pocket and the side chains involved in binding are strikingly similar (see Figure 5). All sterols have similarities in orientations, binding affinities, and expression profiles (see Table 3 and Figure 5). 17βHSDs are known to exhibit variability in substrate binding (Blomquist et al., 2006; Lathe and Kotelevtsev, 2014).

Data herein provides new information about the diversity of ligands, including xenobiotics that bind to NPC2, PTCH3, SMO, and the 17βHSDs. The differential expression of these sterol-binding proteins suggests that these proteins handle natural sterols and xenobiotics differently. Laboratory-induced exposure to cholesterol, E2, and T requires modulating transcript copy numbers to handle the environmental flux of natural ligands. However, laboratory-induced exposures of BBP and BP-3 appear to require transcriptional replacement for proteins after binding xenobiotics. All GOIs except SMO exhibit up-regulated expression in the BP-3 and BBP treatments (see Figures 4E,F).

Hedgehog (HH) signaling is a highly conserved developmental pathway that relies on cholesterol interacting with specific sterol-binding proteins. Endocrine disruption is known to interfere with HH signaling (Johansson and Svingen, 2020). Developmentally-regulated signaling pathways involving apoptosis, embryogenesis, cell proliferation, and development of diseases present in higher animals are also found in cnidarians including Exaiptasia (Technau et al., 2005; Matus et al., 2008; Sinkovics, 2015; Bucher et al., 2016). Four GOIs (PTCH3, SMO, DHH, GLI2) can be linked to HH signaling. E2 acts as a direct HH pathway antagonist (Lipinski and Bushman, 2010) and the E2 expression profile in Figure 3 reaffirms this inhibitory effect. Recognized inhibitors of HH signaling also include the well-characterized steroidal alkaloid cyclopamine (Chen 2002). The down-regulated expression profiles of PTCH3 in the E2, T, and cholesterol treatments and the down-regulated expression of SMO in the T and cholesterol treatments suggests all three sterols exhibit similar inhibitory effects (see Figure 3). The HH pathway is also responsive to xenobiotics as well as sterols (Johansson and Svingen, 2020). Phthalates are known to impact DHH and PTCH expression, and BP-3 is considered estrogenic (Thibaut and Porte, 2004; Borch et al., 2006; Downs et al., 2016; Hong et al., 2016; LaPlante et al., 2018; Matouskova et al., 2020). Although PCTH3 and DHH were down-regulated in response to sterols, PTCH3, DHH, and GLI2 were up-regulated in response to BP-3 and BBP (see Figure 3), suggesting opposing effects of sterols and xenobiotics on HH signaling, presumably with sterols exerting antagonists effects and xenobiotics exerting agonist effects. Regardless, the expression profiles in Figure 3 highlight the inherent vulnerability of the HH pathway to diverse chemical signals.

Gametogenesis is associated with HH signaling. HH signaling molecules are considered stem-cell factors (Zhang and Kalderon, 2001). Evidence exists in other metazoans that secreted molecules in the HH signaling pathway are one of the primary regulators of stem cells which are necessary precursors of oogenesis (Zhang and Kalderon, 2001). Shikina and Chang (2016) identifies cnidarian germline stem cells as originating in endodermal mesentery. HH signaling GOIs (PTCH3, SMO, DHH) are known to be restricted to the endoderm of male and female gonophores (Sanders et al., 2014). Germline stem cells are the essential cells responsible for this process in various metazoans. In Drosophila, somatic cells in the ovary of cannot proliferate normally in the absence of Smoothened activity (Zhang and Kalderon, 2000), and excessive HH signaling can result from aberrant expansion of stem cell pools (Zhang and Kalderon, 2001). Results show all exogenous treatments (except E2) caused down-regulation of SMO (see Figure 4). With the exception of PTCH3, results also show that HH signaling GOIs in the E2 treatment are not significantly expressed (see Figure 3).

Exaiptasia reproductive strategies include both sexual and asexual reproduction (Clayton, 1985; Grawunder et al., 2015). Subsets of anemones used for each treatment were taken from the same population as the control anemones. Although not confirmed, it is highly likely that all anemones used in this investigation were clonal. Detecting the presence of sex-specific tissues enhances the ability to interpret anemone responses to the EDCs. Although the sex of the colony was not known, expression profiles of the GOIs could provide clues about the sex of these animals. Vitellogenin is a female-specific protein commonly found across the animal kingdom and is the most abundant protein associated with oocyte development in cnidarians (Shikina et al., 2013; Levitan et al., 2015). Hormonal and/or environmental cues are known to stimulate release of oocytes (Levitan et al., 2015). Exogenous E2 is known to influence pedal laceration and egg bundle formation in cnidarians (Tarrant et al., 2004; Thorn et al., 2015). Blue-light cues have also been used to induce Aiptasia sp gametogenesis in a laboratory setting (Grawunder et al., 2015).

This investigation also did not attempt to identify the sex of the anemones nor their reproductive status, per se. Hermaphroditism and even Trioecy has been observed in closely related species (Schlesinger et al., 2010; Armoza-Zvuloni et al., 2014). A small amount of data from a previous study suggests male tissues from hermaphroditic corals also express VTG (Hayakawa et al., 2005). Results reveal VTG expression in the control treatment suggesting anemones are female and already transcribing VTG (see Figure 1). Decreases in vitellogenin in females can be useful for identifying EDCs (Ankley et al., 2002). If these anemones were male, then up-regulation of vitellogenin would have been expected. Up-regulation of vitellogenin in males is widely recognized as the prominent estrogenic effect in many aquatic animals exposed to exogenous E2 (Denslow et al., 1999; Matozzo et al., 2008). These results cannot confirm nor reject the possibility that male reproductive tissues were also present.

A second piece of evidence suggest that these anemones were female was the expression of 17βHSD14 in E2 which was significantly up-regulated (p < 0.0001) compared to the T, cholesterol, BBP, and BP-3 treatments. Female anemones perform oogenesis which is stimulated by the presence of E2 so when exogenous E2 is introduced, additional 17βHSD14 is required to reduce E2 potency by transformation into E1. In this investigation, VTG expression was crucial for differentiating E2 exposure from other EDC exposures. The expression profiles of 17βHSD14 and 17βHSD12 can be indirect measures of cellular demands for sex hormones. When exposed to exogenous sex hormones, transcription of 17βHSD14 is up-regulated while 17βHSD12 is not significantly different from the control. These results suggest these anemones were females responding to the EDCs in this investigation. Future investigations using these GOIs, coupled with lab-controlled methodologies to induce gametogenesis in recognized clonal lines of anemones should provide greater resolution in how these organisms respond to EDCs.