The gut mycobiota of rural and urban individuals is shaped by geography.

BACKGROUND: Understanding the structure and drivers of gut microbiota remains a major ecological endeavour. Recent studies have shown that several factors including diet, lifestyle and geography may substantially shape the human gut microbiota. However, most of these studies have focused on the more abundant bacterial component and comparatively less is known regarding fungi in the human gut. This knowledge deficit is especially true for rural and urban African populations. Therefore, we assessed the structure and drivers of rural and urban gut mycobiota.
RESULTS: Our participants (n = 100) were balanced by geography and sex. The mycobiota of these geographically separated cohorts was characterized using amplicon analysis of the Internal Transcribed Spacer (ITS) gene. We further assessed biomarker species specific to rural and urban cohorts. In addition to phyla which have been shown to be ubiquitous constituents of gut microbiota, Pichia were key constituents of the mycobiota. We found that geographic location was a major driver of gut mycobiota. Other factors such as smoking where also determined gut mycobiota albeit to a lower extent, as explained by the small proportion of total variation. Linear discriminant and the linear discriminant analysis effect size analysis revealed several distinct urban and rural biomarkers.
CONCLUSIONS: Together, our analysis reveals distinct community structure in urban and rural South African individuals. Geography was shown to be a key driver of rural and urban gut mycobiota.

Background

By comparison to prokaryotes (bacteria and archaea), eukaryotes are considered part of the rare “biosphere” of the gut [1, 2]. Despite their low abundances, fungi play significant roles in host physiology [2-5]. Recent studies have shown that the gut fungal community composition is less stable over time, compared to bacterial communities [4, 6, 7]. These studies suggest that the gut mycobiota is more variable than bacterial communities, and may be influenced substantially by environmental factors [3, 7]. Despite evidence confirming the gut microbiota is diverse and interacts with the host immune system [3, 8, 9], knowledge regarding the community structure of healthy human gut mycobiota remains scant.

Most studies have focused on the potential roles played by the mycobiota in the aetiology of gut diseases [10-12]. These studies have provided crucial insights on the role of the mycobiota as a potential driver of immunological disorders and as opportunistic pathogens in immunocompromised hosts [13]. Further, dysbiosis of gut mycobiota has been linked to obesity, colorectal cancer and Inflammatory Bowel Diseases (IBDs) [12, 14, 15]. Decreased abundances of Saccharomyces cerevisiae and higher proportions of Candida albicans were found in IBD patients compared to healthy controls. A recent study showed that Crohn’s disease-specific gut environments may select for fungi to the detriment of bacteria suggesting disease-specific inter-kingdom network alterations in IBD [12]. Yet, despite these effects, there remains a clear deficit in knowledge regarding the precise role played by the gut mycobiota in disease prevention. Relatedly, the factors which drive the diversity and community structure of gut mycobiome remain underexplored. Assessing the influence of environmental factors on the gut mycobiome across a wider group of participants is crucial for determining the effects on host-microbiota dynamics and health.

Several studies have evaluated the effects of age [16-18], gender [17], diet [19], diabetes and obesity [15, 20, 21], anorexia nervosa [22], differences across body sites [23, 24] and geographical locations [6, 25, 26] on mycobiome composition and diversity. Yet, these studies are mostly disease centric or focussed on Asian [26] and/or Western populations [6, 7, 19]. To our knowledge, only one study has investigated the gut eukaryotic diversity of African individuals [27]. Although these studies improved our understanding of the mycobiome, there may be several confounding factors such as genetic differences. These differences make it difficult to assess, for instance, the effects of living in urban or rural areas on the microbiome. The effects of diet, geographic locality and lifestyle, on the gut microbiome are often assumed but rarely examined. Where these relationships are assessed, the majority of studies have primarily focused on the ecologically abundant bacteria [28, 29] with assertions that their patterns will likely hold for other taxa, including mycobiomes.

