The Epidermal Microbiome Within an Aggregation of Leopard Sharks (Triakis semifasciata) Has Taxonomic Flexibility with Gene Functional Stability Across Three Time-points

The epidermis of Chondrichthyan fishes consists of dermal denticles with production of minimal but protein-rich mucus that collectively, influence the attachment and biofilm development of microbes, facilitating a unique epidermal microbiome. Here, we use metagenomics to provide the taxonomic and functional characterization of the epidermal microbiome of the Triakis semifasciata (leopard shark) at three time-points collected across 4 years to identify links between microbial groups and host metabolism. Our aims include (1) describing the variation of microbiome taxa over time and identifying recurrent microbiome members (present across all time-points); (2) investigating the relationship between the recurrent and flexible taxa (those which are not found consistently across time-points); (3) describing the functional compositions of the microbiome which may suggest links with the host metabolism; and (4) identifying whether metabolic processes are shared across microbial genera or are unique to specific taxa. Microbial members of the microbiome showed high similarity between all individuals (Bray–Curtis similarity index = 82.7, where 0 = no overlap, 100 = total overlap) with the relative abundance of those members varying across sampling time-points, suggesting flexibility of taxa in the microbiome. One hundred and eighty-eight genera were identified as recurrent, including Pseudomonas, Erythrobacter, Alcanivorax, Marinobacter, and Sphingopxis being consistently abundant across time-points, while Limnobacter and Xyella exhibited switching patterns with high relative abundance in 2013, Sphingobium and Sphingomona in 2015, and Altermonas, Leeuwenhoekiella, Gramella, and Maribacter in 2017. Of the 188 genera identified as recurrent, the top 19 relatively abundant genera formed three recurrent groups. The microbiome also displayed high functional similarity between individuals (Bray–Curtis similarity index = 97.6) with gene function composition remaining consistent across all time-points. These results show that while the presence of microbial genera exhibits consistency across time-points, their abundances do fluctuate. Microbial functions however remain stable across time-points; thus, we suggest the leopard shark microbiomes exhibit functional redundancy. We show coexistence of microbes hosted in elasmobranch microbiomes that encode genes involved in utilizing nitrogen, but not fixing nitrogen, degrading urea, and resistant to heavy metal. Supplementary Information The online version contains supplementary material available at 10.1007/s00248-022-01969-y.


