Exposure to global change and microplastics elicits an immune response in an endangered coral

Global change is increasing seawater temperatures and decreasing oceanic pH, driving declines of coral reefs globally. Coral ecosystems are also impacted by local stressors, including microplastics, which are ubiquitous on reefs. While the independent effects of these global and local stressors are well-documented, their interactions remain less explored. Here, we examine the independent and combined effects of global change (ocean warming and acidification) and microplastics exposures on gene expression (GE) and microbial community composition in the endangered coral Acropora cervicornis. Nine genotypes were fragmented and maintained in one of four experimental treatments: 1) ambient conditions (ambient seawater, no microplastics; AMB); 2) microplastics treatment (ambient seawater, microplastics; MP); 3) global change conditions (warm and acidic conditions, no microplastics; OAW); and 4) multistressor treatment (warm and acidic conditions with microplastics; OAW+MP) for 22 days, after which corals were sampled for genome-wide GE profiling and ITS2 and 16S metabarcoding. Overall A. cervicornis GE responses to all treatments were subtle; however, corals in the multistressor treatment exhibited the strongest GE responses, and genes associated with innate immunity were overrepresented in this treatment. ITS2 analyses confirmed that all coral were associated with Symbiodinium ‘fitti’ and 16S analyses revealed similar microbiomes dominated by the bacterial associate Aquarickettsia, suggesting that these A. cervicornis fragments exhibited remarkably low variability in algal and bacterial community compositions. Future work should focus on functional differences across microbiomes, especially Aquarickettsia and viruses, in these responses. Overall, results suggest that when local stressors are coupled with global change, these interacting stressors present unique challenges to this endangered coral species.


