Transcriptomes and genomes of Symbiodiniaceae and symbiotic anthozoan species were obtained from publicly available sources (Symbiodiniaceae—Table 1, anthozoans—Table 2). For anthozoan species, we selected the highest-quality gene sets which also covered a range of hexacorals and Exaiptasia diaphana (n = 21 datasets). Due to the reduced number of available datasets for Symbiodiniaceae, all available assemblies (n = 28) were analysed, regardless of their quality. This spanned six genera, with most covering associations with hexacorals. BUSCO analysis, a metric to measure genome/transcriptome completeness through the identification of near-ubiquitous single-copy orthologs, was used to quantify the quality of these datasets (Manni et al., 2021). The range of BUSCO scores within the cnidarian datasets using the eukaryote database (v5) was between 78 and 93.8% complete, with an average of 87.8% (Figure 2). BUSCO completeness scores of the Symbiodiniaceae datasets using the eukaryote database (v5) ranged from 0.8 to 80.8%, with an average score of 49.2% (Figure 3). There was a particular focus throughout this analysis on the sea anemone Exaiptasia diaphana (“Aiptasia”), a widely used model organism for cnidarian-dinoflagellate symbiosis that can be maintained in an aposymbiotic state and re-colonised with multiple Symbiodiniaceae species (Weis et al., 2008).

A literature search was used to identify conserved domains of the PI kinases and phosphatases. The conserved domains were used to search the genomic and transcriptomic resources using the “aphid” R-package (Wilkinson, 2018). Specifically, 17 Hidden Markov Models (HMMs) that provided complete coverage of the PI signalling pathway were used as queries to identify proteins of interest (Table 3). Resulting proteins with a score >20 (an observed cut-off in quality and accuracy below this score) were run through the InterProScan plugin in Geneious (v2019.2.3) to identify additional domains. To verify the presence or absence of PI pathway genes, we also used: (1) a custom BLAST analysis against a database of the 28 Symbiodiniaceae and 21 cnidarian species (E-value cut-off 10^–2); (2) motif identification through thorough literature searches; (3) protein alignments using Muscle in Geneious; and (4) similarity scores produced during alignments. Moreover, genes identified above were assigned to gene families, or orthogroups, using OrthoFinder 2.5.2 (Emms and Kelly, 2015). In complex cases, such as in the distinction between phosphatidylinositol 3- and 4-kinases, where the same domain is used to extract distinct target proteins of interest, maximum-likelihood phylogenetic trees were reconstructed using the conserved protein domains with RAxML (v4.0 in Geneious, Protein model: Gamma Blosum 62) (Supplementary Figure 1).

OrthoFinder was used for ortholog assignment of cnidarians and Symbiodiniaceae separately (Emms and Kelly, 2015). Four different analyses were carried out: comparison within the available Symbiodiniaceae datasets; within selected symbiotic cnidarians; within symbiotic and non-symbiotic cnidarians; and within the Symbiodiniaceae plus non-symbiotic unicellular organisms (for symbiotic/non-symbiotic comparison). Protein groups were further delineated by similarity to each other using BLAST, which was especially useful when analysing Symbiodiniaceae datasets, and when attempting to identify differences in the identifiable pathways between symbiotic and non-symbiotic cnidarians/Symbiodiniaceae.

