16S rRNA amplicon sequencing reveals a polymicrobial nature of complicated claw horn disruption lesions and interdigital phlegmon in dairy cattle

Lameness represents an intractable problem for the dairy industry. Complicated claw horn disruption lesions, interdigital hyperplasia, and interdigital phlegmon are important lameness causing foot lesions. Their aetiology is multifactorial, but infectious processes are likely implicated in disease pathogenesis. Our aim was to investigate the bacterial profiles of these lesions using 16S rRNA gene sequencing of samples obtained from 51 cattle across ten farms in the UK. In this study, interdigital hyperplasia, interdigital hyperplasia with signs of interdigital dermatitis, interdigital phlegmon, complicated sole ulcers, complicated toe ulcers lesions, and complicated white line lesions were investigated; corresponding healthy skin control samples were also analysed. All diseased tissues displayed reduced microbial richness and diversity (as described by Chao1, Shannon, and Simpson alpha-diversity indices) compared to their healthy skin control samples. Our results confirm the association of Treponema spp with some of these disorders. Other anaerobic bacteria including Fusobacterium spp., Fastidiosipila spp. and Porphyromonas spp. were implicated in the aetiology of all these lesions with the exception of interdigital hyperplasia. Complicated claw horn disruption lesions, and interdigital phlegmon were found to have similar bacterial profiles. Such sharing of bacterial genera suggests many of the infectious agents detected in these foot lesions are acting opportunistically; this finding could contribute towards future treatment and control strategies.


Results
Swab samples were collected from 10 dairy farms in Cheshire and North Wales, UK. The numbers of each lameness causing lesion obtained from each different farm are described in Table 1. All collected samples were eventually submitted for sequencing.
After quality filtering processes, a total of 47,661,917 sequences were used for further analyses (Mean = 433,290, SD = 58,425 sequences per sample). Chao1, Shannon, and Simpson alpha-diversity indices for different lesion groups and their control samples are shown in Table 2. Beta-diversity was calculated through weighted UniFrac distances, and non-metric multidimensional scaling (NMDS) values were charted in 3D scatterplots. Bacterial composition of the samples did not indicate any clear clustering between farms (Fig. 1). On the other hand, there was a clear distinction between healthy skin control samples and IP, SU, TN, and WLD lesions. IH and IIH samples were however grouped together with the control samples (Fig. 2).
Relative abundances of the 15 most prevalent phyla are shown in Fig. 3. Main phyla for all the samples were Firmicutes followed by Bacteroidetes. Spirochaetes and Fusobacteria phyla were found in increased relative abundances in SU, TN, and WLD lesions compared to their control samples. In Fig. 4, relative abundances of the 15 most prevalent genera are shown.
The response screening analysis results provided a more comprehensive analysis of the differences at genus level between lesion and control samples; these results are charted in bubble plots. Fusobacterium spp., Porphyromonas spp., Helcococcus spp., Parvimonas spp., and Peptostreptococcus spp. were found to be significantly more prevalent in IP lesion samples compared to their healthy skin control samples. On the other hand, the prevalence of Corynebacterium spp., Acinetobacter spp., Christensenellaceae R-7 group, Ruminococcaceae UCG-005, Rikenellaceae RC9 gut group, and other genera of the family Corynebacteriaceae were significantly higher in healthy skin control samples compared to IP lesion samples (Fig. 5).
Analysis of weighted UniFrac distances suggested that the microbial communities of complicated CHDL and IP lesions were similar. For this reason, response screening analysis was also performed for all the complicated CHDL and IP lesions together comparing them to all their healthy skin control samples. Several genera showed significant differences between lesions and control samples, therefore only the genera with mean relative abundance higher than 0.01 were charted in order to decrease complexity. Porphyromonas spp., Fusobacterium spp., Treponema spp., Clostridiales Family XI, Fastidiosipila spp., Peptoniphilus spp., Ruminococcaceae UCG-014,  uncultured Bacteroidetes bacterium, and other genera of the family Porphyromonadaceae were found to be significantly more prevalent in lesion samples compared to healthy skin control samples. On the other hand, the prevalence of Corynebacterium spp., Atopostipes spp., Acinetobacter spp., Ruminococcaceae UCG-005, Ruminococcaceae UCG-010, Christensenellaceae R-7 group, Eubacterium coprostanoligenes group, Rikenellaceae RC9 gut group, and other genera of the family Lachnospiraceae were significantly higher in healthy skin control samples compared to lesion samples ( Fig. 9). In IH lesions, Erysipelothrix spp., Guggenheimella spp., Peptococcus spp., Petrimonas spp., Clostridiales Family XI, and other genera of the family Porphyromonadaceae were significantly more prevalent compared to their healthy skin control samples. However, their relative abundance was found to be low. The prevalence of Alistipes spp., Bacteroides spp., Christensenellaceae R-7 group, Eubacterium coprostanoligenes, and Clostridiales Family XIII AD3011 group were significantly higher in healthy skin control samples compared to IH lesion samples. Lastly, Peptoniphilus spp. were significantly more prevalent in IIH lesions than their healthy control samples.  Besides the lesion and healthy skin control samples, the ZymoBIOMICS ™ Microbial Community DNA Standard was used as a mock microbial community. This community comprises eight microbial and two fungal strains; seven out of eight bacterial strains were successfully assigned at species level, and in their expected relative  abundances. One bacterium was correctly identified at the genus level. Only 0.05% of the sequences from the positive control sample were assigned to other unexpected bacteria. Despite being in extremely low concentrations, and not visible on the agarose gel, three PCR negative controls were also sequenced. Two of these negative controls yielded less than 1,000 sequences, and one of them yielded 20,951 sequences; less than 5% of the average number of sequences obtained from the lesion and healthy skin samples. In addition to these, sterile swabs were processed as negative DNA extraction controls. Since their concentrations were very low and they were invisible in the gel, these samples were not further analysed.

