Multi-Tiered Assessment of Gene Expression Provides Evidence for Mechanisms That Underlie Risk for Type 2 Diabetes

Introduction

Type 2 diabetes (T2D) affects approximately 34 million Americans, and an additional 88 million are estimated to be at high risk of developing the disease. ¹ T2D is a significant public health problem due to its increasing prevalence, chronic nature, economic burden, lifestyle-related risk factors, health disparities, long-term consequences, and the need for effective prevention strategies. ¹ One of the primary challenges in preventing and treating T2D is its multifactorial etiology.² Several genetic risk factors for T2D have been reported, however modifiable lifestyle and behavioural characteristics are also important risk factors.^3,4 Recent studies have identified at least 4–5 T2D subtypes that differ in phenotypic characteristics, associated inflammatory biomarkers, and disease outcomes.^5–8

The prevalence of different subtypes of T2D may have associations with genetic risk factors, but it is widely acknowledged that environmental factors, such as lifestyle and behavior, have a significant impact on the development of the condition.⁹ Behavioral interventions may not have the same effect across all subtypes due to these genetic and environmental variations.⁹ Biomarkers that capture the interplay between genes and the environment could be valuable in identifying individuals at risk for T2D and related conditions earlier, allowing for targeted prevention strategies and more effective behavioral interventions. Assessing the transcriptome, which represents the entire set of expressed genes at a given time, can offer insights into how an individual’s genetic makeup responds to environmental factors. The transcriptome is dynamic and can change over time in response to external influences, making it a potentially useful tool for evaluating the combined influence of genetic and environmental risk factors for T2D. By analyzing the transcriptomic profiles, researchers can not only study individual genes but also examine the coordinated activity of groups of genes involved in specific biological processes. These gene co-expression network analysis methods help identify biological pathways and their genetic signatures, which can aid in understanding the underlying mechanisms of the disease.¹⁰

The purpose of this study was to characterize the differential gene expression in biospecimens from underactive adults who participated in a 12-week behavioral intervention for T2D risk reduction in a previously completed randomized clinical trial. Differentially expressed genes were further examined within the context of annotated biological pathways. The goal of this research is to identify the specific genes and biological pathways that contribute to the risk of T2D based on the combination of their genetic characteristics within the context of extrinsic environmental factors.

Materials and Methods

Study Sample and Design

The study sample included a subset of participants from the previously completed Practicing Restorative Yoga versus Stretching for the Metabolic Syndrome (PRYSMS) clinical trial (clinicaltrials.gov identifier NCT01024816), which tested the effects of restorative yoga versus active stretching on fasting blood glucose (FBG) in overweight adults at risk for T2D.¹¹ The trial showed that participants in the restorative yoga group had a 6mg/dL decrease in FBG compared to the stretching group. However, stretching did have a positive effect on FBG for some participants in the stretching group (n = 3; 30%), and therefore the study sample was combined to look at relationships between gene expression and biological pathways in individuals who were responsive to behavioral intervention overall. Participants in the PRYSMS study were recruited from the San Francisco and San Diego areas and met the International Diabetes Federation criteria for metabolic syndrome.¹² The Institutional Review Board at the University of California, San Francisco approved the PRYSMS trial, and all participants provided informed consent in accordance with the Declaration of Helsinki. Only the subset of participants who provided consent to have biospecimens banked for future research were included in the study described in this manuscript. Additional inclusion criteria included availability of banked plasma specimens from the baseline study visit with availability of at least two additional follow-up timepoints (n = 52). Exclusion criteria from the PRYSMS trial are described in detail elsewhere.¹² Briefly, participants were excluded if FBG was ≥126 mg/dl, HbA1c was ≥7.0%, fasting triglycerides were ≥300 mg/dl, weight was ≥400 pounds, participants had neurological conditions that limited their mobility, a hospitalization for coronary heart disease within the past 6 months, current pregnancy or lactation, history of bariatric surgery, substance abuse, and use of medications affecting metabolic factors.¹² Participant demographics, lifestyle characteristics, and medical history were documented at baseline.

Clinical data were collected at baseline, 3-, 6-, 9-, and 12-months. Participant weight was measured on a standard balance beam scale and height using a stadiometer. Waist circumference was measured using a Gulick II tape spring-tension measure at the site of maximum circumference midway between the lower ribs and the anterior superior iliac spine. The mean of two waist circumference measurements was calculated. Body mass index (BMI) was calculated as weight in kilograms divided by height in meter squared.

