BACH2 inhibition reverses β cell failure in type 2 diabetes models

Type 2 diabetes (T2D) is characterized by pancreatic β cell function and mass abnormalities (1). The clinical features of β cell failure include limited adaptive potential of β cell mass (2, 3), increased glucagon tone (4), rapid progression in response to incipient hyperglycemia (5), and limited response to treatment (6, 7). In this regard, whether β cell failure is reversible remains arguably the key question bearing on the issue of disease modification (8). Understanding islet cell dysfunction is critical to design mechanism-based, durable, and safe disease-modifying interventions (9). Yet, the limited ability to access the endocrine pancreas in vivo has hindered our ability to functionally interrogate islet cells. To wit, despite intensive, multicenter efforts to standardize procurement and analysis of human islets, there is no consensus on a diabetic gene expression signature (10), and mechanistic understanding of β cell failure pathways is still lacking. Indeed, while these studies have demonstrated islet cell gene expression heterogeneity (11–14), clear definition of a diabetic islet cell signature remains elusive.

Single-cell-based expression profiling can reveal disease determinants, but its ability to detect subtle differences, such as those between a diabetic and nondiabetic cell type (15), is hampered by the fact that differential gene expression may not be assessed for up to 80% of the genes (gene dropout), due to low sequencing depth (16). To address this problem, we have developed the metaVIPER algorithm — an extension to single cells of the VIPER algorithm (17) — which can accurately measure the activity of approximately 6,500 regulatory and signaling proteins in single cells, based on the enrichment of their transcriptional targets (regulon) in differentially expressed genes (18). metaVIPER can accurately identify rare cell subpopulations that otherwise escape detection by expression-based analyses, by measuring protein activity in single cells (19).

In this work, we applied metaVIPER to identify master regulatory (MR) proteins controlling the potentially noisy gene expression signature of diabetic islet cells. In a critical validation of the computational analysis, we show that single-cell-based, CRISPR-mediated inhibition of metaVIPER-inferred MR proteins, such as BACH2, can effectively reverse T2D-specific β cell features. Moreover, immunohistopathology of diabetic pancreata demonstrates an enrichment of BACH2-expressing cells, as predicted by the metaVIPER algorithm. Finally, we show that treatment with a BACH inhibitor improves disease phenotypes in diabetic mice and insulin secretion in diabetic islets of mice and humans, thus providing proof of concept of the potential translational value of the findings. Thus, single-cell protein-activity-based analyses can discover diabetes-causing phenotypes, identify their MR proteins, and guide their pharmacological targeting. This work highlights the clinical potential of disease modification as a conceptual advance in diabetes research.

