Disentangling the complex microbial community of coral reefs using standardized Autonomous Reef Monitoring Structures (ARMS).

Autonomous Reef Monitoring Structures (ARMS) have been applied worldwide to describe eukaryotic cryptic reef fauna. Conversely, bacterial communities, which are critical components of coral reef ecosystem functioning, remain largely overlooked. Here we deployed 56 ARMS across the 2,000-km spread of the Red Sea to assay biodiversity, composition and inferred underlying functions of coral reef-associated bacterial communities via 16S rRNA gene sequencing. We found that bacterial community structure and diversity aligned with environmental differences. Indeed, sea surface temperature and macroalgae cover were key in explaining bacterial relative abundance. Importantly, taxonomic and functional alpha diversity decreased under more extreme environmental conditions (e.g., higher temperatures) in the southern Red Sea. This may imply a link between bacterial community diversity and functional capabilities, with implications for conservation management. Our study demonstrates the utility of ARMS to investigate the response of coral reef-associated bacterial communities to environmental change.

INTRODUCTION

Microbial communities of coral reef ecosystems are responsible for fundamental ecological processes such as photosynthesis and nitrogen fixation (Charpy, Casareto, Langlade, & Suzuki, 2012; Diaz & Rützler, 2001). Indeed, microbes associated with coral reefs putatively provide a substantial proportion of the reef nitrogen budget (Shashar, Cohen, Loya, & Sar, 1994). Furthermore, microbial biofilms have been implicated in facilitating larval settlement (Lau, Thiyagarajan, Cheung, & Qian, 2005; Wieczorek & Todd, 1998) as well as protecting macro‐organisms from biofouling (Armstrong, Yan, Boyd, Wright, & Burgess, 2001). This illustrates that the microbial community is critical for the resilience of coral reefs (Ainsworth, Thurber, & Gates, 2010; Garren & Azam, 2012). However, to date, studies of the microbial communities associated with these ecosystems have mainly focused on bacteria associated with corals or sponges (Lee et al., 2011; Sunagawa, Woodley, & Medina, 2010; Webster & Taylor, 2012; Ziegler et al., 2019, 2016; Ziegler, Seneca, Yum, Palumbi, & Voolstra, 2017), coral mucus (Frade, Roll, Bergauer, & Herndl, 2016; Ritchie, 2006), or other reef invertebrates (Bourne et al., 2013; Roik et al., 2016). Some studies have also concentrated on assessing the bacterial diversity in adjacent sediments (Rusch, Hannides, & Gaidos, 2009; Uthicke & McGuire, 2007) or the surrounding water column (Kelly et al., 2014; Somboonna et al., 2017). Yet, bacterial communities of coral reef frameworks—i.e., reef‐colonizing surfaces, cracks, and cavities of the limestone framework and associated biota (Gaidos, Rusch, & Ilardo, 2011; Schöttner et al., 2011)—remain understudied (Egan, Thomas, & Kjelleberg, 2008).

Assuming that these habitats represent the majority of niches within the complexity of coral reefs, understanding the composition and role of bacterial communities within the coral reef framework is crucial (Hernandez‐Agreda, Gates, & Ainsworth, 2017). This is especially true in an era where anthropogenic and natural pressures are having stronger and more frequent impacts on coral reefs (Hughes et al., 2017). Thus, understanding how environmental perturbations will affect the functioning of coral reefs and the underlying microbial processes is vital for the management of these iconic ecosystems (Garren & Azam, 2012). Environmental gradients along the latitudinal extent of the Red Sea provide a natural laboratory to investigate how bacterial communities associated with coral reefs respond to environmental changes (Berumen et al.,2019; Neave, Apprill, Aeby, Miyake, & Voolstra, 2019). However, progress has been partly hindered by the lack of a standardized and nondestructive tool to comprehensively investigate the complex microbial communities within coral reefs.

