Physiological and Transcriptomic Variability Indicative of Differences in Key Functions Within a Single Coral Colony

Polyps in different locations on individual stony coral colonies experience variation in numerous environmental conditions including flow and light, potentially leading to transcriptional and physiological differences across the colony. Here, we describe high-resolution tissue and skeleton measurements and differential gene expression from multiple locations within a single colony of Stylophora pistillata, aiming to relate these to environmental gradients across the coral colony. We observed broad transcriptional responses in both the host and photosymbiont in response to height above the substrate, cardinal direction, and, most strongly, location along the branch axis. Specifically, several key physiological processes in the host appear more active toward branch tips including several metabolic pathways, toxin production for prey capture or defense, and biomolecular mechanisms of biomineralization. Further, the increase in gene expression related to these processes toward branch tips is conserved between S. pistillata and Acropora spp. The photosymbiont appears to respond transcriptionally to relative light intensity along the branch and due to cardinal direction. These differential responses were observed across the colony despite its genetic homogeneity and likely inter-polyp communication. While not a classical division of labor, each part of the colony appears to have distinct functional roles related to polyps’ differential exposure to environmental conditions.


INTRODUCTION
Coral reefs remain one of the most biodiverse ecosystems on Earth (Knowlton et al., 2010), despite facing a suite of anthropogenic stressors with both local and global causes (Hughes, 1994;Hoegh-Guldberg et al., 2019). Most stony corals (order Scleractinia), the sessile builders of the reef system, grow as colonies of thousands of genetically identical polyps formed by budding, with colonies ranging from centimeters to meters in diameter. The reef-building corals are also host to intracellular symbiotic dinoflagellates whose photosynthetic activity restricts the coral holobiont to the top ∼100 m of the ocean and fulfills much of the host's carbon nutritional requirements . Nevertheless, active feeding on particles from the water column plays a significant role as well (Martinez et al., 2020), with active predation utilizing toxins produced and delivered by cnidocytes (recently reviewed by Schmidt et al., 2019). These fixed carbon sources ultimately power the production of the external calcium carbonate skeleton that forms the reef and persists after the animal has died.
While the individual polyps in most colonial corals are genetically identical, each polyp inhabits a distinct location within the colony, which can be morphologically complex (e.g., in the case of branching corals). Thus, it has been suggested that individual polyps within a coral colony are differentially exposed to environmental conditions such as water flow and light availability (Carpenter and Patterson, 2007;Chang et al., 2009), which have distinct effects on coral physiology. These micro-environmental conditions may play a role in feeding capability , reproductive success (Mass et al., 2011), photosynthesis (Dennison and Barnes, 1988;Mass et al., 2007), calcification rate (Dennison and Barnes, 1988), and waste removal (Mass et al., 2010). For instance, reductions in flow between the colony branches can trap food particles, making it easier for hunting tentacles of the polyps to efficiently pull them from the water column (Chang et al., 2009), and leading to downstream effects on metabolism (Gladfelter et al., 1989). Similarly, differences in light exposure due to the coral's branching leads to self-shading that can dramatically reduce light availability (Anthony et al., 2005), potentially leading to finescale spatial heterogeneity in photosynthetic parameters between individual polyps as well as between polyps and the coenosarc (Ralph et al., 2002). We hypothesize that such physiological differences will be underpinned by transcriptional differences. Plasticity in host gene expression has been observed at coral branch tips versus bases across multiple colonies (Hemond et al., 2014). However, it has not previously been linked, at high spatial resolution, to tissue and skeleton morphological measurements indicative of relevant environmental parameters and physiology within a single coral colony.
Here we combined high-resolution transcriptomic analysis with field measurements and tissue and skeleton physiological assays to describe intra-colony plasticity in the Indo-Pacific stony coral, Stylophora pistillata. We sought to answer two questions: first, are there differences in function by groups of polyps at different locations throughout the colony; and secondly, can we attribute those differences to specific location characteristics that are likely associated with varying microenvironments across the colony? Our whole-genome expression and tissue and skeletal morphology data are from a single coral colony at three locations along each branch taken from a different location within the colony. This high-resolution dataset enables us to resolve variability along three main axes: tip-junction-base, top to bottom, and the cardinal direction (Figure 1). Each of these axes is associated with changes in gene expression and tissue properties or skeleton morphology, which we suggest are due in part to light and current direction and strength. It appears that polyps in the interior of the colony are exposed to relatively stagnant conditions, whereas the colony periphery likely experiences higher flow, light exposure, chances for predation, and requirements for competition and defense leading to the dynamic tissue and skeletal morphology and transcriptional responses observed here. Our findings have implications for the interpretation of data from branches collected from different locations across multiple coral colonies. Specifically, when fragments are taken from multiple colonies for replication, considerations of the inherent differences across a colony should be included in sampling design. Our findings also shed new light on the potential drivers of plasticity of corals' responses to environmental stress and disease, and should be considered when reconstructing paleophysiology and past microenvironmental conditions from well preserved fossil coral colonies.

Sample Collection
Ex situ measurements were conducted on a 15 cm diameter colony of S. pistillata collected during the morning in April 2016 at 8 m depth under a special permit from the Israel Nature and Parks Authority in the coral nursery of the Inter-University Institute for Marine Science (IUI), Eilat, Israel using SCUBA ( Figure 1A). Just prior to colony collection, we took 70 photographs of the colony, ensuring full 360 • coverage. We then created a three-dimensional model of the colony from the photographs following the methods of Lavy et al. (2015). Briefly, we uploaded the photographs to 123D Catch online software (Autodesk, Inc.), performed image corrections, paired photographs at 10 • -20 • differences, and constructed a united XYZ coordinate mesh for all photographs. The generated mesh was visualized in MeshLab (Cignoni et al., 2008) to produce Figures 1B,C, and inset. After photographing it, the colony was brought to the lab in a tank filled with seawater and divided into individual branches and branch segments within 20 min of collection. Nine branches (numbered 1-9 in Figures 1B,C and Supplementary Material) were taken from three horizontal concentric rings (one "top" branch, four "middle" branches, and four "bottom" branches) and classified according to their relative height within the coral colony (Figures 1B,C). The eight middle and bottom branches also represent duplicates from each cardinal direction (N, S, E, and W). All branches were further divided into sections according to their position on the branch (tip, junction, and base) and flash-frozen in TRI-Reagent (Sigma) (Figure 1 inset) for differential gene expression analysis. Five additional branches (numbered 20-24 in Figures 1B,C and Supplementary Material) distributed around the colony were chosen for examinations of the tissue and skeleton, fragmented into tip, junction, and base segments, and then kept at -80 • C.