Here, we applied amplicon sequencing of the fungal internal transcribed spacer (ITS) of the rRNA genes on samples collected from individuals living in urban and rural areas in Africa. We provide the first insights regarding the drivers of mycobiome community structure and potential biomarkers specific to individuals from urban and rural locations. Previous studies of the gut mycobiome have primarily investigated small groups with fewer than 20 individuals [25, 30, 31] with very few studies investigating larger groups [6, 7]. This study represents the first analysis of the faecal mycobiota in a large group of healthy sub-Saharan individuals (100 volunteers). Furthermore, this is the first study which compares the composition and diversity of the gut mycobiome of geographically separated non-western individuals with the same ethnicity. We further explored potential biomarker taxa in urban and rural individuals and explore how these taxa vary between the two areas. Using extensive predictor variables collected from participants, we show that geography and lifestyle structure the gut mycobiome of rural and urban South African individuals.

Results

Similarities and differences between urban and rural individuals

We assessed the faecal mycobiota of South African adults living in rural (n = 50) and urban (n = 50) locations by assessing stool samples (see details regarding sample recruitment in Methods). We recruited an equal number of male and female volunteers. The volunteers from the rural areas were from two villages in the Limpopo province of South Africa. These villages are located roughly 500 km from the urban site in Pretoria (Fig. 1a). To gain insights regarding predictive variables, which may shape the gut microbiome, detailed questionnaires were distributed to all volunteers (Additional file 1 Questionnaire). The volunteers from Ha-Ravele and Tshikombani villages (population size of roughly 200,000, representing the rural participants) were on average 24 years (mean ± 6.3). Volunteers from Pretoria (total population of approximately 2.1 million) were on average 31 years (mean ± 9.1). The mean age of the participants was 27 years (mean ± 7.9) across all samples. The average height and weight of the participants was 1.64 m (mean ± 0.1) and 69.8 kg (mean ± 17.6), respectively. The average body mass index (BMI) of all participants was 26.02 kg/m2 (mean ± 6.4), resulting in a group of participants classified as overweight and obese, less than 15.9% of participants were smokers.

Fig. 1

