The Gut Microbiota Profile According to Glycemic Control in Type 1 Diabetes Patients Treated with Personal Insulin Pumps

Recently, several studies explored associations between type 1 diabetes (T1DM) and microbiota. The aim of our study was to assess the colonic microbiota structure according to the metabolic control in T1DM patients treated with insulin pumps. We studied 89 T1DM patients (50.6% women) at the median age of 25 (IQR, 22–29) years. Pielou’s evenness (p = 0.02), and Shannon’s (p = 0.04) and Simpson’s diversity indexes (p = 0.01), were higher in patients with glycosylated hemoglobin (HbA1c) ≥ 53 mmol/mol (7%). There were no differences in beta diversity between groups. A linear discriminant analysis effect size (LEfSe) algorithm showed that one family (Ruminococcaceae) was enriched in patients with HbA1c < 53 mmol/mol, whereas one family (Streptococcaceae) and four species (Ruminococcus torques, unclassified species of Lactococcus, Eubacteroim dolichum, and Coprobacillus cateniformis) were enriched in patients with HbA1c ≥ 53 mmol/mol. We found that at class level, the following pathways according to Kyoto Encyclopedia of Genes and Genomes were enriched in patients with HbA1c < 53 mmol/mol: bacterial motility proteins, secretion system, bacterial secretion system, ribosome biogenesis, translation proteins, and lipid biosynthesis, whereas in patients with HbA1c ≥ 53 mmol/mol, the galactose metabolism, oxidative phosphorylation, phosphotransferase system, fructose, and mannose metabolism were enriched. Observed differences in alpha diversity, metabolic pathways, and associations between bacteria and HbA1c in colonic flora need further investigation.


Introduction
Type 1 diabetes mellitus (T1DM) is an autoimmune disease that results from the destruction of the pancreatic beta cells producing insulin [1]. Because of that, patients with T1DM require lifelong insulin therapy that mimics the basal and mealtime release of insulin by the pancreas [1]. Intensive insulin therapy involves multiple insulin injections (MDI) or continuous subcutaneous insulin infusion (CSII) [1].
The pathogenesis of T1DM is not fully understood [1]. The genetic-specific HLA predisposition, and epigenetic and environmental factors such as diet or viral infections, are reported to contribute to the development of the disease [1,2]. It has been postulated that one of the potential factors contributing to the development of T1DM or triggering the disease could be dysbiosis of interstitial bacterial flora [3,4]. The underlying mechanisms are unknown. It has been proposed, for instance, that alternations in gut microbiota can cause a decrease in bacteria producing butyrate, a short chain fatty acid (SCFA), and as a consequence, can cause an increase in permeability and autoimmunity for T1DM [5]. Recently, several studies have shown the relationships between beta-cell autoimmunity, occurrence of T1DM, and microbial dysbiosis [6][7][8][9][10][11][12][13][14].
The potential association between gut dysbiosis and T1DM may be, however, not only casual, but it could potentially affect glycemic control in the already developed disease. Metabolic endotoxemia is described as a condition without acute infection when there is an elevated lipopolysaccharide (LPS) level, a main component of the Gram negative bacteria membrane [15]. Increased fasting and postprandial LPS levels are observed in both T1DM and T2DM [15]. Some data, mostly based on type 2 diabetes mellitus (T2DM), indicate that derangements in the composition of the microbiota may play a role in the development of "metabolic endotoxemia" leading to hyperglycemia based on immune independent mechanisms [16,17].
A recent meta-analysis has shown that the modification of gut microbiome following probiotics usage may improve glycemic status in T2DM [18]. Since those mechanisms and interventions could be also applicable to T1DM, one could speculate that there could also be a relationship between gut microbiota structure and glycemic status in T1DM individuals.
It has been reported that probiotics/prebiotics can decrease glucose level, increase SCFAs, diminish the inflammatory pathway, and reduced systematic endotoxemia in T1DM animal studies [2]. It is hypothesized that prebiotics via changing intestinal microbiota could lead to better glycemic control in T1DM patients [19].
Given reported improvement in glucose control after changes in microbiota in T1DM animal studies, we hypothesized that colonic flora differs according to the metabolic control in T1DM patients.
The aim of our study was to compare gut microbiota structure and their metabolic pathways between T1DM patients meeting the current criteria of glycemic control [20,21] with glycosylated hemoglobin (HbA1c) below 53 mmol/mol (7%) and patients with HbA1c equal to or greater than 53 mmol/mol (7%).

