First Insights into the Gut Microbiota of Mexican Patients with Celiac Disease and Non-Celiac Gluten Sensitivity

Gluten-related disorders (GRDs) are common chronic enteropathies and increasing evidence suggests an involvement of the gut microbiota. We examined the gut microbiota in Mexican people afflicted with GRDs. Ultra-high-throughput 16S marker sequencing was used to deeply describe the duodenal and fecal microbiota of patients with celiac disease (CD, n = 6), non-celiac gluten sensitivity (NCGS, n = 12), and healthy subjects (n = 12) from our local area. Additionally, we also investigated the changes in gut microbiota after four weeks on a gluten-free diet (GFD) in a subset of patients from whom paired samples were available. Despite a high inter-individual variability, significant differences in various microbial populations were identified. The linear discriminant analysis (LDA) effect size (LEfSe) method revealed that the genus Actinobacillus and the family Ruminococcaceae were higher in the duodenal and fecal microbiota of NCGS patients, respectively, while Novispirillum was higher in the duodenum of CD patients (p < 0.05, LDA score > 3.5). Interestingly, paired samples from NCGS patients showed a significant difference in duodenal Pseudomonas between the baseline period (median: 1.3%; min/max: 0.47–6.8%) and the period after four weeks on GFD (14.8%; 2.3–38.5%, p < 0.01, Wilcoxon signed-rank test). These results encourage more research on GRDs in México.


Introduction
The gut microbiota is comprised of thousands of microbial species that vary widely among individuals [1] and also over time within the same individual due to environmental factors such as dietary patterns [2]. The gut microbiota helps modulate the immune system [3] and has been associated with diseases related to the alimentary tract such as obesity and inflammatory bowel diseases [4]. Given the close relationship between the immune system and microorganisms inside the gut, it is believed that most disorders of the digestive tract bear some relationship with the gut microbiota although a cause-effect relationship can hardly be established [5].

Bioinformatics
The open-source bioinformatics pipeline Quantitative Insights into Microbial Ecology [30] v.1.8 was used for most of the core analyses. Operational Taxonomic Units (OTUs) were chosen using two approaches. First, using the pick_open_reference_otus.py accordingly to the suggestions by Rideout et al. [31]. This approach is capable of detecting OTUs that are not necessarily represented in the reference databases. Further taxonomic and diversity analyses were performed using all OTUs (i.e., the full OTU table) and a filtered OTU table (OTUs with <0.005% of all sequences were removed as suggested by Navas et al. [32]. Second, using the pick_closed_reference_otus.py to then be able to use the OTU table for the prediction of functional metagenome using Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) [33]. The GreenGenes database [34] at 97% similarity was used as the reference 16S database. All sequence and metadata information are publicly available (NCBI, PRJNA401920).

Statistical Analysis
A chi-squared test was used to compare the frequencies (e.g., the proportion of women, number of patients showing a clinical improvement) and the non-parametric Kruskal-Wallis test was used for comparison of health parameters (e.g., blood parameters) and microbial groups. The linear discriminant analysis (LDA) effect size (LEfSe, [35]) was used to determine the organisms that explain the differences in microbiota. Please note that in LEfSe, the idea is that the significant biomarkers (in this case microbial phylogroups) are ranked based on the effect size (the magnitude of the variation) rather than on the statistical significance. When comparing two sets of data (e.g., before and after GFD), the Wilcoxon signed-rank test or the Mann-Whitney test were used. The unique fraction metric (UniFrac) was used to measure the phylogenetic distance among taxa [36]. Both weighted and unweighted UniFrac were calculated and analyzed using Principal Coordinate Analyses (PCoA) because they can lead to different insights into the factors that shape the composition of bacterial communities [37,38]. The ANOSIM and Adonis tests were used to determine whether the grouping of samples by a given category (e.g., health status) is statistically significant based on the UniFrac distances. Two age groups (young < 35 years; old > 35 years) and two body mass index (BMI) groups (low < 24.5; high > 24.5) were created to evaluate the potential contribution of these factors to the similarity of bacterial communities. STAMP [39] was used to analyze PICRUSt data using non-parametric tests.

