Urbanization comprehensively impairs biological rhythms in coral holobionts.

Coral reefs are in global decline due to climate change and anthropogenic influences (Hughes et al., Conservation Biology, 27: 261-269, 2013). Near coastal cities or other densely populated areas, coral reefs face a range of additional challenges. While considerable progress has been made in understanding coral responses to acute individual stressors (Dominoni et al., Nature Ecology & Evolution, 4: 502-511, 2020), the impacts of chronic exposure to varying combinations of sensory pollutants are largely unknown. To investigate the impacts of urban proximity on corals, we conducted a year-long in-natura study-incorporating sampling at diel, monthly, and seasonal time points-in which we compared corals from an urban area to corals from a proximal non-urban area. Here we reveal that despite appearing relatively healthy, natural biorhythms and environmental sensory systems were extensively disturbed in corals from the urban environment. Transcriptomic data indicated poor symbiont performance, disturbance to gametogenic cycles, and loss or shifted seasonality of vital biological processes. Altered seasonality patterns were also observed in the microbiomes of the urban coral population, signifying the impact of urbanization on the holobiont, rather than the coral host alone. These results should raise alarm regarding the largely unknown long-term impacts of sensory pollution on the resilience and survival of coral reefs close to coastal communities.

INTRODUCTION

Human activities have impaired normal ecosystem functioning across most of the Earth's surface (Dominoni et al., 2020), highlighted by the loss of biodiversity, even in remote areas (Motesharrei et al., 2016; Vitousek et al., 1997). Coastal ecosystems proximal to urban centers in the tropics are likely to be among the most vulnerable, as they face not only climate change but unpredictable and fluctuating nutrient (nitrogen and phosphorus) and xenobiotic (hormones and other organic contaminants) local pollution, besides chronic exposure to noise and light pollution (Duarte et al., 2021; Heery et al., 2018). Thus, urbanization of coastal areas near coral reefs is a global issue and is gaining momentum; not only has there been continuous growth of cities such as Jakarta, Singapore, and Hong Kong, but also major new developments have occurred or are planned that are likely to directly impact coral reefs in the near future. The population of the Chinese coastal city of Shenzhen has grown from <1 million in 1990 to >12.5 million in 2021, further expansion being predicted (https://www.macrotrends.net/cities/20667/shenzhen/population), and the rapid development of the Jeddah Corniche on the Red Sea coast of Saudi Arabia is of particular concern given the unique nature of the adjacent reef system (Kleinhaus et al., 2020).

Human activity has impacted animal and marine habitats in almost every conceivable way. This includes urbanization, buildings, lights at night, chemicals from industry or farming, tourism etc. which are known as anthropogenic or sensory pollutants. Sensory pollutants can mask environmental cues, interfere with the cellular processing of information, or alter cue perception leading to distracted responses by the organism (Dominoni et al., 2020). Therefore, formerly dependable cues may no longer be reliable in environments altered by humans. These anthropogenic stimuli can decrease animal survival and reproductive success and may ultimately alter populations and ecological communities. To understand and mitigate the effect of these stimuli, it is crucial to study the impact underlying the sensory reception of these pollutants, on marine habitats like coral reefs (Halfwerk & Slabbekoorn, 2015). Although light pollution has been shown to disrupt the timing of coral spawning (Ayalon et al., 2020; Loya, 2004), with potentially devastating consequences, the broader impacts of urbanization on corals are unknown, as are the underlying molecular mechanisms.

Coral reefs, which support the highest concentrations of marine biodiversity and provide essential ecosystem services to millions of people, are among the most impacted coastal ecosystems in response to human act (Hughes et al., 2013). One consequence of increased human activity near coastlines is that the community structure of many reefs has changed—in extreme cases to dominance by algae or other taxa rather than corals (Guest et al., 2016; McManus & Polsenberg, 2004). Coral reefs exposed to anthropogenic stress typically exhibit lower structural complexity, are dominated by “stress‐resistant,” generalist coral species, and exhibit decreased coral cover, all of which compromise ecosystem function (Brandl et al., 2019; Heery et al., 2018). Even though some corals have survived the selective pressures of “city life,” growing for years under chronic anthropogenic stress, they can face sudden fluctuations in pollutant levels to which corals are not adapted, thus having little chance to acclimatize.

Urban corals can be distinguished from corals in otherwise degraded reef ecosystems by differences in physiology, cellular processes, and growth characteristics (Heery et al., 2018; Nyström et al., 2000). Here we determined the effects of proximity to an urban environment on coral biorhythms over diel, monthly, and annual cycles, resolved in‐natura, using nitrogen and carbon stable isotopes, physiological monitoring, profiling gene expression in the coral host, and analyses of microbiome assemblage patterns. Our comparison between colonies from an urban area to a non‐urban area (Figure S1a) reveals the extent to which coral biorhythms are disrupted and modified by anthropogenic influences.

