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 (Reddy et al., 2014).

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 subjected to DNA isolation using the QIAmp DNA stool mini kit (Qiagen, Valencia, CA). DNA samples were measured on a NanoDrop^TM 2000 spectrophotometer (Thermoscientific, USA) and Qubit® 2.0 Fluorometer (Life Technologies Corporation, Grand Island, NY, USA), to assess DNA quantity as previously described by Habimana et al. (2018). Thereafter, libraries were generated using the Nextera XT DNA Library and index kit. The quality and amount of libraries produced were evaluated using the Agilent Bioanalyzer TapeStation 2200 (Agilent Technologies, Santa Clara, USA) with Agilent High Sensitivity DNA Kit (Agilent Technologies), while the quantification was achieved using Qubit^® fluorometer (Life Technologies). Quality filtered libraries were sequenced on the Illumina Miseq platform at the Biosciences eastern and central Africa - International Livestock Research Institute, (BecA ILRI) Hub laboratory in Nairobi.

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, non-redundant 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. One-way 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-zinc-cadmium 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).

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 non-agricultural 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 5-triphosphate 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 (Reddy et al., 2014). In both the SEED  Subsystem  Levels 1 and 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.

The datasets generated and/or analysed during the current study are available in the https://www.mg-rast.org/linkin.cgi?project=mgp79210, using the following MG-RAST IDs; mgm4736079.3, mgm4752069.3, mgm4748083.3, mgm4752189.3, mgm4751799.3, mgm4735643.3

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.

The authors have not declared any conflict of interests.

The BecA-ILRI Hub supported this work through the Africa Biosciences Challenge Fund (ABCF) program. The ABCF program is funded by the Australian Department for Foreign Affairs and Trade (DFAT) through the BecA-CSIRO partnership; the Syngenta Foundation for Sustainable Agriculture (SFSA); the Bill & Melinda Gates Foundation (BMGF); the UK Department for International Development (DFID) and the Swedish International Development Cooperation Agency (Sida).