Gene expression alterations from reversible to irreversible stages during coral metamorphosis

Abstract

For corals, metamorphosis from planktonic larvae to sedentary polyps is an important life event, as it determines the environment in which they live for a lifetime. Although previous studies on the reef-building coral Acropora have clarified a critical time point during metamorphosis when cells are committed to their fates, as defined by an inability to revert back to their previous states as swimming larvae (here referred to as the “point of no return”), the molecular mechanisms of this commitment to a fate remain unclear. To address this issue, we analyzed the transcriptomic changes before and after the point of no return by inducing metamorphosis of Acropora tenuis with Hym-248, a metamorphosis-inducing neuropeptide. Gene Ontology and pathway enrichment analysis of the 5893 differentially expressed genes revealed that G protein-coupled receptors (GPCRs) were enriched, including GABA receptor and Frizzled gene subfamilies, which showed characteristic temporal expression patterns. The GPCRs were then classified by comparison with those of Homo sapiens, Nematostella vectensis and Platynereis dumerilii. Classification of the differentially expressed genes into modules based on expression patterns showed that some modules with large fluctuations after the point of no return were biased toward functions such as protein metabolism and transport. This result suggests that in precommitted larvae, different types of GPCR genes function to ensure a proper environment, whereas in committed larvae, intracellular protein transport and proteolysis may cause a loss of the reversibility of metamorphosis as a result of cell differentiation.

Background

Planula larvae of corals determine their positions of attachment while dispersing in ocean currents, which consequently influences the limits of their distribution in the ocean [1]. Various factors, such as sensory system stimulation [2], surface structure [3], and chemical substances [4,5,6,7], are involved in the transition from swimming larvae to sessile polyps. When larvae encounter signals that are restricted to small areas of substratum, they cease to swim, rest on the substratum on their aboral end and assume a round shape (Supplementary Fig. 1). If the larvae receive sufficient stimuli from the substratum, they undergo metamorphosis, and stable attachment to the substratum is realized as a result of irreversible cell differentiation processes.

In the genus Acropora, which is used as a model for developmental studies of reef-building corals, cues from crustose coralline algae (CCA) covered by their surface biofilms [8,9,10] have been shown to affect the transition of the planula larvae. Settlement cues from CCA and its surface biofilm are thought to be received by planulae, presumably via sensory neurons [11], and converted into signals to undergo metamorphosis. Alternatively, Hym-248, a GLWamide neuropeptide found in freshwater Hydra [12,13,14,15,16,17], is unique in that, unlike endogenous neuropeptides examined thus far, it can experimentally trigger metamorphosis of acroporids, presumably by mimicking a native internal signaling ligand [12, 15]. As Hym-248 induces metamorphosis without larvae needing to rest on substrata, this exogenous molecule may bypass physiological changes that precede the induction of metamorphosis [18, 19]. Thus, Hym-248-induced metamorphosis is an ideal model system for examining the complex signaling processes and early stages of metamorphosis.

Metamorphosing larvae can revert back to swimming and seek other stopping places even when apparently favorable environmental signals are present (Supplementary Fig. 1), suggesting that sufficient metamorphosis stimuli beyond a ‘threshold’ must be received for the progression of morphogenesis. Similarly, although treating larvae of some coral species (e.g., Acropora tenuis) with Hym-248 causes a morphological change from the swimming to the benthic form, a previous study showed that planula larvae returned to the swimming form if Hym-248 was removed within 4 h at 26 °C [14]. This result suggests that the metamorphosis process can be divided into at least two phases: one is a reversible phase where the level of metamorphosis stimuli is below the threshold, and the other is an irreversible phase where the larval cells are committed to certain cell fates and further progress in metamorphosis after sufficient stimuli are received. Here, we call the critical dividing point the “point of no return” (PNR).

Although gene expression profiles from ESTs and RNA-seq have been analyzed in many studies to search for genes involved in planular metamorphosis [18,19,20,21,22], to our knowledge, no studies have focused on irreversible metamorphosis around the PNR. To better understand the process of metamorphosis of the planula, we focused on the early stages of metamorphosis and identified the metabolic pathways potentially involved in determining reversible metamorphosis.