Introduction
The epidermis is the largest organ in the body, and serves as the first line of defense against environmental influences [70]. Along with providing broad protections, the skin of the organism interacts with surrounding microbes and a specific microbiome occurs on different areas of the body and varies between species [1,28,32,92,98]. In the marine environment, the skin of elasmobranchs and teleost fishes have microbiomes that are distinctive to the surrounding water Michael P. Doane and Colton J. Johnson are co-first authors. [26,58,81,89] and a comparison across seven elasmobranchs, including thresher sharks (Alopias vulpinus), whale sharks (Rhincodon typus), leopard sharks (Triakis semifasciata), blacktip reef sharks (Carcharhinus melanopterus), nurse shark (Ginglymostoma cirratum), lemon shark (Negaprion brevirostris), round rays (Urobatis helleri), and southern stingrays (Hypanus americanus), confirmed that the elasmobranch epidermal microbiomes are species-specific [9,26,28,82]. However, the selection and maintenance of the microbiome on sharks remains an outstanding question, and may be the result of several processes, including, but not limited to, (1) biophysical properties of the epidermis, (2) epidermal secretions of the shark, and (3) interspecific interactions between the microbes.
The biophysical properties of the chondrichthyan fish epidermis are uniquely covered with dermal denticles that improve hydro-dynamism [53]. Dermal denticles have ridges and troughs arranged in overlapping patterns that alter the hydrodynamic properties of water close to the epidermis [90,99,103] and this biophysical property potentially affects microbial growth. In synthetic models mimicking shark skin, two bacteria, Escherichia coli and Staphylococcus aureus, showed differential growth patterns compared to a smooth surface. E. coli had high attachment rates to the structured surface, whereas S. aureus was not able to attach [13], therefore, suggesting that only specific types of microbes are able to attach to the shark skin. Further, modelling of microbial community dynamics showed that high surface structural complexity (like the skin of sharks) can reduce competition and increase coexistence between microbes [61]. Therefore, the high micro-structured nature of the shark skin surface may produce similar dynamics in the microbiome. The epidermis of Chondrichthyan fish is characterized by a lower level of secretory cells compared with Actinopterygian fishes [70] and the product is a thin heterogeneous outer mucous layer (5-8 μm thickness; pH between 6 and 7). The mucus contains high levels of proteins, with only a few proteins being identified [95] and has moderate to high disulfide concentrations that are suspected to provide mechanical stabilization of the mucus coating [71]. In addition to the physical structure of the dermal denticles, the epidermis is covered in protrusions that may stabilize the thin mucus layer [70] and microvilli-like apical protrusions, which may act as chemoreceptors [80], contributing to the highly structured nature of the skin surface.
The metabolism of elasmobranchs relies on ketone bodies, fatty acids, carbohydrates, and amino acids for energy production [88]. Amino acids are ketogenic precursors, oxidative fuel, and nitrogen suppliers [88]. Nitrogen metabolism supports the production of high amount of urea and trimethylamine N-oxide (TMAO) for osmoregulation [6,88]. Exocytotic activity is observable where the contents of the mucous vesicles below the apical membrane are transported onto the epidermal surface or into the thin mucus layer [70], suggesting potential processes for the products of host metabolism to be available to the microbiome. Sharks bioaccumulate toxins and heavy metals from lower trophic levels [33,67], absorb trace metals into their skin at a faster rate than teleost fishes [42][43][44][45][46], and have been shown to exhibit maternal off-loading to their offspring in mucus surrounding the developing embryo (Lyons and Lowe [63]). Therefore, metabolically derived compounds and harmful contaminants are being leached from elasmobranchs into their mucus and may affect the function of the skin microbes. For instance, the microbiome signature from thresher sharks (Alopias vulinpus) suggests potential heavy metal concentration at the skin interface signified by high abundance of heavy metal resistance genes [26].
Interactions between the microbes on the surface and those from the surrounding water are also expected to affect the characteristics of the microbiome. The low levels of mucus secreted by the sharks and minimal microbial biofilm development that occurred on the synthetic models of the shark skin, compared with the pronounced biofilm development on the smooth surface [13], suggest that microbial growth on the epidermis is suppressed, and the turnover of microbial communities may be high. Conversely, stability of a microbiome is important for the host's health and functioning, and rapid or unexpected variations in microbiome taxonomy can lead to a state of dysbiosis, resulting in a decline in organismal health [20,31,58]. Therefore, an understanding of how the microbiome varies through time is important, but there are limited investigations of microbiome stability over time in the marine environment. A few examples include the epidermal microbiome of humpback whales (Megaptera novaeangliae) and Atlantic cod (Gadus morhua) which varied seasonally [8,101], and similar seasonal changes occurred in the microbiome of the stony coral, Montastrea faveolata [51]. While the respiratory microbiome of bottlenose dolphin species (Tursiops truncatus and Tursiops aduncus) maintained taxonomic stability over 2 months, individual dolphins had significantly different microbiome taxonomic structures [56]. In an analysis over 8 months, captive dolphins had an individual signature of their lung microbiome and shared between 17 and 41% of the microbiome with other individuals [96]. In contrast, 77.6% of the microbes in the epidermal microbiome of thresher sharks (A. vulpinus) were shared across individual sharks [26], however, this study was only conducted at one time-point. Therefore, whether the epidermal microbiomes of marine vertebrate hosts are temporally stable or flexible remains an outstanding question in science.
Here we analyze both the taxonomic make-up and functional capabilities of the leopard shark (Triakis semifasciata) epidermal microbiome over three time-points measured across 4 years using the following objectives: (1) to describe the variation of microbiome taxa over time and identify those members which are recurrent (present on all individuals across all time-points); (2) to investigate the relationship between the recurrent and flexible taxa (those which are not found consistently across time-points); (3) to describe the functional compositions of the microbiome which may suggest links with the host metabolism; and (4) to identify whether the metabolic processes are shared across microbial genera or are unique to specific taxa.