Oxygen Measurements
Dissolved oxygen (DO) in seawater between branches of five S. pistillata colonies located in the IUI coral nursery was measured in situ at 5-10 m water depth during the same morning as the removal from the coral nursery of the colony intended for high spatial resolution examination. This nursery is a pair of raised 4 × 3 m frame structures (part of one can be seen behind the coral colony in Figure 1A) ranging in depth from 5 to 10 m, so that conditions experienced by all well-spaced colonies are quite similar across the nursery. Measurements were taken as up to five radial transects per colony, from near the core of each colony, mid-way from the core to the edge of the colony, and between the tips at the edge of the colony (i.e., center, mid-way, and periphery), using 1.1 (tip diameter) × 40 mm (length) Unisense oxygen microsensors (version OX-N) mounted on a Unisense UnderWater Meter System. Prior to the measurements, the sensors were calibrated using 0% oxygen (100 mg Na 2 SO 3 /100 ml water) and 100% oxygen (oxygen-bubbled water) calibration points. Due to the small depth range over which the sampled colonies were mounted on the coral nursery, a constant temperature was assumed. The microsensors were moved through each colony at a controlled pace in order to keep the intra-colony water as undisturbed as possible and to prevent exchange of interbranch water within the colony and with the surrounding water. Each measurement was taken for a duration of 2 min to allow the oxygen microsensor readings to stabilize. All measurements for each transect were normalized to the outermost reading for that transect (Supplementary Table 1). DO concentrations across one representative transect were calculated according to the following equation to establish that such values were normoxic: DO µmol L = DO sat µmol L (mV sat − mV 0 ) ) × in situ mV measurement.

Microbial Counts
Microbial counts outside and within the coral colony were performed at the same time as DO measurements on samples collected in situ via radial transects on five different colonies growing at the same site. We conducted preliminary experiments using small aliquots of concentrated fluorescein dye injected to water between branches of S. pistillata skeletons, the same size as that used in our single colony assays, in a static laboratory aquarium to test the efficacy of extracting water from between branches and at colony peripheries with minimal intra-colonial water mixing. Briefly, we injected the dye on one side of a branch and then used a syringe to carefully extract various volumes of water, from 0.4 to 12 ml in intervals of 0.4 ml, from the other side of the branch, with five replicate experiments per extract volume tested. We then measured the fluorescence of the extracted water, subtracted the fluorescence background of the surrounding water, and represented the results as a percentage of the fluorescence from the injected fluorescein dye aliquot. This showed that we could extract up to 1.6 ml of seawater from between coral branches and from colony peripheries without pulling in water from the immediate environment of the surrounding branches (Supplementary Figure 1). Water flow during our sampling dive was noticeably very low. We therefore felt confident extracting 1 ml of seawater from between the branches near the core of the colonies ("in" position) and the tips of the branches at the peripheries ("out" position) of the in situ colonies. Upon completion of the dive, all water samples were immediately carried to the lab and fixed in cryovials in a final concentration of 0.0625% glutaraldehyde, incubated in the dark for 10 min, and then flash-frozen at -80 • C until analysis by flow cytometry.
Just prior to analysis, microbial count samples were thawed at room temperature and diluted 1:1 in ultra-pure water, and 2 µm-diameter fluorescent beads were added to each sample as an internal standard (Polysciences, Warminster, PA, United States). The samples were analyzed on a FACSCanto TM II Flow Cytometry Analyzer Systems (BD Biosciences). We first examined the natural fluorescence of the cells (chlorophyll and phycoerythrin pigments). We then stained the cells with SYBR Green I (Molecular Probes/ThermoFisher) according to the manufacturer's instructions and counted the total microbial population as well as cyanobacterial sub-populations. Data were acquired and processed with FlowJo software. Flow rates of the instrument were determined several times during each run and the average value for a sample-free test run was used for calculating cells per ml.

Tissue Removal
Tissue was removed from five frozen coral skeleton branches, each previously fragmented into tip, junction, and base segments, by airbrush with phosphate buffered saline (PBS) and homogenized using an electrical homogenizer (MRC, HOG-160-1/2) for 10 s. A sub-sample of each fragment's homogenate was stored for cell counts, chlorophyll extraction, and total protein measurements and the rest of the homogenate was centrifuged for 5 min at 5,000 rpm. The supernatant was retained for host protein and hemolysis assay. The pellet was observed under the microscope to quantify photosynthetic symbionts. All quantities were adjusted to fragment extraction volume and normalized to cm 2 of the fragment's surface area.

Chlorophyll Extraction
Two ml of tissue homogenate was filtered on a Whatman GF/C filter and incubated with 1 ml 90% acetone for 2 h at 4 • C. After incubation, the filter was manually homogenized and the solution was filtered through a 0.22 µm syringe filter into a glass cuvette. Spectrophotometric measurements were conducted on a NanoDrop (Thermo-Fisher) and chlorophyll a concentrations were calculated from light absorbance results based on the following equation (Ritchie, 2008): chl-a µg ml = −0.3319 (ABS630) − 1.7485 (ABS 647 nm ) + 11.9442 (ABS 664 nm ) − 1.4306 (ABS 691 nm ) .

Nematocysts and Photosynthetic Symbiont Counts
100 µl of tissue homogenate was used for nematocyst (supernatant) and photosymbiont (pellet) counts by hemocytometer (BOECO, Germany); cells were counted in five randomly selected fields per fragment (1 mm 2 each) on a Nikon Eclipse Ti-S Inverted Microscope System. Nematocysts were counted under white light with a contrast filter whereas photosynthetic symbionts were viewed by fluorescence at 440 nm excitation and 590 nm emission. NIS ELEMENTS (Nikon) software was used for the cell counts with automatic counting settings limited to cell diameters smaller than 15 µm and circularity set to >0.5 but <1.

Total Protein and Host Protein
We quantified homogenate and host total protein; sonication (Ultrasonic Atomizer Probe Sonicator) was used for further extraction of symbiont protein. Protein concentration of each sample was measured by bicinchoninic acid (BCA) assay (Pierce) against a bovine serum albumin standard curve at an absorbance of 540 nm.