Despite being sessile organisms lacking specialized sensory organs, corals can sense chemical or physical environmental cues (Armoza‐Zvuloni et al., 2016; Paul & Puglisi, 2004) via complex repertoires of receptors that respond to the external cues by triggering signal‐transduction pathways to initiate specific biological processes (Armoza‐Zvuloni et al., 2016; Levy et al., 2007). The ability of these sensory systems to detect natural environmental cycles has been refined over millions of years of evolution. In contrast, anthropogenic influences are recent and can interfere with these mechanisms through which corals synchronize with the environment (Ayalon et al., 2019, 2020; Loya, 2004; Rosenberg, Doniger, & Levy, 2019; Shlesinger & Loya, 2019). For example, recent works (Ayalon et al., 2020; Shlesinger & Loya, 2019) clearly implicated light pollution in the impairment of coral gametogenesis and spawning synchrony. This study focuses on the common coral species Acropora eurystoma from the Gulf of Aqaba in the northern Red Sea to determine the impacts of urbanization on biological rhythms in corals. The fringing reefs in the Gulf of Aqaba are located at an unusually close distance, a few meters, from the shore (Loya, 2004). They are, therefore, particularly exposed to the impact of the surrounding urban environment (Ayalon et al., 2019; Loya, 1972, 2004; Loya et al., 2004; Rosenberg, Doniger, & Levy, 2019; Shlesinger & Loya, 2019). The General Circulation Model in the Gulf of Aqaba points on north to south current while a later model claim that coral larvae connectivity between “source” and “sink” in the Gulf can range between 9 ± 13 km (Berenshtein, 2018). Extrapolating to a global scale, the impacts of chronic exposure to sensory pollutants will likely further decrease the resilience of coral reefs but are not considered in current projections of the future of coastal coral reefs.

METHODS

Field observation and sample collection

Ten mature colonies of A. eurystoma were tagged in January 2016 in two separate environmental areas in the Gulf of Aqaba in the Red Sea at 5‐m depth. The urban zone is defined as heavily light‐polluted area in the northern part of the Gulf (29°32.678″ N, 34°58.204″ E) where we sampled five colonies. The second location, with five colonies and non‐urban conditions, is a dark at night area in the south of the Gulf (Figure 1a; 29°30.251″ N, 34°55.211″ E). Both areas were chosen based on prior light measurements made by Tamir et al., (2017) and represent two environments that are different in their urbanic threshold, while the coral reef biodiversity remains approximately the same. Although the two areas are only six km apart, the human and urbanization impact on the water quality and nutrients are well noticed (Figure S1b). Our sampling points were planned to cover different variables (Table S1):

FIGURE 1

Environmental and physiological parameters of the non‐urban and urban coral colonies. (a) Physiological assay parameters compared between the areas across the seasons, shown are mean (horizontal bars). Range of δ13C (b) and δ15N (c) values in algae (Symbiodiniaceae) and coral tissue (Host) in both areas across the seasons, shown are mean (horizontal bars). (d) Non‐metric multidimensional scaling (NMDS) with area x condition effect (preforming multivariate analysis of variance test), vectors showing response between the areas. p > .05, *p ≤ .05, **p ≤ .01, ***p ≤ .001, ****p ≤ .0001

Seasonal changes—sampling at all four seasons of the year (based on sea temperature change) was chosen to represent the yearly variations in gene expression and physiology parameters. Therefore, we sampled during February (winter), June (spring), August (summer), and November (fall).

Moon phase—at every month of sampling we sampled at two opposite moon phases to present the changes in gene expression that vary between moon phases. We sampled during the full and new moon of every one of the four seasons.