Epidermal Microbiome Sampling
Leopard sharks (Triakis Semifasciata) form consistent aggregations in late summer, early fall at near-shore locations [76,77] along the coast of California making this species a good candidate to investigate temporal microbiome dynamics of a wild animal. Sharks in this study were sampled from La Jolla Cove, CA, USA (32° 51′ 12″ N, 117° 15′ 48″ W). Sampled individuals were all females, reflecting previously described sex-specific aggregation behavior exhibited by T. semifasciata along the southern California coast [76,77]. It is suspected that females aggregate in this location because they are pregnant [77], however, pregnancy status cannot be detected visually and requires the use of speciesspecific hormone analysis and ultra-sound equipment [21], which were beyond the scope of this study. All individuals were > 0.9 m, which is the estimated reproductive age for T. semifasciata [19]. We targeted at least three sharks per year and captured seven sharks in 2013, three in 2015, and eight in 2017; thus, the study was conducted at three time-points collected over a 4-year period. While T. semifasciata form aggregations, it is not always the same individual sharks that return each year [76]. In addition, T. semifasciata are small, and their fins are not strong enough to support a satellite tracker that could be used to pinpoint the position of an individual for recapture. Even if tagging was possible, seeing a tagged shark from the surface and recatching the tagged individual is unlikely. Therefore, recapture of the same individual shark across multiple years was not feasible and we opted to measure the microbiome across different sharks captured across the 4 years. The sharks were caught with a handline on baited barbless circle hooks and brought onboard a small skiff for processing using a scoop net; all sharks were released after processing. Microbiome samples from the sharks were collected using a microbial collection tool, which flushes the epidermis with ~ 250 ml of sterile seawater and dislodges the microbes. The microbial slurry is collected into the back end of the microbial tool, as we have conducted previously [26-28, 47, 55]. Each epidermal microbiome sample was taken from the dorso-lateral region in line with the dorsal fin and above the lateral line. The collected microbiome samples were filtered through a 0.22µm Sterivex filter (Millipore, Inc.) and stored in a − 20 °C freezer until processing.

DNA Extraction and Sequencing
Microbiome DNA was extracted using a cell lysis and spincolumn filtration method by Macherey Nagel Tissue Kit. The purified DNA was prepared for random shot-gun sequencing using the standard protocol from Swift 2S Plus Kit (Swift Biosciences). The DNA library was sequenced using a Mi-Seq Illumina sequencer with Mi-Seq v3 Reagent Kit. Samples were bar-coded and the T. semifasciata microbiomes were mixed in with a range of microbiome samples (e.g., kelp, fish, and seagrass microbiome samples) and run on several sequencing runs by the undergraduates in San Diego State University ecological metagenomics class [30].

Metagenome Annotation
The sequenced metagenomes were processed through PRIN-SEQ ± ± software for quality control, removing sequences less than 70 bp, those that had a single N (ambiguous base call), or a quality score below 25. MID tags were also removed with PRINSEQ ++ [11,86]. The resulting two fasta files were paired using FLASH software [64]. After pairing, metagenomes were uploaded into MG-RAST Version 4.0.3, which provided taxonomic and functional gene annotations of the metagenomes using standard cut-offs [5,69]. The sequences with the highest bit score were reported as we have done previously [39,73,75,97]. The taxonomic annotation outputs were filtered in MG-RAST to only include those from Bacteria and Archaea domains. Viruses and eukaryotes were not used in the analysis because they were underrepresented in the metagenome. The functional gene characterization was obtained from the SEED platform [79], similar to previous analysis [24,25]. Metagenomes were compared using proportional abundance, which is preferred to rarefaction [10,68,83].