Hemolysis Assay
Hemolytic assays were performed against O Rh positive human blood cells obtained from the Rambam/Yoseftal Hospital Blood Bank, as described previously (Primor and Zlotkin, 1975). Briefly, 2 ml of whole blood was diluted to a final volume of 15 ml in PBS (pH 7.4) and centrifuged at 3,000g for 5 min. The supernatant was removed, and the process was repeated with the pelleted erythrocytes until the supernatant was clear. Washed erythrocytes were then resuspended in PBS to a final concentration of 20% (v/v). For each hemolysis assay, 160 µl of homogenate supernatant was incubated with 40 µl of washed erythrocyte suspension (4% erythrocytes v/v) at 37 • C for 30 min in a water bath. At the end of the incubation, 400 µl of PBS was added, and the assays were centrifuged at 3,000g for 3 min. The supernatant fluid containing the hemoglobin released from lysed erythrocytes was transferred to 96-well microplates and the absorbance at 540 nm was determined on a spectrophotometric microplate reader (Perkin-Elmer). In addition, each experiment was normalized to a positive (100% hemolysis) and negative control (0% hemolysis) by incubating erythrocytes with DDW and PBS alone, respectively. HU 50 was defined as the amount of homogenate required to cause 50% hemolysis (Bartosz et al., 2008) in a dilution series of the coral extract.

Fragment Surface Area
Surface areas were estimated using the aluminum foil method (Marsh, 1970). Briefly, aluminum foil was molded over the skeleton of each fragment (without measuring skeleton connective areas between the divided fragments), carefully removed, and weighed. This was repeated three times on every fragment and the surface area was estimated using a standard curve of the derived relationship between foil area and weight.

Skeleton Cleaning for Assessments of Surface Area and Corallite Measurements
After tissue extraction, all fragments were cleaned of organic residues overnight in 3% sodium hypochlorite; the skeletons were then washed in deionized water and dried at 55 • C.

Skeleton Micromorphological Analysis
The distance between neighboring polyps and the surface area of the polyp was measured using a Nikon binocular microscope and analyzed using NIS-Elements (Nikon) software. For each coral fragment, at least four corallites located in a position parallel to the focal plane were chosen. Major and minor axes and circularity were measured for each corallite and the ratio of the two axes (major divided by minor) was calculated; a ratio of 1 indicates a circle and the larger the ratio, the more elliptical the shape. Surface area of each measured corallite was calculated from the two-measured axes, using the following equation: S = p × r major axes × r (minor axes) . Finally, distances between neighboring corallites were also measured to evaluate polyp density in each area. Corallites were qualitatively analyzed by scanning electron microscopy.

RNA Extraction, Processing, and Sequencing
Nine branches, each segmented into three portions along their axes, yielded 36 transcriptome libraries across the threedimensional space of the coral colony (Figure 1, colony locations in Supplementary Table 3, read count statistics in Supplementary Table 4). Total RNA was extracted from the holobiont in fragments stored in TRI-Reagent (Sigma) following the manufacturer's protocol, with some modification at the homogenization step. Briefly, samples frozen in TRI-Reagent were heated and centrifuged, and then bromochloropropane, at a ratio of 1:10, was added to the samples for separation. After incubation at room temperature for 10 min and centrifugation, RNA was purified from the clear phase using a Purelink RNA Mini Kit (Ambion) according to the manufacturer's protocol. The RNA was washed in 70% ethanol and then on-column DNase digestion was performed using a Qiagen RNase-free DNase Kit according to the manufacturer's protocol. RNA-Seq libraries were prepared using an in-house protocol at the Weizmann Institute of Science. Briefly, the polyA fraction (mRNA) was purified from 500 ng of total RNA per sample followed by fragmentation and generation of double-stranded cDNA. Then, end repair, A base addition, adapter ligation and PCR amplification steps were performed. Libraries were evaluated by Qubit (Thermo Fisher) and TapeStation (Agilent). Sequencing libraries were constructed with barcodes to allow multiplexing of all samples to be run in each lane. Approximately 578 million total high-quality 125 bp paired end reads (15.45 ± 1.5 million paired reads per sample) were sequenced on an Illumina HiSeq 2500 across three different lanes (i.e., each sample run in triplicate to remove batch effects).

Host and Symbiont RNA-Seq Differential Expression Analysis
Standard RNA-Seq quality filtering was conducted as described previously . Briefly, RNA-Seq reads were adapter-trimmed using cutadapt 1.15 1 , and then low-quality regions were removed with Trimmomatic 0.3 (Love et al., 2014). In order to detect correct taxonomic classifications within our transcriptome data, high-quality reads were mapped to all available NCBI and Reefgenomics [reefgenomics.org; ] genome-based proteome databases of Symbiodiniaceae species Symbiodinium microadraticum, Cladocopium goreaui, and Fugacium kawagutii [formerly Symbiodinium spp. clades A, C1, and F, respectively (LaJeunesse et al., 2018)], and Cnidaria as well as selected stramenopiles/alveolates/Rhizaria and Metazoa databases using Diamond (Buchfink et al., 2015). Top hits almost exclusively belonged to robust corals (mainly S. pistillata), and S. microadriaticum [formerly "clade A" (LaJeunesse et al., 2018)]. These results are expected since S. microadraticum was found to be a predominant symbiont in shallow water in Eilat . For symbiont RNA-Seq mapping, we further aligned the reads to the merged host genome assembly (NCBI GCA_002571385.1) and the symbiont genome assembly (NCBI GCA_001939145.1) using STAR (Dobin et al., 2013;Aranda et al., 2016;Voolstra et al., 2017). Across all branch fragments, 72-80% (∼11-16 million reads) were concordantly mapped to the host genome, and 2-10% (∼0.5-1.5 million reads) were concordantly mapped to the symbiont genome. In total, 83-84% of the reads were aligned to either host or symbiont genomes. Differential expression (DE) analysis was conducted using Bioconductor DEseq2 (Love et al., 2014), separately for the host and the symbiont genes, using a DEseq2 generalized linear model. Specifically, we tested: (1) the additive effect of branch position, ring height and cardinal direction; (2) the additive and interaction effect of branch position and ring height (excluding cardinal direction factor; Supplementary Tables 5, 6, 7). Note that the first model cannot include the four top branch samples (two tips, junction, and base), since the single top branch is represented by an incomparable cardinal direction ("center"); as such, these samples were excluded from the analysis of expression of putative toxins and known biomineralization-related genes. For the symbiont analysis, only the first model was used. Although DESeq2 is designed to analyze samples of variable total read counts, in order to avoid read normalization problems, we only selected symbiont genes whose 25% quantile of read count was at least 12 reads (with 12,735 symbiont genes remaining); further, only symbiont samples with >300 k total reads were retained. After filtering, the correct symbiont read normalization was indicated by a typical close-to-symmetric distribution of fold-changes around the DESeq2 MA plot horizontal axis (Supplementary Figure 2). NMDS analysis was conducted using metaDMS in the R Vegan package based on log 10 FPM values of all expressed genes. For the symbiont NMDS, since read counts were relatively low, pairs of tips from the same branch are represented together based on the sum of read counts for both tips; however, symbiont DESeq2 data were not pooled for differential gene expression analysis.