Study Setting and Eligibility
The study was performed from 2016 to 2018. Ninety-four consecutive patients with T1DM treated with insulin pumps at our Outpatients Diabetic Clinic were enrolled.
Since the main goal of the study was to assess the relationship between microbiota structure and HbA1c per se, we aimed to perform our analysis on as homogenous a group as possible.
Patients with T1DM who use an insulin pump and confirmed readiness to cooperate with the research center were eligible for the study.
The exclusion criteria were chosen carefully and included clinical conditions and patients' behavior, which could have a significant impact on microbiota status irrespective of glycemic control. These included the following: • confirmed infection of the gastrointestinal tract or using probiotics or taking antibiotics for up to 30 days before delivery of a stool sample, • chronic inflammatory bowel disease of unknown etiology, active cancer (especially of the gastrointestinal tract), • immunodeficiency, • the presence of advanced late complications of diabetes, and • low and very high carbohydrate consumption defined as below 100 or over 400 g of carbohydrates verified by questionnaire and insulin pump downloads.
The study received approval from the Jagiellonian University Ethics Committee (number KBET/256/B/2014, date 27.11.2014). All participants provided written informed consent in accordance with the Declaration of Helsinki.

Study Investigations
We obtained reports from personal insulin pumps (14 days). The reports provided the information about average glucose level along with standard deviation (SD); basal, bolus, and total daily insulin dose; the amount of glucose measurements per day; and the average amount of carbohydrates eaten per day. The patients were asked to complete the survey regarding eating habits. HbA1c was assessed (high-performance liquid chromatography using the Variant II Turbo analyzer, Hercules, CA, USA, or HbA1c, Quo Lab HbA1c, EKF Diagnostics) and stool samples were collected from all patients.

DNA Isolation and 16S Metagenomic Sequencing
Bacterial DNA was isolated using Genomic Mini AX Stool Spin (A&A Biotechnology, Gdynia, Poland), with several modifications, as described in our previous publication [22].
Libraries were prepared according to the Illumina 16S Metagenomic Sequencing Library Preparation protocol. After that, the libraries at a concentration 10 pM with 20% PhiX spike-in control were sequenced on Illumina MiSeq (Illumina, Inc., San Diego, CA, USA) using the V3 sequencing kit (300 bp paired-end reads) [22].

Sequencing Data Analysis
Samples were processed and analyzed using the Quantitative Insights Into Microbial Ecology 2 (QIIME2, version 2019.7) [23] custom pipeline. Briefly, the quality of demultiplexed paired-end reads from MiSeq (2 × 300 bp) was evaluated and the reads were trimmed to remove primers and poor-quality bases with cutadapt [24]. Trimmed pair-end sequences were denoized and merged with DADA2 [25] to generate amplicon sequence variants (ASVs). Next, reference-based chimera filtering was queried against the reference database (Greengenes version 13.8) at 99% similarity with vsearch [26].
Low-abundance ASVs were eliminated when they appeared in less than three samples or the number of counts across all samples was <5. The SATé-enabled phylogenetic placement (SEPP) algorithm was used to build the tree for phylogenetic diversity computation [27]. QIIME2 diversity core-metrics-phylogenetic analyses were used to compute alpha and beta diversity values, and rarefaction curve analysis was used to estimate the completeness of microbial community sampling. We also computed default alpha and beta diversity metrics and generated principal coordinates analysis (PCoA) plots for each of the beta diversity metrics using the EMPeror [28]. The alpha and beta diversity indices of the groups were compared using QIIME2 longitudinal pairwise-differences plugin using the t-test. Correlations with alpha and beta diversity indices were calculated with the QIIME2 plugins using Spearman correlation. Generated ASVs were assigned to taxonomy using a naive Bayes classifier [29] that was pre-trained on the v3-v4 rRNA regions in Greengenes version 13.8, at 99% similarity. The bacterial composition was analyzed at phylum, class, order, family, genus, and species levels. Differential abundance between groups at each taxonomic level was tested using analysis of the composition of microbiomes (ANCOM) [30]. Additionally, we used 'MaAsLin2 package [31] in R statistical environment to verify the multivariable association between clinical metadata and relative abundance of taxa at different taxonomical levels. The regression coefficient (coef) in MaAsLin2 is the effect size, which represents the rate of change in abundance of taxa per 1 score of measurements. Both, ANCOM and MaAsLin2 control for FDR was at level 0.05. To find taxa and pathways that explain differences between our group, we used LEfSe (Linear Discriminant Analysis Effect Size) [32], which is an algorithm for high-dimensional biomarker discovery and a tool to identify genomic features, such as taxa or pathways. The LEfSe was used with the default parameters (p < 0.05 and Linear Discriminant Analysis [LDA] score >2.0). First, we performed differential abundance on species taxonomical level, and then we predicted KEGG (Kyoto Encyclopedia of Genes and Genomes) [33] based on functional profiles using PICRUSt (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) software v.1.1.4 [34]. Since PICRUSTt requires an operational taxonomic unit (OTU)_ table 'closed-reference' picked against Greengenes database, we clustered our ASVs sequences with Greengenes 13.8 database, at 99% similarity, using the vsearch cluster-features-closed-reference [26] plugin in QIIME2 [23].