Quantitative Evaluation of Leopard Shark Microbiomes
To identify whether the taxonomy of the T. semifasciata epidermal microbiome remained consistent or varied over time, a PERMANOVA on the relative abundance of each taxon level (from phyla to genera) was conducted across years using the entire dataset [2,83]. PERMANOVA is designed for non-parametric data, particularly when there is a larger number of variables compared with samples, unequal group sample size, and where there is no requirement for multivariate normalization because normalization is met by the permutation function [2,3]. All data were fourth root transformed [55], which balances the effects of a community structured on a few abundant species and a community structured on all species, and thereby influenced by the occurrence of the rarest taxa [16,17]. A similarity of percentages (SIMPER) analysis was used to identify the taxa contributing to the dissimilarity between years [14]. A principal coordinate analysis was performed on the taxonomic data to visualize the variation in the structure of the microbiomes over the years using a Bray-Curtis similarity, 100 = similar, 0 = no overlap. Last, a PERMDISP analysis compared the distribution of microbial taxa across each year around the year's group centroid [4,15]. These analyses were conducted initially on the entire dataset and then on the recurrent taxa.
We have identified two groups of microbial genera, including recurrent microbes (or core) defined as those genera that were present on every individual shark, on a presence/absence basis and non-recurrent microbes, defined as all other genera that were occasional members of the microbiome. The 19 most abundant recurrent (selected because they had a mean of > 1% of sequences across all metagenomes) and non-recurrent microbial taxa (abundances calculated within the recurrent/non-recurrent datasets, separately) were compared with a Bray-Curtis similarity analysis to identified groups of genera that had a similar proportional abundance across years, visualized with a dendrogram. A Pearson correlation was conducted to identify positive or negative relationships between these taxonomic groupings, displayed as a heatmap. The output of these two analyses was combined to show how the microbial genera related to one another. The analysis identified negative and positive correlations between the microbial groups found in the T. semifasciata shark epidermal microbiome.
We used metagenomics to describe the abundance of genes found by microbiome as a proxy for gene expression: although metagenomics does not measure which functional genes are being expressed at the point the sample was taken, it measures which functional genes are important for the bacteria in that environment [18,24,25]. There is a high level of correlation between the metagenomes and metatranscriptomes [34], where the abundance of a gene in metagenomes is a predictor of its expression level in the metatranscriptome and areas where the two analyses vary are associated with short-term changes in expression rather than bacterial functions that are under strong selective pressure and are well adapted to their environment [35,36,65]. For the leopard shark analysis, we used metagenomes to provide a functional potential, because of the stability of sampling DNA compared with mRNA in the marine setting. To identify whether the functional potential of the T. semifasciata epidermal microbiome remains stable or varies temporally, a PERMANOVA was conducted to compare variation in functional components across years. The SEED subsystem provides a hierarchal description in four levels including level 1 -major metabolic processes (such as carbohydrate metabolism); level 2 -metabolic pathways (such as di-and oligosaccharides); level 3 -gene clusters (such as 2-O-alpha-mannosyl-D-glycerate utilization); and level 4 -gene functions (such as mannosyl-D-glycerate utilization repressor MngR) [79] and all levels were tested. The SIMPER analysis identified potential functions that contributed to the dissimilarity between years. A principal coordinate analysis was performed to visualize the variation in the functional potential of the microbiomes over the years and a PERMDISP analysis compared the distribution of functional potential of T. semifasciata epidermal microbiome across years. All analyses were performed in PRIMER Version 6.1.15 and PERMANOVA ± Version 1.0.5 from PRIMER-E (2012). GraphPad Prism Version 8.3.1 was used to visualize the data.
We identified the functional potential of the 17 most proportionally abundant genera using the Krona plug-in [78] in MG-RAST [5,69]. Only 17 genera were reported, because they occurred in sufficient abundance to explore the functions of the genera. The sequences were further investigated using Kyoto Encyclopedia of Genes and Genomes (KEGG) Orthology (KO) database's KEGG mapper plug-in on MG-RAST [48,49], which provided identification of metabolic and biochemical pathways. The proportion of sequences within each genus for each functional group was calculated as the number of sequences within a genus divided by the total number of sequences within the function across the 17 genera. The proportion of sequences associated with each function in each genus was categorized into one of four quartiles: > 75%, 50-75%, 25-50%, 1-25%, and 0% and visualized using a heatmap. The heatmap was used to visualize the connection between the microbes and the sharks over time and identified whether there was switching of specific functional roles between genera in the microbiome.
The SIMPER analysis identified a similarity coefficient of 83.2 between 2013 and 2015, 80.4 for 2013 and 2017, and 84.6 for 2015 and 2017 (Bray-Curtis similarity index, 100 = similar, 0 = no overlap). Genera of the Gammaproteobacteria class contributed the most to the dissimilarity between the years, followed by genera in the Alphaproteobacteria and Betaproteobacteria, and the Flavobacteria class from the Bacteroidetes phylum ( Table 2). The variation in the relative abundance of Alcanivorax and Alteromonas contributed the most to the dissimilarity across years.

Temporally Recurrent Microbes
Microbial genera present in all 18 metagenomes were categorized as "recurrent" which included 188 of 597 genera. From the 188 recurrent genera, we determined how stable the relative abundance was across the three sampled time-points and identified a significant

