Skin mucus in fish is the first barrier between the organism and the environment but the role of skin mucus in protecting fish against pathogens is not well understood. During copulation in sharks, the male bites the female generating wounds, which are then highly likely to become infected by opportunistic bacteria from the water or from the male shark’s mouth. Describing the microbial component of epithelial mucus may allow future understanding of this first line of defense in sharks. In this study, we analyzed mucus and skin samples obtained from 19 individuals of two shark species and a stingray: the nurse shark (Ginglymostoma cirratum), the lemon shark (Negaprion brevirostris) and the southern stingray (Hypanus americanus).   Total DNA was extracted from all samples, and the bacterial 16S rRNA gene (region V3-V4) was amplified and sequenced on the Ion Torrent Platform.  Bacterial diversity (order) was higher in skin and mucus than in water.  Order composition was more similar between the two shark species.  Alpha-diversities (Shannon and Simpson) for OTUs (clusters of sequences defined by a 97% identity threshold for the16S rRNA gene) were high and there were non-significant differences between elasmobranch species or types of samples. We found orders of potentially pathogenic bacteria in water samples collected from the area where the animals were found, such as Pasteurellales (i.e. genus Pasteurella spp. and Haemophilus spp.) and Oceanospirillales (i.e. genus Halomonas spp.) but these were not found in the skin or mucus samples from any species.  Some bacterial orders, such as Flavobacteriales, Vibrionales (i.e. genus Pseudoalteromonas), Lactobacillales and Bacillales were found only in mucus and skin samples. However, in a co-occurrence analyses, no significant relationship was found among these orders (strength less than 0.6, p-value > 0.01) but significant relationships were found among the order Trembayales, Fusobacteriales, and some previously described marine environmental Bacteria and Archaea, including Elusimicrobiales, Thermoproteales, Deinococcales and Desulfarculales.  This is the first study focusing on elasmobranch microbial communities. The functional role and the benefits of these bacteria still needs understanding as well as the potential changes to microbial communities as a result of changing environmental conditions.

Sample collection

Mucus and skin tissue samples were collected from 19 apparently healthy individuals (no visible wounds, normal swimming activity); 14 of them from animals captured in Bimini, Bahamas (25°43¢59 N, 79°14¢60 W): four corresponded to juvenile nurse sharks (Ginglymostoma cirratum), six to juvenile lemon sharks (Negaprion brevirostris), and four to adult southern stingrays (Hypanus americanus). Samples from an additional five adult nurse sharks were collected at Oceanario from Islas del Rosario (CEINER), in the Colombian Caribbean (10°10¢30 N, 75°45¢00 W).  For each individual, we obtained a sample of skin tissue and mucus, following sampling protocols approved by the Animal Care Committee of Universidad de los Andes (CICUAL) (Bogota, Colombia). The skin tissue sample was cut, using a sterile blade for each specimen, from the posterior part of the dorsal fin (1 cm³ or less) and the mucus from the skin surface, using a sterile 1.5 ml microcentrifuge tube to scrape the skin surface, ideally filling at least half of the tube.  Animals were manipulated for approximately 5 minutes and immediately released.  A water sample was also collected in sterile 15 ml tube from the sampling location of each individual.  Thus, three samples were associated with each individual, for a total of 57 samples.  The individuals were captured and raised slightly above the surface of the water, so that the samples could be taken outside the water, while the animal could continue breathing. Skin samples were preserved in ethanol 90%. All samples were maintained at 4 ºC for less than one week, until processing.

DNA Extraction and PCR amplification

DNA was extracted from the entire sample collected for all samples. The Tissue and Cells DNA Isolation Kit (MoBio Laboratories, Inc.) was used, following the manufacturer instructions. Water samples were filtered through a 0.8 mm cellulose nitrate filter before DNA extraction.  The primers 515f and 806r were used in order to amplify the region V4 from the bacterial and archaea 16S rRNA gene using the primers 515F (5´-GTGCCAGCMGCCGCGGTAA-3´) and 806R (5´ggactahvgggtwtctaat-3´) (Caporaso et al., 2010). PCR amplification conditions were as follows: an initial denaturation at 94 °C for 3 minutes, followed by 35 cycles of denaturing at 94 °C for 45 seconds, annealing for 45 seconds at 50 °C and extension for 45 seconds at 72 °C, followed by a final extension of 20 minutes at 72 °C. A negative PCR control was always included to reduce the chance of contaminant amplification.  Successful amplification was confirmed on 1 % agarose gel.  

Ion torrent library preparation, quantification and sequencing