Discussion
We have investigated the bacterial composition of complicated CHDL, IP, and IH samples and compared them to their healthy skin control samples. 16S rRNA gene sequencing revealed that all lesions are of polymicrobial nature rather than being associated with single taxa. Bacterial profiles of the healthy skin control samples were significantly different from lesion samples (with the exception of IH and IIH samples). In addition, healthy skin samples displayed an increased diversity compared to samples obtained from lesions. Recent studies on bovine mastitis have shown similar findings regarding microbial diversity of diseased versus healthy control samples 28 . As shown in previous studies, Treponema spp. appear to play an important role in the aetiology of some of these lesions 12,29 ; however, other anaerobic bacteria such as Fusobacterium necrophorum, Fastidiosipila spp. and Porphyromonas spp. were also found to be highly prevalent in most of the studied lesions. Complicated CHDL were previously shown to be associated with DD Treponema spp. using species-specific PCR primers 12 . Here, we have used universal primers that allowed us to detect other bacteria that could potentially be associated with the progression of these lesions. Our results indicated that Treponema spp. were statistically significantly more prevalent in complicated SU lesions comparing to their corresponding healthy skin control samples. This was not the case for complicated WLD lesions. Treponema spp. were prevalent only in a few of the WLD lesions; their relative abundance across all the WLD samples was not statistically significantly different to their corresponding healthy skin control samples. Results associated with IIH, IP and TN should be treated cautiously since we only managed to obtain a small number of samples from these lesions.
The cross sectional design of our study only allows us to describe a snapshot of the differences in microbiota profiles and therefore we cannot make assumptions regarding the importance of different taxa at different time points of disease progression. In some of our cases, a specific pathogen (e.g. Fusobacterium spp. in IP) could have been primarily responsible for the lesion which was then also colonized by other opportunistic pathogens. CHDLs are considered to be of non -infectious aetiology and what we describe here is most likely the secondary invasion of the exposed corium by a number of different opportunistic bacteria. Further, larger scale, longitudinal studies could better elucidate these diseases' aetiopathogenesis. A shotgun metagenomics approach could also be employed and would allow for a more in-depth investigation of the studied lesions' microbial communities.
Spirochetes of the genus Treponema were previously described as the predominant bacteria in lameness associated DD lesions [30][31][32] . In this study, they were found to have significantly higher relative abundance in SU and TN samples compared to their healthy skin control samples. T. medium, Treponema phylotype 18, and other Treponema spp. were significantly prevalent in SU and TN lesions. Treponema phylotype 18 was shown to be sharing 95% sequence identity with recognized Treponema species (T. putidum ATCC 700334, T. pedis T3552B, and T. denticola ATCC 35405) 33 . Toe necrosis lesions were also populated by T. denticola, and Treponema canine oral taxon 233. The detrimental effects of T. medium, and T. denticola were previously described in human oral lesions 34,35 , bovine digital dermatitis 9 , and contagious ovine digital dermatitis (CODD) 36 . However, it should be noted that Treponema spp. were only significantly prevalent in a small number of the TN and WLD samples within this study; in many of them we were not able to detect Treponema spp. sequences ( Supplementary  Figures 3, 4) and this would suggest that they are acting opportunistically in these cases that can also be complicated by other pathogens.
Porphyromonas spp. were formerly isolated from different human and animal infections 37,38 . In our study, all lesions were mainly dominated by P. levii which has the ability to synthesize anti-IgG 2 protease, and reduce macrophage activity 39 . P. levii is a well known opportunistic pathogen and is also commonly found in the vaginal discharge of metritic cows 40 . Complicated SU, TN, and WLD samples were also harbouring Odoribacter denticanis, which belongs to the same family (Porphyromonoadaceae) and was also previously associated with DD 41 . Fusobacterium spp., another well known opportunistic anaerobic pathogen, were previously found to be associated with lameness, particularly with DD and IP lesions 42,43 . Our results confirm its potentially important role in the development of IP but also indicated the significant presence of Fusobacterium necrophorum in all the other studied lesions. Invasion and colonization of tissues by Fusobacterium necrophorum is mediated by its virulence factors such as endotoxin, leukotoxin, and secreted proteases [44][45][46] . Its role in the aetiopathogenesis of dairy cattle metritis is also well known 47,48 .
Several types of bacteria in the Clostridiales order of the Firmicutes phylum were shown to have significant prevalence in both lesions and healthy skin control samples. One of these; Fastidiosipila spp. was previously reported in human osteitis 49 , and was found here in higher relative abundance in SU, TN, and WLD lesions (compared to healthy skin samples). Therefore, Fastidiosipila spp. could have a potentially important role in the aetiology of these lesions. Moreover, Clostridiales Family XI was significantly more abundant in SU, and WLD lesions (comparing to respective control samples); Helcococcus spp., and Parvimonas spp. from the same family were significantly more abundant in IP lesions. Murdochiella spp. were significantly more prevalent in WLD lesions and also belong to the Clostridiales family which is known to be associated with many human and animal diseases 50,51 . Complicated CHDL and IP samples displayed similar microbiota profiles. Therefore, all these lesions were analysed together and compared to all control samples. Similarly to individual analysis of these lesions, the prevalence of Fusobacterium necrophorum, Porphyromonas spp., Fastidiosipila spp., and Treponema spp. was significantly higher in lesion samples compared to the prevalence of these bacteria in healthy skin control samples.
Conversely, our weighted UniFrac analysis suggested that IH and IIH samples were grouped together with healthy skin control samples. Therefore, these lesions might not be a result of bacterial infection. Alternatively, the causative bacteria in these lesions might be mainly located deep inside the lesions, and thus biopsy samples would have been more appropriate in order to identify them. Viral infections could also be implicated. Erysipelothrix spp. observed in these lesions, are known to cause erysipelas in pigs and are associated with acute septicaemia, cutaneous lesions, abortions, or chronic infection causing endocarditis and arthritis 52 . Moreover, Guggenheimella spp. were previously reported in DD lesions 53 , but for other lesions in this study (IP, TN, WLD samples), Guggenheimella spp. were observed in higher relative abundances in healthy skin control samples. A notable presence of Porphyromonas spp., and Treponema spp. was also observed in IH and IIH samples, but the difference between lesions and healthy skin control samples was not statistically significant.
In conclusion, we have characterised the bacterial profiles of six different lameness causing lesions using 16S rRNA gene amplicon sequencing and have used multivariate analysis approaches to analyse our data and described significant differences between lesions and their control samples. Our results showed that most of these lesions were associated with a range of pathogens, which are most likely acting opportunistically. Porphyromonas spp. were more prevalent in IP, SU, and WLD lesions, Fusobacterium spp. were more prevalent in IP, and SU lesions, and Treponema spp. were more prevalent in SU samples (compared to their respective healthy skin control samples). Fastidiosipila spp., a pathogen not previously associated with lameness causing lesions in cattle, showed a noteworthy prevalence in SU, TN, and WLD lesions. Our findings could contribute towards future treatment and control strategies for these disorders.