Relationship Between Recurrent and Non-recurrent Genera
To investigate the relationship between the recurrent and non-recurrent genera (top 19 most relatively abundant genera not occurring across all sampling points), we first clustered the genera together with respect to their relative abundance across time-points based on Bray-Curtis similarity. Clustering identified three groups from the 19 dominant recurrent microbes labeled A, B, and C and five groups of the 19 most relative abundance non-recurrent microbes labeled A, B, C, D, and E (Fig. 4). The recurrent group A consisted of Pseudomonas, Erythrobacter, Alcanivoriax, and Marinobacter, which showed an even proportional abundance across the three time-points except Erythrobacter and Marinobacter which had higher proportional abundance in 2013. The recurrent group B consisted of Leeuwenhoekiella, Fig. 2 The microbial genera of the Triakis semifasciata epidermal microbiome clustered by year, analysed on proportional sequence distribution of all genera Altermonas, Gramella, Zunongwangia, Pseudoaltermonas, Flabacterium, and Maribacter which had higher proportional abundance in 2017. The recurrent group C consisted of six genera including Burkolderia, Novoshingobium, and Sphingopyxis which had even distribution across the years; Rugeria, which had proportional higher abundance in 2013 and 2017; and Sphingobium and Sphingomonas, which had higher proportional abundance in 2015. The non-recurrent genera in group A included Propionibacterium, Pirellula, and Plesiocystis and had lower proportional abundance in 2017. The non-recurrent microbial group B consisted of Aromatoleum and Chromobacterium and had low but even proportional abundance across all years. The non-recurrent microbial group C only consisted of Dokdonia, which displayed higher proportional abundance in 2017. The nonrecurrent group D included Acidiphillium, Azorhizobium, which had low but even proportional abundance across years, and Gluconobacter and Sphingobacterium which showed higher proportional abundance in 2017. The nonrecurrent microbial group E, including Leadbetterella with highest proportional abundance in 2015 and Klebsiella, Neisseria, and Thauara with even proportional abundance and Oligotropha, Frankia, and Leadbetterella displayed its highest proportional abundance in 2015. The last group included recurrent Xylella and Limnobacter and non-recurrent Endoritia and Providencia genera, therefore forming the mixed cluster on the dendrogram (Fig. 4).
We then investigated the correlative relationship of recurrent and non-recurrent groups (Fig. 4). The heatmap showed that Erythrobacter and Marinobacter (recurrent group A) were negatively correlated with all microbes in recurrent group B. Whereas Pseudomonas and Alcanivorax (recurrent group A) showed positive correlations with recurrent group B. Recurrent group A, however, was positivity correlated with non-recurrent groups A and B. In contrast, recurrent group B was negatively correlated with non-recurrent groups A and B but strongly associated with non-recurrent group C, Dokdonia. Recurrent group C was positively correlated with non-recurrent groups A, B, D, and E and the mixed cluster (Fig. 4). These associations between taxa will be investigated further in the functional analysis to identify metabolic processes that may drive the microbiome composition.

Functional Potential of the Shark Epidermal Microbiome
Annotated microbial sequences were categorized into 26 major metabolic processes (SEED level 1 subsystems), 153 metabolic pathways (level 2), 949 functional gene clusters (level 3), and 5,358 functional genes (level 4). The functional profiles were highly stable across the three time-points as evidenced by the high functional similarity between individual sharks (mean 97.7, Bray-Curtis similarity index). There was no significant difference between time-points in the functional potential of the microbiome across the three hierarchical levels (PERMANOVA: major metabolic processes, pseudo-F df = 2, 15 = 1.8, P(perm) > 0.05; metabolic pathways, pseudo-F df = 2, 15 = 0.525, P(perm) > 0.5; functional gene clusters, pseudo-F df = 2, 15 = 1.182, P(perm) > 0.1). The similarity of the level 1 functions across time-points was shown by the single dense cluster on the principal coordinate analysis that explained 64.4% of the variation in the first two axes (Fig. 5).

Microbial Gene Function Across Genera
The 17 dominant recurrent genera were used to investigate whether functional genes were equally dispersed across the genera or confined to a small subset of microbes (Fig. 6). Pseudomonas had the highest proportion of sequences annotated to all metabolic processes, including those involved in nitrate and nitrite ammonification and denitrification. Pseudomonas and Alcanivorax also conducted ammonia assimilation. These two genera were positively correlated ( Fig. 6) with most other recurrent genera, suggesting nitrogen metabolism within this group may facilitate survival of other T. semifasciata epidermal microbiome members. In comparison, Marinobacter also specialized in these two functions plus urea decomposition and was negatively correlated Pseudomonas and Alcanivorax. Ruegeria and Novosphingobium were the only other genera that had a high proportional abundance of genes associated with urea decomposition, which may have caused the negative correlation with recurrent group B. The genera within recurrent group B had a high proportion of genes associated with carbon metabolism, with Gramella having a high proportion of genes within TCA cycle, glycolysis and serine cycle.
Maribacter has a high proportion of genes in glycolysis, denitrification, and fatty acid biosynthesis. Many genera, including Pseudomonas, Alcanivorax, Erythrobacter, Leeuwenhoekiella, Altermonas, Gramella, and Flavobacterium, had protein biosynthesis and degradation, and fatty acid biosynthesis genes, which are associated with the shark ketosis dominated metabolism. Genera within recurrent group C, including Ruegeria, Sphingopyxis, and Novosphingobium, had a high proportion of genes associated with denitrification and fatty acid biosynthesis. All 17 genera, except one, had sequences that matched to copper, cobalt, and cadmium resistance and copper homoeostasis genes, suggesting the importance of heavy metal tolerance or manipulation for microbes living on the surface of sharks.