Statistical Analysis
Continuous variables are expressed as median (interquartile range (IQR)) or mean (± SD). Categorical variables are presented as count (percentages). The normality of continuous variables was verified using the Shapiro-Wilk test. The equality of variances was checked by the Levene's test. Groups were compared using Student's t-test, Welch t-test, Mann-Whitney U test, and Fisher's exact test, appropriately. Calculations were performed using Statistica 13 software (StatSoft Inc., Tulsa, OK, USA) and package R ver 3.6.2 [35].

Study Population
From all 94 stool samples, 5 samples were excluded due to low read count after DADA2 denoising, chimera removal, and filtering. Therefore, the final study population consisted of 89 patients (50.6% women). The median age of the cohort was 25 (22)(23)(24)(25)(26)(27)(28)(29) years. The patients were divided into two groups according to the HbA1c level: patients with HbA1c below 53 mmol/mol (7%) (n = 43) and patients with HbA1c equal to or greater than 53 mmol/mol (7%) (n = 46). There were no differences in age, body mass index (BMI), and duration of diabetes between both groups (Table 1). Patients with HbA1c below 53 mmol/mol (7%) more often measured their glucose level and declared higher consumptions of carbohydrates per day. There were no differences in daily insulin dose or in dose of insulin per kilogram body weight between groups (Table 1).
Twenty patients suffered from hypothyroidism, two from celiac disease, two from vitiligo, one patient had autoimmune hepatitis, one patient had Addison-Biermer anemia, and four patients had hypertension ( Table 1). The groups did not differ in hypothyroidism diagnosis. The groups did not differ in general eating habits (Table 2).

16S rRNA Sequencing Analysis
After DADA2 denoising, chimera removal, and filtering of 94 patients, 5 samples were excluded due to low read count. Sequencing analysis for the 89 patients in the final cohort provided a median of 14,865 (IQR, 11,067-17,897). The best sample contained 33,933 reads, while the worst contained 6267 reads. DADA2 pipeline resulted in 480, features (ASVs), with total frequency of 1,340,777 features. Based on rarefaction curves, the diversity analysis was calculated with a minimum number of 6267 sequences per samples, which corresponded to the minimum frequency.
There were the following baseline bacterial profiles with an abundance of >1% at the phylum level in patients with HbA1c below 53 mmol/mol (7%) and patients with HbA1c  HbA1c level (as a continuous variable) was not associated with beta diversity indices, Bray-Curtis dissimilarity (R = 0.03, p = 0.48), Jaccard distance (R = 0.04, p = 0.39), or unweighted (R = 0.05, p = 0.28) and weighted UniFrac (R = 0.05, p = 0.32). Average glucose level, standard deviation of average glucose level, and insulin per kilogram body weight were also not correlated with any of beta diversity indices (data not shown).

Differential Abundance of Microbial Taxa
There were no differences in ANCOM results for phylum, class, order, family, genus, or species between both groups. When adjusted for sex, age, BMI, diabetes duration in years, and total daily insulin in units, still no differences on tested taxonomical level were found between both HbA1c groups. Then, we decided to use MaAsLin2 to verify the multivariable association between clinical metadata, such as sex, age, BMI, diabetes duration in years, total daily insulin, daily carbs, and HbA1c in a linear model. We performed anal- Core microbiome analysis at species level showed that both subgroups shared 129 taxa, three were unique to patients with HbA1c < 53 mmol/mol (7%) and five were unique to subjects with HbA1c ≥ 53 mmol/mol (7%). Detailed taxonomic information about overlaps between groups is available in supplementary data (Supplement 1).
The groups did not differ in Firmicutes/Bacteroidia (F/B) ratio (p = 0.59). We observed borderline significance in F/B ratio between men and women (p = 0.08), although there were no differences in the F/B ratio between HbA1c groups in females (p = 0.39), and in men (p = 0.98). The F/B ratio was not correlated with HbA1c (p = 0.33) or mean glucose level (p = 0.28). There was no association between F/B ratio and dose of insulin per kilogram in patients adjusted for HbA1c (p = 0.77) or mean glucose level (p = 0.78).

