The microbial biosphere of the coral Acropora cervicornis in Northeastern Puerto Rico

Coral samples

Coral sampling was conducted within the La Cordillera Natural Reserve (LCNR), specifically at the northern shores of Isla Palomino, which is located 3.2 km from the northeastern coast of Puerto Rico (18°21′10.8″N 65°34′24.4″W, Fig. 1). The sampling site is a low topographic relief, fringing reef with relatively clear waters year-around and moderate to high wave energy depending on the water depth. Due to the reef’s northeastern orientation, the site is exposed to easterly trade winds. In addition, given the absence of perennial freshwater streams in Palomino Island, concomitant with the distance from the Fajardo watershed, (the nearest freshwater system is ∼6 km west of the site), coastal waters are mostly free of terrestrial sediments and thus have a horizontal water visibility well over 10 m. The corals assemblage at the collection site is dominated by gorgonian corals such as Gorgonia ventalina and Pseudopterogorgia aerosa and small colonies of the scleractinian corals Orbicella annularis, Acropora palmata, Porites astreoides and spread Acropora cervicornis clusters composed of several individuals. For a further description of the study area please see Hernandez-Delgado et al. (2006), Mercado-Molina et al. (2015) and Ruiz-Diaz et al. (2016).

Within this site, six A. cervicornis colonies were collected in August of 2015. Three of them were collected at a depth of 1.5 m (hereafter “shallow” samples) and three from 11 m depth (hereafter “deep” samples). Shallow samples were collected from the reef crest and were at least 5–6 m apart from each other. Corals of the reef crest are nearly continuously exposed to relatively high water motion, and reef crest temperature averages ∼29 °C. Solar radiation was ∼11,203.55 Lux (SI unit of illuminance = 1 lumen/m2) at the time of sampling. Deep samples were collected from at least 8–10 m apart within the back-reef zone. This zone is characterized by lower water motion and wave action, and average temperature is ∼28 °C. Solar radiation was ∼3,429.36 Lux at the time of collection. Temperature and light intensities were estimated by placing one Hobo Pendant temperature/light data logger 64k-UA-002-64 (Onset Company) per depth at the study site at the time of collection. Although we do not have multiple log readings per site, a previous study conducted by the current authors in the same collection sites found that the temperature and light were significantly different between the two depths (Ruiz-Diaz et al., 2016). Collected fragments were 6 cm in length from the tip of the colonies. Each fragment was individually placed in a sterile vial while underwater and the vials containing the coral fragments were put in dry ice once they reached the surface. Samples were stored at −80 °C until processing. Sampling was approved by the Department of Natural and Environmental Resources of Puerto Rico permit number DRNA: 2016-IC-175 issued to Carlos Toledo-Hernández.

DNA extraction

Coral samples were prepared for DNA extraction by scraping the surfaces of the tips (top 1.5 cm) of the branches (see Fig. 1) with sterile scalpel blades. This resulted in ∼400 mg of mucus, tissue and some skeletal material for each samples. The scraped material was subsequently processed using the PowerSoil DNA isolation kit (MO BIO, Carlsbad, CA, USA) following the manufacturer’s specifications. To increase DNA yield, a homogenization step was performed by mixing the lysate with beads using a PowerLyzer™ 24 Bench Top Bead-Based Homogenizer (MO BIO, Carlsbad, CA, USA) for 2 min at 2,000 rpm. Additionally, a second DNA extraction was done using the resultant pellet formed after the first homogenization step of the DNA extraction. Both DNA extractions (from the tissue and the tissue fragments in the bead tube) were pooled. Genomic DNA quality control was assessed using agarose gel electrophoresis with DNA standards of known molecular weight and concentration yielding 20–30 ng gDNA per sample. No further DNA purification steps were performed as PowerSoil DNA isolation kits are known to be effective at removing PCR inhibitors (Santos et al., 2012).

16S rRNA gene PCR and amplicon deposition

