Metagenomic Shotgun Sequencing Reveals Specific Human Gut Microbiota Associated with Insulin Resistance and Body Fat Distribution in Saudi Women

(1) Background: Gut microbiota dysbiosis may lead to diseases such as insulin resistance and obesity. We aimed to investigate the relationship between insulin resistance, body fat distribution, and gut microbiota composition. (2) Methods: The present study included 92 Saudi women (18–25 years) with obesity (body mass index (BMI) ≥ 30 kg/m2, n = 44) and with normal weight (BMI 18.50–24.99 kg/m2, n = 48). Body composition indices, biochemical data, and stool samples were collected. The whole-genome shotgun sequencing technique was used to analyze the gut microbiota. Participants were divided into subgroups stratified by the homeostatic model assessment for insulin resistance (HOMA-IR) and other adiposity indices. (3) Results: HOMA-IR was inversely correlated with Actinobacteria (r = −0.31, p = 0.003), fasting blood glucose was inversely correlated with Bifidobacterium kashiwanohense (r = −0.22, p = 0.03), and insulin was inversely correlated with Bifidobacterium adolescentis (r = −0.22, p = 0.04). There were significant differences in α- and β-diversities in those with high HOMA-IR and waist–hip ratio (WHR) compared to low HOMA-IR and WHR (p = 0.02, 0.03, respectively). (4) Conclusions: Our findings highlight the relationship between specific gut microbiota at different taxonomic levels and measures of glycemic control in Saudi Arabian women. Future studies are required to determine the role of the identified strains in the development of insulin resistance.


Introduction
Globally, the prevalence of type 2 diabetes mellitus (T2DM) mellitus is increasing at an alarming rate [1]. This is largely triggered by the rising rates of insulin resistance [2] and the obesity pandemic [3]. As of 2017, T2DM has affected an estimated 462 million people around the world and has led to over 1 million deaths per year, with rates expected to rise [1]. A systematic review and meta-analysis of data from 2000 to 2018 reported a combined prevalence of T2DM of about 15% in Middle Eastern adults [4]. Saudi Arabia is no exception to this crisis, with an estimated national prevalence of 18% among adults aged 20-79 years in 2019 [5]. Regional [6] and nationwide [7,8] studies have also revealed a higher prevalence of T2DM among women than among men. Concurrently, there is a marked increase in the prevalence of obesity globally, with 13% of adults suffering from obesity in 2016 [9]. In the Middle East, in 2016, 50% and 24% of women were estimated to be overweight or obese, respectively [10,11]. The Saudi Arabian Health Interview Survey 2019 also confirmed this finding, revealing a higher prevalence of obesity among women than men (22% vs. 19%, respectively) [12].

