The skin microbiome of the common thresher shark (Alopias vulpinus) has low taxonomic and gene function &bgr;‐diversity

The health of sharks, like all organisms, is linked to their microbiome. At the skin interface, sharks have dermal denticles that protrude above the mucus, which may affect the types of microbes that occur here. We characterized the microbiome from the skin of the common thresher shark (Alopias vulpinus) to investigate the structure and composition of the skin microbiome. On average 618 812 (80.9% ± S.D. 0.44%) reads per metagenomic library contained open reading frames; of those, between 7.6% and 12.8% matched known protein sequences. Genera distinguishing the A. vulpinus microbiome from the water column included, Pseudoalteromonas (12.8% ± 4.7 of sequences), Erythrobacter (5. 3% ± 0.5) and Idiomarina (4.2% ± 1.2) and distinguishing gene pathways included, cobalt, zinc and cadmium resistance (2.2% ± 0.1); iron acquisition (1.2% ± 0.1) and ton/tol transport (1.3% ± 0.08). Taxonomic community overlap (100 - dissimilarity index) was greater in the skin microbiome (77.6), relative to the water column microbiome (70.6) and a reference host-associated microbiome (algae: 71.5). We conclude the A. vulpinus skin microbiome is influenced by filtering processes, including biochemical and biophysical components of the shark skin and result in a structured microbiome.