The methods described above were used to identify proteins involved in PIP-binding within datasets of the target species E. diaphana and Breviolum minutum, as well as in datasets of selected species spanning a variety of symbiotic/non-symbiotic symbiont and cnidarian types (Table 4). E. diaphana and B. minutum were utilised due to their common association in nature and the use of them as a model organism (Weis et al., 2008). Previous literature on PIP-binding proteins were used to identify regions of conservation (HMMs) (Cernikova et al., 2019, 2020). Within each of these databases, 14 HMMs were queried. Proteins with at least one of the required PIP-binding domains were included in the gene enrichment analysis (Table 5). The PIP-binding proteins for each species were run through Pannzer2 to ascertain Gene Ontology (GO) terms (Törönen et al., 2018). The “TopGO” R-package (Alexa and Rahnenfuhrer, 2016) was used for the gene enrichment analysis to compare the identified PIP-binding proteins to the entire GO repertoire of biological processes, molecular functions, and cellular compartments within each target species. The GO terms with adjusted p-value < 0.05 were further condensed using Revigo semantic similarity clustering (Supek et al., 2011) (Settings: Small (0.5) similarity, p-value, whole Uniprot database, “SimRel”: a semantic similarity measure). A literature search was then completed to find a link for the significant biological terms to either a mutualistic or parasitic interaction.

Across the symbiotic anthozoans analysed there was a complete and conserved PI pathway, analogous to that found in multicellular eukaryotic model organisms (Michell, 2011; Balla, 2013; Figure 2). PI3- and PI4- kinases constitute enzymatic families in higher eukaryotes that are subdivided into multiple classes. Members from each group of the kinases as well as members of the 3- and 4-phosphatases were also identified across the surveyed anthozoans. Phylogenetic trees constructed for the PI3- and PI4-kinases identified each kinase subclass, suggesting the conservation of these kinases within Anthozoa (Supplementary Figure 1). Conversely, 5-phosphatases were not conserved throughout the symbiotic anthozoans, with the loss of INPP5B and SHIP1 in all taxa, while others such as SKIP, PIPP, INPP5E and inositol 1,4,5-trisphosphate 5-phosphatase (INPP5A) were absent from various lineages (Figure 2).

We did not identify any direct evidence of altered regulation of anthozoan host PI-kinases and -phosphatases in response to symbiosis in published datasets. Yet, the mechanism of activation of these enzymes—by phosphorylation—means that more commonly used “omics” technologies (transcriptomics, proteomics) may not identify altered activity in the PI signalling pathway. Interestingly, however, sodium and proton myo-inositol co-transporters are upregulated in symbiotic Aiptasia at the transcriptome (Lehnert et al., 2014) and proteome levels (Sproles et al., 2018b). Myo-inositol is the biologically active precursor to the PI pathway that plays a central role in signal transmission (Figure 1). Metabolomic analysis also found the accumulation of inositol in thermally stressed anemones, indicating that its clearance is important for maintaining the cnidarian-dinoflagellate symbiosis (Hillyer et al., 2016). Furthermore, proteomic analysis identified down-regulation of multiple sugar-inositol transporters, the proton myo-inositol co-transporter (HMIT), the myo-inositol/Na + co-transporter (SMIT1), and the myo-/chiro-inositol/glucose transporter (SMIT2) when the host was inoculated with the heat tolerant, but physiologically sub-optimal (Sproles et al., 2018b) symbiont Durusdinium trenchii (Matthews et al., 2017). Taken together, this suggests that the efficient transport of myo-inositol is essential for the establishment and maintenance of an optimal symbiosis.