Subjects
A total of six patients with CD, twelve patients with NCGS and twelve control subjects were successfully enrolled over the six months enrollment period (Table 1). Please note that not all samples were obtained from all subjects mainly because of the lack of compliance, especially with the submission of stool samples. Among the CD patients, one had a Marsh I classification, two had Marsh II and three had Marsh IIIa. The impact of these varying baseline scores on clinical development and gut microbiota is uncertain but something to look for in future studies. The history of CD among relatives was more common in CD patients, CD patients had lower BMIs and hemoglobin levels and higher intraepithelial lymphocyte counts ( Table 2). There was no difference between the CD and NCGS patients at baseline with regards to abdominal pain and bloating ( Table 2).  Yes  19  CD  25  23  Woman  Yes  NA  20  CD  47  25  Woman  Yes  NA  23  CD  73  21  Woman  Yes  Only on GFD  1  NCGS  23  28  Woman  Yes  Yes  3  NCGS  21  24  Woman  Only baseline  NA  5  NCGS  24  25  Woman  Yes  Only on GFD  6  NCGS  23  29  Woman  Yes  Yes  7  NCGS  22  25  Woman  Yes  NA  8  NCGS  24  27  Woman  Yes  NA  10  NCGS  27  23  Man  Yes  Yes  11  NCGS  23  29  Man  Yes  Only baseline  13  NCGS  37  31  Woman  Yes  Yes  17  NCGS  59  19  Woman  Yes  Yes  22  NCGS  34  26  Woman  Only baseline  NA  24  NCGS  38  24  Woman  Yes  Only baseline  2  Control  23  33  Man  Only baseline  NA  4  Control  24  33  Man  Only baseline  NA  12  Control  23  24  Woman  Only baseline  NA  14  Control  25  23  Woman  Only baseline  NA  15  Control  26  21  Woman  Yes  Yes  21  Control  24  29  Man  Yes  NA  25  Control  45  28  Woman  Yes  Only baseline  26  Control  64  24  Man  Yes  Yes  27  Control  23  25  Man  Yes  Only baseline  28  Control  39  25  Man  Yes  NA  29  Control  58  26  Woman  Yes  NA  30  Control  42  27  Woman Yes NA

16S Sequencing and Taxonomic Classification of Sequence Reads
A total of 2.3 million (biopsies, n = 30) and 1.5 million (fecal, n = 14, many patients refused to provide a stool sample) good-quality 16S rDNA sequences (median length: 300 base pairs) were obtained from the baseline samples and used for OTU picking and further analyses. A total of 32,800 OTUs were originally detected using the open OTU picking approach (unfiltered OTU table); the removal of low-abundant OTUs (i.e., OTUs with <0.005% of total reads) yielded 975 and 916 OTUs (only~3% of all original OTUs) in biopsy and fecal samples, respectively. It is outside the scope of this current work to discuss the consequences of removing low abundant OTUs but please be aware that the so-called rare microbes may in fact be keystone species regulating the function of different microbial environments, including host-associated microbiomes [40].

Microbiota in Duodenal Biopsy Samples at Baseline
Overall 16S reads were assigned to a total of 27 phyla in all samples but only five phyla (Proteobacteria, Firmicutes, Actinobacteria, Bacteroidetes, and Fusobacteria) comprised the vast majority (>90%) of reads in most samples (Figure 1), as shown elsewhere. At the phylum level, there was a significantly lower abundance of Bacteroidetes (p = 0.022, Kruskal-Wallis test) and Fusobacteria (p = 0.052) in duodenal biopsies from CD patients (n = 30, Figure 2). This lower abundance of Bacteroidetes and Fusobacteria in CD patients was also true when analyzing the duodenal microbiota of women only (n = 22, p = 0.028 and p = 0.067, respectively).

