In light of the chronic stress and mass mortality reef-building corals face under climate change, it is critical to understand the processes driving reef persistence and replenishment, including coral reproduction and development. Here, we quantified gene expression and sensitivity to ocean acidification across a set of developmental stages in the rice coral, Montipora capitata. Embryos and swimming larvae were exposed to pH treatments of 7.8 (ambient), 7.6 (low) and 7.3 (extreme low) from fertilization to 9 days post-fertilization. Embryo and larval volume, and stage-specific gene expression were compared between treatments to determine the effects of acidified seawater on early development. Although there was no measurable size differentiation between pH treatments at the fertilized egg and prawn chip (9 h post-fertilization) stages, early gastrulae and larvae raised in reduced pH treatments were significantly smaller than those raised in ambient seawater, suggesting an energetic cost to developing under low pH. However, no differentially expressed genes were found until the swimming larval stage. Notably, gene expression patterns of larvae developing at pH 7.8 and pH 7.3 were more similar than those of larvae developing at pH 7.6. Larvae from pH 7.6 showed upregulation of genes involved in cell division, regulation of transcription, lipid metabolism and response to oxidative stress in comparison to the other two treatments. Although low pH appears to increase energetic demands and trigger oxidative stress in larvae, the developmental process is robust to this at a molecular level, with the swimming larval stage reached in all pH treatments.

As anthropogenic carbon dioxide (CO2) levels continue to rise and contribute to global climate change (Bindoff et al., 2019), approximately one-third of atmospheric CO2 is also being absorbed by the ocean (Sabine et al., 2004). Increasing concentrations of CO2 in the ocean leads to ocean acidification (OA) through a shift in the chemical equilibrium of carbon in seawater (Sabine et al., 2004) that reduces ocean pH and the availability of carbonate ions (Caldeira and Wickett, 2003) for calcium carbonate formation. The absorption of atmospheric CO2 by the ocean has already driven a 0.1 unit decrease in pH in the ocean's surface layer (where gas exchange occurs between the atmosphere and the ocean), which is projected to decrease by an additional 0.036–0.042 or 0.287–0.29 units by the year 2100 under the Intergovernmental Panel on Climate Change 2019 high and low mitigation scenarios, Representative Concentration Pathways 2.6 and 8.5, respectively (Bindoff et al., 2019).

A growing body of studies has reported negative effects of reduced seawater pH on marine organisms. OA can negatively impact biomineralizing marine taxa that build calcium carbonate shells and skeletons (Doney et al., 2009; Kroeker et al., 2010) because of the lower concentration of carbonate ions available for calcium carbonate formation (Caldeira and Wickett, 2003). Coral reefs are one marine system threatened by OA as reef corals are ecosystem engineers, building calcium carbonate skeletons that generate three-dimensional complexity and thus create habitat for an entire ecosystem. Although studies investigating the effects of OA on marine calcifiers have primarily examined its impact on calcification, less is known about the impact of OA on embryonic development. Sensitivity to environmental stressors during early development can present a bottleneck for species persistence in a changing ocean (Byrne, 2012; Przeslawski et al., 2015). These early life stages are crucial for the resilience of marine organisms given their role in repopulating and replenishing existing populations (Adjeroud et al., 2017; Ritson-Williams and Arnold, 2009).

