We used data from The Maastricht Study, a prospectively designed, population-based observational cohort study. The rationale and methodology have been described previously. [9] In brief, the study focuses on the etiology, pathophysiology, complications and comorbidities of type 2 diabetes mellitus and is characterized by an extensive phenotyping approach. Eligible for participation were all individuals aged between 40 and 75 years and living in the southern part of the Netherlands. Participants were recruited through mass media campaigns and from the municipal registries and the regional Diabetes Patient Registry via mailings. Recruitment was stratified according to known type 2 diabetes status, with an oversampling of individuals with type 2 diabetes, for reasons of efficiency. [9] The present report includes cross-sectional data of 3,451 participants who completed the baseline survey between November 2010 and September 2013. The examinations of each participant were performed within a time window of three months. The study has been approved by the institutional medical ethics committee (NL31329.068.10) and the Minister of Health, Welfare and Sports of the Netherlands (Permit 131088-105234-PG). All participants gave written informed consent. [9].

Heart rate variability was assessed from ECGs, as reported previously. [10] All ECGs were recorded by use of a 12-lead Holter system (Fysiologic ECG Services, Amsterdam, the Netherlands) over a 24-h period. During recording time, participants were asked to follow their normal daily activities, except that they were asked not to take a shower or a bath. Recordings were analyzed with proprietary Holter Analysis Software at Fysiologic ECG Services (Fysiologic ECG Services, Amsterdam, the Netherlands) with an algorithm that excluded non-sinus cardiac cycles (e.g., artifacts and premature/ectopic beats), validated by manual inspection afterward. The software from Fysiologic ECG Services provided the intervals between the individual R waves of sinus beats (i.e. interbeat intervals, in milliseconds). From the obtained interbeat intervals, HRV was calculated by use of the publicly available free GNU Octave software, [11] according to the standard time- and frequency- domain measures defined by the (recently updated) recommendations of the Task Force document on HRV. [6, 12] After exclusion of non-sinus cardiac cycles, the minimum duration of the recording for HRV analysis was 18 h. We calculated the following time domain measures: the standard deviation (SD) of all normal-to-normal (NN) intervals (SDNN, in milliseconds [ms]); the SD of the averages of NN intervals in all 5-min segments of the entire recording (SDANN, in ms); the square root of the mean of the sum of squares of differences between adjacent NN intervals (RMSSD, in ms); the mean of the SDs of all NN intervals for all 5-min segments of the entire recording (SDNN index, in ms); the number of pairs of adjacent NN intervals differing by 50 ms in the entire recording (NN50 count, number); and the NN50 count divided by the total number of all NN intervals (pNN50, a percentage). Then, we calculated the following frequency domain measures using the Fast Fourier Transform: the variance of all NN intervals ≤ 0.4 Hz (total power [TP], in ms squared); power in the ultralow-frequency range (ULF, in ms squared; ≤0.003 Hz); power in the very-low-frequency range (VLF, in ms squared; 0.003–0.04 Hz); power in the low-frequency range (LF, in ms squared; 0.04–0.15 Hz); and the power in the high-frequency range (HF, in ms squared; 0.15–0.4 Hz). Individual z-scores were calculated for the time- and frequency- domain measures. We calculated the overall time- domain variable composite z-score using the following formula: ([SDNN + RMSSD + SDANN + SDNN index + pNN50]/5); and calculated overall frequency- domain variable using the following formula: ([TP + ULF + VLF + LF + HF]/5). Before we composed the composite z-score we checked whether associations of individual time- and frequency-domain HRV indices with individual measures of beta cell response were directionally consistent (in the complete study population), and this was the case.

After an overnight fast, all participants (except those who used insulin or had a fasting plasma glucose concentration above 11.0 mmol/L) underwent a standard 2-hour 75 g OGTT. Venous blood samples were collected before and at 15, 30, 45, 60, 90, and 120 min post oral glucose load intake. We estimated beta cell response using the following parameters: C-peptidogenic index t0-30, overall insulin secretion (C-peptide AUC/ glucose AUC [CPAUC/GAUC]), beta cell glucose sensitivity, beta cell potentiation factor, and beta cell rate sensitivity. We used formula-based methods and mathematical modeling, as previously described. [5, 13] First, C-peptidogenic index t0-30, which reflects early insulin secretion after a glucose stimulus, was calculated as the difference in C-peptide level between baseline and 30 min post glucose load intake (ΔCP30) divided by the difference in glucose level between baseline and 30 min post glucose load intake (ΔG30) [5, 13]. Second, overall insulin secretion, an analogous index reflecting the overall OGTT response, was calculated as the area under the curve (AUC) of C-peptide divided by the AUC of glucose (CPAUC/GAUC) [5, 13]. To calculate AUC, C-peptide or glucose were plotted against time. Third, beta-cell function parameters were calculated using a previously developed model, [5] which describes the relationship between insulin secretion and glucose concentration by means of a dose-response function relating the two variables and an early secretion component. The dose-response is characterized by its average slope, termed glucose sensitivity, and early secretion by a parameter denoted as rate sensitivity, a marker of early phase insulin release. The dose-response function is modulated by a time-varying potentiation factor, which accounts for effects of sustained hyperglycemia and incretins. The potentiation factor excursion was calculated as the ratio between the values at the end of the 2-h OGTT and at baseline.

