Bacterial diversity and the antimicrobial resistome in the southwestern highlands of Saudi Arabia.

Soil is a reservoir of microbial diversity and the most supportive habitat for acquiring and transmitting antimicrobial resistance. Resistance transfer usually occurs from animal to soil and vice versa, and it may ultimately appear in clinical pathogens. In this study, the southwestern highlands of Saudi Arabia were studied to assess the bacterial diversity and antimicrobial resistance that could be affected by the continuous development of tourism in the region. Such effects could have a long-lasting impact on the local environment and community. Culture-dependent, quantitative polymerase chain reaction (qPCR), and shotgun sequencing-based metagenomic approaches were used to evaluate the diversity, functional capabilities, and antimicrobial resistance of bacteria isolated from collected soil samples. Bacterial communities in the southwestern highlands were mainly composed of Proteobacteria, Bacteroidetes, and Actinobacteria. A total of 102 antimicrobial resistance genes (ARGs) and variants were identified in the soil microbiota and were mainly associated with multidrug resistance, followed by macrolide, tetracycline, glycopeptide, bacitracin, and beta-lactam antibiotic resistance. The mechanisms of resistance included efflux, antibiotic target alteration, and antibiotic inactivation. qPCR confirmed the detection of 18 clinically important ARGs. In addition, half of the 49 identified isolates were phenotypically resistant to at least one of the 15 antibiotics tested. Overall, ARGs and indicator genes of anthropogenic activities (human-mitochondrial [hmt] gene and integron-integrase [int1]) were found in relatively lower abundance. Along with a high diversity of bacterial communities, variation was observed in the relative abundance of bacterial taxa among sampling sites in the southwestern highlands of Saudi Arabia.

Introduction

Saudi Arabia represents approximately 80% of the Arabian Peninsula (Almazroui et al., 2012) and is mainly composed of desert. Unlike the majority of the country, which receives little precipitation and experiences high temperatures, the southwestern highlands are famous for their diverse vegetation (El-Juhany and Aref, 2013). Moreover, as the average rise in global temperature increases, Saudi Arabia is relatively more vulnerable than other countries to changes in climate and higher temperatures, particularly in the western highlands owing to their comparatively high moisture (Almazroui, 2020). Overall, the impacts of global warming on microbial community structure and subsequently the microbe-assisted decomposition processes are higher in alpine, arctic, and highland regions than in other regions (Bardgett et al., 2008, Li et al., 2016, Frindte et al., 2019, Smith et al., 2019). Our previous study demonstrated that changes in soil microbial communities vary along the elevation gradients of the southwestern highlands of Saudi Arabia (Yasir et al., 2015). The dominant phyla found in that study included Proteobacteria, Actinobacteria, and Acidobacteria (Yasir et al., 2015). So far, the functional capabilities and antimicrobial resistome of the microbial communities from the southwestern highlands of Saudi Arabia have not been investigated.

Soil is considered a hotspot for antimicrobial-producing bacteria (Forsberg et al., 2014). Bacteria that produce antibiotics develop resistance mechanisms to evade the effect of self-antibiotics (Martínez, 2008, Forsberg et al., 2014). In addition, neighboring bacteria also undergo adaptation processes and also develop resistance (Martínez, 2008, Forsberg et al., 2014). With a declining trend in the discovery of novel antibiotics, researchers have focused on exploring the natural antibiotic resistome (Martínez, 2008). These efforts will lead to understanding the ecology and evolution of antibiotic resistance in nonclinical settings (Martínez, 2008, Nesme and Simonet, 2015). Such investigations of reservoirs of both known and novel antimicrobial resistance mechanisms could subsequently add to antibiotic resistance management (Raymond, 2019). A study from the Netherlands compared pre-antibiotic and post-antibiotic soil samples and illustrated a relatively greater abundance of resistance genes in the post-antibiotic samples (Knapp et al., 2010). Moreover, antibiotic-resistant bacterial strains have been identified in entirely isolated caves that are known to completely lack any human antibiotic use (Bhullar et al., 2012).