Discussion
The T. semifasciata epidermal microbiome was consistent between individual sharks with many genera shared across the three time-points spanning 4 years, one of the longest studies of a marine vertebrate microbiome to date. A recurrent set of 188 microbial genera were present on all individuals. However, there was flexibility in the proportional makeup of the microbiome genera across years, but the microbiome variability was consistent across years, suggesting similar ecological dynamics are occurring from year to year, despite some switching in dominate microbial genera. Regardless of the proportional abundance of each microbial genera, the functional potential of the microbiome was constant, suggesting that the epidermal microbes were responding to the unique environment of the T. semifasciata epidermis and that functional redundancy across microbial genera may allow the proportional switching of the genera observed across time-points.

Leopard Shark Microbiome Exhibits Specificity and Flexibility
The T. semifasciata epidermal microbiome has a high relative abundance of Pseudomonas, Erythrobacter, Alcanivorax, and Marinobacter; whereas the thresher shark (Alopias vulpinus) epidermal microbiome had a high relative abundance of Pseudoalteromonas, Erythrobacter, Idiomarina, and Limnobacter [26], and the Black-tip reef sharks (Carcharhinus melanopterus) epidermal microbiome a high relative abundance of Psychrobacter, Pseudoalteromonas, Rhodobacter, and Alteromonas [82]. While some genera were shared across shark species, the proportional abundance differed, suggesting that each elasmobranch microbiome has distinct characteristics which selects for differing microbiome members.
The highly shared microbiome across individual sharks within a time-point and flexible communities across timepoints coupled with evidence that shark species microbiomes are distinct may be caused by the dermal denticles [13,23] and metabolites secreted into the skin mucus as a result of the sharks metabolism. Dermal denticles alter the hydrodynamic properties of water close to the epidermis, reducing drag on the shark [91] and micro-organism growth [99,103]. A study that tested the recruitment of two lab microbial species found that only one species was able to successfully recruit, and biofilms did not develop [13,50]. These features are consistent across the microbiome of four shark species [26,28]. The epidermal microbiome of embryonic skates showed a dramatic change at the time of dermal denticle hardening, going from similar across all body parts to a unique microbiome on the epidermis [72], confirming the effect of the dermal denticles in microbiome structure.
Modelling of microbial community dynamics on surfaces with high structural complexity, like the conditions on a shark, suggests there is lower interspecific competition that promotes limited extinction, and results in a chaotic spatial distribution pattern of microbial species [61]. An expected outcome in this situation is a high level of shared microbes across individuals, but with varying abundance patterns. The T. semifasciata epidermal microbiome supports this expectation where 188 genera coexisted on all sharks across the three time-points. The flexibility in the relative abundance of those genera suggests that the epidermis creates multiple niches driven by the physical structure of the surface. Similarly, the skin microbiome of three lizard species was more diverse and had a larger core number of OTUs than that of six frog species, which may be associated with the structured nature of the rough surface created by the lizard species skin relative to the mucus covered frog epidermis and lower levels of antimicrobial peptides on the lizard skin [98]. The high diversity of the T. semifasciata microbiome may be further enhanced through potential interactions among groups of co-occurrence genera which drive niche partitioning and resource sharing (Fig. 4). For example, Pseudomonas and Alcanivorax showed positive correlations with recurrent group B members, including Leeuwenhoekiella, Altermonas, Pseudoaltermonas, Flabacterium, and Maribacter, suggesting that these two genera facilitate one another in some capacity through shared resources. Conversely, Erythrobacter and Marinobacter were negatively correlated with all microbes in recurrent group B, suggesting strong competition; this group was however positively correlated with non-recurrent groups A and B, suggesting transient microbes add to the microbiome's growth at different years.