Introduction
Anthropogenic global change represents one of the greatest scientific challenges of our time. As atmospheric concentrations of carbon dioxide (CO 2 ) continue to increase, the world's ecosystems face unprecedented challenges from the resulting warming, sea level rise, and extreme weather events (Rosenzweig et al., 2008). The effects of global change are being documented across ecosystems, from declines in terrestrial plant diversity (Harrison, 2020) to shifts in species distributions (Cristofari et al., 2018). However, marine environments are faced with additional threats not present in terrestrial ecosystems (e.g., sealevel rise and ocean acidification), and these threats are especially pressing for tropical coral reefs (Hoegh-Guldberg and Bruno, 2010).
Coral reef ecosystems are experiencing major declines under global change stressors, especially ocean warming that causes the breakdown of the symbiosis between the coral host and dinoflagellate algae [family Symbiodiniaceae; (Muller-Parker et al., 2015;LaJeunesse et al., 2018)]. These declines are particularly evident across the Caribbean where reefs have experienced significant coral loss since the late 20th century, largely due to major disease outbreaks, overfishing leading to macroalgal blooms, and thermal stress events (Contreras-Silva et al., 2020;Bove et al., 2022b;Randazzo-Eisemann et al., 2022). In addition to warming, ocean acidification represents a significant challenge to reef-building corals that can lead to reduced growth rates, dissolution of existing reef framework, and reduced holobiont (animal host, algal symbionts, bacteria, viruses, etc.) physiology and metabolic rates (Anthony et al., 2008;Cornwall et al., 2021). While responses of corals under ocean acidification and warming can vary both within and across species , it is clear these global stressors represent immediate threats to the future of coral reefs.
Along with changing oceanic conditions, coral reefs are susceptible to pollution via terrestrial input, which includes microplastics (Hall et al., 2015;Soares et al., 2020;Hankins et al., 2021;Oldenburg et al., 2021). Microplastics are small plastics, such as beads or fibers, smaller than 5 mm, that come from a variety of sources, including synthetic clothing and personal care products, that enter the waterways via runoff and wastewater treatment (Ivar do Sul and Costa, 2014). Microplastics are ubiquitous across coastal ecosystems, and pose risk to ecosystems via adsorbed chemical pollutants and novel microbiomes (Browne et al., 2011;Rochman et al., 2019). Microplastic particles have been found within coral gastrovascular cavities, suggesting corals ingest plastics (Hall et al., 2015;Hankins et al., 2018;Rotjan et al., 2019). Active ingestion of microplastic particles can reduce coral growth (Reichert et al., 2019;Hankins et al., 2021;Huang et al., 2021) directly via blockages of coral digestive cavities (Allen et al., 2017). Further, interaction with microplastics may impact the energetic budget of the coral, altering polyp behavior and feeding, impeding gamete fertilization success, and affecting the holobiont immune system and coral disease prevalence and susceptibility (Tang et al., 2018;Berry et al., 2019;Rotjan et al., 2019;Hankins et al., 2021;Huang et al., 2021). Passive surface interactions/surface adhesion with the plastics (Reichert et al., 2018) can cause abrasion and injury, increasing coral susceptibility to infection and disease (Page and Willis, 2008;Lamb et al., 2016). While the impacts of microplastics on corals are an active area of research, little is known about how this emerging stressor interacts with other stressors, especially ocean warming and acidification.
Of the work to date examining global change stressors and microplastics, results are variable and species dependent (Axworthy and Padilla-Gamiño, 2019;Huang et al., 2021). For example, recent work found that warming led to consistent reductions in fitness-related traits in corals, while microplastics resulted in mixed responses (Reichert et al., 2021). Similarly, another study found that microplastics had no effect on Symbiodiniaceae cell density under ambient or elevated temperatures (Plafcan and Stallings, 2022), suggesting that microplastics may represent a minor stressor when compared to warming. Little is known about the combined effects of ocean warming and microplastics on corals and even less work has explored how ocean acidification may interact with microplastic pollution to impact corals. Previous work in other marine invertebrates reports impaired immunity in adult mussels (Huang et al., 2022) and altered larval development in urchins (Bertucci et al., 2022) in response to microplastics and ocean acidification, suggesting that the interaction of these stressors may exacerbate coral stress.
One way to gain a more mechanistic understanding of how an organism is responding to a stressor is to profile their gene expression in response to multiple stressors [reviewed in Rivera et al., (2021)]. Currently, gene expression responses of corals to the combined effects of global change and microplastics remain underexplored. Corals have been shown to shuffle their algal symbionts in response to a variety of stressors (Ros et al., 2021;Rodriguez-Casariego et al., 2022), representing a potential acclimation strategy for corals under stress. In addition to their obligate, intracellular photosymbionts (Symbiodiniaceae), corals maintain diverse but specific microbiomes consisting of bacteria, archaea, fungi, viruses, and protists (Bourne et al., 2016;van Oppen and Blackall, 2019). Research has demonstrated that coral health and survival is mediated by their microbiomes (Glasl et al., 2016;Ziegler et al., 2017;Ricci et al., 2019), and that exposure to stressors is reflected in dysbiosis, or a microbiome-wide disturbance/shift in taxonomic composition (Zaneveld et al., 2016;Apprill, 2017;McDevitt-Irwin et al., 2017). Laboratory experiments have demonstrated the capacity of microplastic particles to transmit bacteria into corals via particle ingestion (Rotjan et al., 2019). Moreover, microplastics are known to harbor microbial biofilms that are taxonomically distinct [the "plastisphere", (Zettler et al., 2013;Amaral-Zettler et al., 2020)] from those found in naturally occurring particles suspended in seawater, and bacteria associated with coral diseases have been detected on microplastic pieces (Goldstein et al., 2014;Feng et al., 2020). However, the implications of exposure to microplastics and microplastics-associated microbes on the composition and activity of the coral microbiome, including vectoring pathogens to corals, is not yet well understood.
Understanding how all of the members of the coral holobiont respond in highly endangered coral species, such as Acropora cervicornis, remains a top research priority given the current rates of coral decline (Miller et al., 2002;Schutte et al., 2010;Contreras-Silva et al., 2020). Caribbean A. cervicornis populations are experiencing dramatic declines and are the focus of major restoration efforts (Schopmeyer et al., 2017;Ware et al., 2020). Disease susceptibility in A. cervicornis appears to be mediated by host genotype, environment, and the microbiome (Klinges et al., 2020;Williams et al., 2022). The recently described "Candidatus Aquarickettsia rohweri" (hereafter, Aquarickettsia) (Klinges et al., 2019) is extremely widespread across Caribbean acroporids and can represent as much as 99% of the detected sequences in 16S-based studies (Klinges et al., 2020). Aquarickettsia appears to be parasitic to acroporids and is likely an indicator and determinant of disease susceptibility. Conversely, taxa from the genus Endozoicomonas appear to be significantly and directly correlated with disease resistance and resilience (Chu and Vollmer, 2016). After bleaching/thermal stress, Aquarickettsia abundance declines dramatically in the host, suggesting that disease is the indirect result of Aquarickettsia depletion of holobiont nutritional deficiencies, providing a niche for opportunistic pathogens during/following thermal stress (Klinges et al., 2020). Understanding how the microbiomes of these important reefbuilding corals are involved in holobiont response to stressors is an important consideration for informed conservation, restoration, and management (Ware et al., 2020).
Here, we reared fragments of Acropora cervicornis under projected global change (ocean acidification and warming), microplastic pollution, and their interaction for 22 days. We assessed responses of the holobiont to these stressors via coral gene expression profiling, Symbiodiniaceae community ITS2 m e t a b a r c o d i n g , a n d m i c r o b i o m e c o m m u n i t y 1 6 S metabarcoding. Overall, we hypothesized that the combined stressor treatment (ocean acidification, warming, and microplastics) would elicit the strongest responses by all holobiont members given that these stressors may interact to cause cellular disruption and limit heterotrophy through microplastic ingestion.