Diurnal changes—at every sampling day (full moon and the new moon of all four seasons, a total of eight sampling days) we sampled at 12:00 p.m. and 21:30 p.m. to see changes in gene expression varying across a solar day. These two sampling points across the day represent the strongest illumination from the sun (12:00 p.m.) and the moon (21:30 p.m.). Sampling was conducted by SCUBA diving. At each sampling point a small branch from each colony, measuring an average of five cm in length, using pliers, was collected. Each sample was inserted into a marked tube filled with seawater. After all the colonies were sampled, the tubes were transferred to shore and placed in a small piece of aluminum foil with a tag containing the sample ID, time, and colony. The small aluminum bundle was snap‐frozen in liquid nitrogen and transferred to a −80 freezer. In total, there were 16 sampling points across the year taking place simultaneously at both areas, sampling fragments from five colonies per area. The environmental monitoring data was retrieved from The Israel National Monitoring Program at the Gulf of Eilat (http://www.iui‐eilat.ac.il/NMP).

Physiology

Coral fragments from the full moon of each month of sampling (n = 20 samples per area) were tested for protein concentration, algae density, and chlorophyll concentration to assess coral health and yearly changes in these parameters. Tissue was removed from frozen coral fragments using an airbrush and ice‐cold filtered (0.22 μM) seawater (FSW). Skeletons were retained for surface area determination using the wax dip technique (Stimson & Kinzie, 1991). Tissue samples were homogenized for 30 s using an electrical tissue homogenizer (Kinematica Polytron™ PT2100 Benchtop Homogenizer). A sub‐sample (100 µl) of the supernatant was taken to determine host protein concentration by a colorimetric method using a multi‐scan spectrum spectrophotometer (595 nm, 450 nm, Biotek HT Synergy plate reader) and bovine serum albumin as a standard (Quick Start Bradford Protein Assay, BIO‐RAD). Protein concentration was used as a biomass and normalization index for the coral fragments. Samples were centrifuged to separate host and symbiont tissues. Further centrifugation and washing with filtered seawater were performed to isolate symbiont cells for cell counts (hemocytometer) and photosynthetic pigment extraction. Pigments were extracted for 24 h in 90% acetone at 4°C in the dark and Chlorophyll (Chl) a and c 2 concentrations were measured spectrometrically at 630, 664, and 750 nm with a Multiskan spectrum microplate spectrophotometer (Thermo Fisher Scientific). Chlorophyll concentration was determined as previously described and normalized to algae cells and surface area. Algae were counted with a hemocytometer under a microscope and normalized to the coral surface area. Each result is an average of five fragments per treatment. Physiology assay results from both treatments were compared to each other in each parameter using the r package “unpaired two‐sample t‐test” to represent statistical relevance. A comparison of both areas across all seasons was analyzed using one‐way ANOVA and Tukey HSD.

Stable isotopes

The isotopic measurements were made at the stable isotopes laboratory in the Department of Earth and Planetary Sciences, the Weizmann Institute of Science, Israel. The remaining homogenate from the physiology assays (n = 20 samples per area) was centrifuged to separate the algae from the host tissue at 3000 g for 5 min at 4°C. Supernatant containing host tissue was transferred to a new falcon tube. The remaining pellet containing algae was re‐suspended in 1 ml of double‐distilled water (DDW) and transferred to a new Eppendorf tube. Host tissue samples were centrifuged for a second time at 3000 g for 5 min at 4°C to eliminate all algae tissue, followed by a centrifuge at maximum speed (13,500 rpm) for 10 min at 4°C. Remaining supernatant was discarded and host pellet was re‐suspended in 1 ml DDW and transferred to a new Eppendorf tube. To all samples (host and algae) 0.1 ml of HCL (1 N) was added, short vortex and left for 10 min at room temperature. Samples were centrifuged at maximum speed (13,500 rpm) for 5 min at 4°C. The supernatant was discarded, and pellets were re‐suspended in 1 ml DDW. Samples were centrifuged at maximum speed (13,500 rpm) for 5 min at 4°C, the supernatant was discarded and re‐suspended in 1 ml DDW for two more times. After the final centrifuge, the supernatant was carefully removed using a small pipet, and samples were lyophilized. For Isotope analysis, a small fraction (190–270 µg) from each sample (n = 40) was weighed and placed in a small tin cap. Dried tissue and algae were analyzed using an elemental analyzer (CE 1110) interfaced to the MAT 252 mass spectrometer. Long‐term precision of working standards for ™13C is 0.05‰ and for ™15N is 0.1‰ relative to V‐PDB and air, respectively (±1σ SD). Presented statistical relevance was performed by using one‐way and two‐way ANOVA coupled with Tukey HSD.

RNA extraction and RNA‐Seq library preparation

A total of three samples per area from every sampling point was used for RNA analysis (n = 96 samples). Total RNA was extracted using TRIzol reagent (Invitrogen) and a modified version of the manufacturer's protocol that included an additional chloroform extraction and Lithium Chloride precipitation overnight. This method was chosen as we targeted coral host RNA. Purified RNA samples were analyzed using a NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific) to assess RNA quantity and a 2100 Bioanalyzer (Agilent) to assess RNA quality (RIN >8.5). One and a half microgram RNA from each of the 96 samples was sent for sequencing at the Technion Genome Center in Haifa, Israel. RNA samples were prepared using the Illumina TruSeq RNA Library Preparation Kit v2, according to the manufacturer's protocol. Libraries from each sampling point ran on three lanes of an Illumina HiSeq2500 using the multiplexing strategy of the TruSeq protocol. The protocol starts with polyA selection that results in mRNA selection only. On average 35 million read pairs, 100 bp long, were obtained for each sample.

Transcriptome assembly

de‐novo transcriptome assembly was based on the paired‐end (PE) sequences of all 96 samples. Sequences were first trimmed from Illumina adapters using cutadapt (Martin, 2011) and then quality‐filtered using Trimmomatic (Bolger et al., 2014). The resulting sequences were assembled in Trinity 2.4 (Grabherr et al., 2011), which included 1,065,622 Trinity transcripts (Total assembly length = 995,608,831 bp; N50 = 1785 bp). Quantitation of transcript abundances was done using Kallisto‐0.44.0, and a further quantification at the gene‐level was done using tximport (Soneson et al., 2016) on R3.4.2. Among all Trinity genes, Kallisto/tximport pipline identified a total of 191,577 genes expressed with >2 reads per Trinity gene cutoff. The expressed Trinity transcripts were mapped to SwissProt database provided by Trinotate and confirmed Metazoan Trinity genes (100,568 genes, total assembly length = 128,952,937 bp; N50 = 2174 bp) were retained for subsequent analyses. We further selected for Trinity genes that were aligned to protein sequences curated in the UniProt SwissProt protein database using BlastX with a e‐value threshold of ≤1e−6. This analysis included 10,000 Trinity transcripts were retained and used for analysis. In order to assess representation of the symbiont expression in the data, we examined mapping of the sequencing data against the draft genome assembly of Cladocopium goreaui (Clade C; http://symbs.reefgenomics.org/download/) using the STAR aligner (v2.7.8a) for two randomly selected samples. Mapping percentages were 12.6% and 16.5% corresponding to 1.8 and 2.6 million read pairs, respectively, which we considered insufficient for differential expression analysis.

Differential expression analysis

Differential expression analysis was conducted based on transcript‐level quantification using r version 3.4.2 DESeq2 (Love et al., 2014). Urban and non‐urban subsets were first analyzed separately using DESeq2 GLM method, considering the additive effect of the factors: coral colony, season, moon‐phase, and diel cycle. We subsequently compared the two subsets (urban and non‐urban corals) to each other considering the additive effect of area, season, diel cycle, and the interaction effect of the area and season. When comparing the two subsets to each other the effect of the colony cannot be calculated since all of the colonies were specific to only one area. Genes were considered differentially expressed only when having a Benjamini–Hochberg adjusted p‐value < .05. Trinity contigs were mapped to GO (Gene Ontology [GO]) terms and NCBI accessions using Blast2GO version 5, based on Blastx against NCBI‐nr Metazoa subset, and protein InterPro domain search (prosite.expasy.org). GO functional enrichment analysis was conducted using GOSeq (Young et al., 2010). GO terms were considered enriched where Benjamini–Hochberg adjusted p‐values < .05. To cluster enriched GO terms based on shared Trinity contigs, the Jaccard index was calculated between GO terms (only differentially expressed Trinity contigs were included). The resulting matrix was used for calculation of hierarchical tree using “hclust” command using UPGMA algorithm in r3.4.2. To resolve clusters of GO terms the tree was cut at three different scales of branch heights (0.20% height, 40% height, and 60% height) using r package “dendextend” function “cutree,” and then clusters were denoted based on functional similarity among GO terms.

16S rRNA library preparation and sequencing

A total of 10 µg of total RNA from each of the 96 samples, subjected to RNA‐Seq, was dried in GenTegra‐RNA 0.5 ml screwcap microtubes following manufacturer's instructions and sent for 16S amplicon library preparation at the University of Konstanz, Germany. Samples were recovered according to the manufacturer's protocol and residual DNA contamination was removed using the RNase‐Free DNase Set kit (Qiagen). Synthesis of cDNA was carried out using the High‐Capacity cDNA Reverse Transcription Kit (Thermo Fisher Scientific) following manufacturer's instructions. Amplification of the 16S rRNA was done using the Qiagen Multiplex PCR kit, with 1 μl cDNA as a template and using primers 784F (5′‐AGGATTAGATACCCTGGTA‐3′) and 1061R (5′‐CRRCACGAGCTGACGAC‐3′) with unique 8‐mer barcodes at their 5′ ends at a final primer concentration of 0.5 μM in a reaction volume of 10 μl (Andersson et al., 2008; Bayer et al., 2013). In addition, we did a null template (no DNA input) negative control reaction to assess for PCR reagent contamination. Thermal cycler conditions were 95°C for 15 min, 27 cycles of 95°C for 30 s, 55°C for 90 s, and 72°C for 30 s, followed by a final extension step of 72°C at 10 min. To confirm successful amplification, 2 µl of each PCR product were run on a 1% agarose gel. Samples were cleaned using ExoProStar 1‐step (GE Healthcare) and normalized using the SequalPrep Normalization Plate Kit (Thermo Fisher Scientific). Sample number 94 (fall, day, new moon, Kisosk) did not yield any amplification. The 95 coral samples plus 1 negative control were then pooled in two 1.5‐ml Eppendorf tubes (46 samples per tube, 10 μl per sample) and concentrated using a speedvac (Concentrator plus, Eppendorf). Quantification was done using Qubit (Qubit dsDNA High Sensitivity Assay Kit; Invitrogen). Samples were PE sequenced (2 × 250bp) on the NovaSeq 6000 platform at the Novogene Sequencing Centre (Cambridge, UK).

Bacterial community analysis

Demultiplexed and adaptor‐free sequences were used to infer amplicon sequence variants (ASVs) using DADA2 v.1.18.0 (Callahan et al., 2016). Reads were truncated at the 3′ end of forward and reverse reads at 200 base pairs and reads with an expected error >2 or with ambiguous bases were discarded (Supporting Information data file—bacterial ASV inference). After merging reads, ASVs were checked for chimeras and subsequently annotated using the SILVA database (version 138; Quast et al., 2013; Supporting Information data file—bacterial ASV abundance). Relative abundances of the most abundant 10 families were generated using ggplot2 v.3.3.3 (Wickham, 2011). Principal component analyses were plotted using Euclidian distances of centered log‐ratio (CLR) transformed ASV counts using Phyloseq v.1.34.0 (McMurdie & Holmes, 2013). Differences between bacterial communities were determined using PERMANOVA implemented in Vegan v.2.5 (Dixon, 2003) and pairwise PERMANOVAs were done using Pairwise.Adonis v.0.4 (Supporting Information data file—bacterial beta diversity). Differential abundance analysis of bacterial ASVs was done using ANCOM‐BC v.1.0.2 (Lin & Peddada, 2020). The false discovery rate method was used to correct p‐values for multiple comparisons. An ASV was considered differentially abundant if adjusted p‐values were <.05 between comparisons tables—ANCOM output and differentially abundant ASVs. Venn diagrams of differentially abundant ASVs were generated using VennDiagram v.1.0.12 (Chen & Boutros, 2011) and heatmaps of CLR‐transformed counts were plotted using pheatmap v 1.0.12 (available from https://github.com/raivokolde/pheatmap). We further looked for pathways involved in nitrogen and carbon fixation as well as photosynthesis in the different bacterial groups. Presence and abundance of gene orthologues from the Kyoto Encyclopedia of Genes and Genomes (KOs) were predicted from 16S rRNA amplicon data using PICRUSt2 (Douglas et al., 2020). Abundance of KEGG pathways between urban and non‐urban sites was compared using ANCOM‐BC as described above.

RESULTS

Coral physiology and isotopic signatures

The analyses presented here are based on a year‐long, in‐natura, experiment initiated in 2016 and involved sampling A. eurystoma colonies at 5 m depth from different environmental areas showing diverse characteristics at the northern end of the Gulf of Aqaba. The urban area is exposed to salinity fluctuations, periodic nutrient (N and P) enrichment and high human activity, inflicting light and noise pollution (Figure S1b; Armoza‐Zvuloni et al., 2016; Loya, 2004; Tamir et al., 2017). The non‐urban area is approximately 6 km south and is largely devoid of tourist and recreational activities. Sampling took place four times a year (winter, spring, summer, and fall) aligned with the seawater temperature cycle, twice a month (full and new moon), and twice a day (noon and 21:00). Physiological assay results (Figure 1a) showed no significant differences in algal symbiont cell density or total protein concentration between the two coral communities across all four seasons (p = .087). However, chlorophyll c 2 and chlorophyll a levels differed significantly between seasons and areas (p < .001), with urban samples having continuously higher chlorophyll concentrations compared to samples taken from the non‐urban area. Isotopic signatures measured in both host and symbiont tissue (p < .001) revealed higher carbon (δ13C) and nitrogen (δ15N) isotope signature values in samples from the urban area (Figure 1b,c, respectively). Non‐metric multidimensional scaling (Figure 1d) plots summarize the variables measured.

Diel, monthly and seasonal temporal gene expression

The physiological impacts of urban proximity were investigated across four seasons, two moon phases, and at two times of the day for microbiome profiling and host gene expression analyses. As is evident in Figure 2, the diel, monthly, and annual cycle gene expression patterns showed higher similarity in samples from the same area (n = 10,000 genes assessed), highlighting the dominance of sampling environment as a factor influencing gene expression patterns in each case (Figure 2). Differential expression analysis (Supporting Information data file, “gene expression” tab) disclosed that in samples from the non‐urban area ~7% of the genes (701 genes) responded to the diel cycle and changed significantly between night and day while only ~4% (406 genes) of the genes from the urban area samples responded to the diel cycle (p adj < .05; Figure 2a). We found that the expression of ~5% of the genes (477 genes) in the non‐urban and ~0.9% of the genes (89 genes) in urban corals changed according to the monthly cycle of the lunar phase (p adj < .05; Figure 2b). The largest transcriptomic response at both areas was associated with seasonality changes. Seasonal analysis was performed by comparing each season to the other three seasons separately (resulting in six comparisons) to reveal unique gene enrichments exclusive to only one season. An average of ~18% of genes (1771 genes) in non‐urban samples and ~14% of genes (1384 genes) in urban samples were differentially expressed when comparing two different seasons (p adj < .05; Figure 2c–h).

FIGURE 2

Scatter plot highlighting the transcriptomic response between non‐urban and urban corals across diel, monthly, and seasonal cycles. The log2 fold change in the non‐urban corals is represented on the x‐axis. The y‐axis shows the log2 fold change in the urban corals. Each gene is marked as a dot; unique differentially expressed genes of urban corals are highlighted in dark purple, unique differentially expressed genes for non‐urban corals highlighted in dark green, and significantly differentially expressed genes in both areas: in light blue when expressed in the opposite manner, in yellow if expressed in the same manner in both areas and the non‐significant gene were labeled as grey dots. Gene comparison is between the areas at the different cycles. (a) diel cycle, (b) monthly cycle and (c–h) seasonal cycle

Our data revealed that most diel processes, including developmental and growth processes, signal transduction, and response to outer stimulation, were enriched during the night and observed only in corals from the non‐urban area (Figure 3a; Figure S2a). By contrast, diel processes associated with endoplasmic reticulum (ER) stress were enriched only in corals from the urban area. Enrichment of genes associated with cell proliferation, transcription, and nitrogen metabolic processes on a diel cycle were observed in corals from both environments. Non‐urban corals followed diel cycles of gene expression resulting in developmental processes, signal transduction, and cellular growth enriched during the night, as previously documented in other corals (Levy et al., 2011; Maor‐Landaw & Levy, 2016; Oren et al., 2015). In urban corals, the process of protein mannosylation was enriched during the night, whereas lipid metabolism was enriched during the day (Figure 3a). The temporal uncoupling of these two processes can lead to ER stress, thereby engaging the unfolded protein response resulting in reduced global protein synthesis and increased apoptosis (Hetz, 2012; Høyer‐Hansen & Jäättelä, 2007; Schröder & Kaufman, 2005). Redox cycling between day and night in symbiotic corals is extreme because the algal symbionts impose photosynthetic rhythms driving diurnal cycles of hyperoxia (day) and anoxia (night) on the underlying circadian clock of the coral host, resulting in complex patterns of gene expression (Alderdice et al., 2021; Levy et al., 2011; Sorek et al., 2018). Although the putative clock genes Cry1 and Clock exhibited similar expression in day/night comparisons in corals from both areas, the diel patterns of gene expression differed between the coral groups. Corals from the non‐urban area presented upregulation of HIFα and several of its targets, including the glycolytic enzyme phosphofructokinase, during the night, following the predicted response to hypoxia stress (Alderdice et al., 2021). In contrast, we found no evidence of oxygen tension cycling in urban corals.

FIGURE 3

Gene Ontology enrichment analysis and underlying biological processes based on significant differentially expressed genes in each cycle. (a) Scatter plot representation of biological processes enriched in the diel cycle and (b) between moon phases. Scatter represent the relationship between biological processes and the Log2 fold‐change of each process. The number of genes clustered to each mother term process is represented through circle size

During the lunar cycle, a variety of biological processes, including cellular growth, organization, adhesion, polarity, movement, migration, oogenesis, protein de‐phosphorylation, actin and myosin movement, developmental processes, and cytoskeleton organization were enriched in non‐urban corals, whereas only genes involved in cytoskeleton organization were enriched in urban corals (Figure 3b; Figure S2b).

Elucidating seasonal patterns, GO enrichment analysis revealed common patterns at both coral groups for several biological processes, including cell morphogenesis, cellular communication, differentiation, and localization, signaling and developmental processes, intracellular signaling, amino acid metabolic/catabolic processes, RNA processing, cellular respiration, and cellular response to stimuli (Figure 4a; Figure S2c). However, seasonal shifting was observed for other biological processes, such as transcription regulation, immune response, Ras signaling, and RNA splicing, which were found enriched at both areas, but during different seasons. In addition, several biological processes were enriched exclusively in urban or non‐urban corals (Figure 4a). Urban corals exhibited enrichment of genes associated with the heme biosynthetic process, ATP metabolic process, and cell redox homeostasis. Non‐urban corals were enriched for sexual reproductive processes, the endoplasmic‐reticulum‐associated protein degradation (ERAD) pathway, chromatin assembly, cellular biosynthetic processes, cellular localization, sensory system development, and calcium ion signaling (Figure 4a). Corals, like other animals, typically display seasonal changes in reproductive processes, growth and DNA repair, as longer days bring the potential for higher levels of DNA damage. GO term analysis (Figure 4a) showed that the seasonal gene expression patterns associated with reproduction, cell cycle events, biosynthetic processes, and chromatin assembly were disrupted only at the urban area. For example, disruption of the reproduction process can be seen by the expression of the major egg protein vitellogenin (VIT and VIT6) that was altered earlier at the urban area than at the non‐urban. In addition, DM‐domain genes (i.e., doublesex/mab‐3 related genes) have been implicated in sex‐determination processes across the animal kingdom (Matson & Zarkower, 2012). Consistent with a role in reproduction, we identified a homolog of human DMRT3 that changed in expression in winter/spring and spring/summer comparisons at the non‐urban area, whereas at the urban area, no differential expression was apparent. Genes involved in RNA splicing, including homologs of Dicer, DCL1 (Dicer‐like 1), argonaute, MET25, and LIN41 also displayed seasonally‐shifted expression between the urban and non‐urban areas. Coral skeletal growth differs over seasons (Barnes & Lough, 1989; Goreau & Yonge, 1977; Guzmán & Tudhope, 1998; Roik et al., 2016) and many of the genes implicated in calcification showed altered expression at the urban area, including the known (Ramos‐Silva et al., 2013) skeletal organic matrix proteins hephaestin‐like, SAAR2, USOM1, USOM5, USOM6, USOM7, and galaxin. Season patterns of regulation of the ERAD system (seen at the urban but not at the non‐urban area) presumably reflect responses to higher stress exposure resulting in increased cellular damage.

FIGURE 4

Temporal seasonal biological processes enrichment pattern at each sampling area during the sampling period

Biorhythmicity of the coral microbiome

Evident disturbances in the microbiome were observed at the urban area (n = 18,897 ASVs were considered for PERMANOVA). Seasonal effects accounted for a greater proportion of the total variation (8.3%) than area differences (3.9% of the total variation), with moon phase and diel cycles accounting significantly less (Figure 5a,b). Microbiomes exhibited stronger seasonality patterns during both full and new moons at the non‐urban area, evidenced by the generally higher dissimilarities between seasons (Figure 5c). Differences between full and new moon microbiomes were consistently larger in the non‐urban compared with the urban area corals across seasons (Figure 5d). In contrast to moon cycles, no overall changes in the bacterial community were observed between day and night at any of the areas or seasons. We aggregated the abundances of differentially abundant ASVs to the family level to elucidate taxonomic patterns over diel and moon cycles. Most of the bacterial families that responded to diel and moon cycles were common between non‐urban and urban area corals; there was a temporal asynchrony in the response from each area. For example, members of the bacterial family Rhodobacteraceae peaked during the full moon and day in winter at the non‐urban area. In contrast, at the urban area, they peaked during the summer full moon and day (Figure 5e,f). The seasonal patterns based on differentially abundant ASVs were more defined in the non‐urban area than the urban area. Moreover, only a few ASVs were consistently enriched in full moon and day samples across all seasons (Figure S3). At the non‐urban area, ASV2433 (Unclassified Rickettsiales) was consistently enriched during the full moon, and ASV0120 (Chlamydiales—Simkaniaceae) and ASV0185 (Unclassified Rickettsiales) were consistently enriched during the day. In the non‐urban area, ASV0163 (Rhizobiales—Beijerinckiaceae), and ASV0237 (Vibrionales—Vibrionaceae) were consistently enriched during the day. When investigating the trophic level of corals from both areas by analyzing the microbial community and looking for pathways affiliated in nitrogen and carbon fixation and photosynthesis we found that in the non‐urban samples genes involved in N fixation and photosynthesis (Figure S4a–c) were more abundant than in the urban samples. In addition, ABC transporter genes (Figure S4d–f) were more lavish in urban samples, which could indicate toward different metabolic pathways preferred by corals at each site. Moreover, our predicted functional profiles based on bacterial composition showed that genes involved in nitrogen and carbon fixation, as well as photosynthesis, were more abundant in the non‐urban site (Figure S4a–c), whereas genes encoding ABC transporters were more abundant in the urban site (Figure S4d–f).

FIGURE 5

Coral microbiome activity changes over diel cycles and moon phases. (a) Principal component analysis (PCA) of bacterial communities over seasons and areas. PCA is based on CLR‐transformed amplicon sequence variant (ASV) abundances using Euclidean distances. (b) Overall factors that influence the bacterial community composition denoted by PERMANOVA R 2 values. Significance values indicate *p ≤ .05, **p ≤ .01, and ***p ≤ .001. (c) Seasonal bacterial community dissimilarities between urban and non‐urban areas. (d) Community dissimilarities between new and full moons per season and area. Dissimilarities were calculated as Bray Curtis distances. Differentially abundant ASVs clustered by family for (e) moon phases and (f) diel cycles. ASV counts were binned into families and transformed using the CLR transformation. Z‐scores were applied on families across seasons (rows), and hierarchical clustering was done using the “complete” method

DISCUSSION

Our year‐long in‐natura experiment comparing corals from urban and non‐urban areas indicate that normal diel, monthly, and annual biorhythms of corals are considerably disrupted by the urban conditions in the Gulf of Aqaba. The Gulf of Aqaba is an ideal study site because corals exhibit remarkable thermal tolerance, facilitating the determination of the effects of urbanization, without having to take into consideration the confounding impact of ongoing and increasing thermal stress by ocean warming (Savary et al., 2021; Voolstra et al., 2021). At the same time, the Gulf of Aqaba is a very small and constrained water body, which made the incorporation of more study sites (ideally, replicated urbanized and non‐urbanized areas) impossible. Rather, the urban and non‐urban study areas were only 6 km apart. It is noteworthy, that this spatial proximity renders a potential influence of genetic isolation or divergence that could arguably contribute to the observed differences unlikely. In this regard, the apparent chronic disruption of coral physiology is particularly troubling, because it directly aligns with coastal urbanization, that is, the presence of humans in the study area.

Despite the presence of large colonies of A. eurystoma of healthy appearance at similar densities in both areas, the comprehensive analyses conducted here suggest that major differences are apparent at physiological and metabolic level (Figure 6). The enrichment of genes involved in N fixation and photosynthesis in the non‐urban samples suggest that diazothrophs and phototrophs were more abundant at this area than in the urban one. The signature of “heavier” carbon and nitrogen isotope ratio in both coral tissue and symbiont samples from the urban area, together with the abundance of ABC transporter genes, supports our notion that there is an anthropogenic local disturbance of eutrophication (Figure S1b). This probably effects the metabolism and photosynthesis performances which impact the isotopic fractionation and photopigment synthesis (Ferrier‐Pagès & Leal, 2019; Muscatine, 1994; Muscatine & Cernichiari, 1969; Muscatine et al., 1989; Rädecker et al., 2015; Wall et al., 2019; Figure 1b,c; Figures S1b and S4d–f). The variations in δ13C values between the urban and non‐urban samples could indicate changes in the biomass composition (protein: lipid: carbohydrate ratios) of corals from each sampling area. In addition, the measured δ15N values could reflect on the temporal variability in nitrogen sources in both areas affecting the symbiont nitrogen demands (Levy et al., 2010; Wall et al., 2019).

FIGURE 6

Conceptual illustration of urbanization effects on coral reefs. Environmental conditions are represented by arrows pointed from the outside towards the coral. Biological outputs are represented by arrows pointed from the coral towards to outside

Artificial light at night (ALAN) is likely to play a significant role in disrupting the normal diel, lunar, and annual cycles of physiology and gene expression. Direct effects of ALAN on the circadian clock have been documented in a range of organisms, including corals (Ayalon et al., 2019; Davis et al., 2001; Rich & Longcore, 2013; Rosenberg, Doniger, Harii, et al., 2019; Rosenberg, Doniger, & Levy, 2019). The disrupted lunar rhythms of gene expression in the urban corals documented here are consistent with ALAN interfering with moonlight‐sensing systems as observed in other marine invertebrates (Kronfeld‐Schor et al., 2013). As seen by the monthly cycle of enriched processes in the non‐urban corals, most processes respond to the full moon where illumination is the strongest but, even at the full moon, the intensity of artificial light penetrates the water column at the urban area is greater (Tamir et al., 2017). The monthly enriched processes found in the non‐urban corals represent features of cellular organization, tightly connected to the actin and myosin filaments, enriched during the full moon (correlated to previously published work (Rosenberg et al., 2017). The reproductive processes, which are all aligned with the moon phase, are also abolished in the urban corals. Our results support the notion of moon light regulating expression patterns of clock genes in corals (Brady et al., 2016; Hoadley et al., 2011; Levy et al., 2007; Reitzel et al., 2013; Rosenberg, Doniger, Harii, et al., 2019) and the destructive effect of artificial light, over riding the moon light, causing delayed gametogenesis and loss of synchrony in gamete release as observed in many coral species across the globe (Ayalon et al., 2020; Jokiel et al., 1985; Kaniewska et al., 2015; Van Woesik et al., 2006). Our molecular data for these two cycles (diel and moon phase) emphasize the ability of light pollution to override the natural light/dark cycle and moon light, as well as masking biological processes associated with these cycles effecting corals in urban areas. Additionally, recent work (Lin et al., 2021) showed that dim light during the night suppressed spawning in the coral Dipsastraea speciosa. Importantly, this later study showed that the period of darkness between sunset and moonrise is essential to trigger synchronized mass spawning.

Normal seasonal rhythms of gene expression were also severely disrupted in the urban area, as the GO term enrichment analyses clearly indicate (Figure 4a). Seasonal GO term enrichment of sexual reproductive and cell cycle processes was only observed at the non‐urban area or in some cases, seasonally shifted between the two areas (e.g., RNA splicing). In the GO term analysis, the few processes that were only enriched in the urban area most likely reflect stress responses. Those included heme biosynthesis which is a biomarker for evaluating contamination in marine environments (Bogovski et al., 1998; Hongo et al., 2017), cell redox homeostasis which refers to the capacity of cells to continuously deal with challenges brought by different stressors, metabolic or environmental (Ursini et al., 2016), and ATP metabolic process, which is known indicator for the presence of contaminants (Kroll, 2009).

Diel and seasonal light regimes drive respiration and primary production, and in turn, are fundamental determinants of nutrient transformation and heterotrophic microbial diversity. We hypothesize that ALAN will affect these ecological processes and, therefore, diminish the microbial biorhythmicity as reefs are impacted by anthropogenic perturbation. Correspondingly, we found that the biorhythmicity of coral‐associated bacteria was strongly dependent not only on the season but also on the sampling area (i.e., reefs with different degrees of anthropogenic perturbation; Roder et al., 2015; Ziegler et al., 2019). Moreover, microbial community composition differed between moon phases in urban and non‐urban areas. The bacterial taxa that significantly responded to moon phases were common to both areas but fluctuated asynchronously between new and full moon at different seasons. In addition, lower bacterial beta diversity across seasons in the urban area indicated a less pronounced seasonality compared to the non‐urban area. Our results suggest that anthropogenic‐derived ALAN affects the biorhythmicity of coral microbiomes due to the potential demotion of light as an essential seasonal cue, as shown for microbial communities in sediment over time (Hölker et al., 2015). Coral‐associated bacteria are primarily heterotrophic but have been previously shown to respond to diel fluctuations driven by primary producer‐derived dissolved organic carbon (Kelly et al., 2019). Similarly, our results show that higher light availability states (i.e., full moon and daytime) were associated with a stronger microbial response, evidenced by a higher number of enriched bacterial taxa (ASVs). However, consistent with previous work (Baquiran et al., 2020; Silveira et al., 2017), we found no overall changes in the bacterial community composition (activity) between diel cycles. Instead, we identified taxa from bacterial families such as Rhodobacteraceae and Sphingomonadaceae oscillating between day and night, which putatively have food webs tightly coupled with Symbiodiniaceae and/or phytoplankton‐derived organic matter. Notably, our approach to assess community changes based on the active microbiome (cDNA‐based) may have helped resolve such differences. According to the predictions of bacterial functional signatures, photoautotrophs and diazotrophs were more abundant in urban samples than in non‐urban samples dominated by heterotrophs. This opens the possibility that bacterial communities in urban sites prefer autotrophic metabolism contrary to the heterotrophic metabolism preferred in urban sites. However, the accuracy of functional predictions for non‐human samples are limited, and the derived hypotheses require further investigation. Overall, our data suggest seasonal microbiome variation is affected by urbanization, indicating the impact of urbanization on the holobiont in general and not only the coral host.

Contemporary coral reef ecosystems are thought to have evolved in the last 45–50 million years (Close et al., 2020). During this period, reefs have experienced a wide range of natural disturbance regimes differing in magnitude, duration, and frequency, to which these complex ecosystems and their reef‐building corals have adapted and evolved (Buddemeier & Smith, 1999; Hatcher, 1997; Nyström et al., 2000). These natural disturbance regimes have led to the high species diversity, complex community structure and dynamics characteristic of pre‐industrial coral reefs (Pandolfi, 1999). Conversely, human‐induced disturbances often happen in a more persistent manner and occur at frequencies that prohibit adaptation, acclimatization, or recovery (Connell, 1997). In the longer run, even low levels of chronic stress can have severe impacts on coral reef ecosystems, causing decreased reproduction and growth rates, and compromising coral immunity (Richmond, 1993). In this presented work we have focused on interpreting the differences between urban and non‐urban sites mainly in regard to light pollution, nutrients, and eutrophication. This is since we have active monitoring data from the Israeli National Monitoring Program, regrading those stressors and knowledge of how they can impact the lifestyle of coral reefs from physiological and molecular aspects. However, we are not ignoring the fact that other sensory pollutants, such as chemical and hormonal pollution, which are not monitored, might add another layer of stress to the system. Given the critical importance of the coral holobiont to the fabric of coral reef ecosystems, the impact of increasing urbanization on coral biorhythms adds a further level of threat to an already compromised reef ecosystem unaccounted for in current projections of reef loss (Anthony et al., 2008; Fitt & Warner, 1995; Jokiel & Coles, 1990; Negri et al., 2005; Rj, 1997).

Finally, the coral reefs in the northern part of the Red Sea are considered as coral refuge from climate change and ocean acidification (Krueger et al., 2017). The increased economic interest, future development planned along the Red Sea coastlines, which are still not heavily populated, will eventually expose the Red Sea fringing reefs to human‐based disturbances in addition to the environmental threats (Fine et al., 2019; Loya, 2004); therefore, we hope our study can serve as warning to the potential sensory pollutants chronic disturbances can impair coral reefs.

CONFLICT OF INTEREST

Authors declare that they have no competing interests.

AUTHOR CONTRIBUTIONS

Conceptualization: YR, OL; Methodology: YR, OL, NSB, RY, AC, CRV; Theoretical analysis: ML, SA, AC, YR; Visualization: ML, SA, AC, YR; Funding acquisition: OL; Project administration: OL, YR, NSB; Writing—original draft: YR, OL, CRV, DJM; Writing—review & editing: YR, OL, NSB, ML, RY, AS, SA, GHP, AC, CRV, DJM.

Supplementary Material

Click here for additional data file.