FBG was measured using an automated analyser with an immobilized enzyme biosensor (YSI 2300 STAT Plus, YSI Life Sciences, Yellow Springs, OH). Total cholesterol, triglycerides and HDL-cholesterol were measured by enzymatic calorimetric methods (Quest Diagnostics, San Jose, CA), and LDL-cholesterol was calculated using the Friedewald equation.¹³ Blood used for banking of plasma was collected by venipuncture into vacutainers containing the preservative EDTA, centrifuged at 4°C to separate plasma from cellular blood components, and stored at −80°C.

For the current analysis and as we have done previously, participants were defined as “responders” or “non-responders”.¹⁴ Participants with at least a 6mg/dL decrease in FBG after 12-months were classified as responders, while participants with <6mg/dL change in FBG at 12 months were classified as non-responders. The approximate physiologic range in FBG for ambulatory adults without type 2 diabetes is 50–126 mg/dL. The mean change in FBG in the intervention group in the PRYSMS trial was 6 mg/dL, which represents an 8% change within this approximate physiologic range. As evidenced by the landmark Diabetes Prevention Program trial, which also identified a 6 mg/dL decrease in FBG, this magnitude of change is clinically meaningful for decreasing incidence of T2D.⁶

Molecular Data Collection

Blood was collected in PAXgene vacutainers containing reagents to lyse cells and stabilize RNA molecules according to standard protocol (Qiagen, Hilden, Germany). Vacutainers were stored at −80°C.

RNA Extraction

RNA isolation was completed using the PAXgene blood RNAeasy kit according to standard protocol (Qiagen, Hilden, Germany). The quantity of RNA was checked by Nanodrop (Thermo Fisher). The integrity of total RNA will be checked on Fragment Analyzer (Agilent, Cat. No. DNF-472), only RNA with RQN number of above 6 was used for library construction.

Library Preparation

A starting quantity of 100ng of total RNA was used according to vendor instructions with Universal plus mRNA with Nu Quant (TECAN, Cat. No. 0520), and the final library PCR amplification was >17 cycles. After library completion, individual libraries were pooled equally by volume and quantified using Fragment Analyzer (Agilent, Cat. No. DNF-474).

RNAseq

Quantified library pools were diluted to 1nM, loaded onto a flow cell (Illumina, Cat. No. FC-420-1001) and sequenced on an Illumina MiniSeq to check for quality of reads. Individual libraries were normalized according to MiniSeq output reads, specifically by the percentage of reads mapping to protein coding genes, loaded onto a flow cell (Illumina, Cat. No. 20028313) and sequenced in four lanes of an Illumina NovaSeq6000.

Statistical Analysis

Descriptive statistics were calculated to examine and evaluate the demographic and clinical characteristics of participants. Comparisons between the baseline characteristics and clinical characteristics of responders and non-responders were assessed using independent groups Student’s t-test for continuous variables and chi-squared tests for categorical variables (R, 2019).

Data Preprocessing

PRSYMS RNA seq Fastq files were transferred to the New York University Greene high performance computational cluster. The fastq files were assessed for quality, and adapters were removed using Fastp (version 0.20.1). The cleaned sequences were mapped to the human genome (build 102) using HISAT2 (version 2.2.1) resulting in Sequence Alignment Map (SAM) files. SAM files were converted into BAM (Binary Alignment Map) files, which were then assembled into transcripts using StringTie (version 2.1.4). A gene count matrix was created using the transcript results and PrepDE python script.

Differential Gene Expression Analysis

Differential expression analysis was performed using the Deseq2 R package (version 1.32) for the responder variable to assess for differences in gene expression between responders vs non-responders. Gene counts generated in an earlier data pre-processing step and the responder phenotypes were converted into Deseq2 objects. Genes that could be annotated to an Entrez ID and with a minimum 10 counts in ≥90% of the sample were retained for the analysis. Using the surrogate variable analysis (SVA) package (version 3.40), batch correction was performed to identify variation not attributed to the phenotype or covariates (a total of 12 surrogate variables were identified and included in the final Deseq2 models). Differentially expressed genes were identified using the filtered SVA-corrected expression set. Genes that had an unadjusted p-value <0.05 were considered significant. The Deseq2 package calculates p-values using the Wald test and adjustment for multiple comparisons was done using the Benjamini-Hochberg false discovery rate method.¹⁵