Microbiota in Duodenal Biopsy Samples at Baseline
Overall 16S reads were assigned to a total of 27 phyla in all samples but only five phyla (Proteobacteria, Firmicutes, Actinobacteria, Bacteroidetes, and Fusobacteria) comprised the vast majority (>90%) of reads in most samples (Figure 1), as shown elsewhere. At the phylum level, there was a significantly lower abundance of Bacteroidetes (p = 0.022, Kruskal-Wallis test) and Fusobacteria (p = 0.052) in duodenal biopsies from CD patients (n = 30, Figure 2). This lower abundance of Bacteroidetes and Fusobacteria in CD patients was also true when analyzing the duodenal microbiota of women only (n = 22, p = 0.028 and p = 0.067, respectively).   LEfSe analysis confirmed the finding of statistically significant differences in various bacterial groups among the three groups of subjects at the baseline ( Figure 3). For instance, there was a higher abundance of Actinobacillus (Gammaproteobacteria), Finegoldia (Clostridia), and the phylum TM7 in NCGS patients, while Sphingobacterium (Bacteroidetes) was higher in the healthy subjects ( Figure 3). The separate LEfSe analysis of samples from women confirmed the higher abundance of TM7 in NCGS patients and Sphingobacterium in healthy subjects and also revealed significant differences in various other bacterial groups (e.g., women with CD were deprived of Campylobacterales, Paraprevotellaceae, and Fusobacteriaceae; see Supplementary Figure S1). The health status of the patients was not related to significant differences in any index of richness or diversity with the exception of Shannon diversity index (lower in CD patients; Table 3 and Supplementary Table S1). This overall lack of difference in alpha diversity was also true when only analyzing samples from women (n = 22).  LEfSe analysis confirmed the finding of statistically significant differences in various bacterial groups among the three groups of subjects at the baseline ( Figure 3). For instance, there was a higher abundance of Actinobacillus (Gammaproteobacteria), Finegoldia (Clostridia), and the phylum TM7 in NCGS patients, while Sphingobacterium (Bacteroidetes) was higher in the healthy subjects ( Figure 3). The separate LEfSe analysis of samples from women confirmed the higher abundance of TM7 in NCGS patients and Sphingobacterium in healthy subjects and also revealed significant differences in various other bacterial groups (e.g., women with CD were deprived of Campylobacterales, Paraprevotellaceae, and Fusobacteriaceae; see Supplementary Figure S1). The health status of the patients was not related to significant differences in any index of richness or diversity with the exception of Shannon diversity index (lower in CD patients; Table 3 and Supplementary Table S1). This overall lack of difference in alpha diversity was also true when only analyzing samples from women (n = 22). LEfSe analysis confirmed the finding of statistically significant differences in various bacterial groups among the three groups of subjects at the baseline ( Figure 3). For instance, there was a higher abundance of Actinobacillus (Gammaproteobacteria), Finegoldia (Clostridia), and the phylum TM7 in NCGS patients, while Sphingobacterium (Bacteroidetes) was higher in the healthy subjects ( Figure 3). The separate LEfSe analysis of samples from women confirmed the higher abundance of TM7 in NCGS patients and Sphingobacterium in healthy subjects and also revealed significant differences in various other bacterial groups (e.g., women with CD were deprived of Campylobacterales, Paraprevotellaceae, and Fusobacteriaceae; see Supplementary Figure S1). The health status of the patients was not related to significant differences in any index of richness or diversity with the exception of Shannon diversity index (lower in CD patients; Table 3 and Supplementary Table S1). This overall lack of difference in alpha diversity was also true when only analyzing samples from women (n = 22).   The differences in duodenal microbiota at the phylum (e.g., Bacteroidetes) and lower taxonomic levels (e.g., Actinobacillus) were not enough to differentiate the bacterial communities as a whole, as evaluated by the PCoA plots of weighted and unweighted UniFrac distances (Supplementary Figure S2) and this was also true when only analyzing the samples from women. The PICRUSt predicted metabolic features with the lowest uncorrected p values were flavonoid biosynthesis, dioxin degradation, and riboflavin metabolism ( Figure 4) but after the Bonferroni correction, there was no significant difference in any metabolic feature.  The differences in duodenal microbiota at the phylum (e.g., Bacteroidetes) and lower taxonomic levels (e.g., Actinobacillus) were not enough to differentiate the bacterial communities as a whole, as evaluated by the PCoA plots of weighted and unweighted UniFrac distances (Supplementary Figure  S2) and this was also true when only analyzing the samples from women. The PICRUSt predicted metabolic features with the lowest uncorrected p values were flavonoid biosynthesis, dioxin degradation, and riboflavin metabolism ( Figure 4) but after the Bonferroni correction, there was no significant difference in any metabolic feature. To summarize the results for the baseline duodenal microbiota, we found significant differences in the relative abundance of several bacterial groups but these differences were not enough to modify the diversity parameters (with the exception of Shannon diversity indexes) or predicted metabolic features.