Materials and methods

Collection and maintenance of Acropora planula larvae

Colonies of A. tenuis were collected with permission from Okinawa Prefecture (#31–65) at Sesoko Island (26° 37′41″N, 127°51′38″E, Okinawa, Japan). Corals were kept in tanks with running natural seawater and under partially shaded natural light at Sesoko Marine Station (University of Ryukyus, Okinawa, Japan). After spawning on June 7, 2020, bundles of aposymbiotic gametes from ten colonies were mixed, and the resulting planula larvae were maintained in plastic bowls at approximately 1000 larvae/L in 10 μm filtered natural seawater (FNSW) or artificial sea water (ASW) using Instant Ocean sea salt (Aquarium systems); the water was exchanged daily.

Induction of metamorphosis and sampling for RNA-seq

At 17 days postfertilization (dpf), larvae were induced to metamorphose in 6-well plates with 200 individuals/3.5 ml ASW supplemented with 16 μM Hym-248 neuropeptide (NH2-EPLPIGLW-CoNH2, synthesized by Scrum Inc., Japan). The timing of metamorphosis induction and the concentration of planula were determined based on previous studies [23, 24], respectively. After the addition of Hym-248, the morphology was observed using a stereomicroscope (OLYMPUS SZX-FOF, OLYMPUS DP 27, Olympus, Japan) at 0 (before treatment), 0.5, 1, 2, 4, 6, and 24 h. Samples for RNA extraction were obtained at 0, 1, and 6 h. For RNA extraction, 20 individuals were soaked in 1.8 ml of RNAlater (Thermo Fisher Scientific, Massachusetts, USA) and stored at − 80 °C until RNA extraction.

RNA extraction and sequencing

RNA was extracted from one single larva preserved in RNAlater (4 biological replicates for each time point, 0, 1 and 6 h). After blotting to remove excess RNAlater, RNA extraction was performed according to the protocol of the Maxwell® RSC Plant RNA Kit (Promega, Wisconsin, USA). The quality and quantity of RNA were verified using an Agilent RNA 6000 Pico Kit on an Agilent Bioanalyzer (Agilent Technologies, California, USA) and a Nanodrop spectrophotometer (Thermo Fisher Scientific), respectively. We confirmed that the Bioanalyzer RIN values ranged from 6 to 9.6 and that all the total RNA spectra were quite similar and suitable for RNA-seq samples. Ten nanograms of total RNA was used for cDNA amplification with the SMART-Seq® v4 Ultra® Low Input RNA Kit for Sequencing (Takara Bio, Japan). Total cDNA samples from planulae were subjected to library preparation using the NEB Next Ultra RNA Library Prep Kit (New England Biolabs, Ipswich, MA, USA) according to the manufacturer’s protocol (NEB #E7530). These mRNA libraries were sequenced on an Illumina NovaSeq6000 (S2 flow cell) in dual flow cell mode with 150-mer paired-end sequences (Filgen Inc., Japan).

Transcriptome analysis

A total of 12 libraries (four larva at each of 0, 1 and 6 h) were obtained and filtered using fastp (with options -q 30 -u 30), and paired output reads were used for analysis. The reads from each library were mapped onto the A. tenuis genome sequences [25] using STAR 2.7 [26] with the default setting. The read count data, transcripts per million (TPM), were calculated using RSEM [27].

To obtain a visual overview of the effect of Hym-248 treatment time on the global gene expression patterns, principal component analysis (PCA) was performed on the TPM values using the function “prcomp” in R. Genes harboring 5 or more mean counts across samples were included in the PCA. Differentially expressed genes (DEGs) among the three groups were detected by a likelihood ratio test using the R package “edgeR” [28] with the count data as input. To define DEGs or genes showing statistically significant differences in expression among the three groups, a false discovery rate (FDR), or q-value, of 0.01 was used as a cutoff.

GO term and KEGG pathway enrichment analysis

InterProScan 5.48 [29] was used to annotate A. tenuis protein sequences enriched with Gene Ontology (GO) terms, resulting in 14,482 genes linked to GO terms. GO term enrichment analysis was performed using the “goseq” package in R, which can compensate for the gene length bias and calculate p values for each GO category using the Wallenius approximation to test over- and underrepresentation among the differentially expressed genes [30]. KAAS [31] was used to annotate the A. tenuis protein sequences for Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis, resulting in 11,283 genes possessing KO annotation. The resulting KO annotations were mapped to the KEGG pathway using the KEGG mapper [32]. Enrichment analysis of the KEGG pathway was performed using WebGestalt [33] based on the obtained KO annotation-KEGG pathway information, with the option ‘weighted set cover’ for redundancy reduction.

Identification and characterization of A. tenuis GPCRs

To search for the GPCR homologs in A. tenuis, a domain search of Pfam [34] was carried out using InterProScan 5.48 [29] for all peptide sequences of A. tenuis, and the genes associated with 45 types of Pfam accessions belonging to the ‘Family A G protein-coupled receptor-like superfamily GPCR_A (CL0192)’ were extracted. The obtained TPM values were normalized by the ‘trimmed mean of M values’ (TMM) method using the R package “edgeR”. The z scores of each normalized TPM value for each gene were calculated in R (Supplementary material 1). The heatmap of the z scores was created using the function “heatmap” in R.

For clustering, the protein sequences of the GPCR homologs of A. tenuis and the protein sequences belonging to the ‘Family A G protein-coupled receptor-like superfamily’ of Homo sapiens, Nematostella vectensis, and Platynereis dumerilii (Platynereis data contained only Family 1 and “others”) were analyzed in CLANS [35] (Supplementary Data 1, 2, 3, 4). All Homo sapiens GPCR protein sequences associated with each Pfam accession (reviewed in (Swiss-Prot), using a manually annotated filter) from the UniProt database were used. Nematostella GPCR protein sequences were identified based on the Pfam domain annotations found in a previous study [36]. Based on the GPCR protein sequences found in a previous study [37], Platynereis GPCR protein sequences were classified into each family by the Pfam domain search in InterProScan 5.48 [29]. CLANS input files were created using the MPI Bioinformatics Toolkit [38] (Scoring Matrix BLOSUM 62, E-value is 1e − 40). CLANS was performed approximately 10,000 times and was drawn with Cluster in 2D.

Analysis of gene modules with time-dependent expression patterns

The TPM data normalized by the TMM method were also used for weighted gene coexpression network analysis (WGCNA [39]) to identify coregulated groups of genes (modules). The signed adjacency matrix was calculated using a soft threshold power of 12 and a minimum module size of 30. Modules were characterized by GO enrichment performed using the “goseq” package in R [30]. The z scores calculated from the normalized TPM values for each gene were used to draw line charts for each module using the R package “ggplot2” [40]. The hub gene of each module was searched using the ChooseTopHubInEachModule of WGCNA. PANTHER annotations output by InterProScan 5.48 [29] were used for analysis.

Results

Morphological and gene expression changes during metamorphosis after the addition of Hym-248

To confirm the effect of Hym-248 on metamorphosis, morphological changes in A. tenuis planula larvae were characterized after the addition of Hym-248 (Fig. 1A). Morphological changes were detectable within 30 min: larvae became rounded and then gradually flattened. At 6 h, they became almost flat, and after 24 h, they showed the form of a juvenile polyp with early skeleton formation. Such morphological changes were in good agreement with those observed in previous studies [41, 42]. After Hym-248 treatment, sampling was carried out at 0 h (before treatment), 1 h (before PNR) and 6 h (after PNR) to obtain transcriptomic data by RNA-seq.

Fig. 1

Global pattern of gene expression in samples used for RNA-seq analysis. A Morphological changes in Acropora tenuis after the addition of Hym-248. The samples at 0 (i.e., before the addition of Hym-248) h, 1 h, and 6 h were subjected to RNA-seq analysis. Scale bar = 100 μm. B PCA of gene expression patterns of each RNA-seq sample. C GO enrichment analysis of DEGs based on three categories: biological process (BP), cellular component (CC), and molecular function (MF)

In this study, RNA-seq produced 536,582,100 reads from 12 samples. The total number of mapped reads of the quadruplicates to the reference sequence was 85,433,460 reads for 0 h (pretreated), 74,442,951 for 1 h, and 101,027,843 for 6 h of treatment with Hym-248. The count values of the genes in the transcriptomic dataset were obtained under three conditions (0, 1 and 6 h). Global gene expression patterns were visualized by PCA, and groupwise separation was confirmed (Fig. 1B).

Here, 5893 genes were identified as candidate DEGs related to the metamorphosis process as their expression levels fluctuated significantly among the three groups. GO and KEGG pathway enrichment analyses were performed. GO enrichment analysis showed that terms associated with GPCRs were enriched in both ‘biological process (BP)’ and ‘molecular function (MF)’ categories (BP, G protein-coupled receptor signaling pathway; MF, G protein-coupled receptor activity) (Fig. 1C). In the KEGG pathway analysis, four pathways, including neuroactive ligand–receptor interactions, were enriched (Table 1). These results suggest that the expression of GPCRs, which are involved in neuronal activity, is significantly altered in the early stages of metamorphosis in A. tenuis.

Table 1 Results of KEGG pathway analysis. All genes: Number of A. tenuis genes in the pathway. DEGs: Number of DEGs in the pathway

Expression pattern changes of A. tenuis GPCR genes during early metamorphosis

To further characterize the GPCR genes involved in early metamorphosis, all homologs were identified from the A. tenuis genome and classified into Families 1, 2, 3, and “others” based on the Pfam domain definitions. Approximately 20 to 45% of the genes in each family (approximately 25% of the total genes) were identified as DEGs without any apparent familywise bias (Table 2). To visualize the patterns in each GPCR family, heatmaps using z scores of gene expression levels were created, and subfamily level groups were defined based on the expression patterns (Fig. 2A, 3A, C, E). Most GPCRs are orphan receptors for which no ligand has been identified in cnidarians. Therefore, to predict the functions of A. tenuis GPCRs from the expression pattern, the GPCR genes were clustered by CLANS, with the H. sapiens GPCR sequences used as references. The results showed that although some genes were homologous, there were many A. tenuis GPCRs with low homology to H. sapiens GPCRs (Fig. 2B, Fig. 3B, D, F).

Table 2 The composition of Acropora tenuis GPCR families. DEGs were detected by multigroup comparison using 0, 1, and 6 h samples

Fig. 2

DEGs in GPCR Family 1 of Acropora tenuis. A Heatmap of gene expression (z score). DEGs were divided into six groups (1a to 1f) according to their expression patterns. B Cluster analysis with Homo sapiens ligand–receptor annotations. Rhombus indicates DEGs of A. tenuis. Rhombus colors correspond to the six DEG groups in Fig. 2A

Fig. 3

DEGs in GPCR Families 2, 3, and others of Acropora tenuis. A, C, E. Heatmap of the gene expression (z score) of Families 2 (A), 3 (C), and others (E). The DEGs were divided into groups according to their expression patterns. B, D, F. CLANS cluster analysis of Families 2 (B), 3 (D), and others (F) with Homo sapiens ligand–receptor annotations. Rhombus colors correspond to the DEG groups in Families 2 (A), 3 (C), and others (E)

DEGs classified as Family 1 were divided into six groups, 1a-1f, based on their expression patterns (Fig. 2A). From the CLANS results, members in each group of Family 1 were distributed throughout the homologous clusters, and no biased GPCR functions were predicted for each group (Fig. 2B). DEGs classified into Family 2 were divided into three groups (2a, 2b, and 2c) based on their expression patterns (Fig. 3A). Regardless of their expression pattern, some members showed high homology to human adhesion GPCRs, but others showed no homology to any human GPCRs (Fig. 3B). DEGs classified as Family 3 were divided into four groups, 3a-3d, based on their expression patterns (Fig. 3C). Group members with increased expression at 1 h or 6 h (groups 3c and 3d) were biased toward metabotropic glutamate receptors (Fig. 3D). In contrast, group 3a, with the highest expression level at 0 and 1 h, was biased toward a cluster of GPCRs acting as gamma aminobutyric acid (GABA) receptors. DEGs classified as “others” were divided into four groups, 4a-4d, based on their expression patterns (Fig. 3E). Group members highly homologous to proteins classified as Frizzled belonged to group 4d; their expression level decreased 1 h after treatment with Hym-248 and recovered at 6 h (Fig. 3F).

Characteristics of gene modules involved in the early phase of A. tenuis metamorphosis

To further characterize the gene expression of all DEGs, we conducted WGCNA to classify them into six modules (MD1–6) (Fig. 4A). In addition, the function of each module was characterized by GO analysis (Fig. 4B, Supplementary Tables 1, 2, 3), and hub genes for each module were detected (Table 3). In modules MD2, 3 and 4, no enriched GO terms were detected. In module MD1, GO terms relevant to transcriptional regulation, GPCRs, and transmembrane transport were enriched (Supplementary Table 1). Many of the MD1 genes (1049 out of 1744 genes) showed the highest expression at 1 h. In module MD5, in which the gene expression level showed a large change at 1 h after the addition of Hym-248 and then recovered, the term for transcriptional regulation was enriched, but potential downstream pathways were not found in the enrichment analyses. In module MD6, GO terms associated with redox reactions, protein metabolism (proteolysis, ubiquitin-dependent protein catabolic process, vesicle−mediated transport, and protein transport) and bioluminescence were enriched. In module MD6, eleven genes associated with the GO term ‘bioluminescence’ were all annotated as ‘green fluorescent protein’ by the Pfam domain search analysis. Most of the MD6 genes were expressed at low levels at 0 h and increased at 6 h (6 h z score > 0 in 1812 genes, < 0 in 1190 genes), and the hub gene of this module was RasGAP-like protein 1.

Fig. 4

Modular analysis of the expression patterns of the DEGs by WGCNA in the three-group comparison. A The dendrogram of the modules and the expression pattern of the genes in each module. The number in parentheses indicates the number of genes clustered into each module. The average values of the z scores of gene expression levels over time are shown by a line plot. B GO enrichment analysis for each module. Only BP categories are shown. MD2, 3, and 4 have no enriched terms

Table 3 Hub genes of each module detected by WGCNA

Discussion

It is critical to delineate the coral metamorphosis process at the molecular level to further understand how larval dispersal can be regulated, which has a high impact on the diversity and stability of coral reef ecosystems in tropical and subtropical oceans. Previous studies have identified many metamorphosis-related candidate genes by analyzing changes in gene expression at many developmental stages [20] before and after metamorphosis [18, 19] and focusing on metamorphosis competence [21]. In this study, by analyzing the gene expression changes at a fine time resolution before and after the PNR, candidate molecular mechanisms involved in the reversibility of metamorphosis were identified. These results suggest that alterations in environmental signal perception through GPCRs and proteolytic regulation of cellular components can lead to drastic changes that make the metamorphosis process irreversible. These results also show that the gene expression of green fluorescence proteins may be under strict control after the PNR, implying a link to understudied physiological and ecological roles of fluorescence in A. tenuis at early developmental stages [43, 44].

A. tenuis differentially regulates GPCR transcripts during early metamorphosis

In this study, a group of genes (GPCRs) were found to have altered gene expression patterns at the early stage of metamorphosis, and the GPCRs with altered expression accounted for approximately 25% of the GPCR homologs of A. tenuis. Previous studies have suggested that GPCRs are involved in the metamorphosis of planula larvae in A. millepora [21]. This study suggests that the involvement of GPCRs in metamorphosis is conserved in the genus Acropora and provides the first comprehensive analysis of the expression patterns of GPCRs at the gene family and subfamily levels.

GABA receptors and Frizzled GPCRs showed dynamic expression patterns in this study. The expression level of the GABA receptor decreased within 1 h after the addition of Hym-248 (Fig. 3C, D). It has been shown that the addition of GABA_B antagonists to seawater does not significantly affect metamorphosis itself [21], and GABA receptors are expected to play a role in searching for a substratum before metamorphosis [18] in A. millepora. These results suggest that GABA receptors may be involved in the search behavior and that this ability is lost as metamorphosis proceeds because the search behavior of the planula is no longer necessarily activated. Since GABA_B receptors are known to stimulate cell differentiation [45], some GABA receptors may be activated at a later developmental stage.

The expression of Frizzled was suppressed 1 h after the addition of Hym-248 (Fig. 3E, F). Frizzled is a receptor involved in the reception of Wnt signaling, and binding of Wnt to Frizzled initiates transcriptional regulation via β-catenin. In Hydra and Hydractinia, Wnt is expressed in the oral tip and functions as a head organizer [46, 47]. In corals, Wnt may also drive the differentiation of the oral part of the developing polyp during metamorphosis, and changes in Frizzled expression may affect this process. In addition, Ras GTPase Activating Protein (RasGAP)-like protein 1 was detected as a hub gene of MD6, as classified by WGCNA (Table 3), suggesting that Frizzled-Wnt signaling may be mediated by RasGAP after the PNR.

Protein-level metabolism is involved in the regulation of metamorphosis

Network analyses identified the gene modules responsible for metabolic pathways potentially involved in the reversibility of metamorphosis (Fig. 5). The results revealed that some metabolic pathways were important throughout the early metamorphosis process, while others made a switch in regulation around the PNR. Transcriptional regulation and receptor signaling were identified as functional characteristics of the precommitted phase, but no specific downstream functions were predicted (see MD1 and, 5 in Fig. 4). In addition to transcriptional regulation, the committed phase after the PNR was characterized by GO terms related to the ubiquitin–proteasome system associated with the genes encoding 26S proteasome and ubiquitin-related enzyme subunits (MD6 in Fig. 4, Supplementary Table 3). This result is consistent with the fact that metamorphosis becomes irreversible after the PNR: compared to the transcription-level regulation detected before the PNR, it would be far more difficult to restore a system to the original form after it has been modified or degraded by proteolytic regulation.

Fig. 5

Metabolic pathways expected to be important in the early metamorphosis induced by Hym-248. A model illustrating that in early metamorphosis (between 0 and 6 h; regardless of the PNR) in A. tenuis, GPCR expression repertoire replacement, transcription, and transmembrane transport play important roles. After the PNR, the onset of protein transport and proteolysis functions become more relevant to the loss of metamorphosis reversibility. The colors of bands (blue and purple) correspond to the colors of the modules (MD1 and MD6, respectively) where the relevant GO terms are enriched, as shown in Fig. 4

Regulation of proteolysis is important in bacterial development [48] and may play an important role in the progression of coral metamorphosis. During the metamorphosis of Acropora, tissue differentiation proceeds rapidly (e.g., formation of the solid, single-layer endoderm and oral organs, including tentacles, mesenteries, and pharynx) [42]. Apoptosis and cell division are known to occur actively after the PNR in the reef building coral Stylophora [49]. This result suggests that active tissue reconstitution, including the degradation of existing tissues and the differentiation of new tissues, is preceded by the regulation of the expression of genes related to protein degradation. Although the targets of proteasomal degradation remain unclear, proteolytic regulation may be involved in the downstream differentiation cascades, leading to a complete loss of reversible metamorphosis. Genetic engineering of coral larvae [50] to suppress specific protein degradation processes is an ideal platform to substantiate this hypothesis.

Conclusions

This study identified key metabolic pathways in the initial stages, especially before and after the PNR, of the metamorphosis of the reef-building coral A. tenuis. Comprehensive gene expression pattern analysis suggested that GPCRs with various functions can play important roles in the different phases of coral metamorphosis. In addition, modular analysis of the gene expression fluctuations identified molecular mechanisms associated with the loss of the reversibility of metamorphosis, indicating the onset of proteolytic regulation in the early phases of metamorphosis.