The V4 hypervariable region of the 16S ribosomal RNA was amplified by PCR using the universal bacterial and archaeal primers: 515F (5′ GTGCCAGCMGCCGCGGTAA 3′) and 806R (5′ GGACTACHVGGGTWTCTAAT 3′) as used with the Earth Microbiome Project (Caporaso et al., 2012). Amplification conditions were 1 cycle of 94  °C for 3 min, 35 cycles of 94 °C for 45 s, 50 °C for 60 s, 72 °C for 90 s and a final extension of 72 °C for 10 min. The six amplicons of ∼300 bp were barcoded to allow for sample multiplexing and paired-end sequenced in the Illumina MiSeq platform at the Sequencing and Genotyping Facility of the University of Puerto Rico. The resulting demultiplexed raw sequences per sample, as well as the 803 16S rRNA gene sequence representatives of the operational taxonomic units (OTUs) per sample were deposited in the NCBI BioProject ID PRJNA379103 with SRA accession SRP102061.

Community profiling and bioinformatics

Demultiplexed reads underwent quality control using QIIME (Kuczynski et al., 2012), selecting those reads with Phred scores > 20 (99% confidence) and lengths > 200 bp which were searched for chimeras with the usearch61 hierarchical clustering method (Edgar, 2010). Sequences were binned into OTUs in QIIME using de novo OTU assignment methods (thus the “de novo” ids for OTUs), with a 97% sequence similarity with Greengenes core representative sequences Gg_13_8_99.taxonomy (McDonald et al., 2011). The algorithm chooses an OTU “centroid sequence” to be the representative sequence for each OTU. Sequence alignment was done using the Python nearest alignment space termination tool (PyNAST), and taxonomy assignment was done with the uclust consensus taxonomy assigner using QIIME’s default settings (Kuczynski et al., 2012). Chloroplast and mitochondria OTUs were removed from downstream analyses using the script filter_taxa_from_otu_table available in QIIME (Kuczynski et al., 2012). Additionally, stringent OTU filtering included the removal of singletons and OTUs with less than two sequences per sample to eliminate overestimation caused by sequencing artifacts. Data analyses including diversity estimates and taxonomic composition of the coral samples was done using QIIME 1.9.1 (Kuczynski et al., 2012) after data was subsampled with a rarefaction level of ∼28,000 sequences per sample, to mitigate bias in the analyses due to differences in sampling depth. Alpha diversity was calculated using the phylogenetic diversity metric of Faith (PD) to assess community diversity of the samples. Abundance was not taken into account, but rather the branch lengths of the phylogenies connecting all species to each community (Faith, 1992). Significant differences in the rarefaction curves were calculated with a non-parametric two-sample t-test using 999 Monte Carlo permutations using the QIIME script compare_alpha_diversity.py. Beta and alpha diversity analyses, as well as core microbiome analyses, were done through QIIME (Kuczynski et al., 2012). Taxa summaries were built by modifying the QIIME L2 and L6 taxonomy tables with the R package reshape2 (Wickham, 2007).

To explain differences among microbial communities inhabiting shallow and deep corals, we used principal coordinates analysis (PCoA) on UniFrac distances, a beta-diversity measure that uses phylogenetic information to compare samples (Lozupone, Hamady & Knight, 2006). Statistical analyses on beta diversity were made using ANOSIM, a non-parametric statistical test that compares ranked beta diversity distances between different group depths found in the mapping file and calculates a p-value based on the unweighted Unifrac table used to generate the 3D PCoA plots (Kuczynski et al., 2012). The test was done using the script compare_categories.py in QIIME (Kuczynski et al., 2012) with the Unifrac distance matrix as the input file. Unweighted UniFrac PCoA biplots were visualized in the EMPeror Visualization Program (Vazquez-Baeza et al., 2013). Additionally, beta diversity analyses using non-metric multidimensional scaling (nMDS) ordinations of Bray Curtis dissimilarity was done using distance metrics computed from the rarefied OTU table and the metadata. We used nMDS ordination, achieved by the metaMDS wrapper function from the vegan package in R (Oksanen et al., 2008). The ordination was applied such that the data was scaled down to two dimensions.