Geographic locations and diversity estimates (a) The three sampling locations in Gauteng (Pretoria) and Limpopo (Ha-ravele and Tshikombani) provinces of South Africa (b) The differences in mycobiota species richness between the two locations, gender and age group and, (c) The relative abundances of taxa at phylum and class levels within each location. The abundance of each taxon was calculated as the percentage of sequences per gender (RF = Rural female, RM = Rural male, UF = urban female and UM = Urban male) from each location for a given microbial group. The group designated as ‘Unknown’ encompasses unclassified sequences together with classes representing > 0.1% of the total sequences. The bar size represents the relative abundance of specific taxa in the particular group, with colours referring to taxa according to the legend. The map was sourced from d-maps.com (https://d-maps.com/carte.php?num_car=23735&lang=en) and manually edited to indicate the study locations

Amplicon sequence data from 95 volunteers (samples from 5 rural volunteers were excluded due to low quality reads) generated 5,936,454 raw reads. Of these, 1,636,180 reads were clustered into OTUs at 3% divergence (97% similarity) and 1911 OTUs were taxonomically classified. The resulting accumulation curves showed reasonable sequence saturation at a regional level (Additional file 3 Fig. S1).

A higher proportion of fungal OTUs were unique to location with urban and rural samples accounting for 47.9 and 45.3% of reads, respectively (Additional file 4 Fig. S2). Fungal species richness was higher in the stool samples of urban volunteers compared to rural volunteers (Fig. 1b). However, there were no significant differences in species richness based on location (W = 915, p-value = 0.118) and sample type (Kruskal-Wallis chi-squared = 5.103, df = 3, p-value = 0.164). A significant difference was detected in species richness between the three age groups (Kruskal-Wallis chi-squared = 12.215, df = 2, p-value = 0.002).

Two ubiquitous fungal phyla in urban and rural locations

Overall, four distinct fungal phyla were detected in urban and rural gut mycobiomes, based on sequences with relative abundances above 0.1% (Fig. 1c). The majority of sequences were assigned to members of the phyla Ascomycota and Basidiomycota, that constituted 80.7 and 6.1% of the total relative abundance, respectively. Unknown sequences constituted 12.9% of the total relative abundance of the mycobiome. The relative abundance of Ascomycota was almost equal between the urban (40.1%) and rural (39.8%) location. Whereas, the relative abundance of Basidiomycota was higher in the urban (4.4%) compared to the rural (1.7%) location. Nonetheless, no significant differences were detected in taxa relative abundance between location (W = 8720, p-value = 0.3826), sample type (Kruskal-Wallis chi-squared = 0.9454, df = 3, p-value = 0.8145) and age group (Kruskal-Wallis chi-squared = 0.29477, df = 2, p-value = 0.863) at phylum level. In total, 8 distinct fungal classes were identified with Saccharomycetes constituting the majority of sequences (52.6%) followed by Dothideomycetes (20.9%), Eurotiomycetes (4.3%), Sordariomycetes (4%) and Tremellomycetes (3.98%). The relative abundance of Saccharomycetes was higher in the rural (26.8%) compared to the urban (25.8%) location. Whereas, Dothideomycetes relative abundance was higher in the urban (14.8%) compared to the rural (6.1%) location. However, no significant differences were detected in taxa relative abundance between location (W = 3480, p-value = 1), sample type (Kruskal-Wallis chi-squared = 2.7506, df = 3, p-value = 0.4317) and age group (Kruskal-Wallis chi-squared = 0.55361, df = 2, p-value = 0.7582). Pichia dominated our participants (18.9% of total relative abundance), followed by Candida (18.2%) and Cladosporium (6.1%). Whereas, unknown fungal genera accounted for 18.1% of the relative abundance. The relative abundance of Pichia was higher in the rural (12.7%) compared to the urban (6.2%) location. Whereas, the relative abundance of Candida and Cladosporium was higher in the urban (10.9 and 4.2%, respectively) compared to the rural (7.3 and 2%, respectively) location. However, no significant difference was found between taxa abundance at the genus level for the gut mycobiota of rural and urban participants based on location (W = 4417, p-value = 0.9936), sample type (Kruskal-Wallis chi-squared = 1.1117, df = 3, p-value = 0.7742) and age group (Kruskal-Wallis chi-squared = 0.51268, df = 2, p-value = 0.7739).

To assess the distribution and contribution of taxa in a given sample to the overall community composition, we assessed the local contribution to beta diversity (LCTBD). In line with findings from alpha diversity (observed species richness) analyses, we found that samples from urban volunteers (greater species richness) had a more significant contribution to the overall community diversity (p-Value < 0.05). Samples with high LCTBD had a high abundance of Basidiomycota and other unknown taxa. In contrast, only two samples from the rural location contributed more significantly to overall community diversity (Additional file 5 Fig. S3).

Distinct mycobiota among urban and rural volunteers unrelated to gender

Differences in the fungal community structure between the rural and urban localities were visualized in an non-metric multi-dimensional scaling (NMDS) plot (Fig. 2a). Urban and rural samples formed distict clusters [permutational multivariate analysis of variance (PERMANOVA) (R = 0.070; p-Value = 0.0001), ANOSIM (R = 0.43, p-Value = 0.001) and ADONIS (R = 0.07034 p-Value = 0.0001)]. However, male and female samples did not cluster separately. Pairwise analysis using PERMANOVA showed that there was no significant difference between gender within each location (R = 0.018; R = 0.023; respectively and p-Value > 0.4 for both). Nevertheless, there was a significant difference between the gut mycobiota of female and male participants between the two locations (R < 0.074 for all; p-Value = 0.001 for all).

Fig. 2

Overview of mycobiota structure and significant environmental drivers (a) The non-metric multidimensional scaling (NMDS) plot based on Bray–Curtis dissimilarity and, (b) Redundancy analysis (RDA) showing community structure in response to four selective variables. The filled shapes reflect fungal community composition in the different locations, with colours referring to location and the different explanatory variables according to the legend

Ecological drivers of gut mycobiota

Redundancy analysis (RDA) was performed to determine which predictor variables significantly explained the variation in fungal community composition (Fig. 2b). Four predictive variables were significant (r > 0.2; p-Value < 0.05) drivers of community composition and structure. Variation partitioning analysis demonstrated that only 3% community variations were explained by these four predictive variables (Additional file 2 Partition of variance in RDA). Predictive variables which included; breastfeeding, smoking, mode of birth and location; all of which significantly influenced the fungal community composition.

We conducted correlation analyses to explore the relationships among dominant gut species. Our results showed a few highly positive correlations in the rural participants: between Mucor and Dipodascus, Mucor and Naganishia, Clavispora and Lentendrea, and between Udeniomyces and Lentendrea (Fig. 3). Whereas, the strongest negative correlation was found between Dipodascus with Trichoderma, Dipodascus with Ascotricha and Dipodascus with Chalastospora. Within the urban cohort, Xeromyces and Agaricus, Diutina and Clavispora, and Dekkera and Diutina exhibited the strongest positive correlations (Fig. 3). The strongest negative correlations were detected between Clavispora with Filobasidium, and with Verrucaria and Malassezia.

Fig. 3

Correlations occurring between fungal taxa in (a) rural and (b) urban fungal mycobiota with P < 0.05 after FDR adjustment. Red squares represent significant negative correlations and blue squares represent significant positive correlations. The darker colours represent stronger correlations and non-significant correlations have been excluded from the plot

Biomarker taxa

Linear discriminant analysis (LDA) and the linear discriminant analysis effect size (LEfSe) [32] test for biomarkers was used to detect taxa that showed the strongest effect on group differentiation (Fig. 4a). OTU level analysis uncovered 14 urban-associated species from 10 genera. Whereas, 17 rural-associated species from 11 genera, were detected as possible biomarkers. The most abundant rural-associated biomarker genera were Hypopichia and Dipodascus, with species Hypopichia burtonii and Dipodascus geotrichum being the most abundant (Fig. 4b). The urban-associated biomarkers were dominated by the class Tremellomycetes and genera Dekkera and Hannaella. Species Dekkera bruxellensis and Hannaella sinensis dominated the urban-associated biomarkers.

Fig. 4

The results of Linear discriminant analysis (LDA) effect size (LefSe) analysis of rural and urban gut mycobiota (a) The cladogram shows the output of the LEfSe algorithm, which identifies taxonomically consistent differences between rural (Ha-ravele and Tshikombani villages) and urban (Pretoria) fungal community members, respectively. Taxa with nonsignificant differences are represented as yellow circles and the diameter of the circle is proportional to relative abundances (b) The histogram of the LDA scores was computed for differentially abundant taxa between the rural and urban gut mycobiota. The bar size represents the effect of the size of specific taxa in the particular group at species level

Discussion

The results from this study suggest that the gut mycobiome of the South African population is structured by geography and lifestyle. This finding is supported by the clustering of a large proportion of the fungal OTU’s into discrete rural and urban groups within the Venn diagram. Only a small percentage of OTUs were shared between the two populations, which may suggest that factors such as the environment, age and diet may play a role in shaping the differences in OTU clustering. These results were further corroborated by NMDS and PERMANOVA analyses, which not only show that both populations cluster distinctly and share just a few taxa, but that they have diverse fungal community composition consistent with rural and urban locality. Redundancy analysis (RDA) also predicted that the measured variables account for only a small proportion of the variability (R2 = 0.2) in fungal community composition. Indicating that other undetermined factors may be driving differences in fungal community composition between the two groups.

Several studies have investigated the healthy human mycobiome [6, 7, 19, 25, 30, 33]. In these studies, geography was not considered as a potential factor structuring the gut mycobiota. For instance, previous studies found no association between host phenotypic characteristics with mycobiome profile [6]. Nash et al. (2017) also suggests that diet, the environment, diurnal cycles, and host genetics may substantially influence the human gut mycobiome. However, the finding that the majority of the variation could not be explained by their metadata does suggest that other environmental factors, such as geography, may contribute to structuring the human mycobiome [6].

Our study provides the first results showing the importance of geography in African populations. Geographic locality may be associated with different environmental factors, such as different climatic regimes, which may effect structural changes in the mycobiota. For example, climate significantly influences vegetation and farming practices and leads to region specific diets. These region-specific diets may ultimately influence the gut mycobiota. This is a reasonable prediction given previous findings showing that fungi have climate dependent biogeographic patterns [34, 35]. These patterns are likely to determine the type of fungi individuals may be exposed to, which may in turn impact the colonization of fungi in the human gut. The most abundant rural-associated biomarker species found in this study, Dipodascus geotrichum, is ubiquitous in nature [36] whereas, Hypopichia burtonii is commonly isolated from corn, wheat, and rice [37]. The urban-associated biomarkers were dominated by the species Dekkera bruxellensis, which are commonly isolated from fermented food such as wine, beer, feta cheese and sour dough [38-40]. In contrast, Hannaella sinensis is commonly isolated from plants and soil [41, 42]. The staple diet of the rural South African population primarily consists of a corn-based porridge (called ‘pap’). It is therefore not uncommon for a fungal species commonly isolated from corn to be a dominant biomarker for the rural population. Conversely, the urban population diet was more diverse and included fermented foods such as wine, sour dough bread and feta cheese, which are commonly available in supermarkets. Thus, the species Dekkera bruxellensis was identified as a dominant biomarker in the urban population.

In addition to geographic location, we found that smoking, mode of birth and breastfeeding significantly influenced gut fungal communities. Several studies have previously reported that these factors may significantly influence the initial colonization, subsequent composition and structure of bacterial members of the human gut microbiome [28, 43–45]. Suhr et al. (2016) and Hallen-Adams et al. (2017) investigated the gut mycobiome of two cohorts that were exclusively on a vegetarian or a western diet. These studies found that the distribution of fungi differed considerably between the two cohorts [7, 46]. Plant-associated Fusarium, Malassezia, Penicillium and Aspergillus species were detected at higher abundances within the vegetarian cohort, compared to the cohort on a conventional diet. The finding that smoking affected fungal community composition and structure is supported by several studies [47-49]. The approximately 4000 chemical compounds produced by cigarettes have been shown to alter the composition of the gut microbiome [47, 49–52]. The reported increase of Clostridia induced by smoking in murine models has also been indirectly confirmed in humans where an increased rate of C. difficile infection was greater in former and current smokers compared to never smokers [51]. Moreover, the abundance of the fungus Candida tropicalis has also been reported to be significantly higher in C. difficile infection patients compared to healthy individuals [53]. The abundance of C. tropicalis has also been detected to be positively correlated with levels of anti-Saccharomyces cerevisiae antibodies (ASCA) [53]. In our study C. tropicalis was detected to be higher in individuals who smoke compared to non-smokers whereas the inverse was true for S. cerevisiae. These findings may confirm the antagonistic association between the species C. tropicalis and S. cerevisiae, as previously reported by Hoarau et al. (2016) [53].

Most studies have identified the genera Candida, Saccharomyces, Malassezia and Aspergillus as the three most abundant in the gut of healthy individuals [6, 7, 25]. To the best of our knowledge, our study is the first to report Pichia as one of the top four (Pichia, Candida, Cladosporium and Paraconiothyrium) most abundant genera in the human gut mycobiome. This may be due to several factors including differences in cohort characteristics (e.g., geographical location, diet, genetic predisposition and climate). Pichia have been identified as both constituent members of the human oral [54, 55] and gut microbiome [33]. Mukherjee detected a 1:1 abundance ratio in the oral mycobiome of individuals when Candida and Pichia were present together [55]. Pichia was also observed to have an antagonistic effect against Candida, Fusarium and Aspergillus [55].

The yeast genera, Pichia, Candida and Cladosporium, dominated the South African gut mycobiome. Our findings agree with previous studies which show that members of the Aspergillus, Candida, Debaryomyces, Malassezia, Penicillium, Pichia, and Saccharomyces genera were the most recurrent and/or dominant fungal genera [33, 46, 56]. In contrast to previous findings, our data indicate higher relative abundances of Cladosporium, detection of Mucor and the absence or low abundance of genera such as Cyberlindnera, and Galactomyces [6, 19, 57]. Previous studies found that the gut mycobiome of a cohort from Houston, Texas, was dominated by Saccharomyces, Malassezia and Candida [6]. By contrast, the genus Malassezia was not detected in the gut mycobiome of a Pennsylvania cohort, which was instead dominated by the genera Saccharomyces and Candida [19]. Differences in study methodologies may be a source of these conflicting findings [6]. One study amplified the Internal Transcribed Spacer 2 (ITS2) region of the fungal rRNA gene [6], and the second amplified the ITS1 region [19]. Studies similar to the work by Gardes et al. (1993) and White et al. (1990), where ITS1F and ITS2 primer sets were used to amplify the ITS2 region, did not detect Malassezia [58, 59]. The second reason for the observed differences has been attributed to differences in cohort characteristics, such as diet and/or geographical location. Strati (2016) and Raimondi’s (2019) investigating cohorts in Italy, detected same dominant fungal genera [31, 57], and the investigation of cohorts in two different states in the USA observed different results [6, 45]. We used ITS1 and ITS4 in this study and found that the genera Pichia, Candida and Cladosporium dominated the urban cohort, whereas genera Pichia, Candida and Aspergillus dominated the rural cohort. The dominant taxa identified in urban and rural locations further support our assertion that geographic location plays a major role in the observed differences.

Candida albicans was the most dominant taxon in our cohort and is frequently reported as the most abundant Candida species in both diseased [60] and healthy individuals [61]. Candida spp. not only colonize the gut [19, 33] but several other body sites, including the oral cavity [54, 62], vagina [63], and skin [64, 65]. However, Candida are autochthonous to the mammalian digestive tract and species including Candida albicans, C. tropicalis, C. parapsilosis, and C. glabrata may grow and colonize at 37 °C [7]. A review of the literature suggest that C. albicans carriage in healthy individuals ranges from 30 to 60% [66] and that living mammals are considered a niche for these species as they are not found in significant concentrations in soil, food or air [67, 68]. Raimondi et al., (2019) reported that C. albicans was frequently detected and dominated the cultivable mycobiota of different faecal samples [31].

Conclusions

This study provides the first insight into the importance of geography and lifestyle factors on the gut mycobiome in rural and urban locations in Africa. We found that fungi in the gut display distinct patterns consistent with geographic locality. Redundancy analysis showed that several lifestyle factors were major drivers explaining the distinct community structure. The results of biomarker analysis revealed several ecologically important fungal taxa, which were unique to individuals from urban and rural areas. The finding that certain taxa may be biomarker species have potential consequences for certain groups including immunocompromised individuals living in rural and urban locations. Increases in the abundances of these taxon may lead to deleterious effects on the health of these groups. Such findings provide a valid basis for the development of novel therapeutics or preventative measures reliant on modulating the gut mycobiome.

Methods

Participant enrolment criteria for urban and rural areas

Volunteers were recruited from two rural locations and one urban location. For rural volunteers, we recruited individuals following traditional diets, with generally low levels of processed foods. Urban cohorts reported mixed diets and increased consumption of processed foods. Volunteers from the Ha-Ravele (females; n = 1 and males; n = 15) and Tshikombani (females; n = 24 and males; n = 10) villages located in the Vhembe District of the Limpopo Province comprised the rural cohort. Both villages are approximately 391 km and 439 km, respectively, from the closest city (Pretoria). This city, in the Gauteng province of South Africa, served as the urban sampling area. In total, 100 stool samples were collected from healthy volunteers. These samples were equally divided between gender and locality [i.e. rural (25 males and 25 females) and urban (25 males and 25 females)]. Participants were categorized by age into young adults (ages 18–27 years; n = 61), middle-aged adults (ages 28–37 years, n = 28), and older adults (aged older than 37 years, n = 11). The height of the participants was measured in meters, weight in kilograms and the BMI was calculated using the formula BMI = kg/m2 [69] where kg is a person’s weight in kilograms and m2 is their height in meters squared. Participants were categorized by BMI into Underweight = BMI < 18.5 (n = 10), Normal weight BMI = 18.5–24.9 (n = 46), Overweight BMI = 25–29.9 (n = 22) and Obese = BMI of 30 or greater (n = 22) [69]. Self-stool collection kits were provided to all volunteers (Easy Sampler® Stool collection Kit, Hounisen Lab Equipment A/S, Skanderborg, Denmark).

Inclusion and exclusion criteria

The participants were all healthy adults age 18–50 years. Volunteers reporting antibiotic use/other treatments within 6 months prior to participating in the study and sample collection were excluded from the study. Similarly, individuals who had been diagnosed with any inflammatory-related bowel diseases or gastrointestinal diseases within 6 months prior to sample collection were excluded from the study.

DNA extraction

DNA was isolated using the PowerSoil DNA Isolation Kit (MO BIO Laboratories Inc., Carlsbad, CA) following the manufacturer’s specifications with minor modifications. Briefly, approximately 0.25 g of stool sample was transferred into the Power-Bead tubes using a sterile disposable wooden spatula (Lasec Laboratories, RSA). The sample was homogenized by gently vortexing the tubes for 10 s before adding 60 μL of the lysis buffer. This was then incubated for 30 min. at 55 °C prior to centrifugation at room temperature for 30 s at 10,000 x g. The supernatant from this step was transferred to sterile 2 mL tubes and 250 μL of inhibitor removal reagent was added to this. The samples were incubated on ice for 5 min., thereafter approximately 1.2 mL of binding buffer was added. Next, 70% ethanol (500 μL) was added and the contents precipitated by centrifugation at room temperature for 60 s at 10,000 x g. The DNA was eluted with 100 μL filter-sterilised autoclaved Millipore water and quantified using the NanoDrop™ 2000/2000c Spectrophotometer (Thermo Scientific, Waltham, MA, USA). The quality of isolated DNA was confirmed by agarose gel electrophoresis, on 1% (w/v) agarose gel in 1 X TAE buffer (0.2% [w/v] Tris, 0.5% [v/v] acetic acid, 1% [v/v] 5 M EDTA [pH 8]) at 90 Volts for 45 min. in a BioRad Sub-Cell® GT gel electrophoresis system with gel red visualising agent. The gel was visualised using the BioRad Gel Doc system and viewed with a UV Trans-illuminator.

ITS gene region amplification, sequencing and data processing

The internal transcribed spacer (ITS) region (420 to 825 bp) was amplified using fungal-specific primers [70]: ITS1F (5′-CTTGGTCATTTAGAGGAAGTAA-3′) [58] and ITS4 (5′-TCCTCCGCTTATTGATATGC-3′) [59] with barcode inserted on the forward primer. Briefly, the HotStarTaq Plus Master Mix Kit (Qiagen, USA) was used for the PCR amplification reaction (94 °C for 3 min., followed by 30 cycles of 94 °C for 30 s, 53 °C for 40 s, 72 °C for 1. min and final elongation step at 72 °C for 5 min.). The PCR products were checked in 2% agarose gel to determine the success of amplification and the relative intensity of bands. Amplicons from different samples were pooled to equal proportions based on their molecular weight and DNA concentrations. The pooled DNA was purified of short fragments using Agencourt Ampure beads (Agencourt Bioscience Corporation, USA). Then the pooled and purified PCR product was used to prepare Illumina library. Paired end 2 × 250 bp sequencing was performed on an Illumina MiSeq instrument (Illumina Inc., San Diego, CA, USA) at Mr. DNA (Shallowater, TX, USA).

The resultant data were analysed using the Quantitative Insights into Microbial Ecology (QIIME2) software version 2018.8.0 [71]. Demultiplexed sequences were merged and assessed for quality. Sequences shorter than 200 bp, with quality scores below 25, containing more than two ambiguous characters or more than one mismatch to the sample-specific barcode or the primer sequences, were excluded from further downstream analyses. Sequences were denoised, chimeric sequence removed and operation taxonomic units (OTUs) were defined by clustering at 3% divergence (97% similarity) using USEARCH v11 [72]. Taxonomies were assigned to each OTU using the UNITE (release 7_99) databases for fungi [73]. Singletons were excluded, and each sample was randomly subsampled (rarefied) to the same number of sequences per sample (17980).

Statistical analyses

All statistical analyses were performed in R version 3.5.1 using R studio [74, 75]. Alpha diversity (observed richness), together with rarefaction curves were calculated and visualized using the R packages “phyloseq” and “ggplot”. First, the Shapiro-Wilk’s test was used to determine whether the data had a normal distribution [76]. Subsequently, the unpaired two-sample Wilcoxon rank sum test [77, 78] was applied to determine significant differences between the alpha diversity indices using the R packages “dplyr” version 0.4.3 and the “ggpubr” version 0.1.8 [79, 80]. In these analyses, the rural or urban location was specified as a random factor.

The R packages “phyloseq” [81] and “microbiomeseq” [82] were used to calculate and visualize relative taxa abundance at phylum and class level. OTU abundance was transformed to relative abundance and taxa with relative abundance less than 0.1% were removed. The Wilcoxon rank sum test was applied to determine significant differences between taxa relative abundance in the urban and rural samples. Whereas, the Kruskal-Wallis test [81, 83] was applied to determine significant differences in taxa relative abundance between the four sample types (rural female, rural male, urban female and urban male).

The LCBD was calculated according to [84]. The LCBD describes the degree of uniqueness of a given sample in relation to the overall community composition. The taxa abundance was normalized to obtain the proportion of most abundant taxa per sample. Location was used as the grouping variable and the Hellinger method [85] was used for the dissimilarity coefficients calculation.

Pairwise similarities among samples were calculated using the Bray–Curtis index of similarity. The resulting matrix was represented visually in a nonmetric multidimensional scaling (NMDS) plot to observe community structure. Using the vegan package [86], a permutational multivariate analysis of variance (PERMANOVA) [87] based on 9999 permutations of the data, was performed to test whether differences between sample groupings in the NMDS ordinations were statistically significant. Microbial community similarities and the homogeneity of dispersion between the rural and urban sample groups were tested using the ANOSIM and ADONIS tests, respectively [88, 89].

The effect of the different recorded environmental factors on fungal community composition and structure was determined through redundancy analysis (RDA). The contribution of highly correlating OTUs (p-Value < 0.05) with redundancy axes was identified using the envfit function from the R package vegan [86].

Fungal-fungal relationships were interrogated using SparCC [90]. Correlation was based on measuring the linear relationship between log transformed abundances. First, data were filtered to remove OTUs that had less than 2 reads on average. SparCC was used to generate true correlation coefficients from which pseudo p-values were calculated. The calculate pseudo p-values were false discovery rate (FDR) adjusted [91] and the correlation matrix was visualized using the “corrplot” function [92] in R.

Potential biomarker taxa which differed in abundance and occurrence between the two geographic groups were detected by linear discriminant analysis (LDA) effect size (LEfSe) [32]. The LEfSe was calculated using the online Galaxy web application [93] with the Huttenhower lab’s tool (https://galaxyproject.org/learn/visualization/custom/lefse/). First the nonparametric factorial Kruskal–Wallis sum rank tests (alpha = 0.01) was used to detect differential abundant features (at genera, family, class and phylum level) within the two geographic locations (rural and urban). The phylogenetic consistency was then tested using the pairwise Wilcoxon rank-sum tests (alpha = 0.01). Finally, the effect size of each differentially abundant feature was estimated using the LDA. The all-against-all classes were compared (most stringent) and a linear discriminant analysis score value of 2.0 was chosen as threshold for discriminative features.

Additional file 1. Questionnaire Details of the questionnaire provided to participants prior to enrolment in the study. The questionnaire details essential required information, clinical information, voluntary dietary information and questions regarding data sharing.

Additional file 2. Results from the partition of variance analysis in RDA

Additional file 3: Fig. S1. Rarefaction plot showing sequencing coverage. The estimated average sequence coverage of high-quality paired end reads after quality control assessed using Nonpareil (in alignment mode).

Additional file 4: Fig. S2. Venn diagram showing the unique and shared phylotypes for samples collected from urban and rural participants.

Additional file 5: Fig. S3. Taxa abundance data was normalised to obtain the proportion of most abundant taxa per sample. The diameter of the points at the bottom of the plot corresponds to the magnitude of the LCBD value for a particular sample. The bars correspond to taxa that are most abundant with the top taxa sharing a bigger portion of the bar for each sample.