Branch Genetics and SNPs Analysis
Single nucleotide polymorphisms (SNPs) analysis was conducted for each sample using the Broad Institute-recommended RNA-Seq SNPs practice 2 . The pipeline includes four steps: STAR reads mapping, a pre-processing step aimed at removing various alignment biases, variant calling using GATK version 3.5, and finally variant filtration and annotation (DePristo et al., 2011). In order to exclude variant-call biases due to changes in transcript abundance between samples, we considered only SNP loci with coverage greater than 30 RNA-Seq reads in all tested samples. Identity By Decent (IBD) analysis was conducted using SNPRealte in R.

Host and Symbiont Functional Enrichment Analysis
Biological terms were assigned to genes based on Uniprot S. pistillata 3 , KEGG S. pistillata 4 , S. pistillata Trinotate annotations (Bryant et al., 2017), and Uniprot S. microadriaticum databases. Enrichment analysis was conducted in Bioconductor GOSeq (Young et al., 2012), which corrects for enrichment biases associated with correlations between gene size and DE significance. We also searched for functional enrichment using the score-based tool GSEA (Subramanian et al., 2005;Sergushichev, 2016), which may be more sensitive than p-value cutoff-based searches (such as GOSeq) when small non-significant changes in relatively large groups of genes are expected. In GSEA, we used log 2 fold changes as scores. Because many enriched terms are functionally related, we hierarchically clustered the biological terms based on pairwise distances between groups of genes: where a, b are two sets of gene ids, and xa and xb are DE genes ids from a and b, respectively. Terms trees were constructed using Bioconductor-ggtree (Yu et al., 2018). Terms enriched in ring height comparisons to the top of the colony should be viewed with caution as only one branch was taken from that classification. Only significant enrichment results are reported.

Statistical Analyses
Statistical tests of all tissue, skeletal, and environmental data described above were performed in RStudio (2019). If a parameter displayed both homogeneity of variance (Fligner-Killeen test) and normal distribution (Shapiro-Wilke test), a one-way ANOVA (chlorophyll/total protein, chlorophyll/photosymbiont cell, photosymbiont/surface area, total protein/surface area, host protein/surface area, nematocysts/total protein, nematocysts/surface area, corallite diameter, and corallite distance) was performed, with TukeyHSD post hoc analysis for parameters with significant differences between locations. Parameters with non-homogenous variance and/or non-normal distribution were analyzed using Wilcoxon's rank sum test (water column microbes, DO) or the Kruskal-Wallis rank sum test (chlorophyll/surface area, photosymbiont/total protein, hemolytic units/total protein, hemolytic units/surface area, and corallite axis ratio), with Dunn's test for post hoc analysis for parameters with significant differences between locations. Environmental parameters were pooled between all colonies and tested between measured locations as described above. Coral and photosymbiont tissue and skeleton morphological parameters were tested between branch axis locations across the colony (i.e., five replicates each of junctions and bases, and 10 replicates of tips); however, there was not sufficient replication at cardinal directions or ring height for statistical evaluation of parameters between these locational classifications.

Differences in Water DO and Microbial Counts Across Colony Structure
Because branching corals' three-dimensional structure can lead to different micro-environments within the colony (Kaniewska et al., 2008;Chang et al., 2009), we first aimed to characterize indicators of the environmental conditions experienced by polyps at different locations within the colony that could potentially affect their gene expression. Rather than making direct measurements of flow within S. pistillata colonies at the IUI reef, we instead measured several within-colony water parameters affected by flow. We focused on dissolved oxygen because its concentration can vary within coral colonies due to oxygen consumption by the coral holobiont (nighttime) (Haas et al., 2013) and microbes in the water (Wild et al., 2010), oxygen production by photosynthetic microbes of the coral holobiont (daytime) (Mass et al., 2010;Haas et al., 2013), and transport of water across the colony under higher flow conditions (Chang et al., 2009), affecting coral physiology (reviewed by Nakamura, 2010). Likewise, microbial counts represent an easy to quantify variable that can be used as a proxy of water stability, with lower flow conditions potentially supporting the retention of microbes actively growing between branches (Schiller and Herndl, 1989;Ritchie, 2006), or those chemotaxing toward the coral tissue (Garren et al., 2014). We measured higher oxygen concentrations in water collected close to the core (center; ∼250 µmol/L) of the corals, compared to samples collected at colony peripheries between branch tips (∼230 µmol/L; p < 0.05, Wilcoxon rank sum test; Figure 2A and Supplementary Tables 1,2), with locations at colony peripheries more likely to experience mixing with the overlying ocean water. Although there was variability in DO signal between colonies (Supplementary Figure 3), the "center" versus "periphery" differences held across all colonies. We also observed higher microbe densities in water collected from the colony interiors versus the immediately surrounding seawater (p < 0.05, Wilcoxon rank sum test; Figure 2B and Supplementary Tables 1,2). Both of these observations may be due to lower rates of water exchange at branch bases due to decreased water flow within the colony, as suggested by Chang et al. (2009). Lower mass transport has been observed within the less densely branched Pocillopora eydouxi whose overall colony structure is similar to the S. pistillata studied here ( Figure 1A) compared to the more densely branched P. meandrina (Hossain and Staples, 2020). These differences further support the notion that the environment differs along the coral branches from colony exterior to interior. Some polyps likely inhabit more or less favorable areas, which we observe being reflected in transcriptomic and tissue and skeletal differences across the colony, ultimately leading to, for instance, previously observed lower calcification rates (Dennison and Barnes, 1988) and higher pH in the diffusive boundary layer adjacent to coral tissue (Chan et al., 2016) under reduced flow conditions, and optimal prey capture in colonial polyps at moderate flow rates (Wijgerde et al., 2012).