OA can affect multiple early developmental processes in corals and other marine invertebrates. Specifically, lower pH impacts sperm motility (Albright, 2011; Havenhand et al., 2008; Morita et al., 2010), fertilization success (Albright, 2011; Albright and Mason, 2013; Albright et al., 2010; Byrne, 2012; Han et al., 2021; Havenhand et al., 2008; Kurihara and Shirayama, 2004), larval settlement (Albright, 2011; Albright and Langdon, 2011; Albright et al., 2010) and post-settlement growth (Albright and Langdon, 2011; Albright et al., 2008, 2010; Scucchia et al., 2021a). In corals, decreased fertilization efficiency has been observed under low pH, which may be explained by decreased sperm flagellar motility and increased polyspermy (Albright and Mason, 2013; Albright et al., 2010). Low pH is also correlated with metabolic depression in many coral species (Moya et al., 2012, 2015; Nakamura et al., 2011; Putnam et al., 2013), which may explain observations of delayed development, smaller embryo/larval volume and reduced post-settlement growth in corals developing in acidified conditions (Albright and Langdon, 2011; Nakamura et al., 2011; Rivest and Hofmann, 2014; Yuan et al., 2018). The increased energetic costs of maintaining homeostatic processes under disrupted acid–base balances and ionic gradients can additionally slow development and reduce growth of marine invertebrate larvae and recruits (Albright and Langdon, 2011; Byrne, 2012; Przeslawski et al., 2015; Timmins-Schiffman et al., 2013). For example, low pH is correlated with delayed cleavage (Foo et al., 2012), smaller larvae (Kelly et al., 2013) and depressed metabolism (O'Donnell et al., 2010) during the early development of sea urchins. Developmental delays (Barton et al., 2012; Timmins-Schiffman et al., 2013) and decreases in fertilization rate (Barros et al., 2013; Parker et al., 2009, 2010), hatching success (Barros et al., 2013; Gazeau et al., 2010; Van Colen et al., 2012), larval survival (Barros et al., 2013; Van Colen et al., 2012) and larval size (Barros et al., 2013; Barton et al., 2012; Van Colen et al., 2012) have also been observed in bivalves under low pH. These compounding effects of low pH and co-occurring stressors can present a major life history bottleneck, leading to downstream effects on the resilience of marine invertebrate species (Byrne, 2012; Byrne and Przeslawski, 2013).

Most studies examining the effects of OA on coral development thus far have focused primarily on fertilization (Albright and Mason, 2013; Albright et al., 2010), settlement (Albright, 2011; Albright and Langdon, 2011; Albright et al., 2010) and post-settlement growth (Albright and Langdon, 2011; Albright et al., 2008, 2010; Scucchia et al., 2021a), with less of a focus on embryonic development. In many species, embryogenesis represents a critical window in development during which organisms are particularly sensitive to environmental fluctuations outside of the expected range (Byrne, 2012; Byrne et al., 2009; Dahlke et al., 2020; Ericson et al., 2010, 2012; Foo et al., 2012; Przeslawski et al., 2015). Several studies have implicated the loss of maternal defenses during the maternal-to-zygotic transition (MZT) as playing a role in increased environmental sensitivity during this period (Dahlke et al., 2017; Ericson et al., 2012; Foo et al., 2012; Przeslawski et al., 2015). The MZT describes the developmental reprogramming of gene expression from maternal to zygotic regulation (Tadros and Lipshitz, 2009; Vastenhouw et al., 2019). During the MZT, the maternally provided gene products that regulate early development and protect embryos from expected stressors are degraded as zygotic transcription activates and intensifies (Tadros and Lipshitz, 2009; Vastenhouw et al., 2019). Higher mortality in response to environmental perturbations has been observed during this transition in fish gastrulae (Dahlke et al., 2017, 2020), sea urchin blastulae/gastrulae (Ericson et al., 2012; Foo et al., 2012) and several other marine invertebrates (Przeslawski et al., 2015).

Maternal mRNA provisioning and the timing of the MZT in reef-building corals, thus far, have only been investigated in the Hawaiian rice coral, Montipora capitata (Chille et al., 2021). Under ambient conditions, the maternal mRNA complement in M. capitata eggs primarily provide housekeeping functions, including cell cycle, biosynthesis, transcription, signaling and protein processing (Chille et al., 2021; Van Etten et al., 2020). However, DNA repair mechanisms, which may protect early embryos from ultraviolet (UV) radiation and reactive oxygen species (ROS), are also over-represented in the M. capitata maternal mRNA complement of eggs (Chille et al., 2021). As of yet, the maternal mRNA complement in reef-building corals has not been assessed under stressful conditions. However, based on timing of the MZT in M. capitata (Chille et al., 2021), and general developmental gene expression studies in Acropora millepora (Strader et al., 2018), Acropora digitifera (Reyes-Bermudez et al., 2016) and Nematostella vectensis (Helm et al., 2013), coral embryos may be most sensitive to stress during early gastrulation (between 9 and 14 h post-fertilization; hpf), when a large portion of maternal transcripts are degraded and zygotic transcription begins.

Previous transcriptomic examinations of pH stress during coral embryonic development have been completed on aposymbiotic Acropora spp. (Strader et al., 2018) and on post-settlement stages after the completion of the MZT (Moya et al., 2012; Scucchia et al., 2021a). Here, we used M. capitata, a hermaphroditic broadcast-spawning coral in the family Acroporidae with vertical symbiont transmission (Heyward, 1986; Padilla-Gamiño and Gates, 2012), to examine the embryonic sensitivity to OA in the context of the MZT. In this study, we examined the physiological and transcriptomic outcomes of early coral development in five life stages from single-celled embryo to swimming larva under pH 7.8 (ambient), and two low pH conditions (7.6 and 7.3). This time series allows us to determine sensitivity to, and phenotypic consequences of, OA during and after the MZT.

Specimen collection

Montipora capitata Dana 1846 egg–sperm bundles were collected by scooping from the water's surface under the Department of Land & Natural Resources Special Activity Permit 2018-50 from the reef adjacent to the Hawaiʻi Institute of Marine Biology, Coconut Island, Kāneʻohe Bay, Oʻahu Hawaiʻi (21°25'58.0″N, 157°47'24.9″W), at ∼21:00 on 13 June 2018. The gametes collected therefore likely represent a mix of genotypes from an unknown number of parents, characteristic of a typical mass spawning event for Kāneʻohe Bay Reef 1, Moku o Loʻe (numbered as designated by the State of Hawaiʻi, Division of Aquatic Resources). After collection, egg–sperm bundles were allowed to break apart. The eggs are positively buoyant and the sperm negatively buoyant. Thus, aliquots of concentrated eggs were separated from the sperm via pipetting into a new tube and were then rinsed three times with 0.2 μmol l−1 filtered seawater, snap-frozen and stored at −80°C in ∼300 µl aliquots. The remainder of the live eggs were split evenly between nine replicate 1.6 l Brine Shrimp Direct hatchery conicals (7.5 ml of concentrated eggs and 25 ml of homogeneous sperm mixture per 1.6 l container; Ogden, UT, USA).

Experimental set-up

Fertilization and early development occurred in the nine 1.6 l hatchery conicals set to pH 7.8 (ambient), pH 7.6 and pH 7.3 (n=3 conical tanks treatment; Fig. 1A). At 22 hpf, the embryos were transferred to replicate 118 ml bins with 152 µm plankton mesh bottoms floating in 74 l glass tanks (n=3 bins per tank, 2 replicate tanks per treatment; Fig. 1B). These pH treatments were chosen to reflect current ambient conditions occurring in Kāneʻohe Bay, Hawaiʻi (Drupp et al., 2011), and projected future conditions for coastal embayments (Jury et al., 2013; Shaw et al., 2016), which, owing to benthic feedbacks and longer residence time, tend to have lower mean pH and higher natural diel fluctuations. All treatments were controlled by a pH-stat feedback system in a header tank, with one flow-through header tank per treatment. Header tank conditions were continuously monitored using a Neptune Apex Reef Aquarium set-up (Energy Bar 832; Controller Base Unit; PM1 pH/ORP Probe Module; Extended Life Temperature Probe; pH Probe), and 99.99% food-grade CO2 was injected into the header tanks with gas flow solenoids (Milwaukee, MA, USA) using inline Venturi injectors (Forfuture-go G1/2 Garden Irrigation Device Venturi Fertilizer Injector) in a pressure-driven line connected to a circulating pump (Pondmaster Pond-mag Magnetic Drive Water Pump Model 5) in each header tank (Fig. 1A). Water flow was delivered to each replicate conical by 1.9 l h−1 Pressure Compensating Drippers, for a turnover rate of once every ∼50 min. The average light intensity (mean±s.e.m.=96±8.823 µmol m−2 s−1; n=24), was measured using a handheld PAR sensor (Underwater Quantum Flux Apogee instruments, model MQ-510; accuracy=±4%).

Fig. 1.

Experimental summary. (A) Experimental design and (B) sampling timeline, including example photographs of each life stage. Early life stages were exposed to treatment immediately after fertilization. Analyses to assess the morphological effects of low pH were performed on unfertilized eggs (0 hpf; n=50 ambient), fertilized eggs (1 hpf, n=75 replicates per treatment), cleavage (4 hpf; n=3 replicate samples per treatment), early gastrulae (14 hpf, n=75 ambient and extreme low, n=50 low) and swimming larvae (planula, 9 dpf, n=50 replicates per treatment). Gene expression analyses were performed on fertilized eggs (1 hpf), cleavage (4 hpf), prawn chip (9 hpf), early gastrulae (14 hpf) and swimming larvae (planula, 9 dpf) exposed to pH 7.8 (ambient), pH 7.6 (low) and pH 7.3 (extreme low) [n=3 replicates per life stage per treatment, except where n=2 replicate samples per treatment (ambient 1 hpf, extreme low 4 hpf, extreme low 9 hpf, and extreme low 14 hpf)]. Embryos were moved from treatment conicals to tanks at 22 hpf.

Fig. 1.

Experimental summary. (A) Experimental design and (B) sampling timeline, including example photographs of each life stage. Early life stages were exposed to treatment immediately after fertilization. Analyses to assess the morphological effects of low pH were performed on unfertilized eggs (0 hpf; n=50 ambient), fertilized eggs (1 hpf, n=75 replicates per treatment), cleavage (4 hpf; n=3 replicate samples per treatment), early gastrulae (14 hpf, n=75 ambient and extreme low, n=50 low) and swimming larvae (planula, 9 dpf, n=50 replicates per treatment). Gene expression analyses were performed on fertilized eggs (1 hpf), cleavage (4 hpf), prawn chip (9 hpf), early gastrulae (14 hpf) and swimming larvae (planula, 9 dpf) exposed to pH 7.8 (ambient), pH 7.6 (low) and pH 7.3 (extreme low) [n=3 replicates per life stage per treatment, except where n=2 replicate samples per treatment (ambient 1 hpf, extreme low 4 hpf, extreme low 9 hpf, and extreme low 14 hpf)]. Embryos were moved from treatment conicals to tanks at 22 hpf.

Temperature, salinity and pH were measured three times daily in each treatment tank using a handheld digital thermometer (Fisherbrand Traceable Platinum Ultra-Accurate Digital Thermometer, accuracy=±0.05°C, resolution=0.001°C) and a portable multiparameter meter (Thermo Fisher Scientific Orion Star A series A325; accuracy=±0.2 mV, 0.5% of psu reading, resolution=0.1 mV, 0.01 psu) with pH and conductivity probes (Mettler Toledo InLab Expert Pro pH probe 51343101; Orion DuraProbe 4-Electrode Conductivity Cell Model 013010MD), respectively (Table 1). Temperature and pH probes were calibrated using Tris (Dickson Laboratory Tris Batch 27, Bottles 70, 75, 167, 245 and 277) standard calibrations. Total alkalinity (µmol kg−1 seawater) was measured from water samples daily using an automated titrator (Mettler Toledo T50) and hydrochloric acid titrant (Dickson Laboratory Titrant A3), with certified reference material (Dickson Laboratory CO2 CRM Batch 149) used as a quality control standard. Carbonate chemistry parameters (pH, CO2, pCO2, HCO3, CO3, dissolved inorganic carbon, total alkalinity and aragonite saturation) were calculated from total alkalinity measurements using seacarb v3.0.11 (https://CRAN.R-project.org/package=seacarb) in RStudio (Table 1; Table S1) using Kf (equilibrium constant of hydrogen fluoride, mol kg−1; Kf=[H+][F]/[HF]) from Perez and Fraga (1987), K1 (first dissociation constant of carbonic acid, mol kg−1; K1=[H+][HCO3]/[CO2]) and K2 (second dissociation constant of carbonic acid, mol kg−1; K2=[H+][CO32−]/[HCO3]) from Lueker et al. (2000) and Ks (stability constant of hydrogen sulfate, mol kg−1; Ks=[H+]F[SO42−]/[HSO4]) from Dickson (1990). Seawater pH was calculated on the total scale (Dickson et al., 2007; SOP6a), using the gas constant (8.31447215 J K−1 mol−1) and the Faraday constant (96485.339924 coulombs mol−1).

Table 1.

Descriptive statistics for temperature, salinity and carbonate chemistry across the 9 days of exposure (means±s.e.m.)

Descriptive statistics for temperature, salinity and carbonate chemistry across the 9 days of exposure (means±s.e.m.)
Descriptive statistics for temperature, salinity and carbonate chemistry across the 9 days of exposure (means±s.e.m.)

Sample collection and preservation

Physiological and molecular samples were taken immediately upon breakup of egg–sperm bundles and at five time points following transfer to experimental conditions (Fig. 1). Sampling times represent visually distinct developmental stages including fertilized embryo (1 hpf), cleavage (4 hpf), prawn chip (9 hpf), early gastrulation (14 hpf) and swimming larva (planula, 9 days post-fertilization, dpf). Physiological samples (∼100–200 embryos or larvae per replicate, n=3 replicate samples per treatment, for a total of 45 samples) were preserved in 500 µl 25% zinc buffered formalin in 0.2 µm filtered seawater, and stored at 4°C. Molecular samples [∼300 µl of concentrated embryos or larvae and surrounding seawater per replicate, n=3 replicate samples per treatment, except where n=2 replicate samples per treatment (ambient 1 hpf, extreme low 4 hpf, extreme low 9 hpf and extreme low 14 hpf), for a total of 41 samples] were snap-frozen in liquid nitrogen and stored at −80°C. In total, aliquots were estimated to contain 4500 to 6500 embryos or larvae based on calculations using average volume at each life stage (see Results, ‘Embryo development and growth’).

Embryo development and growth

Stereomicroscope images of fixed samples were analyzed for blastomere number at 4 hpf (cleavage stage), and for maximum length and perpendicular width of the unfertilized eggs, gastrulae and planulae using ImageJ (Schneider et al., 2012). Cell division stage was assessed in cleavage-stage samples with 82 to 282 embryos per sample (n=3 replicate samples per treatment). Cleavage-stage embryos were classified as ‘one cell’, ‘two cells’, ‘four cells’ or ‘greater than four cells’. The percentage of embryos at each stage of cell division was then calculated for each sample to examine potential variation in developmental timing between treatments. The volume of unfertilized eggs (n=50, ambient), fertilized embryos (n=75 replicates per treatment), early gastrulae (n=75 ambient and extreme low, n=50 low) and planulae (n=50 replicates per treatment) was calculated using the equation for volume of an oblate ellipse (V=4/3πab2, in which V is equal to volume, a is the half of maximum width and b is the half maximum length) to examine differences in size between treatments through early development. The volume of cleavage and prawn chip stages was not calculated owing to their irregular shape.

All statistical analyses for morphological data were performed in RStudio (v1.3.959; https://www.rstudio.com/), using R version 4.0.2 (https://www.r-project.org/). A beta regression model (formula=proportion∼cleavage stage+treatment+cleavage stage:treatment) was used to analyze differences in the proportion of cells at each cleavage stage and treatment using the betareg R package (v3.1-4; Cribari-Neto and Zeileis, 2010). Differences between cleavage stages and treatments were then computed using the joint_tests function from the emmeans R package (v1.4.4; https://CRAN.R-project.org/package=emmeans), which effectively runs the beta regression model as a type III ANOVA. Additionally, a one-way nested ANOVA analysis tested the effect of tank on embryo and planula size wherein tank ID was nested within treatment. After determining that tank effects were non-significant (P>0.05), one-way ANOVA analysis was used to test for the effect of treatment on fertilized embryo, gastrula and planula volume. Post hoc Tukey HSD tests were conducted when the effect of treatment was significant (P<0.05). Data were visually examined for normal distribution and equal variance. The dependent variable was square-root transformed for the gastrulation and planula life stages in order to meet statistical assumptions prior to ANOVA analysis. Data points for these life stages were back-transformed for visualization in Fig. 2.

Fig. 2.

Developmental age and size of Montipora capitata. (A) Proportion of embryos at each cleavage stage at 4 hpf (n=3 replicate samples per treatment), and volume (mm3) of (B) unfertilized eggs (0 hpf; n=50 ambient), (C) fertilized eggs (1 hpf, n=75 replicates per treatment), (E) early gastrulae (14 hpf, n=75 ambient and extreme low, n=50 low) and (E) swimming larvae (planula, 9 dpf, n=50 replicates per treatment). Prawn chip stages were not analyzed given the fold in the morphology, which generates added assumptions for volume calculations (see Fig. 1 prawn chip picture).

Fig. 2.

Developmental age and size of Montipora capitata. (A) Proportion of embryos at each cleavage stage at 4 hpf (n=3 replicate samples per treatment), and volume (mm3) of (B) unfertilized eggs (0 hpf; n=50 ambient), (C) fertilized eggs (1 hpf, n=75 replicates per treatment), (E) early gastrulae (14 hpf, n=75 ambient and extreme low, n=50 low) and (E) swimming larvae (planula, 9 dpf, n=50 replicates per treatment). Prawn chip stages were not analyzed given the fold in the morphology, which generates added assumptions for volume calculations (see Fig. 1 prawn chip picture).

RNA extraction, sequencing and processing

To extract RNA, molecular samples were digested at 55°C in 300 µl of DNA/RNA Shield (Zymo Research, Irvine, CA, USA) for 2–3.5 h and centrifuged at 2200 RCF for 1 min to separate the remaining solids. Total RNA was extracted from each supernatant with the Zymo Quick-DNA/RNA™ Miniprep Plus Kit (Zymo Research) following the manufacturer's protocol for tissue samples. RNA was quantified (ng μl−1) with a Thermo Fisher Qubit Fluorometer and quality was assessed using the Eukaryotic RNA analysis on the Agilent TapeStation 4200 System. Total RNA samples were then sent to Genewiz (South Plainfield, NJ, USA) for library preparation and sequencing. cDNA libraries were constructed following the TruSeq Stranded mRNA Sample Preparation protocol (Illumina, San Diego, CA, USA) and sequenced on a HiSeq instrument targeting 15 million reads per sample. Quality trimming, alignment and assembly were conducted as described in Chille et al. (2021). In short, sequences were filtered for quality by applying a 5 bp sliding window to remove reads with an average quality score less than 20, and retain sequences with quality scores greater than or equal to 20 in at least 90% of bases and sequence length greater than or equal to 100 bases. Trimmed reads were then aligned to an M. capitata reference genome assembly (Shumaker et al., 2019) using HISAT2 in the stranded paired-end mode (Kim et al., 2019), and assembled with StringTie in the stranded setting (v2.1; Pertea et al., 2015). Finally, the StringTie prepDE python script (Pertea et al., 2015) was used to generate a single gene count matrix containing the gene count data for each sample (Pertea et al., 2015).

Host gene expression analyses

Experiment-wide patterns in gene expression were assessed with a principal components analysis (PCA). Prior to analysis, low-coverage genes were removed from the gene count matrix using the Genefilter pOverA function (v1.70.0; http://www.bioconductor.org/packages/release/bioc/html/genefilter.html). Here, genes that had a count of 10 or less (A<10, where A=counts) in fewer than two out of 41 samples (p=0.0488, where p=proportion of samples) were removed from the gene count matrix. A proportion of p=0.0488 was chosen under the assumption that owing to the MZT, a gene count of zero likely represents a true zero because some genes in M. capitata embryos may not be activated until the onset of gastrulation or later (Chille et al., 2021). Then, counts were transformed with the DESeq2 (v1.28.1) variance stabilizing transformation (VST) function (Love et al., 2014a) after confirming that all size factors were less than 4. Finally, the VST-transformed counts were used for PCA to visualize expression results. The DESeq2 plotPCA function, with default parameters, was used to generate a PCA of all samples based on sample-to-sample distance. All gene expression analyses were performed in RStudio (v1.3.959), using R version 4.0.2.

PCA and differential gene expression analysis with DESeq2 were used to assess within-stage differences in expression owing to treatment. As described above, genes with low overall expression were filtered using the Genefilter pOverA filter function (Gentleman et al., 2020). After testing a range of proportions from p=0.33 to p=0.9, a more conservative value of p=0.875 was chosen to prevent bias owing to high duplication levels in the first four life stages, where zeros likely represent true zeros compared with the planula samples, which had lower duplication. As such, we retained genes with greater than 10 counts in at least seven out of eight samples (pOverA, p=0.875, A=10) for further analysis. The resulting genes were then VST-transformed using the DESeq2 VST-transformation function (Love et al., 2014b). PCA was conducted with the DESeq2 plotPCA function, with default parameters (Love et al., 2014b), using VST-transformed counts to calculate sample-to-sample distance. Differential expression analysis between pH treatments was then performed for each life stage using the filtered counts (pOverA, p=0.875, A=10). Pairwise differences in gene expression between treatments were estimated using the Wald model in DESeq2 (Love et al., 2014b). As DESeq2 results showed zero differentially expressed genes (DEGs; Padj<0.05) in all treatment comparisons prior to the planula stage, only genes expressed at the planula stage were used for further analysis. To visualize the DEGs that were shared and distinct for each treatment comparison within the planula stage, the results were summarized as a Venn diagram (Fig. 4A) using the VennDiagram R package v1.6.20 (Chen and Boutros, 2011). Overall expression patterns among the planula DEG set were visualized as a heatmap using the ComplexHeatmap package v2.5.4 (Gu et al., 2016), with clustering based on expression relative to the row mean in the VST-transformed count matrix. Clusters in the heatmap are based on k-means clustering using the R package NbClust v3.0 (Malika et al., 2014) to calculate the optimal number of clusters using 30 indices.

Host functional enrichment analyses

Gene Ontology (GO) and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analyses were performed on the planula DEG set to determine which biological processes, molecular functions and pathways are over-represented in each DEG k-means cluster. Firstly, functional annotation of the M. capitata genome was retrieved from Chille et al. (2021). KEGG annotation was additionally completed using KofamScan (v1.3.0; Aramaki et al., 2020) to retrieve KEGG terms in the eukaryote profiles database (v. 23 April 2021). GO and KEGG enrichment analyses were performed using the default ‘Wallenius’ model on the DEGs with log2FoldChange>|1| in each k-means cluster using the R package Goseq (v1.40.0; Young et al., 2010). Subsequently, slim categories for the enriched GO terms were obtained using the goSlim function in the R package GSEABase, using the generic goSlim obo database as reference (http://current.geneontology.org/ontology/subsets/goslim_generic.obo; v1.2; Ashburner et al., 2000). Finally, enriched biological process and molecular function terms with numInCat>5 and their associated slim categories were visualized as a heatmap using the ggplot2 geom_tile function (Wickham, 2011).

Symbiont transcriptome assembly and gene expression analysis

To assess correlation between pH treatment and gene expression of Symbiodinaceae within the planula host, the symbiont transcriptomes were extracted from a de novo-assembled holobiont sequence library and analysed for differential expression between treatments. The holobiont transcriptome was assembled de novo using Trinity v2.12.0 (Grabherr et al., 2011) with the stranded ‘RF’ option. Following quality trimming (see RNA extraction, sequencing and processing), the nine planula cDNA libraries (n=3 replicate samples per treatment) were used to construct the holobiont transcriptome. The three adult M. capitata libraries from Chille et al. (2021) were additionally used for de novo assembly in order to establish a symbiont-rich library. Following assembly, completeness was assessed using BUSCO v4.0.6 (Simão et al., 2015) to detect the presence of single-copy orthologs universal to Alveolata (for the symbionts) and metazoans (for the host). The host and symbiont transcriptomes were then identified from the holobiont assembly using psytrans (https://github.com/sylvainforet/psytrans) with Cladocopium C1 (Liew et al., 2016; Voolstra et al., 2015) and M. capitata (Shumaker et al., 2019) predicted proteins as references. Following host–symbiont separation, the symbiont reads were used as a reference file for the Trinity align_assemble perl script (Grabherr et al., 2011) using RSEM v1.3.3 (Li and Dewey, 2011) for abundance estimation and Bowtie v-1.2.1.1 as the alignment method (Langmead, 2010). Finally, the symbiont gene counts were obtained using the Trinity abundance_estimates_to_matrix perl script (Grabherr et al., 2011). Differential gene expression of the planula endosymbionts was carried out using the methods described for the host early life stages (see Host gene expression analyses). Finally, the identities of the symbiont DEGs were determined through a local alignment search of the corresponding sequences against the NCBI non-redundant database using Blastx v2.11.0 (Altschul et al., 1997) with default parameters.

Embryo development and growth

Immediately after release, prior to treatment exposure, eggs were 0.040±0.023 mm3 (n=50; Fig. 2A). Following treatment exposure, the proportion of embryos at each cleavage stage at 4 hpf and embryo size in the fertilized egg (n=75 replicates per treatment; 1 hpf), early gastrula (n=75 replicates per treatment ambient and extreme low, n=50 low; 14 hpf) and planula (n=50 replicates per treatment; 9 dpf) samples was assessed to characterize morphological effects of developing in ambient (7.8), low (7.6) and extreme low pH (7.3) conditions. At the cleavage stage, across all treatments combined, undivided cells were on average more frequent (44.77±21.49%; mean±s.d.) than those at the two-cell (26.51±12.63%), four-cell (18.37±11.23%) or greater-than-four-cell stage (10.34±13.40%). Beta regression analysis of the proportion of cells at each cleavage stage and treatment (proportion∼cleavage stage+treatment+cleavage stage:treatment; Fig. 2B) shows that, although there were significant differences (P<0.0001, F=14.144) in proportions of two-, four- and greater-than-four-cell stages overall, the interaction between cleavage stage and treatment were not significant (n=82–282; P=0.0765, F=1.902) with the given sample size. Additionally, one-way ANOVA showed there were no significant differences between pH treatments in embryo size in the fertilized eggs (P=0.2560, F=1.082; Fig. 2C). However, at both the early gastrula (P=0.0003, F=8.535; Fig. 2D) and planula stages (P=0.0001, F=9.716; Fig. 2E), embryos developing under ambient pH conditions were significantly larger than those developing under pH 7.6 and 7.3 conditions.

RNA sequencing, quality control, and mapping

TruSeq Illumina Stranded Paired-End (PE) sequencing resulted in 935,000,000 PE (19,392,156 mean reads per sample), with 694,800,000 PE reads (14,241,176 mean reads per sample) remaining after quality filtering and adapter removal. On average, 79.28±1.95% (mean±s.e.m. reads per sample) aligned to the reference M. capitata genome (Shumaker et al., 2019). Assembly quality assessment with GFFcompare showed that 63,227 genes queried from the mapped GTF files were matched to the 63,227 genes in the reference genome (Shumaker et al., 2019), with 40,883 matching intron chains, 63,218 matching loci and zero missing elements.

Host gene expression patterns

After filtering to remove genes with fewer than 10 counts in ∼4.88% of samples (i.e. the proportion representing each treatment in a given life stage), 30,271 genes remained to characterize global gene expression patterns. A PCA (Fig. 3A) showed that sample-to-sample variation was primarily related to differences in life stage, with PC1 (83%) and PC2 (12%) explaining the majority of variation between life stages.

Fig. 3.

Principal components analyses (PCA) of sample-to-sample distances based on VST-transformed gene expression. (A) PCA of all life stages [n=3 replicates per life stage per treatment, except where n=2 replicate samples per treatment (ambient 1 hpf, extreme low 4 hpf, extreme low 9 hpf, and extreme low 14 hpf)] based on sample-to-sample distance computed from all genes passing a low counts filter, wherein a gene must have a count of 10 or greater in at least 2 out of the 41 samples (pOverA ∼0.05, 10). (B) PCA of the sample-to-sample distance of planula samples only (n=3 replicates per treatment; pOverA ∼0.875, 10).

Fig. 3.

Principal components analyses (PCA) of sample-to-sample distances based on VST-transformed gene expression. (A) PCA of all life stages [n=3 replicates per life stage per treatment, except where n=2 replicate samples per treatment (ambient 1 hpf, extreme low 4 hpf, extreme low 9 hpf, and extreme low 14 hpf)] based on sample-to-sample distance computed from all genes passing a low counts filter, wherein a gene must have a count of 10 or greater in at least 2 out of the 41 samples (pOverA ∼0.05, 10). (B) PCA of the sample-to-sample distance of planula samples only (n=3 replicates per treatment; pOverA ∼0.875, 10).

PCAs computed from variation in expression of genes with greater than 10 counts in seven of eight fertilized egg (10,751 genes; Fig. S1A), cleavage (12,837 genes; Fig. S1B), prawn chip (16,392 genes; Fig. S1C), early gastrula (17,669 genes; Fig. S1D) and planula (23,225 genes) samples showed no differentiation between treatments until the planula stage (Fig. 3B). At the planula stage, treatment groups separated strongly along PC1 and PC2, which explained 47% and 25% of the variation between samples, respectively. Notably, ambient pH (7.8) and extreme low pH (7.3) samples clustered more similarly on PC1 than ambient and low pH (7.6) samples. However, with the exception of one ambient sample, ambient and low pH samples clustered closer than ambient and extreme low pH samples on PC2 (Fig. 3B).

Pairwise treatment comparisons within each life stage using DESeq2 (Table S2) showed zero DEGs between treatments until the planula stage. Differentiation of samples from planula raised in different treatments were supported by 5406 DEGs between planula samples raised under the two low pH conditions relative to the ambient pH condition. Mirroring the clustering of treatments along PC1 in Fig. 3B, a Venn diagram (Fig. 4A) of shared and unique DEGs in each treatment comparison showed 2202 DEGs between planula samples at the low and ambient pH treatments, 1078 of which were also differentially expressed between the extreme low and low treatments. Only 66 DEGs exist between planula samples at the extreme low and ambient pH treatments. Planula DEGs were grouped into two k-means clusters based on the consensus of 11 of 30 clustering indices (Fig. 4B), with genes in cluster 1 upregulated in the low pH treatment compared with the ambient and extreme low pH treatments, and genes in cluster 2 exhibiting the opposite pattern (Fig. 4B).

Fig. 4.

Differential gene expression. (A) Venn diagram and (B) heatmap of all genes that are differentially expressed between treatments at the planula stage. Clusters in the heatmap are based on k-means clustering using the R program NBclust to calculate the optimal number of clusters using 30 indices. Heatmap colors for each gene are based on difference in expression compared with average VST across all samples, where red cells are more highly expressed compared with the row mean and blue cells more lowly expressed compared with the row mean.

Fig. 4.

Differential gene expression. (A) Venn diagram and (B) heatmap of all genes that are differentially expressed between treatments at the planula stage. Clusters in the heatmap are based on k-means clustering using the R program NBclust to calculate the optimal number of clusters using 30 indices. Heatmap colors for each gene are based on difference in expression compared with average VST across all samples, where red cells are more highly expressed compared with the row mean and blue cells more lowly expressed compared with the row mean.

Host functional enrichment

After filtering the DEG set to retain genes with log2FoldChange>|1|, 344 DEGs remained in cluster 1 (upregulated under low pH), and 256 remained in cluster 2 (downregulated under low pH) for GO enrichment analysis. GO enrichment analysis revealed 39 biological processes and 31 molecular functions over-represented in cluster 1. Biological processes overrepresented in cluster 1 (Fig. 5A; Table S3) are primarily related to seven slim terms: biosynthetic process, cellular component assembly, cellular nitrogen compound metabolic process, immune system process, lipid metabolic process, response to stress and transport. Likewise, molecular functions overrepresented in cluster 1 (Fig. 5B; Table S3) can be summarized by eight slim terms: ATPase activity, cytoskeletal protein binding, DNA binding, ion binding, kinase activity, oxidoreductase activity, structural molecule activity and transferase activity. Ion binding, in particular, includes several offspring terms under the GO category calcium ion binding. In cluster 2, GO enrichment revealed five overrepresented biological processes (Fig. 5A; Table S3) represented by the slim terms biosynthetic process, immune system process and transport, as well as 10 overrepresented molecular functions (Fig. 5B; Table S3) represented by the slim terms ATPase activity, nucleotidyltransferase activity, peptidase activity and transmembrane transporter activity. The GO terms enriched in clusters 1 and 2 are supported by significant overrepresentation of 93 KEGG terms related to them (Table S3). This set of significantly overrepresented KEGG terms includes peroxidase (K19511), Amt family ammonium transporter (K03320), cytosolic phospholipase A2 (K16342) and calreticulin (K08057) in cluster 1, and E3 ubiquitin-protein ligase (K22754) and [histone H4]-lysine20 N-methyltransferase SETD8 (K11428) in cluster 2.

Fig. 5.

Gene ontology enrichment analysis. (A) Biological process and (B) molecular function terms that are enriched in planulae developing in the low treatment compared with the extreme low and ambient treatments. Up: terms associated with DEGs (Padj<0.05, log2FoldChange>|1|, numInCat>5) from cluster 1 (C1); down: terms associated with DEGs (Padj<0.05, log2FoldChange>|1|, numInCat>5) from cluster 2 (C2).

Fig. 5.

Gene ontology enrichment analysis. (A) Biological process and (B) molecular function terms that are enriched in planulae developing in the low treatment compared with the extreme low and ambient treatments. Up: terms associated with DEGs (Padj<0.05, log2FoldChange>|1|, numInCat>5) from cluster 1 (C1); down: terms associated with DEGs (Padj<0.05, log2FoldChange>|1|, numInCat>5) from cluster 2 (C2).

Symbiont transcriptome assembly and gene expression

The M. capitata planula and adult holobiont de novo assembly was constructed from 698,861,565 bases, yielding 991,330 transcripts and 656,603 genes. The N50 (measure of the average length of scaffolds in a genome assembly) of the holobiont assembly was 1102, with an average contig length of 704.97 bases and a GC content (proportion of guanine/cytosine residues in a genome assembly) of 43.7%. A BUSCO search of the holobiont assembly against the metazoan and Alveolata databases for the host and symbiont showed high completeness (90.1% Alveolata and 97.5% metazoan). In total, 83.37% (826,513) of transcripts from the holobiont assembly were assigned to the host and 16.63% (164,817) of transcripts were assigned to the symbiont. Finally, using the separated assemblies as reference transcriptomes for Trinity, 4.77±0.43% (mean±s.e.m.) of planula sequences aligned to the symbiont, yielding 133,811 symbiont transcripts for differential expression analysis.

After filtering to remove reads with fewer than 10 counts in seven out of eight samples, 1365 genes remained for PCA and differential expression analyses. A PCA (Fig. S2A) of the symbiont samples showed large variation in gene expression between the three low pH (7.6) samples along the PC1 (28%) and PC2 (14%) axes, but the samples from the ambient pH (7.8) and extreme low pH (7.3) each clustered by treatment. Differential gene expression analysis with DESeq2 (Table S4) revealed 12 DEGs between the ambient and low treatments, 17 DEGs between the ambient and extreme low treatments, and 0 DEGs between the low and extreme low treatments. Ten symbiont genes upregulated in the low treatment compared with the ambient aligned to the histone H3 protein as well as actin and tubulin proteins, whereas the two downregulated genes aligned to the pre-mRNA-splicing factor ATP-dependent RNA helicase and a cell wall-associated hydrolase (Table S4). However, the 17 symbiont genes upregulated in extreme low treatment (Table S4) compared with the ambient were primarily associated genes regulating rDNA expression, a mucin-like protein, cytochrome P450, COPB protein, as well as several hypothetical proteins.

Limited research has been conducted so far on the effects of OA on coral embryonic development, with a primary focus on physiological processes (Chua et al., 2013a,b), although, little is known about physiological changes occurring at the transcriptomic level. Here, we showed that although acidic seawater altered the morphological development of M. capitata embryos, it did not drive any substantial change in gene expression until after the MZT. Early gastrulae and planulae developing at pH 7.6 (low) and pH 7.3 (extreme low) were significantly smaller than those developing at pH 7.8 (ambient), suggesting an energetic cost to developing under low pH. This is also supported by functional enrichment of the genes associated with cell division, regulation of transcription, lipid metabolism and oxidative stress, which are significantly upregulated in swimming larvae (9 dpf) exposed to pH 7.6. We additionally showed that transcriptomically, larvae developing at pH 7.8 and pH 7.3 were more similar than those developing at pH 7.6, and hypothesize that this may be due to high CO2 stimulation of symbiont photosynthesis at pH 7.3, which may alleviate the effects of OA on host energetics.

The MZT under ocean acidification

Although the morphological effects of low pH were apparent by the early gastrula stage (22 hpf; Fig. 2D,E), pH did not drive any substantial changes in gene expression until the larval stage, i.e. after the zygotic transcriptome has assumed autonomy (Fig. S1; Chille et al., 2021). These results suggest that the MZT in M. capitata is robust to low pH and that maternal gene products (i.e. mRNA, proteins, transcription factors) may play a role in buffering early embryos from environmental stress (Hamdoun and Epel, 2007). Although our study does not include survivorship data, this finding is supported by other studies on OA stress during the embryonic development of fish (Dahlke et al., 2017, 2020) and sea urchins (Ericson et al., 2012; Foo et al., 2012) that show that although early development is robust to OA, embryos are more vulnerable after the MZT. In M. capitata, early development is likely to be influenced by maternal transcripts until 22–36 hpf, when the zygotic transcriptome exhibits the capacity for responding to stimuli (Chille et al., 2021). Therefore, the minimal transcriptomic response to pH prior to the swimming larval stage is unsurprising, given the embryo's limited capacity for gene expression regulation prior to late gastrulation. Seeing as a transcriptomic response to decreased pH was only apparent at the planula stage (Figs 3B and 4), this suggests that either the MZT is robust to pH stress, or the outcomes of developing under low pH may be duration dependent. However, even though the MZT does appear to be robust to these conditions, maintaining homeostasis under low pH comes with an energetic cost, as indicated by the smaller sizes of embryos developing at low pH. As such, although the embryonic development of M. capitata may be robust to OA, the energetic cost of developing under low pH may have latent effects on larval competence and recruitment (Albright, 2011).

Outcomes of developing in low pH

In our study, low pH (pH 7.6) elicited a robust transcriptomic response in planula that implicates several processes affecting energy metabolism and larval size (Fig. 2E), including depressed metabolism, cellular stress response and biomineralization. We detail these responses below.

Acid–base parameters such as pH and pCO2 are generally considered to lead to depressed metabolism through cellular acidosis, leading to alterations in cell size (Kaniewska et al., 2012; Reipschläger and Pörtner, 1996). Here, the enrichment of terms related to ATP synthesis in the genes downregulated at low pH (pH 7.6; Table S3), and lipid metabolism in the genes upregulated at low pH (pH 7.6; Table S3) suggest alteration of larval metabolism at reduced pH. Metabolism is highly interconnected with both growth and cell cycle progression (Björklund, 2019). Additionally, mitochondria play a crucial role in setting growth rates through energy generation and lipid biosynthesis (Miettinen and Björklund, 2017). In this context, alteration in metabolism at low pH might cause changes in cell size. Here, the enrichment of terms related to cell growth and cytoskeleton in the genes upregulated at low pH (pH 7.6; Table S3) support changes in cell cycle and cell size with acidification. Cell growth-related genes such as JNK cascade and GTPase activity have been shown to change with acidification (Liew et al., 2018), and acidic cellular pH was demonstrated to alter processes such as cytoskeletal integrity (Tominaga and Barber, 1998). Such alterations in cytoskeletal integrity could explain the observed smaller size of planula that developed under reduced pH conditions.

Chronic exposure to low pH (pH 7.6) during early development also appears to elicit the onset of cellular stress and antioxidant responses. For example, GO terms enriched in planula developing at pH 7.6, including oxidative stress, immune response, superoxide metabolic process, oxidation–reduction process and microtubule-based process (Table S3), are linked to protein domains that protect the cell from oxidative stress through oxidoreductase activity (e.g. superoxide dismutase-like, peroxidasin-like and FAD-linked oxidase; DeSalvo et al., 2008, 2010; Fernandes et al., 2007; Kaniewska et al., 2012; Lesser and Farrell, 2004; Martínez, 2007; Sampayo et al., 2003; Zhao et al., 2016) and the removal of damaged cells (e.g. Toll/interleukin-1 receptor, TNF receptor-associated factor 4; Fuchs and Steller, 2015; Linkermann and Green, 2014; Poquita-Du et al., 2019; Todgham and Hofmann, 2009; van de Water et al., 2015). Upregulation of cellular stress and antioxidant responses is consistent with stress response mechanisms documented in adult and larval corals upon exposure to multiple stressors, including low pH (Kaniewska et al., 2012), UV radiation (Lesser and Farrell, 2004), and thermal (DeSalvo et al., 2008, 2010; Polato et al., 2013; Rosic et al., 2014) and nutrient stress (Rosic et al., 2014). Here, cellular stress and antioxidant responses may have been induced from macromolecular damage owing to ROS released from the rapid peroxidation of planula lipid stores for energy use (Polato et al., 2013). Notably, however, the onset of cellular stress response may also come at an energetic cost to the planula (Sokolova et al., 2012) and lead to smaller larvae by depleting maternal lipid stores.

Several studies have found that actively calcifying larvae may be more susceptible to OA compared with non-calcifying embryos and larvae (Byrne, 2012; Przeslawski et al., 2015) because biomineralization is more chemically and energetically difficult under low pH (Ries, 2011). In the brooding coral Stylophora pistillata, calcium carbonate deposition is initiated during the pre-settled larval phase (Akiva et al., 2018). Although mineralization was not specifically examined in our study, as 9-day-old larvae, it is likely that the planulae here were at the calcifying stage. In fact, GO enrichment analysis shows that genes associated with calcium ion binding and calcification, such as calreticulin, are upregulated in larvae exposed to pH 7.6 (Table S3). These genes provide structural and cell adhesion functions that aid in the biomineralization process (Mass et al., 2016) and may be upregulated to compensate for the chemical shift in carbon equilibrium in the seawater (Comeau et al., 2018; DeCarlo et al., 2018; Zhao et al., 2016). However, upregulation of these genes in response to OA may increase energetic demand, as observed in some adult corals (Davies et al., 2016; Vidal-Dupiol et al., 2013; Scucchia et al., 2021b). This increased demand could potentiate the rapid depletion of larval lipid stores if the energy derived from their endosymbionts is insufficient to keep up with demand (Alamaru et al., 2009; Ben-David-Zaslow and Benayahu, 2000; Harii et al., 2010), ultimately decreasing larval size.

Symbiont stimulation under high CO2

Predicted future increases of CO2 in seawater have been shown to stimulate the photosynthetic activity of the algal symbiont (Krief et al., 2010; Reynaud et al., 2003). Such stimulation can enhance host fitness under acidification (Guillermic et al., 2021). In the present study, the enrichment of terms related to coral oxidative stress in larvae developing under pH 7.6 (low pH; Table S3) compared with larvae developing under pH 7.3 (extreme low pH) suggests low pH (i.e. increased availability of CO2 in seawater) stimulates photosynthesis, thus capacity for electron transport, and therefore potentially minimizes the generation of reactive oxygen species. Such a stimulation of photosynthesis during exposure to low pH was also observed in the endosymbionts of S. pistillata primary polyps and adults (Scucchia et al., 2021a,b). It has been hypothesized that at higher concentrations of CO2, the stimulation of symbiont photosynthesis could prevent acidosis of the host cells, thus helping with preserving cellular acid–base homeostasis (Gibbin et al., 2014), and ultimately alleviating the effects of OA on the host fitness. The increase in translocated photosynthates provided by the endosymbiont under this scenario could additionally supplement the larvae's diet (Harii et al., 2010), thereby decreasing consumption of endogenous lipids and the associated influx of ROS. The differential expression of housekeeping functions, including tubulin, rRNA promoter binding proteins and rDNA transcriptional regulator 15, in our symbiont data (Table S3) suggest to us that there is not enough endosymbiont expression coverage in the dataset to have the power to test this hypothesis. However, it is possible that owing to the potential for post-transcriptional regulation of gene expression in the Symbiodiniaceae through the use of trans-splicing of spliced leader sequences (Zhang et al., 2007; Lin et al., 2010; Lin, 2011), small RNAs and microRNAs (Baumgarten et al., 2013; Lin et al., 2015), and/or RNA editing (Liew et al., 2017), differential expression may not necessarily be observed.

Conclusions

In this study, we examined developmental sensitivity to, and phenotypic consequences of, OA during and after the MZT in the rice coral, M. capitata. Although on a transcriptomic level, no differential gene expression was observed in the host until after the MZT, phenotypically, early gastrulae (14 hpf) and swimming larvae developing in the pH 7.6 and pH 7.3 treatments were significantly smaller than those in the control pH 7.8 treatment. This supports the idea that the development of M. capitata embryos is resilient to low pH conditions, even those substantially lower than what currently occurs in Kāneʻohe Bay. Although the MZT appears to be robust to OA, the energetic strain on embryos and larvae from maintaining homeostasis under low pH could create negative carryover effects for the energetically costly process of metamorphosis (Edmunds et al., 2013) and may make recruits more susceptible to competition and predation (Anlauf et al., 2011; Pechenik et al., 2002; Thiyagarajan et al., 2007). Lowered fitness in these early life stages could create an ontogenetic bottleneck that affects the repopulation and replenishment of coral reefs and contributes to their continued decline (Hughes and Tanner, 2000). Of note, however, the reductions in body size detected so far in coral larvae represent the effects of future OA conditions on present-day populations without consideration of potential moderating effects of rapid adaptation or acclimatization, or prior parental exposure to OA, which may be pivotal to include in forecasting the impacts of climate change on marine organisms (Kelly et al., 2013; Putnam, 2021).

We would like to extend our gratitude to the faculty and staff of the Hawai'i Institute of Marine Biology and the University of Rhode Island Computing for their technical support.

Author contributions

Conceptualization: T.M., H.M.P.; Methodology: E.L.S., E.E.C.; Formal analysis: E.E.C., H.M.P.; Investigation: E.E.C., E.L.S., V.S., M.N., M.O.S., H.M.P.; Resources: H.M.P.; Data curation: E.E.C.; Writing - original draft: E.E.C., E.L.S., F.S., H.M.P.; Writing - review & editing: E.E.C., E.L.S., F.S., V.S., M.N., M.O.S., T.M., H.M.P.; Visualization: E.E.C., H.M.P.; Supervision: H.M.P.; Funding acquisition: T.M., H.M.P.

Funding

This work was funded by United States-Israel Binational Science Foundation grant 2016321 to H.M.P. and T.M. This work was also partially supported by the USDA National Institute of Food and Agriculture, Hatch Formula project accession number 1017848.

Data availability

All raw data can be accessed on the NCBI SRA repository under the BioProject accession IDs PRJNA731596 (adult) and PRJNA616341 (all other time points). Additionally, all scripts used for bioinformatic and statistical analyses are available on the GitHub repository, Mcapitata_OA_Developmental_Gene_Expression_Timeseries (https://github.com/echille/Mcapitata_OA_Developmental_Gene_Expression_Timeseries). All data and scripts are available from Zenodo: doi:10.5281/zenodo.7072713.

Adjeroud
,
M.
,
Kayal
,
M.
and
Penin
,
L.
(
2017
).
Importance of recruitment processes in the dynamics and resilience of coral reef assemblages
. In
Marine Animal Forests: The Ecology of Benthic Biodiversity Hotspots
(ed.
S.
Rossi
,
L.
Bramanti
,
A.
Gori
and
C.
Orejas
), pp.
549
-
569
.
Cham
:
Springer International Publishing
.
Akiva
,
A.
,
Neder
,
M.
,
Kahil
,
K.
,
Gavriel
,
R.
,
Pinkas
,
I.
,
Goobes
,
G.
and
Mass
,
T.
(
2018
).
Minerals in the pre-settled coral Stylophora pistillata crystallize via protein and ion changes
.
Nat. Commun.
9
,
1880
.
Alamaru
,
A.
,
Yam
,
R.
,
Shemesh
,
A.
and
Loya
,
Y.
(
2009
).
Trophic biology of Stylophora pistillata larvae: evidence from stable isotope analysis
.
Mar. Ecol. Prog. Ser.
383
,
85
-
94
.
Albright
,
R.
(
2011
).
Reviewing the effects of ocean acidification on sexual reproduction and early life history stages of reef-building corals
.
J. Mar. Biol.
2011
,
1
-
14
.
Albright
,
R.
and
Langdon
,
C.
(
2011
).
Ocean acidification impacts multiple early life history processes of the Caribbean coral Porites astreoides: ocean acidification impacts coral recruitment
.
Glob. Change Biol.
17
,
2478
-
2487
.
Albright
,
R.
and
Mason
,
B.
(
2013
).
Projected near-future levels of temperature and pCO2 reduce coral fertilization success
.
PLoS One
8
,
e56468
.
Albright
,
R.
,
Mason
,
B.
and
Langdon
,
C.
(
2008
).
Effect of aragonite saturation state on settlement and post-settlement growth of Porites astreoides larvae
.
Coral Reefs
27
,
485
-
490
.
Albright
,
R.
,
Mason
,
B.
,
Miller
,
M.
and
Langdon
,
C.
(
2010
).
Ocean acidification compromises recruitment success of the threatened Caribbean coral Acropora palmata
.
Proc. Natl. Acad. Sci. USA
107
,
20400
-
20404
.
Altschul
,
S. F.
,
Madden
,
T. L.
,
Schäffer
,
A. A.
,
Zhang
,
J.
,
Zhang
,
Z.
,
Miller
,
W.
and
Lipman
,
D. J.
(
1997
).
Gapped BLAST and PSI-BLAST: a new generation of protein database search programs
.
Nucleic Acids Res.
25
,
3389
-
3402
.
Anlauf
,
H.
,
D'Croz
,
L.
and
O'Dea
,
A.
(
2011
).
A corrosive concoction: the combined effects of ocean warming and acidification on the early growth of a stony coral are multiplicative
.
J. Exp. Mar. Biol. Ecol.
397
,
13
-
20
.
Aramaki
,
T.
,
Blanc-Mathieu
,
R.
,
Endo
,
H.
,
Ohkubo
,
K.
,
Kanehisa
,
M.
,
Goto
,
S.
and
Ogata
,
H.
(
2020
).
KofamKOALA: KEGG Ortholog assignment based on profile HMM and adaptive score threshold
.
Bioinformatics
36
,
2251
-
2252
.
Ashburner
,
M.
,
Ball
,
C. A.
,
Blake
,
J. A.
,
Botstein
,
D.
,
Butler
,
H.
,
Cherry
,
J. M.
,
Davis
,
A. P.
,
Dolinski
,
K.
,
Dwight
,
S. S.
,
Eppig
,
J. T.
et al. 
(
2000
).
Gene ontology: tool for the unification of biology. The Gene Ontology Consortium
.
Nat. Genet.
25
,
25
-
29
.
Barros
,
P.
,
Sobral
,
P.
,
Range
,
P.
,
Chícharo
,
L.
,
Matias
,
D.
(
2013
).
Effects of sea-water acidification on fertilization and larval development of the oyster Crassostrea gigas
.
J. Exp. Mar. Biol. Ecol
440
,
200
-
206
.
Barton
,
A.
,
Hales
,
B.
,
Waldbusser
,
G. G.
,
Langdon
,
C.
and
Feely
,
R. A.
(
2012
).
The Pacific oyster, Crassostrea gigas, shows negative correlation to naturally elevated carbon dioxide levels: implications for near-term ocean acidification effects
.
Limnol. Oceanogr.
57
,
698
-
710
.
Baumgarten
,
S.
,
Bayer
,
T.
,
Aranda
,
M.
,
Liew
,
Y. J.
,
Carr
,
A.
,
Micklem
,
G.
and
Voolstra
,
C. R.
(
2013
).
Integrating microRNA and mRNA expression profiling in Symbiodinium microadriaticum, a dinoflagellate symbiont of reef-building corals
.
BMC Genom.
14
,
704
.
Ben-David-Zaslow
,
R.
and
Benayahu
,
Y.
(
2000
).
Biochemical composition, metabolism, and amino acid transport in planula-larvae of the soft coral Heteroxenia fuscescens
.
J. Exp. Zool.
287
,
401
-
412
.
Bindoff
,
N. L.
,
Cheung
,
W. W. L.
,
Kairo
,
J. G.
,
Arístegui
,
J.
,
Guinder
,
V. A.
,
Hallberg
,
R.
,
Hilmi
,
N. J. M.
,
Jiao
,
N.
,
Karim
,
M. S.
,
Levin
,
L.
et al. 
(
2019
).
Changing ocean, marine ecosystems, and dependent communities
. In IPCC Special Report on the Ocean and Cryosphere in a Changing Climate (ed.
H.-O.
Pörtner
,
D. C.
Roberts
,
V.
Masson-Delmotte
,
P.
Zhai
,
M.
Tignor
,
E.
Poloczanska
,
K.
Mintenbeck
,
A.
Alegría
,
M.
Nicolai
and
A.
Okem
et al. 
), pp.
477
-
587
.
Cambridge
:
Cambridge University Press
.
Björklund
,
M.
(
2019
).
Cell size homeostasis: metabolic control of growth and cell division
.
Biochim. Biophys. Acta Mol. Cell Res.
1866
,
409
-
417
.
Byrne
,
M.
(
2012
).
Global change ecotoxicology: identification of early life history bottlenecks in marine invertebrates, variable species responses and variable experimental approaches
.
Mar. Environ. Res.
76
,
3
-
15
.
Byrne
,
M.
and
Przeslawski
,
R.
(
2013
).
Multistressor impacts of warming and acidification of the ocean on marine invertebrates’ life histories
.
Integr. Comp. Biol.
53
,
582
-
596
.
Byrne
,
M.
,
Ho
,
M.
,
Selvakumaraswamy
,
P.
,
Nguyen
,
H. D.
,
Dworjanyn
,
S. A.
and
Davis
,
A. R.
(
2009
).
Temperature, but not pH, compromises sea urchin fertilization and early development under near-future climate change scenarios
.
Proc. Biol. Sci.
276
,
1883
-
1888
.
Caldeira
,
K.
and
Wickett
,
M. E.
(
2003
).
Oceanography: anthropogenic carbon and ocean pH
.
Nature
425
,
365
.
Chen
,
H.
and
Boutros
,
P. C.
(
2011
).
VennDiagram: a package for the generation of highly-customizable Venn and Euler diagrams in R
.
BMC Bioinform.
12
,
35
.
Chille
,
E. E.
,
Strand
,
E.
,
Neder
,
M.
,
Schmidt
,
V.
,
Sherman
,
M.
,
Putnam
,
H. M.
and
Mass
,
T.
(
2021
).
Developmental series of gene expression clarifies maternal mRNA provisioning and maternal-to-zygotic transition in the reef-building coral Montipora capitata
.
BMC Genom.
22
,
815
.
Chua
,
C. M.
,
Leggat
,
W.
,
Moya
,
A.
and
Baird
,
A. H.
(
2013a
).
Temperature affects the early life history stages of corals more than near future ocean acidification
.
Mar. Ecol. Prog. Ser.
475
,
85
-
92
.
Chua
,
C. M.
,
Leggat
,
W.
,
Moya
,
A.
and
Baird
,
A. H.
(
2013b
).
Near-future reductions in pH will have no consistent ecological effects on the early life-history stages of reef corals
.
Mar. Ecol. Prog. Ser.
486
,
143
-
151
.
Comeau
,
S.
,
Cornwall
,
C. E.
,
DeCarlo
,
T. M.
,
Krieger
,
E.
and
McCulloch
,
M. T.
(
2018
).
Similar controls on calcification under ocean acidification across unrelated coral reef taxa
.
Glob. Chang. Biol
24
,
4857
-
4868
.
Cribari-Neto
,
F.
and
Zeileis
,
A.
(
2010
).
Beta regression in R
.
J. Stat. Softw.
34
,
1
-
24
.
Dahlke
,
F. T.
,
Leo
,
E.
,
Mark
,
F. C.
,
Pörtner
,
H.-O.
,
Bickmeyer
,
U.
,
Frickenhaus
,
S.
and
Storch
,
D.
(
2017
).
Effects of ocean acidification increase embryonic sensitivity to thermal extremes in Atlantic cod, Gadus morhua
.
Glob. Change Biol.
23
,
1499
-
1510
.
Dahlke
,
F. T.
,
Lucassen
,
M.
,
Bickmeyer
,
U.
,
Wohlrab
,
S.
,
Puvanendran
,
V.
,
Mortensen
,
A.
,
Chierici
,
M.
,
Pörtner
,
H.-O.
and
Storch
,
D.
(
2020
).
Fish embryo vulnerability to combined acidification and warming coincides with a low capacity for homeostatic regulation
.
J. Exp. Biol.
223
,
jeb212589
.
Davies
,
S. W.
,
Marchetti
,
A.
,
Ries
,
J. B.
and
Castillo
,
K. D.
(
2016
).
Thermal and pCO2 stress elicit divergent transcriptomic responses in a resilient coral
.
Front. Mar. Sci.
3
,
112
.
DeCarlo
,
T. M.
,
Comeau
,
S.
,
Cornwall
,
C. E.
and
McCulloch
,
M. T.
(
2018
).
Coral resistance to ocean acidification linked to increased calcium at the site of calcification
.
Proc. Biol. Sci.
285
,
20180564
.
DeSalvo
,
M. K.
,
Voolstra
,
C. R.
,
Sunagawa
,
S.
,
Schwarz
,
J. A.
,
Stillman
,
J. H.
,
Coffroth
,
M. A.
,
Szmant
,
A. M.
and
Medina
,
M.
(
2008
).
Differential gene expression during thermal stress and bleaching in the Caribbean coral Montastraea faveolata
.
Mol. Ecol.
17
,
3952
-
3971
.
DeSalvo
,
M. K.
,
Sunagawa
,
S.
,
Voolstra
,
C. R.
and
Medina
,
M.
(
2010
).
Transcriptomic responses to heat stress and bleaching in the elkhorn coral Acropora palmata
.
Mar. Ecol. Prog. Ser.
402
,
97
-
113
.
Dickson
,
A. G.
(
1990
).
Standard potential of the reaction: AgCl(s)+12H2(g)=Ag(s)+HCl(aq), and the standard acidity constant of the ion HSO4 in synthetic sea water from 273.15 to 318.15 K
.
J. Chem. Thermodyn.
22
,
113
-
127
.
Dickson
,
A. G.
,
Sabine
,
C. L.
and
Christian
,
J. R.
(
2007
).
Guide to Best Practices for Ocean CO2 Measurements
.
North Pacific Marine Science Organization
.
Doney
,
S. C.
,
Fabry
,
V. J.
,
Feely
,
R. A.
and
Kleypas
,
J. A.
(
2009
).
Ocean acidification: the other CO2 problem
.
Ann. Rev. Mar. Sci.
1
,
169
-
192
.
Drupp
,
P.
,
De Carlo
,
E. H.
,
Mackenzie
,
F. T.
,
Bienfang
,
P.
and
Sabine
,
C. L.
(
2011
).
Nutrient inputs, phytoplankton response, and CO2 variations in a semi-enclosed subtropical embayment, Kaneohe Bay, Hawaii
.
Aquat. Geochem.
17
,
473
-
498
.
Edmunds
,
P. J.
,
Cumbo
,
V. R.
and
Fan
,
T.-Y.
(
2013
).
Metabolic costs of larval settlement and metamorphosis in the coral Seriatopora caliendrum under ambient and elevated pCO2
.
J. Exp. Mar. Bio. Ecol.
443
,
33
-
38
.
Ericson
,
J. A.
,
Lamare
,
M. D.
,
Morley
,
S. A.
and
Barker
,
M. F.
(
2010
).
The response of two ecologically important Antarctic invertebrates (Sterechinus neumayeri and Parborlasia corrugatus) to reduced seawater pH: effects on fertilisation and embryonic development
.
Mar. Biol.
157
,
2689
-
2702
.
Ericson
,
J. A.
,
Ho
,
M. A.
,
Miskelly
,
A.
,
King
,
C. K.
,
Virtue
,
P.
and
Byrne
,
M.
(
2012
).
Combined effects of two ocean change stressors, warming and acidification, on fertilization and early development of the Antarctic echinoid Sterechinus neumayeri
.
Polar Biol.
35
,
1027
-
1034
.
Fernandes
,
P. N.
,
Mannarino
,
S. C.
,
Silva
,
C. G.
,
Pereira
,
M. D.
,
Panek
,
A. D.
and
Eleutherio
,
E. C. A.
(
2007
).
Oxidative stress response in eukaryotes: effect of glutathione, superoxide dismutase and catalase on adaptation to peroxide and menadione stresses in Saccharomyces cerevisiae
.
Redox Rep.
12
,
236
-
244
.
Foo
,
S. A.
,
Dworjanyn
,
S. A.
,
Poore
,
A. G. B.
and
Byrne
,
M.
(
2012
).
Adaptive capacity of the habitat modifying sea urchin Centrostephanus rodgersii to ocean warming and ocean acidification: performance of early embryos
.
PLoS One
7
,
e42497
.
Fuchs
,
Y.
and
Steller
,
H.
(
2015
).
Live to die another way: modes of programmed cell death and the signals emanating from dying cells
.
Nat. Rev. Mol. Cell Biol.
16
,
329
-
344
.
Gazeau
,
F.
,
Gattuso
,
J.-P.
,
Dawber
,
C.
,
Pronker
,
A. E.
,
Peene
,
F.
,
Peene
,
J.
,
Heip
,
C. H. R.
and
Middelburg
,
J. J.
(
2010
).
Effect of ocean acidification on the early life stages of the blue mussel Mytilus edulis
.
Biogeosciences
7
,
2051
-
2060
.
Gibbin
,
E. M.
,
Putnam
,
H. M.
,
Davy
,
S. K.
and
Gates
,
R. D.
(
2014
).
Intracellular pH and its response to CO2-driven seawater acidification in symbiotic versus non-symbiotic coral cells
.
J. Exp. Biol.
217
,
1963
-
1969
.
Grabherr
,
M. G.
,
Haas
,
B. J.
,
Yassour
,
M.
,
Levin
,
J. Z.
,
Thompson
,
D. A.
,
Amit
,
I.
,
Adiconis
,
X.
,
Fan
,
L.
,
Raychowdhury
,
R.
,
Zeng
,
Q.
et al. 
(
2011
).
Full-length transcriptome assembly from RNA-Seq data without a reference genome
.
Nat. Biotechnol.
29
,
644
-
652
.
Gu
,
Z.
,
Eils
,
R.
and
Schlesner
,
M.
(
2016
).
Complex heatmaps reveal patterns and correlations in multidimensional genomic data
.
Bioinformatics
32
,
2847
-
2849
.
Guillermic
,
M.
,
Cameron
,
L. P.
,
De Corte
,
I.
,
Misra
,
S.
,
Bijma
,
J.
,
de Beer
,
D.
,
Reymond
,
C. E.
,
Westphal
,
H.
,
Ries
,
J. B.
and
Eagle
,
R. A.
(
2021
).
Thermal stress reduces pocilloporid coral resilience to ocean acidification by impairing control over calcifying fluid chemistry
.
Sci. Adv.
7
,
eaba9958
.
Hamdoun
,
A.
and
Epel
,
D.
(
2007
).
Embryo stability and vulnerability in an always changing world
.
Proc. Natl. Acad. Sci. USA
104
,
1745
-
1750
.
Han
,
Y.
,
Shi
,
W.
,
Tang
,
Y.
,
Zhao
,
X.
,
Du
,
X.
,
Sun
,
S.
,
Zhou
,
W.
and
Liu
,
G.
(
2021
).
Ocean acidification increases polyspermy of a broadcast spawning bivalve species by hampering membrane depolarization and cortical granule exocytosis
.
Aquat. Toxicol.
231
,
105740
.
Harii
,
S.
,
Yamamoto
,
M.
and
Hoegh-Guldberg
,
O.
(
2010
).
The relative contribution of dinoflagellate photosynthesis and stored lipids to the survivorship of symbiotic larvae of the reef-building corals
.
Mar. Biol.
157
,
1215
-
1224
.
Havenhand
,
J. N.
,
Buttler
,
F.-R.
,
Thorndyke
,
M. C.
and
Williamson
,
J. E.
(
2008
).
Near-future levels of ocean acidification reduce fertilization success in a sea urchin
.
Curr. Biol.
18
,
R651
-
R652
.
Helm
,
R. R.
,
Siebert
,
S.
,
Tulin
,
S.
,
Smith
,
J.
and
Dunn
,
C. W.
(
2013
).
Characterization of differential transcript abundance through time during Nematostella vectensis development
.
BMC Genom.
14
,
266
.
Heyward
,
A. J.
(
1986
).
Sexual reproduction in five species of the coral Montipora
.
Hawaii Inst. Mar. Biol. Tech. Rep.
37
,
170
-
178
.
Hughes
,
T. P.
and
Tanner
,
J. E.
(
2000
).
Recruitment failure, life histories, and long-term decline of Caribbean corals
.
Ecology
81
,
2250
-
2263
.
Jury
,
C. P.
,
Thomas
,
F. I. M.
,
Atkinson
,
M. J.
and
Toonen
,
R. J.
(
2013
).
Buffer capacity, ecosystem feedbacks, and seawater chemistry under global change
.
Water
5
,
1303
-
1325
.
Kaniewska
,
P.
,
Campbell
,
P. R.
,
Kline
,
D. I.
,
Rodriguez-Lanetty
,
M.
,
Miller
,
D. J.
,
Dove
,
S.
and
Hoegh-Guldberg
,
O.
(
2012
).
Major cellular and physiological impacts of ocean acidification on a reef building coral
.
PLoS One
7
,
e34659
.
Kelly
,
M. W.
,
Padilla-Gamiño
,
J. L.
and
Hofmann
,
G. E.
(
2013
).
Natural variation and the capacity to adapt to ocean acidification in the keystone sea urchin Strongylocentrotus purpuratus
.
Glob. Change Biol.
19
,
2536
-
2546
.
Kim
,
D.
,
Paggi
,
J. M.
,
Park
,
C.
,
Bennett
,
C.
and
Salzberg
,
S. L.
(
2019
).
Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype
.
Nat. Biotechnol.
37
,
907
-
915
.
Krief
,
S.
,
Hendy
,
E. J.
,
Fine
,
M.
,
Yam
,
R.
,
Meibom
,
A.
,
Foster
,
G. L.
and
Shemesh
,
A.
(
2010
).
Physiological and isotopic responses of scleractinian corals to ocean acidification
.
Geochim. Cosmochim. Acta
74
,
4988
-
5001
.
Kroeker
,
K. J.
,
Kordas
,
R. L.
,
Crim
,
R. N.
and
Singh
,
G. G.
(
2010
).
Meta-analysis reveals negative yet variable effects of ocean acidification on marine organisms
.
Ecol. Lett.
13
,
1419
-
1434
.
Kurihara
,
H.
and
Shirayama
,
Y.
(
2004
).
Effects of increased atmospheric CO2 on sea urchin early development
.
Mar. Ecol. Prog. Ser.
274
,
161
-
169
.
Langmead
,
B.
(
2010
).
Aligning short sequencing reads with Bowtie
.
Curr. Protoc. Bioinform.
Chapter 11
,
Unit 11.7
.
Lesser
,
M. P.
and
Farrell
,
J. H.
(
2004
).
Exposure to solar radiation increases damage to both host tissues and algal symbionts of corals during thermal stress
.
Coral Reefs
23
,
367
-
377
.
Li
,
B.
and
Dewey
,
C. N.
(
2011
).
RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome
.
BMC Bioinform.
12
,
323
.
Lin
,
S.
(
2011
).
Genomic understanding of dinoflagellates
.
Res. Microbiol.
162
,
551
-
569
.
Lin
,
S.
,
Zhang
,
H.
,
Zhuang
,
Y.
,
Tran
,
B.
and
Gill
,
J.
(
2010
).
Spliced leader-based metatranscriptomic analyses lead to recognition of hidden genomic features in dinoflagellates
.
Proc. Natl. Acad. Sci USA
107
,
20033
-
20038
.
Lin
,
S.
,
Cheng
,
S.
,
Song
,
B.
,
Zhong
,
X.
,
Lin
,
X.
,
Li
,
W.
,
Li
,
L.
,
Zhang
,
Y.
,
Zhang
,
H.
,
Ji
,
Z.
et al. 
(
2015
).
The Symbiodinium kawagutii genome illuminates dinoflagellate gene expression and coral symbiosis
.
Science
350
,
691
-
694
.
Liew
,
Y. J.
,
Aranda
,
M.
and
Voolstra
,
C. R.
(
2016
).
Reefgenomics.Org – a repository for marine genomics data
.
Database
2016
,
baw152
.
Liew
,
Y. J.
,
Li
,
Y.
,
Baumgarten
,
S.
,
Voolstra
,
C. R.
and
Aranda
,
M.
(
2017
).
Condition-specific RNA editing in the coral symbiont Symbiodinium microadriaticum
.
PLoS Genet.
13
,
e1006619
.
Liew
,
Y. J.
,
Zoccola
,
D.
,
Li
,
Y.
,
Tambutté
,
E.
,
Venn
,
A. A.
,
Michell
,
C. T.
,
Cui
,
G.
,
Deutekom
,
E. S.
,
Kaandorp
,
J. A.
,
Voolstra
,
C. R.
et al. 
(
2018
).
Epigenome-associated phenotypic acclimatization to ocean acidification in a reef-building coral
.
Sci. Adv.
4
,
eaar8028
.
Linkermann
,
A.
and
Green
,
D. R.
(
2014
).
Necroptosis
.
N. Engl. J. Med.
370
,
455
-
465
.
Love
,
M.
,
Huber
,
W.
and
Anders
,
S.
(
2014a
).
Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2
.
Genome Biol.
15
,
550
.
Love
,
M.
,
Anders
,
S.
and
Huber
,
W.
(
2014b
).
Differential analysis of count data – the DESeq2 package
.
Genome Biol.
15
,
10
-
1186
.
Lueker
,
T. J.
,
Dickson
,
A. G.
and
Keeling
,
C. D.
(
2000
).
Ocean pCO2 calculated from dissolved inorganic carbon, alkalinity, and equations for K1 and K2: validation based on laboratory measurements of CO2 in gas and seawater at equilibrium
.
Mar. Chem.
70
,
105
-
119
.
Malika
,
C.
,
Ghazzali
,
N.
,
Boiteau
,
V.
and
Niknafs
,
A.
(
2014
).
NbClust: an R package for determining the relevant number of clusters in a data Set
.
J. Stat. Softw.
61
,
1
-
36
.
Martínez
,
R.
(
2007
).
Effects of ultraviolet radiation on protein content, respiratory electron transport system (ETS) activity and superoxide dismutase (SOD) activity of Antarctic plankton
.
Polar Biol.
30
,
1159
-
1172
.
Mass
,
T.
,
Putnam
,
H. M.
,
Drake
,
J. L.
,
Zelzion
,
E.
,
Gates
,
R. D.
,
Bhattacharya
,
D.
and
Falkowski
,
P. G.
(
2016
).
Temporal and spatial expression patterns of biomineralization proteins during early development in the stony coral Pocillopora damicornis
.
Proc. Biol. Sci.
283
,
20160322
.
Miettinen
,
T. P.
and
Björklund
,
M.
(
2017
).
Mitochondrial function and cell size: an allometric relationship
.
Trends Cell Biol.
27
,
393
-
402
.
Morita
,
M.
,
Suwa
,
R.
,
Iguchi
,
A.
,
Nakamura
,
M.
,
Shimada
,
K.
,
Sakai
,
K.
and
Suzuki
,
A.
(
2010
).
Ocean acidification reduces sperm flagellar motility in broadcast spawning reef invertebrates
.
Zygote
18
,
103
-
107
.
Moya
,
A.
,
Huisman
,
L.
,
Ball
,
E. E.
,
Hayward
,
D. C.
,
Grasso
,
L. C.
,
Chua
,
C. M.
,
Woo
,
H. N.
,
Gattuso
,
J.-P.
,
Forêt
,
S.
and
Miller
,
D. J.
(
2012
).
Whole transcriptome analysis of the coral Acropora millepora reveals complex responses to CO2-driven acidification during the initiation of calcification
.
Mol. Ecol.
21
,
2440
-
2454
.
Moya
,
A.
,
Huisman
,
L.
,
Forêt
,
S.
,
Gattuso
,
J.-P.
,
Hayward
,
D. C.
,
Ball
,
E. E.
and
Miller
,
D. J.
(
2015
).
Rapid acclimation of juvenile corals to CO2-mediated acidification by upregulation of heat shock protein and Bcl-2 genes
.
Mol. Ecol.
24
,
438
-
452
.
Nakamura
,
M.
,
Ohki
,
S.
,
Suzuki
,
A.
and
Sakai
,
K.
(
2011
).
Coral larvae under ocean acidification: survival, metabolism, and metamorphosis
.
PLoS One
6
,
e14521
.
O'Donnell
,
M. J.
,
Todgham
,
A. E.
,
Sewell
,
M. A.
,
Hammond
,
L. M.
,
Ruggiero
,
K.
,
Fangue
,
N. A.
,
Zippay
,
M. L.
and
Hofmann
,
G. E.
(
2010
).
Ocean acidification alters skeletogenesis and gene expression in larval sea urchins
.
Mar. Ecol. Prog. Ser.
398
,
157
-
171
.
Padilla-Gamiño
,
J. L.
and
Gates
,
R. D.
(
2012
).
Spawning dynamics in the Hawaiian reef-building coral Montipora capitata
.
Mar. Ecol. Prog. Ser.
449
,
145
-
160
.
Parker
,
L. M.
,
Ross
,
P. M.
and
O'Connor
,
W. A.
(
2009
).
The effect of ocean acidification and temperature on the fertilization and embryonic development of the Sydney rock oyster Saccostrea glomerata (Gould 1850)
.
Glob. Change Biol.
15
,
2123
-
2136
.
Parker
,
L. M.
,
Ross
,
P. M.
and
O'Connor
,
W. A.
(
2010
).
Comparing the effect of elevated pCO2 and temperature on the fertilization and early development of two species of oysters
.
Mar. Biol.
157
,
2435
-
2452
.
Pechenik
,
J. A.
,
Jarrett
,
J. N.
and
Rooney
,
J.
(
2002
).
Relationships between larval nutritional experience, larval growth rates, juvenile growth rates, and juvenile feeding rates in the prosobranch gastropod Crepidula fornicata
.
J. Exp. Mar. Biol. Ecol.
280
,
63
-
78
.
Perez
,
F. F.
and
Fraga
,
F.
(
1987
).
Association constant of fluoride and hydrogen ions in seawater
.
Mar. Chem.
21
,
161
-
168
.
Pertea
,
M.
,
Pertea
,
G. M.
,
Antonescu
,
C. M.
,
Chang
,
T.-C.
,
Mendell
,
J. T.
and
Salzberg
,
S. L.
(
2015
).
StringTie enables improved reconstruction of a transcriptome from RNA-seq reads
.
Nat. Biotechnol.
33
,
290
-
295
.
Polato
,
N. R.
,
Altman
,
N. S.
and
Baums
,
I. B.
(
2013
).
Variation in the transcriptional response of threatened coral larvae to elevated temperatures
.
Mol. Ecol.
22
,
1366
-
1382
.
Poquita-Du
,
R. C.
,
Huang
,
D.
,
Chou
,
L. M.
,
Mrinalini
, and
Todd
,
P. A.
(
2019
).
Short term exposure to heat and sediment triggers changes in coral gene expression and photo-physiological performance
.
Front. Mar. Sci.
6
,
121
.
Przeslawski
,
R.
,
Byrne
,
M.
and
Mellin
,
C.
(
2015
).
A review and meta-analysis of the effects of multiple abiotic stressors on marine embryos and larvae
.
Glob. Change Biol.
21
,
2122
-
2140
.
Putnam
,
H. M.
(
2021
).
Avenues of reef-building coral acclimatization in response to rapid environmental change
.
J. Exp. Biol.
224
,
jeb239319
.
Putnam
,
H. M.
,
Mayfield
,
A. B.
,
Fan
,
T. Y.
,
Chen
,
C. S.
and
Gates
,
R. D.
(
2013
).
The physiological and molecular responses of larvae from the reef-building coral Pocillopora damicornis exposed to near-future increases in temperature and pCO2
.
Mar. Biol.
160
,
2157
-
2173
.
Reipschläger
,
A.
and
Pörtner
,
H. O.
(
1996
).
Metabolic depression during environmental stress: the role of extracellular versus intracellular pH in Sipunculus nudus
.
J. Exp. Biol.
199
,
1801
-
1807
.
Reyes-Bermudez
,
A.
,
Villar-Briones
,
A.
,
Ramirez-Portilla
,
C.
,
Hidaka
,
M.
and
Mikheyev
,
A. S.
(
2016
).
Developmental progression in the coral Acropora digitifera is controlled by differential expression of distinct regulatory gene networks
.
Genome Biol. Evol.
8
,
851
-
870
.
Reynaud
,
S.
,
Leclercq
,
N.
,
Romaine-Lioud
,
S.
,
Ferrier-Pagés
,
C.
,
Jaubert
,
J.
and
Gattuso
,
J.-P.
(
2003
).
Interacting effects of CO2 partial pressure and temperature on photosynthesis and calcification in a scleractinian coral: effects of pCO2 and temperature on coral
.
Glob. Change Biol.
9
,
1660
-
1668
.
Ries
,
J. B.
(
2011
).
A physicochemical framework for interpreting the biological calcification response to CO2-induced ocean acidification
.
Geochim. Cosmochim. Acta
75
,
4053
-
4064
.
Ritson-Williams
,
R.
and
Arnold
,
S. N.
(
2009
).
New perspectives on ecological mechanisms affecting coral recruitment on reefs
.
Smithson. Contrib. Mar. Sci.
38
,
437
.
Rivest
,
E. B.
and
Hofmann
,
G. E.
(
2014
).
Responses of the metabolism of the larvae of Pocillopora damicornis to ocean acidification and warming
.
PLoS One
9
,
e96172
.
Rosic
,
N.
,
Kaniewska
,
P.
,
Chan
,
C.-K. K.
,
Ling
,
E. Y. S.
,
Edwards
,
D.
,
Dove
,
S.
and
Hoegh-Guldberg
,
O.
(
2014
).
Early transcriptional changes in the reef-building coral Acropora aspera in response to thermal and nutrient stress
.
BMC Genom.
15
,
1052
.
Sabine
,
C. L.
,
Feely
,
R. A.
,
Gruber
,
N.
,
Key
,
R. M.
,
Lee
,
K.
,
Bullister
,
J. L.
,
Wanninkhof
,
R.
,
Wong
,
C. S.
,
Wallace
,
D. W. R.
,
Tilbrook
,
B.
et al. 
(
2004
).
The oceanic sink for anthropogenic CO2
.
Science
305
,
367
-
371
.
Sampayo
,
J. N.
,
Olsen
,
A.
and
Lithgow
,
G. J.
(
2003
).
Oxidative stress in Caenorhabditis elegans: protective effects of superoxide dismutase/catalase mimetics
.
Aging Cell
2
,
319
-
326
.
Schneider
,
C. A.
,
Rasband
,
W. S.
and
Eliceiri
,
K. W.
(
2012
).
NIH Image to ImageJ: 25 years of image analysis
.
Nat. Methods
9
,
671
-
675
.
Scucchia
,
F.
,
Malik
,
A.
,
Zaslansky
,
P.
,
Putnam
,
H. M.
and
Mass
,
T.
(
2021a
).
Combined responses of primary coral polyps and their algal endosymbionts to decreasing seawater pH
.
Proc. Biol. Sci.
288
,
20210328
.
Scucchia
,
F.
,
Malik
,
A.
,
Putnam
,
H. M.
and
Mass
,
T.
(
2021b
).
Genetic and physiological traits conferring tolerance to ocean acidification in mesophotic corals
.
Global Change Biology
27
,
5276
-
5294
.
Shaw
,
E. C.
,
Hamylton
,
S. M.
and
Phinn
,
S. R.
(
2016
).
Incorporating benthic community changes into hydrochemical-based projections of coral reef calcium carbonate production under ocean acidification
.
Coral Reefs
35
,
739
-
750
.
Shumaker
,
A.
,
Putnam
,
H. M.
,
Qiu
,
H.
,
Price
,
D. C.
,
Zelzion
,
E.
,
Harel
,
A.
,
Wagner
,
N. E.
,
Gates
,
R. D.
,
Yoon
,
H. S.
and
Bhattacharya
,
D.
(
2019
).
Genome analysis of the rice coral Montipora capitata
.
Sci. Rep.
9
,
2571
.
Simão
,
F. A.
,
Waterhouse
,
R. M.
,
Ioannidis
,
P.
,
Kriventseva
,
E. V.
and
Zdobnov
,
E. M.
(
2015
).
BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs
.
Bioinformatics
31
,
3210
-
3212
.
Sokolova
,
I. M.
,
Frederich
,
M.
,
Bagwe
,
R.
,
Lannig
,
G.
and
Sukhotin
,
A. A.
(
2012
).
Energy homeostasis as an integrative tool for assessing limits of environmental stress tolerance in aquatic invertebrates
.
Mar. Environ. Res.
79
,
1
-
15
.
Strader
,
M. E.
,
Aglyamova
,
G. V.
and
Matz
,
M. V.
(
2018
).
Molecular characterization of larval development from fertilization to metamorphosis in a reef-building coral
.
BMC Genom.
19
,
17
.
Tadros
,
W.
and
Lipshitz
,
H. D.
(
2009
).
The maternal-to-zygotic transition: a play in two acts
.
Development
136
,
3033
-
3042
.
Thiyagarajan
,
V.
,
Pechenik
,
J. A.
,
Gosselin
,
L. A.
and
Qian
,
P. Y.
(
2007
).
Juvenile growth in barnacles: combined effect of delayed metamorphosis and sub-lethal exposure of cyprids to low-salinity stress
.
Mar. Ecol. Prog. Ser.
344
,
173
-
184
.
Timmins-Schiffman
,
E.
,
O'Donnell
,
M. J.
,
Friedman
,
C. S.
and
Roberts
,
S. B.
(
2013
).
Elevated pCO2 causes developmental delay in early larval Pacific oysters, Crassostrea gigas
.
Mar. Biol.
160
,
1973
-
1982
.
Todgham
,
A. E.
and
Hofmann
,
G. E.
(
2009
).
Transcriptomic response of sea urchin larvae Strongylocentrotus purpuratus to CO2-driven seawater acidification
.
J. Exp. Biol.
212
,
2579
-
2594
.
Tominaga
,
T.
and
Barber
,
D. L.
(
1998
).
Na–H exchange acts downstream of RhoA to regulate integrin-induced cell adhesion and spreading
.
Mol. Biol. Cell
9
,
2287
-
2303
.
Van Colen
,
C.
,
Debusschere
,
E.
,
Braeckman
,
U.
,
Van Gansbeke
,
D.
and
Vincx
,
M.
(
2012
).
The early life history of the clam Macoma balthica in a high CO2 world
.
PLoS One
7
,
e44655
.
van de Water
,
J. A. J. M.
,
Ainsworth
,
T. D.
,
Leggat
,
W.
,
Bourne
,
D. G.
,
Willis
,
B. L.
and
van Oppen
,
M. J. H.
(
2015
).
The coral immune response facilitates protection against microbes during tissue regeneration
.
Mol. Ecol.
24
,
3390
-
3404
.
Van Etten
,
J.
,
Shumaker
,
A.
,
Mass
,
T.
,
Putnam
,
H. M.
and
Bhattacharya
,
D.
(
2020
).
Transcriptome analysis provides a blueprint of coral egg and sperm functions
.
PeerJ
8
,
e9739
.
Vastenhouw
,
N. L.
,
Cao
,
W. X.
and
Lipshitz
,
H. D.
(
2019
).
The maternal-to-zygotic transition revisited
.
Development
146
,
dev161471
.
Vidal-Dupiol
,
J.
,
Zoccola
,
D.
,
Tambutté
,
E.
,
Grunau
,
C.
,
Cosseau
,
C.
,
Smith
,
K. M.
,
Freitag
,
M.
,
Dheilly
,
N. M.
,
Allemand
,
D.
and
Tambutté
,
S.
(
2013
).
Genes related to ion-transport and energy production are upregulated in response to CO2-driven pH decrease in corals: new insights from transcriptome analysis
.
PLoS One
8
,
e58652
.
Voolstra
,
C.
,
Miller
,
D.
,
Ragan
,
M.
,
Hoffmann
,
A.
,
Hoegh-Guldberg
,
O.
,
Bourne
,
D.
,
Ball
,
E.
,
Ying
,
H.
,
Foret
,
S.
,
Takahashi
,
S.
et al. 
(
2015
).
The ReFuGe 2020 Consortium—using “omics” approaches to explore the adaptability and resilience of coral holobionts to environmental change
.
Front. Mar. Sci.
2
,
68
.
Wickham
,
H.
(
2011
).
ggplot2: ggplot2
.
Wiley Interdiscip. Rev. Comput. Stat.
3
,
180
-
185
.
Young
,
M. D.
,
Wakefield
,
M. J.
,
Smyth
,
G. K.
and
Oshlack
,
A.
(
2010
).
Gene ontology analysis for RNA-seq: accounting for selection bias
.
Genome Biol.
11
,
R14
.
Yuan
,
X.
,
Yuan
,
T.
,
Huang
,
H.
,
Jiang
,
L.
,
Zhou
,
W.
and
Liu
,
S.
(
2018
).
Elevated CO2 delays the early development of scleractinian coral Acropora gemmifera
.
Sci. Rep.
8
,
2787
.
Zhang
,
H.
,
Hou
,
Y.
,
Miranda
,
L.
,
Campbell
,
D. A.
,
Sturm
,
N. R.
,
Gaasterland
,
T.
and
Lin
,
S.
(
2007
).
Spliced leader RNA trans-splicing in dinoflagellates
.
Proc. Natl. Acad. Sci. USA
104
,
4618
-
4623
.
Zhao
,
X.
,
Ren
,
X.
,
Zhu
,
R.
,
Luo
,
Z.
and
Ren
,
B.
(
2016
).
Zinc oxide nanoparticles induce oxidative DNA damage and ROS-triggered mitochondria-mediated apoptosis in zebrafish embryos
.
Aquat. Toxicol.
180
,
56
-
70
.

Competing interests

The authors declare no competing or financial interests.

Supplementary information