Analogous to the PI pathway found within Apicomplexa (Wengelnik et al., 2018; Cernikova et al., 2019; Nakada-Tsukui et al., 2019; Cestari and Stuart, 2020), Symbiodiniaceae displayed no genes of the 4-phosphatase group (Figure 3). However, conversely to apicomplexans, Symbiodiniaceae showed a greatly increased diversity both of and within the 5-phosphatases, with multiple 5-phosphatase homologs with dual domain structure, identified within the majority of Symbiodiniaceae species (Supplementary Tables 1, 2). Anthozoan genomic and transcriptomic data showed an average of four and three proteins with the PIPK (used for identification of PIP5K proteins) domain, respectively. However, Symbiodiniaceae genomic and transcriptomic datasets both showed an average of 15 proteins containing the PIPK domain, despite the reduced BUSCO scores. As the lowest BUSCO score for anthozoans was 78%, we used that as a threshold to look at comparative numbers of Symbiodiniaceae PIPK domain-containing proteins, and the average proteins further increased to 21. Of particular note was the presence of a dual domain protein within 15 of the Symbiodiniaceae datasets, with both the PIPK domain and also the PI-PLC domain. Searches of the BLAST database showed no results for similar proteins in other species. The range of PI kinases discovered was diverse, with members from each kinase group identified in all but one Symbiodiniaceae species (Figure 3). The exception was “F. kawagutii #1” —a dataset of Fugacium kawagutii—that only returned a 0.8% complete BUSCO score. However, the lack of functional annotation of proteins from other dinoflagellates makes identification difficult. We relied primarily on sequence homology to the characterised 5-phosphatases (n = 10) and other well-annotated results achieved through NCBI BLAST analysis. Caution should be exercised when examining gene absence in this group because of the incompleteness of some of the datasets used (BUSCO scores found in Figure 3). Despite the sparseness of the available datasets, we were able to identify a completely functional PI pathway across the analysed Symbiodiniaceae species.

Transcriptomic analysis of four species spanning a diversity of four Symbiodiniaceae genera [Symbiodinium (ITS strain A2), Breviolum (ITS strain B2), Cladocopium (ITS strain C1), and Durusdinium (ITS strain D1)] revealed that six cellular pathways are common to all four genera, including the PI signalling pathway and the inositol phosphate metabolism pathway, highlighting the potential global significance of these pathways in Symbiodiniaceae (Rosic et al., 2015). Notable genes encoding phosphatidylinositol 4-phosphate 5-kinases (PIP5K) and phosphatidylinositol 4-kinase (PI4K) were identified across all four Symbiodiniaceae species analysed. PIP4K synthesises PI 4-phosphate (PI4P), which is vital for the regulation of chloroplast division (Bovet et al., 2001; Okazaki et al., 2015). The ubiquity of PIP4K across the Symbiodiniaceae datasets analysed here is therefore unsurprising, with exceptions likely reflecting the quality of available resources. Meanwhile, PIP5K is the main route for the production of phosphatidylinositol 4,5-bisphosphate, and has cellular roles in actin dynamics, endocytosis, exocytosis, and focal adhesion assembly (Tan et al., 2015; Nakada-Tsukui et al., 2019). PIP5K also plays an important role in endosomal trafficking, with the specific isoform PIPKIyi2 being localised to recycling endosomes to facilitate transmembrane protein recycling (Thapa et al., 2012). Exocytosis and endocytosis are important functions in signalling processes within symbiosis (Weis, 2008; Oakley et al., 2016; Rosset et al., 2021) and hence, the down-regulation of this protein when the alga is in culture could reflect an increased demand for vesicle transport after colonisation of a host. Meanwhile PIP4K synthesises PI 4-phosphate (PI4P), which is vital for the regulation of chloroplast division (Bovet et al., 2001; Okazaki et al., 2015). The ubiquity of PIP4K and PIP5K across the Symbiodiniaceae datasets analysed here is therefore unsurprising, with exceptions likely reflecting the quality of available resources.

A transcriptomic study of B. minutum demonstrated a significant down-regulation of PIP5K activity when the alga was in culture compared to in hospite, providing further support for the importance of PIP5K in the functionality of symbiosis (Maor-Landaw et al., 2020). A similar transcriptomic analysis explored the differences between the transcriptome of the thermally tolerant D. trenchii when free-living and in hospite in Aiptasia (Bellantuono et al., 2018). Our targeted analysis of the dataset from Bellantuono et al. (2018) revealed that many of the PI pathway-related genes were significantly differentially expressed between the symbiont when in hospite and in culture (Supplementary Figures 2A, B). Notably there was a large downregulation of PIP5K “2” and a large upregulation of PIP5K “9” and “7,” among other types, when the symbiont was in symbiosis compared to free-living. These observations of up-regulated PIP5K expression in symbiosis could reflect an increased demand for vesicle transport in symbiosis. It is also possible that, as observed in other host-microbe interactions (Baumgartner et al., 2000; Drecktrah et al., 2004; Leirião et al., 2005; Ham et al., 2011; Ebrahimzadeh et al., 2019), Symbiodiniaceae modifies PI-kinase expression as a molecular strategy to manipulate host PIP synthesis to facilitate host membrane dynamics and regulate cellular mechanisms essential for symbiosis persistence.