Microbiota in Fecal Samples at the Baseline
Only 14 samples were available for the analysis of fecal microbiota at the baseline. Fecal samples showed an unexpected high abundance of Firmicutes (~85%) and a low abundance of Bacteroidetes (~1%) regardless of the disease status ( Figure 1). Despite the low number of samples, there was a clear To summarize the results for the baseline duodenal microbiota, we found significant differences in the relative abundance of several bacterial groups but these differences were not enough to modify the diversity parameters (with the exception of Shannon diversity indexes) or predicted metabolic features.

Microbiota in Fecal Samples at the Baseline
Only 14 samples were available for the analysis of fecal microbiota at the baseline. Fecal samples showed an unexpected high abundance of Firmicutes (~85%) and a low abundance of Bacteroidetes (~1%) regardless of the disease status ( Figure 1). Despite the low number of samples, there was a clear higher abundance of fecal Ruminococcaceae in NCGS patients, and this difference was significant according to the LEfSe analysis (Supplementary Figure S3).

Effect of GFD on Duodenal Microbiota
All subjects had a GFD adherence above 90%. Sixty-seven percent of CD patients (4/6) and ninety percent of NCGS patients (9/10) reported a global improvement of symptoms after four weeks on GFD but this difference was not significant (p = 0.247, chi-squared test). There was also no difference in any other clinical or physiological parameter with the exception of abdominal pain (lower during GFD in CD patients, Table 4). An additional 2.9 million sequences (1.8 million from a total of 24 biopsy samples, and 1.1 million from a total of 12 fecal samples) were obtained from subjects on GFD. Paired samples were not obtained for all subjects mainly because of the lack of compliance, especially with the submission of stool samples (Table 1). Despite an apparently clear distinctive abundance and distribution of phyla in duodenum between the periods with and without dietary gluten ( Figure 5), there was no significant difference in the abundance of any taxa between the two periods of time (p > 0.1), likely due to the high inter-individual variability. Additional analyses of relative abundances in paired samples revealed that each group of patients (controls, CD, and NCGS) displayed a distinctive variation over-time after consuming the GFD for four weeks (for example, most NCGS patients displayed little change after the GFD period, Figure 6). Interestingly, and despite a relatively more stable division-wide composition, 9 out of 10 paired samples of patients with NCGS showed an increase in the duodenal Pseudomonas on the GFD (Figure 7, baseline median: 1.3%, min/max: 0.47-6.8%; median after four weeks on GFD: 14.8%, 2.3-38.5%, p < 0.01, Wilcoxon signed-rank test, only subject 7 showed a decrease in this group, from 4.3% to 2.5%). This difference in most individuals was specific for Pseudomonas and not for other members of the duodenal microbiota (Figure 7). In contrast, only half of the paired samples (3 out of 6) from CD patients showed increases in Pseudomonas but these increases were so pronounced that they also affected median values (Figure 7). Additional analyses revealed that the 16S sequences from Pseudomonas were not different among the groups of subjects (see Figure S4 in "3.1 Pseudomonas in duodenum" in the Supplementary Information), thus suggesting that taxonomically similar Pseudomonas populations react differently in the presence of similar environmental conditions, in this case, in the absence of dietary gluten. This is particularly relevant in a context of the ecological significance of microdiversity [41]. was specific for Pseudomonas and not for other members of the duodenal microbiota (Figure 7). In contrast, only half of the paired samples (3 out of 6) from CD patients showed increases in Pseudomonas but these increases were so pronounced that they also affected median values ( Figure  7). Additional analyses revealed that the 16S sequences from Pseudomonas were not different among the groups of subjects (see Figure S4 in "3.1 Pseudomonas in duodenum" in the Supplementary Information), thus suggesting that taxonomically similar Pseudomonas populations react differently in the presence of similar environmental conditions, in this case, in the absence of dietary gluten. This is particularly relevant in a context of the ecological significance of microdiversity [41].   was specific for Pseudomonas and not for other members of the duodenal microbiota ( Figure 7). In contrast, only half of the paired samples (3 out of 6) from CD patients showed increases in Pseudomonas but these increases were so pronounced that they also affected median values ( Figure  7). Additional analyses revealed that the 16S sequences from Pseudomonas were not different among the groups of subjects (see Figure S4 in "3.1 Pseudomonas in duodenum" in the Supplementary Information), thus suggesting that taxonomically similar Pseudomonas populations react differently in the presence of similar environmental conditions, in this case, in the absence of dietary gluten. This is particularly relevant in a context of the ecological significance of microdiversity [41].    LEfSe analysis of the taxa at the genus level confirmed the results on Pseudomonas and showed that other Proteobacteria (e.g., Stenophomonas and Novosphingobium) were significantly more abundant on the GFD, while Actinomycetaceae was lower before the GFD (Figure 8). On the other hand, LEfSe did not reveal any taxa significantly associated with either a specific health status or with diet as class and health status as a subclass. Interestingly, LEfSe analysis revealed that Brevundimonas (a very low abundant group, <0.5% of all reads) was significantly enriched in CD patients when analyzing health status as class and diet as a subclass. This result was mainly due to a higher abundance of Brevundimonas in CD patients on the GFD. Other factors such as age, sex, or BMI were not significantly associated with the abundance of any bacterial taxa accordingly to the LEfSe analysis; however, this lack of significance must be taken cautiously because of the low sample size in each subgroup of patients.
There was no significant difference in the bacterial richness or diversity in the duodenum when comparing the period at baseline and after four weeks on GFD or health status using either the full OTU table (Table 3) or the filtered OTU table (Table S1). The ANOSIM and Adonis tests revealed interesting results to the factors associated with the differences in microbial composition among the samples based on UniFrac distances (Table 5 and Supplementary Table S2). For example, the diet factor almost reached a level of significance when analyzing the weighted UniFrac distances (Table 5  and Table S2). Additionally, the grouping of duodenal samples based on health status was found to be statistically significant when using unweighted UniFrac distances and almost reached a level of significance when using weighted UniFrac distances (Adonis test, Table 5 and Table S2). The age of the patients also seemed to contribute to the separation of duodenal communities, especially when using the filtered OTU table (Table S2). These results were supported by significance in the ANOSIM test but the associated R values were very low (R < 0.10), indicating that the clustering of samples was relatively weak (Table 5 and Table S2). Please note that the analysis of both the full and the filtered Figure 7. The box plots showing the relative abundance of 16S rDNA reads corresponding to the three most abundant bacterial groups in the duodenum at the genus level. * Significantly higher compared to the period of gluten consumption (p < 0.05, Wilcoxon signed-rank test). Please note that 90% (9 out of 10) and 50% (3 out of 6) of paired samples of patients with NCGS and CD (respectively) showed an increase in Pseudomonas (see main text for more details on this). CD: celiac disease, Ctrl: controls, NCGS: non-celiac gluten sensitivity, GFD: gluten-free diet.
LEfSe analysis of the taxa at the genus level confirmed the results on Pseudomonas and showed that other Proteobacteria (e.g., Stenophomonas and Novosphingobium) were significantly more abundant on the GFD, while Actinomycetaceae was lower before the GFD (Figure 8). On the other hand, LEfSe did not reveal any taxa significantly associated with either a specific health status or with diet as class and health status as a subclass. Interestingly, LEfSe analysis revealed that Brevundimonas (a very low abundant group, <0.5% of all reads) was significantly enriched in CD patients when analyzing health status as class and diet as a subclass. This result was mainly due to a higher abundance of Brevundimonas in CD patients on the GFD. Other factors such as age, sex, or BMI were not significantly associated with the abundance of any bacterial taxa accordingly to the LEfSe analysis; however, this lack of significance must be taken cautiously because of the low sample size in each subgroup of patients.
There was no significant difference in the bacterial richness or diversity in the duodenum when comparing the period at baseline and after four weeks on GFD or health status using either the full OTU table (Table 3) or the filtered OTU table (Table S1). The ANOSIM and Adonis tests revealed interesting results to the factors associated with the differences in microbial composition among the samples based on UniFrac distances (Table 5 and Supplementary Table S2). For example, the diet factor almost reached a level of significance when analyzing the weighted UniFrac distances (Table 5  and Table S2). Additionally, the grouping of duodenal samples based on health status was found to be statistically significant when using unweighted UniFrac distances and almost reached a level of significance when using weighted UniFrac distances (Adonis test, Table 5 and Table S2). The age of the patients also seemed to contribute to the separation of duodenal communities, especially when using the filtered OTU table (Table S2). These results were supported by significance in the ANOSIM test but the associated R values were very low (R < 0.10), indicating that the clustering of samples was relatively weak (Table 5 and Table S2). Please note that the analysis of both the full and the filtered OTU table revealed similar results, thus suggesting that low-abundant OTUs did not play an important role in the separation or lack thereof of communities.    The variables included Diet (gluten/GFD), Disease (Control/CD/NCGS), Group (six groups of samples comprise this category: Control, CD and NCGS before and after four weeks of GFD), Age (Young/Old), BMI (high/low). p values that reached statistical significance (p < 0.05) or were close to reaching significance (p < 0.1) are highlighted in bold for better visualization.
Considering the differences in relative abundance among the different bacterial groups of the duodenal microbiota (e.g., the higher Pseudomonas on GFD in CD and NCGS patients), we hypothesized that the beta diversity analyses for different bacterial populations may offer clues regarding the effect on different factors such as the diet and health status. Therefore, we used the Adonis and ANOSIM tests to compare UniFrac distances for either all non-Proteobacteria OTUs and Pseudomonas OTUs only (Supplementary Table S3). Despite obtaining low R values, thus suggesting the weak clustering of communities, this additional analysis revealed that the effect of these factors   p values that reached statistical significance (p < 0.05) or were close to reaching significance (p < 0.1) are highlighted in bold for better visualization.
Considering the differences in relative abundance among the different bacterial groups of the duodenal microbiota (e.g., the higher Pseudomonas on GFD in CD and NCGS patients), we hypothesized that the beta diversity analyses for different bacterial populations may offer clues regarding the effect on different factors such as the diet and health status. Therefore, we used the Adonis and ANOSIM tests to compare UniFrac distances for either all non-Proteobacteria OTUs and Pseudomonas OTUs only (Supplementary Table S3). Despite obtaining low R values, thus suggesting the weak clustering of communities, this additional analysis revealed that the effect of these factors is different in distinct populations of microbes. For example, the effect of the BMI was stronger for Pseudomonas populations (Table S3).