Differential Abundance of Microbial Taxa
There were no differences in ANCOM results for phylum, class, order, family, genus, or species between both groups. When adjusted for sex, age, BMI, diabetes duration in years, and total daily insulin in units, still no differences on tested taxonomical level were found between both HbA1c groups. Then, we decided to use MaAsLin2 to verify the multivariable association between clinical metadata, such as sex, age, BMI, diabetes duration in years, total daily insulin, daily carbs, and HbA1c in a linear model. We performed analysis on phylum, class, order, family, genus, and species level, but did not find any significant association (all FDR > 0.05) between tested clinical data and relative abundance of gut microbiota in T1DM patients. Since there were no differences in differential abundance, the changes in the microbial community may be more discrete; therefore, we performed additional analysis with the LEfSe algorithm to determine the bacteria that might be associated with differences in HbA1c (Figure 4). In patients with HbA1c < 53 mmol/mol (7%), one family (Ruminococcaceae) was enriched. One family (Streptococcaceae) and four species (Ruminococcus torques, unclassified species of Lactococcus, Eubacteroim dolichum, and Coprobacillus cateniformis) were enriched in patients with HbA1c ≥ 53 mmol/mol (7%).
Microorganisms 2021, 9, x FOR PEER REVIEW 9 of 14 gut microbiota in T1DM patients. Since there were no differences in differential abundance, the changes in the microbial community may be more discrete; therefore, we performed additional analysis with the LEfSe algorithm to determine the bacteria that might be associated with differences in HbA1c (Figure 4). In patients with HbA1c < 53 mmol/mol (7%), one family (Ruminococcaceae) was enriched. One family (Streptococcaceae) and four species (Ruminococcus torques, unclassified species of Lactococcus, Eubacteroim dolichum, and Coprobacillus cateniformis) were enriched in patients with HbA1c ≥ 53 mmol/mol (7%). Additionally, in MaAsLin2 HbA1c (as a continuous variable), average glucose level or standard deviation of average glucose level were not correlated with any OTU at phylum, class, order, family, genus, or species level. The dose of insulin per kilogram body weight was positively correlated with Bacteroides uniformis (coef = 0.008, FDR = 0.037), and we found that insulin per kilogram body weight after adjustment for age, sex, and standard deviation of average glucose level was negatively correlated with phyla Firmicutes (coef= −0.017, FDR = 0.043).

Figure 4.
The linear discriminant analysis effect size (LEfSe) method showed the significantly different taxa of patients with type 1 diabetes (T1DM) associated with glycosylated hemoglobin (HbA1c) < 53 mmol/mol (7%) (green) and ≥ 53 mmol/mol (7%) (red). The taxa with significantly different abundances among the groups were presented as barplot (a) with linear discriminant analysis (LDA) score, and on the clarogram (b). From inside to outside, the dots denote the kingdom, phylum, class, order, family, genus, and species.
Additionally, in MaAsLin2 HbA1c (as a continuous variable), average glucose level or standard deviation of average glucose level were not correlated with any OTU at phylum, class, order, family, genus, or species level. The dose of insulin per kilogram body weight was positively correlated with Bacteroides uniformis (coef = 0.008, FDR = 0.037), and we found that insulin per kilogram body weight after adjustment for age, sex, and standard deviation of average glucose level was negatively correlated with phyla Firmicutes (coef= −0.017, FDR = 0.043).

Functional Profiles of Gut Microbiota
Using LEfSe, we also explored KEGG pathways that were differentially enriched between HbA1c < 53 mmol/mol (7%) and HbA1c ≥ 53 mmol/mol (7%). At phylum level, we observed only four pathways, and interestingly, the group with higher HbA1c was enriched in the pathway associated with Energy Metabolism and Carbohydrate Metabolism (Figure 5a). At class level, we observed 10 metabolic pathways differentially enriched between patients in our groups (Figure 5b). In patients with HbA1c < 53 mmol/mol (7%) the following pathways were enriched: bacterial motility proteins, secretion system, bacterial secretion system, ribosome biogenesis, and translation proteins lipid biosynthesis, whereas in patients with HbA1c ≥ 53 mmol/mol (7%), the galactose metabolism, oxidative phosphorylation, phosphotransferase system, fructose, and mannose metabolism were enriched. The linear discriminant analysis effect size (LEfSe) method showed the significantly different taxa of patients with type 1 diabetes (T1DM) associated with glycosylated hemoglobin (HbA1c) < 53 mmol/mol (7%) (green) and ≥53 mmol/mol (7%) (red). The taxa with significantly different abundances among the groups were presented as barplot (a) with linear discriminant analysis (LDA) score, and on the clarogram (b). From inside to outside, the dots denote the kingdom, phylum, class, order, family, genus, and species.