We used two publicly available datasets to perform validation of the differentially expressed transcripts identified in this study. The first (GSE193436) was sequenced on an Illumina HighSeq 4000 using standard Illumina protocols and filtration.¹⁶ SRA files were downloaded and converted to FASTQ using fasterq-dump (version: 3.0.5) after which the same differential gene expression analysis as described above was applied. This study sample included overweight individuals with (n = 6) and without (n = 6) T2D recruited in China. The second dataset (GSE101931) used the Illumina HumanHT-12 V4.0 expression bead chip array.¹⁷ Log2 normalization and probes that had a median expression in at least 90% of samples were retained. SVA was performed (version 3.40) and included in the subsequent analysis. Differential expression of the filtered normalized data was performed using Limma (version 3.54.0). This study sample included Black individuals, overweight and obese individuals with prediabetes (n = 5) and controls (n = 5) recruited in the US.

Pathway Analysis

Over-representation analysis was performed for the SVA corrected differentially expressed genes using clusterprofiler packages (version 4.0.5) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. Genes that had an unadjusted p-value <0.05 were considered significant. The background signal was defined as all filtered and annotated genes. Dotplot and Gene-Concept network plots were generated using Clusterprofiler.

Weighted Gene Co-Expression Network Analysis

SVA-corrected normalized RNAseq counts were used as the input variables. Standard deviations across the participants for all the genes were calculated. WGCNA (version 1.71) was performed on the genes that had a standard deviation ≥0.3. Modules significantly associated with phenotypes of interest (ie, responder and body mass index (BMI)) were retained. Significant genes were used to perform over representation analysis using clusterprofiler (version 4.0.5) and pathway annotation using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database.¹⁸ Genes that had an unadjusted p-value <0.05 were considered significant. The background signal was defined as all filtered and annotated genes. To identify potential mechanistic interactions between the genes selected for this module, protein–protein interactions were evaluated using STRINGdb (version 2.10.0 R package). Interactions were filtered to include human protein–protein interactions with direct connections.

Pathway Network Analysis

Ingenuity Pathway Analysis (IPA) path explorer (version 01–21-03) was used to identify direct and indirect interactions, in any direction, between genes, proteins, and their products (eg, enzymes, transcription factors, peptidases, etc.) that were identified in the differential expression analysis to be overrepresented in KEGG pathways (p < 0.05). (QIAGEN Inc., https://www.qiagenbioinformatics.com/products/ingenuity- pathway-analysis) Filters in the Ingenuity expert information tool included mammalian, non-chemical, and non-miR interactions.

Results

Participant Characteristics

A total of 52 participants from the PRYSMS trial had complete gene expression data and were included in this study. Participant characteristics are shown in Table 1. Overall, the mean age of the sample was 54 ± 7 years, participants primarily identified as female (70%, n = 36) and White race (73%, n = 38), and were overweight (BMI 35.5 ± 6.8 kg/m²). The only significant difference between the groups was the baseline FBG, which was higher in the responders (n = 10; 112 ± 9 mg/dL) compared to non-responders (n = 42; 100 ± 10 mg/dL) (p < 0.05). Because both groups fall in the prediabetes range, this was not considered a covariate in subsequent analyses.

Table 1 Demographic and Clinical Characteristics by Group

Table 2 KEGG Pathway Themes

Gene Expression

Of a total of 60,675 annotated genes, 17,294 met the quality control thresholds of having ≥10 reads in at least 90% (n = 47) of the sample. There were 570 transcripts that were differentially expressed in responders compared to non-responders (p<0.05)), however none of these remained significant after adjustment for multiple comparisons (q < 0.05) (Supplemental Table 1). The expression of the 10 transcripts with the highest magnitude difference in expression levels between responders and non-responders is shown in Figure 1 (Panel A). While not statistically significant in this sample, Panel B in Figure 1 shows the distribution of gene counts between the responders and non-responders with clear differences between the median values. We then leveraged publicly available external datasets to assess for differential expression of genes that were significant prior to adjustment for multiple comparisons and found that, excluding duplicate probes on the array, were 99 transcripts validated by the external dataset that also used RNAseq methods, and 20 transcripts were further validated in the external dataset that used array-based methods (Supplemental Table 2).