Effect of GFD on Fecal Microbiota
Gluten proteins are not completely digested in the small intestine and several members of the fecal microbiota have the capacity to metabolize gluten [42]; therefore, the removal of gluten from the diet may also affect the distal gut microbiota. In this study, however, both diet and health status were not associated with differences in fecal bacterial richness or diversity (Table 3 and Supplementary Table S1). LEfSe analysis did not find any indication to suggest a difference in fecal microbial communities according to diet as the class or diet as a class and health status as a subclass. Interestingly, the LEfSe approach revealed a diverse group of microorganisms that were significantly enriched in each of the disease states when using health status as the main class ( Figure 9). The family Veillonellaceae, which was found to be lower in the feces of healthy subjects on GFD [21] and contains sulfite reducer members [43], was included in this group (higher in CD patients, Figure 9). The analysis of health status as class and diet as subclass revealed that Proteobacteria (in general, without an indication of a particular taxon within) was more abundant in CD patients. Beta-diversity analyses of UniFrac distances showed a significant grouping of fecal samples accordingly to BMIs and this relationship was also independent of low-abundant OTUs (Adonis test, Table 5 and Table S2).
Nutrients 2018, 10, x FOR PEER REVIEW 13 of 19 is different in distinct populations of microbes. For example, the effect of the BMI was stronger for Pseudomonas populations (Table S3).