Single-cell protein activity analysis of pancreatic islet cells. To capture T2D-specific cellular states, we used single-cell RNA sequencing (scRNA-Seq) profiles from 4 nondiabetic controls (ND) and 6 T2D donors (Supplemental Table 1; supplemental material available online with this article; https://doi.org/10.1172/JCI153876DS1) to build islet-specific transcription regulatory networks and to identify the transcription factors (TFs) and cofactors (co-TFs) that control their differential gene expression signatures (Figure 1A). To measure differential TF activity in individual cells, we first used ARACNe — a highly validated reverse-engineering algorithm (20) — to build tissue-specific regulatory networks from single cells, as shown in Ding et al. (18). Use of single-cell data greatly exceeded ARACNe’s requirements for large data sets (n ≥ 100 samples), resulting in a comprehensive, first-of-a-kind pancreatic islet cell regulatory network, assembled from 6,137 scRNA-Seq profiles from T2D or ND islets. The resulting network is composed of regulons for 1,813 TFs, 969 co-TFs, and 3,370 signal transduction (ST) proteins. Since direct physical regulation of gene expression signatures is implemented by TFs and co-TFs, we used the network to measure the differential activity of TF and co-TF proteins in each single cell, compared with the average of all single cells, using metaVIPER (ref. 18 and Figure 1A).

Figure 1

metaViper analysis buffers donor-to-donor variability. (A) Schematic workflow of islet-specific regulatory network generation (ARACNe) and protein activity analyses (metaVIPER) at the single-cell level from sc-RNA-Seq data. (B) Single cells from human ND and T2D islets were projected onto 2D t-SNE space based on protein activity inferred from islet-specific regulatory networks. Each dot represents a single cell, color coded according to donor. (C) INS (insulin) and GCG (glucagon) mRNA expression plotted in t-SNE at the single-cell level. RPM, reads per million mapped reads. (D) Computationally inferred MAFA and IRX2 protein activity plotted in t-SNE at the single-cell level.

t-SNE (t-distributed stochastic neighbor embedding) analysis based on metaVIPER-measured protein activity effectively removed most donor-to-donor batch effects, resulting in 2 distinct, yet donor-independent, cell clusters (Figure 1B and Supplemental Figure 1). The 2 clusters did not represent the 2 major canonical islet populations (β and α cells), as shown by integrative assessment of INS (insulin) and GCG (glucagon) mRNA levels (Figure 1C), as well as MAFA and IRX2 protein activity (Figure 1D). Yet, INS and GCG mRNA expression aligned perfectly with metaVIPER-measured MAFA and IRX2 activity, respectively, confirming metaVIPER’s accuracy in marker activity measurement.

Statistically significant superiority of protein activity versus gene-expression-based clustering was recently demonstrated (21), including at the single-cell level (22). Thus, to further cluster islet cells, including both ND and T2D cells, we performed unsupervised hierarchical clustering analysis using iterative clustering (iterClust) algorithms (23) with the full, metaVIPER-inferred protein activity matrix. The analysis yielded 11 clusters (Figure 2A and Supplemental Figure 2), which were further characterized by assessing established lineage and functional markers, including (a) hormone mRNA expression, (b) activity of TFs characteristic of either β or α cell identity (24), (c) activity of metabolic inflexibility/stress-response drivers (25), and (d) activity of endocrine progenitor and stem-like-cell TFs (ref. 26 and Supplemental Table 2). The latter were chosen based on experimental evidence from animal studies of their role in disease progression (25, 27, 28).

Figure 2

iterClust classifies ND and T2D islet cells into different biological states. (A) iterClust analyses performed using ND and T2D islet cells. The resulting cluster architecture was visualized as a heatmap. Each subgroup is color coded. Each bar denotes a single cell. Black bars represent T2D cells and white bars ND cells. KRT19, AMY2A, PPY, SST, GCG, and INS mRNA expression is plotted at the single-cell level. metaViper-inferred protein activity for β or α cell factors, metabolic inflexibility, endocrine progenitor, and stemness markers is plotted at the single-cell level. RPM, reads per million mapped reads. (B) Violin plots showing the distribution of cells in each cluster based on integrated activity of β cell factors (top), α cell factors (middle), or stemness markers (bottom). (C) 3D plot showing integrated β cell factor, α cell factor, and stemness activity on the x, y, and z axes, respectively, at the single-cell level. (D) 3D plot as in C but based on the average cell behaviors of each cluster. Red arrows indicate T2D-enriched clusters, MI+² (T2D-β-like) and MI–⁴ and MI–⁵ (both T2D-α-like).

Protein-based cluster analysis identified 2 main clusters, characterized by differential activity of metabolic inflexibility/stress-response markers (25, 27) (MI+ and MI–) — a stage in the progression of β cell failure (25) — including PPARα and PPARγ (Figure 2A and Supplemental Figure 3). Expression of the corresponding genes was virtually undetectable in most cells (Supplemental Figure 4), thus confirming metaVIPER’s superior sensitivity and complementarity, compared with mRNA measurements (Supplemental Figure 5). Regulatory proteins, including TFs and co-TFs, are often expressed at low levels in terminally differentiated cells (29). More importantly, their mRNA expression is a poor predictor of their activity, due to complex posttranscriptional and posttranslational events. scRNA-Seq data have low sequencing depth, and as a result most genes (generally >80%) cannot be detected by even a single read (dropouts; ref. 18). metaVIPER effectively overcomes this limitation, because the differential activity of a protein is measured from the differential expression of 50 or more of its transcriptional targets, akin to a highly multiplexed gene reporter assay, thus allowing reliable protein activity assessment even from low-depth profiles (18, 22).

MI+ cells showed clear evidence of a metabolic stress response (30), denoted by high FOXO1, HIF1α, HSF1, and TP53 activity (Supplemental Figure 3, A and B); endocrine progenitor features, denoted by high RFX6 (31, 32) and RFX7 activity (Figure 2A and Supplemental Figure 3C); and cell stemness, denoted by high NANOG, MYCL, and POU5F1 activity in a subset of RFX6/7-positive MI+ cells (Figure 2A and Supplemental Figure 3D). Thus, the 2 main functional clusters in ND and T2D islets are defined by differential activity of metabolic inflexibility/stress-response and endocrine progenitor/stem-cell-like markers (25, 27).

Iterative clustering further divided the MI+ cluster into 5 distinct subclusters — labeled MI+¹ to MI+⁵ — and the MI– cluster into 6 distinct subclusters — labeled MI–¹ to MI–⁶. However, MI–⁶ appeared to contain KRT19-positive pancreatic ductal cells and was thus excluded from further analyses. Four clusters (MI+², MI+⁴, MI–¹, and MI–²) displayed β cell features. Specifically, MI–¹ contained the healthiest β cells, as indicated by the highest levels of INS expression and β cell–specific TF activity, as well as by inactive metabolic inflexibility/stress response, progenitor/stem cell, and α cell identity (Figure 2, A–D); MI–² had similar, yet less pronounced features, while MI+² and MI+⁴ differed from the other 2 clusters based on activation of metabolic inflexibility/stress-response and progenitor/stem cell features (Figure 2, A–D, and Supplemental Videos 1 and 2).

In contrast, 6 clusters (MI+¹, MI+³, MI+^5, MI–³, MI–⁴, and MI–⁵) had α cell or α cell–like features. Specifically, MI+⁵ showed the strongest α cell identity, with high GCG expression and IRX2 activity, and slight activation of metabolic inflexibility/stress response (Figure 2, A–D); MI+³ had similar features but stronger activation of metabolic inflexibility/stress response, while MI+¹, also showed high progenitor/stem cell activities (Figure 2, A–D, and Supplemental Videos 1 and 2). Finally, MI–⁵, MI–⁴, and MI–³ showed gradually decreasing α cell activity, despite preserved GCG expression and absence of metabolic inflexibility/stress response or progenitor cell features (Figure 2, A–D).

Single-cell protein activity analysis reveals T2D-enriched islet subpopulations. Next, we assessed whether the 10 molecularly distinct clusters were differentially represented in T2D versus ND cells and evaluated the statistical significance of their association (Figure 3A). Indeed, t-SNE plots revealed highly distinctive cluster representation patterns and differential subpopulation representation (Supplemental Figures 1 and 2). Notably, cluster analysis was performed without separating ND and T2D cells and using a common, islet-specific regulatory network, thus allowing for the unbiased, unsupervised determination of T2D- versus ND-specific cells across detected clusters.

Figure 3

T2D-enriched subclusters in human islets. (A) Bar plots presenting the percentage of ND or T2D cells in each subcluster. A dashed line represents the proportion of cells from ND or T2D islets analyzed by scRNA-Seq. P values were derived from Fisher’s exact test. (B) Single cells from human ND and T2D islets were projected onto 2D t-SNE space based on protein activity inferred from islet-specific regulatory networks. ND cells are shown in gray and T2D cells are shown in red as background. T2D cells in T2D-enriched clusters, MI+² (T2D-β-like) and MI–⁴ and MI–⁵ (both T2D-α-like), were color coded as indicated.

The analysis revealed that ND islets were highly enriched in “healthy” β cells (MI–¹), and α-like cells characterized by weak α cell markers but inactive metabolic inflexibility/stress response and progenitor/stem cell features (MI–³) (P = 3.4 × 10^–8 and 8.2 × 10^–10, respectively, by Fisher’s exact test). In contrast, T2D islets were highly enriched in 3 clusters: MI+², MI–⁴, and MI–⁵ (P = 9.6 × 10^–13, 7.9 × 10^–71, and 1.0 × 10^–80, respectively; Figure 3, A and B). Among these, MI+² is composed of β cells with an active metabolic inflexibility/stress response and progenitor/stem-cell-like features (Figure 2, A–D, and Supplemental Videos 1 and 2). These cells are distinguished from the β cell cluster MI–² by their higher MAFA and metabolic inflexibility/stress-response activities (Figure 2, A and B). The biphasic distribution of MAFA activity among different β cell clusters is reminiscent of the rodent islet stress response associated with diabetes progression (30). However, this also heralds the onset of metabolic inflexibility (PPARα/γ; ref. 25) and dedifferentiation (27), as indicated by RFX6 activation (31). MI–⁴ and MI–⁵ combine weak α-like features with a muted stress response. However, MI–⁵ showed elevated β cell and stemness TF activities compared with MI–⁴ (Figure 2, B and D, and Supplemental Videos 1 and 2), consistent with a putative β- to α-like-cell reprogramming. Although we cannot establish directionality of the transition based on these data alone, functional testing described below supports the β- to α-like model. In summary, this analysis identified endocrine cell subpopulations enriched in T2D. They include β cells characterized by metabolic inflexibility/stress response and progenitor/stemness signatures (MI+²), as well as cells with α-like features (MI–⁴ and MI–⁵), consistent with cell reprogramming.

Different driver proteins elicit distinct T2D cell state transitions. To identify candidate MR proteins that mechanistically regulate the transcriptional identity of cells in T2D-enriched clusters, we assessed the most activated and inactivated proteins (candidate MR proteins), based on metaVIPER analysis of the 2 most representative T2D-enriched clusters, MI+² (T2D-β like) and MI–⁵ (T2D-α like), versus the healthiest ND β and α cells (MI–¹ and MI+⁵, respectively; Supplemental Data 1).

Based on their activity in T2D-enriched clusters and association with T2D susceptibility loci (Supplemental Table 3), we selected 15 candidate MR proteins for experimental validation using gain-of-function (GOF) assays in ND islets (Supplemental Table 4). For this purpose, we modified the Perturb-Seq protocol (33, 34) (single-cell GOF sequencing [scGOF-Seq]) (Figure 4A), in which the transduced gene can be mapped to an individual cell using a DNA barcode, allowing its identity and the effect of the gene manipulation to be read out by scRNA-Seq, for further metaVIPER analysis (Supplemental Figure 6 and Supplemental Table 3). We transduced islets with pooled adenoviruses (n = 3 in each pool) to simultaneously screen multiple candidates and test potential epistatic interactions among the cotransduced TFs.

Figure 4

scGOF-Seq experiments in human islets. (A) Schematic drawing of scGOF-Seq plasmids and experimental procedure. ND islets (ND5) were used for this experiment. (B) Violin plots showing the distribution of cells following transduction with each individual candidate or combination thereof analyzed according to T2D-β-like signature, an integrated value of RFX6, RFX7, FOXO1, PPARα, PPARγ, RB1, POU5F1, NANOG, and MYCL protein activities. Nontransduced and BFP-transduced ND islets serve as negative controls. (C) Bar plots showing the normalized proportion of islet cells with a positive T2D-β-like signature (activity > 0) in each scGOF-Seq condition over nontransduced control. (D) Violin plots showing cells with a T2D-α-like signature, which is an integrated value of IRX2, ZNHIT1, ZFPL1, PAX6, and DRAP1. (E) Bar plots showing the normalized proportion of islet cells with a T2D-α-like signature as in B. P values in C and E were derived from Fisher’s exact test.

We measured the MR candidates’ ability to reprogram ND cells toward a T2D transcriptional identity, as represented by the MI+² (β-like with metabolic inflexibility/stress response and progenitor/stem-like features) and MI–⁵ (putatively converted, α-like) subtypes (T2D-β-like and T2D-α-like in Figure 2D). As an integrated metric of the MI+² (T2D-β-like) state, we used a protein activity signature composed of RFX6, RFX7, FOXO1, PPARα, PPARγ, RB1, POU5F1, NANOG, and MYCL — henceforth referred to as T2D-β-like signature. As a metric of the MI–⁵ (T2D-α-like) state, we used a protein activity signature composed of IRX2, ZNHIT1, ZFPL1, PAX6, and DRAP1 — henceforth referred to as T2D-α-like signature (Figure 4). Results are presented as violin plots illustrating the shift from the basal (nontransduced or BFP-transduced) to the scGOF-Seq cell population, and as bar graphs (Figure 4 and Supplemental Figure 7). For some candidate MR combinations, the number of transduced cells was too low to allow statistical evaluation.

The strongest emergence of the T2D-β-like signature was observed in cells transduced with the BACH2/TSZH2 and RARB/GAS7/ZNF385D combinations, as well as with ZRANB3/MYT1L and CUX2/RFX7 (Figure 4, B and C). Next, we asked which TFs can drive conversion to a T2D-enriched α-like state. AFF3 showed the strongest effect, followed by the BACH2/TSHZ2 and ZRANB3/MYT1L combinations (Figure 4, D and E; see also Supplemental Figure 7). These data suggest that active reprogramming to an α-like-cell state is characteristic of T2D.

To determine the mechanisms driving cell fate transitions to a T2D-β-like signature, we analyzed TFs and co-TFs whose activity changed in coordinated fashion upon conversion. We found that BACH2, FOXO1, MYTL1, NFATC3, RFX7, and TCF4 were coregulated in all conditions associated with a T2D-β-like endpoint signature (Supplemental Figure 8). While FOXO1 or TCF4 ectopic expression did not drive the T2D-β-like signature, their activity significantly increased upon BACH2- or ZRANB3-dependent conversion, suggesting that these 6 TFs form a hierarchical module whose coordinated activity drives T2D-β-like conversion. This observation is consistent with RNA-Seq analyses of β cell–specific FoxO1-knockout mice, which showed decreased Bach2 mRNA levels (25, 28). We confirmed the data in a replication experiment using islets from a different donor (Supplemental Figure 9), obtained on different dates with 2 independent adenovirus infection experiments. We present the 2 experiments separately (Figure 4 for ND5 and Supplemental Figure 9 for ND6; Supplemental Table 4) to circumvent batch effects that may arise from a merged analysis. Each experiment was internally controlled using negative control guide RNAs (gRNAs), and the two are statistically independent.

Reversibility of diabetic islet cell signature by gene perturbation in T2D islets. Next, we asked whether targeted, CRISPR-mediated ablation of selected MR activity could reverse disease features in T2D islets (Supplemental Table 4), using a modified Perturb-Seq approach (ref. 34 and Supplemental Figure 10A). We generated individual adenoviruses encoding 4 different gRNAs for each of the top 3 candidate MRs emerging from the GOF assays, including AFF3, BACH2, and CUX2. Following transduction, we measured the extent of cotransduction with each gRNA by decoding gRNA barcodes (Supplemental Table 5) and plotted the effect of different gRNA combinations for each candidate (MR) as SET1 to SET3, where each SET corresponds to a different gRNA for each gene (Supplemental Figure 10B). We measured mRNA levels and computed protein activity signatures to identify sets that reproducibly reprogrammed cells to an ND-like state (Figure 5A, blue color).

Figure 5

Perturb-Seq experiments in human T2D islets. (A) Heatmap showing protein activity and mRNA expression of master regulators (MRs) in Perturb-Seq analyses. Black arrows indicate the strongest knockout efficiency in each Perturb-Seq condition. T2D islets (T2D7) were used for this experiment. (B) Bar plots showing the normalized proportion of islet cells with a positive T2D-α-like signature, an integrated value of IRX2, ZNHIT1, ZFPL1, PAX6, and DRAP1 (activity > 0) in each Perturb-Seq condition. (C) Bar plots showing the normalized proportion of islet cells with a T2D-β-like signature, which is an integrated value of RFX6, RFX7, FOXO1, PPARα, PPARγ, RB1, POU5F1, NANOG, and MYCL protein activities.

AFF3 SET2 decreased AFF3 mRNA abundance and converted T2D-α-like cells back to an ND-like state but had no effect on the T2D-β-like signature (Figure 5, B and C). BACH2 SET1 decreased BACH2 mRNA and reversed both T2D-α-like and T2D-β-like signatures toward an ND state. Finally, CUX2 SET2 inhibited CUX2 mRNA and reversed T2D-α-like and T2D-β-like signatures. In each instance, the extent of knockdown correlated with the effects (Figure 5, B and C). These assays strongly corroborate the findings from the scGOF studies.

We also separately analyzed restoration of β and α cell features using INS and GCG mRNA as an indirect measure (Supplemental Figure 11A), as well as MAFA and IRX2 activities (Supplemental Figure 11, B–E). This subanalysis was consistent with the findings in Figure 5.

Functional effects of MRs on fate-converted cells. We next tested whether converting ND β cells (Supplemental Table 4) to a T2D-enriched signature impairs β cell function, as assessed by measuring intracellular Ca²⁺ flux in response to glucose in single β cells identified by the INS reporter. For these experiments, we selected BACH2 (due to its effects on T2D-β-like signatures and genetic link to diabetes susceptibility; ref. 35), AFF3 (due to its effects on T2D α-like cell conversion), and TCF4 (as a negative control). We gated β cells by cotransducing RIP-zsGreen together with adenovirus expressing each candidate TF (36) to monitor Ca²⁺ flux (Figure 6A). β Cells transduced with BFP showed a brisk Ca²⁺ response to glucose or KCl (Figure 6B). The glucose response was blunted in cells transduced with AFF3 or BACH2, while the KCl response was intact, indicating that AFF3 or BACH2 GOF impairs β cell glucose sensing but do not alter membrane depolarization (Figure 6, C and D). This is consistent with the loss of β cell features observed in scGOF-Seq experiments (Figure 4D). In contrast, TCF4 GOF had no statistically significant effect on Ca²⁺ influx compared to Tag-BFP (Figure 6E), consistent with the scGOF-Seq results. These data further validate our approach and strengthen the observations on AFF3 and BACH2.

Figure 6

Single–β cell calcium microfluorimetry. (A) Schematic drawing of the single-cell Ca²⁺ imaging procedure. (B–E) Representative traces of Ca²⁺ flux measured by Rhod-2 loading in Ad-BFP– (B), Ad-AFF3– (C), Ad-BACH2– (D), or Ad-TCF4–transduced (E) primary human β cells from 3 technical repeats using 2 ND donors (ND7 and ND8). Red arrows indicate the timing of addition of 16.8 mM glucose, and black arrows indicate addition of 40 mM KCl. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001 by ANOVA with Dunnett’s method.

Increased BACH2 expression in human T2D pancreata. To test the clinical relevance of these observations, we performed immunofluorescence analysis of pancreatic samples obtained from ND and T2D patients (Supplemental Table 6) with anti-BACH2 antibodies. BACH2-immunoreactive cells are rare in ND islets, while T2D pancreata show strong BACH2 immunoreactivity both in α and β cells (Figure 7, A and B). Quantitative analysis of the results showed a 4-fold increase in BACH2-positive β cells and a 3-fold increase in BACH2-positive α cells in T2D (Figure 7, C and D, and Supplemental Table 6). Also, the percentage of BACH2-positive cells among bihormonal cells (INS⁺GCG⁺) increased significantly in T2D islets (22% ± 5.6% vs. 80% ± 4.6%, ND vs. T2D, P < 0.0001; Figure 7E and Supplemental Table 6). These data are consistent with metaVIPER, showing activation of the BACH2 program in association with dedifferentiation and β to α transition.

Figure 7

BACH2 immunostaining in human pancreata. (A and B) Coimmunostaining of ACH2 (green), insulin (white), glucagon (red), and DAPI (blue) in human ND (A) or T2D pancreata (B). INS⁺GCG⁺ bihormonal cells expressing BACH2 are indicated by the pink arrows. Orange arrows indicate BACH2-positive β cells and white arrows indicate BACH2-positive α cells. Scale bars: 20 μm. (C–E) Ratio of BACH2-positive cells to INS⁺GCG^– β cells (C), GCG⁺INS^– α cells (D), and INS⁺GCG⁺ bihormonal cells (E). All data are expressed as mean ± SEM. **P < 0.01, ***P < 0.001 by 2-tailed, unpaired Student’s t test performed between the 2 groups (n = 9 per group).

A BACH2 inhibitor improves glycemia in diabetic mice. Finally, we tested whether inhibition of BACH2 using BACH2 inhibitors affects diabetes in experimental animals. Bach2 mRNA expression was significantly increased in diabetic mouse β cells compared with controls (Supplemental Figure 12A), consistent with human data (Figure 7). To test the effect of BACH2 inhibition, we administered compound 8 (Supplemental Figure 12B and ref. 37), a dual BACH2/BACH1 inhibitor (Supplemental Figure 12, C and D) to db/db mice for 2 weeks and monitored glucose levels throughout. Starting on day 7, we observed a glucose-lowering effect (Figure 8A) as well as decreased fasting glucose in compound 8–treated animals in a dose-dependent manner (Figure 8B). Whereas plasma insulin levels dropped in vehicle-treated mice, consistent with progressive β cell failure, this trend was reversed in compound 8–treated mice, so that insulin levels gradually increased, consistent with improved β cell function upon BACH2 inhibition (Figure 8C). Moreover, oral glucose-tolerance tests showed that compound 8 treatment improved glucose tolerance in db/db mice compared with vehicle-treated controls (Figure 8D and Supplemental Figure 12, E and F). It is possible that some of these effects are related to peripheral actions of compound 8. Thus, to analyze the direct effect on islets, we examined the effect of compound 8 on diabetic islets isolated from db/db mice or obtained from diabetic organ donors (Supplemental Table 4). In both instances, treatment with compound 8 resulted in a 50% increase in glucose-stimulated insulin secretion and total insulin content, consistent with restoration of β cell function (Figure 8, E and F, and Supplemental Figure 13). These data are consistent with the Perturb-Seq data in human islets (Figure 5) and support the conclusion that BACH2 inhibition can reverse β cell failure.

Figure 8

The effect of BACH2 inhibitor treatment in diabetic mouse models. (A) Nonfasting blood glucose levels in male db/db mice treated with compound 8 or vehicle control throughout 2 weeks of treatment. (B) Overnight-fasted glucose levels in db/db mice after 2 weeks of compound 8 treatment. (C) Plasma insulin levels in db/db mice treated with compound 8 as in A. (D) Oral glucose-tolerance test (OGTT) in db/db mice treated with compound 8 or vehicle control (n = 8 mice per group). (E) Insulin secretion at 2.8 mM or 16.8 mM glucose in islets isolated from db/db mice following 24 hours of treatment with 1 μM compound 8 (n = 10 per group). (F) Insulin secretion as in E but in hT2D islets (n = 10 per group, T2D8 and T2D9). All data are expressed as mean ± SEM. *P < 0.05; **P < 0.01; ***P < 0.001 by ANOVA with Dunnett’s method performed between the 2 groups.

Clinical progression of T2D is characterized by rapid degradation of β cell function (5) that becomes refractory to treatment, requiring combination therapy and resulting in life-threatening complications. GLP-1–based agents have β cell protective effects (38), but do not appear to significantly delay monotherapy failure, or reverse established β cell failure (39, 40).

Through identification and functional characterization of diabetes-enriched islet cell populations, we provide proof of principle of disease reversibility at the cellular level and validate it by cellular studies in primary β cells of diabetics, immunohistochemistry, and animal experimentation. These observations dovetail with the clinical reversibility of β cell failure in the early phases of T2D (1), or in response to diet (41, 42), but also add a pharmacological dimension to the concept. Compound 8 was originally developed as a BACH1 inhibitor to treat multiple sclerosis but our data showed that it inhibits both BACH1 and BACH2 (Supplemental Figure 12). BACH1 and BACH2 are highly homologous, with significant similarity in the bZIP and BTB domains. However, their expression pattern is different. While BACH1 is expressed ubiquitously, BACH2 is expressed in in B cells, T cells, alveolar macrophages, and neural cells (43). Our study suggests that BACH1/2 inhibitors (44, 45) may be repurposed in exploratory animal studies to evaluate their effectiveness against β cell failure. BACH2 expression has been linked to genetic susceptibility to type 1 diabetes (35), and changes to its chromatin structure have been reported in diabetic islets (46). In addition, BACH2 is regulated by the Akt/mTOR pathway and has emerged from a previous analysis of β cell dedifferentiation (28). However, it has not been identified as a mechanistic determinant of T2D cell states. BACH2 binds to the Maf response element (MARE) with small MAFs, and inhibits genes involved in antioxidant responses, which are activated by NRF2 when BACH2 is absent (43). Therefore, inhibition of BACH2 in T2D can increase NRF2-dependent antioxidant response genes, decreasing cell damage and increasing β cell function.

Unbiased analysis of cell transition states identified 3 salient features of diabetic cell populations: (a) α cell–like features, (b) metabolic inflexibility/stress response, and (c) endocrine progenitor/stem-cell-like features. A striking feature of the MR analysis is the existence of intermediate states, reflecting continuity between the main β and α cell identity, as detected by identification of α cell subpopulations with increased β cell signature and vice versa. This can be reconciled with the observation that the latter can be reprogrammed into surrogate β cells (47), as well as with the overall plasticity of endocrine islet cell fate. Within this distribution, healthy β cells (cluster MI–¹) are enriched in ND, while “stressed” β cells (i.e., β cells with markers of metabolic inflexibility/stress response/progenitor/stemness) are enriched in T2D (cluster MI+²). In the latter cluster, activation of the 2 nuclear hormone receptors PPARα and -γ may appear counterintuitive, since they oversee different lipid metabolic functions, but is in fact consistent with the observation that, as islets fail, synthetic and oxidative branches of lipid metabolism are coactivated (25). We have suggested that activation of PPARγ is in this context compensatory by directing lipids toward synthesis as opposed to oxidation. This provides a potential explanation for the longstanding clinical observation that PPARγ agonists have beneficial effects on diabetes prevention and β cell function (48). Nevertheless, given the side effects associated with either PPARγ agonism (49, 50) or PPARα antagonism (51), we do not believe that these data can be exploited pharmacologically. Rather, they provide unbiased validation of the metaVIPER approach by linking it to mechanisms of disease identified in experimental models (25, 52).

Within “stressed” clusters (MI+), T2D islets are enriched in stressed β cells (MI+²), whereas ND islets are enriched in stressed α-like cells (MI+³). In this regard, there appears to be an inverse correlation between GCG mRNA and activation of the IRX2-dependent network on the one hand, and metabolic inflexibility/stress response on the other. Cells in which these features are more pronounced prevail in ND islets, whereas cells with more muted α-like characteristics are enriched in T2D (MI–³ vs. MI–⁴ and MI–⁵), which is somewhat reminiscent of mixed-feature cells identified previously (53). What can this mean? One potential interpretation is that α-like cells displaying a more robust stress response are more resilient, and hence capable of reversal to β cells (if indeed this is their origin), whereas those lacking a stress response are transitioning toward a functionally quiescent state. To the extent that these processes are reversible, at least up to a point in the progression of the disease (54), there are likely to be subpools of cells undergoing β to α and α to β cell transitions. However, the clinically salient features of decreased insulin production and increased glucagon production make the β to α cell transition a key disease process.

Our functional studies show that analysis of regulatory networks, reverse engineered from single cells of pancreatic islets from T2D and ND individuals, can pinpoint drivers of different aspects of the cellular pathology, and, consistent with our findings in other contexts (55), favor a model in which a small subset of TFs can drive different aspects of diabetic islet pathology. This regulatory model and associated protein activity signatures differ from those identified thus far and could not have been predicted based on mRNA expression, GWAS, or epigenetic analyses (56). The results that emerge are that BACH2 strongly affects dedifferentiation, whereas AFF3 has the strongest effect on driving the α-like cell phenotype.

In conclusion, the key findings of this work are the functional assessment of unbiased scRNA-Seq data, their functional validation in a disease-relevant context, and the identification of actionable targets for intervention that address the underlying pathophysiology.

Human islet procurement. Human ND and T2D islets were from the NIH’s Integrated Islet Distribution Program (IIDP). Upon arrival, human islets were plated at a density of 10,000 islet equivalents per 10-cm nontreated tissue culture dish (Corning, 430591) into 10 mL of islet culture medium (Prodo Labs, PIM(S), cPIM-CS001GMP), supplemented with 5 mL PIM(G) Glutamine/Glutathione (Prodo Labs, PIM-G001GMP), 5% PIM(ABS) Human AB Serum (Prodo Labs, PIM-ABS001GMP), and triple antibiotics, PIM(3×), which includes ciprofloxacin (Corning, 61-277-RF, 10 mg/L), gentamicin (Sigma-Aldrich, G1272, 10 mg/L), and amphotericin B (Omega, FG-70, 2500 μg/L). Islets were cultured for no longer than 1 week after arrival and medium was replaced every 2 days. On the day of scRNA-Seq, islets were collected, washed in phosphate-buffered saline (PBS) once, and dispersed into single cells by mechanical shaking at 37°C using 0.05% trypsin (Gibco, 25300054). For scGOF-Seq or Perturb-Seq, adenovirus-transduced cells were further gated for PacBlue-positive cells at an excitation wave length of 401 nm, and collected at 452 nm.

scRNA-Seq library preparation using the Fluidigm C1 800 platform. Sorted cells were suspended in C1 Cell Suspension Reagent (Fluidigm) and loaded onto each inlet of the C1 high-throughput integrated fluidic circuit (HT IFC). The number of cells captured in each chamber was visualized and noted using a phase contrast microscope. Only chambers with single-cell capture were used for analysis. Cells were lysed, and mRNA reverse transcribed and PCR amplified using C1 Single-cell Auto Prep IFC (Fluidigm, protocol 101-4964). The quality and yield of cDNA were determined by Agilent Bioanalyzer using Agilent High Sensitivity DNA Chip. Libraries for sequencing were prepared using Nextera XT DNA library preparation kit (Illumina FC-131-1096) and sequenced with 50 paired-end cycles on an Illumina HiSeq 2500. Each library pool was sequenced in 2 lanes of the Illumina HiSeq 2500.

10× Genomics platform scRNA-Seq library preparation. Sorted cells for scGOF-Seq or Perturb-Seq were treated with a Chromium Single Cell 3′ Library and Gel Bead Kit v2 (PN-120237), Chromium Single Cell 3′ Chip Kit v2 (PN-120236), and Chromium i7 Multiplex Kit (PN-120262), and analyzed with a 10× Genomics Chromium for Single-Cell Library Preparation Instrument, per the manufacturer’s specifications and 150-bp, paired-end sequenced on a HiSeq 4000 to a depth of 90,000 unique molecular identifiers (UMIs) per cell. UMI counts for each cellular barcode were quantified and used to estimate the number of cells successfully captured and sequenced. The Cell Ranger Single Cell Software suite was used for demultiplexing, barcode processing, alignment, and initial clustering of the raw scRNA-Seq profiles.

We used the Chromium instrument and the Single Cell 3′ Reagent kit (v1) to prepare individually barcoded scRNA-Seq libraries following the manufacturer’s protocol (10× Genomics). For QC and to quantify the library concentration we used both the BioAnalyzer (Agilent BioAnalyzer High Sensitivity Kit) and qPCR (Kapa Quantification kit for Illumina Libraries). Sequencing with dual indexing was conducted on an Illumina NextSeq machine using the 150-cycle High Output kit. Sample demultiplexing, barcode processing, and single-cell 3′ gene counting was performed with the Cell Ranger Single Cell Software Suite CR2.0.1. Each droplet partition’s contents were tagged with a UMI — a barcode encoded as the second read of each sequenced read pair. We followed the 10× Single Cell 3′ Reagent Kit’s v2 protocol as written, using 12 cycles for cDNA amplification and 12 cycles for sample index PCR. Samples were sequenced to a depth of approximately 400 million reads per sample on a NovaSeq 6000 (R1 = 26 bp, R(i) = 8 bp, R2 = 91 bp).

Plasmid construction. We synthesized open reading frames (ORFs) of each scGOF-Seq candidate with a Tag-BFP and an 18-nt unique barcode (Supplemental Table 2) (Qinglan Biotech) and cloned them into the pENTR2b vector using KpnI and EcoRV (AFF3, CUX2, FOXO1, GAS7, TSHZ2, and ZFN385D), BamHI and NotI (BACH2, BNC2, EBF1, RARB, RFX7, and TCF4), SalI and NotI (MYT1L and NFATC3), or KasI and NotI (ZRANB3). For modified Perturb-Seq, we synthesized the DNA fragment containing the hU6 promoter and gRNA scaffold followed by CMV promoter and Streptococcus pyogenes Cas9 ORF with a Tag-BFP. We also inserted NotI and HindIII enzyme sites between Tag-BFP and the bGH poly(A) signal (Genewiz). This synthesized DNA fragment was cloned into the pENTR2b vector using KpnI and EcoRV (pENTR2B-Cas9). gRNA was cloned into the pENTR2B-Cas9 vector using BsmbI (pENTR2B-gRNA-Cas9), and corresponding guide barcode was sequentially cloned into the pENTR2B-gRNA-Cas9 vector using NotI and HindIII (Supplemental Table 3).

Adenovirus generation. Recombinant adenoviruses for scGOF-Seq were generated using the pAd/CMV/V5-DEST Gateway recombination system (Life Technologies) after cloning the full-length cDNA into the pENTR vector. Individual adenoviruses were packaged and amplified in HEK-293A cells, and then pools of 3 P1 adenoviruses as detailed below were expanded into high-titer virus. Pool 1: ZNF385D, RARB, and GAS7; Pool 2: EBF1, FOXO1, and TCF4; Pool 3: BACH2, TSHZ2, and NFATC3; Pool 4: ZRANB3, BNYC2, and MYT1L; and Pool 5: AFF3, RFX7, and CUX2. Adenoviruses were purified by PD-10 column (GE Healthcare, 17085101). Titers were determined by plaque assay (PFU). Each virus pool was transduced into HCT116 cells and expression was analyzed by qPCR. FOXO1, TCF4, NFACT3, RFX7, and AFF3 were amplified individually from P1 adenovirus. Adenovirus for modified Perturb-Seq were generated using the pAd/PL-DEST Gateway recombination system (Life Technologies). Individual adenoviruses were packaged and amplified in HEK-293A cells, and then pools of 4 P1 adenoviruses encoding gRNA to target the same TF were expanded into high-titer virus.

Adenovirus transduction. Two hundred to 300 human ND islets were handpicked for each transduction condition and placed in 5 mL round-bottom polystyrene test tubes. Thereafter, islets were washed and incubated with 100 μL serum-free islet culture medium containing 1 mM EGTA and transduced at an MOI of 20. After a 5- to 6-hour transduction, 1 mL of complete islet culture medium with 5% PIM(ABS) Human AB Serum (Prodo Labs, PIM-ABS001GMP) was added overnight. Islets were then transferred to 60-mm nontreated tissue culture dishes (Thermo Fisher Scientific, FB0875713A), and the medium was replaced with fresh islet culture medium every 2 days for 7 days for scGOF-Seq or Perturb-Seq experiments, and 3 days for single-cell intracellular calcium microfluorimetry. The 2 ND donors for scGOF-Seq were procured on different dates and adenoviral transduction resulted in different sets of transcription factors being transduced. For this reason, we analyzed the 2 data sets separately in Figure 4 and Supplemental Figure 9 for ND5 and ND6 (Supplemental Table 4), respectively, lest batch effects yield a biased analysis.

Single-cell intracellular calcium microfluorimetry. Similarly sized human islets from ND donors were handpicked and transduced with adenovirus expressing each candidate cDNA. The day after transduction, islets were partially dispersed using 0.05% trypsin for 5 minutes at 37°C, and then plated on 35-mm glass-bottom dishes with 10-mm microwells (In vitro Scientific, D35-10-0-N) precoated with fibronectin (Sigma-Aldrich, F1141). Cells were washed with islet media and allowed to rest for 2 additional days. On the third day, each plate was washed with KRBH buffer (10 mM HEPES pH 7.4, 140 mM NaCl, 1.5 mM CaCl₂, 3.6 mM KCl, 0.5 mM NaH₂PO₄, 0.5 mM MgSO₄, 2 mM NaHCO₃, and 0.1 % BSA) and incubated with 2.8 mM glucose–containing KRBH buffer for 30 minutes, and then loaded in the dark with 5 μM Rhod-2, AM (Thermo Fisher Scientific, R1244) in KRBH buffer. Cells were washed and transferred into a perifusion chamber placed in the light path of a Zeiss Axiovert fluorescence microscope, and perifused with low glucose (2.8 mM), high glucose (16.8 mM), or KCl (40 mM) in KRBH buffer. β Cells were excited by a Lambda DG-4 150 watt xenon light source (Sutter), using alternating wavelengths of 340 and 380 nm at 0.5-second intervals, and imaged at 510 nm. For each data set, regions of interest corresponding to the locations of 10 to 20 individual cells were selected and images were recorded using an AxioCam camera controlled by Stallion SB.4.1.0 PC software (Intelligent Imaging Innovations). Single-cell intracellular Ca²⁺ mobilization data are plotted as a function of time.

Immunohistochemical and morphometric analyses. Paraffin sections of pancreas far from the margin of pancreatectomy were collected from our previous research (57). In brief, all cases with partial pancreatectomy performed in Ruijin Hospital between 2013 to 2019 were assessed. Cases with a malignant tumor were excluded and then 9 cases of T2D patients and 9 age- and BMI-matched ND subjects were finally included in this study. Detailed information and clinical characteristics for each patient are listed in Supplemental Table 6.

We fixed and processed tissue for immunohistochemistry as previously described (57). We performed immunofluorescence assays using 5-μm-thick paraffin sections obtained and processed as described previously (57). All slides were treated by tissue antigen recovery to improve the fluorescent immune detection of various proteins. Slides were incubated at 4°C overnight with a cocktail of primary antibodies diluted in the blocking reagent for chromogenic and fluorescent immunohistochemical assays. Primary antibodies were prepared at the following dilutions: guinea pig anti-insulin (1:800, DAKO, A056401–2), mouse anti-glucagon (1:400, Abcam, ab10988), and rabbit anti-BACH2 (1:100, Sigma, SAB2108650). Detection was performed using secondary antibodies conjugated to Alexa Fluor 488, 594, and 647 (Jackson ImmunoResearch or Life Technologies). Nuclei were stained with DAPI (Vectashield, Vector Labs) as a marker. Images were captured with an Olympus Microscope or Zeiss LSC 710 confocal microscope. We processed the quantification in a blinded fashion using the cell counter of ImageJ software (NIH) to analyze individual cells located throughout the whole cross section. This analysis scores numbers of positive cells for each marker and calculates the number of cells marked by different signals. To process cell counting, we examined insulin-positive cells (271 ± 63 vs. 207 ± 34 in ND vs. T2D; P = not significant) and glucagon-positive cells (124 ± 31 vs. 205 ± 39 in (ND vs. T2D; P = not significant) per group.

Quantitative analyses of single-cell gene expression. For scRNA-Seq profiles of primary donors (GSE98887), raw sequencing reads of single cells were aligned to hg19 reference genome by Bowtie2-2.2.6 (58). Aligned reads were sorted and indexed by samtools-1.2 (59). Count matrices were measured with R packages GenomicFeatures (1.24.5; ref. 60), GenomicAlignments (1.8.4; ref. 60), and TxDb.Hsapiens.UCSC.hg19.knownGene (3.2.2). Cells with fewer than 100,000 reads were filtered out as low-quality cells. For Perturb-Seq (GSE137766) and scGOF-Seq (GSE136887) profiles, UMI matrices generation, including quality filtering, were performed using the standard 10× Genomics Cell Ranger pipeline (v2.1.1) with hg38 reference genome. Barcode matrices were generated using the standard scPLATE-Seq data processing pipeline, which is available at https://github.com/califano-lab/scPLATE-Seq Specifically, alignment and mapping were performed by STAR (2.5.2a) with flag --outFilterMultimapNmax 1 to exclude multimapping reads for accurate barcode quantification (61). QC reports of cells quality filtering, including distributions of transcriptome complexity (number of genes) and sequencing depth (number of reads/UMI) are included in Supplemental Figures 14–16.

Regulatory networks and transcriptional regulator activity analysis. Islet-specific regulatory networks were reverse engineered by ARACNe (20) on an individual patient basis. ARACNe was run with 200 bootstrap iterations using 1,813 transcription factors (genes annotated in Gene Ontology molecular function database as GO:0003700, “transcription factor activity” or as GO:0003677, “DNA binding” and GO:0030528, “transcription regulator activity” or as GO:00034677 and GO: 0045449, “regulation of transcription”); 969 transcriptional cofactors (a manually curated list, not overlapping with the transcription factor list, built upon genes annotated as GO:0003712, “transcription cofactor activity” or GO:0030528 or GO:0045449); and 3,370 signaling pathway related genes (annotated in GO Biological Process database as GO:0007165 “signal transduction” and in GO cellular component database as GO:0005622, “intracellular”, or GO:0005886, “plasma membrane”). Parameters were set to 0 data processing inequality (DPI) tolerance and mutual information (MI) P value (using MI computed by permuting the original data set as null model) threshold of 1 × 10^–8. Protein activity profiles were then generated by metaVIPER by integrating across all 11 donor-specific regulatory networks (18).

Dimension reduction and clustering analysis. Dimension reduction was done using both gene expression and metaVIPER-inferred activity profiles. For gene expression, cells were first projected to principal component space using Principal Component Analysis (PCA), further projected to 2D t-SNE space. The t-SNE function in the CRAN R package t-SNE-0.1.3 was used for t-SNE analysis. Specifically, 1 – r (cell-wise Pearson correlation in principal component space) was used as distance matrix. For metaVIPER-inferred activity, cells were projected to 2D t-SNE space. The t-SNE function in CRAN R package t-SNE-0.1.3 was used for t-SNE analysis. Specifically, cell-wise activity dissimilarity was used as distance matrix (17).

Clustering analysis was done using the iterClust iterative clustering analysis framework (23). The iterative clustering analysis was performed using the iterClust function in Bioconductor R package iterClust-1.0.2, at metaVIPER-inferred activity level. In total, 3 iterations were done, separating cell populations at different metabolic stress states and cell types sequentially. Using scRNA-Seq data from cells with more than 0.1 million mapped reads, we first projected single cells onto on 2D space with t-SNE.

ELISA for HMOX1. HMOX1 protein induction as a measure of inhibition of the BACH1/2 cellular activity was quantified in cell lines expressing different combinations of these transcription factors. HMOX1 protein levels were determined using the DuoSet IC Human Total HO-1/HMOX-1 kit (R&D Systems, DYC3776). Cells were cultured according to the manufacturer’s protocol and seeded (5,000 cells/well) in 96-well poly-D-lysine–coated plates. Two days after seeding, cells were treated with different concentrations of compound 8 in culture media for 6 hours and then lysed. Lysates were stored at –80°C overnight. ELISA plates were coated overnight with the capture antibody and the next day lysates were tested in the HMOX1 ELISA assay following the manufacturer’s protocol. Optical density of plates was read at 450 nm on a Cytation microplate reader (BioTek).

In vivo mouse studies with BACH inhibitors. Five-week-old male C57BLKS/J db/db mice were purchased from The Jackson Laboratory. All mice were fed Lab Diet 5008, maintained under approved Animal Care and Use protocols for Lilly Research Laboratories, and housed under a 12-hour light/12-hour dark cycle. At 6 weeks of age, compounds were administered orally (10 mL/kg) daily for 15 days in 1% hydroxyethylcellulose, 0.25% Tween 80, and 0.05% antifoam. Nonfasted glucose and insulin were assessed throughout the study. On day 13, after an overnight fast, oral glucose-tolerance tests were performed with a bolus of glucose (0.3 gm/kg).

Glucose-stimulated insulin secretion assays. Islets were preincubated in KRBH buffer for 1 hour at 2.8 mM glucose, followed by incubation in KRBH at 2.8 or 16.7 mM glucose for 1 hour at 37°C. At the end, we collected islets by centrifugation and assayed the supernatant for insulin secretion or lysates for insulin content by ELISA (Mercodia, 10-1247-01). The insulin levels were normalized to protein concentration.

Statistics. Figure 3A and Supplemental Figure 3, C and D quantified the alteration of cell population proportions as opposed to the corresponding control groups, thus P values from Fisher’s exact test (FET) were used to evaluate the significance level. For example, to test whether cluster A1 in Figure 3A is enriched for ND or T2D cells, we constructed a 2 × 2 table (ND%, T2D%) × (A1 cells, all cells). Such table was then used for FET, and the yielded P value of 0.64 indicted that cluster A1 is not significantly enriched for ND or T2D cells. The same FET was performed on other clusters in Figure 3A, as well as cell populations presented in Figure 4, C and E, and Supplemental Figures 7–9. For Figure 8, ANOVA with Dunnett’s method was performed between the 2 groups. A P value of less than 0.05 was considered significant.

Study approval. Animal studies were performed under approved Animal Care and Use protocols by Lilly Research Laboratories. Human studies were approved by the Institutional Review Board of the Ruijin Hospital affiliated to Shanghai Jiao Tong University School of Medicine and were in accordance with the principles of the Declaration of Helsinki II.

Data availability. scRNA-Seq data using the Fluidigm C1 system for all donors have been deposited in the NCBI’s Gene Expression Omnibus (GEO) database (GEO GSE98887). scGOF-Seq and Perturb-Seq data using the 10× Genomics Chromium platform have been deposited in the GEO database under accession numbers GSE136887 and GSE137766, respectively.

J Son and HD contributed equally to designing and performing experiments, to analyzing data, and writing the manuscript. TBF, AME, JLG, SKS, and ZL performed BACH2 inhibitor experiments. J Sun and QW performed immunostaining experiments. AC and DA designed experiments, oversaw research, and wrote the manuscript.