Figure 1 Panel A shows the heat map for the 10 transcripts that showed the highest magnitude fold-change between responders and non-responders. Individual participants are represented by columns; those indicated by purple boxes are the non-responders; those indicated by yellow boxes are the responders. Rows represent individual gene transcripts. Intensity of the colors blue or red indicates the magnitude of the z-score for fold-change. Panel B shows box and whisker plots for the distribution of gene counts for the 10 transcripts with the greatest magnitude differential expression between responders and non-responders. The center thick black bar shows the median, the upper and lower borders of the boxes represent the 75th and 25th percentiles, and individual dots are outliers. Non-responders are shown in purple; responders are shown in yellow.

We also identified transcripts for 33 microRNAs and found three that were significantly differentially expressed in responders compared to non-responders (ie, miR-181, miR-3667, miR-378; p < 0.05); however, none of these remained significant after adjustment for multiple comparisons (q < 0.05).

Pathway Analysis

We identified KEGG pathways that contained over-representation of genes that were differentially expressed between responders and non-responders to behavioral interventions (p < 0.05). The top 10 pathways are shown in Figure 2. In the WGCNA of the responder versus non-responder groups, no significant modules were identified. Because overweight and obesity are significant risk factors for T2D, we also evaluated BMI to include in the WGCNA analysis. For BMI, one significant module was identified (r = 0.35, p < 0.05), containing 99 genes. In an over-representation analysis, 64 KEGG pathways were identified to contain at least one of the 99 genes within this module (top 10 shown in Figure 3). Of the 99 genes in this module, 79 were annotated in the StringDB protein interaction database.¹⁹ A greater number of interactions between these gene products were observed (ie, 427) than expected (ie, 108; p < 0.0001). Among the observed interactions, immune function was a common theme.

Figure 2 The top 10 pathways that contained over-representation of the differentially expressed transcripts related to responders versus non-responders to behavioral interventions are shown. The size of the dot represents the number of transcripts within the pathway that were over-represented relative to the total number of transcripts included in the analyses. The color of the dot shows the magnitude of the false discovery rate-adjusted p-value.

Figure 3 The top 10 pathways that contained over-representation of the differentially expressed transcripts contained within the module derived from weighted gene correlation network analysis for body mass index are shown. The size of the dot represents the number of transcripts within the pathway that were over-represented relative to the total number of transcripts included in the analyses. The color of the dot shows the magnitude of the false discovery rate-adjusted p-value.

Discussion

The purpose of this study was to assess for differential expression of genes between individuals who did versus did not respond to behavioral interventions to decrease risk for T2D. T2D is a complex condition that occurs as a result of numerous underlying mechanisms, which can vary between subgroups of individuals, as previously described.⁵ The potential applications of this research are improved understanding of the specific mechanisms underlying risk for T2D for these subgroups of individuals, which can then inform optimal treatment strategies. This analysis identified 570 transcripts that were differentially expressed between study participants at risk for T2D who responded versus those who did not respond to the interventions. In order to more accurately assess the complex coordinated biological activity that underlies the development of risk for T2D, building on the initial analysis of differentially expressed single genes, we also determined annotated biological pathways that contained the genes observed to be differentially expressed. Previously observed biological pathways were validated, and new pathways were identified.

Among the top ten KEGG pathways targeted by differentially expressed genes (Figure 2), there is a primary theme of neurodegeneration and protein-misfolding disorders, including Huntington's disease, amyotrophic lateral sclerosis, prion disease, and Parkinson’s disease pathways. The nervous system relies on a delicate balance of chemical factors, and the T2D risk factors collectively known as metabolic syndrome can disrupt this balance, linking nervous system damage and neurodegenerative disorders to T2D risk.²⁰ The mechanism of this linkage is complex, including endothelial damage, inflammation, and impaired autophagy.²⁰ Chronic physiological stress-related hormones and inflammation can impede autophagy, and thus stress-reducing lifestyle interventions may enhance autophagy and reduce inflammation.^21,22 Individuals in our sample who responded to stress-reducing lifestyle interventions with decreased FBG may have more resilient nervous system pathways that maintain sufficient autophagy and mitigate inflammation to preserve function. As with most compensatory mechanisms, time of exposure may affect cells’ resiliency. Because of the gradual nature of T2D risk onset, it is difficult to precisely measure the duration of time that individuals in our study have been exposed to the biochemical environment that increases risk to the nervous system cells; thus, we are not able to control for exposure time in our sample.