Effect of GFD on Fecal Microbiota
Gluten proteins are not completely digested in the small intestine and several members of the fecal microbiota have the capacity to metabolize gluten [42]; therefore, the removal of gluten from the diet may also affect the distal gut microbiota. In this study, however, both diet and health status were not associated with differences in fecal bacterial richness or diversity (Table 3 and Supplementary  Table S1). LEfSe analysis did not find any indication to suggest a difference in fecal microbial communities according to diet as the class or diet as a class and health status as a subclass. Interestingly, the LEfSe approach revealed a diverse group of microorganisms that were significantly enriched in each of the disease states when using health status as the main class ( Figure 9). The family Veillonellaceae, which was found to be lower in the feces of healthy subjects on GFD [21] and contains sulfite reducer members [43], was included in this group (higher in CD patients, Figure 9). The analysis of health status as class and diet as subclass revealed that Proteobacteria (in general, without an indication of a particular taxon within) was more abundant in CD patients. Beta-diversity analyses of UniFrac distances showed a significant grouping of fecal samples accordingly to BMIs and this relationship was also independent of low-abundant OTUs (Adonis test, Table 5 and Table S2).

Effect of GFD on the Predicted Functional Profile
The closed OTU picking approach yielded a total of 4958 OTUs in biopsies and fecal samples. PICRUSt revealed no significant difference in the predicted functional profile of duodenal or fecal microbiota accordingly to diet or health status. Interesting results were found (for fecal samples only) when analyzing the group factor (six groups, control, CD and NCGS patients before and on GFD). For example, the proportions of genes related to the propanoate metabolism were higher in CD patients on a GFD (see "3.2 Predicted functional profile" in Supplementary Information) but caution must be exerted because of the low sample size in each subgroup of patients.