Functional Profiles of Gut Microbiota
Using LEfSe, we also explored KEGG pathways that were differentially enriched between HbA1c < 53 mmol/mol (7%) and HbA1c ≥ 53 mmol/mol (7%). At phylum level, we observed only four pathways, and interestingly, the group with higher HbA1c was enriched in the pathway associated with Energy Metabolism and Carbohydrate Metabolism (Figure 5a). At class level, we observed 10 metabolic pathways differentially enriched between patients in our groups (Figure 5b). In patients with HbA1c < 53 mmol/mol (7%) the following pathways were enriched: bacterial motility proteins, secretion system, bacterial secretion system, ribosome biogenesis, and translation proteins lipid biosynthesis, whereas in patients with HbA1c ≥ 53 mmol/mol (7%), the galactose metabolism, oxidative phosphorylation, phosphotransferase system, fructose, and mannose metabolism were enriched.

Discussion
To our knowledge, we report for the first time differentially enriched KEGG metabolic pathways according to HbA1c among a group of patients with T1DM in the Polish cohort.
Surprisingly, in our study, we noticed slightly greater alpha diversity in T1DM patients with HbA1c ≥ 53 mmol/mol (7%) compared to patients with T1DM and <53 mmol/mol (7%). We also found a positive correlation between alpha diversity and HbA1c in the T1DM cohort. Moreover, we observed some significantly different taxa using the LEfSe algorithm between groups. Our study demonstrates differences in colon microbiota regarding the level of metabolic control of diabetes in T1DM patients.
First of all, we would like to underline that our cohort of young adults with T1DM was well-characterized and homogenous. All patients were treated with personal insulin pumps, with downloading of the devices allowed to obtain precise information concerning insulin dosing. Study participants were free from advanced complications of diabetes, quite homogenous with regard to diet, and the majority of them were free from comorbidities that could affect study results (e.g., celiac disease, n = 2). This minimizes the effect of confounding factors on the relationship of single variable, HbA1c, and the microbiota profile.
In our study, we did not observe differences in ANCOM results for phylum, class, order, family, genus, or species between individuals with HbA1c level below 53 mmol/mol (7%) and individuals with HbA1c equal to or greater than 53 mmol/mol (7%).