Individuals diagnosed with Huntington’s disease have seven times greater risk for T2D,²³ and individuals with T2D are at greater risk for Parkinson’s disease.²⁴ In contrast, individuals with T2D are at decreased risk for amyotrophic lateral sclerosis.²⁵ As with these neurodegenerative conditions, T2D has been characterized as a protein misfolding condition when accumulation of the islet amyloid polypeptide (IAPP) hormone results in pancreatic β-cell loss and disease progression.²⁶ Polymorphisms in the IAPP gene cause increased amyloid formation of this protein in the pancreas that begins to accumulate before clinically detectable disease and is associated with early onset T2D.²⁶ Preliminary evidence has suggested the possibility that T2D may be transmitted as an infectious disease via the IAPP protein in the same manner that prions are the vector for transmission of bovine spongiform encephalopathy (ie, mad cow disease).²⁶ Prion diseases in humans are most often developed sporadically, but they may also be genetic or acquired via prion-infected tissue (eg, blood transfusions, grafts, human-derived hormones, or rarely from ingesting infected animals), or contaminated surgical instruments.²⁷ This line of inquiry about a prion-based T2D transmission needs extensive development beyond the current state of evidence, but if true, it may mean that responders to behavioral interventions have an enhanced ability to reject infectious prions, that the prions are dormant or present in smaller amounts that are not yet impacting function, or that responders may be experiencing a different mechanism of T2D risk than the non-responders (ie, prion exposure-related dysfunction vs non-prion-related dysfunction).

While this study did not directly measure IAPP expression, our network analysis (Figure 4) indicated that signal transducer and activator of transcription 3 (STAT3), farnesyl-diphosphate farnesyltransferase 1 (FDFT1), and cytochrome P450 family 51 subfamily A member 1 (CYP51A1), which were differentially expressed between responders and non-responders in this dataset (p < 0.05; q >0.05), or their products, have indirect connections to IAPP. Stat3 has been described as a key mediator in the injury responses in islets transgenic for human IAPP and from individuals with T2D.²⁸ Further, Stat3 activates cytokine signaling-3 (Socs3), an inhibitor of leptin signaling in the brain, with greater time-sensitive activity in the evening before and during T2D observed in mouse models.²⁹ In another mouse model study, high levels of Stat3 led to hyperglycemia, hyperlipidemia, and decreased liver glycogen.³⁰

Figure 4 Solid lines show direct interactions. Dashed lines show indirect interactions. Arrows show activation. Blunted lines show inhibition. Vertical diamonds are enzymes. Horizontal diamonds are peptidases. Vertical ovals are transmembrane receptors. Horizontal ovals are transcription regulators. Circles are other. Double circles are complexes. Vertical rectangles are G-protein coupled receptors. Squares are cytokines. Trapezoids are transporters. Inverted triangles are kinases. Molecules outlined in purple are gene transcripts or their products that were located in the top 10 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways that contained overrepresented genes from differential expression analysis.

High glucose levels led to the upregulation of FDFT1 in a rat model for diabetes, which contributed to altered β-cell differentiation.³¹ Both FDFT1 and CYP51A1 are related to BMI and T2D via increased intracellular cholesterol biosynthesis.^32,33 CYP51A1 may be involved in neurodegenerative disorders by interacting with the Sigma 1 receptor on the endoplasmic reticulum membrane, which under physiological stress may lead to protein misfolding and accumulation.³⁴ Although STAT3, FDFT1, and CYP51A1 are associated with T2D pathology, IAPP, and neurodegenerative diseases, potential mechanisms that connect these pathways are undefined.