Discussion
Increasing evidence suggests a role of the gut microbiota in the onset and clinical development of GRDs but this phenomenon has been mostly studied in Europe. This study sheds light for the first time into the complex host-microbiota interactions in control subjects and patients with CD and NCGS from México. Additionally, this study offers relevant clues regarding the potential effect of GFD on health and gut microbiota.

Effect of GFD on the Predicted Functional Profile
The closed OTU picking approach yielded a total of 4958 OTUs in biopsies and fecal samples. PICRUSt revealed no significant difference in the predicted functional profile of duodenal or fecal microbiota accordingly to diet or health status. Interesting results were found (for fecal samples only) when analyzing the group factor (six groups, control, CD and NCGS patients before and on GFD). For example, the proportions of genes related to the propanoate metabolism were higher in CD patients on a GFD (see "3.2 Predicted functional profile" in Supplementary Information) but caution must be exerted because of the low sample size in each subgroup of patients.

Discussion
Increasing evidence suggests a role of the gut microbiota in the onset and clinical development of GRDs but this phenomenon has been mostly studied in Europe. This study sheds light for the first time into the complex host-microbiota interactions in control subjects and patients with CD and NCGS from México. Additionally, this study offers relevant clues regarding the potential effect of GFD on health and gut microbiota.
The gluten metabolism is an interesting physiological phenomenon and growing evidence suggests a strong involvement of the gut microbiota [44]. However, each individual carries a highly specific group of microorganisms even at the strain level [1], and therefore such an involvement must be highly individualized. More importantly, the response of these unique communities to environmental factors (e.g., dietary changes, antibiotic administration) is also unique and may never return to the exact same baseline state before the challenge [45]. Finally, the region where the individuals live is an important factor, in fact, one study showed important interactions between the patients' geographical location and the clinical and microbiological manifestations of inflammatory bowel disease [16]. In this study, for example, our results are unlikely to apply to patients with GRDs from other cities, even within the same state of Veracruz.
From a clinical perspective, four weeks on GFD often improves symptoms and the quality of life in patients with CD or NCGS and this paper shows that this period of time was also enough to change the gut microbiota in our group of subjects, for example, duodenal Pseudomonas in NCGS patients. In contrast, Tjellström et al. [46] showed that fecal short-chain fatty acids output (a direct result of microbial activity) in CD patients with more than one year on GFD was significantly different compared to the output in CD patients with less than one year on GFD and CD patients at the presentation, thus suggesting that a long period of time on GFD may be necessary to fully re-establish the functioning of the gut microbial ecosystem in some patients. It has also been shown that a subgroup of patients does not respond positively even while adhering to a strict GFD and that these patients seem to harbor distinctive microbiota [47]. Here we showed that each individual carries a highly specific gut microbial composition, that the microbiota is different between healthy subjects and people with GRDs, and that this microbiota can experience variation due to the removal of gluten. It is important to note that this change also varies widely among individuals (the most significant and consistent change was associated with duodenal Pseudomonas in NCGS patients but every individual showed a unique increase or decrease in the abundance of these and other microorganisms).
The (unexpected) finding of higher abundance of Pseudomonas in some patients during GFD deserves special attention. For instance, whether the increase in Pseudomonas is beneficial or not to the integrity of the duodenal mucosa is uncertain. Clinicians often associate Pseudomonas with diseases because of the pathogenic nature of some strains of P. aeruginosa and other species. However, Pseudomonas is a highly heterogenic bacterial genus that includes thousands of non-pathogenic, highly divergent strains inhabiting a wide variety of environments [48]. Unfortunately, very few studies have paid attention to native gut-associated Pseudomonas [49][50][51][52]. The finding that a GFD is associated with a higher abundance of Pseudomonas in the duodenum could be explained using at least two hypotheses. First, gluten may lead to a given immunological status in the mucosa that interferes negatively with the presence of autochthonous Pseudomonas, thus explaining the lower abundance at baseline. Second, some members of Pseudomonas may act as a protective microbe and its low abundance may prompt a more sensitive state to dietary allergens. This is supported by the relatively lower abundance of Pseudomonas in CD and NCGS patients before the GFD (Figure 7).
The possibility that some members of Pseudomonas can act as protective agents suggest that some strains of Pseudomonas may even be considered as probiotics for patients with some GRDs. Interestingly, Gao et al. [53] showed that Pseudomonas and other bacteria were reduced in cancerous tissues compared to adjacent non-cancerous tissues, thus suggesting a protective role in the gut mucosa. Wei et al. [54] identified an interesting aciduric gluten-degrading enzyme from P. aeruginosa with a therapeutic potential for CD; yet this does not explain whether GFD would lead to a higher or lower abundance of gluten-degrading Pseudomonas (we reasoned that gluten-degrading Pseudomonas populations would grow preferentially only if gluten consumption offers a selective advantage). One study showed higher abundances of Pseudomonas in the duodenum of adult CD patients on a GFD compared to controls but this finding was not discussed at all [55]. This current study also suggests that other non-Pseudomonas Proteobacteria (e.g., Stenophomonas) deserve attention in terms of gluten degradation and gut health.
This study also shows that the health status in terms of gluten sensitivity may be related to differences in the distal digestive microbiota. For example, this study showed a higher abundance of Ruminococcaceae in the fecal microbiota of NCGS patients. Additionally, Veillonellaceae, a pro-inflammatory taxon that has been shown to be increased in patients with inflammatory bowel disease and inflammatory bowel syndrome [56][57][58], was shown to be increased in fecal samples from CD patients. This adds valuable information to a growing literature showing that the distal microbiota is also worth looking at in gluten-related disorders [59].
This study has limitations that are relevant to future studies. First, this and other studies lack a large enough sample size to generalize phenomena and even with bigger samples sizes the results cannot be extrapolated from one population to others [17]. Second, gluten-free diets vary widely around the world and these may or may not lead to a microbial state more similar to healthy controls [13]. Third, 16S sequencing does not inform about the microbe-immune system interaction at the cell level. In this regard, De Palma et al. [60] showed interesting differences in IgA-coated fecal bacteria in treated and untreated CD patients, thus suggesting that a simple molecular characterization of microbes is not enough to fully capture the complex relationship. Fourth, one study showed that serum concentrations of short-chain fatty acids were similar in the control and CD patients; however, the authors found an interesting difference between genders [61]. This is particularly important because the reasons explaining the differences between genders with regards to the clinical presentation and severity of GRDs and other autoimmune disorders have not been fully clarified. One hypothesis suggests that infections can induce autoimmune diseases [62]. Finally, we only looked at the bacterial microbiota here, yet non-bacterial organisms (e.g., yeasts) may play a role in these disorders [63].

Conclusions
In summary, this study generates valuable preliminary data about the relationship between the gut microbiota and gluten-related disorders in Mexican people. Interestingly, the four-week consumption of GFD was associated with an increased abundance of Pseudomonas in duodenal biopsies of patients with these disorders, particularly in NCGS patients. This change was noticed despite a general lack of differences in richness or diversity. Pseudomonas comprises strains with gluten-degrading capabilities that deserves more attention. It is our hope that these results can contribute to starting to visualize alternatives for the more effective treatment of afflicted patients in our area.