Location-Specific Intracolonial Variability of Gene Expression
To assess how the differences in environmental conditions resulted in different tissue and skeleton properties and gene expression across a colony, immediately after measuring DO and collecting microbial samples we removed a single S. pistillata colony from the IUI reef. We then examined its branches both at different cardinal directions and heights above the substrate (i.e., ring height), with each branch further segmented along its axis from base to tip (Figure 1 inset and Supplementary  Table 3). The identified patterns in gene expression and tissue and skeleton properties therefore represent responses to microenvironments experienced by the colony at the specific time of collection and depth at which the colony was growing. The S. pistillata colony examined here was genetically homogenous (not chimeric, Supplementary Figure 4), and thus any changes in gene expression between regions were likely a response of the polyps to their local micro-environmental conditions rather than genotypic differences. While our study was limited to a single colony to achieve high spatial resolution of data, many similar genetic patterns are conserved between two major stony coral clades (see below). Additionally, differences were observed in hemolytic potential between branch tips and  Table 7) and marked by dashed and solid ellipses, respectively. bases, which were seen also in two additional colonies from an independent study (Ben-Ari et al., 2018). Taken together, it is clear that there are transcriptional and physiological differences between various locations within the coral colony, and that many of these differences are likely to be common across branching colonial corals. The expression of >10,000 of the ∼24,000 tested host genes was affected by at least one of the tested factors: location at branch, ring height, and cardinal direction (Supplementary Tables 4, 5, 6). We observed a clear differentiation in host gene expression between branch positions and from the top to the bottom of the colony (Figures 2C, 3 and Supplementary Table 5), although comparisons to the lone top branch should be viewed with caution. Combined, we observed an interaction effect between factors (Supplementary Tables 4, 6), as, for instance, groups of genes involved in signaling pathways important in development that were significantly enriched in branch tips versus bases as compared to middle versus bottom interactions ( Figure 3B and Supplementary Figure 5). Comparing colony interior with periphery, host gene expression patterns are more similar between junction and base branch positions compared to tips, although not the same at junctions and bases ( Figure 2C). While roughly half of all tested genes showed differential expression along the branch axis or by cardinal direction or ring height, the only indicators of physiological differences along the branch axis were hemolytic potential and skeleton morphology (described in detail below). Elucidating the microenvironmental drivers of these genetic and physiological differences across a colony may aid in, for example, understanding spatial differences of polyps' interactions with various infectious organisms that target the oral ectoderm (Gibbin et al., 2018) leading to diseases such as skeletal band eroding disease (Winkler et al., 2004), which has been observed to commence at colony interiors. Other previously reported response differences across colonies that may be due to microenvironmental variations include bleaching that commences at colony perimeters in Oculina (Shenkar et al., 2005), greater storage of lipids and total fatty acids toward colony interiors in Acropora (Conlan et al., 2018), and differing relative abundance of various photosymbiont species in regions that later displayed differential bleaching susceptibility (Kemp et al., 2014).
In contrast to the host, of ∼12,000 tested photosymbiont genes, only ∼150 were differentially expressed (DE) (Supplementary Table 7). The host versus symbiont differences are partly attributed to technical differences in statistical power due to lower symbiont read counts, and likely represent an underestimation of symbiont DE genes. DE gene patterns of  Table 9), with shaded regions representing genes showing the same DE patterns and noting "up" and "down" regulated genes. Enrichment patterns of clustered host KEGG (S. pistillata; B) and photosymbiont GO terms (UniProt S. microadriaticum; C) across factor comparisons. Host terms enriched in both S. pistillata (here) and two species of Acropora (Hemond et al., 2014) for tip versus base comparisons, are shown as green text. Hierarchical tree clustering reflects functional similarity between terms. Color-coded heatmap scores equal to log 10 (adjusted p-value) of enriched terms multiplied by 1 and -1 for up-and down-regulated genes, respectively (for clarity, values <-2 and >2 are shown as 2 and -2).
the photosymbionts are not as clearly delineated as for the host ( Figure 2D); nevertheless, the expression patterns differ between cardinal directions (p < 0.0007, ADONIS test, Supplementary  Table 8), likely reflecting the different light regimes experienced on different sides of the colony, visible as shading on the north (right) side of the colony and no shading on the south (left) side of the colony in Figure 1A, although we did not measure light levels around the colony. The ultimate drivers of the differences observed for groups of polyps in locations differing by ring height and cardinal direction remain to be resolved. Further, the observed differences across a single colony could have implications for studies using branches collected from different locations on separate colonies. While intracolonial differences that may diminish observable transcriptional and physiological differences between treatments-if chosen branches are from differing parts of selected colonies-may lessen over multi-year studies (Palumbi et al., 2014), epigenetic responses that differ across a colony may persist across generations (Liew et al., 2020).

