Shotgun Metagenomic Sequencing Reveals Functional Genes and Microbiome Associated with Bovine Digital Dermatitis

Metagenomic methods amplifying 16S ribosomal RNA genes have been used to describe the microbial diversity of healthy skin and lesion stages of bovine digital dermatitis (DD) and to detect critical pathogens involved with disease pathogenesis. In this study, we characterized the microbiome and for the first time, the composition of functional genes of healthy skin (HS), active (ADD) and inactive (IDD) lesion stages using a whole-genome shotgun approach. Metagenomic sequences were annotated using MG-RAST pipeline. Six phyla were identified as the most abundant. Firmicutes and Actinobacteria were the predominant bacterial phyla in the microbiome of HS, while Spirochetes, Bacteroidetes and Proteobacteria were highly abundant in ADD and IDD. T. denticola-like, T. vincentii-like and T. phagedenis-like constituted the most abundant species in ADD and IDD. Recruitment plots comparing sequences from HS, ADD and IDD samples to the genomes of specific Treponema spp., supported the presence of T. denticola and T. vincentii in ADD and IDD. Comparison of the functional composition of HS to ADD and IDD identified a significant difference in genes associated with motility/chemotaxis and iron acquisition/metabolism. We also provide evidence that the microbiome of ADD and IDD compared to that of HS had significantly higher abundance of genes associated with resistance to copper and zinc, which are commonly used in footbaths to prevent and control DD. In conclusion, the results from this study provide new insights into the HS, ADD and IDD microbiomes, improve our understanding of the disease pathogenesis and generate unprecedented knowledge regarding the functional genetic composition of the digital dermatitis microbiome.


Introduction
Bovine digital dermatitis (DD) is a prevalent disease that affects cattle worldwide [1][2][3]. It is associated with reduced milk production, decreased reproductive performance, high treatment cost and increased risk of culling, resulting in significant economic losses [4][5][6][7]. Current healthy skin (HS), which were harvested using a 0.6 cm diameter punch biopsy instrument (Biopsy Punch, Miltex Inc., PA). The punch biopsy was performed in the center of the DD lesion. Healthy skin samples were collected from the same hoof affected by DD, 2 cm from the center of the DD lesions. For all biopsy samples, the biopsy area was cleaned and disinfected with clorhexidine (Clorhexidine 2%; VetOne, Boise, ID) and a local anesthesia was performed by 3 subcutaneous injections of 2 mL of lidocaine (Lidocaine 2%; VetOne, Boise, ID). Following collections, samples were placed in sterile 2.0 ml microcentrifuge tubes, transported on ice and stored at -80°C until further analysis.

Whole genome-Shotgun sequencing
An aliquot of each extracted sample was normalized to 0.2 ng/μl. Once the normalization was done the samples could be used as an input to the Nextera XT DNA Sample Prep Kit (Illumina Inc. San Diego, CA). Tagmentation of samples was done using 1 ng of template, as directed by manufacturer. Following tagmentation, PCR amplification was done according to manufacturer's instructions using a unique combination of barcode primers (provided by manufacture) for each of the 16 samples allowing samples to be multiplexed. Following amplification, each DNA library had short DNA fragments removed using AMPure XP bead purification and then normalized through Library Normalization beads/additives. In addition, to prepare for cluster generation and sequencing, equal volumes of normalized libraries were combined, diluted in hybridization buffer and heat denatured, according to Nextera XT protocol. Finally, Pair-end sequencing was performed using the MiSeq Reagent Kit v3 (600-cycles) through the Illumina MiSeq platform.

