Energetics, but not development, is impacted in coral embryos exposed to ocean acidification

As anthropogenic carbon dioxide (CO₂) levels continue to rise and contribute to global climate change (Bindoff et al., 2019), approximately one-third of atmospheric CO₂ is also being absorbed by the ocean (Sabine et al., 2004). Increasing concentrations of CO₂ 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 CO₂ 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.

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⁻¹ 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).

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 CO₂ 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⁻¹ 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⁻² s⁻¹; n=24), was measured using a handheld PAR sensor (Underwater Quantum Flux Apogee instruments, model MQ-510; accuracy=±4%).

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⁻¹ 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 CO₂ CRM Batch 149) used as a quality control standard. Carbonate chemistry parameters (pH, CO₂, pCO₂, HCO₃, CO₃, 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 K_f (equilibrium constant of hydrogen fluoride, mol kg⁻¹; K_f=[H⁺][F⁻]/[HF]) from Perez and Fraga (1987), K₁ (first dissociation constant of carbonic acid, mol kg⁻¹; K₁=[H⁺][HCO₃⁻]/[CO₂]) and K₂ (second dissociation constant of carbonic acid, mol kg⁻¹; K₂=[H⁺][CO₃²⁻]/[HCO₃⁻]) from Lueker et al. (2000) and K_s (stability constant of hydrogen sulfate, mol kg⁻¹; K_s=[H⁺]_F[SO₄²⁻]/[HSO₄⁻]) from Dickson (1990). Seawater pH was calculated on the total scale (Dickson et al., 2007; SOP6a), using the gas constant (8.31447215 J K⁻¹ mol⁻¹) and the Faraday constant (96485.339924 coulombs mol⁻¹).

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’).

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πab², 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.

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⁻¹) 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).

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; P_adj<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.

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 log₂FoldChange>|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).

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.

Immediately after release, prior to treatment exposure, eggs were 0.040±0.023 mm³ (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.

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.

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.

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).

After filtering the DEG set to retain genes with log₂FoldChange>|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.

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 CO₂ stimulation of symbiont photosynthesis at pH 7.3, which may alleviate the effects of OA on host energetics.

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).

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 pCO₂ 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.

Predicted future increases of CO₂ 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 CO₂ 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 CO₂, 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.

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.