Study design

This study was carried out within three sub-cohorts of the Rotterdam Study (RS), a prospective cohort study of adult aged 45 years and older living in the well-defined district of Ommoord in Rotterdam, the Netherlands. A detailed description of the Rotterdam Study methodology is described elsewhere [12]. Briefly, recruitment of participants for the first sub-cohort (RS-I) started in the period of 1989–1993 among inhabitants aged ≥ 55 years (n = 7983). In 2000–2001, the study was extended with a second sub-cohort (RS-II) of new individuals (n = 3011) who had become 55 years of age or moved into the study area after 1990. In 2006–2008, a third sub-cohort (RS-III) was recruited with new individuals aged 45 years and older (n = 3932). By the end of 2008, the overall study population contained 14,926 participants. Upon entering the study, participants underwent home interviews and a series of examinations in our research center every 3–5 years.

The Rotterdam Study has been approved by the institutional review board (Medical Ethics Committee) of Erasmus Medical Center and by the review board of The Netherlands Ministry of Health, Welfare and Sports. The approval has been renewed every 5 years. All participants gave informed consent.

Population for current analyses

For the current study, we used data from all three sub-cohorts (Fig. 1). Of the 14,926 participants, we excluded those without valid dietary data (no dietary data (n = 5141) or unreliable dietary intake according to a trained nutritionist or an estimated energy intake of < 500 or > 5000 kcal/day (n = 84) [13]) at baseline (RS-I-1: 1989–1993, RS-II-1: 2000–2001, RS-III-1: 2006–2008), and those without diabetes information or with prevalent T2D at baseline (n = 2903), leaving 6798 participants included as main population for analysis.

Fig. 1 Flow diagram of participant selection Full size image

From this group of 6798 participants, 6514 participants had data on HOMA-IR before onset of T2D and were included in the longitudinal HOMA-IR analyses. For the analyses on prediabetes risk, we excluded those with prevalent prediabetes at baseline (n = 1005) or without follow-up of prediabetes (n = 25), leaving 5768 participants. In the analyses assessing risk of T2D, we excluded participants without follow-up of T2D (n = 28), leaving 6770 participants. The flow-diagram of the included participants is presented in Fig. 1.

Dietary assessment

Dietary intake was assessed at baseline in all three sub-cohorts using semi-quantitative food-frequency questionnaires (FFQ) as described in more detail elsewhere [13]. We used an FFQ with 170 food items to assess dietary intake at baseline of RS-I (1989–1993) and RS-II (2000–2001) [14]; and at baseline of RS-III (2006–2008) we used an FFQ with 389 food items [15]. The 170-item FFQ was validated in a subsample of the Rotterdam Study (n = 80) against fifteen 24-h food records and four 24 h urinary urea excretion samples [14]; and the 389-item FFQ was previously validated in other Dutch population against measurement of biomarkers, against a 9-day dietary record, and against a 4 week dietary history [16]. In general, the validation studies demonstrated that the FFQs were able to adequately rank participants according to their intake [13]. Food intake data were converted to energy and nutrient intake based on Dutch Food Composition tables (NEVO).

Plant-based dietary index

We constructed an overall plant-based dietary index, which was a modified version of two previously created indices [11, 17]. More specifically, our index is similar to the “provegetarian food pattern” of Martínez-Gonzáles et al. [17] and to the “overall plant-based diet index” of Satija et al. [11], but was adapted to include slightly different types and numbers of food categories.

First, the food items as measured by the FFQs were divided into 23 food categories (Supplemental Table 1), on the basis of the main food groups in the Dutch diet and the Dutch food-based dietary guidelines [18, 19]. Twelve of the categories were plant-based and eleven were animal-based. Food items that were not clearly animal-based or plant-based, such as pizza, as well dietary supplements, were not included in the food categories for the index.

Dietary intake for each of the 23 food categories (g/day) was calculated for each participant. Subsequently, for each category, the intake was divided into cohort-specific quintiles. Each quintile was assigned a value between 0 and 4. For the twelve plant-based food categories, consumption within the highest quintile was scored a 4, consumption within the second highest quintile was scored a 3, and so on, ending with consumption within the lowest quintile receiving a score of 0. The eleven animal-based food categories were scored reversely: consumption within the highest quintile was scored a 0 consumption within the second highest quintile was scored a 1, ending with consumption within the lowest quintile receiving a score of 4. Furthermore, we ensured that all participants with null-consumption were given the score belonging to the lowest quintile by re-scoring when necessary.

Finally, these category quintile-scores were added up for per participant to create their overall score on the plant-based dietary index. The resulting index yielded a score for each participant that measured adherence to a plant-based versus animal-based diet on a continuous scale, with a lowest possible score of 0 (low adherence to a plant-based diet) and a highest possible score of 92 (high adherence: high plant-based and low animal-based). Information on intake of each food category across quintiles of scores on the plant-based dietary index is shown in Supplemental Table 2.

Assessment of insulin resistance

Fasting blood samples were collected at RS-I (RS-I-3: 1997–1999, RS-I-5: 2009–2010), RS-II (RS-II-1: 2000–2001, RS-II-3: 2010–2011), and RS-III (RS-III-1: 2006–2008, RS-III-2: 2011–2012). Glucose levels were examined with the glucose hexokinase method. Serum insulin was measured by electro chemiluminescence immunoassay technology. Insulin resistance was calculated using the homeostasis model assessment of insulin resistance (HOMA-IR). The following formula was used: fasting insulin (mU/L) × fasting glucose (mmol/L)/22.5.