Bioinformatics and statistical analysis
Raw data files from shotgun sequencing were de-multiplexed and converted to fastq using Casava v.1.8.2 (Illumina, Inc, San Diego, CA, USA). Fastq files were concatenated and uploaded to the MG-RAST [21] server for further analysis, using two annotation source data, SEED and M5NR. In MG-RAST, sequences were subjected to quality control, which includes dereplication (remove artificial sequences produced by sequencing artifacts), removing host specific species sequences (B. Taurus, UMD v3.0), ambiguous base filtering (removing sequences with >5 ambiguous base pairs) and a length filtering (removing sequences with a length of >2 standard deviation from the mean). Organism abundance was analyzed using a "Best Hit Classification" approach with a maximum e-value of 1 x 10 −5 , minimum identity of 60% and a minimum alignment length of 15 measured in amino acids for proteins and base pairs for RNA databases. A functional abundance analysis was performed using "Hierarchical Classification", with a maximum e-value of 1 x 10 −5 , minimum identity of 60% and a minimum alignment length of 15 measured in amino acids for proteins and base pairs for RNA databases. One sample with less than 1,000 sequences from IDD was excluded, yielding 15 samples for downstream analyses. To reduce the impact of experimental noise/error, the normalized data option of MG-RAST was used. Data underwent a normalization procedure, all values have undergone a log2-based transformation (log2 (x + 1)) followed by standardization within each sample and linear scaling (across all samples). Details can be found in the MG-RAST manual (ftp://ftp.metagenomics.anl.gov/data/manual/mg-rast-manual.pdf).
Phyla taxa abundance, normalized abundance of Treponema spp., flagellar motility, flagellum proteins, bacterial chemotaxis, genus distribution (unclassified, derived from Eukaryota), Candidatus Aemobophilus asiaticus and distribution of resistance to antibiotics and toxic compounds from HS, ADD and IDD were extracted from MG-RAST and analyzed using ANOVA in conjunction with Tukey test for multiple comparisons using JMP Pro 11(SAS Institute Inc., NC). A series of multivariable screening analyses using JMP Pro 11 were performed to determine which flagella protein were most important in differentiating the HS microbiome from the microbiomes of ADD and IDD. The false discovery rate (FDR) [22] was used to correct for multiple comparisons. Heatmap for phyla and function distribution from HS, ADD and IDD, and circular tree of Treponema strains from ADD and IDD lesion states were performed using MG-RAST. Results are presented as least square means follow by the standard error of the mean unless stated otherwise. Several

Sequencing data
Whole-genome sequencing was employed to compare HS, ADD and IDD from 16 dairy cows affected by DD. Table 1 illustrates the sequencing data obtained.

Normalized abundances
Comparison of the 12 most abundant phyla between HS, ADD and IDD are illustrated in Fig 1. Chordata, Firmicutes, Actinobacteria, Arthropoda, and Echinodermata were identified in HS as the most abundant phyla with Chordata being the most dominant phylum. Statistically significant differences at the phylum level were not identified between ADD and IDD and Chordata, Bacteroidetes, Firmicutes, Spirochaetes and Protobacteria were the most abundant phyla. The phylum Chordata, most likely bovine DNA, was more abundant (P<0.05) in HS when compared to ADD and IDD and the phylum Spirochaetes was significantly more abundant (P<0.05) in IDD when compared to HS (Fig 1). We also evaluated the phylum distribution of HS, ADD and IDD samples using a Heatmap (Fig 2). Spirochaetes and Bacteroidetes were highly prevalent in ADD and IDD samples followed by Proteobacteria and Firmicutes (Fig 2). The most dominant phyla in HS samples include Chordata, Actinobacteria and Firmicutes ( Fig  2). Analysis of the normalized abundance of Treponema spp. from HS, ADD and IDD demonstrated that T. phagedenis, T. denticola, T. vincentii and T. pallidum were the most abundant Treponema spp. and were significantly more abundant (P<0.05) in ADD and IDD compared to HS samples (Fig 3). Evaluation of the Treponema spp. across ADD and IDD identified a higher abundance (P<0.05) of T. maltophilum in IDD compared to ADD (Fig 3). The heatmap showing the functional composition of HS, ADD and IDD samples is illustrated in   Differences in functional categories were associated with motility/chemotaxis, respiration, iron acquisition, phosphorus metabolism, cell division and cell cycle and regulation and cell signaling (Fig 4).
To assess putative virulence factors associated with Treponema spp. and DD, the normalized abundance of genes associated with flagellar motility, bacterial chemotaxis and flagellar proteins from HS, ADD and IDD was evaluated and is illustrated in Figs 5 and S1 and S1 Table, respectively. The majority of genes with a role in flagella structure and synthesis in Treponema spp. were overrepresented in ADD and IDD when compared to HS. Two screening analyses were performed to identify the most prevalent and dominant flagella proteins for HS to ADD and IDD samples. The results indicated that 11 genes were distinct (P<0.00005) and highly prevalent in ADD and IDD compared with HS (Fig 6).
Comparison of normalized abundance of distribution of resistance to antibiotics and toxic compounds for HS, ADD and IDD is presented in Fig 7. It is evident that ADD and IDD have increased antibiotic resistance compared to that of healthy samples. Comparison of the antibiotic resistance profiles for IDD and ADD shows a similar trend for both samples.
The normalized abundance of genus distribution (unclassified, derived from Eukaryota) was compared for HS, ADD and IDD. The data indicates the presence of Acanthamoebas (free-living amoebas) in ADD and IDD (S2 Fig). Comparison of the normalized abundance of Candidatus Aemobophilus asiaticus demonstrated the presence of this bacterium in ADD (0.35 ± 0.12) and IDD (0.44 ± 0.14) but the bacterium was completely absent in HS.

Discussion
Metagenomic methods amplifying 16S ribosomal RNA genes (rRNA) have been used to describe the bacterial and archaeal phylogenetic composition of healthy and lesion stages of DD [15,16,19]. However, there is a lack of information regarding the functional composition of the metagenome associated with the disease. Whole genome shotgun sequencing of DNA is another approach to study microbial composition and diversity. It has previously been demonstrated that while data generated from 16S rRNA and shotgun sequencing are distinct [23,24], there is a huge overlap. Moreover, PCR bias is largely reduced in whole genome sequencing. The present study used shotgun metagenomics and the MiSeq Illumina sequencing platform to describe the microbial diversity. MG-RAST metagenomic analysis was used to assess gene function in HBS, ADD and IDD of dairy cows affected by DD. Bacteroidetes, Firmicutes, Spirochaetes, Proteobacteria and Actinobacteria were the most abundant bacterial phyla across all three samples types (HS, ADD, and IDD). The most abundant bacterial phyla identified in ADD and IDD samples were Spirochetes followed by Bacteroidetes, Firmicutes and Proteobacteria, which is in agreement with previous studies that have identified these phyla as the most abundant in DD lesions [15,16,19]. The two most abundant bacterial phyla in HS samples were Firmicutes and Actinobacteria, which is in line with data generated from other studies, identifying these three phyla as the most common phyla in healthy samples [15][16][17]. From the data obtained in our study it can be concluded that there is a distinct difference be tween the phyla types associated with DD samples compared to that of healthy samples.
DD was first perceived to be a polymicrobial disease, however more comprehensive studies have identified Treponema spp. as key components in the microbial etiology of the disease [18,19,25]. In our work we identified specific species of Treponema namely T. denticola-like,  T. phagedenis-like and T. vincentii-like as overrepresented in ADD and IDD compared to HS. This finding is consistent with work conducted by Evans et al. (2008) [25], where they identified T. phagedenis-like, T. medium/T. vincentii-like and T. putidum/T. denticola-like in bovine DD samples in abundances of 98%, 96.1% and 74.5%, respectively [25]. In the present study, we performed recruitment plot analyses where metagenomic sequences from HS, ADD and IDD samples were aligned against the genomes of T. denticola ATCC 35405, T. vincentii ATCC 35580 and T. pallidum subsp. pallidum str. Nichols. These recruitment plots determined that ADD and IDD sequences mapped with a high degree of similarity to T. denticola ATCC 35405 and T. vincentii ATCC 35580, signifying the presence of these specific Treponema spp. in diseased samples. Conversely, we established that a low number of sequences mapped to T. pallidum subsp. pallidum str. Nichols suggesting that our Treponema isolate does not belong to T. pallidum subsp. pallidum str. Nichols. Recruitment plot analysis was not performed for T. phagedenis and T. medium because their complete genomes were not fully sequenced at the time of data analysis. Circular tree analysis of the Treponema strains associated with ADD and IDD identified a high number of sequence reads for these specific Treponema spp. A small number of sequence reads were observed for other Treponema spp. previously associated with DD, including T. maltophilum, T. medium, T. pedis and T. paruluiscuniculi [16][17][18]25]. However, this low association may potentially be due to the small sample size used in this study.
Differences in functional features across the microbiomes of HS, ADD and IDD were identified through our MG-RAST analyses. We predominantly focused on further evaluating the most significantly altered functional features. Treponemas are characterized as Gram-negative, motile, anaerobic bacteria and are a member of the Spirochetes phylum. Typically Spirochetes have a unique structure consisting of a periplasmic flagellum located between the inner and the outer membrane, thus distinguishing these from other bacteria. Depending on the extracellular environmental stimuli, chemotaxis enables Spirochetes to move towards specific niches for colonization [26,27]. The ability to penetrate viscous media such as the skin tissue is an important virulence factor associated with pathogenic Spirochetes [28][29][30]. Lux et al. (2001) [31] established an in vitro tissue penetration assay using wild-type T. denticola, T. pallidum and T. phagedenis and demonstrated the ability of T. denticola and T. pallidum to penetrate tissue. However T. phagedenis was unable to overcome the epithelial barrier, suggesting that during the penetration process some Treponema spp. facilitate the entry of other Treponemes [31]. An association between Treponema spp. and deeper layers of skin affected by DD has been identified, thus further supporting the ability of Treponema spp. to penetrate skin tissue [16,18].
To gain a deeper insight into the virulence factors employed by DD-associated Treponema spp., comparison of the functional composition of HS compared to ADD and IDD was completed. Motility and chemotaxis were determined to be dominant functional features. A large portion of genes identified in ADD and IDD were related to flagella structure and to genes previously identified as playing a role in flagella synthesis in Treponema spp. Screening analysis to detect highly prevalent flagella genes in HS, ADD and IDD identified eleven genes, which did not exhibit a distinct difference between ADD and IDD. Although the function of FlaA sheath flagella protein is not well understood, it has been suggested that this protein increases rigidity of the flagellum and may permit more efficient swimming in viscous environments like the skin tissue [32][33][34]. Our results showed a high prevalence of FlaA sheath protein in ADD and IDD compared to HS samples.
Chemotaxis signals enable bacteria to move directly in response to environmental stimuli through a network of complex signal transduction pathways [26,35]. Sim et al., (2005) [35] constructed a model of T. denticola chemotaxis signaling pathways and demonstrated the role of specific proteins involved in chemotaxis signaling, specifically the role of CheA, CheW, CheY, Che X, Che R and chemoreceptors [35]. When the normalized abundance of bacterial chemotaxis genes in HS, ADD and IDD was analyzed, all genes were highly present in ADD and IDD compared to HS. Our findings strongly support a role for motility and chemotaxis in the disease process of DD. Another significant functional feature detected in ADD and IDD was iron acquisition and metabolism. Iron is a crucial component for growth and development of almost all bacteria [36][37][38][39]. Vertebrates are devoid of free iron; subsequently pathogenic bacteria employ multiple systems to sense and sequester iron in an iron-limited milieu and transport this across the outer membrane into the cell [40]. Our results indicate that during infection, Treponema spp. express genes possibly related to surface receptors. A putative role for these genes would be in iron acquisition, however further work to characterize these is required.
Foot bathing using disinfectant or antibiotic solution has been the most common method routinely used to prevent and control DD. However, there is controversy regarding foot bathing strategies and their success. Speijers et al. (2010) [41] compared the efficacy of several foot bathing methods and solutions in the treatment of DD. They demonstrated that 5% of CuSO4 per week was the optimum treatment [41]. The relative efficacy of a commercial disinfectant agent used twice a week compared to using 5% formalin and 10% CuSO4 was assessed and it was established that the commercial product had greater overall efficacy [10]. However, Laven and Hunt (2002) [8] reported that hoof bathing daily for seven days with 5% formalin and 10% CuSO4 was effective in controlling DD [8]. Although it appears that CuSO4 is effective, most published studies do not include a negative control group [42]. Here we have identified a higher normalized abundance in the expression of genes associated with resistance to antibiotic and toxic compounds in ADD and IDD, namely copper and zinc. The chronicity of this disease has been explained as a consequence of unsuccessful immune response and the ability of Treponemes to change from spiral to encysted form in deeper layers of the skin [38,[43][44][45][46]. This resistance may represent a mechanism by which Treponema spp. can persist in lesions. However, further research is needed to demonstrate the efficacy of CuSO4. Moreover, resistance was not observed against antibiotics such as oxitetracycline and lincomycin HCl, which are commonly used to topically treat DD [47].
A novel bacterium Candiatus Aemobophilus asiatiaticus has been identified in DD lesion samples and it is thought to play an important role in the pathogenesis of the disease [16]. It has also been reported that this bacterium has a symbiotic relationship with amoebas [48,49] and thus a second possibility is that Candiatus Aemobophilus asiatiaticus is an indirect predictor of amoebas [16]. Although in the present study we demonstrated the presence of this bacterium in ADD and IDD, we confirm the presence of Acanthamoebas in ADD and IDD supporting the hypothesis that Candiatus Aemobophilus asiatiaticus indirectly indicates the existence of free-living amoebas without necessarily being involved in DD pathogenesis. Amplicon sequencing of highly conserved genes, such as segments of the 16S rRNA gene, have been used in order to determine which organisms were present in different stages of DD lesions [15,16,18,19]. These studies, using a large sample size, identified differences among the microbiome of ADD and IDD. However, the present study did not identify differences among the microbiome of ADD and IDD and are likely attributed to the restricted number of samples and the experimental method used. Despite these differences, the technique used in this study revealed the potential metabolic profiles of a microbial community. How this potential is translated to functional activity remains unknown and further work needs to be completed to identify which genes are differentially expressed in ADD and IDD measured by metatranscriptome.

Conclusion
Here, we report the first shotgun metagenomic study of the microbiome of bovine digital dermatitis lesions. The MiSeq Illumina sequencing platform was used to assess the microbial diversity of DD and T. denticola-like, T. phagedenis-like and T. vincentii-like were identified as the dominant strains present in ADD and IDD. MG-RAST metagenomic analysis identified virulence factors including motility/chemotaxis and iron uptake. Additionally, an increase in the abundance of copper and zinc resistance genes has been reported, suggesting that further research is necessary to optimize the foot bathing technique for the control and treatment of DD. Overall, our findings provide significant information to better understand disease pathogenesis as well as to settle preventative and efficient treatments against microbes involved in DD pathogenesis.