Patterns in Gene Expression Are Conserved Between Coral Taxa
A previous study, analyzing the differences in gene expression between branch tips and bases in Acropora spp. (Hemond et al., 2014), provided an opportunity to identify conserved intracolonial patterns in gene expression between these two coral taxa, which likely diverged at least 200 million years ago (Park et al., 2012). We identified ∼9,600 putative orthologous expressed genes between S. pistillata and two species of Acropora spp. (Hemond et al., 2014), based on reciprocal BLAST. Interestingly, from a total of 952 genes that were DE in Acropora spp. along the branch axis, ∼70% were also DE in S. pistillata (672), and ∼90% of these exhibited the same fold change trend (p-value < 2.2e-16, Fisher's exact test; Figure 3A and Supplementary Table 9). Overall, these results indicate that micro-environmental dependence of gene differential expression at tips versus bases reflects an evolutionarily conserved rather than a species-specific pattern. Importantly, these results were observed in highly divergent species that are members of the robust (here) and complex (Hemond et al., 2014) scleractinian clades. Additionally, the conserved gene expression patterns were observed between taxa which possess physiologically unique axial polyps (Acropora; (Wallace, 1978(Wallace, , 1985 and a genus which does not (Stylophora). Our findings therefore support that the occurrence of intracolonial differential responses to microenvironments in general, and up-regulation of genes for specific processes such as protein production and energy generation in particular ( Figure 3B, green font), can likely be generalized across branching coral taxa. Further study is needed, however, to confirm our findings from this currently small dataset.

Gene Expression Patterns Suggest Differences in Host Versus Photosymbiont Metabolism
We next inferred potential broad-scale changes in cellular and metabolic functions by identifying KEGG and Gene Ontology (GO) enrichment patterns associated with micro-environmental changes in both the host (Figure 3B and Supplementary Tables 10,11) and the photosymbiont (Figure 3C and Supplementary Table 12). As expected, broad clusters of host KEGG terms and photosymbiont GO terms were DE (Figures 3B,C), mirroring the pattern observed for all DE genes within the host and symbiont (Figures 2C,D), with FIGURE 4 | Toxicity potential. Nematocyst abundance (A) and hemolytic units per area (B), and gene expression patterns of putative toxins as log 2 FC between branch tips (positive values) and bases (negative values) (C), and log 10 FPKM along S. pistillata branch axis (D). Only significantly DE toxin genes are shown (p < 0.05). The two most abundant transcripts, a cytolysin annotated as stichotoxin and an uncharacterized gene, for which reciprocal BLAST hits suggest it is a SCRiP (Supplementary Table 14), are noted with a red star (D). many more terms DE in the coral host along the branch axis compared to the ring height or cardinal direction (Figure 3B), and more terms differed along cardinal direction for the photosymbiont (Figure 3C).
Most DE host genes grouped by KEGG term tended to be up-regulated in branch tips relative to the junctions or bases ( Figure 3B and Supplementary Tables 10, 11). Specifically, many terms related to cell metabolism were enriched in the tips versus the bases, including carbohydrate metabolism (e.g., glycolysis and the TCA cycle, oxidative phosphorylation, mitochondrial biosynthesis), amino acid biosynthesis and metabolism, and lipid metabolism ( Figure 3B). We further tested which functional groups were enriched with conserved tip versus base DE gene patterns between coral taxa (Figures 3A,B). As indicated by green KEGG term text in Figure 3B, many fundamental processes such as metabolism show similar DE trends and enrichment between these two genera (Supplementary Tables 9, 10), indicating a conserved response between the taxa. "Zooming in" from the KEGG terms to the underlying genes, we note that genes related to oxidative phosphorylation, such as ATP synthase, cytochrome c oxidase, co-enzyme Q, and NADH dehydrogenase tend to be up-regulated at branch tips across the species. In addition, multiple genes for succinate dehydrogenase as well as other individual genes comprising the citric acid cycle, were also more highly expressed at branch tips than junctions or bases in S. pistillata and many of these showed similar regulation patterns in the Acropora spp. Finally, most of the >100 DE genes with KEGG terms for mitochondrial biogenesis were up-regulated at tips compared with junctions and bases in S. pistillata and many of their likely homologs were also up-regulated at branch tips in Acropora spp. The increased rates in physiological processes involved in energy production at branch tips, as suggested by our gene expression analysis, have been observed in the branching coral A. palmata as increased respiration and mitotic index (Gladfelter et al., 1989).
Intriguingly, the patterns of gene expression observed in the host were not all recapitulated in the photosymbiont. Although some symbiont GO terms associated with photosynthesis were enriched in the tips versus bases in the symbionts (e.g., those involved in light-harvesting and in the structure and function of the thylakoid), likely in response to higher light intensities at tips compared to bases Kaniewska et al., 2008), others were not, including terms involved with carbon fixation (Figure 3C). Additionally, no differing patterns in chlorophyll and symbiont density per both total protein and area were found along branch axis (Supplementary Tables 2,13 and Supplementary Figure 6), although we acknowledge that these are limited measures associated with photophysiological processes such as photosynthetic rate  and quantum efficiency of Photosystem II (Cohen and Dubinsky, 2015). Indeed, previous studies have shown that the products of photosynthesis can be translocated from branch bases to the tips in branching corals (Rinkevich and Loya, 1983). Rather than exhibiting tip-to-base functional differences, many of the terms associated with symbiont photosynthesis (both light and dark reactions) were found to differ along cardinal direction, being enriched in south-facing branches (Figure 3C), possibly due to the direction of the sun. As such, we anticipate that future studies will find differences in gene expression in the symbiont across a diel cycle, with, for example potentially decreased differences between north-and south-facing branches during nighttime hours. As noted above, these enrichment patterns are likely an underestimation. Furthermore, correlations between mRNA abundance and physiological patterns may not be high due to translation level effects, translocation of host/symbiont products, food particle capture by the host, etc., which may lead to differences in symbiont transcriptional control.

Developmental Pathways Are DE Along the Branch Axis
In addition to pathways related to metabolism, several KEGG terms associated with development and morphological plasticity in colony architecture were DE between branch tips and junctions or bases. For example, the Wnt signaling pathway was enriched at branch tips versus the bases in both S. pistillata ( Figure 3B) and Acropora spp. (Hemond et al., 2014). This pathway, which transmits information from outside to inside cells, is associated with cell differentiation and axis patterning in cnidarians (reviewed in Lee et al., 2006;Duffy et al., 2010) and bilaterians (Siegfried and Perrimon, 1994). Similarly, the TGFβ signaling pathway was also enriched in S. pistillata branch tips (Figure 3B), although this was not observed in Acropora spp. (Hemond et al., 2014); this pathway has been implicated for a role in development and calcification in corals (Gutner-Hoch et al., 2017). At the individual gene level, we observed that fibroblast growth factors, which could aid in development of limbs (reviewed in Xu et al., 1999) and in production of nerve cells (reviewed in Guillemot and Zimmer, 2011), were significantly up-regulated at branch tips. Half of the genes involved in mitochondrial biosynthesis were also up-regulated at branch tips. Further, two different bone morphogenetic proteins (BMPs), which may play a role in axial patterning in cnidarian polyps (Reinhardt et al., 2004) and in budding (Surekha et al., 2020), were up-regulated at branch tips here. Finally, several retinol dehydrogenases, proteins associated with photoreceptors in higher organisms (Albalat, 2012), were up-regulated at branch tips, where light levels tend to be higher, an observation shared with Acropora spp. (Hemond et al., 2014).
In addition to the plasticity described above, gamete production has also been shown to vary along branches axes across several coral taxa with gametogenesis typically occurring preferentially deeper in the colony (Stimson, 1978). We observed the highest expression of several testis-specific genes, such as testis-specific chromodomain protein Y 1-like and testis-specific protein kinases at S. pistillata branch junctions (Supplementary Table 4). In humans, these genes are associated with chromosome Y (Lahn and Page, 1999) and/or may be involved in sperm structure (Walden and Cowan, 1993). Interestingly, we observed up-regulation of octopamine receptors at branch tips, suggestive of higher maturation of ovaries and testes in that location (Chiu et al., 2020).