Assessment of prediabetes and type 2 diabetes

Information on prediabetes and T2D was collected from general practitioners’ records, pharmacies’ databases, and follow-up examinations in our research center. Data of prediabetes and T2D in our analyses were collected until January 1, 2012. Prediabetes and T2D were identified according to WHO criteria: prediabetes was defined as a fasting blood glucose concentration of > 6.0 and < 7.0 mmol/L, or a non-fasting blood glucose concentration of > 7.7 mmol/L and < 11.1 mmol/L; T2D was defined as a fasting blood glucose concentration of ≥ 7.0 mmol/L, a non-fasting blood glucose concentration of ≥ 11.1 mmol/L (when fasting samples were unavailable), or the use of blood glucose-lowering drugs or dietary treatment and registration of the diagnosis diabetes. All possible cases of prediabetes and T2D were formally judged by two independently working study physicians or, in case of disagreement, by an endocrinologist [20].

Assessment of covariates

Information on age, sex, smoking status, educational level, medication use, food supplement use, and family history of diabetes, was obtained from questionnaires at baseline. Information on physical activity was obtained using the adapted version of the Zutphen Physical Activity Questionnaire at RS-I-3 and RS-II-1, and using the LASA Physical Activity Questionnaire at RS-III-1. Physical activities were weighted according to intensity with Metabolic Equivalent of Task (MET), from the Compendium of Physical Activities version 2011. To account for differences between the two questionnaires, questionnaire-specific z-scores of MET-hours per week were calculated. At our research center at baseline, body weight was measured using a digital scale and body height was measured using a stadiometer, while participants wore light clothing and no shoes, and BMI was calculated (kg/m2). Information on hypertension, hypercholesterolemia, coronary heart disease (CHD), cancers, and stroke was obtained from general practitioners, pharmacies’ databases, Nationwide Medical Register, or follow-up examinations in our research center.

Data analyses

To obtain a normal distribution for HOMA-IR, we applied a natural-log transformation. Non-linearity of associations of score on the plant-based dietary index with all outcomes were explored using natural cubic splines (degrees of freedom = 3). As no indications for non-linear associations for the main models were found, all primary analyses were performed using models assuming linearity. We examined the association between score on the plant-based dietary index with longitudinal HOMA-IR using linear mixed models, with a random-effects structure including a random intercept and slope (for time of repeated measurements of HOMA-IR). We examined the association between score on the plant-based dietary index and risk of prediabetes and risk of T2D using Cox proportional-hazards regressions. Hazard ratios (HRs) and regression coefficients (βs) were presented per 10 units higher score on the plant-based dietary index, along with the corresponding 95% confidence intervals (CIs). All analyses were performed in participants of the three sub-cohorts combined and in the three sub-cohorts separately.

All analyses were adjusted for energy intake, age, sex and RS sub-cohort in model 1, and for the analyses of longitudinal HOMA-IR we additionally adjusted for the time of repeated measurements of HOMA-IR. In model 2, we additionally adjusted for smoking status, educational level, physical activity, food supplement use, and family history of diabetes. Baseline BMI was added to model 3 to examine its potential mediating effect.

We examined effect modification by including interactions of the plant-based index with age, sex, or BMI for all outcomes in model 2.

Several sensitivity analyses were performed based on model 2. First, to check if the associations were driven by any specific components of the plant-based dietary index, we repeated our main analyses by excluding each one of the 23 components from the plant-based dietary index one by one at a time, and additionally adjusting for the excluded component. Second, to check if the associations were mainly driven by plant-based beverages combined, we examined the associations by excluding all plant-based beverages combined (category “coffee and tea”, category “alcoholic beverages”, and category “sugary beverages”) from the plant-based dietary index at a time, and additionally adjusting for them. Third, we examined the associations by excluding less healthy plant-based foods combined (category “sweets”, category “sugary beverages”, category “potatoes”, and category “refined grains”) from the plant-based dietary index at a time, and additionally adjusting for them. To further examine whether these less healthy plant foods contributed to the association of the plant-based dietary index; we created a less healthy plant foods score, for which, positive scores were given to these four types of less healthy plant-based food groups; and reverse scores were given to healthy plant food groups and animal food groups [21]. Fourth, to examine if potential associations of the plant-based dietary score with outcomes were independent of overall quality of the diet based on adherence to dietary guidelines, we examined the correlation between the plant-based dietary score and the dietary guidelines score; and we repeated analyses with additional adjustment for dietary guidelines score. Fifth, we additionally adjusted for hypertension and hypercholesterolemia. Sixth, we excluded the participants with chronic diseases at baseline, such as participants with CHD, cancers, or stroke, to exclude the possibility of a significant change of diet and life style at follow-up. Last, we excluded the participants who developed prediabetes and T2D in the first 2 years of follow-up in the analyses for risk of prediabetes and T2D, respectively.

Missing values on covariates (ranging from 0.3 to 3.9%) were accounted for using multiple imputations (n = 10 imputations). We used SPSS version 21 (IBM Corp., Armonk, NY, USA) and R version 3.1.2 (R Foundation for Statistical Computing, Vienna, Austria) to perform these analyses.