From the 57 samples, 32 were used to construct libraries (Supplementary Table 1). Samples were chosen depending on their final DNA concentration, once the PCR products were cleaned using magnetic beads and run on a 1.5% agarose gel.   Only the samples that had a clear strong band were used for library construction. Two libraries, each with 16 barcodes, were prepared using the protocol Ion Xpress™ Plus gDNA Fragment Library Preparation (Life Technologies). Libraries were quantified with the Qubit kit.  Templates were prepared following the Ion PGM™ Template OT2 200 Kit (Life Technologies) protocols. Libraries were prepared for sequencing using the protocol Ion PGM™ Sequencing 200 Kit v2 (Life Technologies). Libraries were pooled to equimolar concentration and loaded on two Ion 316 chips and sequenced in the Ion Torrent PGM (Life Technologies).

16S rRNA datasets used in this manuscript with accompanying metadata has been submitted to Dryad.

Bioinformatic and statistical analyses

Sequences were separated by barcodes directly by the Ion Torrent PGM and saved by the ion reporter in different files; sequence quality was analyzed using FastQC (Andrews, 2014). The file format was changed from Fastq to Fasta.  Demultiplexing was conducted by comparing the mapping file of the chip with the files containing the sequences. For the core diversity analysis, Qiime2 (Bolyen et al. 2019) was used via command line using the moving pictures tutorial as reference. The files were imported as “MultiplexedSingleEndBarcodeInSequence” and demultiplexed using “cutadapt”, eliminating sequences shorter than 50 bp. The sequences went through DADA2 (Callahan et al. 2016) for quality control to delete sequences with lower qscore than 20 and then the remaining sequences were aligned de novo with align-to-tree-mafft-fastree. In parallel, the sequences were clustered into OTUs used to perform non phylogenetic analysis. The rooted tree obtained with fasttree2 (Price et al. 2010) was used to perform an alpha rarefaction with a 1000 sequence depth. For taxonomic assignment, analyses were performed on the Galaxy online platform (Afgan et al. 2016) following one amplicon data workflow on Mothur v.1.28.0 (Schloss et al., 2009). This workflow started by merging all read files into group files. Group files were identified as samples from each of the three elasmobranch species and also as type of sample (skin, mucus or water). The next step of the workflow identified unique sequences and generated a file with these sequences and a second file in which the number of each unique representative sequence was kept. Following this, reads were filtered based on quality and length. Parameters to remove low quality sequences (quality control) was for those with less than 20 Phred score and shorter than 50 bp. (minimum length) followed by a step to remove poorly aligned sequences and chimeric sequences. Finally, reads were clustered based on their degree of similarity, with a minimum of 97% identity threshold and aligned to the Silva V4 reference database (Quast et al., 2013), followed by a classification step into taxonomic categories (order, family, genus and species when possible). 

Rstudio version 1.1.463 (R Development Core Team, 2010) was used (Wickham, 2009) for alpha (a) diversity analyses (Simpson and Shannon) (package Vegan) (Oksanen et al. 2015) which were conducted for OTUs, using the number of OTUs per sample, comparing among species (N. brevirostris, G. cirratum juveniles and adults and H. americanus) and among sample types (tissue, mucus and water). In this analysis, OTUs with less than 0.2% presence were not included. A Shapiro-Wilk normality test was conducted to evaluate normality among samples belonging to each elasmobranch species (Supplementary Table 1), including the additional category of adults and juvenile for nurse sharks, or to each category of sample type before performing any statistical tests.  Since results fell outside the normality assumption, a Kruskal-Wallis test was used to evaluate whether a diversity was significantly different among elasmobranch species or among sample type.  To estimate beta (b) diversity (Bray-Curtis dissimilarity index and a Principal Component Analysis (PCA) the taxonomic category “order” was used.  Venn diagrams (package DVenn) were used to visualize shared orders among elasmobranch species and among sample types.

In order to find co-occurrence between different bacterial and/or Archaea orders a correlation matrix was created in R using the Spearman´s co-efficient as in Ju et al (2013). Correlations had to be stronger than 0.6 with a p-value < 0.01 to be considered to have a significant co-occurrence with other orders. All orders, including those with less than 0.2% presence were included in the co-occurrence analysis. A chord plot was created to visualize the relations between the different orders.

The raw dataset is available to the public on Figshare.

Caballero, Susana; Galeano, Ana Maria; Lozano, Juan Diego; Vives, Martha (2020): Description of the microbiota in epidermal mucus and skin of sharks (Ginglymostoma cirratum and Negaprion brevirostris) and one stingray (Hypanus americanus).