Metagenomic assessment of the rumen resistome, mobilome and stress response genes in smallholder dairy cattle in Kenya

Smallholder dairy cattle rumen microbiotas are subjected to a wide range of antimicrobials as well as sudden fluctuations in diets. As such, they develop an enormous reservoir of resistant genes, mobilome and stress response genes. However, information on metagenomic reactions to such dietary variations, especially for cattle reared in the tropics, remains largely unexplored. This meta-analysis was conducted to assess if antibiotic and toxic compound resistance genes (ARGs), stress response genes and bacterial phages, prophages and transposable element genes were present, and to what extent, in three dairy cattle genotypes (Friesian, FriesianXJersey crossbreed, Jersey) reared in a farm that practiced judicious use of antimicrobials. Potential bacterial hosts to these genes were also explored. The rumen metagenomes generated from Next Generation Sequencing (NGS) technology were analyzed using MG-RAST. According to the results stress reaction, resistance and mobilome genes were present in similar amounts in all the three genotypes. Cobalt-zinc-cadmium resistance, fluoroquinolone resistance, methicillin resistance in Staphylococci and multidrug resistance efflux pumps were the most abundant resistant genes and were spread across 20, 24, 16 and 21 bacterial classes, respectively. Bacteria in charge of phage integration and excision, phages replication and phage packaging were mostly allocated to the phyla Firmicutes, Bacteroides and Proteobacteria. Within the stress response genes, metagenomic assembly-based host-tracking analysis identified the extended heat shock dnaK gene cluster as the most abundant genes, while Bacteroides and Clostridium were the principal bacterial hosts. The results show that even with proper use of antimicrobials, the cattle rumen contained an immense distribution of responses to stress, ARGs and mobilome genes distributed in a vast assemblage of hosts. There is also a high correlation between these three functional groups.


INTRODUCTION
Ruminants are a mammalian group that includes domestic cattle, sheep, and goats. The importance of domesticated *Corresponding author. E-mail: mkibegwa@gmail.com Tel: +254723382047.
Author(s) agree that this article remain permanently open access under the terms of the Creative Commons Attribution License 4.0 International License ruminants is derived from their capability to change forages into high-quality, high-protein products for human consumption, through rumen fermentation (Ross et al., 2012). The rumen is the ruminant stomach's first chamber, and it has been bestowed upon by nature to harbor a large amount of symbiotic microflora and fauna . Smallholder dairy cattle are subjected to a range of diets which mainly consist of lignocellulosic by-products such as cereal straws. These plant products are digested by cattle through microbial processes in the rumen (Kong et al., 2010;Morgavi et al., 2013;Patel et al., 2014). The most important microbes in the rumen are bacteria; this is because they are the predominant group. Additionally, the rumen contains an assortment of archaea, fungi and protozoa (Hespell et al., 1997).
Colonization of the rumen by microbial communities begins within the first 24 h of birth (Jami et al., 2013;Li et al., 2012) and continues to increase in mass and density as the cow is exposed to different diets. However, as the rumen microbes increase, there is also a corresponding increase of the resistance, mobile and stress response genes in the cattle digestive tract. Of more importance are the antimicrobial-resistant (AMR) genes and bacteria that the cattle harbor and that have the potential to spread to humans who consume products from these animals. AMR genes and bacteria are mainly acquired from phylogenetically distant microbiomes in the animals environment (soil, food, water, other animals and or humans) (Westphal-Settele et al., 2018). This genetic transfer is facilitated through horizontal gene transfer of the mobilomes (mobile genetic elements) (Forsberg et al., 2012;Wallau et al., 2018). Besides AMR acquisition, AMR genes have been shown to develop in animals. This has been mainly linked to the antimicrobial usage in the production cycle (Cameron and McAllister, 2016;Penders et al., 2013). Nevertheless, AMR has been reported in non-medicated animals (Chambers et al., 2015;Poole, 2012). This appearance of AMR genes has largely been associated with microbial responses to stresses in dietary changes (Auffret et al., 2017). Nutritional modifications have been shown to trigger a "bloom" of particular microbial communities or increase the abundance of other stress-response genes in the gut microbial population (Keto-Timonen et al., 2016;Shin et al., 2015). It is therefore important to understand the diversity of these resistances and stress response genes in smallholder dairy cattle which are subjected to constant variations in diets. Additionally, since the gastrointestinal tract is an open system which comes across numerous bacteria every day, (Baquero, 2012), an understanding of the mobile genes in this system becomes imperative.
For many years, conventional culture-dependent methods have been used to assess rumen microbial population structure. However, these techniques provide a restricted and biased image of the existing rumen microbial community (de Menezes et al., 2011;Hess et al., 2011;Jami and Mizrahi, 2012). These challenges have been overcome with the advent of metagenomics, a next-generation sequencing (NGS) technique that isolates DNA from the whole community irrespective of individual microbial culturing conditions. Using metagenomics, previous studies have assessed functional selections of all genes that confer resistance, including putative or precursor resistance genes (resistome) (Berendonk et al., 2015;Wright, 2007). The resistome has been postulated to be the cause recognized techniques of resistance, such as modifying an antibiotic target, enabling the cell to efflux antibiotic compound or producing an enzyme capable of disabling active compounds (Gomez-Alvarez et al., 2012). Additionally, applications of metagenomics has revealed a complex network of genetic exchange between bacterial pathogens and environmental reservoirs in antibiotic resistance and stress genes studies .
It is therefore imperative that the relationship between the antimicrobial resistance genes, mobile genes and stress response genes be properly investigated especially for animals where proper antimicrobial use is practiced. In addition, breed differences in these genes should be assessed especially at the smallholder farm level. Given the foregoing, the aim of this study was to provide a description of the phylogenetic and functional potential of rumen resistome, mobilome, and stress genes in smallholder dairy cattle reared in the tropics as well as the host bacteria associated with these genes.