Discussion
To our knowledge, we report for the first time differentially enriched KEGG metabolic pathways according to HbA1c among a group of patients with T1DM in the Polish cohort.
Surprisingly, in our study, we noticed slightly greater alpha diversity in T1DM patients with HbA1c ≥ 53 mmol/mol (7%) compared to patients with T1DM and <53 mmol/mol (7%). We also found a positive correlation between alpha diversity and HbA1c in the T1DM cohort. Moreover, we observed some significantly different taxa using the LEfSe algorithm between groups. Our study demonstrates differences in colon microbiota regarding the level of metabolic control of diabetes in T1DM patients.
First of all, we would like to underline that our cohort of young adults with T1DM was well-characterized and homogenous. All patients were treated with personal insulin pumps, with downloading of the devices allowed to obtain precise information concerning insulin dosing. Study participants were free from advanced complications of diabetes, quite homogenous with regard to diet, and the majority of them were free from comorbidities that could affect study results (e.g., celiac disease, n = 2). This minimizes the effect of confounding factors on the relationship of single variable, HbA1c, and the microbiota profile.
In our study, we did not observe differences in ANCOM results for phylum, class, order, family, genus, or species between individuals with HbA1c level below 53 mmol/mol (7%) and individuals with HbA1c equal to or greater than 53 mmol/mol (7%). However, we noticed differences between groups in some alpha diversity indices. We observed higher alpha diversity in patients with higher HbA1c and a positive correlation between HbA1c and alpha diversity. A previous study showed a decrease in alpha diversity in patients after the seroconversion and before T1DM development, and the decrease in alpha diversity was also observed between T1DM patients in comparison to controls [36]. Also, patients with newly diagnosed T2DM had decreased diversity compared to healthy control group [37]. There are inconclusive results whether obesity is associated with changes in alpha diversity [38], as an increased diversity in obese patients compared to controls was also reported [39].
We did not find differences in F/B ratio between our two groups of patients. In some previous cohort studies, it has been demonstrated that patients with T1DM had a decreased F/B ratio [40,41]. However, we found that insulin per kilogram body weight after adjustment for age, sex, and standard deviation of average glucose level was negatively correlated with phyla Firmicutes. That could indicate the potential role of hipo/hyperinsulinemia itself on microbiota composition.
We did not notice any correlation between HbA1c level and single microbal OTU. We found that the family Ruminococcaceae was enriched in a group of patients with HbA1c < 53 mmol/mol (7%) based on LEfSe algorithm, which could be partially in agreement with the a study of Hang et al., who observed negative correlation between HbA1c and both Ruminococcaceae and its genus Faecalibacterium [40]. Previously, Qi et al. observed a positive correlation between relative abundance of Blautia and HbA1c among newly diagnosed children with T1DM and healthy controls [11], and Fassatoui et al. noticed a negative correlation between HbA1c and amount of Akkermansia muciniphia [41].
As we did not perform the shotgun sequencing, we predicted KEGG pathways based on functional profiles using PICRUSt. We found that microbiota profile in patients with HbA1c < 53 mmol/mol (7%) were characterized by more abundant metabolic pathways related to bacterial motility proteins, secretion system, bacterial secretion system, ribosome biogenesis, translation proteins, and lipid biosynthesis proteins, whereas in patients with HbA1c ≥ 53 mmol/mol (7%), they were related to the galactose metabolism, oxidative phosphorylation, phosphotransferase system, fructose, and mannose metabolism. The lipid biosynthesis proteins include the acetyl-coenzyme A (CoA) pathway, which is the main carbohydrate-driven pathway [42]. The conversion of acetyl-CoA to butyrate is essential to bacteria growth [42]. Butyrate is the main source of energy for colon cells [42][43][44]. In our study, we only determined that the lipid biosynthesis pathway was enriched in gut microbiota from patients with HbA1c < 53 mmol/mol (7%) compared to individuals with HbA1c ≥ 53 mmol/mol (7%). Our finding that family Ruminococcaceae was enriched in a group of patients with HbA1c < 53 mmol/mol (7%) could support the thesis that this pathway could be enriched in this group as some butyrate-producers belong to Ruminococcaceae.
It should be highlighted that associations should be considered as hypothesis-generating findings that should be investigated in further studies. Our observations indicate the importance of exploring the metabolic pathways of gut microbiota according to the metabolic control, and the importance of further studies on associations between metabolic control and microbiota in patients with T1DM.
One could speculate that increased glucose level favors bacteria with some metabolic pathways are responsible for carbohydrates metabolism. Also, poor metabolic control can play a role in shifting for some metabolic pathways. There may also be, independently, a two-way relationship-worse metabolic control that may change microbiota structure and that may hinder a return to better metabolic control, or that microbiota dysbiosis can worsen the metabolic control in patients with T1DM. Further studies are needed to explore it, because we were not able to answer this question fully as it was not the aim of our study.
Some study limitations should be acknowledged. The study population was relatively small. We did not obtain the information about the exact food eaten two weeks before the sample collection; however, the patients completed a survey on commonly consumed products (data not shown). We collected only one stool sample from every patient in a single point time; therefore, we cannot exclude casual changes in time. We did not obtain information about the patients' antibodies profile at the time of diagnosis and the time of the study. We did not assess HLA alleles in patients. We performed 16S sequencing; however, shotgun sequencing could give information about viruses and fungi, as well as direct information about which metabolic pathways are encoded. In our study, we could only predict KEGG pathways based on functional profiles using PICRUSt. We also had a lack of a healthy control group, which could have made a study more comprehensive. It would also be of value to investigate this issue in the population of T2DM since gut microbiota might modulate glucose homeostasis via different mechanisms in T1DM and T2DM.

Conclusions
We observed differences in metabolic pathways and alpha diversity in colonic flora according to HbA1c among patients with T1DM. Further studies should be performed to assess the association between metabolic control and microbiota and their metabolic pathways, also regarding autoimmunity profile and HLA alleles, as well as the possible association between the microbiota profile and insulin resistance.