In this study, we used a shotgun metagenome-based approach to demonstrate the functional capabilities and antimicrobial resistance of apparently pristine environments in the southwestern highlands of Saudi Arabia. We then used a culture-based approach to evaluate the antimicrobial resistance in bacterial isolates. The study was intended to reveal the structure of the antimicrobial resistome in highlands with limited anthropogenic activities documented based on visual observation, and was quantified by indicator genes of anthropogenic activities. Exploring the functional capabilities of the resident microbes will help uncover the interaction between climate change and the ecology of resident microorganisms.

Materials and methods

Sampling

The southwestern highlands constitute a geographically unique region in Saudi Arabia that is covered by natural vegetation and is well known as a plant biodiversity hotspot on the Arabian Peninsula (El-Juhany and Aref, 2013). Annual rainfall in this region is in the range of 300–500 mm, and summer temperature rise to 30 °C and above, which is comparatively higher than other highlands around the world (El-Juhany and Aref, 2013). Soil samples were collected from three sites in this study. The SWH6 samples were collected from Sabt Al Alayah in the vicinity of Shibanah Wild Park. Plants at the site included Olea europaea, Tarchonanthus camphorates, and bushes (Table S1). The SWH8 and SWH11 samples were collected from mountainous locations with scattered vegetation (Table S1). All the sites appeared to have limited anthropogenic activities. At each sampling site, quadrates of 3.5 × 3.5 m were selected (Yasir et al., 2015), and the top 3 cm of surface soil was removed. Five replicate samples were collected with a sterilized spatula from the four corners and the middle of each quadrate, at a depth of 10 cm, and the samples were sieved through 2-mm mesh. Samples from each spot were mixed in a sterilized bag to obtain a homogenized sample and frozen at − 80 °C for further processing.

Shotgun sequencing and data analysis

Total genomic DNA was extracted from the 300 mg of each homogenized soil sample using the PowerSoil DNA Isolation Kit (MO BIO, USA) following the manufacture protocol. DNA concentration was measured with Qubit dsDNA BR assay using Qubit fluorometer (Invitrogen, USA). The DNA libraries of ∼ 400 bp insert size were prepared using Nextera XT DNA library preparation kit (Illumina, Inc. USA). Sequencing was performed with 2 × 300 bp chemistry on a MiSeq platform (Illumina, Inc.) using a V3 kit (Illumina, Inc.). The raw sequence reads were uploaded to Metagenomic Rapid Annotations using Subsystems Technology (MG-RAST) and were analyzed using the default settings (Meyer et al., 2008). The classification of bacterial communities at different taxonomic levels from samples metagenomes was performed by annotating the ORFs against the NCBI reference sequence (RefSeq) database in MG-RAST. For functional analysis, the ORFs were mapped to Clusters of Orthologous Group of Proteins (COGs) and KEGG Orthology (KO) databases. Antimicrobial resistance genes (ARGs) were identified from the metagenomes using an artificial intelligence-based DeepARG tool with default settings (Arango-Argoty et al., 2018).

Quantitative polymerase chain reaction of antimicrobial resistance genes

Quantification of 22 clinically important ARGs (Table S2) was performed from the total genomic DNA of the southwestern highland samples with qPCR using an Applied Biosystems® 7500 instrument (Thermo Fisher Scientific, USA) as described previously (Tan et al., 2018, Ullah et al., 2019). Briefly, a no-template negative control was included for each primer set. Human mitochondrial (hmt) and class 1 integron-integrase (int1) genes were used as indicators of anthropogenic activities with each studied sample (Tan et al., 2018). Bacterial 16S rRNA gene quantification was performed to normalize the resistance gene copy numbers per sample to the total bacterial community in each sample (Tan et al., 2018). To prepare the standard curve, the tetB gene was cloned in Escherichia coli DH5α using pGEM®-T Easy Vector Systems (Promega, USA). The recombinant plasmid was extracted and quantified. The quantity of DNA was calculated for 101 to 104 copies of the gene in g/µl of DNA, and dilutions were prepared accordingly (Tan et al., 2018). qPCR was performed in triplicate for each copy number dilution.