Over the last decade, coral reef scientists investigating changes in reef biodiversity have begun to focus on looking at small organisms that are missed during traditional biodiversity reef surveys (e.g., photo transect and visual censuses). Fostered by a need for a more comprehensive understanding of reef biodiversity in an era of accelerating global climate change and coral reef degradation, a team of scientists developed the Autonomous Reef Monitoring Structures (ARMS) (https://www.pifsc.noaa.gov/cred/survey_methods/arms/overview.php). ARMS are artificial structures using standardized materials and provide an ideal opportunity to study the small organisms within the reef habitat in a nondestructive way. Since standard protocols were developed, covering steps from deployment to sample collection and processing, the ARMS standardized framework ensures the samples obtained are comparable across various spatial scales. ARMS is now a validated tool to describe the biodiversity of the eukaryotic cryptic fauna associated with coral reefs (Al‐Rshaidat et al., 2016; Carvalho et al., 2019; Leray & Knowlton, 2015; Pearman, Anlauf, Irigoien, & Carvalho, 2016; Pearman et al., 2018; Plaisance, Caley, Brainard, & Knowlton, 2011; Ransome et al., 2017). More widely, ARMS have been applied to study the biodiversity in hard bottom marine habitats of the Atlantic, Indian and Pacific Oceans, as well as in the Black, Mediterranean and Red Seas (David et al., 2019). The global extent of deployments using a standardized tool will allow for a better understanding of marine biodiversity and will be particularly helpful to unveil how overlooked yet critical components of coral reefs respond to environmental changes from local to global scales. Despite the increasingly broad use of ARMS to assess reef biodiversity worldwide, they have not yet been utilized for the assessment of microbial diversity.

Here we sought to investigate the applicability of ARMS in elucidating compositional changes of bacterial communities within coral reef frameworks in different latitudinal regions. To do this, we conducted a survey using a set of 56 ARMS deployed in 19 coral reefs of the Red Sea spanning a latitudinal gradient of 16° and a distance of approximately 2,000 km. We hypothesized that the broad latitudinal gradient and associated differences in environmental variables might reveal drivers of bacterial community structure and function. We demonstrate that, in combination with high‐throughput sequencing techniques, ARMS allow the assessment of a broad range of bacterial communities (from surface‐associated bacteria to those in symbiotic relationships with sessile eukaryotic organisms) within coral reef ecosystems. Our approach argues for a more comprehensive evaluation of the complex microbiological community within coral reefs through a nondestructive and standardized methodology. Applying standardized approaches across variable spatial scales will contribute to a better understanding of the forces that govern the structure and function of microbial communities on coral reefs.

MATERIAL AND METHODS

Study area and environmental characterization

The study area encompassed the latitudinal extent of the Saudi Arabian Red Sea coastline (2,000 km). The Red Sea is a semi‐enclosed basin with only a single natural entrance at its southern end and has pronounced environmental gradients (Raitsos et al., 2011; Sofianos & Johns, 2002), being characterized by decreasing salinity and increasing temperature gradients from north to the south (Sofianos & Johns, 2002). In addition, inflow of Gulf of Aden water into the southern Red Sea brings high levels of nutrients and subsequently high chlorophyll a (Chla) concentrations. As the water progresses north, nutrient levels decline and the central and northern regions are classified as oligotrophic with periodic increases in nutrients due to seasonal physical oceanographic processes, particularly in the north (Acker, Leptoukh, Shen, Zhu, & Kempler, 2008; Kheireddine et al., 2017). For this study, based on a long‐term analysis (i.e., a 10‐year high‐resolution data set) of changes in sea surface temperature (SST) and productivity (Chla as a proxy) (Raitsos, Pradhan, Brewin, Stenchikov, & Hoteit, 2013), along with genetic connectivity boundaries (Nanninga, Saenz‐Agudelo, Manica, & Berumen, 2014), three hydrogeographical regions (North, Central and South) were considered.

During the period of deployment, monthly averages of Chla, SST, particulate organic carbon (POC) and photosynthetically active radiation (PAR) from the MODIS A satellite system were retrieved for each investigated reef at 4‐km resolution using the obpgcrawler (Tupper, 2018) package in R (R Core Team, 2018). Monthly records were averaged to allow for a comparative characterization of the environmental conditions where the ARMS were deployed. Benthic reef community surveys (20 m × 1 m transects at 8–10 m depth) were conducted at the time of ARMS deployment. In total, triplicate transects (separated by 5 m) were undertaken at each reef with photographs taken every 2 m as described by Ellis et al. (2017). Identification of the benthic groups (hard coral, soft coral, turf algae, macro algae and abiotic) along each transect utilized the software coral point count with excel extensions for 48 randomly distributed points over a 1‐m2 frame.

ARMS deployment–recovery and sample processing

Triplicate ARMS (spaced 2–5 m apart) were deployed at 19 reefs along the Saudi Arabian Red Sea coast (Table S1). Each ARMS unit was submerged for 2–3 years at 10 m depth on the hard framework of the reef (details of deployment, recovery and depth are provided in Table S1). Five reefs were investigated in the southern region, eight in the central region and six in the northern region (Table S1). Recovery and laboratory processing of the ARMS were undertaken following the protocols described by Leray and Knowlton (2015). Briefly, ARMS were retrieved from the reef sites and temporarily stored in 100‐μm filtered water from the reef during transportation to the laboratory. Samples were immediately processed upon arrival on shore.

ARMS units comprise a variety of niches (e.g., high light, low light and various flow regimes) and feature nine stacked PVC layers to mimic the complex coral reef framework. For this study, the settled community attached to the nine plates (top and bottom surfaces) was of interest. Plates were brushed to remove any mobile eukaryotic organisms and placed in 0.2‐μm filtered seawater before the top and bottom of each plate was scraped, combined and blended in a food processor. This resulted in a bulk homogenized sample containing all organisms which settled on the plate. A subsample of this bulk community (~10 g) was stored in dimethyl sulphoxide until DNA extraction was undertaken. This approach was performed to yield an integrative representation of the multiple reef niches and to allow for a nondestructive method to sample the bacterial communities of coral reef frameworks. Ten grams of the blended material was used for DNA extraction. The bacterial community, including surface‐associated bacteria, and those in association with sessile eukaryotic organisms, was then analysed by amplifying an ~430‐bp fragment of the v3 and v4 region of the 16S rRNA gene by polymerase chain reaction (PCR) (see Table S2 for details on primers and for PCR conditions see Klindworth et al., 2013). Approximately 5 ng of DNA was run in triplicate PCRs (25 cycles) using Invitrogen Taq Polymerase with negative PCR controls being included. The triplicate PCR products were pooled and 20 µl of the combined PCR products was cleaned and normalized using SequelPrep Normalization plates (ThermoFisher Scientific) resulting in a concentration of ~1 ng/μl. A second round of PCR amplification of eight cycles using KAPA 2 × HiFi Hot Start ReadyMix was undertaken to add Illumina Nextera tags following the manufacturer's recommendations. A second clean up and normalization step was done using SequelPrep Normalization plates (ThermoFisher Scientific). The Illumina MiSeq sequencing platform (v3 chemistry) at the King Abdullah University of Science and Technology (KAUST) Bioscience Core Laboratory (BCL) was used to determine sequence reads (2 × 300 bp). Raw reads were deposited at the NCBI Short Red Archive under the project accession PRJNA479721.

Amplicon sequence variant (ASV) inference and taxonomic assignments

Raw reads were processed, subsequent to primers being removed with cutadapt (Martin, 2011), using the dada2 package (Callahan et al., 2016) within R. Briefly, reads were trimmed based on quality, and filtered with a maxEE (maximum number of “expected errors”) of 2 for forward reads and 6 for reverse reads (reads not reaching this threshold were discarded). dada2 first constructs a parametric error matrix (based on the first 108 bp in the data set), then the sequences are dereplicated and sequence variants for the forward and reverse reads are inferred based on the derived error profiles from the samples. Singletons observed in the inference step are discarded. Subsequently, paired‐end reads were merged with a maximum mismatch of 1 bp and a required minimum overlap of 10 bp. Forward and reverse reads that did not merge were not included in further analysis. Chimeras were removed using the function remove BimeraDenovo. The resulting chimera‐checked, merged ASVs were used for taxonomic classification using the SILVA version 128 database (Pruesse et al., 2007) within the dada2 package, which is based on the rdp classifier (Wang, Garrity, Tiedje, & Cole, 2007) with a bootstrap of 50. The results were parsed into a table using the phyloseq package, and sequences assigned as eukaryotes, chloroplasts and mitochondria were removed. Samples were subsampled to an even depth of 25,100 sequences for comparison.

Functional characterization

Functional characterization was undertaken by performing a pathway analysis. To achieve this, the reference sequence of each ASV present in each reef was obtained and replicated to match the relative abundance of the ASV in the ASV table. To make metabolic inferences regarding the community, phylogenetic placement was undertaken using the resulting fasta files with the paprica (Bowman & Ducklow, 2015) pipeline. This places the query sequences into a reference phylogenetic tree produced from the 16S rRNA gene of sequenced genomes. Each node has a consensus genome derived from the genomes comprising the node. Paprica conducts a functional inference for each sequence by reconstructing the probable metabolism, based on phylogenetic comparisons to the published genomes deposited in MetaCyc (Caspi et al., 2007). MetaCyc is a highly curated database containing metabolic pathways that have been experimentally validated and reported in the scientific literature (Caspi et al., 2007). This analysis resulted in the relative abundance of pathways assigned per reef.

Shifts in bacterial communities and functions across the latitudinal gradient

Community composition was first calculated for the whole basin in phyloseq and was visualized using the sunburst package (Santesmasses, 2018). Subsequently, regional differences in the average number of ASVs and pathways were tested by analysis of variance (ANOVA) using a log transformation of the observed ASV numbers (two factors: Region [three levels] and Reefs [19 levels] nested within Region). For the pathway analysis, there was a single factor (Region [three levels]). Shared ASVs were calculated amongst regions in phyloseq and depicted using the venndiagram package (Chen & Boutros, 2011). Composition at the family level (agglomeration of ASVs belonging to the same family) was analysed for each reef separately (mean contribution of ASVs across the triplicate ARMS at each reef) with the positioning of the reefs provided—using the packages marmap (Pante & Simon‐Bouhet, 2013) and ggplot2 (Wickham, 2016). Differences in the contribution of bacterial families to the three hydrogeographical regions were similarly tested with ANOVA. Mantel correlations (Pearson) were used to test for differences in the distance matrices for the ASVs and functions.

Underlying processes influencing bacterial community assembly

To identify processes that contribute to the assembly of bacterial communities across the Red Sea gradient, we used ecological community modelling frameworks developed and updated by Stegen and colleagues (Stegen et al., 2013; Stegen, Lin, Fredrickson, & Konopka, 2015). This framework assesses the relative extent to which the assembly of ecological communities is governed by three processes as defined by Stegen et al. (2015): “(a) selection resulting from among‐taxa differences in performance; (b) dispersal resulting from organismal movement; and (c) ecological drift resulting from stochastic changes in population sizes.” This approach calculates the mean‐nearest‐taxon‐distance (βMNTD) metric to assess pairwise phylogenetic turnover between communities (Fine & Kembel, 2011). To undertake this, a phylogenetic tree was produced incorporating all ASVs with a mean abundance >0.001%. This metric quantifies the phylogenetic distance between each ASV in one community and its closest relative in a second community (Fine & Kembel, 2011). A null distribution of βMNTD values (repetitions = 999) was calculated using the R package picante (Kembel et al., 2010) under the assumption that ecological selection was not the primary cause of differences in pairs of communities. The beta‐nearest taxon index (βNTI) was evaluated by comparing the observed βMNTD values to the mean of a null distribution of βMNTD values and normalizing by its standard deviation (Stegen, Lin, Konopka, & Fredrickson, 2012). Selective pressures are indicated in this analysis by deviations away from the null βMNTD distribution. If environmental conditions are similar across samples, selective pressure will be consistent, which results in low levels of changes in the community. This is defined by having a βNTI value <−2. On the other hand, βNTI > 2 indicated variable selection. This is where changes in environmental conditions and among‐taxa fitness differences result in higher than expected pairwise differences in communities. Those pairwise comparisons that did not deviate from the null βMNTD distribution indicated that selective pressures were weak and were thus assessed for the influence of the stochastic measures, dispersal limitation and homogenizing dispersal. This was undertaken by calculating the Raup–Crick metric adapted to account for species relative abundances (RCBray). Again, a null distribution was calculated and the observed value was compared against this null distribution and standardized to values between −1 and +1 (Stegen et al., 2015). If selective pressures are low then large differences in the community are primarily due to dispersal limitation, as communities drift apart compositionally and is indicated by RCBray values >0.95. Finally, high rates of dispersal between a given pair of communities can be the primary cause of community similarity and this is called homogenizing dispersal and is denoted by RCBray values <−0.95. Values >−0.95 but <0.95 were interpreted as having no dominant assembly process.

Generalized linear models (GLMs) were used to investigate the effects of the environmental variables on the relative abundance patterns of prominent bacterial families. A stepwise selection (direction = both) was undertaken in the R package mass (Venables & Ripley, 2002) to find the variables that resulted in the model with the lowest Akaike information criterion (AIC). The predictors from the best model were then assessed for their relative importance with the package relaimpo (Grömping, 2006) using the lmg method (Lindeman, 1980). Pearson's correlations of these families against the environmental factors were undertaken. Principal coordinates analysis was run on both the ASV and the functional data sets (see above for how functions were assigned using paprica) based on Bray–Curtis dissimilarity matrices in R. Regional differences in the Bray–Curtis dissimilarity matrices, for both ASVs and functional pathways, were tested using the adonis function in vegan (Oksanen et al., 2007) (one factor: Region; three levels: North, Central, South). Bacterial indicator species and pathways for the regions were determined using the R package indicspecies (De Cáceres & Jansen, 2015). Indicator species were designated if they were significant for the region (p < .05). For further analysis, indicator ASVs had to have a relative abundance of 1% in at least two reefs. The resulting dendrogram and heatmap were plotted using the package cowplot (Wilke, 2016). Log values of the relative abundance of indicator ASVs and pathways were tested for correlation against environmental parameters using Pearson's correlations. Canonical correspondence analysis (CCA) was performed using the ordinate function in phyloseq for those ASVs belonging to the family Rhodobacteraceae and plotted in ggplot2.

RESULTS

Compositional changes of bacterial communities associated with coral reef frameworks highlight regional differences

From sequencing of the 16S rRNA gene associated with ARMS structures, we observed a total of 24,385 ASVs across 19 reefs of the Red Sea spanning 16° of latitude. Proteobacteria accounted for the highest proportion of ASVs (59%) and sequences (70%). The Proteobacteria predominantly comprised the classes Alphaproteobacteria (38% of all sequences) and Gammaproteobacteria (21% of all sequences) with the family Rhodobacteraeceae (class Alphaproteobacteria, order Rhodobacterales; 17% of all sequences) being especially prevalent (Figure S1). Bacteroidetes were the second largest phylum in terms of ASV numbers (12% of ASVs), but they only contributed 4% of the total sequences. Although accounting for a lower number of ASVs, in terms of number of sequences, Chloroflexi (5%), Acidobacteria (4%) and Cyanobacteria (4%) showed a comparable contribution.

Significant regional differences were observed in the community structure at both the ASV and the family levels (adonis F = 2.53; p < .001 for ASV and F = 2.64; p < .001 for family). In general, Rhodobacteraceae was the most predominant family (Figure 1), ranging from 7% to 34% of sequences per reef. Nitrosomonadaceace (F = 15.653; p = .001) showed a statistically higher relative abundance in the central region, compared with the south. Rhodobacteraceae (F = 5.23; p = .37) and Rhodospirillaceae (F = 4.55; p = .05), which contributed substantially to the Red Sea community, showed a significant interaction between site and region.

Figure 1

Bacterial community composition (family level) across the latitudinal spread of the Red Sea (19 reefs, 2,000 km) denoting three environmentally distinct regions (North, Central, South). The most abundant families (top 10) in terms of sequence number for each region are shown. Families are represented by different colours and the size of the bubble indicates the family's proportional relative abundance in the reef. Each line of bubbles represents a reef and the specific reef is denoted by the arrow. The stacked bar chart indicates the contribution of the families depicted on a regional scale. Short codes (where appropriate) for reef names are: AFL, Al Fahal; AMF, Abu Madafi; ASA, Abu Shoosha; ASF, Abu Shootaf; CG, Coastguard; NAL, North Al Lith; SAL, South Al Lith; WS, Whale Shark. Family I refers to a cyanobacterial family

The vast majority of the ASVs observed (~85%) showed restricted spatial ranges, being found in two or fewer reefs (Figure S2). However, these ASVs only accounted for 17.5% of the sequences. Only 3.7% of the ASVs were observed in all three regions (i.e., North, Central, South). However, they accounted for 50.9% of all sequences. Significant differences in the number of ASVs were found amongst reefs (nested within region) (ANOVA: F = 5.626; p < .001) and amongst regions (ANOVA: F = 72.830; p < .001). The highest number of ASVs was observed in the northern region (12,901 ASVs) compared with the 4,892 ASVs found in the south. A similar pattern was detected for the average number of ASVs per reef—i.e., a significantly higher number of ASVs (mean = 3,391) was observed in the northern region compared to the central (mean = 1,941 ASVs; p = .03) and southern regions (mean = 1,195 ASVs; p = .002) (Figure S3b). The northern region showed the highest number of ASVs restricted only to a single region. The central region showed the highest number of shared ASVs, especially with the northern region (Figure S3a).

Differences in bacterial community composition aligned weakly with differences in underlying functional capability along the latitudinal gradient (Mantel test R = .326; p < .003). Pathway profiles (the composition of pathways for each reef) were statistically different (F = 3.288; p = .001) amongst regions. In total, 697 pathways were detected, 655 of which were shared amongst all regions, with the vast majority being found in all reefs (Figure S2 and Figure S3c). There was a significant difference (ANOVA; F = 6.04; p = .011) in pathway numbers with a higher diversity in the northern region compared to the south (p = .008) (Figure S3d).

Distribution of bacterial community members determined by selective processes

Applying the ecological modelling framework of Stegen et al. (2015), we identified variable selection (i.e., the shifting of selective pressure resulting from differences in environmental conditions) to be the main factor determining the compositional change in the bacterial community (accounting for 67% of pairwise comparisons between reefs) (Figure 2; βNTI results in Table S3). However, homogeneous selection, in which environmental pressures are constant, was only observed in 6% of pairwise comparisons, mainly being present in comparisons between northern reefs. Homogenizing dispersal (i.e., areas where high dispersal rates lead to low levels of spatial compositional changes) and ecological drift, due to dispersal limitation, contributed 17% and 3%, respectively, to compositional changes (Figure 2).

Figure 2

Assembly of the bacterial community showing the proportion of pairwise comparisons attributed to the different contributing processes. DL, dispersal limitation; HD, homogenizing dispersal; HS, homogeneous selection; NDP, no dominant process; VS, variable Selection

Because spatially divergent selection pressures (from changing environmental conditions) were shown to be relevant in determining the distribution of community members, we investigated regional‐scale differences in environmental and biological components. As expected, regional differences were observed for SST (χ 2 = 15.819, p < .001) and Chla (χ 2 = 8.66, p = .013), with highest values being present in the south. However, the percentage cover of the main benthic components (hard coral, soft coral, macroalgae and turf algae) did not show significant regional differences.

Investigation of how much of the variance in relative abundance of the dominant bacterial families was explained by environmental variables indicated that the best generalized models for each family had R 2 values ranging between 0.29 and 1 (Table S4). The main explanatory variable for the changes in relative abundance of the families Rhodobacteraceae, Phyllobacteriaceae, Flammeovirgaceae, Sva0996 and TK85 was SST (Figure 3). The latter three families were negatively correlated with SST while Phyllobacteriaceae had a positive relationship. The percentage cover of macroalgae was the main explanatory variable for changes in the relative abundance of Nitrosomonadaceae, SAR116 and Pseudoalteromonadace. Nitrosomonadaceae and Sva0996 had a significant negative relationship with macroalgae cover while Pseudoalteromonadaceae was positively correlated with macroalgae cover. Proportional cover of turf algae and soft corals were the dominant predictor variables for two families each (turf: Rhodospirillaceae and Planctomyceaceae; soft coral: Flavobacteriaceae and JTB255). Although the percentage of the reef assigned as abiotic (i.e., sand, rock, rubble) was not the main explanatory variable for the bacterial families studied, it did contribute substantially in explaining the variance of the bacterial reef community (Figure 3). The family Rhodobacteraceae (the most abundant family) revealed differential trends in the distribution of the most abundant ASVs. Two ASVs of the genus Roseovarius showed high relative abundance in the northern region, while the four most abundant ASVs of Rubribacterium were associated with the central region (Figure S4). Notably, Rubribacterium correlated negatively with SST, Chla, and the percentage cover of macroalgae, but positively with soft coral and turf algae percentage cover (Table S5). Similarly, Roseovarius was negatively correlated with Chla while Ruegeria was positively correlated with SST, Chla and macroalgae (Table S5).

Figure 3

Environmental factors contributing to bacterial community differences across reefs of the Red Sea. The top 10 families from each region as depicted in Figure 1 are shown. Stepwise selection of variables was undertaken for each family and the subset of variables present in the best GLM (lowest AIC result) was selected. For each family the relative importance (sum of variables will equal 1) of the subset of variables was calculated in relaimpo using the lmg method. Environmental variables: SST = sea surface temperature; Chla = chlorophyll a; HC = % cover of hard coral; SC = % cover of soft coral; Turf = % cover of turf algae; MA = % cover of macroalgae. Abiotic = % cover of abiotic elements. The bottom panel indicates the relative proportion of average environmental parameters in each region

The indval analysis revealed a total of 1,007 significant ASVs as putative indicator taxa for the regions, with 921 being indicators of the northern region. The vast majority of these indicator taxa had a low relative abundance and restricted spatial distribution. A total of 14 ASVs (Table S6) occurred in at least two reefs with a mean relative abundance greater than 1%. Abundant indicator ASVs represented a broad range of seven families with three ASVs being characteristic of the southern region, three of the northern region and eight of the central region (Figure S5a; Table S6). The response of those characteristic ASVs to environmental variables revealed that the two indicator ASVs of the family Nitrosomonadaceae were highly correlated with soft corals (Table S6). Indicator ASVs of the families JTB23_unclassified and FamilyI were negatively correlated with SST. Two SAR116 indicator ASVs and the Phyllobacteraceae ASV showed positive relationships with SST. Significant indicator pathways were found for reefs in the northern (five pathways; Figure S5b) and southern (one pathway) regions. The indicator pathways (five) of the northern region were negatively correlated with SST (Table S7).

DISCUSSION

Here we present a large‐scale survey using a standardized approach to the study of coral reef bacterial communities along 16° of latitude across strong environmental gradients. Bulk analysis of the scraped component of the ARMS enabled us to resolve for the first time the bacterial communities colonizing hard substrates as well those associated with sessile eukaryotes inhabiting the coral reef framework.

The main contributors to the reef framework bacterial communities were the Proteobacteria, with the classes Alphaproteobacteria and Gammaproteobacteria being especially prevalent. The predominance of these bacterial classes agrees with other studies investigating bacterial diversity within sediments of coral reefs (Rusch et al., 2009; Uthicke & McGuire, 2007) or associated with holobionts of coral or sponges in reef environments (Olson & Kellogg, 2010; Rusch et al., 2009). However, variations in the relative abundance of these classes across the reefs were observed. Alphaproteobacteria have been shown to be major components of coral reef invertebrates lacking photosymbionts, such as particular sponges and ascidians (Bourne et al., 2013; Erwin, Pineda, Webster, Turon, & Lopez‐Legentil, 2014; Evans, Erwin, Shenkar, & López‐Legentil, 2017; Webster & Hill, 2001), while Rhodobacteraceae have been shown to be important components of marine biofilms (Dang & Lovell, 2000; Jones, Cottrell, Kirchman, & Dexter, 2007; Witt, Wild, Anthony, Diaz‐Pulido, & Uthicke, 2011). In contrast, Gammaproteobacteria, the second most abundant group in the current study, have been mainly associated with organisms harbouring photosynthetic symbionts (Bourne et al., 2013; Bourne & Munn, 2005; Cárdenas, Rodriguez‐R, Pizarro, Cadavid, & Arévalo‐Ferro, 2012; Hernandez‐Agreda et al., 2017). Other classes typical of microbial communities associated with coral reefs (e.g., Cyanobacteria, Bacteroidetes, Chloroflexi, Acidiomicrobia, Actinobacteria, Planctomycetes) contributed to a diverse consortium of bacteria associated with the benthic reef habitat (Gaidos et al., 2011; Kelly et al., 2014).

Our data show that coral reefs along the investigated latitudinal gradient harbour distinct bacterial communities in terms of structure and composition. For example, communities associated with the northern region showed significantly higher diversity compared with the southern region, including a large number of ASVs that were indicative of the northern region. However, the vast majority of these ASVs were of low relative abundance and frequency. Although it is acknowledged that the inference of functional pathways based on 16S rRNA gene data is biased towards the availability of sequenced genomes (Bowman & Ducklow, 2015), we have demonstrated that functional capabilities of reef‐associated bacteria are regionally structured across the latitudinal gradient. However, a substantial portion of pathways were conserved across regions. Interestingly, contrasting patterns were observed between the frequency of occurrence of ASVs and pathways, with ASVs generally being restricted to a single reef while pathways were found across the Red Sea. This suggests that despite changes in the community, functional capabilities are present throughout. This is further corroborated by the significant positive correlation between the community structure and relative abundance of pathways. These results further indicate the coexistence of organisms with the Red Sea community that share specific functions, but differ in other ecological requirements, as has been previously suggested for the planktonic community (Galand, Pereira, Hochart, Auguet, & Debroas, 2018). As bacteria are vital in the functioning of coral reef ecosystems (Ainsworth et al., 2010), understanding the degree of functional redundancy is critical to better understand how species loss will ultimately affect coral reefs. Here, we found a correlation between ASV diversity and functional diversity. Indeed, the southern Red Sea exhibited the lowest number of bacterial taxa and the lowest level of functional diversity. The relationship between taxonomic and functional diversity has been suggested to correlate if there is an increase in the coverage of niche space as species richness increases (Díaz & Cabido, 2001). Although this relationship suggests that conservation management decisions targeting high‐diversity areas would preserve the highest amount of diversity, care needs to be taken to ensure that low‐diversity regions do not contain unique functions.

The Red Sea presents a unique combination of environmental variables that have previously been shown to shape the planktonic community (Kürten et al., 2016, 2014; Ngugi, Antunes, Brune, & Stingl, 2012; Pearman et al., 2017; Pearman, Kurten, Sarma, Jones, & Carvalho, 2016) and the metazoan component of the reef cryptobiome (Carvalho et al., 2019). However, knowledge of the factors structuring the composition and function of bacterial communities within coral reefs remains limited (Roik et al., 2016; Neave et al., 2019). Interestingly, despite the 2,000 km latitudinal range investigated, drift as a result of dispersal limitation was not a substantial contributor to the structure of the bacterial communities. Furthermore, within the Red Sea selective processes appeared dominant. Variable selection contributed substantially to the structuring of bacterial communities; this trend has been previously shown for prokaryotes in a variety of environments (Dumbrell, Nelson, Helgason, Dytham, & Fitter, 2010; Graham et al., 2017; Tripathi et al., 2018). Homogeneous selection, although not a substantial contributor to determining microbial communities across the Red Sea, did contribute in the northern Red Sea where it was the predominant structuring process. This finding may suggest that the contribution of the two types of selection (homogeneous and variable) is partially scale‐dependent. Notably, from our results, homogeneous selection may be more dominant within a region (North) whereas variable selection is more prevalent between regions, where the differences in environmental gradients are strongest.

Importantly, SST was identified as a key driver for various bacterial families in the reef community, in line with its critical role for the microbial communities comprising the coral microbiome (Grottoli et al., 2018), and as reported in global microbial plankton studies (Sunagawa et al., 2015). Various bacterial families showed correlations with temperature (e.g., Planctomycetaceae, SAR116 and Phyllobacteriaceae having a positive relationship, and Flammeovirgaceae having a negative correlation). Understanding how critical components of coral reef‐associated bacteria respond to temperature gradients will be vital as SST is predicted to rise with climate change. Chla concentrations were highest in the southern region, indicative of this region being more productive due to nutrient input from the Gulf Aden (Kheireddine et al., 2017; Qurban, Wafar, Jyothibabu, & Manikandan, 2017). The increased nutrient enrichment occurring in the south could explain the reduced bacterial diversity found in this region (Qurban et al., 2017). A lower diversity with increasing productivity has been previously observed in other regions and habitats (Milici et al., 2016; Wemheuer et al., 2014). Further studies evaluating the importance and role of environmental factors across environmental gradients, as well as investigating seasonal differences, will provide a better understanding of the drivers of microbial reef community assembly and function.

The relative proportion of the macrobenthic community within a reef (i.e., hard corals, soft corals, turf and macroalgae) has been postulated to influence microbial community composition (Kelly et al., 2014). Indeed, for a variety of bacterial families we found that macroalgae cover had a high relative importance in explaining the variance in relative abundance. The shift from a dominance of calcifying to noncalcifying benthic groups (turf algae, soft corals, fleshy macroalgae) has repercussions for the resilience of coral reefs (Elmhirst, Connolly, & Hughes, 2009) with distinct bacterial taxa being associated with coral and macroalgae (Barott et al., 2011). For example, Pseudoalteromonadaceae responded positively to the percentage of macroalgal cover, while other abundant families such as Nitrosomonadaceae and Sva0996, as well as the genus Rubribacterium showed the opposite trend. These differential responses of microbial groups to one or several potential indicators of reef degradation (in this example, the percentage of macroalgae cover) may lead to changes in the functional profile of the bacterial community. Further analysis of the composition of the eukaryotic composition settling on the plates would enable a more fine‐scale determination of the effect of eukaryotic organisms on the formation of prokaryotic communities on reef hard substrates.

The impact of climate change on marine ecosystems makes studying the effects of ocean warming, ocean acidification, changes in oxygen and nutrient saturation, etc., a priority. The latitudinal gradient of the Red Sea provides a natural laboratory to investigate how bacterial communities associated with coral reefs respond to changes in environmental variables and especially SST. The observed SST gradient allowed us to reveal differences in the composition and function of the bacterial community from the cooler northern regions, typical of other current global ocean conditions, to the south where temperatures 2°C above those observed in the north may be seen as representative of future conditions under global warming predictions (Hughes et al., 2017).

ARMS were originally designed to target the invertebrate fauna within a coral reef. The standardized design of the structure, incorporating a variety of niches (e.g., high light, low light and various flow regimes), provides an ideal methodology for investigating large‐scale patterns in bacterial diversity. The use of artificial substrates is likely to have biases in the resulting bacterial community as selective pressures from the macro‐organisms will be missing (Sweet, Croquer, & Bythell, 2011). However, the use of standardized materials will enable a comparative approach to be undertaken which is likely to be representative of changes occurring in the reef community. Furthermore, while ARMS do provide a variety of niches for bacteria to colonize, it must be acknowledged that other environments within the reef habitat or system, such as reef sediment, microbiomes of the mature reef framework (such as corals and sponges) and the water column, cannot be assessed using ARMS methodologies and because they target different habitats will give complementary information on the processes occurring within a reef.

Currently, conservation measures are mainly based on the responses of macro‐organisms such as changes in coral coverage and fish abundance and/or biomass (Hill & Wilkinson, 2004), but these are likely to respond to the presence and absence of prokaryotic organisms, which are vital players to reef functioning (Glasl, Bourne, Frade, & Webster, 2018). Our study highlights that a standardized methodology, which can be applied on a global scale, can provide insight into coral reef microbial communities (free‐living, benthic and organism‐associated) and determine differences along environmental gradients. Thus, we advocate that ARMS can provide standardized comparative data on a global scale and thus should be included in integrative conservation management decisions that incorporate information from prokaryotes to macro‐organisms.

CONFLICT OF INTEREST

The authors declare no competing interests.

AUTHOR CONTRIBUTIONS

S.C. designed the experiment. J.K.P., E.A., R.V., H.A. and S.C. undertook sampling and processing of the samples. J.K.P. and E.A. worked on the molecular samples and the bioinformatics. J.K.P., E.A., C.R.V. and S.C. wrote the paper while all authors commented on various drafts.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.