nutrients Article Gut Microbiota Bacterial Species Associated with Mediterranean Diet-Related Food Groups in a Northern Spanish Population Carles Rosés 1 , Amanda Cuevas-Sierra 2 , Salvador Quintana 3 , José I. Riezu-Boj 2,4 , J. Alfredo Martínez 2,4,5 , Fermín I. Milagro 2,4,5 and Anna Barceló 1, * 1 Servei de Genòmica i Bioinformàtica, Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain; carles.roses@uab.cat 2 Center for Nutrition Research, Department of Nutrition, Food Sciences and Physiology, University of Navarra, 31008 Pamplona, Spain; acuevas.1@alumni.unav.es (A.C.-S.); jiriezu@unav.es (J.I.R.-B.); jalfmtz@unav.es (J.A.M.); fmilagro@unav.es (F.I.M.) 3 Independent Researcher, 08021 Barcelona, Spain; squintana@mutuaterrassa.es 4 Navarra Institute for Health Research (IdISNA), 31008 Pamplona, Spain 5 Centro de Investigación Biomédica en Red de la Fisiopatología de la Obesidad y Nutrición (CIBERobn), Instituto de Salud Carlos III, 28029 Madrid, Spain * Correspondence: anna.barcelo@uab.cat Abstract: The MD (Mediterranean diet) is recognized as one of the healthiest diets worldwide and is associated with the prevention of cardiovascular and metabolic diseases. Dietary habits are considered one of the strongest modulators of gut microbiota, which seem to play a significant role   in health status of the host. The purpose of the present study was to evaluate interactive associations between gut microbiota composition and habitual dietary intake in 360 Spanish adults from the Citation: Rosés, C.; Cuevas-Sierra, Obekit cohort (normal weight, overweight, and obese participants). Dietary intake and adherence to A.; Quintana, S.; Riezu-Boj, J.I.; the MD tests were administered and fecal samples were collected from each participant. Fecal 16S Martínez, J.A.; Milagro, F.I.; Barceló, A. Gut Microbiota Bacterial Species rRNA (ribosomal Ribonucleic Acid) gene sequencing was performed and checked against the dietary Associated with Mediterranean habits. MetagenomeSeq was the statistical tool applied to analyze data at the species taxonomic Diet-Related Food Groups in a level. Results from this study identified several beneficial bacteria that were more abundant in Northern Spanish Population. the individuals with higher adherence to the MD. Bifidobacterium animalis was the species with the Nutrients 2021, 13, 636. https:// strongest association with the MD. Some SCFA (Short Chain Fatty Acids) -producing bacteria were doi.org/10.3390/nu13020636 also associated with MD. In conclusion, this study showed that MD, fiber, legumes, vegetable, fruit, and nut intake are associated with an increase in butyrate-producing taxa such as Roseburia faecis, Academic Editor: Jose Lara Ruminococcus bromii, and Oscillospira (Flavonifractor) plautii. Received: 17 January 2021 Accepted: 13 February 2021 Keywords: Bifidobacterium animalis; gut microbiota; short-chain fatty acids; obesity; butyrate Published: 16 February 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in 1. Introduction published maps and institutional affil- iations. Gut microbiota status has an impact on the health and disease of the host [1]. Dietary habits are considered one of the strongest modulators of gut microbiota. However, it is not clear how the two relate to each other or which gut microbiota profile to pursue in order to ensure good health. Currently, the majority of the world’s population is exposed to Western diets, characterized by a high intake of saturated and omega-6 fatty acids, Copyright: © 2021 by the authors. Licensee MDPI, Basel, Switzerland. reduced omega-3 and fiber intake, an overuse of salt, and too much refined sugar and This article is an open access article processed food [2]. Along with a sedentary lifestyle, these factors are increasing the preva- distributed under the terms and lence of obesity worldwide, with half of the world’s population now considered to be conditions of the Creative Commons overweight [3]. As a result of this lifestyle, serious conditions can appear and in a variety Attribution (CC BY) license (https:// of ways. Hypertrophied adipocytes release inflammatory molecules (i.e., interleukins and creativecommons.org/licenses/by/ tumor necrosis factor), which can act as false alarms in the immune system, causing the en- 4.0/). tire immune system to reduce its sensitivity, such that the response to a real condition may Nutrients 2021, 13, 636. https://doi.org/10.3390/nu13020636 https://www.mdpi.com/journal/nutrients Nutrients 2021, 13, 636 2 of 17 be delayed [4]. This can aid the development of several inflammation-related disorders such as metabolic syndrome, cardiovascular disease, colorectal cancer, and neurodegenera- tive diseases [5,6]. Most of these disorders have also been associated with alterations in microbiota composition in humans, especially those with reduced bacterial richness and diversity [7]. These changes have been related to disturbed gut barrier functions, increased gut permeability, and increased plasma concentrations of lipopolysaccharide (LPS) and other bacterial by-products, which cause low-grade inflammation that, again, triggers the development of insulin resistance, obesity, metabolic syndrome, colorectal cancer [8], and autoimmune disorders such as Crohn disease, ulcerative colitis, and allergies [7]. In this context, the MD (Mediterranean diet) is recognized as one of the healthiest diets worldwide as it contains to a high proportion of fiber, mono- and poly-unsaturated fatty acids, antioxidants, and polyphenols, present in vegetables, fruits, pulses, and extra virgin olive oil (EVOO), which are all strongly associated with a reduced risk of devel- oping non-communicable diseases related to Western diet and lifestyle [9,10]. In a study with overweight, obese participants with lifestyle risk factors for metabolic disease, an isocaloric MD intervention reduced their blood cholesterol and caused multiple changes in their microbiome and metabolome, thus improving their metabolomic health [11]. A 12-month-long MD intervention with elderly subjects showed a taxa enrichment that was associated with lower frailty and improved cognitive functions, but negatively correlated with inflammatory markers, thus promoting healthier aging [12]. Carbohydrates and fiber present in the MD are fermented by gut microbiota, through which large quantities of biologically active metabolites such as SCFAs (short-chain fatty acids) are produced [13]. Also, greater concentrations of phenolic metabolites are excreted in feces when there are high concentrations of bioactive compounds coming from polyphenols and fiber [14]. Fur- thermore, de Filippis et al. reported that a high adherence to a MD rich in plant foods beneficially impacts gut microbiota and the associated metabolome [15]. Therefore, we would expect a modulation of the gut microbiota to be one of the positive health effects of the MD [13]. Mitsou et al. found a positive correlation with gastrointestinal symptoms, fe- cal moisture, total bacteria, and Bifidobacteria, but a reduced representation of Lactobacilli and butyrate-producing bacteria induced by fast food consumption [16]. The main objective of the present study was to relate the adherence to the MD to specific metagenomic traits, focusing on those bacterial taxa that are more abundant in individuals with a high adherence to the MD. In addition, a specific food group consump- tion assessment was carried out to increase our knowledge of the impact of diet on gut microbiota composition. The focus was on the bacteria that are most closely associated with a high adherence to the MD and how these bacteria are influenced by specific food groups characteristic of the Mediterranean pattern. 2. Material and Methods 2.1. Participants This cross-sectional study enrolled 360 Spanish adults of self-reported European ances- try (251 females and 109 males) with ages ranging from 45.0 ± 10.5 years old. Participants were recruited at the Center for Nutrition Research of the University of Navarra, Spain, and took part in the Obekit study. Major exclusion criteria included a history of diabetes mellitus, cardiovascular disease and hypertension, pregnant or lactating women, and current use of lipid-lowering drugs. Patients with a diagnosis of primary hyperlipidemia were also excluded. This investigation followed the ethical principles for medical research in humans from the 2013 Helsinki Declaration [17]. The research protocol (ref. 132/2015) was approved by the Research Ethics Committee of the University of Navarra. Written informed consent from each participant was obtained before the inclusion in the study. The characteristics of the Obekit research project, including study design and registration, have been reported elsewhere [18]. Inclusion criteria were body mass index between 25 and 40, a physical examination and assessment of vital signs considered as normal or clinically insignificant by the researcher. Nutrients 2021, 13, 636 3 of 17 Also, in the case of individuals with chronic, stable-dose drug treatment during the last three previous months and at baseline, the investigators assessed possible inclusion. 2.2. Anthropometric and Biochemical Measurements Anthropometric measurements such as height (cm) and body weight (kg) were col- lected in the fasting state by trained nutritionists following validated procedures [19]. Body mass index (BMI) was calculated as the ratio between weight and squared height (kg/m2 ), and, according to the WHO (World Health Organization) standards, the volunteers were classified as normal weight when BMI was 18.5–24.9 (n = 64), overweight when BMI was 25.0–29.9 (n = 115), and obese when BMI > 30.0 (n = 181). Venous blood samples were drawn by venipuncture in a clinical setting after an overnight fast. Blood biochemistry (glucose, total cholesterol (TC), high-density lipoprotein cholesterol (HDL), and triglycerides) was analyzed with a Pentra C200 clinical chemistry analyzer (HORIBA Medical, Madrid, Spain) and suitable kits provided by the company. Low-density lipoprotein cholesterol (LDL-c) was estimated using the Friedewald equation: LDL-c = TC-HDL-(triglycerides/5). Insulin was measured using a specific enzyme-linked immunosorbent assay (Mercodia, Uppsala, Sweden) read with an automated analyzer system (Triturus, Grifols, Barcelona, Spain). The homeostatic model assessment for insulin resistance (HOMA-IR) was calculated as follows: fasting insulin (microU/L) × fasting glucose (nmol/L)/22.5. 2.3. Dietary Estimation Habitual dietary intake at baseline was collected with a validated food frequency questionnaire that included 137 food items with corresponding portion sizes [20]. All participants enrolled were asked to provide information about the number of times they had consumed each food item during the previous year according to four frequency categories: daily, weekly, monthly, or never. Total energy (kcal) and macronutrient intakes (%) were determined with ad hoc software and the information available from valid Spanish food composition tables [21]. Among all this information, six groups were defined as the most representative of the MD: fiber, legumes, vegetables, fruit, olive oil, and nuts. Consumption of these groups was represented in g per day and these scores were transformed into tertiles as cut-off points. The second tertile was removed and first and third were compared. A 14-item questionnaire, PREDIMED (Prevención con Dieta Mediterránea) validated test, was also used in this study to appraise the adherence of participants to the MD [22]. The MD score ranged from 0 (minimal adherence) to 14 (maximal adherence). Ranges of scores were defined again by tertiles: first tertile (0–6 score) for low, second tertile (7–8) was removed, and third tertile (9–14) for high adhesion to the MD. 2.4. Faecal Sample Collection and DNA Extraction Volunteers self-collected fecal samples at baseline using OMNIgene.GUT kits from DNA Genotek (Ottawa, ON, Canada), according to the standard instructions provided by the company. The isolation of DNA from fecal samples was performed with QIAamp® DNA kit (Qiagen, Hilden, Germany), following the manufacturer’s protocol. 2.5. Metagenomic Data: Library Preparation Bacterial DNA sequencing was performed by the Servei de Genòmica i Bioinformàtica from the Universitat Autònoma de Barcelona (Bellaterra, Cerdanyola del Vallés, Spain). Metagenomics studies were performed by analyzing the variable regions V3–V4 of the prokaryotic 16S rRNA (ribosomal Ribonucleic Acid) gene sequences, which gives 460 bp amplicons in a two-round PCR protocol. In a first step, PCR is used to amplify a template out of a DNA sample using specific primers with overhang adapters attached that flank regions of interest. The full-length primer sequences, using standard IUPAC (International Union of Pure and Applied Chem- istry) nucleotide codes, to follow the protocol targeting this region were: Forward Primer: Nutrients 2021, 13, 636 4 of 17 50 TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCCTACGGGNGGCWGCAG and Reverse Primer: 50 GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGACTACHV GGGTATCTAATCC. PCR was performed in a thermal cycler using the following condi- tions: 95 ◦ C for 3 min, 25 cycles of (95 ◦ C for 30 s, 55 ◦ C for 30 s, and 72 ◦ C for 30 s), and 72 ◦ C for 5 min. To verify that the specific primers had been correctly attached to the samples, 1 µL of the PCR product was checked on a Bioanalyzer DNA 1000 chip (Agilent Technologies, Santa Clara, CA, USA). The expected size on a Bioanalyzer is ~550 bp. In a second step and using a limited-cycle PCR, sequencing adapters, and dual indices barcodes, Nextera® XT DNA Index Kit, FC-131-1002 (Illumina, San Diego, CA, USA), were added to the amplicon, which allows up to 96 libraries for sequencing on the MiSeq sequencer with the MiSeq® Reagent Kit v3 (600 cycle) MS-102-3003 to be pooled together. PCR was performed in a thermal cycler using the following conditions: 95 ◦ C for 3 min, eight cycles of (95 ◦ C for 30 s, 55 ◦ C for 30 s, and 72 ◦ C for 30 s), and 72 ◦ C for 5 min. Subsequently, the Index PCR ran a second Bioanalyzer DNA 1000 chip to validate the library. The expected size was ~630 bp. The next step consisted of the quantification of the libraries using a fluorometric quantification and dilution of the samples before pooling all samples. Finally, paired-end sequencing was performed on a MiSeq platform (Illumina) with a 600 cycles Miseq run [23] and with 20 pM sample and a minimum of 20% PhiX. The mean reads obtained were 164,387. Only samples with more than 40,000 reads were used for further analysis. All the sequencing data were deposited by the authors in SRA (Sequence Read Archive) and the accession key has been included in the text (PRJNA623853). 2.6. Metagenomics Data: Analysis and Processing The 16S rRNA gene sequences obtained were filtered following the quality criteria of the OTUs (operational taxonomic units) processing pipeline LotuS (release 1.58) [24]. This pipeline includes UPARSE (Highly accurate OTU sequences from microbial amplicon reads) de novo sequence clustering and removal of chimeric sequences and phix contaminants for the identification of OTUs and their abundance matrix generation [25,26]. Taxonomy was assigned using HITdb (Highly scalable Relational Database), achieving up to species sensitivity. BLAST (Basic Local Alignment Search Tool) was used when HITdb failed to reach a homology higher than 97% [27,28]. Thus, OTUs with a similarity of 97% or more were referred to species. However, OTUs that did not reach this percentage of similarity were checked and updated using the Basic Local Alignment Search Tool (BLASTn) to compare with the 16S rRNA gene sequences for bacteria and archaea database of GenBank of National Center for Biotechnology information in order to find an assignment to a species. These sequences in which the BLASTn tool found a new assignment were indicated using the GenBank access number and the percentage of homology following the species name. The abundance matrices were first filtered and then normalized in R/Bioconductor at each classification level: OTU, species, genus, family, order, class, and phylum. This study focused mainly on the species level. Briefly, taxa with a less than 10% frequency in our population were removed from the analysis, and a global normalization was performed using the library size as a correcting factor and log2 data transformation [29]. 2.7. Richness and Evenness Richness was defined as the total of species. Evenness was calculated using Pielou’s evenness index according to the following formula: J’ = H/ln(S). Alpha diversity was assessed using the Shannon index. Beta diversity was calcu- lated using the Bray–Curtis index, PERMANOVA (Permutational Multivariate Analysis of Variance) statistical method, and NMDS (Non-metric Multidimensional Scaling) as ordination methods. Nutrients 2021, 13, 636 5 of 17 All the sequencing data were deposited by the authors in the Sequence Read Archive (SRA) from NCBI (National Center for Biotechnology Information), with PRJNA623853 as accession key. 2.8. Statistical Analysis The microbiome Analyst tool [30] was used for statistical differences in microbiota pro- files between groups (tertiles) through a Zero-inflated Gaussian approach of Metagenome- Seq and using the cumulative sum scaling (CSS) normalization. 2.9. Prediction of Functional Potential of Gut Microbiota Computational prediction of the functional capabilities using data from 16S rRNA metagenomics was performed using the Tax4fun tool from MicrobiomeAnalyst [30]. KEGG (Kyoto Encyclopedia of Genes and Genomes) Orthologus (KO) provided by Tax4Fun was comparatively analyzed using the Shotgun Data Profiling section in the MicrobiomeAnalyst. A total of 2094 KO low-abundance features were removed based on prevalence (<20% of prevalence in samples) and 336 low-variance features based on inter-quantile range, with a total of 3018 features remaining after filtering. The EdgeR statistical analysis compared differentially expressed KO between groups (low and high adherence to MD) and allowed to visualization of the results within KEGG metabolic networks, along with pathway analysis between groups. 3. Results 3.1. Participant Characteristics Baseline characteristics of the population that participated in this study separated by adherence to the MD and BMI (normal weight, overweight, and obese, according to the World Health Organization criteria [31]) are shown in Table 1, including age, an- thropometric measures, and biochemical and dietary data. Additional data on the study population (separated by sex and adiposity status) are shown as supplementary material (Supplementary Table S1). Table 1. Baseline characteristics of the population separated by adherence to the MD (Mediter- ranean diet). High Adherence Low Adherence Variables (3rd Tertile) (1st Tertile) p Value (n = 94) (n = 128) Age (y) 47.3 ± 1.2 41.7 ± 0.9 0.020 * BMI 27.7 ± 0.6 29.9 ± 0.6 0.002 * Glucose (mg/dL) 95 ± 2 93 ± 1 0.394 Total cholesterol (mg/dL) 209 ± 4 208 ± 4 0.822 HDL (mg/dL) 59 ± 1 57 ± 1 0.261 LDL (mg/dL) 133 ± 4 132 ± 3 0.904 Triglycerides (mg/dL) 87 ± 5 93 ± 4 0.336 HOMA-IR 1.4 ± 0.15 1.6 ± 0.1 0.011 * Carbohydrate intake (%) 54.5 ± 0.7 54.7 ± 0.6 0.135 Protein intake (%) 21.8 ± 0.4 22.2 ± 0.4 0.286 Fat intake (%) 23.6 ± 0.5 23.0 ± 0.4 0.085 Fiber intake (g/day) 35.9 ± 1.6 24.5 ± 0.7 0.000 * Energy intake (kcal/day) 2932 ± 106 2872 ± 84 0.680 The p value was calculated depending on the distribution of the variables by T-Student or Mann–Whitney test. * The score is significant at the 0.05 level. The high-adherence group was made up of older people, with a significantly lower BMI and less resistance to insulin. Despite the total energy intake being very close, the data in Table 1 indicate that the third tertile values were healthier. Concerning dietary composition, the third tertile was characterized by a higher intake of fiber. Nutrients 2021, 13, 636 6 of 17 3.2. Microbiota Composition: MD Adherence MD tertiles 1 and 3 were compared through metagenomeSeq analysis. Significant differences appeared when comparing both tertiles (FDR (False Discovery Rate) < 0.05). Species shown in Table 2 were strongly influenced by the MD score. Participants with a higher adherence to the MD were represented in the third tertile while those who are far from the MD model were in the first tertile. Table 2. Bacterial species with a significant relation with adherence to the MD (FDR (False discovery rate) < 0.05) by metagenomeSeq test. High Adherence (3rd Tertile) Low Adherence (1st Tertile) SPECIES FDR SPECIES FDR OTU100|NN = Eubacterium saphenum Bifidobacterium animalis 1.21 × 10−7 4.44 × 10−5 GU427005|D = 91 OTU375|NN = Succinivibrio Bacteroides cellulosilyticus 4.47 × 10−7 0.0001 dextrinosolvens Y17600|D = 97 OTU946|NN = Paraprevotella clara OTU759|NN = Gordonibacter pamelaeae 1.72 × 10−5 0.0005 AB331896|D = 86.8 AB566419|D = 87.6 OTU1682|NN = Oscillibacter valericigenes OTU11|NN = Butyricicoccus pullicaecorum 3.42 × 10−5 0.0002 AB238598|D = 91.1 EU410376|D = 89 OTU1065|NN = Oscillospira (Flavonifractor) 3.42 × 10−5 Christensenella minuta 0.0020 plautii Y18187|D = 86.6 OTU1173|NN = Roseburia faecis 0.0008 Parabacteroides goldsteinii 0.0073 AY804149|D = 94.9 OTU1517|NN = Catabacter hongkongensis OTU1625|NN = Anaerotruncus colihominis 0.0008 0.0120 AB671763|D = 87 DQ002932|D = 89 OTU1296|NN = Ruminococcus bromii 0.0120 Alistipes timonensis 0.0155 DQ882649|D = 92.3 Erysipelatoclostridium ramosum 0.0176 Prevotella corporis 0.0192 OTU521|NN = Papillibacter cinnamivorans 0.0463 AF167711|D = 89 This work focused on the high-adherence species and their distribution, with box plots (Figure 1). All box plots represent those species with significant differences between high- and low-adherence tertiles to the respective dietary pattern or food group. According to the species shown in Table 2, the relationship between these bacteria and the intake of certain foods (g of food per day) was analyzed. The food groups chosen for this approach were those most relevant in the MD: legumes, vegetables, fruit, olive oil, nuts, and also total fiber. 3.3. Microbiota Composition: Fibre Intake It has been proposed that a substantial part of the beneficial effects of the MD could be attributed to a high intake of fiber. Species upregulated in the high MD adherence group were analyzed to see which of them were related to fiber intake (FDR < 0.05) (tertiles 1 and 3). The most abundant species in the tertile of higher fiber intake (third tertile) are shown in Table 3 and Figure 2. Nutrients 2021, 13, 636 Nutrients 2021, 13, x FOR PEER REVIEW 8 of 19 7 of 17 OTU1682|NN=Oscillibacter valerici- OTU1065|NN=Oscillospira (Flavonifractor) OTU1173|NN=Roseburia fae- Bifidobacterium animalis genes plautii cis OTU1296|NN=Ruminococcus Erysipelatoclostridium ramosum OTU521|NN=Papillibacter cinnamivorans Bacteroides cellulosilyticus bromii OTU946|NN=Paraprevotella clara OTU1517|NN=Catabacter hongkongensis Figure 1. Figure 1. species Bacterial Bacterial species that were that were significantly significantly more in more abundant abundant in with the group the group with high adherence high adherence to MD (FDRto