Experimental design and rumen sampling
The experimental animals were maintained at the University of Nairobi's Faculty of Veterinary Medicine farm. Six 3-4 years old healthy dairy cattle were used in this study. The animals were from three genotypes; Friesian (Fri), Jersey (Jer) and Friesian X Jersey cross (FriXJer), with two animals per genotype. Animals were purposefully selected based on their breed, body condition, medical history, and stage of lactation. Prior to the study, all animals were maintained on pasture, which was predominantly kikuyu grass (Pennisetum clandestinum) with daily supplementation with dairy meal (a formulated concentrated diet during milking). During the dry seasons, the animals were also offered rhodes grass (Chloris gayana) hay, Napier grass (Pennisetum purpureum), and maize (Zea mays) stover or silage. Additionally, mineral supplement and fresh clean water were given to the cows ad libitum. On the sampling day, individual rumen liquor samples were collected at approximately 0900 h using a flexible stomach tube as previously suggested by Lodge-Ivey et al. (2009). Samples were promptly put on ice, transported to the laboratory and stored at -20°C before metagenome analysis.

DNA extraction, library preparation and Illumina Miseq Sequencing
A homogenized subsample of the rumen liquor samples was

Bioinformatics analysis
Data analyses were performed using the publicly available Metagenome Rapid Annotation Subsystem Technology (MG-RAST) pipeline. The raw sequence information, together with their related quality results (FASTQ format), were used for the optional first quality control (QC) filter to remove duplicate and poor-quality reads. For functional and diversity assessment, the reads that passed the quality filters were subjected to the M5NR database (M5 non-redundant protein database, http://tools.metagenomics.anl.gov/m5nr/) applying an e-value threshold 1.0e-05. The M5NR is a single, searchable, nonredundant database that contains protein sequences and annotations from various sources and their associated tools. In addition, SEED subsystems (Overbeek et al., 2005) (http://www.theseed.org/wiki/Home_of_the_SEED), was used for the functional hierarchical classification applying an e-value threshold of 1.0e-05 (Edwards et al., 2006). The difference in the gene content of resistance to antibiotics and toxic compound, phages and prophages and stress response between different cattle genotypes was quantified using one-way analysis of variance (ANOVA) with Bonferroni adjustment using Genstat version 14 software (Payne, 2011). Significance was established at P ≤ 0.05. To identify which organisms were associated with the genes allocated to each subcategory, sequence alignment using BLAT incorporated in the MG-RAST database was undertaken.

Ethics statement
This study was approved and performed following the University of Nairobi's Faculty of Veterinary Medicine Animal Care and Use Committee (ACUC) guidelines. Animals were handled by experienced animal health professionals to minimize discomfort and injury.

Taxonomic characterization of the sequencing samples
In this study, genetic diversity and functional capacity of the small holder cattle rumen microbiota were characterized through metagenomic sequences. Supplementary Table S1 summarizes the metagenome information. Illumina Miseq sequencing resulted in a total of 4.85 million reads from all the samples, with different reads length (90-94 bp). At the domain level, data was dominated by bacteria followed by eukaryotes, archaea and viruses (Table 1). Within the eukaryotic sequences, Fungi, Metazoa and Viridiplantae were the most abundant. We hypothesize that this could have been as a result of plant DNA contamination. The phylum level breakdown in our data showed that Bacteroidetes, Firmicutes and Proteobacteria predominated in all the genotypes ( Figure 1).

Predicted gene functions of the rumen metagenomes
Sequences annotations using the SEED interface demonstrated the existence of functional protein encoding genes (PEGs). All the metagenomes showed that PEGs belonging to four subsystems namely, carbohydrates (14.47 to 14.87%), protein metabolism (13.7 to 14.38%), clustering-based subsystems (12.68 to 12.81%) and amino acids derivatives (10.17 to 10.57%), were most abundant. The functional classification by the subsystem database also indicated the presence of mobilome, resistome and stress genes as shown in Figure 2. Focusing at specific metabolic pathways, reads assignment to virulence, disease and defence were 1.86, 1.9 and 1.73% in Friesian, FriesianXJersey cross and Jersey, respectively. Genes that coded for Phages, Prophages, Transposable elements and Plasmids accounted for 1.8% in Friesian, 1.93% in FriesianXJersey cross and 1.63% in Jersey. Stress response genes were 1.81, 1.92, 1.63 respectively in the Friesian, FriesianXJersey cross and Jersey cattle genotypes. Oneway ANOVA analysis for the three functional systems, revealed that there were no significant differences within genotypes; virulence, disease and defence (P = 0.803), Phages, Prophages, Transposable elements (P = 0.283) and Stress response (P = 0.355). Pearson correlation analysis of the resistance genes, mobilome and stress genes indicated a high positive correlation between these functional groups (Table 2).

Resistome analyses
MG-RAST classification of virulence, disease and defence identified seven classes. The most abundant cluster among the seven was resistance to antibiotics and toxic compounds (79.05 to 80.70%). About 12.4 to 14.5% of the genes fell in the category of poorly characterized genes associates with resistance called NULL. Virulence, disease and defence associated proteins involved in adhesion to host cells were in the range of 2.59 to 3.57%. Also, 1.22 to 2.61% of genes were associated with host detection. The other classes included genes responsible for ribosomally synthesized bacteriocins, invasion and intracellular resistance and toxins and superantigens (Figure 3). The statistical analysis using One-Way ANOVA showed that only copper homeostasis: copper tolerance was significantly different between all the three genotypes (Table 2). Within the resistance to antibiotics and toxic compounds, the four most abundant subgroups were resistance to fluoroquinolones (28.45 to 29.81%), cobalt-zinc-cadmium resistance (25.75 to 27.12%), multidrug resistance efflux   pumps (14.91 to 16.32%) and methicillin resistance in Staphylococci (6.35 to 8.39%) ( Table 3). MG-RAST -BLAT integration revealed that the distribution of specific RATC gene categories is non-random among bacterial taxa. The fluoroquinolone resistance, cobalt-zinccadmium resistance, multidrug resistance efflux pumps and Methicillin resistance genes were generally scattered across 24, 20, 21 and 14 bacterial classes, respectively. The predominant phylum in all three genotypes was bacteroidetes followed by firmicutes and proteobacteria (Table 4).

Gene assignments to stress response
In the category of stress response genes, heat shock, oxidative stress and poorly characterized genes associates with stress were predominant in all samples. Conversely, cold shock, detoxification and desiccation stress were the least abundant (Figure 4). Amongst the most principal heat shock system, heat shock dnaK gene cluster extended (33.99 -37.77%), rubrerythrin (9.26 -10.51%), regulation of oxidative stress response (8.28 -9.89%), sigmaB stress response regulation (5.65 -7.32%) and oxidative stress (6.24 -6.39%) were predominant in the all treatments samples. No significant differences were observed in all stress response genes between all the three genotypes on One-Way ANOVA (Table 5). Phyla and class-wise affiliation of stress response genes are given in Table 6.

Functional classification of phages, prophages, transposable elements and plasmids
Within the mobilome, phages and prophages, transposable elements and pathogenicity islands had the highest proportions that ranged from; 59.94 to 65.95%, 18.74 to 21.48% and 14.57 to 18.46% respectively. On the other hand, the least amounts were those of gene transfer agent (GTA) (0 to 0.14%), plasmid related functions (0.05 to 0.07%) and poorly characterized functions (0 to 0.16%) ( Figure 5). r1t-like streptococcal phages, phage integration and excision, phage replication and phage packaging machinery were the four major subclasses in phages and prophages. The four subclasses accounted for 74-77% in all the metagenomes. Apart from this subclasses, MG-RAST also identified 13 other subclasses (Table 7). On ANOVA, Phage entry and exit (P = 0.023) and Phage nin genes -N-independent survival (P = 0.031) were significantly different in the three genotypes. In both of this function, the FriesianXJersey crossbreed had significantly higher proportions than the other two genotypes. Bacilli were the chief taxa that contributed to the expression of the genes responsible for r1t-like streptococcal phages and phage packaging machinery while Bacteroidia was more in phage integration and excision and phage replication. Worth noting were the uniquely high proportions of Bacteroidia and Alphaproteobacteria in r1t-like streptococcal phages and phage packaging machinery respectively within the Jersey genotype (Table 8).

DISCUSSION
This research shows that shotgun next generation sequencing can also be used to identify relayed genes from smallholder cattle rumen metagenomes for resistance to antibiotics and toxic compound, phages and prophages and stress response. The shotgun technique outlined for deriving rumen microbiome profiles enables comparison of samples based on the entire population. Sample analysis on MG-RAST identified bacteria was predominant, followed by eukaryotes and viruses as the domains in our samples. This finding was in agreement with a previous study on cattle rumen microbiota (Pitta et al., 2016). At the phylum level, the sequences distribution indicates that the five most predominant clusters belonged to Bacteroidetes, Firmicutes, Proteobacteria, Actinobacteria and Fibrobacteres in all the cattle genotypes. This result was in congruent with findings from an earlier research by Jose et al. (2017a, b) who assesses the rumen microbial and carbohydrate-active enzymes profile in Indian crossbred cattle. Similar to previous studies by Kalyuzhnaya et al. (2008) in cattle rumen and Singh et al. (2012) in buffalo rumen metagenome, the Level 1 subsystem classification showed that the largest percentage of gene fragments allocated to known features in all samples was correlated with clustering subsystems, carbohydrates and protein metabolisms. Genes in the clustering-based subsystems    group are regularly observed together in various species for which particular features are not yet known (Durso et al., 2011). As observed in Figure 2, significant proportions of genes were assigned to the mobilome, resistome and stress genes in the SEED Subsystem Level 2 subcategory. Observation of antimicrobial resistance genes in this study serves to reinforce the theory proposed by Auffret et al. (2017), that the presence of AMR genes is largely linked to microbial stress response to dietary changes, since the animals in the study were reared in a farm with proper use of antimicrobials. Even more important, were the high correlations observed between these three functional groups (mobilome, resistome and stress genes). This result supports the hypothesis suggested by Reddy et al. (2014), which indicates that a complex network exists in the study of antibiotic resistance and stress genes, between bacterial pathogens and environmental reservoirs. Specifically, within the RATC, the most frequently occurring functional groups were fluoroquinolone resistance genes, multidrug resistance efflux pumps and Methicillin resistance in Staphylococci. Similar observations were reported by Durso et al. (2012) in cattle feces microbiome within agricultural and nonagricultural metagenomes and by Cardoso et al. (2012) in snail metagenome. The finding of RATC genes in these cattle that had been reared in a farm with judicious use of antimicrobials further affirmed previous studies that suggested antibiotic resistance genes and antibiotic resistance occurs as an ancient phenomenon in a variety of human-impacted and pristine habitats (Kim and Bae, 2011;Reyes et al., 2010). Seed subsystem composition of phages indicated the predominance of pathogenicity islands that were majorly from phages replication, phage packaging and rlt streptococcal phage genes. This finding indicates that the expression of phage encoding genes in cattle rumen were a reflection of the induction of prophages in rumen bacteria. Genes associated with integrases and islands of pathogenicity were also identified similar to a previous research by Hacker and Kaper (2000). Genes involved in adapting to stress reactions were present in all metagenomes. These comprised of the genes encoding for Chaperone protein (DnaK), Chaperone protein (DnaJ) and nucleoside 5triphosphate RdgB. These genes play an important role in adaptations to psychrophilic lifestyles (Rodrigues and Tiedje, 2008). Contrary to previous studies, in this study, the cattle rumen metagenomes had an enhanced representation of the Sigma B gene, a general stress regulon which induces a wide range of genes in response to various stressful conditions including heat, acid, salt and starvation (Höper et al., 2005). Matches allocated to the genes of more-constituent proteins connected with cold adaptation chaperones DnaK and DnaJ were also of significant interest due to their high abundance in the rumen microbial communities. These genes have been shown to be induced in bacteria when exposed to cold temperatures (Cavicchioli, 2006). We theorize that the presence of these genes could be due to the cold climate in which the animals were reared.
The findings of this study also showed that Clostridia and Bacilli were the predominant bacterial taxa actively playing a part in the resistome, mobilome and stress reactions in the current research. This was in agreement with previous studies on the Indian Buffalo rumen . In both the SEED Subsystem Levels 1 and Percent of genes within each metagenome that are assigned to the listed taxa.
2, the resistance genes, phage and prophage genes and stress responses had statistically no significance, thus giving an indication that the genotype of animals did not affect these microbial functions. This pilot study elucidates the resistome, mobilome and stress responses in Friesian, Jersey and FriesianXJersey crossbreed cows reared in the tropics. Additionally, the study demonstrates that despite the judicious use of antimicrobials, the rumen microbiota in the study animals had a vast assemblage of antimicrobial resistance genes. The study also showed that there exists a correlation between stress genes, antimicrobial resistance genes and mobile genes. Moreover, we found no support for the hypothesis that the resistome, mobilome and stress responses were correlated with the animal breed with minimal differences observed in gene abundance between the three cattle genotypes. However, validating these results using more cows per breed and increasing the sampling repetitions in future studies could undoubtedly augment our understanding of the disparity between individual cows and breeds. In addition, comparing the results with farms where there is indiscriminate use of antimicrobials could present better insights into these microbial profiles.

Authors' contributions
FK, RB, CG, EM, FM designed the research project and helped prepare the manuscript. FK, FM conducted the bioinformatics and statistical analysis. FK, EM, and FM performed the laboratory analyses. All authors read and approved the final manuscript.