Functions for Life on a Shark Skin
We investigated the gene functions present in the groups of microbial genera that may explain the positive or negative relationships between genera in Fig. 4. In addition, because host metabolism and microbial functions are linked [59,62], we identified functional genes in each genus that may be linked to the host metabolism and classify three broad functional groups, including generalists (recurrent group A), opportunists (recurrent group B) and specialists (recurrent group C) based on the presence of potential genes in these microbial genera. The group A recurrent microbes are typical marine species that have highly flexible genomes and growth strategies, including heterotrophic and photoheterotrophic growth. These organisms can utilize many of the nutrients anticipated to be present on the shark epidermis or excreted in the mucus, including proteins, fatty acids, and disulfides [70]. Pseudomonas, for example, is capable of mixotrophy [38] and the positive correlation of Pseudomonas with many other microbes (Fig. 4) may have been enhanced by traits, such as phosphate solubilization, and ammonia production [37,102]. Alcanivorax species, a highly abundant recurrent group A microbe, has the genes for reducing alkane hydrocarbons, and possessed many enzymes for β-oxidation of fatty acids, the glyoxylate bypass, andgluconeogenesis. These gene functions carried by Alcanivorax suggest they break down the complex hydrocarbons and fatty acid-rich metabolites and may facilitate other member of the microbiome by providing simpler metabolic products that can be utilized. They also include enzymes for synthesis of riboflavin and unsaturated fatty acids and cardiolipin [84]. Cardiolipin is a diphosphatidylglycerol lipid [85] and is more highly saturated in sharks than other vertebrates [87].
Erythrobacter species, a highly abundance group A microbe that was negatively correlated with many other microbes, are mostly aerobic anoxygenic phototrophic bacteria, containing bacteriochlorophyll a, but lack the genes of autotrophic CO 2 fixation pathway, thus photoheterotrophic metabolism requiring a supply of organic substrates. Several Erythrobacter species have genes encoding enzymes for glycolysis and the tricarboxylic acid cycle but lack genes for nitrogenase or nitrate reductase thus are reliant on other forms of organic nitrogen [52]. Erythrobacter carry heavy metal resistance genes, including resistance to lead, cadmium, zinc, mercury, nickel, cobalt, and arsenicals [104]. The behavior of the leopard sharks, being in shallow water with higher light conditions, may have provided the Erythrobacter with a competitive advantage over other microbial groups.
The last group A recurrent microbe was Marinobacter and these organisms are opportunistic generalists that switch rapidly between lithoheterotrophy to heterotrophy, in both anaerobic and aerobic conditions, in response to nutrient pulses [40]. Marinabacter can respire inorganic compounds that are usually found in metal rich environments [40], and we show they have many genes associated with metal resistance and transport.
The recurrent group B microbes, the opportunists, were negatively correlated with most other groups and had a high relative abundance of genes associated with carbohydrate metabolic processes. Gramella species degrades high molecular weight organic matter and encoding hydrolytic enzymes predicting a preference for polymeric carbon sources and have range of gliding motility genes that provide for surface adhesion [7], features that may enable the organism to live on the surface of elasmobranchs. However, these microbes had fewer genes in ammonia assimilation, promoting their coexistence with the more abundant species. Leeuwenhoekiella, a marine Flavobacteria, has heterotrophic growth, is pigmented, and often found in association with photosynthetic microbes [94], such as those described in the recurrent group A genera. Alteromonas are non-phototrophic, heterotrophic, with flexible genome for taking advantage of influx of nutrients, particularly polymers, and some degrade polycyclic aromatic hydrocarbons [66]. Alteromonas have high flexibility in the metal resistant genes, such as the copper-zinc-cobalt genes and efflux pumps [60] and these were highly abundant in this genus in the T. semifasciata microbiome. Pseudoalteromonas and Alteromonas utilized extracellular hydrolysis as the major decomposition pathway of peptides and released fragments of amino acids into the surrounding environment [57] and may have enabled them to outcompete Marinobacter which is unable to conduct extracellular hydrolysis processes [40]. However, these free amino acids may be used by other microbes, such as the non-recurrent group C and D microbes.
The group C recurrent microbes, the specialists, had low relative abundance overall but showed high relative abundance of a few metabolic pathways (Fig. 5). Ruegeria had a high proportion of sequences associated with urea metabolism, compared with its relatively low proportional abundance. Ruegeria are able to use trimethylamine (TMA) and trimethylamine N-oxide (TMAO) as an energy source to produce intracellular ATP [54]. TMAO and TMA are used by elasmobranchs for osmoregulation and TMAO-associated genes were identified in the microbiome. Sphingopyxis also had low relative abundance in the microbiome, utilizing carbon and nitrogen compounds as substrates for growth but are generally slow growing [74]. Serine, a pathway that was overrepresented in Sphingopyxis, can serve as a nitrogen source for the growth of some species, but not as a carbon source [100]. Sphingomonas are closely related to Sphingoyxis and both produce exopolysaccharides (EPS) [12]. The production of these EPS compounds may promote biofilm development on the shark skin. Therefore, each recurrent microbial group was providing a different set of functional genes that enhanced coexistence and reflected the products that may be associated with the shark mucus.

Gene Functions Linked to Host Metabolic Processes
The stability of the functional potential of the epidermal microbiome suggests the elasmobranch physiology and dermal denticle topography may both be drivers of microbiome structure, and we have developed a model to link key microbiome functions with shark metabolism (Fig. 7). Elasmobranchs rely on protein as the primary energy source and use ketogenic pathways that produce high ketone bodies and hydrocarbons [88]. The functions that were present in metagenomes suggests that the microbes were responding to the presence of these compounds in the minimal mucus that covers the sharks. For example, Alcanivorax and Marinobacter carried the genes to degrade hydrocarbons. There