The same analysis was done to understand the depth of diversity of the dominant Rickettsiales OTUs. A species table was prepared, and OTUs were filtered to retain only Rickettsiales OTUs; then, the resulting OTUs were compared between the two sampling depths. To analyze species composition similarities across sites a Principal Coordinate Analysis (PCoA) ordination was used. The PCoA was constructed using a UniFrac distance matrix and visualized through QIIME. Additionally, we used a log-likelihood ratio to test which OTUs changed significantly in relative abundance between the two depths. This test compares the ratios of the OTU frequencies in the sample groups to an “extrinsic hypothesis” that assumes that all sample groups have equal OTU frequencies, thus revealing which taxa have significantly different OTU frequencies.

An alpha rarefaction plot was built for richness estimations using the Chao 1 values. These values represent the estimated true species richness of a sample and are calculated through the workflow script for performing alpha rarefaction in QIIME that in turn implements the Chao 1 abundance-based estimator (Chao, 1987). An alpha rarefaction plot was then built for the Chao 1 richness estimations between microbial communities in deep and shallow coral samples. Additionally, we used a non-parametric two-sample t-test statistical with Monte Carlo permutations to compare the richness curves between shallow and deep group samples.

OTUs that changed significantly between the two sample categories were found using the script group_significance.py through the implementation of a G-test (log-likelihood ratio test) that compares the frequency of OTUs across all samples and finds those that are significantly different between both groups. Those significantly different OTUs (p < 0.05) were grouped into an OTU table that underwent DESeq2 negative binomial Wald normalization for visualization purposes as the numbers of individuals per sample greatly varied. This normalization step was implemented in QIIME using the script normalize_table.py. The normalization in QIIME uses a variance stabilization transformation function (VST) that is applied to the count data, as described in https://www.rdocumentation.org/packages/DESeq/versions/1.24.0/topics/varianceStabilizingTransformation. The heatmap showing the taxa that significantly differed in abundance between depths (p < 0.05) was built using the heatmap.3 function in R (Zhao et al., 2014). Significant taxa were further highlighted using boxplots made with the vegan package in R (Oksanen et al., 2008). A detailed repository and tutorial with all the necessary intermediary files and scripts to serve as a guide to generate the plots and perform the analyses used in this publication is available on GitHub: https://github.com/meglab2017/The-microbial-biosphere-of-the-coral-Acropora-cervicornis-in-Northeastern-Puerto-Rico.

Alpha and beta diversity estimates of corals at different depths

Coral individuals collected at 1.5 m (shallow samples 1, 2 and 3 hereafter) received more solar radiation compared to those at 11 m (deep samples 1, 2 and 3 hereafter). The total number of raw reads was 715,825 for the six samples. A total of 173,137 sequences were used for analysis, in which the average ± standard deviation of the number of reads was 29,135 ± 177.4 for the shallow samples and 28,578 ± 368 for deep samples. These sequences correspond to a subsampling to a rarefaction level of 28,000 sequences, without a replacement, thus guaranteeing an unbiased analysis (Table 1). The binning of the 173,137 high quality reads resulted in 803 operational taxonomic units (OTUs) (Table 1, Table S1).

nMDS ordination based on the relative dissimilarities of the samples (Bray Curtis) shows that shallow samples have higher dispersion and are separated from the deep samples by axis 1 (Fig. 2A). PCoA revealed that the first principal component, sample origin, represented the highest variance (PC1 = 46.16%, Fig. 2B) indicating that bacterial communities from coral samples partition according to sampling depth. Nonetheless, the non-parametric two-sample t-test ANOSIM revealed non-significant differences (Rstat = 0.59; p = 0.221) thus demonstrating that the similarity between groups is not significantly greater than the similarity within each of the sample groups. In addition, alpha diversity measures revealed no significant differences in phylogenetic diversity (t-stat = 0.548; p = 0.891, Fig. 2C). When visualizing the most abundant bacterial taxa associated with each sample, biplot indicated that Rickettsiales, Serratia marcescens and Lactococcus were most likely to be found in shallow coral samples, while Prevotella, Lactobacillus and Campylobacterales were more likely associated with the deep samples (Fig. 2B). Taken together, these results indicate that the species diversity did not vary significantly within each sampling depth and that there were only subtle differences in diversity between shallow and deep samples.

Taxonomic profiles of the A. cervicornis microbiome at different depths