Bacterial culturing and antimicrobial susceptibility screening

We used low-nutrient medium for culture-dependent analysis. The medium was composed of R3A broth at 1:3 dilution and actinomyces broth at 1:3 dilution supplemented with 1.5% agar following the procedure described previously (Yasir, 2018). Five grams of each sediment sample was resuspended in 45 ml sterilized normal saline. The homogenous suspension of each sample was serially diluted (101 to 106), and a 100-μl aliquot from each dilution was spread on media agar plates (Yasir et al., 2019). The plates were incubated at 30 °C for 72 h, and colony-forming units (CFUs) were counted. The colonies were purified by sub-culturing and stored in 15% glycerol at − 80 °C (Farman et al., 2019).

The purified colonies were identified by matrix-assisted laser desorption ionization-time of flight mass spectrometry using the VITEK MS (bioMérieux, France) system described previously (Farman et al., 2019). The calibration was performed using E. coli ATCC 25,922 strain to validate the run. The unidentified strains by VITEK MS were processed for 16S rRNA gene sequencing as described previously (Yasir et al., 2009) using universal primers 27F and 1492R. The isolates were identified from 16S rRNA gene sequence identity with closely related type strains from the blast at EzBioCloud Database (https://www.ezbiocloud.net/). The cut-off value was set to ≥ 99% for bacterial isolate identification. Antibiotic sensitivity susceptibility testing was accomplished using Kirby - Bauer disk diffusion method (Hudzicki, 2009). The antibiotics were used in the following concentration; carbenicillin (100 mcg), amoxicillin/clavulanic acid (30 mcg), imipenem (10 mcg), meropenem (10 mcg), cefotaxime (30 mcg), cefepime (30 mcg), streptomycin (10 mcg), gentamycin (10 mcg), kanamycin (30 mcg), amikacin (30 mcg), azithromycin (15 mcg), trimethoprim/sulfamethoxazole (25 mcg), tetracycline (30 mcg), chloramphenicol (30 mcg), and ciprofloxacin (5 mcg). The results were interpreted according to the French Society of Microbiology standard for Gram-positive and Gram-negative bacteria (SFM Antibiogram Committee, 2003).

Statistical analysis

Statistical analysis of the qPCR results was performed using the nonparametric equivalent of ANOVA. p-Values ≤ 0.05 were considered significant. Statistical analysis was performed using SPSS 22 software (IBM, USA), and the heat maps were constructed with an online GENE-E tool (https://software.broadinstitute.org/morpheus/).

Results

Metagenomic analysis of bacterial diversity

A total of 104,336 sequence reads containing 24,356,528 bp were obtained from SWH6; 287,202 sequence reads containing 65,017,973 bp were obtained from SWH8; and 406,376 sequence reads containing 94,861,195 bp were obtained from SWH11. Dereplication identified 157 sequences as artificial duplicate reads in SWH6, 532 sequences in SWH8, and 644 in SWH11. Of the sequences analyzed, 12.7% in SWH6, 6.7% in SWH8, and 14.2 % in SWH11 failed to pass the quality control of the MG-RAST pipeline. The mean length of the reads in the three metagenomes was 231 ± 67 bp. The mean G + C content was 59 ± 10% in SWH6, 63 ± 9% in SWH8, and 64 ± 8% in SWH11. The total number of predicted protein features in SWH6 was 90,280; 262,775 in SWH8; and 344,573 in SHW11. The predicted rRNA features were 303 in SWH6, 897 in SWH8, and 949 in SWH11. The metagenomic features of the southwestern highland are summarized in Table S3.

The bacterial diversity of the southwestern highlands sites was determined from shotgun sequencing. A total of 27 phyla, including a candidate phylum of Candidatus Poribacteria, and a group of unclassified bacteria were detected in the samples from the studied sites. Nine phyla were found at > 1% relative abundance in each sample, and Proteobacteria was the dominant phylum, composing 43% to 57% relative abundance of the total bacterial community structure (Fig. 1A). Other major phyla included Actinobacteria, Bacteroidetes, Firmicutes, and Verrucomicrobia (Fig. 1A). Except for the rare phylum Dictyoglomi, which was absent from SWH6, all 26 phyla were commonly detected in the analyzed samples (Fig. 1B). At the family level, 204 different families were commonly detected among all samples, but variation was noticed in the relative abundance percentages (Fig. 1B). The five predominant families in SWH6 included Cytophagaceae, Bradyrhizobiaceae, Burkholderiaceae, Flavobacteriaceae, and Comamonadaceae (Fig. 1A). In SWH8, the five most abundant families included Planctomycetaceae, Bradyrhizobiaceae, Streptomycetaceae, Micromonosporaceae, and Mycobacteriaceae (Fig. 1A). The top five abundant families for SWH11 included Bradyrhizobiaceae, Planctomycetaceae, Burkholderiaceae, Solibacteraceae, and Streptomycetaceae (Fig. 1A).

Fig. 1

Bacterial communities’ analysis in the southwestern highland of Saudi Arabia. (A) Percent relative abundance of dominant phyla, families, and genera identified from shotgun sequencing. (B) Venn diagrams analysis for the common and unique phyla, families, and genera among samples. SWH, southwestern highlands.

Diversity analysis at the lower taxonomic level of genus revealed that variation exists among the samples collected at different sites. Among the 532 identified genera, 459 were common among the samples, whereas 13 genera were unique to SWH11, and three genera were unique to each SWH6 and SWH11 (Fig. 1B). A relatively increased number of genera were detected in SWH11 (5 2 6) followed by SWH8 (5 1 1). In SWH6, 467 distinct genera were detected (Fig. 1B). An average number of 16 ± 3 genera were found at ≥ 1% relative abundance in each sample. The genera of Bradyrhizobium (2.2% ± 0.4%), Burkholderia (1.9% ± 0.3%), and Rhodopseudomonas (1.8% ± 0.4%) were commonly found at relatively higher abundance in the sampling sites (Fig. 1A). Among other most abundant genera in SWH6 sample include Pseudomonas (2.1%), Cytophaga (1.6%), and Marivirga (1.5%) (Fig. 1A). Abundant genera in SWH8 included Streptomyces (4.6%), Mycobacterium (2.8%), and Methylobacterium (2.7%) (Fig. 1A). Other abundant genera in SWH11 samples included Streptomyces (3.0%) and Mycobacterium (1.8%) (Fig. 1A). Moreover, we analyzed rare genera detected at < 0.01% relative abundance that generally have a minor representation in biomass but have a major role in microbial diversity and play an important part in biogeochemical cycles. We found 123 rare genera in the western highland metagenomes, including 40 genera from SWH6, 77 from SWH8, and 78 from SWH11. Only 10 genera were common among the samples found at ≤ 0.01% relative abundance in each sampling site.

Metagenomics insights to antimicrobial resistance genes

Metagenome data analysis revealed that various ARGs were present among the southwestern highlands samples. A total of 102 different resistance genes or variants were detected in these samples, and they were previously identified as conferring resistance to 19 classes of antibiotics and included multidrug resistance genes and a group of unclassified genes not assigned to any specific antibiotic class. A relatively increased abundance of ARGs (45.9% ± 3.4%) was found in the southwestern highlands samples producing multidrug resistance. Among antibiotic classes, a relative increase was noticed in resistance genes against macrolide-lincosamide-streptogramin (MLS) (14.7% ± 0.8%), tetracycline (7.6% ± 1.4%), glycopeptide (4.2% ± 1.7%), bacitracin (4.5% ± 0.4%), and beta-lactam (3.9% ± 1.3%) antibiotics (Fig. 2A).

Fig. 2

Resistance against antibiotic classes and antimicrobial resistance genes (ARGs). (A) Percent relative abundance of the antibiotics classes related to identified ARGs retrieved from the metagenomes of southwestern highlands samples. (B) Analysis of common and unique ARGs among the samples. AA, aminoglycoside-aminocoumarin; MLS, macrolide-lincosamide-streptogramin.

Among the 102 ARGs and variants retrieved, 18 were common among the southwestern highland metagenomes (Fig. 2B). Increased diversity of resistance genes was noticed in the SWH11 (n = 80), followed by SWH8 (n = 66 genes; Fig. 2B). In SWH6, a relatively low number of ARGs were retrieved (n = 28). Similarly, the SWH11 sample was found to contain genes for resistance against most of the antibiotic classes (16 antibiotic classes), followed by SWH8 (13 antibiotic classes) and SWH6 (10 antibiotic classes). The abundance ratio of individual ARGs ranged from 0.01 to 0.9 counts/16S rRNA gene in the studied samples (Fig. 3). The ARGs causing multidrug resistance (rpoB2) and resistance against MLS (macB), bacitracin (bcrA, bacA), glycopeptide (vanR), and tetracycline [tetA(48)] were commonly found in the metagenomes at a relatively high range of 0.1 to 0.9 counts/16S rRNA gene (Fig. 3, Fig. S1). The identified ARGs produced antimicrobial resistance mainly through six types of mechanisms. An average increased number of ARGs in western highlands metagenomes were associated with antibiotic efflux (59/102; 35.7 ± 14.9) mainly causing multidrug resistance, followed by antibiotic target alteration (13/102; 8.3 ± 2.9) and antibiotic inactivation (13/102; 6.3 ± 3.8).

Fig. 3

Antimicrobial resistance genes (ARGs) retrieved from metagenomes of southwestern highlands. ARG counts were normalized using 16S rRNA genes. ARG abundance is expressed as copy numbers of ARG per copy of 16S rRNA gene. AE, antibiotic efflux; ATA, antibiotic target alteration; AA, aminoglycoside-aminocoumarin; MLS, macrolide-lincosamide-streptogramin.

qPCR analysis of antimicrobial resistance genes

In the next step, 22 clinically prevalent ARGs and two indicator genes (hmt, int1) were quantified from the total genomic DNA by qPCR. Four genes (qnrS, blaDHA-1, blaGES-1, and tetM) were not detected across the samples, but most of the commonly found ARGs were detected among all samples causing resistance to beta-lactam, aminoglycosides, quinolones, tetracycline, and sulfonamides antibiotics (Fig. 4A). The number of ARGs varied among the samples. Genes encoding for beta-lactam resistance were significantly abundant in the SWH8 sample (blaSHV-1, blaCMY-2, and blaVIM-1) followed by SWH6 (blaSHV-1 and blaAmpC) and SW11 (blaOXA-1). Aminoglycoside resistance gene accC4 in SWH6 (3.6 × 10−3/16S rRNA gene) was ∼ 2-fold higher than in SWH11 (1.9 × 10−5/16S rRNA gene) and was not detected in SWH8 (Fig. 4B). Other aminoglycoside-resistance genes (aac2 and aacC3) were commonly detected in the southwestern highlands metagenomes. Quinolones resistance genes were detected at relatively increased quantities in SWH8 (qnrB, qnrD, and qepA) followed by the SWH11 sample (qnrB, and qnrD). The qnrB and qnrD genes were not detected in the SWH6 sample (Fig. 4B). The tetracycline-resistance genes tetD and tetC were commonly detected in the range of 9.5 × 10−4 to 6.5 × 10−3 and 9.2 × 10−5 to 2.1 × 10−4/16S rRNA gene, respectively, in the studied samples. A relatively low level of tetracycline resistance was observed in the SWH6 sample (Fig. 4B). We analyzed the genes int1 and hmt as proxies for anthropogenic pollution and found values ranging from 1.9 × 10−5 to 2 × 10−4/16S rRNA gene and 2.1 × 10−6 to 1.0 × 10−5/g of soil, respectively. The SWH6 sample had around one-fold higher int1 abundance than the SWH8 and SWH11 sample. The hmt gene quantification revealed around one-fold higher anthropogenic activities in the SWH11 sample than the SWH6 and SWH8 samples (Fig. 4B).

Fig. 4

Quantification of antimicrobial resistance genes (ARGs) in southwestern highlands. (A) Abundance of the resistance genes grouped according to the antibiotics class they conferred resistance. (B) Quantification of ARGs in the respective samples. ARGs and int1gene copy numbers were normalized using 16S rRNA gene copy numbers. White boxes indicate respective gene was not detected. Ab class, antibiotics class.

Phenotypic analysis of antimicrobial resistance

Antimicrobial resistance in the southwestern highlands bacteriome was confirmed by susceptibility analysis of isolates that tested response to 15 antibiotics. Culture-dependent analysis retrieved from 2.2 × 105 to 4.0 × 106 colony-forming units (CFUs) among the samples from the selected sites (Table S1). Considering the morphology, 49 bacteria were identified and classified into 23 distinct bacterial species, mainly from Firmicutes, Proteobacteria, and Actinobacteria. The southwestern highlands samples mainly contained species from Bacillus and Pseudomonas, such as Bacillus simplex, Bacillus megaterium, Pseudomonas fluorescens, and Pseudomonas protegens. Clinically important pathogenic/opportunistic bacteria found in these samples included Klebsiella pneumoniae, Micrococcus luteus, Staphylococcus haemolyticus, Staphylococcus hominis, and Stenotrophomonas maltophilia.

Antimicrobial susceptibility profiling revealed that 49% of the cultured bacteria showed resistance to at least one of the antibiotics tested (Fig. 5). Among gram-positive isolates, five were resistant to cefotaxime from third-generation cephalosporin and azithromycin. Resistance to meropenem was identified in Bacillus simplex and Microbacterium flavescens (Fig. 5). No resistance was observed in gram-positive isolates against streptomycin, gentamycin, and amoxicillin/clavulanic acid. Gram-negative isolates showed a relatively low level of resistance against tested antibiotics (Fig. 5). The K. pneumoniae isolate swh22 was resistant to cefotaxime and amoxicillin/clavulanic acid. Resistance to meropenem was observed in P. protegens isolate swh100.

Fig. 5

Antimicrobial susceptibility of the cultured isolates from southwestern highland of Saudi Arabia against the tested antibiotics. The Y-axis represent the number of isolates. Amox/Clav, amoxicillin/clavulanic acid; Trim/Sulfa, trimethoprim/sulfamethoxazole.

Metagenome functional analysis

Metagenomes from the southwestern highlands were analyzed using Clusters of Orthologous Groups (COG) mapping and KEGG Orthology (KO) mapping. Functional analysis revealed that a major proportion of the metagenome coded for housekeeping functions. According to COG mapping, the open reading frames (ORFs) were mainly associated with metabolism, cellular processes and signaling, and information storage and processing (Fig. S2A). Functional annotation for metabolism was the most dominant (51.8% ± 1%) of the total mapped ORFs in the analyzed samples (Fig. S2A). Of the 1714 annotated functions, 239 were unique to the SWH11 sample, 118 were unique to the SWH8 sample, and 48 were unique to the SWH6 sample (Fig. S2B). All three sample sites shared a total of 846 functions. Additionally, the analysis using the KO database was performed, reflecting relatively similar results. ORFs annotated for metabolism-associated functions were the most predominant (60.5% ±0.1%), followed by genetic information processing (17.9% ± 0.5%) and environmental information processing (15% ± 0.5%) (Fig. S3A).

In total, 149 pathways were found; nine of those were unique in SWH11, five in SWH8, and three in SWH6. Whereas, 114 detected pathways were shared among all three sample types (Fig. S3B). The top KEGG pathways at sub-level (level-2), detected with at least 5% relative abundance, included amino acid metabolism, carbohydrate metabolism, translation, energy metabolism, signal transduction, and metabolism of cofactors and vitamins (Fig. 6A). Further analysis at the subcategory level demonstrated that the predominant pathways (relative abundance range 2% to 8%) were related to genes encoding ABC transporters [PATH:ko02010]; aminoacyl-tRNA biosynthesis [PATH:ko00970]; two-component system [PATH:ko02020]; alanine, aspartate, and glutamate metabolism [PATH:ko00250]; oxidative phosphorylation [PATH:ko00190]; glycine, serine, and threonine metabolism [PATH:ko00260]; valine, leucine, and isoleucine degradation [PATH:ko00280]; purine metabolism [PATH:ko00230]; RNA degradation [PATH:ko03018]; and cysteine and methionine metabolism [PATH:ko00270] (Fig. 6B). All the genes associated with carbohydrates were categorized into 15 KEGG pathways among the three sampling sites (Fig. 6B), with the relative abundance of genes associated with glycolysis/gluconeogenesis [PATH:ko00010] (2.1% ± 0.2%), pyruvate metabolism [PATH:ko00620] (1.8% ± 0.3%), pentose phosphate pathway [PATH:ko00030] (1.4% ± 0.1%), and starch and sucrose metabolism [PATH:ko00500] (1.4% ± 0.1%).

Fig. 6

Percent relative abundance of KEGG Orthology pathways of southwestern highlands metagenomes. (A). Dominant pathways identified at sub-level 2 of relative abundance of ≥ 1%. (B) Functional pathways identified at sub-level 3 of relative abundance of ≥ 1%.

Discussion

The bacterial community of soil is vital for life on Earth, and it plays a significant role in cycling carbon, nitrogen, and organic matter (Jansson and Hofmockel, 2018). Soil is also a supportive environment for the evolution and development of antimicrobial resistance available for exchange with clinical pathogens (Forsberg et al., 2012). Previous studies have identified the exchange of ARGs between environmental bacteria from soil and clinical pathogens and suggest that soil bacteria could be a source of some ARGs identified in health care facilities (Forsberg et al., 2012, Chatterjee et al., 2018). Overall, antimicrobial resistance has been thoroughly investigated in clinical samples and health care facilities from Saudi Arabia and worldwide (Farman et al., 2019, Shah et al., 2019). In contrast, limited studies about the resistome have been published from the natural environment of the Arabian Peninsula. Understanding and measuring the drivers, potential sources, and reservoirs of antibiotic resistance in various environments are of immense importance (Forsberg et al., 2012, Forsberg et al., 2014). This study investigated the bacterial community and associated antimicrobial resistance in the southwestern highlands of Saudi Arabia. This region occupies 4.3% of the total country area and is considered relatively pristine because of limited human activities. However, in the last few years, tourism has increased in this region, which may adversely affect the local environment.

Bacterial diversity from shotgun sequencing revealed Proteobacteria was the predominant phylum among samples from study sites and was relatively abundant in SWH6 and SWH11 sites. Consistent with our findings, Proteobacteria has been reported to be abundant in highlands from other geographical regions and common in different soil microbial communities (Yasir et al., 2015, Menoyo et al., 2017, Yasir, 2018). Zhang et al. (2013) observed that different bacterial phyla respond differently to climate change and found Acidobacteria and Gammaproteobacteria to be more sensitive to climate effects. Comparative analysis of phyla and lower-level taxonomic classifications revealed variation in bacterial community composition and relative abundance of taxa among sampling sites in the southwestern highlands, which was consistent with our previous observations (Yasir et al., 2015). These variations could result from differences in the soil texture, temperature, oxygenation, nutrients, and vegetation between the sites (Zhang et al., 2017). Additionally, we found several pathogens/opportunistic pathogenic bacteria, including K. pneumoniae, M. luteus, S. haemolyticus, S. hominis, and S. maltophilia, in the analyzed sample sites. Our results are in parallel with previous studies, which showed that potentially pathogenic bacteria reside in various environments of Saudi Arabia (Ansari et al., 2015, Ullah et al., 2019). The source of these bacteria may have been animals that were carriers and reservoirs of different pathogenic bacteria (Manyi-Loh et al., 2016).

Resistance to antibiotics is of great concern for both public health and the ecosystem. Environmental microorganisms contribute to the dissemination and spread of such resistance (Forsberg et al., 2012, Forsberg et al., 2014). Previous studies have demonstrated an array of antibiotic resistance potential in soil bacteria (Martínez, 2008, Forsberg et al., 2014). In this study, we identified 49 bacterial isolates from the southwestern highlands of Saudi Arabia based on a culture-dependent approach. Approximately half of the isolates showed resistance to antibiotics that are commonly used in clinical settings. Moreover, multidrug resistance was observed in isolates of B. simplex and M. flavescens against ciprofloxacin, amikacin, cefotaxime, azithromycin, meropenem, and cefepime. Similarly, antimicrobial resistance has been reported along an elevation gradient from Taibai Mountain in China (Peng et al., 2018). In the literature, multidrug resistance has generally been observed in species of Bacillus and Microbacterium but not from B. simplex and M. flavescens, and this finding might be unique to the southwestern highlands of Saudi Arabia (Kim et al., 2009).

Significant studies of the environmental resistome have reported that terrestrial bacteria sourced from anthropogenic activities and indigenous resident microbial flora of soil are the primary sources of ARGs (Martínez, 2008, Knapp et al., 2010, Forsberg et al., 2014). ARGs can be quickly disseminated from commensal microorganisms via various mechanisms, including horizontal gene transfer (Forsberg et al., 2012). A dire need exists to identify the antibiotic resistome and its sources and evolution to mitigate its dissemination from unmanaged ecosystems to the human milieu (Forsberg et al., 2012, Nesme and Simonet, 2015). In this study, we found 102 ARGs and variants in the southwestern highlands of Saudi Arabia. Clinically important ARGs encoding resistance to various antibiotics, including beta-lactam, aminoglycoside, quinolone, tetracycline, and sulfonamide, were confirmed by qPCR. Consistent with previous findings, most ARGs commonly found in environmental bacteria were associated with efflux, which underlies multidrug resistance. However, genes like blaOXA, aacC3, and aacC4 were sporadically found in the samples. Moreover, detecting the high diversity of ARGs and a relatively low concentration of indicator genes (hmt, int1), mainly in the SWH8 sample, suggests that anthropogenic activities were not the only source for the observed ARGs found in the southwestern highlands. However, limited studies have been published regarding the extent of ARGs diversity from highlands (Peng et al., 2018). A discrepancy exists in whether any associations exist between ARGs abundance in such environments and the extent of anthropogenic activities. However, studies from other pristine environments, notably from Antarctica (Van Goethem et al., 2018), have reported detection of resistance genes.

Soil microbial community structure is a highly diverse and integral part of ecosystems that plays a crucial role in organic matter decomposition and nutrient recycling, thus ultimately affected climate change (Davidson and Janssens, 2006). Functional analysis of the southwestern highlands metagenomes revealed that most ORFs were associated with the housekeeping genes and mainly related to amino acid, carbohydrate, and energy metabolism. Soil is the primary carbon reservoir (Heimann and Reichstein, 2008). Any considerable variations in its cycling can significantly affect the levels of atmospheric carbon dioxide. Overall, the functional capabilities of soil microbiota and their interactions with other diverse factors such as vegetation and animals significantly contribute to regulating climate change (Bardgett et al., 2008, Heimann and Reichstein, 2008).

Conclusions

This study explored the bacterial community composition, functional capabilities, and antibiotic resistome in soil of the southwestern highlands of Saudi Arabia. A significant portion of the metagenomes was associated with carbohydrate metabolism, energy metabolism, and environmental information processing that the local environment could affect. The ARGs causing resistance to antibiotics used in health care facilities were detected at relatively low levels compared with ARGs associated with efflux and commonly found in environmental bacteria and pathogenic species. However, further studies are required to investigate the potential risk of disseminating ARGs from southwestern highlands metagenomes to animals or humans that may ultimately compromise the effectiveness of clinically important antibiotics.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.