Ethics and sampling. Ethical approval for the study was granted by the University of Liverpool Research
Ethics Committee (Reference Number: VREC547) all methods were performed in accordance with the relevant guidelines and regulations. As part of this cross-sectional study, a member of the research team accompanied five different professional foot trimmers as they visited clients' farms throughout June and July 2017. Farm visits were organised to coincide with therapeutic foot trimming, as opposed to routine preventative foot trimming visits, to increase success of finding targeted lesions. Lesions targeted within this study included; complicated SU, complicated WLD, TN, IH, IH with signs of interdigital dermatitis, and IP. Cases were defined following the ICAR Claw Health Atlas definitions 54 with the only additional requirement for a CHDL to be enrolled being the presense of obvious signs of infection. Claw horn disruption lesions without signs of infection (not complicated) were not included in this study. The foot trimmer was observed during routine foot trimming and when a targeted lesion was identified, a picture was taken, and a sterile swab used to obtain a sample from within lesion. In the case of IP and IH the skin was cleaned from any gross contamination with the use of a paper towel before sample collection. A second swab sample was taken from the plantar/palmar aspect of the affected foot, targeting the normal skin proximal and adjacent to the interdigital cleft, just above the heel bulbs. The area was cleaned from any gross contamination with the use of a paper towel before sample collection. Once obtained, swabs were placed in sterile tubes and labelled with lesion type, animal id, date, and whether the swab was lesion or control sample. Samples were transported on ice and stored at −80 °C until DNA extraction.  Pooling of PCR amplicons. Amplified libraries were pooled at 8 ng/μl, and 30 ng/μl concentrations for low-concentration and high-concentration samples, respectively. After measuring the concentrations of these two pools, they were mixed to 8.1 ng/μl for each amplicon with maximum volume from negative controls. The final pool was purified with Agencourt AMPure XP beads, and eluted in 30 μl to increase the final concentration.
Sequencing. Concentration and quality of the pooled PCR amplicons was evaluated with Qubit ™ dsDNA HS Assay Kit, and a fragment analyser (High Sensitivity NGS Fragment Analysis Kit, Advanced Analytical Technologies, Inc., Ankeny, IA, USA). The fragment analyser traces indicated that no primer-dimers or other non-target PCR products were present; thus, size selection was not required. A quantitative real-time PCR (qPCR) assay, designed to specifically detect adapter sequences flanking the Illumina libraries, was performed using an Illumina ® KAPA Library Quantification Kit (Kapa Biosystems, Wilmington, USA). This assay was used to quantify the number of DNA templates that had both adaptor sequences on either end and therefore those that would successfully form clusters on a flowcell for sequencing. Each pool of amplicon libraries was sequenced on a lane of the Illumina ® HiSeq2500 platform, in rapid mode with version 2 chemistry using sequencing by synthesis (SBS) technology to generate 2 × 300 bp paired-end reads. 15% PhiX fragment library was added to increase sample diversity.
Bioinformatics. Initial processing and quality assessment of the sequence data used an in-house pipeline.
CASAVA version 1.8.2 (Illumina) was used for base-calling and de-multiplexing of indexed reads to produce 102 samples across the two runs, in FASTQ format. PCR primer sequences and Illumina adapter sequences were trimmed by using Cutadapt version 1.2.1 56 . Sequencing errors were corrected to improve base quality in both forward and reverse reads using the error-corrected module within SPAdes sequence assembler (version 3.1.0) 57 . Read pairs were converted to single sequences that span entire amplicons using PEAR (version 0.9.10) 58 . Sequences with uncalled bases were removed. Size selection was applied to select sequences between 200 bp and 750 bp, thus removing sequences from potential PCR primer dimers or spurious amplification events. Sequences matching PhiX (E-value < 10 −5 ) were filtered out of the dataset.
For each sample, sequences that passed the filters were merged into a single file. This final sequence file with its own metadata file containing description for each sample, was analysed using a custom pipeline based on QIIME 1.9.0 59 . The Silva database (version 123) 60 was used along all the analysis, except for the phylogenetic tree alignment step which was performed using the GreenGenes precomputed 16S rRNA gene tree. The obtained amplicon sequences were sorted in groups to identify the sequence variability in each sample. This step has been performed using SWARM 61 with the strictest parameters (default parameters). Potential chimeric-sequences were discarded. Sequences were then aligned on the identified clusters to calculate the abundance of each cluster using a minimum similarity threshold of 97% for the entire length of the sequence in VSEARCH1.1.3. Taxonomic assignment of each cluster was carried out using QIIME and the RDP classifier 62 to match a representative sequence from each OTU to a sequence from the database. The most abundant sequence within each OTU's cluster was used as the representative sequence.