Coral collection and experimental design
In June 2020, Nova Southeastern University staff collected nine putative Acropora cervicornis genotypes off coastal Fort Lauderdale, FL, USA from the Nova Southeastern University in situ coral nursery (Florida Fish and Wildlife Conservation Commission permit #SAL-19-2200A-SRP). These genotypes were maintained in outdoor tanks for one week in ambient conditions, after which they were transported to the University of North Carolina Wilmington. Fragments were acclimated to recirculating laboratory tank conditions (27.6-28.9˚C; 34-36 ppt salinity) for three months. Two weeks before the start of the experiment, four branches (~5 cm) from each genotype were glued to ceramic plugs, allowed to recover for one week, and then acclimated to individual experimental tanks for one week.
The microplastic treatments were achieved by weighing out UV-fluorescent blue polyethylene microspheres in equal amounts of two different size classes (density: 1.13 g/cc, Cospheric, LLC; diameters: 180-212 µm and 355-425 µm) to obtain a concentration of 1.25 x 10 -4 g ml -1 per tank. This microplastic concentration was selected based on previous studies to represent moderate microplastic exposure and two size classes were used to ensure the corals could ingest the plastics, which would make the treatment more ecologically relevant (Reichert et al., 2018). The microplastics were added to clean tripour containers within each experimental tank to allow the microplastics to acquire a biofilm over 48 hours before a single coral fragment was placed into each tripour container. The microplastics were then gently resuspended within the individual tripours three times a day using a turkey baster. All tripours were treated the same and this process was repeated throughout the experiment to maintain microplastic densities within each replicate.
Water temperature and pH were controlled using Neptune Systems' Apex microcontrollers. Each tank contained a pH probe, temperature probe, and aquarium heater to control individual tank conditions. The OAW and OAW+MP tanks also each contained a tube for bubbling in CO 2 to control the pH in these systems. The Neptune System allowed for precise control of individual tank parameters by modulating the aquarium heaters and the solenoids connected to the CO 2 tubes to achieve desired treatment conditions. Fragments in OAW/OAW+MP treatments were acclimated by starting at 28.6°C and then raising the temperature by 0.5°C each day until 30.5°C was achieved. Similarly, the pH in OAW/OAW +MP treatments started at 8.05 and was reduced by 0.05 pH units each day until 7.94 was reached. Corals were maintained at these experimental conditions for 22 days while monitoring salinity, temperature, pH, and dissolved oxygen daily (Supplemental Figure 1 and Table S1). Ammonia, nitrate, phosphate, and alkalinity were also measured daily in two randomly selected experimental tanks, one ambient and one ocean acidification and warming. Lights were on a 12:12 light: dark schedule with a PAR of 195 mmol m 2 S -1 that is similar to those recorded for waters across southeast Florida (Yates et al., 2019). All fragments were fed Reef-Roids coral food following manufacturer's instructions (2 mL aliquots) three times weekly and 30% water changes were conducted for all tanks the day after feeding. Immediately following the 3-week experiment, tissue samples (1 cm) from each fragment were placed in cryovials, flash frozen in liquid nitrogen, and maintained at -80°C before being transported to Boston University for downstream analyses.

RNA isolation and TagSeq library preparation
Tissue samples were crushed with a razor blade and RNA was isolated using RNAqueous kits (ThermoFisher Scientific) following manufacturer instructions with the exception of eluting in 30 µl of elution buffer. DNA was then removed using DNA-free DNA Removal kits (ThermoFisher Scientific). Of the nine genotypes, the top six ( Figure 1 and Supplemental Table S2) with the highest concentration and quality of RNA were normalized to 15 ng/µl. Normalized samples were sent to UT Austin's Genome Sequencing and Analysis Facility, where TagSeq libraries were prepared following Meyer et al. (2011) and sequenced (single end 100 bp) on a NovaSeq 6000.

Identifying Acropora cervicornis clones
Acropora cervicornis is well-known for its ability to reproduce asexually via fragmentation (Tunnicliffe, 1981). To identify potential A. cervicornis clones in our dataset, we called Experimental design to test the effects of ocean warming, acidification, and microplastics on the endangered coral Acropora cervicornis. Four fragments from each of six A. cervicornis genotypes that were used in downstream sequencing were assigned one of four experimental treatments: ambient conditions (AMB; blue), ocean acidification and warming (OAW; orange), microplastics (MP; dark blue), and OAW and MP (OAW+MP; red). Fragments were maintained individually in isolated tanks not depicted in this diagram. The dark blue rectangle (left) showcases two corals that were identified as clones via SNP calling from TagSeq gene expression data and the^symbol indicates the clone pair that was removed. The asterisk (*) indicates one gene expression library that was removed from downstream analysis due to low read counts.
single nucleotide polymorphisms (SNPs) from TagSeq reads. Briefly, 104.64 million raw reads were generated, with individual library counts ranging from 4.17 to 6.58 million reads per sample (mean = 5.23 million reads) (Supplemental Table S2). Fastx_toolkit removed 5'-Illumina leader sequences and poly (A) + tails. Sequences <20bp in length with <90% of bases having quality cutoff scores <20 were also trimmed. In addition, because degenerate bases were incorporated during cDNA synthesis, PCR duplicates were removed from all libraries. After quality filtering 0.63 to 3.00 million reads remained (mean = 2.07 million) and these resulting quality filtered reads were mapped to the Acropora millepora genome (Fuller et al., 2020) using Bowtie2.2.0 (Langmead and Salzberg, 2012) (Supplemental Table S2). Resulting SAM files were converted to BAM files using samtools (Li et al., 2009). ANGSD (Korneliussen et al., 2014) calculated pairwise identity-by-state (IBS) matrices using a minIndDepth filter of five. IBS matrices were used as input to hclust(), which visualized relatedness of individuals and allowed for the identification and removal of clones. We note here that mapping was also conducted using the Acropora cervicornis transcriptome (Libro et al., 2013) and the A. cervicornis genome (Kitchen et al., 2019), however the transcriptome was found to contain symbiont data and the genome resulted in lower mapping efficiencies relative to the A. millepora genome. Given that results remained the same (e.g., clones were identified and similar overall gene expression patterns were observed), A. millepora has superior annotations, and Cooke et al. (2020) found high synteny among acroporids, we opted to move forward mapping to the A. millepora genome. When clones were detected, only one clonemate was maintained in downstream analyses of gene expression patterns (see Supplemental Table S2).
2.4 Read mapping, differential gene expression and gene ontology enrichment analysis Reads were trimmed as described above and these reads were mapped to the Acropora millepora genome (Fuller et al., 2020) using Bowtie2 (Langmead and Salzberg, 2012). A custom Perl script (samcount.pl) was used to count the number of reads and the resulting raw counts file was then imported into R v. 3.4.2 (R Core Team (2017) for further analyses. Only samples representing clonemates from SNP analyses above and samples with < 400k counts were removed leaving 19 samples in downstream analyses. Differentially expressed genes (DEGs) were identified using DESeq2 v. 1.32.0 in R, with the model: design =~genotype + treatment. A contig was considered significantly differentially expressed if it had an FDR adjusted p-value < 0.1.
Data were rlog-normalized and effects of genotype and treatment on global gene expression profiles were tested using a PERMANOVA [adonis() function; vegan (Oksanen, 2007)] and visualized using principal components analysis (PCA) using Euclidean distances. An additional PCA was computed using only the top 1000 differentially expressed genes for experimental treatment based on raw p-values to further explore sample clustering when only the genes whose profiles were the most divergent were included.
Gene expression plasticity was then calculated using the first two principal component (PC) axes as the distance between a genotype's gene expression in its treatment condition relative to its expression in ambient conditions using a custom function (Bove, 2022a). The effect of treatment on calculated distances (i.e., gene expression plasticity) was assessed using an analysis of variance (function aov) with treatment as a fixed effect. Genotype was not included in the model because it was accounted for in the PC distance calculations. Differences between treatment levels were tested using Tukey's Honest Significant Differences (HSD) tests.
To determine whether global gene expression patterns showed enrichment of different gene ontology (GO) classes, the collection of scripts 'GO_MWU' from Wright et al. (2015) was used (https://github.com/z0on/GO_MWU). Here, the GO database (go.obo v.1.2) was used to test for enrichment of GO terms based on the ranked −log signed p-values of each gene. Gene ontology terms that were over-represented or underrepresented were then visualized in a tree format that groups GO terms with other terms of similar function. All GO enrichment results for all treatment comparisons can be found on the github repository (https://github.com/daviessw/Acer_ OAW-Microplastics); however only results of the double stressor treatment (OAW+MP) relative to ambient conditions (AMB) are described in detail here. Within this comparison, a heatmap of genes (raw p-value 0.10) assigned to GO terms associated with immunity was constructed using the R package pheatmap to illustrate gene expression patterns associated with OAW+MP.

ITS2 and 16S metabarcoding
DNA was extracted using an RNAqueous kit (ThermoFisher Scientific) as described above, except samples were not subjected to the DNA removal step. Internal transcribed spacer region 2 ( I T S 2 ) P C R a m p l i fi c a t i o n w a s p e r f o r m e d u s i n g SYM_VAR_5.8S2 and SYM_VAR_REV primers (Hume et al., 2013;Hume et al., 2015) using the following PCR profile: 26 cycles of 95°C for 40 s, 59°C for 2 min, 72°C for 1 min and a final extension of 72°C for 7 min. A negative control was included and failed to amplify so it was not sequenced. PCR products were cleaned using the GeneJET PCR Purification kit (ThermoFisher Scientific) according to the manufacturer's instructions with the exception of eluting in 30 µl of elution buffer. A second PCR was performed to dual-barcode samples before pooling, which was done based on the visualization of band intensity on a 1% agarose gel. After pooling, the sample was cleaned using the GeneJET PCR Purification kit (ThermoFisher Scientific) according to the manufacturer's instructions with the exception of using 40 ml of elution buffer. 20 ml of the pool was run on a 2% SYBR Green gel, the target band was excised and placed in 30 ml of Milli-Q water overnight at 4°C before submission for sequencing as detailed below. The V4 region of the 16S rRNA gene was amplified via PCR using Hyb515f (Parada et al., 2016) and Hyb806R (Apprill et al., 2015) primers and the following PCR profile: 30 cycles of 95°C for 40 s, 63°C for 2 min, 72°C for 1 min and a final extension of 7 min. The same procedure as described above for the ITS2 samples was then followed, but with the addition of two negative controls using Milli-Q water which were included in sequencing submission. Concentrations of the ITS2 and 16S pools (via DeNovix DS-11+ Spectrophotometer) were used to combine the two pools in a 1:3 ratio, respectively. Libraries were sequenced on Illumina MiSeq (paired-end 250 bp) at Tufts Genomics Core Facility.
Demultiplexed reads were pre-processed using bbmap (Bushnell, 2014) to split ITS2 and 16S reads based on primers, while tossing reads that included neither primer. Resulting ITS2 reads were then analyzed by submitting paired fastq.gz files directly to SymPortal, which identifies specific sets of defining intragenomic ITS2 sequence variants (DIVs) to define ITS2 type profiles that are indicative of genetically differentiated Symbiodiniaceae taxa (Hume et al., 2019).
16S primers were removed using cutadapt (Martin, 2011), then DADA2 (Callahan et al., 2016) was used to conduct quality filtering and inference of 3,493 amplicon sequence variants (ASVs) (see Supplemental Table S3 to track reads lost through filtering). Taxonomy was assigned with DADA2 using the Silva v. 138.1 database (Quast et al., 2013) and National Center for Biotechnology Information's nucleotide database using blast+ (Camacho et al., 2009). ASVs matching mitochondria, chloroplasts, or non-bacterial kingdoms were removed (216 total) and 5 ASVs were removed based on negative controls as contaminants [Decontam; (Davis et al., 2018)]. Cleaned counts were rarefied to 9,409 using vegan (Oksanen, 2007) and trimmed using MCMC.OTU (Wright et al., 2015) to remove ASVs present in less than 0.01% of counts, resulting in 260 ASVs across all 36 samples.
Beta diversity of the 16S trimmed counts based on treatment and genotype was assessed using a PCoA on Bray-Curtis dissimilarity [Phyloseq; (McMurdie and Holmes, 2013)] and a PERMANOVA [vegan; (Oksanen, 2007)]. We then calculated alpha diversity (Shannon index, Simpson's index, ASV richness, and evenness) of the rarefied counts using Phyloseq [function estimate_richness(); (McMurdie and Holmes, 2013)]. Finally, to test for differences in background 16S communities, we removed all ASVs from the dominant genus MD3-55 [Candidatus Aquarickettsia rohweri, referred now as "Aquarickettsia"; (McMurdie and Holmes, 2013;Klinges et al., 2019)] and performed the same beta and alpha diversity assessments. All ITS2 and 16S data, figures, and analyses were completed in R [version 3.6.3; (R Core Team, 2020) and can be found on GitHub ( h t t p s : / / g i t h u b . c o m / s e a b o v e 7 / a c e r _ m i c r o p l a s t i c s ; (Bove, 2022b)].

Identification of coral genotypes
A total of 4,018 SNPs were identified and of the six putative genotypes that were sequenced, one pair of clones was identified (G10a and G12) (Supplemental Figure 2). Thus, genotype G12 samples from each of the four experimental treatments were removed to retain only one representative clone (genotype G10a) in each treatment (Supplemental Table S2).

Double stressor treatment elicits strongest gene expression plasticity
A total of 20 samples were sequenced, which resulted in a total of 104.6 million reads, 41.4 million of which remained after trimming, and 13.9 million of which mapped successfully to the A. millepora reference genome. One sample -genotype G2, from the OAW treatment -was removed from the dataset due to low counts (166,828 total; Supplemental Table S2). Mean counts across samples after outlier removal was 661,358.
Pairwise comparisons between each stressor (MP, OAW, OAW +MP) and ambient conditions (AMB) showcased that individuals in the double stressor (OAW + MP) treatment had the highest number of differentially expressed genes (DEGs) (83 up, 60 down), followed by microplastics (48 up, 9 down) and then OAW (18 up, 18 down) with relatively few genes overlapped between treatments (p FDR <0.10, Supplemental Figure 3). A principal component analysis (PCA) on all rlog-transformed counts found that there was no overall effect of treatment on A. cervicornis gene expression; however, there was a significant effect of genotype on overall expression patterns ( Figure 3A; p GENOTYPE <0.001). A PCA of the top 1,000 DEGs showcased stronger clustering by experimental treatment and the effect of genotype remained significant ( Figure 3B; p GENOTYPE <0.001, p TREATMENT <0.001). This second PCA also showcased that the double stressor (OAW+MP) had the strongest effect on gene expression, which corroborated the numbers of observed DEGs. A gene expression plasticity analysis confirmed this pattern, which showcased a significant effect of treatment on gene expression plasticity (p TREATMENT =0.041), with the OAW+MP treatment exhibiting the highest plasticity ( Figure 3C). However, these differences in plasticity between corals in the OAW+MP treatment relative to corals in the OAW and MP treatments were only marginally significant after multiple test correction (OAW+MP~MP, p=0.07; OAW +MP~OAW, p=0.06).  Relative abundance of major ITS2 types across each sample grouped by treatment. Light grey represents Symbiodinium 'fitti' (A3) and dark grey represents Breviolum minutum (B2).

Corals exhibit mixed functional responses to single stressor treatments
Gene ontology enrichment analysis between microplastics (MP) relative to ambient conditions (AMB) revealed underrepresentation of oxidoreductase (GO: 0016491), suggesting down regulation of stress response in corals in MP treatment. In addition, overrepresentation of GO terms associated with ribosomal processes was observed including structural constituent of ribosome (GO:0003735) and large ribosomal subunit (GO:0015934), suggesting upregulation of terms associated with increased growth under the presence of microplastics.
In response to OAW, we observed an enrichment of regulation of defense response to virus (GO:0050688), terms associated with amino acid catabolism (aromatic amino acid family catabolic process (GO:0009074), cellular amino acid catabolic process (GO:0009063) and alpha-amino acid catabolic process (GO:1901606)), and two terms associated with stress response (chaperone-mediated protein folding (GO:0061077) and unfolded protein binding (GO:0051082)). In addition, we observed downregulation of processes associated with sensory detection including detection of stimulus (GO:0051606), sensory perception (GO:0007600), and cognition (GO:0050890). Lastly, several GO terms associated with synapse were down regulated including synapse (GO:0045202), synaptic membrane (GO:0097060), and cell junction (GO:0030054). All resulting GO trees from these single stressor treatments can be found on the accompanying GitHub repository (https://github.com/ daviessw/Acer_OAW-Microplastics).

Enrichment of immune-related functions under multiple stressors
Gene ontology enrichment analysis between our double stressor (OAW+MP) relative to ambient conditions (AMB) revealed many significantly enriched GO terms within 'Biological processes', of which the OAW+MP-enriched categories were dominated by GO terms associated with immunity (green box, Supplemental Figure 4). These enriched GO terms included: regulation of defense response to virus by host (GO:0050691), activation of innate immune response (GO:0002218), regulation of immune effector process (GO:0002697), NIK/NF-kappaB signaling (GO:0038061), regulation of cytokine production (GO:0001819), suggesting that the combined stressor elicited an upregulation of innate immunity in the coral host. When a heatmap was made of all genes belonging to these 30 GO terms that met our alpha cut-off (unadjusted p-value < 0.10), strong differences in gene expression were observed between corals in OAW+MP relative to those in AMB treatments (Figure 4). Of particular interest, many classic stress response genes were upregulated in OAW+MP relative to AMB Differentially expressed genes (DEGs, unadjusted p-value < 0.10) with annotations associated with immunity gene ontology (GO) terms from Supplemental Figure 4 (green box). Heatmap showing annotated immunity genes where each row is a gene and each column is a unique gene expression sample. The color scale is in log2 (fold change relative to the gene's mean) and genes and samples are clustered hierarchically based on Pearson's correlation of their expression across samples. Colored blocks indicate treatment. Hierarchical clustering of libraries (columns) demonstrates strong differences in gene expression of immunity genes by treatment condition.
including: superoxide dismutase, two heat shock proteins, several ubiquitin genes, tumor necrosis factors/receptors, a proto-oncogene, and several genes associated with apoptosis (bcl2-like 1) and response to cytokines (p38 map kinase).
Beta diversity analyses did not identify any statistical differences across samples based experimental treatment or genotype (Figures 5C, D and Supplemental Table S4). Similarly, alpha diversity was indistinguishable across treatments (Supplemental Figure 6 and Table S5). After removing Aquarickettsia from all samples to assess background bacterial communities, there were still no differences in alpha or beta diversity across treatments or genotypes (Supplemental Figure 7 and Table S6).

FIGURE 5
Acropora cervicornis bacterial (16S) relative abundance across fragments and experimental treatments. Diversity at the phylum level (A) depicts dominance of all samples by taxa from Proteobacteria and this (B) phylum was dominated by the order Rickettsiales in most samples. Bacterial diversity color was assigned alphabetically, not based on abundance. Beta diversity was visualized through (C) multivariate ordination plots (PCoA) of between-sample Bray-Curtis dissimilarity of rarefied ASVs of all taxa and (D) distance to treatment centroids. PCoA ellipses depict 95% confidence intervals and p-values indicate significance of treatment and genotype. Treatment is depicte-d by color in all panels: light blue = ambient conditions (AMB), dark blue = ocean acidification and warming (OAW), orange = microplastics (MP), and red = OAW and MP (OAW+MP).

Microplastics alone do not drive strong gene expression responses in Acropora cervicornis
Tropical coral reef ecosystems are facing extraordinary challenges from both local and global stressors that are altering their function (Eddy et al., 2021). On the local scale, pollution in the form of microplastics may threaten the health and physiology of reef-building corals (Soares et al., 2020;Nanthini devi et al., 2022). While there is evidence that corals ingest microplastics, which can lead to a variety of physiological responses (Hall et al., 2015;Rotjan et al., 2019), here we detected only a muted gene expression response in A. cervicornis exposed to microplastics. In fact, overall gene expression profiles were indistinguishable between the ambient and microplastics treatments, and an effect of microplastics was only apparent when looking at the expression of top 1,000 genes. This result corroborates previous work reporting negligible effects of microplastics on coral physiology, bleaching susceptibility, and mortality (Reichert et al., 2021;Plafcan and Stallings, 2022) and supports the hypothesis that corals that rely more heavily on heterotrophy may be more impacted by microplastic pollution.
Recent studies have considered the physiological impact of microplastics on tropical corals at the molecular level. These studies, conducted on the tropical coral Pocillopora damicornis (Tang et al., 2018) and the habitat-forming octocoral Corallium rubrum (Corinaldesi et al., 2021), report evidence of elevated stress response after exposure to microplastics in as little as a few hours. In contrast, we identified enrichment of GO terms associated with growth (structural constituent of ribosome, large ribosomal subunit) and an underrepresentation of terms associated with stress response (oxidoreductase). Our gene expression results suggest that A. cervicornis is not exhibiting signs of molecular stress after a 3-week exposure to microplastics, and that these corals may actually be growing more compared to those reared under ambient conditions. It should be noted, however, that the enrichment of growth GO terms in response to microplastics may be due to other compensatory mechanisms, such as increasing photophysiology as a way to counteract declining coral host energy reserves (Reichert et al., 2019;Lanctôt et al., 2020). While we did not quantify photophysiology or gene expression of the algal symbionts in this study, it is possible that algae associated with corals in the microplastic treatment may have exhibited enriched photophysiology, which may have led to the enrichment of molecular signatures of growth. However, it is possible that a longer-term exposure to microplastics may have eventually led to signs of physiological stress (Reichert et al., 2019). While these contrasting results may be due to experimental design differences (microplastic concentration and size, type of microplastics, duration, etc.), it is more likely due to species-specific responses (Hankins et al., 2021;Mendrik et al., 2021). Previous work on A. cervicornis suggests that exposure to microplastics in the presence of warming does not impact bleaching, likely because this species may not ingest the microplastics under otherwise ideal conditions (Plafcan and Stallings, 2022). Instead, A. cervicornis relies on photosynthetically-derived carbon from its algal symbionts as the primary source of nutrition rather than heterotrophy, thus avoiding the potential pitfalls of microplastics ingestion observed in other coral species that are largely heterotrophic (Reichert et al., 2019;Rotjan et al., 2019;Hankins et al., 2021;Mendrik et al., 2021). Overall, we recommend that future work incorporate additional physiological assessments along with coral host and algal symbiont molecular responses to disentangle the role that nutritional source (i.e., dependence on heterotrophy vs. autotrophy) serves in coral physiological response to microplastic pollution.
While these findings suggest that microplastics exposure alone does not elicit a severe stress response in A. cervicornis, it is unclear whether prolonged exposures may eventually impact the host (Hankins et al., 2021) and their algal symbionts (Lanctot et al., 2020;Ripken et al., 2020). Further, microplastic pollution is occurring alongside global stressors (Soares et al., 2020), including widespread shifts in the microbial landscape of seawater (Cavicchioli et al., 2019) and increased prevalence of pathogens and disease across the Caribbean (van Woesik and Randall, 2017). Because microplastics have been shown to be a vehicle for microbial transmission into a corals (Rotjan et al., 2019), there is emerging concern for microplastics pollution in combination with increased pathogen exposure that has long been a predicted component of global change (Cavicchioli et al., 2019).

Global change treatment elicits subtle stress signatures
Along with local concerns, global stressors, such as ocean acidification and warming, are contributing to the wide-scale degradation of coral reef ecosystems. Like its response to microplastics, A. cervicornis in the ocean acidification and warming (OAW) treatment exhibited similar overall gene expression profiles to those in the ambient treatment and only gene expression profiles of the top 1,000 genes resulted in a clear response. This subtle response to OAW in A. cervicornis is surprising since this species is known to be sensitive to these global stressors (Enochs et al., 2014;Towle et al., 2015b;Kaufman et al., 2021;Muller et al., 2021). These contrasting results may be due to different experimental design considerations across studies (Bove et al., 2020), especially given that the heat stress temperature employed here (~30.5˚C) was lower than most other warming studies. However, there is also significant genetic variation within A. cervicornis (Drury et al., 2017;Million et al., 2022), thus it is also possible that the genotypes used in this study were particularly resistant. Alternatively, these corals are already exposed to warmer seawater conditions where they were collected compared to other Caribbean reefs (Muñiz-Castillo et al., 2019;Bove et al., 2022b), suggesting the potential for acclimatization to elevated temperatures.
While we only observed minor gene expression differences in A. cervicornis reared in the OAW treatment, these shifts were associated with an enrichment of stress-related GO terms. For example, corals in OAW treatment exhibited enrichment of GO terms linked to protein folding and protein binding, which have been reported previously in heat-stressed corals in situ (Ip et al., 2022). This pattern was accompanied by enrichment of terms associated with amino acid catabolism, suggesting A. cervicornis may have mediated physiological stress through protein catabolism (Davies et al., 2016;Aguilar et al., 2019;Rädecker et al., 2021). While we did not assess host or symbiont physiological traits along with these gene expression profiles, it is likely that the growth and/or energy reserves may have been impacted as a result of treatment conditions Scucchia et al., 2021;Bove et al., 2022a).
Most interestingly, we identified enrichment of GO terms associated with defense against viruses. Previous work has suggested that corals become more susceptible to infection when water temperatures increase (Bruno et al., 2007). Indeed, patterns of Caribbean coral disease have become more prevalent in recent decades (Randazzo-Eisemann et al., 2022) and population declines of A. cervicornis have been specifically linked to increases in disease associated with seawater temperature increases (Muller et al., 2018;Goergen et al., 2019). While the temperature stress applied here was not high enough to induce a classic heat stress response, it suggests that even subtle increases in temperature may make this species of coral more susceptible to infection.

Global change stressors interact with microplastics to invoke an immune response
Ecological stressors can interact to pose further threats to marine organisms that may not be accounted for when assessing responses to stressors independently (Darling and Cote, 2008;Ellis et al., 2019). Indeed, our multistressor treatment (ocean acidification, warming, and microplastics; OAW+MP) resulted in the strongest gene expression response, suggesting that microplastic pollution may interact with ocean acidification and warming to elicit a more severe molecular response. This strong response to multiple stressors is commonly observed in coral (Coles and Jokiel, 1978;Reynaud et al., 2003;Courtial et al., 2017;Muller et al., 2021) and has been reported previously in studies assessing gene expression as well (Ogawa et al., 2013). The observed response of A. cervicornis to the multistressor treatment in this study may be due to the corals experiencing physiological stress from the acidification and warming exposure (i.e., bleaching), which then led to an increase in reliance on heterotrophy to meet their energetic needs (Ferrier-Pagès et al., 2010;Towle et al., 2015a). This increased heterotrophy in turn may have resulted in the corals consuming more microplastics that do not contribute nutritional value to the coral host, thus furthering the energetic deficit and potentially making them more susceptible to harmful pathogens (Mendrik et al., 2021). While it is unlikely that global change and microplastic pollution are directly interacting, the physiological mechanisms employed by corals in response to these single stressors likely lead to a more severe physiological response.
Here, we identified GO terms enriched in the multistressor treatment related to response to viruses and innate immunity, indicating these stressors may invoke an immune response by the coral. Interestingly, a similar response was observed in OAW alone, however, to a lesser extent, suggesting that the presence of microplastics magnified this stress response. While no studies to date have assessed gene expression responses in corals to the stressors together tested here, previous work in the marine mussel, Mytilus coruscus, exposed to ocean acidification and microplastics identified signals of a repressed immunity that may lead to higher disease susceptibility (Huang et al., 2022). Further, gene expression patterns consistent with an upregulation of immune stress response were reported in the brine shrimp Artemia franciscana exposed to microplastics and warming (Han et al., 2021). This is particularly interesting because Artemia is also frequently used as a nutrition source in coral experiments, including microplastic ingestion experiments (Hall et al., 2015;Axworthy and Padilla-Gamiño, 2019;Lanctot et al., 2020), potentially altering the nutritional value of the Artemia and impacting coral health. Overall, these results suggest that the combination of global change stressors and microplastics will likely lead to suppressed immunity in corals in tandem with higher disease transmission associated with warming (Randazzo-Eisemann et al., 2022) and microplastic pollution (Rochman et al., 2019;Rotjan et al., 2019).

Acropora cervicornis microbial communities dominated by Aquarickettsia
Similar microbiome compositions were observed across the four treatment conditions. A total of 31 of the 36 fragments were high-Aquarickettsia (i.e., greater than 50% dominance), with relative abundances ranging from 17.5% up to 98.8% of the total sequence reads. This sequence dominance in A. cervicornis microbiomes is not uncommon and has been documented in this species in the Florida Keys and throughout the Caribbean , although previous work has also showcased variation among individual genotypes (Aguirre et al., 2022). Even when Aquarickettsia-affiliated ASVs were removed from analysis, microbiome composition did not significantly differ across treatments. Importantly, the microbiome composition of these corals before collection from the wild is unknown, therefore we cannot determine if the universally high Aquarickettsia levels were a result of captivity muting microbiome signatures (Galand et al., 2018). Future work leveraging experiments on individuals representing a larger number of host genotypes and previously identified as high-Aquarickettsia and low-Aquarickettsia (e.g., (Aguirre et al., 2022)) may further clarify factors that drive differential coral vulnerability to microplastics, acidification, and warming, especially since low-Aquarickettsia individuals exhibit higher resistance to disease, relative to high-Aquarickettsia individuals .
Genes implicated in the coral immune response were upregulated in response to the OAW and multistressor treatment, suggesting that some sort of external stimulus triggered this response. Microbiome taxonomic profiling did not reveal significant differences in microbiomes across treatments, so we suggest that further research should investigate the transcriptomic and metabolomic responses of the bacterial and viral fractions across these stressors. For example, it is possible that while Aquarickettsia remained dominant, it may have differentially regulated genes in response to OAW conditions and/or OAW+MP. Alternatively, the immune response may reflect host response to physical injury/insult by contact with the plastics (van de Water et al., 2015), but not necessarily microbial infection.
Another potential explanation for the observed immune response is a host response to viruses. Indeed, GO categories consistent with a viral response were enriched in OAW+MP treatments (i.e., cellular response to virus). Coral viromes have been implicated in directing holobiont response to environmental stress (Thurber and Correa, 2011;Leruste et al., 2012;Nguyen-Kim et al., 2015;Thurber et al., 2017). It is thought that viruses may play a role in coral bleaching under thermal stress through host lysis of algal symbiont cells that precedes bleaching signals (Lawrence et al., 2014;Grupstra et al., 2022). Additionally, previous work has demonstrated increased expression of anti-viral transcripts in hosted algal symbionts (Cladocopium, formerly Clade C) in response to thermal stress (Levin et al., 2017). Taken together, it is possible that the A. cervicornis immune response may suggest a virally driven response that precedes the classic bleaching response (Barshis et al., 2013;Dixon et al., 2020).There is also evidence that microplastics and the microbial plastisphere can facilitate survival of certain viruses (Moresco et al., 2021). It is plausible that increased levels of virus-coral interactions could elicit a host immune response, if not also facilitate viral infection of coral cells (Bowley et al., 2021;Loiseau and Sorci, 2022).
We recommend that future work investigating the effects of microplastic pollution on reef-building corals, especially in combination with global change stressors, specifically target the viral dynamics and diversity in the seawater relative to those associated with the coral holobiont. Additionally, assessing the microbiome and virome at several time points throughout stress challenge, coupled with holobiont gene expression and metabolome profiles, will be crucial for understanding the timeline and underlying causes of potential assemblage shifts that inform holobiont phenotypes. Overall, our results highlight the importance of applying multi-omic approaches to better understand the consequences of local and global stressors on reef-building corals.

Author contributions
SS, NF, and SD conceived the research. SS and NF conducted the experiment. KG, NK, AKH, AMH, and CB performed molecular lab work. KG, NK, and SD completed the gene expression analyses. CB, AKH, and AMH completed the 16S and ITS2 analyses. CB and KS interpreted the 16S analyses. CB, KG, and SD wrote the manuscript with help from all co-authors. Funding acquisition, project administration, and resources was led by SD, NF, and SS. All authors contributed to the article and approved the submitted version.

Funding
Funding for this project was facilitated by Boston University (BU) Biology Department start-up funds to SD and by keystone project funds to KG from BU's Kilachand Honors College. CB was funded through a BU Microbiome Initiative Accelerator grant. Funding for the experimental portion was supplied to NF by the UNCW College of Arts and Sciences Research Initiative Award and SS by the Explorers Club. KS was supported by the Institutional Development Award (IDeA) Network for Biomedical Research Excellence from the National Institute of General Medical Sciences of the National Institutes of Health under grant number P20GM103430.