We proceeded to explore the taxonomic profiles of the A. cervicornis-associated microbiome at different depths. At the phyla-level, our community profile analysis showed a total of eight phyla, with a dominance of Proteobacteria (specifically from the order Rickettsiales) at both sampling depths (95%; Fig. 3A, left panel). When these dominant Rickettsiales OTUs were filtered out of our analysis, other taxa from the Proteobacteria phylum dominated at both sampling depths (∼55%) followed by the phyla Firmicutes (∼35%), Bacteroidetes (∼5%) and Actinobacteria (∼2%) (Fig. 3A, right panel). Despite the dominant taxa being shared by corals from both sampling depths, some low dominance phyla appeared exclusively at only one of the sampling depths. Two phyla, Planctomycetes and Nitrospira, were only observed in the shallow samples, whereas Gemmatimonadetes and Verrucomicrobia were uniquely observed in deep samples (Fig. 3A, right panel).

The community profile analysis at the genus level revealed from the rarefied analyses revealed a total of 38 different genus-level OTUs, 19 of them were shared between shallow and deep samples, while 8 and 11 were exclusively isolated from shallow and deep samples respectively (Fig. 3B). Similar to the phylum-level analysis, most of these genus-level OTUs belonged to Proteobacteria and Rickettisales (95%). A total of 725 Rickettsiales-like OTUs were dominant across all samples independently of sampling depth. An analysis of the 78 non-Rickettsiales OTUs (<5%) revealed that the taxa most likely to differ in abundance between depths were low abundance Proteobacteria (Fig. 3B, left panel). In fact, nearly 60% of these were Haemophilus, Acinetobacter, Rhodobacteraceae, Serratia, Pseudomonas, Campylobacterales or other unclassified Proteobacteria (Fig. 3B, left panel). The other dominant phylum in the non-Rickettsiales group was Firmicutes (∼18%), with the most dominant genus being, Lactococcus and Peptostreptococcaceae mostly in the shallow samples while Lactobacillus, Anaerococcus and Bacillus were most dominant in the deep samples.

To gauge the diversity of the Rickettsiales-like OTUs, we performed a community profile analysis comparing Rickettsiales-like OTUs within shallow and deep samples using a rarefaction of 23,472 sequences (Fig. 4A). Out of the 725 Rickettsiales-like OTUs, a single OTU was most dominant across all samples (Fig. 4A, denovo_0). Furthermore, Rickettsiales richness was not significantly different between shallow and deep samples (t-test = 1.16, p = 0.226, Fig. 4B). As shown in Fig. 4B, the curve becomes asymptotic as the OTU number saturates, and each sample depth adds an increasingly smaller number of new OTUs, indicating adequate coverage for the environment. We did find that many Rickettsiales OTUs differ significantly in abundance between shallow and deep samples. For plotting purposes we selected those 44 OTUs whose p < 0.00001 (Fig. 4C).

Differences in the microbiome of A. cervicornis corals inhabiting different depths

Once the microbial taxonomic profiles of A. cervicornis naturally inhabiting different depths was established, we next focused our analysis on those OTUs present in all three samples of one depth and absent in all the samples of the other depth (core depth analyses). This analysis revealed that the taxa shared by all shallow samples were mainly Streptococcus, Haemophillus, Paucibacter and Porphyromonas and in the deep samples the shared taxa were dominated by Campylobacterales, Bradyrhizobium, and Lactobacillus. Overall, 15 OTUs were plotted representing the core taxa, of these only four OTUs were shared between all six coral samples with all other being shared between more than three samples (Fig. 5).

We then proceeded to determine which taxa changed significantly (selected OTUs with p ≤ 0.05) between shallow and deep samples by employing a log-likelihood ratio test. Upon performing this analysis, we found that 38 OTUs varied significantly between shallow and deep samples (Figs. 6A and 6B). Among the taxa that varied significantly between depths we found unclassified Rhodobacteraceae, Lactococcus, Comamonadaceae and Serratia to be more abundant in the shallow samples as compared to deep samples (Figs. 6A and 6B top panel). Conversely, we found a significantly higher abundance of Campylobacterales and Bradyrhizobium, L. plantarum and Erythrobacteraceae in the deeper coral samples (Figs. 6A and 6B). Taken together, our data suggest that only low dominance taxa changed significantly in abundance across depths.