Introduction
Symbiotic interactions have extensive ecological consequences. Multicellular organisms, including diatoms, sponges, humans, amphibians, teleost fishes, brown algae and corals (Rohwer et al., 2002;Turnbaugh et al., 2007;Burke et al., 2011a;Fan et al., 2012;McKenzie et al., 2012;Larsen et al., 2013;Sison-Mangus et al., 2014) host a diverse consortia of microbial symbionts (Bacteria, Archaea, viral like particles and pico-Eukaryota), termed the microbiome (Lederburg and Mccray, 2001;Ursell et al., 2012). The microbiome is defined as the community of microorganisms that associates with a host organism and is hypothesized to have co-evolved with its host (Rosenberg et al., 2007). These communities of microorganisms are critical to the health of the host (Rohwer et al., 2002;Knowlton and Rohwer, 2003;Rosenberg, 2008;Consortium THM., 2012;McFall-Ngai, 2014) being involved in host functions, including nutrient supplementation (LeBlanc et al., 2013), successful development (Hooper et al., 2012), and disease susceptibility (Robinson et al., 2010). The members of the microbiome and their abundances inherently vary through space and time in response to ecological interactions with host and environmental constraints (Van Opstal et al., 2015); thus identifying structural community patterns of the microbiome and important drivers of fluctuation therein are critical to making a microbiomehost health link (Robinson et al., 2010).
Community structure can be defined as the amount of community overlap (b-diversity) from within the microbial community (Whittaker, 1972;Anderson et al., 2011). High b-diversity (low community overlap) across microbiomes of host individuals is expected if microorganisms from a source community experience dispersal limitation or other stochastic processes (random community membership), thus suggesting little influence by the host in selecting for its microbiome. Alternatively, low b-diversity (high community overlap) is the result of ecological interactions such as host selection and within microbiome processes (Robinson et al., 2010). If the host organism has co-evolved with its microbiome, it is reasonable to hypothesize that b-diversity of the taxonomic composition, gene function composition, or both will be low (high community overlap) among individuals from a single host species. Inherent to community structure is the amount of normal variance among metrics used to measure community structure. Deviation from normal microbiome structure can indicate altered host health status. For instance in humans, disease states, including obesity and irritable bowel syndrome, display gut microbiomes with structure that is altered relative to the 'normal' host community (Greenblum et al., 2012), which includes an increase in b-diversity (decrease in community overlap) and greater variance among the individual's microbial communities; thus community structure can be a useful metric in detecting changes to the host's functioning.
In the marine environment, the microbiome of the cutaneous (skin) (Larsen et al., 2001) and gut (Givens et al., 2015) regions of teleost fishes have been investigated. The taxonomic composition of the skin microbiome on fishes is primarily driven by the host fish species (Larsen et al., 2013). However, extreme environmental change does induce fluctuation in the skin microbiome composition, as observed in the Atlantic salmon (Salmo salar) microbiome as the fish migrate from rivers to the ocean (Lokesh and Kiron, 2016). In addition, biogeography has a strong influence on the microbiome composition from Atlantic cod (Gadus morhua) compared with season (Wilson et al., 2008). The microbiome is dynamically influenced by a variety of factors, and the observation that the microbiome is reproducible within designated sampling groups (i.e. within geographic locations or species), suggests intrinsic selective properties structure the microbiome. Host factors selecting for symbiotic skin microorganisms on fishes include the skin mucus (Ingram, 1980;Shephard, 1994;Scardino and de Nys, 2010;Gallo and Nakatsuji, 2011;Koropatkin et al., 2014) and the anti-microbial compounds in the mucus (Ingram, 1980;Subramanian et al., 2008). Skin mucus is a complex substance produced by goblet cells of the epidermis and is composed of mucin backbones which support a diverse array of O-linked glycans (Koropatkin et al., 2014). Skin mucus is common in teleost fishes, however Chondrichthyan fishes (sharks, rays, chimeras) have a thin layer of mucus.
Sharks have few documented infections (Ostrander et al., 2004) and a unique cutaneous structure. In the wild, sharks are commonly observed with open wounds, yet rarely observed as infected; instead these wounds heal quickly (Reif, 1978;Chin et al., 2015), suggesting that despite insult, the surrounding skin retains normal functioning, including the maintenance of a healthy skin microbiome (Reid et al., 2011). The mechanisms for the maintenance of a healthy shark skin microbiome may vary from distantly related teleost species. Sharks, like teleosts support a biologically active mucus layer (Luer et al., 2014;Tsutsui et al., 2014) at the skin interface, however mucus on sharks is inconspicuous relative to that of teleost fishes (Meyer and Seegers, 2012), suggesting sharks have reduced mucus secretion. Another difference between sharks and teleost fishes is that shark skin is textured (microscopic ridging) by dermal denticles that protrude through the epidermal and mucosal layer. A result of the textured surface is reduced fluid friction, which enhances swimming efficiency (Bechert and Bartenwerfer, 1989), and reduced microbial settlement (Sullivan and Regan, 2011;Reddy et al., 2011). The combined effects of the shark skin properties, i.e. the dermal denticles, which reduce fluid friction and the reduced mucus layer, may cause a selective skin microbiome and we will test this by comparing the taxonomic make-up and functional potential of the microbial communities collected from individuals of the same shark species.
We investigated the skin microbiome of the common thresher shark (Alopias vulpinus) which has a constant swimming behaviour, a distinct dermal denticle pattern and the presence of mucus. These features of the shark skin vary from other host organisms studied to date and thus our understanding of the effect of the biophysical and biochemical features of the shark skin on the microbiome are unknown. To date, most microbiome studies have relied on 16S rRNA sequencing, which considers only taxonomic structure. 16S approaches however, suffer biases that limit their comparative value across studies utilizing alternative 16S variable regions (for instance the V2-V4 region with the V5-V7 region). To begin linking the health of sharks to their skin microbiome, we evaluated the microorganisms (Bacteria and Archaea) from the skin of common thresher sharks using random shotgun metagenomics. Shotgun metagenomics provides a description of the taxonomic make-up and potential gene function of the microbial communities: two microbial characteristics which exhibit distinct ecological patterns Haggerty and Dinsdale, 2016). Oh et al. (2014), analysing the microbiome of the human gut confirmed shotgun sequencing to be valid for analysis of the taxonomic composition. Shotgun sequencing has been used to describe a variety of microbial communities including coral reef water Kelly et al., 2014;Haas et al., 2016), coral (Vega Thurber et al., 2009), Stromatolites (Breitbart et al., 2009 and the human body (Oh et al., 2014).  demonstrated the ability of shotgun metagenomics to distinguish between environments by describing the gene functional potential of microbes across nine biomes and Haggerty and Dinsdale (2016) reveal differential patterns in the taxonomic and gene functional potential in free living and particle associated microbial communities in the oceans. For a hostassociated microbiome, it is unknown whether taxonomic and gene functional potential patterns are distinct from each other like those observed for environmental microbial communities (Haggerty and Dinsdale, 2016).
If taxa and function of the host-associated microbiome are distinct, predicting functions from 16S amplicon data using tools like PICRUSt can produce misinformed conclusions. Therefore, the shotgun metagenomic approach provides a more complete picture of non-model hostassociated microbiomes.
We hypothesized that the features of shark skin cause a microbiome that; (i) is reproducible across individuals of the same species, (ii) has taxonomic and gene function composition that varies from water column microbial communities and (iii) is more tightly clustered (i.e. has low b-diversity) compared with microbiomes that are selected by one factor alone (i.e. mucus). We first contrasted the A. vulpinus skin microbiome to a reference coastal marine water microbial community, and then compared the amount of observed variation within the A. vulpinus microbiome to an algal species microbiome.

Results
The microbiomes from the skin surface of six A. vulpinus sharks were compared using shotgun metagenomics. The largest A. vulpinus library contained 1 651 492 reads while the smallest contained 431 043 and the average sequence length across all A. vulpinus libraries was 181 (S.D. 6 6) bps. The number of reads containing protein coding features (ORFs) ranged between 1 366 022 and 353 185 and of those only 7.5% to 12.8% aligned to known proteins (Table 1; Supporting  Information Table S1), suggesting a unique habitat in the shark skin boundary layer with novel microbial organisms and genes. Library size of the water column samples ranged between 204 479 and 56 742 reads and the average read length across all libraries was 367 (S.D. 6 75) base pairs. Proportions of identified protein coding features with alignment to reference protein databases ranged between 22.2% and 61.4%. The majority of aligned reads from both datasets were of bacterial origin (sharks: 98.3% -99.1%; water column: 96.3% and 99.1).
The ten genera that discriminated the A. vulpinus skin microbiome from the southern California water column accounted for a cumulative total of 38.5% of the total dissimilarity (Fig. 2); note that some of the discriminating   Distribution of genera discriminating A. vulpinus (grey) microbiomes from Southern California water column (blue) microbial communities. Contribution percent (Contrib%) is the percent each variable contributes to total variation observed between microbiomes. The total contribution of taxa driving divergence of A. vulpinus microbiome equals 30.3%. A. vulpinus and southern California water column and bar charts are the proportion of each taxonomic unit to its respective metagenomic library. Error bars are calculated as standard deviation. Contribution was calculated using SIMPER. Welch's t-test was used to determine group differences and a Benjamini -Hochberg correction applied to limit the false discovery rate.

Structural community comparison between the A. vulpinus microbiome, algal microbiome and water column microbial communities
Coefficients of similarity (100 -dissimilarity) were calculated to compare the structure of the shark skin microbiome to that of the water column and algal microbiome (Fig. 5A). MDS ordination indicates similar genera profiles across A. vulpinus microbiomes with an average similarity coefficient, which is a measure of b-diversity, of 77.6 6 S.E. 0.65 (100 -dissimilarity index; 0 no overlap, 100 high overlap), relative to the water column microbial community (70.0 6 1.60) and algal microbiome (71.5 6 1.15). The algal samples are examples of hostassociated microbiome that are formed by a single host structuring variable, i.e. mucus. The inclusion of an additional water column microbial community dataset revealed that despite temporal differences (Supporting Information Table S1), water column microbial communities are more similar to each other than to the shark or algae microbiomes (Fig. 5A). The amount of variance within microbial communities of the three locations (shark, water and algae) is different (PERMDISP, F 2, 12 5 5.99, p < 0.05). Pair-wise comparisons indicate the shark group has less variance than the water column microbial communities (t 5 3.19, p < 0.05) and algal microbiome (t 5 3.14, p < 0.05) and this difference is graphically displayed with MDS ordination (Fig. 5A).
The gene pathway composition of the shark skin microbiome is different from both the southern California and Point Loma water column microbial communities (Fig. 5B). Gene pathway overlap in the A. vulpinus microbiome was 87.4 6 0.59, which was similar to U. australis microbiomes (87.5 6 0.65) and greater than the water column microbial communities (82.9 6 1.19) (Fig.  5B). The within group variation among the three microbial community habitats was not different (PERMDISP, F 2, 12 5 1.65).
To investigate similarity of the microbial communities at a nucleotide level, we conducted a reference independent analysis. A pairwise BLASTn was performed to determine if unknown reads were similar across shark metagenomes. This analysis revealed that the A. vulpinus microbiomes had the greatest number of overlapping reads, averaging 62.4 6 1.1% (Fig. 6), compared with the algal microbiome (21.1 6 3.5%) and the water column microbiome community (8.4 6 2.2%).

Discussion
Microbiome community structure on the skin of Alopias vulpinus Here, we have characterized for the first time, the microbial communities from the skin surface of a shark species. We show that the shark skin microbiota is taxonomically and functionally consistent across several A. vulpinus individuals. Sharks shared common microbial genera and functions across individuals, and these microbes vary from those in the water column ( Fig. 5A and B). Taxonomic community overlap (100 -dissimilarity) was greater (lower b-diversity) within the shark microbiome compared with the water microbial communities and algae microbiomes. In addition, reference independent analysis, using a BLASTn pairwise analysis showed the highest sequence overlap occurring within the shark microbiomes compared to within algae microbiomes and water column microbial communities (Fig.  6). Taken together, these results indicate that sharks have a specific skin microbiome. The overlap in bdiversity suggests the shark skin microbiome has greater structure than the algal reference microbiome, a potential indicator that shark skin features have strong filtering properties.
Comparing microbial taxonomic diversity indicates that the shark, water and algal microbial communities share some common community characteristics while also exhibiting differences. Shark microbiomes have greater richness of genera than the water and algal microbiomes, however the shark microbiome has similar evenness (J 0 ). At the gene pathway level, sharks have greater diversity, richness and evenness. This discrepancy in taxonomic and gene function diversity suggests varying mechanisms involved in regulating the skin microbial communities. Higher richness of the shark microbiome compared to the algal microbiome, may reflect the shark skin encountering more microbial taxa due to higher spatial habitat usage, while the similarity in taxonomic evenness and gene function diversity metrics indicate filtering of the microbes. The shark microbiome is distinct from water column microbial community because it has a different ecological guild of microorganisms with different gene functions (Figs 2 and 4). The distinguishing genera in the shark microbiome are those included in the copiotrophic or eutrophic classification (Fierer et al., 2007;Haggerty and Dinsdale, 2016). The accompanying gene functions distinguishing the shark from the water column also correlate with the Fig. 6. Reference independent pairwise alignment analysis indicating the proportion of reads from each library that were identified across other libraries within the same group. Alignment was performed using BLASTn algorithm on unassembled reads. Error bars represent maximum and minimum observed values and box on either side of median are first and third quartile observations. copiotrophic/eutrophic microbial guilds (Haggerty and Dinsdale, 2016). In contrast to the shark, the most discriminating genera from the water column microbial communities are those considered phototrophic or oligotrophic, and gene functions are those from within broad functional groups identified as being positively associated with these ecological guilds. Thus, the types of microorganisms associating with the shark represent those expected to have surface associated life-histories.

Microbial taxa of the shark skin microbiome
The microbiome is challenged by environmental variables, host immune responses and microbial interactions, all of which need to be balanced to maintain the consistent community structure we observed across shark microbiomes. For instance, Pseudoalteromonas, which was overrepresented in the shark microbiome encodes gene functions which may promote healthy microbiomehost interactions (Holmstr€ om et al., 1999). An important feature of stable communities is weak ecological interaction (May, 1972;Coyte and Schluter, 2015). Pseudoalteromonas is known to synthesize antimicrobial compounds which may facilitate this within microbiome competition (Rao et al., 2006). In addition, synthesized anti-microbial compounds from this genera hinder settlement by eukaryotic marine fouling organisms (Holmstr€ om et al., 1999). Erythrobacter was overrepresented in A. vulpinus microbiome and has been reported from coral and algal microbiomes (Burke et al., 2011b;Morrow et al., 2014). Some Erythrobacter spp. perform oxygenic photosynthesis (Kolber et al., 2000) which could supplement a variety of compounds to other members of the A. vulpinus microbiome. Idiomarina spp., synthesize proteins associated with flagellated cell structures and performs heavy metal detoxification (Hou et al., 2004) and this group was overrepresented in the shark microbiome. Flagella, which aid in cell movement and nutrient uptake (Ottemann and Miller, 1997) provide a competitive advantage over other microbial cells. The capacity for Idiomarina to metabolize toxic substances may enhance the success of this group on the skin surface because sharks bio-accumulate heavy metals (Jeffree et al., 2006). The heavy metals may leech onto the skin surface, as occurs in coho salmon (Oncorhynchus kisutch) (Varanasi and Markey, 1977) and interact with the skin microbiome. Marinobacter, was the fourth genus to distinguish A. vulpinus skin from the water column. Members of this genus are known hydrocarbon degraders (M arquez et al., 2005) and produce lipopolysaccharides (K€ upper et al., 2006). Lipopolysaccharides mediate host inflammatory response and in some instances, reducing the inflammatory response; a possible function which would facilitate other microbial members to the shark skin microbiome

Distinguishing potential gene functions important to the shark microbiome
Overrepresented gene functions in the A. vulpinus microbiome included; (i) virulence, disease and defense; (ii) motility and chemotaxis; (iii) iron acquisition and (iv) membrane transport. The most abundant gene pathway overrepresented in the A. vulpinus microbiome was cobalt, zinc and cadmium resistance, a possible indicator of heavy metals that sharks are known accumulate in their tissues, particularly the integumentary tissues (Jeffree et al., 2006). In the microbiome of the alga, Ulva australis, cobalt, zinc and cadmium resistance was also overrepresented compared to the water column microbial community (Burke et al., 2011), but the A. vulpinus microbiome has three times the abundance of this gene pathway relative to the alga microbiome. Iron acquisition was the next greatest contributor to distinguishing the shark microbiome from the communities found in the water. Iron is an essential nutrient for basic metabolic processes but is a limiting resource due to its insolubility at normal physiological pH levels (Sandy and Butler, 2009). Microbiome constituents are hypothesized to exhibit high metabolic activity, competing with one another for limiting resources, including scavenging for iron. The overrepresentation of iron acquisition suggests an intrinsic selection property of the microbiome, whereby individuals with the ability to acquire iron more efficiently are suited to compete in the microbiome arena. Ton/tol system, another nutrient acquisition pathway was sixfold higher in the A. vulpinus microbiome compared to the water column microbial communities. The ton/tol system was shown to be positively correlated to phosphate levels in the water column on coral reefs system in the Line Islands , thus highlighting its importance for the acquisition of limiting nutrients. Another distinguishing potential gene pathways is n-phenylalkanoic acid degradation. Phenylalkanoic acids are natural products produced by Bacteria, Archaea and Eukaryotes, which have anti-microbial properties (Bowman, 2007). The ability of some microbiome members to degrade antimicrobial compounds increases habitable space among the microbiome. The high representation of sequences associated with heavy metal and toxicity reduction strategies and the presence of genera such as Idiomarina suggests that A. vulpinus may be carrying a greater load of heavy metals than previously realized (Lyons and Lowe, 2013) and/or these genes may be important for mediating the settlement of other organisms at the skin interface (Thomason et al., 1994). While the presence of these gene functions suggests a possible involvement in facilitating and providing beneficial services to the host, they may also be involved in biofilm formation. However, if biofilm formation is the primary outcome for the presence of the described gene functions, then taxonomic composition of the biofilm may not exhibit high community structure, as long as those microbial groups have gene functions necessary to produce the biofilm structure. Investigating shark species from several distinct locations would determine whether gene functional composition remains consistent despite location and if taxonomic composition of the skin biofilm is biogeographically dependent. In addition, if the presences of microbial heavy metal genes are linked to metal loads, investigating species of sharks that occupy different trophic positions or locations with varying heavy metal loads will help unravel this.

Sharks and other marine vertebrate microbiomes
Our understanding for the microbiome of other marine vertebrates is limited to 16S rRNA gene marker methods or amplicon sequencing; thus, lacking information for gene repertoires that indicate important gene functions of the microbiome. Larsen et al. (2013) used ribosomal internal spacer analysis (RISA) to compare the microbiome composition of six fish species and contrary to the authors' original prediction, found that the microbial community composition was more similar within individuals of the same fish species despite seasonal changes. Similar to fish, the microbial taxa associated with the epidermal surface of humpback whales (Megaptera novaeangliae) was similar across individuals, despite broad geographic separation of sampled whales (Apprill et al., 2014). The A. vulpinus microbiome had less diversity at the phylum level than either the fish or the humpback whales, being dominated by Proteobacteria, while fish microbial communities had high representation of both Proteobacteria and Firmicutes and humpback whales had high representation of Proteobacteria and Bacteroidetes.
The primary genera present in the shark skin microbiome varied from the microbiomes of fish and humpback whales. Pseudoalteromonas, Erythrobacter, Limnobacter and Idiomarina were abundant genera in the A. vulpinus microbiome, whereas Aeribacillus was reported on all teleost fish species and Pseudomonas was on all but one teleost species (Atlantic croaker; Micropogonias undulates). Janthinobacterium was a dominant genus in the skin microbial community from spotted sea trout (Micropogonias undulates), sand trout (Cynoscion arenarius) and pinfish (Lagodon rhomboids), but in low abundance in the A. vulpinus microbiome. Microbial communities from humpback whales were enriched in Tenecibaculum and Psychrobacter. The Psychrobacter genus was in low abundance on the A. vulpinus skin microbiome and Tenecibaculum was absent. Therefore, the skin of A. vulpinus appears to support a microbiome that is unique compared with other described marine vertebrates.

Proposed model for microbial interaction with the biologically active and topographically heterogeneous surface of shark skin
Our results indicate sharks have a specific microbiome, and the shark microbiome is less variable across individuals compared with the mucus driven microbiome of Ulva australis. Therefore, biophysical and biochemical properties of shark skin are hypothesized to produce the highly structured microbiome observed from the A. vulpinus skin microbiome. We propose a conceptual model that suggests that the combination of three selective features on the shark skin influence microbiome structure. The three selective features are, (i) a micro-patterned surface established by the dermal denticles, (ii) reduced hydrodynamic flow and (iii) antimicrobial properties within the skin mucus. At the epithelium/water interface, marine organisms harbour surfaces covered by secretions that support endogenous immunological proteins to reduce microbial activity (Shephard, 1994;Brown and Bythell, 2005;Meyer et al., 2007). Compared to teleost fishes, most sharks do not have high amounts of mucus at this interface (Meyer and Seegers, 2012) and in addition, scale like projections, called dermal denticles protrude through the epidermis, covering the thin mucus layer. Dermal denticles interact with water flow at the skin interface, reducing fluid friction through shear stress reduction (Bechert and Bartenwerfer, 1989;Bechert et al., 2000). A potential outcome of this interaction is a boundary layers that physically hinders diffusive momentum of particulate matter, or migration by microbial cells to the dermal denticle surfaces (Rusconi and Stocker, 2015). Interestingly, in controlled flow experiments microbial accumulation on surfaces is positively correlated with increasing shear stress. The reduction of shear stress by the dermal denticles may therefore slow the accumulation of water column microbial organisms to the shark skin surface relative to an untextured surface. In addition to creating a strong flow boundary, the surface pattern produced by dermal denticles has a negative influence on biofilm formation (Sullivan and Regan, 2011). Reddy et al. (2011) showed that patterned surfaces resembling shark skin reduced biofilm formation in Escherichia coli. A posteriori analysis of our data reveals that the number of gene functions associated with motility and chemotaxis are enriched in the A. vulpinus microbiome (1.8%; 1505 reads) compared to the U. australis microbiome (0.79%; 1334 reads). There was also a lower representation of quorum sensing pathways in the shark microbiome (0.004%; 3 reads) compared to the U. australis microbiome (0.014%; 24 reads) and conversely, chemotaxis pathways were an order of magnitude greater from the shark microbiome (0.19%; 158 reads) compared to U. australis (0.01%; 18 reads). Quorum sensing (QS) is the recognition of microbial compounds that are involved in biofilm production. However, biofilms exposed to high flow environments have QS functions which are turned off (Kim et al., 2016). If biofilms are mature on sharks, then the need for the QS phenotype of microbial individuals would be low and the genetic representation of genes in the QS pathway may be low in the community. Chemotaxis on the other hand may have higher representation in the shark microbiome because for cells to reach the surface under lower shear stress, the chemotaxis phenotype would be necessary.
In addition to the biophysical constraints of shark skin on the microbiome, the thin mucus layer against the skin surface is a biochemical filter to microorganisms. A Ctype lectin was identified in the skin mucus of the Japanese bullhead shark (Heterodontus japonicas) (Tsutsui et al., 2014) and another anti-microbial peptide called Pentraxin was found from the skin mucus of the common skate (Raja kenojei) (Tsutsui et al., 2009). Though the mucus of sharks has not been characterized, histochemical staining of the thin mucus layer shows a protein rich substrata (Meyer and Seegers, 2012).These observations support the hypothesis that skin mucus of sharks mediates microbial activity. The overrepresentation of reads identified to fatty acid metabolism and nphenylalkanoic acid degradation suggests a possible microbial response to these chemical insults imposed by the host.

Caveats
We recognize possible introduced biases through the analysis of metagenomic libraries constructed from different sequencing platforms; however, there is a paucity of microbiome studies using shotgun metagenomics from non-human hosts (Colston and Jackson, 2016). Several previous studies have incorporated metagenomic libraries sequenced across differing platforms and identified biological variance to be greater than inherent biases generated by sequencing technology Haggerty and Dinsdale, 2016). In addition, sequencing across various platforms generates reads of varying basepair lengths, of which may introduce additional bias (Danhorn et al., 2012). While this is a consideration, biological effects emerged as the drivers of divergent patterns in this study.
The ideal comparative host to a shark host is a marine fish species, however no study is publically available which describes the fish epidermal microbiome using shotgun sequencing (Colston and Jackson, 2016). The algal dataset presented here represents the only microbiome processed using similar DNA isolation methods and sequenced using shotgun metagenomics. In addition, by studying sharks of different trophic positions or individuals from regions influenced differently by coastal run-off, we can determine whether gene repertoires, such as heavy metal tolerance genes are linked with the health state of the host organism or are a normal component of the host's microbiome.
In this study, we describe between 8 and 12% of reads identified as protein coding across A. vulpinus samples with a match to the SEED database. Due to the high level of reads identified as having protein coding motifs without a similarity match to the SEED database, we suggest a high degree of novelty in the shark microbiome. The remaining reads may be novel microbial proteins from yet to be described microorganisms or may be those of phage/prophage origin. To account for the unknown, future studies will elucidate population genomes from the shark microbiome using nucleotide binning approaches. Sharks therefore present an opportunity for novel protein and genome discovery.