Next, we calculated an overall beta cell response composite score (“overall beta cell response”). We standardized individual beta cell response measures (i.e. expressed as z-scores) and calculated the overall beta cell response composite score as follows: (C-peptidogenic index t0-30 + overall insulin secretion + beta cell glucose sensitivity + beta cell potentiation factor + beta cell rate sensitivity)/5. Then, we re-standardized the composite score. Before we composed the composite z-score we checked whether associations of time- and frequency-domain HRV with individual measures of beta cell response were directionally consistent (in the complete study population), and this was the case.

As described previously, [9] we assessed educational level (low, intermediate, high), socio-economic status (income level and occupational status), smoking status (never, former, current), alcohol use (none, low, high), and history of cardiovascular disease by questionnaire; assessed total energy intake with a food frequency questionnaire; [14] assessed lipid-modifying, antihypertensive, and glucose-lowering medication use as part of a medication interview; assessed weight, height, and waist circumference during a physical examination; calculated body-mass index (BMI) based on body weight and height; measured office and 24-hour ambulatory blood pressure; measured total daily physical activity (hours/day) with an accelerometer; measured lipid profile in fasting venous plasma samples; estimated the Matsuda Index, an index of insulin sensitivity, using data from the OGTT; and calculated the estimated glomerular filtration rate (eGFR) based on serum creatinine and cystatin C.

Characteristics of the study population were described as means and standard deviations (SD) for continuous variables or as number and proportions of participants per category for categorical variables (% of study population). We used multivariable linear regression analyses to investigate the associations of the independent variables (time- and frequency- domain HRV) with the dependent variables (overall beta cell response composite score; and C-peptidogenic index, overall insulin secretion, beta cell glucose sensitivity, beta cell potentiation factor, and beta cell rate sensitivity separately). Before we performed these analyses, we checked whether assumptions of linear regression were not violated and this was the case (we checked whether associations were linear, whether the residuals were normally distributed, and whether the homoscedasticity [homogenous distribution of variance] assumption was not violated). We checked these assumptions by visually inspecting relevant plots. These associations were presented as standardized betas with corresponding 95% confidence intervals [CI]). In addition, we inversed HRV so that associations were expressed per SD lower HRV (indicating progressively worse autonomic function). Last, we used complete case analysis.

We first analyzed associations without adjustment (crude model). In model 1, we adjusted for age, sex, and educational status (low, medium, high). We chose these variables because they are key demographic potential confounders. In model 2, we additionally adjusted for Matsuda Index, an index of insulin resistance. [7] We entered Matsuda Index in a separate model because beta-cell function is influenced by insulin resistance. In model 3, we additionally adjusted for cardiovascular risk factors and lifestyle factors. In model 3 A, we additionally adjusted for: BMI, total cholesterol / HDL cholesterol ratio, lipid-modifying medication (yes/no), smoking status (current, former, never), and alcohol consumption status (none, low, high). We adjusted for these potential confounders in a separate model as on one hand these factors may be confounders, but on the other hand these factors are also on the causal pathway, thus adjustment for these factors may be overadjustment (i.e. these factors can induce autonomic dysfunction which can lead to beta cell dysfunction). [15] In model 3B, we additionally adjusted for office systolic blood pressure and use of antihypertensive medication (yes/no). We adjusted for blood pressure in a separate model because blood pressure may be a confounder; a cause of autonomic dysfunction; and a descending proxy of autonomic function. [15, 16]

We tested for interaction by sex and glucose metabolism status to investigate whether associations under study differed by sex (i.e. between men and women) or glucose metabolism status (i.e. between individuals with type 2 diabetes, prediabetes, or normal glucose metabolism). We tested for interaction in the fully adjusted model by including interaction terms with the determinant (e.g. frequency-domain HRV*sex) and all covariates (e.g. age*sex), as previously described. [17].

To assess the robustness of our findings we performed a number of additional analyses. First, we analyzed the associations with the overall beta cell response composite score and individual beta cell response indices of individual measures of HRV (i.e. individual measures that were used to compose time and frequency domain composite scores). Second, and only for associations of time-and frequency-domain HRV with beta cell rate sensitivity as outcome, we performed logistic regression analyses in which beta cell rate sensitivity was categorized into tertiles (we performed this analysis for statistical reasons because the distribution of beta cell rate sensitivity was somewhat skewed and could not be normalized via logarithmical transformations). Third, we repeated the analyses with additional adjustment for total energy intake and physical activity. These potential confounders were not included in the main analyses because data were missing for a relatively large number of participants (up to n = 311 had missing data on one or more of these variables). Fourth, we additionally adjusted for kidney function (eGFR) and history of cardiovascular disease. We adjusted for these covariates in a separate model because they may be confounders but may also (in part) be mediators or descendants of the outcome. Fifth, we replaced waist circumference by BMI; educational status by occupational status or income level; and office systolic blood pressure by office diastolic blood pressure, systolic or diastolic 24-hour ambulatory blood pressure. Last, we analyzed the associations under study stratified by glucose metabolism status to check whether associations of time-and frequency-domain HRV with overall beta cell response were similar over the entire glucose metabolism spectrum (i.e. from normal glucose metabolism to prediabetes and type 2 diabetes).

All analyses were performed with Statistical Package for Social Sciences version 28.0 (IBM SPSS, IBM Corp, Armonk, NY, USA). For all analyses, a P-value < 0.05 was considered statistically significant.