In our previous work, we described the KEGG pathways targeted by microRNAs associated with risk for T2D in this same study sample.³⁵ Table 2 compares those previous results to the findings from the current study. The identified pathways were sorted into themes based on function, including metabolism/inflammation and endocrine/hormone. The study described in this paper aligns with our previous findings in that there were similar inflammation-related pathways identified (ie, coronavirus disease – COVID-19, Chemical carcinogenesis – reactive oxygen species, phagosome). Inflammation is widely known to contribute to, and result from, T2D,³⁶ which is consistent with the recurrent theme of inflammation-related pathways being targeted by differentially expressed genes in individuals at risk for T2D. Our findings suggest that responders to a stress-reducing lifestyle intervention may better attenuate inflammation than non-responders, thus improving glycemic control.³⁷ Additional studies are needed to further refine the evidence about the specific inflammatory pathways and to determine whether modification of the expression of inflammatory genes and consequential perturbation of biological pathways may impact risk for T2D. Administration of synthetic microRNAs or antagomiRs has the potential to modify the expression of individual genes, which may be useful for prevention or therapeutic interventions.

The theme of inflammatory pathway involvement is repeated in our findings from the WGCNA module associated with BMI, along with the addition of immune function pathways. KEGG pathways targeted by genes in the WGCNA module associated with BMI (Figure 3) include viral protein interaction with cytokine and cytokine receptor and melanogenesis. Another common theme for pathways containing overrepresentation of genes in the BMI module is related to immune function, including acute myeloid leukemia, thyroid cancer, basal cell carcinoma, and endometrial cancer. Validating our prior observation that pathways related to endocrine/hormone activity are associated with risk for T2D (Table 2), parathyroid hormone synthesis, secretion and action was one of the top KEGG pathways containing genes within the module associated with BMI. There are two genes from the module associated with BMI that are present in seven of the top ten KEGG pathways. These include lymphoid enhancer-binding factor 1 (LEF1) and transcription factor 7 (TCF7). Both of these genes are involved in innate immune activity. One prior study found that polymorphisms in TCF7 are associated with increased prevalence of metabolic syndrome in a sample of individuals living in Pakistan.³⁸ Further studies are needed to determine whether, similar to inflammatory pathways, there may be immunity-related genes that could be targeted by synthetic microRNAs or antagomiRs to decrease risk for T2D and enhance response to lifestyle interventions.

This was a secondary analysis of a previously completed clinical trial to test the effect of behavioral stress reduction interventions on risk factors for T2D. The characteristics of this subset of the overall trial sample did not differ from the overall sample (data not shown). The study sample was predominantly female and White, which is not representative of the overall United States population. The extent to which these findings can be generalized to other sex or gender and racial or ethnic groups is not known. Prior studies have conflated race and ethnicity with genetic ancestry, which may be one explanation for discrepant findings. Our own prior study highlighted that findings from gene expression studies were inconsistently validated when analyzed by sex or gender and racial or ethnic group.³⁹ Future studies warrant the characterization of genetic ancestry associated with genetic polymorphisms that may impact gene expression in a manner that is distinct from the social construct of race and ethnicity.

Not all participants agreed to storage of biospecimens, and among those who did, not all participants had sufficient biospecimens for inclusion in this study, which resulted in a relatively small sample size that may have impacted the power to detect statistical differences in gene expression. However, as shown in Figure 1 (Panel B), there is a distinct difference in the distribution of the gene counts between the responders and non-responders. This limitation was mitigated by the approach of evaluating biological pathways and networks in addition to single-gene transcripts. Furthermore, we were able to validate our findings using external datasets, even relatively small sample sizes (n = 10, n = 10) that included related but not identical phenotypes.

The purpose of this study was to identify genes and pathways associated with responses to behavioral interventions to decrease risk for T2D. Because single genes are numerous and, individually, capture only a fragment of overall physiological function, there may be limitations to how much information can be gleaned from gene-level assessment, particularly in the setting of a relatively small sample size. Assessment of a higher order of biology, through evaluation of over-representation of differentially expressed genes within annotated biological pathways, may overcome these limitations. Validating prior studies, we identified inflammation, immune function, and metabolic function and common themes for pathways that contained differentially expressed genes in this sample. We also identified a strong theme of neurodegenerative and protein misfolding disease; a small number of recent studies have also characterized T2D as a protein misfolding condition. Additional studies are needed to determine if modulation of expression of genes within these pathways by synthetic microRNAs or antagomiRs could impact risk for T2D and responses to risk-reduction interventions. In addition, further studies are needed to determine whether the observed themes are generalizable to a population with more diverse biological and social characteristics.