The range of proteins within the PI pathway was highly conserved between symbiotic and non-symbiotic cnidarians (Tables 1, 2 and Supplementary Table 3). Both coverage of the pathway and number of proteins within each ortholog were similar, with no significant increases or decreases in gene numbers between cnidarians that associate with symbionts (Supplementary Tables 3, 4).

In contrast, differences were observed between the Symbiodiniaceae and non-symbiotic unicellular algal datasets, including presence/absence of ortholog groups and varying numbers of isoforms/duplicates of individual orthogroups between symbiotic and non-symbiotic datasets. One orthogroup of interest was OG0006447, which was absent in all five of the non-symbiotic species, but present across all except one of the Symbiodiniaceae datasets (“F. kawagutii #1” —previously noted to be largely incomplete based on BUSCO analysis). This orthogroup was analysed further using NCBI BLAST, where it matched consistently and with a high sequence conservation to a phosphatidylinositol 5-phosphate 4-kinase, or PIP5K/PIPK type II protein. This protein has been described across model metazoans but has not yet been described within unicellular organisms; indeed some research has suggested that it is not present in unicellular organisms (Mathre et al., 2019; Raghu, 2021). This protein is responsible for a small proportion of the production of PI 4,5-bisphosphate, however, it is predicted that its primary role is the regulation of PI 5-phosphate (PI5P), which acts as its binding substrate (Rameh et al., 1997; Raghu, 2021). This protein is localised to the plasma membrane, vesicular compartments and nucleus, with roles in PI3K-based regulation of receptor-activated signalling, membrane transport, autophagosome regulation, and the control of PI5P levels in response to cues such as UV irradiation and oxidative stress (Keune et al., 2013; Vicinanza et al., 2015; Wang et al., 2019; Raghu, 2021). Interestingly, the loss of this protein in mouse models leads to the hyperactivation of the innate immune system (Shim et al., 2016). The presence of this protein across all but one of the 28 Symbiodiniaceae datasets in our study is noteworthy, and experimental research to elucidate its role in symbiosis would be valuable.

Another orthogroup (OG0015044) predicted to be a PTEN/TPTE2 3-phosphatase homolog showed a similar pattern. Although there were other orthogroups that showed the presence of PTEN/TPTE2 ubiquitously across symbiotic and non-symbiotic unicellular organisms, this separate orthogroup (OG0015044) indicates an expansion of this protein type within the Symbiodiniaceae datasets. Most of the proteins within this orthogroup showed predominantly PTEN hits when run against the NCBI BLAST database. PTEN is a negative regulator of PI3K and has a variety of roles within cells, including survival, proliferation and regulation of cellular architecture, and it is frequently disrupted in cancer (Song et al., 2012). The presence of an expanded PTEN homolog group is observed in the parasitic protist, Entamoeba histolytica, which has six PTEN homologs, while there is only one present in humans (Nakada-Tsukui et al., 2019). The expansion of the PTEN group within Symbiodiniaceae is indicative of an importance of PTEN functionality, given its role in cell survival and proliferation, and it is therefore another good candidate for future functional studies of the cnidarian-dinoflagellate symbiosis.

A few notable Symbiodiniaceae orthogroups (e.g., OG0012793, OG0014270) were also distinct between Symbiodiniaceae and non-symbiotic unicellular organisms; however, the proteins did not BLAST to any known annotations, or recovered only a few hits, with “unnamed” or “hypothetical protein” responses. When analysed using InterPro scan, these orthogroups contained just one domain, PIPK. However, without targeted experimental work, these protein groups cannot be identified beyond the PIPK group in general.

To further understand the downstream functions of the PI pathway and to unravel the potential role in symbiosis, we expanded the analysis to include enrichment of proteins with the capabilities of binding to members of the PI pathway (termed PIP-binding proteins). All 14 identified domains of interest were found in Aiptasia, whereas four domains were absent in Breviolum (ANTH, FERM, GOLPH3, PTB). Within Aiptasia, 913 PIP-binding proteins were identified, which was 3.12% of the total predicted proteins in the genome, resulting in 306 significantly (p < 0.05) enriched GO terms (Supplementary Table 5). Within Symbiodiniaceae (Breviolum sp. B1), 898 PIP-binding proteins were identified, which was 0.91% of the total predicted proteins in the assembly, resulting in 103 significant GO terms (Supplementary Table 6).

Putative functionalities of the PIP-binding proteins identified in Aiptasia and Breviolum that were predicted to be relevant to symbiosis regulation can be found in Figure 4. Of particular interest, there were seven and six significantly enriched GO terms related to vesicle functioning in Breviolum and Aiptasia, respectively. We hypothesise that cellular processes involving vesicles, such as vesicle-mediated transport, play an important role in the functioning of the cnidarian-dinoflagellate symbiosis.

Enriched GO terms can be found in Supplementary Tables 7–20. The GO term “COPII-coated vesicle budding” was significantly enriched within one of the symbiotic cnidarians (the coral Montipora peltiformis), but more noticeably, all analysed Symbiodiniaceae. One role of COPII vesicles is the transport of proteins from the endoplasmic reticulum, the site of PI synthesis (Agranoff et al., 1958; Sato and Nakano, 2007; Balla, 2013), suggesting an important role of PI in the symbiont. Two GO terms corresponding to immune responses (“innate immune response” and “immune response-activating cell surface receptor”) were found to be significantly enriched across four out of the five Symbiodiniaceae species analysed, whereas these GO terms were not found to be enriched across any of the cnidarians analysed. There is evidence of endosymbiotic Leishmania spp. modulating their host’s immune system through the secretion of exosomes into the host (Silverman and Reiner, 2011). Other research shows Leishmania spp. subverting host immune responses by targeted host pathways such as the PI 3-kinase/Akt signalling cascade to prevent the activation of macrophages (Ruhland et al., 2007; Lambertz et al., 2012). This evidence extends to the cnidarian-dinoflagellate symbiosis, with 16% of differentially expressed genes between symbiotic states being composed predominantly of genes involved in apoptosis and oxidative stress; which are critical functions of immune responses (Kvennefors et al., 2010; Matthews et al., 2017).

Another GO term that regularly appeared as significant is the “positive regulation of TORC1 signalling” in four out of the five analysed Symbiodiniaceae, and “positive regulation of TOR signalling” in three out the five symbiotic cnidarians analysed. TOR (target of rapamycin) signalling activity inhibits autophagy (Juhász et al., 2008), and the presence of significantly enriched GO terms in both the symbiotic cnidarian group and Symbiodiniaceae indicates active inhibition of autophagy across both groups. Transcriptomics supports this finding, with elevated levels of mTORC1 kinase in symbiotic Aiptasia compared to non-symbiotic animals, suggesting a mechanism similar to a feeding response in symbiotic cnidarians (Voss et al., 2019). Up-regulation of TOR-related genes and transcripts during colonisation of plants by mycorrhizae suggests its role in other symbioses (Nanjareddy et al., 2016; Arthikala et al., 2021).

Understanding the implications for the enriched GO terms in the cnidarian-dinoflagellate symbiosis is complex, as just because a function does not appear with a high frequency does not mean its role in the symbiosis is unimportant. Further experimental work, including the inhibition and knock-out of specific genes, will help to elucidate the roles of these highlighted GO terms in the functioning and maintenance of the symbiosis.