Changes in Expression of Genes Potentially Related to Predation Along Branch Axis
As described above, many metabolic processes were enriched in the tips versus the bases, whereas there is little evidence for higher photosynthesis by photosymbionts at branch tips.
Asking whether this could be due to more heterotrophic food input in the tips (e.g., through predation), we analyzed the tissue toxicity (hemolytic activity) and gene expression patterns of putative toxins. Although there was no difference in nematocyst abundance along branch axes (Figure 4A and Supplementary Table 2), tissue hemolytic activity was highest at branch tips ( Figure 4B and Supplementary Table 2), consistent with previous studies in S. pistillata that used branches from two different colonies (Ben-Ari et al., 2018). Thus, the observed branch axis patterning of hemolysis potential is robust across at least three colonies from two independent studies. Using reciprocal BLAST, we identified 22 putative toxin genes (Rachamim et al., 2015), encoding cytolysins, phospholipase A2 toxins, SCRiPs (small cysteine rich peptides) and a potential Kunitz/K channel toxin (Supplementary Table 14). We also compared these genes to the reciprocal BLAST results against the Acropora spp. dataset (Hemond et al., 2014), finding only two orthologous matches, a SCRiP and an uncharacterized protein. Like their S. pistillata counterparts, they were not DE between branch tips and bases (Supplementary Table 15). In our putative toxin dataset, we observed 13 and 11 DE genes between branch tips versus bases or junctions, respectively (Figures 4C,D; Supplementary Table 15 and Supplementary  Figure 7). Up-regulation of genes encoding potential toxins at branch tips has previously been observed in colonial corals (Hemond et al., 2014;Hemond and Vollmer, 2015). The two most abundantly expressed putative toxin genes encode for a cytolysin (up-regulated at tips) and a second SCRiP (up-regulated at bases), with the rest of the toxin genes having much lower overall expression levels (Figures 4C,D and Supplementary  Table 15). The highly expressed cytolysin is of particular interest as it belongs to the actinoporin family, and a hemolytic protein from this family, likely the protein or an isoform of the highly expressed gene here (Supplementary Figure 8), has previously been isolated from S. pistillata (Ben-Ari et al., 2018). Actinoporins are ∼20 kDa water-soluble proteins which are lethal when injected into crustaceans (Giese et al., 1996), and which insert themselves into membranes, particularly those containing sphingomyelins (Alvarez et al., 2009). They have been immunolocalized to tentacle nematocysts as well as to mesenterial filaments (Basulto et al., 2006). SCRiPs are ∼5 kDa peptides which may be anthozoan-specific (Sunagawa et al., 2009;Jouiaei et al., 2015), with evidence that they have neurotoxic effects (Jouiaei et al., 2015), although their function remains to be fully clarified.
The higher hemolytic activity and larger number of toxins more abundantly expressed at branch tips suggest that the tips FIGURE 5 | Skeleton characteristics and biomineralization gene expression. Microstructural differences of corallites and coenosteal spines along the branch axis (A; tip in i, ii; junction in iii, iv; base in v, vi). Corallite spacing (B) and ellipticity (C) along S. pistillata branch axis. Several known biomineralization genes were DE, shown as log 2 FC between branch tips (positive values) and bases (negative values) (D) and log 10 FPKM along the branch axis (E) to clarify genes with high versus low expression. Only significantly DE genes are shown (p < 0.05).
play the dominant role in venom-based competition, defense, and predation, although more work is needed to identify the role of putative toxins such as SCRiPs, which were more abundantly expressed at branch bases, and the potential contribution and roles of photosymbiont toxin production (Nakamura et al., 1993). While it is difficult to disentangle toxin effects on competition and defense from predation, based on higher flow at the edges of colonies (Chang et al., 2009), colonies' edges are likely exposed to higher food loads than are areas closer to the base of the branches. Consumption of these food sources could then be reflected in increased respiration rates as the food is metabolized, as has been observed in other species (Gladfelter et al., 1989) and may yield increased energy production which would require up-regulation of genes involved in this process as we observed here. We also queried genes related to digestion and observed that many peptidases were DE, but not in a clear pattern along the branches with approximately half upregulated at the tips. Similarly, multiple lipases were DE between branch locations, but no general trend was observed. However, pancreatic lipase expression, which has been observed in gastric cirri in a scyphozoan (Steinmetz, 2019), was found to be significantly up-regulated to a high degree at branch tips. We also observed up-regulation of three chitinases in branch tips versus bases, which could support extracellular digestion of arthropod prey (Steinmetz, 2019). In contrast, an amylase which converts starch and glycogen-biomolecules used by the photosymbiont and host, respectively, to store photosynthates (Kopp et al., 2015)-into simple sugars, was down-regulated at the tips versus junctions and bases. Thus, while our results suggest differences along the branch axis of functions related to prey capture and defense, including both toxicity and digestion, further work is needed to fully characterize the relative contributions FIGURE 6 | Micro-environmental variation across a colony leads to physiological and transcriptional differences. Schematic of environmental parameters (left) and resulting differential gene expression patterns and physiological outcomes for key processes focusing on those enhanced at branch tips (right). Higher flow is represented by thicker bands upstream of and above the colony, with vortices of lower flow rates within the colony. Differential lighting across the colony is represented by the shaded tissue on the side of the colony facing away from the sun. We cannot assign gene expression to specific tissue layers or cell types in most cases, except for the likely localization of toxin and biomineralization genes to the oral ectoderm and calicoderm, respectively. of photosynthesized products and predation. Additionally, temporal differences in predation capabilities at both the transcriptional (i.e., genes related to digestion) and physiological (i.e., hemolytic potential) levels across individual colonies should be examined as zooplankton abundance on coral reefs changes on a diel cycle (Yahel et al., 2005). Future studies measuring water flow rates explicitly in conjunction with transcriptional changes across colonies with complex macromorphologies will also clarify precisely how particle capture and digestion are affected by differing flow patterns .