3
are few descriptions of the compounds within the mucus of elasmobranchs, and these included high levels of proteins [95] with moderate to high disulfide concentrations [71]. A recent transcriptome analysis of the integument of cookie cutter sharks identified glyceraldehyde-3-phosphate dehydrogenase and fructose-bisphosphate aldolase A transcripts, both playing key roles in the glycolysis and gluconeogenesis [22] suggesting these products may be available to the microbial community. Protein metabolism along with nitrate and nitrite ammonification and urea decomposition with genes such as trimethylamine N-oxide (TMAO) reductase genes were identified in the microbiomes. The low levels of nitrogen fixation genes and the high proportional abundance of nitrogen metabolism and urea decomposition genes within the microbiome suggest the elasmobranch skin is a nitrogen-rich environment [29,41]. Proteasomes, metallocarboxypeptidase, aminopeptidase genes, protein degradation, and biosynthesis functional genes suggest a propensity of the microbes for protein degradation. An experiment on a human skin modelled showed that mixed microbiome communities had a greater positive effect on the skin function, including the mediation of the epidermal layer thickness, reduction in the number of actively proliferating cells, and increased filaggrin expression [59], suggesting that there is a two-way interaction between the microbes and host metabolic processes. The epidermal microbiome of embryonic skates showed a dramatic change at the time of dermal denticle hardening, going from similar to other body parts prior to denticle formation to unique post denticle formation [72], supporting our hypothesis that the dermal denticles affect microbiome structure. There was no evidence of dentine, hydrolysate, or collagen degradation in the shark skin microbiome, suggesting that the microbes are not breaking down the dermal denticles, which signifies a lack of negative impact on the sharks, at least in the form of degrading the dermal denticles. Therefore, a commensal connection between T. semifasciata and the epidermis microbiome may be established through the passive resource subsidization of these microbiome functions.
Cobalt-zinc-cadmium resistance and ton/tol outer membrane transport system associated with bacterial toxins [93] are highly prevalent functional genes shared across epidermis microbiomes of T. semifasciata (this study), A. vulpinus [26], Rhincodon typus, and Urobatis halleri [28]. Elasmobranchs bioaccumulate heavy metals via their food and shunt these compounds to the epidermis [33,67], and absorb heavy metals from the environment [45,45,46,46]. The epidermal microbiomes are responding to the presence of these heavy metals and future studies should measure heavy metal concentration in the skin in conjunction with the microbiome. The prevalence of heavy metal resistance genes may provide a biomarker of elasmobranchs health and identify those exposed to anthropogenic sources of heavy metals in the environment, providing an early warning sign of declining health.

Conclusions
The highly structured epidermis of elasmobranchs, in this case, T. semifasciata, promotes a diverse microbiome that has a flexible taxonomic makeup where many genera coexist. The microbial community was maintained by functional redundancy driven by the taxonomic flexibility and various microbes "flip-flopping" in abundance with one another, which kept the microbial niches filled across the years. The high relative abundance of genes involved in nitrification and denitrification, urea decomposition, and heavy metal resistance in the epidermis microbiome across all T. semifasciata individuals suggests a connection with the host metabolism potentially through passive subsidization of organic nitrogen, proteins, hydrocarbons, and heavy metals. Our results coupled with other host microbiomes studies, such as the skin surface microbiome of lizards, provide evidence that a structured surface is important for supporting and maintaining diverse skin surface microbiomes. Further efforts to understanding the effects of structure on driving skin microbiome communities will greatly enhance the ability to predict microbiome function of the largest organ in the vertebrate body.
Author Contribution MD collected all field samples, constructed the ideas, and helped write the paper; CJ analyzed the data and wrote the manuscript; EK, AT, AG, and SJ captured and prepared samples and provided feedback on the manuscript; MMM and LL provided sequencing and analysis; AN facilitated the catching of leopard sharks; RD provided the illustration; ED conceptualized the idea, analyzed and wrote the manuscript, and lead the project.
Funding Open Access funding enabled and organized by CAUL and its Member Institutions We acknowledge support from S. Lo and B. Billings. We acknowledge funding support from NSF Division of Undergraduate Education # 1323809 and NSF Division of Molecular and Cellular Science # 1330800.
Data Availability All data is publically available on the SRA database under BioProject PRJNA803589.

Declarations
Ethical Approval and Consent to Participate Animal handling and ethics were reviewed at San Diego State University through IACUC under permit APF #14-05-011D, APF #17-11-010D, APF # 18-05-007D. Sampling was conducted under state permit SCP #12847 and SCP #9893 from the California Department of Fish and Wildlife.

Consent for Publication
Not applicable.

Competing Interests
The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.