Conclusions
Our findings support the hypothesis that sharks have a structured skin microbiome with a unique complement of microbial taxa and potential gene functions. The microbiome is a product of both the host and its environment, and the high structure reported here suggests selective properties are imposed by the shark for the maintenance of a specific microbiome. Future studies should aim to decouple environmental versus host effects of the shark's skin microbiome. For instance, investigation of one shark species across several distinct populations will provide insight to the fidelity of microbiome members and their functions, helping to clarify whether the shark skin is selecting from the environment directly or if historical processes of recruitment and settlement to the skin surface are apparent.

Life history of and sample collection from Alopias vulpinus
Common thresher sharks frequent the northeastern Pacific Ocean along the western coast of North America, where the greatest densities are observed in the California Bight (Hanan et al., 1993); the geographic region extending from Point Conception, California, USA (348 N) to Cabo Colonet, Baja California, Mexico (318 N) (Carlucci et al., 1992). Juvenile common thresher sharks are reported to exhibit site fidelity in this region (Cartamil et al., 2010), but one individual tagged in the California Bight was caught in Bodega Bay, CA (388 33 N) (personal observation, 2015). Thus, juvenile common thresher sharks exhibit large potential ranges over the continental shelf. In this study, six juvenile common thresher sharks (< 300 cm total length) from the California Bight were sampled. Sharks were caught between Malibu, CA, USA (338 58 N, 1188 47 W) andNewport Beach, CA, USA (338 33 N, 1178 52 W) in September 2013 (Supporting Information Table S1) in association with the annual NOAA Southwest Fisheries Science Center juvenile thresher shark cruise. Sharks were captured on long-lines of approximately 2 km length with a soak time of less than one hour. Once caught, A. vulpinus individuals were placed in a cradle mounted along the gunwale of the boat. Microbial samples were collected from the skin surface at the base of the first dorsal fin. Six A. vulpinus individuals (n 5 6) were sampled using a 'supersucker' device (see Supporting Information Fig. S1) . The 'supersucker' is a twoway (blunt-ended) syringe that is prefilled with 0.02mm filtered sea water. The syringe modification allows water to flow over the shark skin surface and be recollected behind the plunger within the syringe body via fittings engineered to link the plunger end of the syringe and chamber. By pressing the chamber of the syringe against the shark, the chamber is isolated from the surrounding environment and filtered sea water is forcefully flowed over the shark surface and recollected via the vacuum created as the plunger is depressed. The microbial and viral free status of the filtered sea water was checked using fluorescent microscopy prior to sampling. Approximately 250 ml of filtered seawater was passed over the skin surface. The sampled water was passed through a 0.22 mm Sterivex filter (Millipore, Inc.), which captures microorganisms including small Eukaryotes, Bacteria, Archaea and some viral-like-particles. Parafilm wrapped Sterivex filters were stored dry at 228C for two days. Once back on land, samples were stored at 2208C until DNA extraction.  Table  S1). The La Jolla water column sample was from a region along the coast adjacent to a sea lion colony. In addition, the sample was taken post rain event. These five water samples are called the southern California (so_Cal) water samples and are well within the documented range of A. vulpinus. Due to temporal discrepancy between the southern California water column samples and shark samples, an additional water column dataset was included from Point Loma off San Diego, CA, USA (32839 0 59.56 00 N, 117814 0 50.62 00 W). The Point Loma water samples were collected at four time points (11/2013, 2/2014, 3/2014 and 11/2014), which coincide with the time of shark sampling. For both water column datasets, approximately 60 l of water was collected just below the ocean surface. Microbial cells from the seawater were concentrated using tangential flow filtration Haas et al., 2014), collected on 0.22 lm Sterivex filters (EMD Millipore, Billerica, MA), and frozen at 2208C until DNA extraction.