Statistical analyses.
To visualize the community composition of each sample, the OTUs abundance table obtained after normalization (at 260,000 sequences per sample) was used to summarise taxon abundance for a given taxonomic rank (from kingdom to species), using QIIME. The OTUs abundance table was also used to investigate the richness and evenness of the samples using the following estimators: total observed sequence variants (i.e. number of OTUs in the sample), Chao1, Shannon, and Simpson indexes. Comparisons of these indexes between the different groups of samples were performed using a series of t-tests.
To consider how the taxa composition changed in relation to groups for each metadata category, the rarefied abundance table was used to build pairwise sample distance matrices, using the Bray-Curtis 63 and the Weighted and Unweighted UniFrac 64 dissimilarity measures. Nonmetric Multidimensional Scaling (NMDS) analysis was performed on the obtained dissimilarity matrices and statistical significance for the obtained dissimilarity matrices calculated using analysis of similarities (ANOSIM). Subsequent P-values were calculated using Student's t-test.
Mean relative abundances of the fifteen most prevalent phyla and genera were charted for each lesion. Log fold changes (Log10) in relative abundance of the genera, with at least 1% abundance in one sample, were calculated for each lesion compared to their respective healthy skin control samples. Robust response screening analysis was performed in JMP Pro 12 (SAS Institute Inc., Cary, NC) in order to evaluate the differences in OTU (genus level assignments) relative abundance between complicated CHDLs, IH, and IP and their corresponding healthy skin control samples. A false discovery rate (FDR) correction was applied and statistical significance was declared at FDR LogWorth of 1.3 (equivalent of a P-value of 0.05). To facilitate data presentation, the genera with a Robust FDR LogWorth of 20 or more were adjusted to 20 (corrected Robust FDR LogWorth). Subsequently, the log fold change was plotted versus the corrected Robust FDR LogWorth value using bubble plot graphs in JMP Pro 12. Genera mean relative abundance defined the bubbles' size, and effect size was indicated by the bubbles' colouring 28 .