Sample Population
Methods for this analytical, observational study have been previously published [37]. In summary, between January 2019 and March 2020, female college students aged >18 years, who were either of normal weight (body mass index (BMI); 18.5-24.9 kg/m 2 ) or obese (BMI ≥ 30 kg/m 2 ) were recruited via emails, the social media network, etc. Subsequently, out of 400 women, 290 were eligible to participate. Vitamins and antibiotic users (n = 88), those with chronic illnesses (n = 20), women who were pregnant (n = 2), those who refused to provide stool samples (n = 193), and participants who had stool samples with DNA concentration below the required level (n = 5), were all excluded. Consequently, 92 female students were included in this study and were invited to visit the study clinic for data collection. Each participant signed a consent form before enrolling in the study. The Institutional Ethics Committee at King Saud University approved this study (KSU-IRB #E-19-3625).

Anthropometric Measurements
Weight was measured twice, in light clothing and without shoes, to the nearest 0.1 kg, immediately after the participants visited the clinic. Height was also measured to the nearest 0.1 cm. BMI was calculated by dividing the weight in kg by the square of height in m 2 . The umbilicus, the lowest rib, and the narrowest waist area were utilized to measure the waist circumference following the World Health Organization standard [9]. At the great trochanter level, the hip circumference was measured using a non-stretchable tape with the legs closely apart. Both readings were recorded to the closest 0.5 cm. If the two measurements differed by >2 cm, a third measurement was made, and the mean of the two most comparative measurements was used. Dividing the mean waist circumference by the mean hip circumference yielded the waist-hip ratio (WHR) which was categorized as follows; low WHR < 0.83 and high WHR ≥ 0.83 [38]. Bioelectrical impedance analysis (BIA) was used to quantify body composition indices (770 BIA; InBody, Seoul, Republic of Korea) [39]. Body fat % (BF%) was considered low or high when BF% ≤ 35% or BF% > 35%, respectively.

Biochemical Measurements
Following an overnight fast of longer than 10 h, blood samples were collected into two 5 mL tubes. Serum was separated and centrifuged within 15 min and sent to the study lab. Fasting blood sugar levels were measured using a biochemical analyzer (Konelab, Espoo, Finland) and insulin levels were measured using a LIAISON XL analyzer (DiaSorin, Saluggia, Italy). The following formula, the homeostatic model assessment of the insulin resistance index (HOMA-IR) [40], was used:

Stool Analysis
Fecal samples were collected aseptically and stored at −80 • C. DNA was extracted from frozen stool and stored at −20 • C for further processing. QIAamp PowerFecal DNA Isolation Kit (Qiagen, Hilden, Germany) was used to extract DNA from 0.25 g aliquots of frozen stool. Following the manufacturer's protocol, it was then eluted in a 100 µL C6 elution buffer. A NanoDrop spectrophotometer (NanoDrop Technologies, Wilmington, USA) was used to measure purity (260/280 ratio) and concentration of the extracted DNA. The DNA libraries were prepared using the Illumina Nextera XT Library Preparation Kit (Illumina, Inc., San Diego, CA, USA), and samples were sequenced using an Illumina sequencer. For quantitative reading, a Qubit ® fluorimeter was used.

Characterization of Microbial Composition
Unrestricted genomic sequencing of all microorganisms present was used to determine the gut microbiota composition; hence, WGS analyses were used [41]. This technique allows for the determination of the gut composition at the major phyla level. Using the CosmosID bioinformatics platform (CosmosID, Inc., Rockville, MD, USA), the unassembled sequencing data were assessed for multi-kingdom microbiome analysis, profiling antibiotic resistance and virulence genes, and quantifying the relative abundance of organisms. With curated genome databases, high-performance data mining, and hundreds of millions of metagenomic sequence reads, this method [42][43][44][45] quickly divides the reads into distinct sequences that produce microbes. Also, the unassembled sequence reads were searched against the curated CosmosID antibiotic resistance and virulence-associated gene databases to identify the community resistome and virulome, referring to the collection of antibiotic resistance and virulence-associated genes in the microbiome, respectively.

Statistical Analyses
IBM SPSS Statistics for Windows (version 24; IBM Corp, Armonk, NY, USA) was used to perform all statistical analyses. The normality of all quantitative variables was tested before the analyses. To determine the impact of insulin resistance, participants were divided into two distinct categories: high and low HOMA-IR, using the median (1.85 mmol/L) as the cut-off point. The independent sample t-test was used to determine the statistical differences between continuous variables, and results are shown as means and standard deviations. Correlation coefficients between gut microbiota and glycemic control were calculated using the Pearson correlation coefficient. A p-value of 0.05 was used to demonstrate statistical significance.
In low and high HOMA-IR groups, heatmaps were created using the pheatmap R package [46]. The α-diversity was calculated using the CosmosID taxonomic analysis software (R software Vegan package, version 2.5-6; Oksanen et al., 2019) to measure the gut microbiota richness. The β-diversity was calculated to determine the gut composition from the species-level relative abundance matrices.
To identify the influence of different adiposity indices, the sample was further divided into subgroups stratified by HOMA-IR and low/high adiposity indices (WHR, BF%, and BMI). Then, αand β-diversity were calculated for all subgroups at the species level. Wilcoxon rank-sum tests were used to calculate the statistical difference between groups for α-diversity by Shannon index, using the ggsignif package for R [47]. For β-diversity, the nonparametric PERMANOVA analysis was used based on the Bray-Curtis distance using Vegan's function adonis2.

Descriptive Characteristics
The current study included a total of 92 participants from Saudi Arabia. Table 1 presents the descriptive traits of participants. Participants known to have high HOMA-IR had higher BMI, BF%, WHR, and WC compared to participants with low HOMA-IR (p-value < 0.0001). Fasting blood glucose (p-value = 0.04), insulin, total cholesterol (p-value = 0.01), LDL (p-value = 0.01), HDL (p-value = 0.05), and triglycerides (p-value < 0.001) were higher in those with high compared to low HOMA-IR. There were no significant differences in the relative abundance of gut microbiota between participants with low and high HOMA-IR ( Table 1). The heatmap analysis presents the composition of the gut microbiota of the participants with high and low HOMA-IR ( Figure 1). Participants with high HOMA-IR had higher Bacteroides vulgatus and Bacteroides sp. 3_1_40A and lower Alistipes putredini and Alistipes shahii compared to those with low HOMA-IR. Biomolecules 2023, 13, x FOR PEER REVIEW 6 of 16

α-and β-Diversity Analysis in Low and High HOMA-IR Groups
The α-diversity, reflected by the Shannon index, was higher in the group with high HOMA-IR compared to the group with low HOMA-IR (5.1 ± 0.5 vs. 5.0 ± 0.5); however, the difference was not statistically significant (p-value = 0.10) (Figure 2A).

α-and β-Diversity Analysis in Low and High HOMA-IR Groups
The α-diversity, reflected by the Shannon index, was higher in the group with high HOMA-IR compared to the group with low HOMA-IR (5.1 ± 0.5 vs. 5.0 ± 0.5); however, the difference was not statistically significant (p-value = 0.10) (Figure 2A).
The β-diversity using the paired Bray-Curtis similarity index indicated marginal significance between those with high compared to low HOMA-IR (p = 0.07) ( Figure 2B).

α-and β-Diversity Analysis for Subgroups Stratified by HOMA-IR and Different Adiposity Indices
The α-diversity was significant in the subgroup with high HOMA-IR and high WHR compared to the subgroup with low HOMA-IR and low WHR only (p = 0.02) ( Figure 3A). No statistically significant differences were observed between other subgroups. Similarly, β-diversity was significant for the subgroup with high HOMA-IR and high WHR compared to those with high HOMA-IR and low WHR (p = 0.03) and for those with high HOMA-IR and high WHR compared to those with low HOMA-IR and low WHR (p = 0.05) ( Figure 3B). The β-diversity using the paired Bray-Curtis similarity index indicated marginal significance between those with high compared to low HOMA-IR (p = 0.07) ( Figure 2B).

α-and β-Diversity Analysis for Subgroups Stratified by HOMA-IR and Different Adiposity Indices
The α-diversity was significant in the subgroup with high HOMA-IR and high WHR compared to the subgroup with low HOMA-IR and low WHR only (p = 0.02) ( Figure 3A). No statistically significant differences were observed between other subgroups. Similarly, β-diversity was significant for the subgroup with high HOMA-IR and high WHR compared to those with high HOMA-IR and low WHR (p = 0.03) and for those with high HOMA-IR and high WHR compared to those with low HOMA-IR and low WHR (p = 0.05) ( Figure 3B).

α-and β-Diversity Analysis for Subgroups Stratified by HOMA-IR and Different Adiposity Indices
The α-diversity was significant in the subgroup with high HOMA-IR and high WHR compared to the subgroup with low HOMA-IR and low WHR only (p = 0.02) ( Figure 3A). No statistically significant differences were observed between other subgroups. Similarly, β-diversity was significant for the subgroup with high HOMA-IR and high WHR compared to those with high HOMA-IR and low WHR (p = 0.03) and for those with high HOMA-IR and high WHR compared to those with low HOMA-IR and low WHR (p = 0.05) ( Figure 3B).  For other subgroups, α-diversity was significant between the subgroup with high HOMA-IR and high BF% compared to low HOMA-IR and high BF% (p = 0.03) ( Figure 4A). The β-diversity was significant for those with high HOMA-IR and high BF% compared to low HOMA-IR and low BF% (p = 0.04) ( Figure 4B). For other subgroups, α-diversity was significant between the subgroup with high HOMA-IR and high BF% compared to low HOMA-IR and high BF% (p = 0.03) ( Figure 4A). The β-diversity was significant for those with high HOMA-IR and high BF% compared to low HOMA-IR and low BF% (p = 0.04) ( Figure 4B).
The α-diversity was significant for the subgroup with obesity and high HOMA-IR in comparison with the subgroup with normal weight and low HOMA-IR (p = 0.03) ( Figure 5A). The β-diversity was significant for the subgroup with obesity and high HOMA-IR compared to the group with obesity and low HOMA-IR (p = 0.05) ( Figure 4B), and for the subgroup with obesity and high HOMA-IR in comparison to the subgroup with normal weight and low HOMA-IR (p = 0.01) ( Figure 5B). and low WHR: (A) α-diversity, (B) β-diversity.
For other subgroups, α-diversity was significant between the subgroup with high HOMA-IR and high BF% compared to low HOMA-IR and high BF% (p = 0.03) ( Figure 4A). The β-diversity was significant for those with high HOMA-IR and high BF% compared to low HOMA-IR and low BF% (p = 0.04) ( Figure 4B).  The α-diversity was significant for the subgroup with obesity and high HOMA-IR in comparison with the subgroup with normal weight and low HOMA-IR (p = 0.03) ( Figure  5A). The β-diversity was significant for the subgroup with obesity and high HOMA-IR compared to the group with obesity and low HOMA-IR (p = 0.05) ( Figure 4B), and for the subgroup with obesity and high HOMA-IR in comparison to the subgroup with normal weight and low HOMA-IR (p = 0.01) ( Figure 5B). The α-diversity was significant for the subgroup with obesity and high HOMA-IR in comparison with the subgroup with normal weight and low HOMA-IR (p = 0.03) ( Figure  5A). The β-diversity was significant for the subgroup with obesity and high HOMA-IR compared to the group with obesity and low HOMA-IR (p = 0.05) ( Figure 4B), and for the subgroup with obesity and high HOMA-IR in comparison to the subgroup with normal weight and low HOMA-IR (p = 0.01) ( Figure 5B).

Discussion
This study aimed to investigate the relationship between insulin resistance, body fat distribution, and gut microbiota in young women in Saudi Arabia using the WGS technique. Results show that specific gut microbiota at different taxonomic levels in Saudi Arabian women are associated with measures of glycemic control, which supports previous hypotheses [48,49]. For instance, some were positively correlated with HOMA-IR (Bacteroides_u_s, Bacteroides thetaiotaomicron) and insulin (Bacteroides thetaiotaomicron), while others were inversely correlated with HOMA-IR (Bifidobacterium adolescentis, Bifidobacterium kashiwanohense), fasting blood glucose (Actinobacteri, Bifidobacterium kashiwanohense), and insulin (Actinobacteria, Bifidobacterium adolescentis, Bifidobacterium kashiwanohense). Although this work investigated markers of diabetes rather than diabetes itself, the relationships found supporting studies showing that gut microflora differs between individuals with and without T2DM [50]. One study looked at this relationship specifically among patients in Saudi Arabia and found that Firmicutes levels were upregulated, and the Firmicutes/Bacteroidetes ratio was elevated among those with T2DM compared to those without [51]. Further, the study revealed that both men and women with T2DM had a higher abundance of prevotellaceae compared to those without [51], in contrast to earlier findings demonstrating that gut microbiota composition varies between genders [52,53]. In comparison to other populations, studies conducted on Chinese and Egyptian populations showed that the composition of gut microbiota is, in fact, associated with T2DM [54,55]. Some evidence showed that gut microbiota is influenced by geographic location and ethnicity [56,57]. For example, a study by He et al. revealed that the composition of the gut microbiota and its association with metabolic disease were strongly dependent on geography [56]. These data collectively suggest the possible role of gut microbiota in T2DM development.
The results can also be compared to prior work that looked at more specific types of bacteria. For instance, when combined with Methanobrevibacter smithii in germ-free mice, Bacteroides thetaiotaomicron increased short-chain fatty acid production [58], which has been demonstrated to lead to increased glycemic control [59]. Earlier work has also shown that the Actinobacteria phyla, Bifidobacteria, was inversely associated with insulin sensitivity [60]. Specifically, in agreement with the present study's results, Bifidobacterium

Discussion
This study aimed to investigate the relationship between insulin resistance, body fat distribution, and gut microbiota in young women in Saudi Arabia using the WGS technique. Results show that specific gut microbiota at different taxonomic levels in Saudi Arabian women are associated with measures of glycemic control, which supports previous hypotheses [48,49]. For instance, some were positively correlated with HOMA-IR (Bacteroides_u_s, Bacteroides thetaiotaomicron) and insulin (Bacteroides thetaiotaomicron), while others were inversely correlated with HOMA-IR (Bifidobacterium adolescentis, Bifidobacterium kashiwanohense), fasting blood glucose (Actinobacteri, Bifidobacterium kashiwanohense), and insulin (Actinobacteria, Bifidobacterium adolescentis, Bifidobacterium kashiwanohense). Although this work investigated markers of diabetes rather than diabetes itself, the relationships found supporting studies showing that gut microflora differs between individuals with and without T2DM [50]. One study looked at this relationship specifically among patients in Saudi Arabia and found that Firmicutes levels were upregulated, and the Firmicutes/Bacteroidetes ratio was elevated among those with T2DM compared to those without [51]. Further, the study revealed that both men and women with T2DM had a higher abundance of prevotellaceae compared to those without [51], in contrast to earlier findings demonstrating that gut microbiota composition varies between genders [52,53]. In comparison to other populations, studies conducted on Chinese and Egyptian populations showed that the composition of gut microbiota is, in fact, associated with T2DM [54,55]. Some evidence showed that gut microbiota is influenced by geographic location and ethnicity [56,57]. For example, a study by He et al. revealed that the composition of the gut microbiota and its association with metabolic disease were strongly dependent on geography [56]. These data collectively suggest the possible role of gut microbiota in T2DM development.
The results can also be compared to prior work that looked at more specific types of bacteria. For instance, when combined with Methanobrevibacter smithii in germ-free mice, Bacteroides thetaiotaomicron increased short-chain fatty acid production [58], which has been demonstrated to lead to increased glycemic control [59]. Earlier work has also shown that the Actinobacteria phyla, Bifidobacteria, was inversely associated with insulin sensitivity [60]. Specifically, in agreement with the present study's results, Bifidobacterium adolescentis has been associated with improved insulin sensitivity in a supplementation study on rats [61]. Further, in a more recent study, Bifidobacterium adolescentis was shown to help regulate blood glucose levels by increasing the amount of short-chain fatty acid-producing flora, alleviating inflammation, and restoring gut microbiota homeostasis in mice [62]. In a randomized, double-blind, cross-over trial investigating the effect of probiotic yogurt on metabolic and inflammatory biomarkers compared to acidified milk [63], Bifidobacterium kashiwanohense was less abundant after 2 weeks of the intervention with probiotic yogurt intake (p < 0.01). Interestingly, in our current work we found that Bifidobacterium kashiwanohense was inversely and significantly correlated with HOMA-IR, fasting blood glucose, and insulin. This finding highlights that Bifidobacterium kashiwanohense may be a potential biomarker of T2DM. Nevertheless, the role of Bifidobacterium kashiwanohense in T2DM development has not been fully identified in the literature, and future investigations are warranted.
When examined independently, there was no significant difference in αand βdiversity between the low and high HOMA-IR groups. However, α-diversity trended higher in the high HOMA-IR group. The previous study conducted in Saudi Arabia of individuals with and without T2DM also found a trend towards increased diversity among participants with T2DM and higher glucose levels [51], a finding similar to other work from the Middle East [64] but different from results from Western populations. For instance, a population-based study that used two European cohorts found that lower gut microbiota diversity was associated with higher levels of HOMA-IR irrespective of BMI [65]. Discrepancies in the findings between the present study and those previously conducted on different cohorts may be attributable to the differing microbial biomarkers between groups that could result from genetics, lifestyle, sex, age, geography-specific microbiota patterns, or measurement techniques [27,[32][33][34]66,67]. Thus, it is imperative to conduct locally representative studies to better understand the relationship between gut microbiota and health measures.
In the present work, adding different adiposity indices to analyses of HOMA-IR levels and diversity allowed for the emergence of significant findings. The nature of these findings appears to depend on the type of diversity and the type of adiposity index examined. For those with a high BF%, the α-diversity differed between those with low and high HOMA-IR. While for some groups, it appears the combination of an adverse HOMA-IR (high) and an adverse adiposity index led to significant differences (i.e., high WHR and obesity for α-diversity and high BF% for β-diversity). For other combinations, β-diversity differed within a HOMA-IR group primarily due to the stratification by a measure of adiposity (i.e., low vs. high WHR for high HOMA-IR and obesity vs. normal weight for low HOMA-IR). These findings suggest the importance of considering adiposity indices in these analyses. A previous study on this same cohort was the first to examine several obesity markers (i.e., WHR, BF%, and obesity) as they relate to gut microbiota among young women in the Middle Eastern region [37]. That study found αand β-diversity was significantly associated with WHR, and β-diversity was significantly associated with BF% [37]. The results from that work highlight the vital relationship between adiposity measures and microbiota diversity in this cohort and help explain the significance that emerged with the addition of adiposity indices in the present work.
Identifying the association between gut microbiota and insulin resistance can help develop strategies to promote insulin sensitivity and manage obesity. For example, gut microbiota signatures can be used to identify individuals at higher risk for diabetes and act as a target for microbial alterations to improve insulin sensitivity. In fact, recent work has demonstrated promising results for treating diabetes with fecal microbiota transplantation in mice [68]. Such efforts are particularly essential given the growing prevalence of T2DM in Saudi Arabia.
A strength of this work is that the study used the high-standard WGS method to determine gut microbiota at the species level, leading to higher sensitivity and resolution than other techniques. Repeated measurement tools were used to increase the accuracy of data collection methods. Furthermore, a combination of tools to measure microbial diversity (richness and variation) was used. On the other hand, this observational study does not allow for conclusions about causation, and there is a possibility of measurement error. Furthermore, the generalizability of these findings is additionally limited since only women participants in a single setting were included. However, selecting a specific cohort addresses the inability to generalize results from other cohorts to this population.
In conclusion, the present study was the first to examine markers of glycemic control, obesity, and gut microbiota among young women in Saudi Arabia. The results support others conducted in Saudi Arabia and Middle Eastern countries highlighting the benefit of using geography-specific cohorts. Lastly, the association between glycemic measures and microbial diversity was further enhanced by the inclusion of measures of adiposity. Our findings highlight the need for further clinical studies to investigate T2DM treatment through manipulation of gut microbiota approaches, such as the use of short-chain fatty acids, probiotics, prebiotics in the diet, fecal bacteria transplantation, or antibiotics. Such approaches will allow meaningful progress in the quest to understand the role of gut microbiota in the development of T2DM and will be an innovative step in the prevention and treatment of T2DM.