Water column sample collection
DNA extraction and metagenomic library preparation of A. vulpinus microbial samples DNA from microorganisms of A. vulpinus were extracted using a modified cell lysis and spin-column filtration method (Macherey-Nagel, Dueren, Germany) with slight modifications (Kelly et al., 2012). Sterivex filters were incubated with 720 ml of T1 buffer and 90 ml of Proteinase-K (2.5 mg/ml) overnight at 558C with rotation. Ends were sealed with luer-locking caps to keep material from leaking. Lysate was then removed with a 3 ml syringe with luer-lock fitting. Extraction was continued per manufacturer recommendations. Purified DNA was prepared for next-generation sequencing per Ion Torrent PGM (Thermo-Fisher Scientific, Carlsbad, CA, USA) protocol. Fragmentation and barcoding were conducted with IonXpress TM Plus Fragmentation Kit and Ion-Xpress TM Barcode Adaptors, respectively. Emulsion PCR was performed using the OneTouch (Thermo-Fisher Scientific, Carlsbad, CA, USA) with 8-cycles, supplemented with IonPGM Template OT2 200 Kit. IonPGM Sequencing Kit 200 v2 and Ion 318 Chip v2 were used for the sequencing step.

DNA extraction and metagenomic library preparation for southern California and Point Loma water column samples
Microbial cells from both water datasets were lysed in 100 ll of SDS (20%) and Proteinase K (2.5 mg/ml) with rotation at 558C for at least 6 h. For the southern California samples, standard phenol-chloroform procedure was used to isolate and purify DNA (Sambrook et al., 1989;. DNA purification for the Point Loma samples was conducted using a spin column purification kit (NucleoSpin Tissue, Macherey-Nagel, D€ uren, Germany), as described above. DNA for each dataset was quantified using a PicoGreen dsDNA assay (Life Technologies, Carlsbad, CA, USA). The southern California samples were sequenced with the Roche 454 GS Junior (Roche Diagnostics, CT), while the Point Loma samples with the MiSeq (Illumina, San Diego, CA), both at San Diego State University.

Library quality control and bioinformatic analysis
Quality control for all libraries (A. vulpinus, southern California and Point Loma) was conducted with PRINSEQ software (Schmieder and Edwards, 2011), removing all artificial duplicate reads and reads shorter than 70 basepairs. Contamination within the shark samples by A. vulpinus nuclear and mitochondrial DNA was checked by aligning reads to the closest sequenced nuclear genome, the elephant shark (Callorhinchus milii) (Venkatesh et al., 2014) and the closest mitochondrial genome, Alopias pelagicus (Chen et al., 2013), with the BLASTn algorithm and default settings (e-value of 10 25 ) (Altschup et al., 1990). Prior to alignment, reads were assembled into contigs using SPADES (Bankevich et al., 2012), with those reads forming each contig being assigned back to their respective metagenome using crAss (Dutilh et al., 2012). Reads that matched DNA of C. milii and A. pelagicus were removed from the microbiome analysis and used in further analysis (Doane et al. in review). High quality, unassembled reads were uploaded into MG-RAST (Meyer et al., 2008), an automated gene calling and protein prediction pipeline. MG-RAST calls taxonomic assignment with BLAT comparisons (Wilke et al., 2012) to the NCBI database and gene functional potential with BLAT comparisons to the SEED protein database (Meyer et al., 2008); each database using a hierarchical classification method. Sequence annotations were conducted using the following parameters: e-value 10 25 , 70% identity and > 30 bp alignment length.

Algal reference samples
A metagenomic dataset from the microbiome of the green alga, Ulva australis, was chosen for analysis because it represents one of a few marine microbiome surface datasets processed in a similar way to the shark microbiome samples (Burke et al., 2011a, b). At present, there are no publicly available shotgun metagenomic samples from fish, other sharks, or marine mammals. Microbial DNA sampling, extraction, processing and sequencing from the surface microbiome of U. australis is described in detail by Burke et al. (2009). Briefly, microbial samples were extracted from the surface of the algal thallus using a bead beating procedure and spin column isolations method. Shotgun sequencing was performed with Sanger sequencing, producing read lengths ranging from 718 to 802 base pairs. All algal microbiome samples are publicly available on the MG-RAST server (Supporting Information Table S1). This data was used to test hypothesis 3. Comparisons between metagenomes sequenced on different sequencers has been found to be valid (Qin et al., 2010;Luo et al., 2012;Dinsdale et al., 2013;Haggerty and Dinsdale, 2016).

Statistical analysis
The annotation of the sequences by MG-RAST provides two sets of data, (i) the number of sequences similar to each taxon from domain to genus and (ii) the number of sequences that are similar to each gene functional potential from major metabolisms (e.g. carbohydrate metabolism) to gene pathways (e.g. serine-glyoxylate cycle). Data were standardized to account for differences in library size, by calculating proportion of each taxa or gene function to the total number of annotated sequences in each library. Arcsine square-root transformation was applied to proportional data prior to statistical analysis. This data normalization was applied due to the heteroscedastic nature of the proportional data (Sokal and Rohlf, 1995). To address our three hypotheses, we have four microbiome datasets all annotated the same way, including (i) shark metagenomes (n 5 6), (ii) water column metagenomes (including southern California: n 5 5; and Point Loma metagenomes: n 5 4) and (iii) the algae metagenomes (n 5 4).
First, we identified whether the shark microbiome was reproducible across individuals by describing the relative differences in the taxa and gene functional potential on the shark skin microbiome (n 5 6) compared to the southern California water column microbial community (n 5 5) at the phylum and genus level for the taxonomy and the major metabolic categories and gene pathways for functional potential.
Second, we tested whether the shark taxonomic and gene function composition varied from water column microbial communities using SIMPER (Similarity Percentage) analysis. SIMPER analysis describes the contribution of each variable to the measured variation between two sample groups (Clarke and Warwick, 2001). The top 15 contributing groups identified in the SIMPER analysis (genera for taxa and gene pathway for function) were tested for differences using a two-sided Welch's t-test. A Benjamini-Hochberg correction (Benjamini and Hochberg, 1995) was applied to p-values to limit false discovery. In addition, distance matrices were generated with Bray-Curtis similarity (100 -dissimilarity score; 0 having no community overlap or high bdiversity, and 100 being a perfect overlap or low b-diversity) on transformed data. MDS ordination was performed on taxonomic and gene function matrices that included A. vulpinus (n 5 6), U. australis (n 5 4), southern California water (n 5 5) and Point Loma (n 5 4) water to visualize group compositional similarities.
Third, we described whether the shark microbiome is more tightly clustered or structured than a microbiome that is selected by one factor alone (i.e. algae mucus) or the water column. To identify the structure in the microbiomes the results of the Bray-Curtis (described above) distance matrices were compared using PERMDISP (permutational dispersion) procedure. PERMDISP was performed to calculate the linear distance of each similarity score within each group to the group centroid, which is a measure of bdiversity and this value was compared across groups in a permutational t-test (Anderson et al., 2006).The datasets used in this analysis include the shark metagenomes (n 5 6), the southern California water column metagenomes (n 5 5) and the algae metagenomes (n 5 4) with the Point Loma water column microbiomes excluded.
The structural differences between the microbiomes could be caused by varying diversity displayed by the microbial communities on each of the organisms and water column. Therefore, we described the diversity of the three datasets using Simpson's diversity (1-D) (Simpson, 1949), Margalef's Index (richness; d) (Margalef, 1958) and Pielou's evenness (J 0 ) (Pielou, 1975). Diversity metrics were calculated on the four datasets (shark microbiome, southern California and Point Loma water, and U. australis microbiome) using the number of sequences that matched to each genus and gene pathway. We also compared the Point Loma water microbial community diversity to the southern California samples to determine whether community diversity of the water microbiome was different. Each diversity metric used standardizes to the total number of annotated reads from each library. SIMPER, PERMDISP, MDS (Dinsdale et al., 2013) and diversity measures were generated using PRIMER/PERMANOVA package (Anderson et al., 2008); and Welch's t-test and p-value correction using STAMP software (Parks and Beiko, 2010).
To further test the three hypotheses and describe the structure in the shark microbiomes, we conducted an annotation independent analysis. This test was performed to assess microbial sequence overlap within metagenomes, from A. vulpinus, southern California water column and Ulva australis using BLASTn analysis. Each metagenomic library was built into a database, totalling 15 databases. Each metagenome was then aligned against each database, resulting in 225 BLASTn alignments. A cut-off e-value of 10 25 was used to identify sequence hits against the databases. This resulted in a matrix with the proportion of libraries overlapping, calculated as the total number of sequences that aligned against the database and standardized to the total number of sequences in that library (Cassman et al., 2012). The libraries aligned against the 'self' database confirmed the validity of this approach verifying the 1:1 match-up between database and library (data not shown).

Supporting Information
Additional Supporting Information may be found in the online version of this article at the publisher's web-site: Table S1. Library metadata by sample. NCB: northern California Bight regions. A. vulpinus (av), U. australis (ua), Point Loma water column (pl), southern California water column (so_cal). Fig. S1. Two-way blunt end syringe used for sampling microbial communities from the skin surface of A. vulpinus individuals. (a) Syringe system loaded with sterile seawater. (b) Syringe in position after sample acquired from A. vulpinus skin surface (water behind the syringe plunger). To obtain sample, syringe is pressed against the skin surface, creating a chamber inside the area labeled (Sison-Mangus et al., 2014). The plunger is depressed (Fan et al., 2012), and water flowed into chamber, over skin and up tubing, running alongside syringe body via tubing at (Turnbaugh et al., 2007). The water is deposited behind the plunger and stored until filtration through 0.22 mm Sterivex filter.