Changes in Skeleton Structure and Expression of Known Biomineralization Genes Along Branch Axis
Shallow water stony corals are most notable for their production of massive reef structures, growing in size by adding mineral and new polyps, particularly for branching corals, at colony tips or edges (D'A and Le Tissier, 1988;Meesters and Bak, 1995). We therefore examined in depth skeletal morphological patterns and related gene expression across the three branch axis locations. Unlike most other parameters indicative of physiological processes across the colony which were not significantly different between branch positions (Supplementary  Tables 2,13), those linked to differences in biomineralization were significantly different along branch axes. Qualitatively, coenosteal spines appear narrowest and thinnest at the tip while base spines exhibit a more distinct microtuberculate texture in their lower portions (Figures 5Aii,iv,vi). Such microtuberculate texture with better differentiation of mineral fiber bundles may suggest higher organic matrix content which separates the mineral phase packages to a greater extent (Mass et al., 2014;Coronado et al., 2019). Corallites were closer together and more elliptically shaped at the tips than at the branch junctions or bases (Figures 5B,C and Supplementary Tables 2,13). Comparable trending toward ellipticity at branch tips has been reported in other Pocilloporids (Schmidt-Roach et al., 2014). Further, genetically identical S. pistillata branches have been shown to have more distantly spaced corallites with thicker theca and longer septa when grown under light spectra biased toward blue wavelengths (Rocha et al., 2014). While we did not measure light spectra in this study, such a spectral shift along a single branch, due to absorbance of red light by water, has been observed in S. pistillata tissue at the colony level on the sides versus at the top of ∼30 cm diameter colonies (Kaniewska et al., 2011) and at the branch level along Acropora millepora branches toward the colony interior (Wangpraseurt et al., 2014). We next examined differential expression of genes encoding proteins sequenced from S. pistillata skeleton and likely directly involved in the calcification mechanism (Puverel et al., 2005;Bertucci et al., 2011;Drake et al., 2013;Mass et al., 2014;Peled et al., 2020;Mummadisetti et al., 2021), as well as a bicarbonate anion transporter observed predominantly in the calicoblastic cell layer (Zoccola et al., 2015; Figure 5E; Supplementary  Table 16 and Supplementary Figure 9). We did not include genes for known skeletal proteins from Acropora millepora and A. digitifera (Ramos-Silva et al., 2013;Takeuchi et al., 2016), as not all orthologs have been found in S. pistillata skeleton although they exist in the genome, or non-orthologous proteins within the same functional group occur in S. pistillata skeleton (Zaquin et al., 2021). Similar to the global gene expression pattern for the host, we observed an up-regulation of biomineralization genes biased toward the tips (Figures 5D,E and Supplementary Table 16). Of the 36 biomineralization-related genes that were DE between branch tips and bases, eight orthologous genes were also DE and in the same direction in Acropora spp. (Hemond et al., 2014;Supplementary Table 17). In S. pistillata, the bicarbonate anion transporter SLCγ 4, a coadhesin-like protein, uncharacterized skeletal organic matrix protein 8 (USOMP8), and a hypothetical protein were up-regulated at branch tips but were expressed in only low amounts. In contrast, a ZP-domain containing protein, carbonic anhydrase (STPCA2), mammalian ependymin-related protein, and coral acid rich protein 3 (CARP3) were moderately up-regulated at branch tips and were very highly expressed overall. These proteins exhibit very different functions with SLCγ 4 and STPCA2 involved in regulating carbonate chemistry within the calcifying space (Zoccola et al., 2015(Zoccola et al., , 2016, while CARP3, which is a highly acidic protein [35% Asp, 15%Glu ] can direct mineralogy toward formation of Mgcalcite and vaterite Laipnik et al., 2020), and is up-regulated in newly-settled coral polyps (Akiva et al., 2018), in colonies growing at shallow depths , and in coral cell cultures reared at increased pCO 2 (Drake et al., 2017). Although the function of USOMP8 is not known, it is structurally similar to a BMP inhibitor . Of the remaining known highly acidic coral skeleton genes, CARPs 1 and 4 [CARP4 is sometimes also called skeletal aspartic acid rich protein 1; SAARP1 (Ramos-Silva et al., 2013)] were up-regulated at branch tips and were moderately expressed; both proteins are biased toward Asp residues  and are likely involved in nucleation initiation in early mineralization zones of growing septa (Mass et al., 2014). In contrast, CARP2 exhibited very low although significantly up-regulated expression at branch tips; while this protein is suggested to have a role in skeleton extension based on its localization to more soluble growth areas (Mass et al., 2014), its down-regulation upon settlement in pocilloporid larvae (Mass et al., 2016;Akiva et al., 2018), association with ACC phases pre-settlement (Akiva et al., 2018), and low overall expression in the present study may instead point to a role in inhibiting aragonite nucleation.
Previous and our present work provide strong genetic support for highest biomineralization activity at branch tips (Goreau, 1959). This faster calcification at the tips may result in the qualitatively thinner and narrower coenosteal spines we observed at branch tips ( Figure 5) that would be thickened by polyps as they age. Energy for higher calcification rates at branch tips could come from two sources. First, there is evidence of translocation of ATP from tissue regions deeper in the colony with relatively high photosynthetic rates (Fang et al., 1989) to branch tips. Secondly, our and previous data suggest that the potential for predation is higher at branch tips, with higher toxin gene expression and hemolytic activity at that location (Figure 4; Hemond et al., 2014;Ben-Ari et al., 2018), and up-regulation of genes related to metabolism and energy production (Figure 3). This further suggests that biomineralization gene expression and skeleton morphology patterns observed at the sub-polyp scale (i.e., centers of calcification in growing tips of septa) are recapitulated at the branch scale .

SUMMARY AND CONCLUSION
Some colonial cnidarians, particularly hydrozoans, exhibit a clear division of labor inside the colony, with functionally and morphologically distinct polyps across the colony (Cartwright et al., 1999). While no obvious division of labor at the polyp scale is observed in corals, our results show, from molecular to tissue and skeleton morphological scales on a single geneticallyidentical colony, a functional specialization by clusters of polyps due to their locations on coral branches and around the colony likely related to the microenvironments they encounter (Rueffler et al., 2012; Figure 6). In simplest terms, polyps at branch tips and toward the top of the colony are more exposed to light and higher water flow than are those found lower in the branches and closer to the substrate with consequent differential expression of a multitude of genes directing many key functions and, at least along the branch axis, differences in toxicity potential and skeleton morphology. Our findings suggest that branch tips may also be sites of higher potential for envenomation for colony defense and predation, as well as possibly for enhanced prey capture. These plastic responses by genetically identical polyps can then be communicated through the coenosarc to the rest of the colony to allow for enhanced resource sharing (Mackie, 1986). Such plasticity may aid in colony recovery after partial mortality events, such as those due to bleaching (Matsuda et al., 2020), and may be integral to adaptation to environmental change (Kenkel and Matz, 2016). Further, the expression patterns underlying the tested micro-environmental differences are found to be evolutionary conserved in the distantly-related stony corals S. pistillata and Acropora spp. Future study of multiple colonies and species will clarify the extent to which the general patterns we observe across the colony studied here are found across differing colony morphologies, diurnal and seasonal cycles, and geographic location. Further, the diversity of gene expression and tissue and skeleton properties across a coral colony needs to be taken into account when designing experiments with coral fragments from multiple colonies.

AUTHOR CONTRIBUTIONS
YP, OY, DT, DS, and TM conceived the study. YP, OY, and ES prepared the samples. JD, AM, YP, OY, ES, and JS conducted formal analysis of the data. All authors wrote the manuscript and approved of this submission.

FUNDING
This work was supported by grant 1239/13 from the Israel Science Foundation and by the associated Grant Number 2236/16 from the Israeli National Center for Personalized Medicine to DS. This work has received funding to TM from the European Research Council under the European Union's Horizon 2020 Research and Innovation Program (Grant Agreement Number 755876). JD was supported by the Zuckerman STEM Leadership Program. Computations presented in this work were performed on the Hive computer cluster at the University of Haifa, which is partly funded by Israel Science Foundation Grant 2155/15.