Research Paper|Volume 15, Issue 2|pp 441—458

Clearance of p16Ink4a-positive cells in a mouse transgenic model does not change β-cell mass and has limited effects on their proliferative capacity

Nadine Bahour¹, Lucia Bleichmar¹, Cristian Abarca¹, Emeline Wilmann¹, Stephanie Sanjines¹, Cristina Aguayo-Mazzucato¹

¹Joslin Diabetes Center, Harvard Medical School, Boston, MA 02215, USA

Received: August 11, 2022Accepted: December 27, 2022Published: January 12, 2023

https://doi.org/10.18632/aging.204483

Copyright: © 2023 Bahour et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

Type 2 diabetes is partly characterized by decreased β-cell mass and function which have been linked to cellular senescence. Despite a low basal proliferative rate of adult β-cells, they can respond to growth stimuli, but this proliferative capacity decreases with age and correlates with increased expression of senescence effector, p16Ink4a. We hypothesized that selective deletion of p16Ink4a-positive cells would enhance the proliferative capacity of the remaining β-cells due to the elimination of the local senescence-associated secretory phenotype (SASP). We aimed to investigate the effects of p16Ink4a-positive cell removal on the mass and proliferative capacity of remaining β-cells using INK-ATTAC mice as a transgenic model of senolysis. Clearance of p16Ink4a positive subpopulation was tested in mice of different ages, males and females, and with two different insulin resistance models: high-fat diet (HFD) and insulin receptor antagonist (S961). Clearance of p16Ink4a-positive cells did not affect the overall β-cell mass. β-cell proliferative capacity negatively correlated with cellular senescence load and clearance of p16Ink4a positive cells in 1-year-old HFD mice improved β-cell function and increased proliferative capacity in a subset of animals. Single-cell sequencing revealed that the targeted p16Ink4a subpopulation of β-cells is non-proliferative and non-SASP producing whereas additional senescent subpopulations remained contributing to continued local SASP secretion. In conclusion, deletion of p16Ink4a cells did not negatively impact beta-cell mass and blood glucose under basal and HFD conditions and proliferation was restored in a subset of HFD mice opening further therapeutic targets in the treatment of diabetes.

Introduction

Type 2 diabetes develops in response to over-nutrition and lack of physical activity in subjects with underlying genetic and acquired predisposition to insulin resistance and β-cell dysfunction. Over time, β-cell compensation for insulin resistance fails, resulting in a progressive decline in function [1–6] and mass [7]. The endocrine pancreas is a slow turnover tissue [8] and β-cell proliferation decreases with age [9]. This age-dependent decline in proliferation is partly due to an increase in senescence marker and mediator p16^Ink4a, a cyclin-dependent kinase inhibitor [10, 11]; mice with an additional copy of p16^Ink4a had a significant decrease in β-cell proliferative capacity [12–15].

Aging and senescence are related but not interchangeable terms. Cellular senescence is a stress response that occurs throughout the lifespan in which cells remain metabolically active with an altered phenotype. Senescent cells accumulate with aging, resulting in changes in structure and function that include irreversible growth arrest, resistance to apoptosis and alterations in gene expression [16]. Senescence also leads to the secretion of an array of cell-specific proteins known as the senescence-associated secretory phenotype (SASP) that can induce dysfunction and entry into senescence of surrounding, healthy cells [17].

Previously, we showed that with age and insulin resistance, β-cell senescence increased [18, 19] and senolysis (selective deletion of senescent cells) improved β-cell function, gene identity, and blood glucose levels [19]. Additionally, it has been shown that senolysis can preserve β-cells in a model of Type 1 Diabetes [20]. INK-ATTAC mice are a transgenic whole-body FLAG-tagged line that allows specific deletion of cells expressing p16^Ink4a upon administration of B/B homodimerizer, a synthetic drug that causes dimerization and activation of the caspase 8 only in p16^Ink4a-positive cells [21, 22]. However, whether the elimination of p16^Ink4a-expressing beta cells negatively impacts β-cell mass is unknown and whether residual cell proliferation can be rescued by removing local SASP remains to be determined.

To this end, we set out to explore the effects of removing p16^Ink4a+ senescent cells on the proliferative capacity and mass of β-cells using INK-ATTAC mice as a transgenic model. We hypothesized that removal of this cell population would decrease overall β-cell mass and rescue the proliferative capacity of the remaining cells due to local SASP elimination.

Results

Removal of p16^Ink4a-expressing cells in non-metabolically challenged middle-aged mice had no effect on β-cell mass while proliferation inversely correlated with senescence load

INK ATTAC mice are a mouse transgenic model where p16^Ink4a cells can be specifically tracked through FLAG staining and removed upon B/B homodimerizer injection (Figure 1A) [21]. Validation of this deletion model in pancreatic islets has been previously published [19] and showed that senescent islet cells were deleted after administration of two 3-day courses of B/B homodimerizer, 14 days apart, in 8–9 months old mice. FLAG immunostaining was used as a senescence marker and treatment with B/B homodimerizer significantly increased the proportion of FLAG-negative (non-senescent) islets [19]. Herein, islets from B/B homodimerizer treated animals were isolated and compared to islets from non-treated animals. A significant decrease of transgene markers eGFP and Caspase8 were detected as well as a decrease in p16^Ink4a mRNA (Figure 1B). It is interesting to note that p16^Ink4a transcript levels decrease by 50% after B/B homodimerizer treatment whereas those of eGfp and Caspase 8 are almost undetectable. This can be due to lower expression of second and third cistrons secondary to internal in-frame AUGs. This phenomenon has been described both in Drosophila [23] and mammalian cells [24]. In this model, FLAG intensity is an indicator of cellular senescence because it is driven by the p16^Ink4a promoter. Therefore, confocal pictures of FLAG-stained islets were taken using the same settings such that differences in intensity reflected differences in protein quantity. Quantification of FLAG intensity in insulin positive islets, showed a significant decrease in pancreas from B/B homodimerizer treated animals when compared to those treated with vehicle (Figure 1C and 1D) consistent with decreased senescence. Additionally, increased apoptotic nuclei were counted after TUNEL staining of pancreatic sections from animals treated with B/B homodimerizer (Figure 1E and Supplementary Figure 1). Since, INK-ATTAC mice are whole body transgenics and to understand the expression of p16^Ink4a in other islet cell-types, analysis of scRNASeq data from mouse islets was used and showed that β-cells are the main islet-cell population expressing Cdkn2a, the gene encoding p16^Ink4a (Figure 1F) both in control and insulin resistance conditions (Supplementary Figure 2). However, other islet cell types can also express the transgene to a lower extent. These results confirm that INK-ATTAC animals are a valid model for studying the effect of removing p16^Ink4a expressing β-cells.

Figure 1.

INK ATTAC mice as a model to remove p16Ink4a⁺ cells. (A) The INK ATTAC transgene is driven by the p16^Ink4a promoter and encodes a Caspase 8 moiety that upon B/B homodimerizer administration leads to specific apoptosis of this cell subpopulation and is compared to animals treated with vehicle (referred to as controls). FLAG tag and eGFP expression can be tracked to measure effectiveness of deletion and are a surrogate marker of senescence load. Diagram modified from [21]. (B) qPCR from islets of animals treated with B/B homodimerizer and compared to control, untreated animals. A significant decrease of the transgene transcripts Casp8 and eGFP and p16^Ink4a was observed after treatment with B/B. (C) Representative confocal picture of islets from B/B treated and untreated animals showing a significant decrease of FLAG staining. (D) Quantification of FLAG intensity using image analysis software (Image J) showed a significant decrease in senescence load in islets from animals treated with B/B homodimerizer; n_control = 5 animals, 290 islets analyzed; n_treated = 8 animals, 518 islets analyzed mean+/−SEM, ^*p = 0.0004 by unpaired t-test. (E) Increased number of apoptotic nuclei per field in pancreas from 8-month-old INK-ATTAC mice; n_control = 5 animals, 201 images analyzed; n_treated = 6 animals, 223 images analyzed; mean+/−SEM, ^*p = 0.02 by unpaired t-test. (F) Heatmap of islet scRNASeq expression of Cdkn2a (encodes p16^Ink4) in different islet cell types reveals enriched expression of Cdkn2a in β-cells. No significant changes of mean of (G) blood glucose levels (mg/dL), (H) body weight (g), (I) pancreatic weight (mg) of 8/9-month-old INK-ATTAC mice, n = 6 per group; and (J) Beta cell mass (mg) of 6-month-old INK-ATTAC mice, n_control = 7 and n_treated = 8.

We used male and female 6–9-month-old INK-ATTAC mice fed a chow diet to evaluate the effects of B/B homodimerizer in animals without a metabolic challenge. Analysis revealed no changes in circulating blood glucose (Figure 1G), body weight (Figure 1H), pancreatic weight (Figure 1I), or β-cell mass (Figure 1J). Transcriptional analysis of insulin-sensitive peripheral tissues (fat, liver, red and white muscle) showed no significant differences in transcription of senescence genes p16^Ink4a (encoded by Cdkn2a) (Supplementary Figure 3A) and p21^Cip1 (encoded by Cdkn1a) (Supplementary Figure 3B).

To quantify the proliferation in the same cohort of middle-aged mice, pancreatic sections were stained for BrdU and insulin. B/B homodimerizer did not have an effect on β-cell proliferation at this age (Figure 2A). Consistent with the literature [13], there was a proportion of mice with undetectable proliferation levels, 26% at 6 months and 33% at 8–9 months (Figure 2B). To further assess remnant proliferation and its relation to cellular senescence, the percentage of animals that maintained proliferative capacity and its correlation to FLAG staining were quantified. A significant direct correlation was found between the proportion of non-senescent FLAG negative islets (negative for FLAG staining) and the percentage of proliferation (positive nuclear staining of BrDU) in 6–9-month-old animals (Figure 2C) supporting the idea that higher loads of senescent cells decrease β-cell basal proliferation while maintaining glucose homeostasis.

Figure 2.

Inverse correlation between cellular senescence load and proliferative capacity. (A) Overall proliferation quantification of 6–9-month-old INK-ATTAC animals, n_control = 14, and n_treated = 13; (B) Percentage of islets with proliferating cells in 6 and 8/9-month-old INK-ATTAC mice measured by the percentage of BrDU⁺ Ins⁺ cells, n_6month = 15, and n_8–9month = 12; (C) Correlation between percentage of non-senescent FLAG⁻ islets (negative for FLAG staining) and proliferation (positive nuclear staining of BrDU) in 6–9-month-old (n = 19) INK-ATTAC mice. Line of best fit is shown along with dotted lines indicating 95% confidence intervals. P-value was calculated using the null hypothesis that the slope of the best-fit line equals 0.

Improved β-cell function with removal of p16^Ink4a expressing cells in a HFD model in 1-year-old animals

In adult animals, β-cell proliferation is induced by increased metabolic demand. Therefore, the effects of removing p16^Ink4a cells in a high-fat diet (HFD) model were evaluated. One-year-old INK-ATTAC animals were maintained on a specific diet for 8 weeks in the following groups: control on a chow diet, HFD only, and HFD with four courses of B/B homodimerizer. Weekly body weight was monitored and showed a significant increase of 19% in both the HFD and HFD with B/B homodimerizer groups (Figure 3A). As previously reported [19], HFD accelerated β-cell senescence, which in this model was indicated as an increased Caspase8 (Figure 3B) and p16^Ink4a (Figure 3C) transcript in islets from animals fed a HFD. Representative confocal images of the pancreas from B/B treated animals showed a decrease in FLAG staining (Figure 3D) when compared with those treated with vehicle. HFD diet increased blood glucose by 27%, these levels were restored in the B/B treated group and comparable to chow animals (Figure 3E and 3F). Insulin tolerance test showed no significant differences among the groups (Supplementary Figure 4) even though results were different between males (Supplementary Figure 4A) and females (Supplementary Figure 4B). Fasting hyperinsulinemia was noted in the HFD group and restored after BB treatment (Figure 3G). Similarly, glucose stimulated insulin secretion (GSIS) in vivo during IPGTT, showed increased basal insulin in the HFD treated group which returned to normal after treatment with BB homodimerizer (Figure 3H). These results show that clearance of p16^Ink4a in a HFD model partially restored glucose homeostasis.

Figure 3.

Specific deletion of p16^Ink4a-expressing cells improved metabolic profile. (A) Bodyweight changes compared to the beginning of the HFD. (B) HFD increased transcription of Caspase 8 by qPCR indicating acceleration of senescence in pancreatic islets. FC (fold change) to C (control). (C) HFD increased transcription of p16Ink4a by qPCR indicating acceleration of senescence in pancreatic islets. (D) Representative confocal images of HFD and HFD+ B/B homodimerizer treated islets. FC (fold change) to C (control). (E) Blood glucose percentage changes compared to the beginning of the treatment and (F) area under the curve for the three groups throughout the treatment period. (G) Fasting insulin (ng/mL) levels collected from tail blood during IPGTT (p = 0.08 by two-tailed unpaired t-test). (H) Glucose stimulated insulin secretion (GSIS) evaluated in vivo during IPGTT revealed restoration of basal insulin secretion after treatment with B/B homodimerizer. 1-year old INK-ATTAC mice; n = 8 control group, n = 8 HFD, n = 9 HFD = BB homodimerizer; males and females. ^*p < 0.05 after t-test.

Removal of p16^Ink4a expressing β-cells induced proliferation in a subset of HFD-treated adult mice

A significant correlation between the proliferation percentage and percentage of non-senescent β-cells was also seen at 12-months of age (Figure 4A) indicating that even at this age, removal of p16^Ink4a+ cells can have an impact on beta-cell proliferation.

Figure 4.

Senolysis restored proliferative capacity in a subset of animals after a HFD metabolic challenge in 1-year old mice. (A) Correlation between percentage of non-senescent FLAG⁻ islets and beta-cell proliferation in 12-month-old animals. Line of best fit is shown along with dotted lines indicating 95% confidence intervals. P-value was calculated using the null hypothesis that the slope of the best-fit line equals 0. (B) Pancreatic weight of the groups; (C) beta-cell mass calculated in (mg). (D) Beta-cell proliferation calculated by %BrDU⁺ Ins⁺ cells. Mean+/− SEM; significance calculated by ordinary one-way ANOVA with Tukey’s multiple comparisons. (E) Comparison between responders and non-responders to BB/homodimerizer treatment in the HFD group. AUC for IPGTT before treatment. (F) Blood glucose percentage changes compared to the beginning of the treatment and AUC (G) between responders and non-responders.

There were no significant differences in pancreatic weight in different treatment groups (Figure 4B). Systematic imaging of pancreatic sections revealed no significant changes in β-cell mass amongst groups (Figure 4C) and overall, clearance of p16^Ink4a expressing cells did not change the rate of β-cell proliferation. However, there was a subset of animals treated with B/B homodimerizer where BrdU incorporation was significantly increased and we termed this group “responders” which represented 44% (Figure 4D) of the total treated population. In the remainder 56% of the animals, proliferation did not increase in response to p16^Ink4a clearance (Figure 4D). Further analysis comparing the non-responders to the responders revealed that the latter had better glucose clearance in the initial IPGTT (Figure 4E) and maintained overall lower blood glucose levels than non-responders throughout the experiment (Figure 4F and 4G). However, no differences were found between the two groups in the following parameters: gender, initial ITT, weight gain, fed or fasting insulin and fasting blood glucose (Supplementary Figure 5). These results suggest that responders had better β-cell function at the beginning of the experiment but further experiments are required to fully understand these differences.

β-cell proliferation in a pharmacological model of insulin resistance

Insulin receptor antagonist S961 was used as an additional model to induce acute insulin resistance in 1.5-year-old male and female mice. S961 induced hyperglycemia (Figure 5A) and hyperinsulinemia (Figure 5B). S961 increased β-cell senescence as seen by an increase in Caspase8 transcript (Figure 5C). Proliferation was significantly induced by S961 even in older animals (Figure 5D, Supplementary Figure 6A) however this proliferative rate was not further increased by deletion if p16^Ink4a cells (Figure 5D, Supplementary Figure 6B) and no changes in beta cell mass were found after treatment (Figure 5E). Removal of p16^Ink4a expressing cells in this group worsened blood glucose levels (Figure 5A) indicating that the functional reserve was adversely affected by the removal of this cell population.

Figure 5.

Senolysis in an acute insulin resistance model, S961, did not affect beta-cell proliferative capacity. (A) Area under the curve of blood glucose levels throughout the treatment period, n_control = 8, n_S961 = 3, n_S961+BB = 5. (B) Fed plasma insulin (ng/mL) levels in control, S961, and S961 + BB; significance calculated by repeated unpaired t-test. (C) qPCR of islets from S961 treated animals increased Caspase 8 transcription consistent with accelerated senescence. (D) Proliferation of beta cells (%) in 8–19-month-old mice. (E) Beta cell mass (mg) 18–19 month-old INK ATTAC mice, male and female. Mean +/− SEM, significance calculated by ordinary one-way ANOVA with Tukey’s multiple comparisons.

The targeted p16^Ink4a subpopulation is non-proliferative and non-SASP producing

To further understand the lack of effects of the removal of p16^Ink4a cells in β-cell proliferation, scRNASeq analysis was performed in previously reported β-cells pooled from control and S961 treated animals [25]. t-SNE analysis revealed β-cell subpopulations that were divided into 6 clusters (Figure 6A). Analysis was concentrated in 4 of these subpopulations (Figure 6B). Two senescent subpopulations: p16Ink4a⁺, p21Cip1⁺, and two non-senescent subpopulations: p16Ink4a⁻/p21Cip1⁻ A and p16Ink4a⁻/p21Cip1⁻ B. In the INK-ATTAC model, only the p16Ink4a⁺ cells were removed (Figure 1B) while p21Cip1⁺ remained unchanged (Figure 6K). scRNASeq allowed further analysis of the proliferative capacity and SASP transcription of these specific cell populations. As previously described, p16Ink4a⁺ cells were non-proliferating (Figure 6C and 6D) whereas p21Cip1⁺ and p16Ink4a⁻/p21Cip1⁻ cells upregulated of proliferative genes (Figure 6C, 6D). We believe that the upregulation of proliferation genes in p21Cip1⁺ represents a population of growth arrest during early senescence [26, 27]. This is further supported by significant upregulation of Cdk inhibitors (Figure 6G and 6H) in this subpopulation. Other cyclins and Id were significantly upregulated in the p16Ink4a⁻/p21Cip1⁻ B subpopulation (Figure 6E and 6F). scRNASeq data was further inquired for the transcriptome of SASP factors known to be upregulated in β-cells (Figure 6I and 6J). These changes were confirmed by performing qPCR on islets from BB homodimerizer treated animals. The main SASP-producing β-cell subpopulation was the p21Cip1⁺ which persisted in the islets from mice treated with BB (Figure 6K).

Figure 6.

p16^Ink4a-expressing cells do not proliferate or secrete SASP compared to other senescent cells. (A) tSNE of beta-cell clusters based on Ins2 expression and four subpopulations highlighted: p16^Ink4a+, p21^Cip1+ (senescent) and p16^Ink4a−/p21^Cip1− A and p16^Ink4a−/p21^Cip1− B (non-senescent). (B) Heat map showing beta cell senescence genes. (C, D) Heat map and violin plot showing beta-cell proliferation genes. (E, F). Heat map and violin plot showing beta-cell expression of cyclins and Id. (G, H) Heat map and violin plot of Cdk-inhibitors; (I, J) Heat map and violin plot of SASP genes. (K) qPCR of islets treated with B/B homodimerizer show no changes in the transcript of p21Cip1 and several transcription factors in the INK ATTAC model. Significance calculated by ordinary one-way ANOVA with Tukey’s multiple comparisons.

These results highlight the existence of different senescence subpopulations of β-cells with varying proliferative and SASP producing characteristics which are differentially targeted in this transgenic model. In this particular case, the removal of p16^Ink4a+ cells did not result in the removal of SASP-producing subpopulations which might account for a lack of effects upon proliferative capacity.

In summary, a transgenic model that targets p16^Ink4+ senescent cells, did not alter β-cell mass and increased proliferation levels only in a subset of mice in response to HFD. This is potentially due to remaining senescent cells which are negative for p16^Ink4a and whose local SASP production continues to impair the proliferative capacity of the remaining cells. Additionally, a direct correlation between β-cell proliferative capacity and non-senescent islet cells was found.

Discussion

Cellular senescence has been studied in the context of type 2 diabetes and shown to play a role in the progression of the disease. Deletion of senescent cells has been associated with improved glucose levels and β-cell function [19]. In this study, the effects of deleting p16^Ink4a-expressing senescent cells on β-cell function, mass, and proliferative capacity were evaluated. Removal of this senescent subpopulation did not affect β-cell mass and increased proliferation only in a subset of animals. However, endogenous proliferative capacity is inversely correlated with β-cell senescence underscoring the relationship between them. scRNASeq revealed that the targeted p16^Ink4a β-cell population was non-proliferative and not SASP producing implying that local SASP production remained and was responsible for the terminal cell cycle arrest of β-cells.

Given the known decrease in β-cell mass at the time of T2D diagnosis, pursuing interventions that eliminate a percentage of the remaining cells is concerning. Therefore, a crucial pending question when considering senolysis as a therapeutic option is whether the β-cell mass will be decreased and therefore worsen the development of the disease. Given that senolysis targets senescent cells, our previous reports have found that middle-aged mice (3–6 m) have less than 5% senescent β-cells whereas this increases to approximately 12% at 2-years of age as measured by β-Gal activity [18]. When considering insulin resistance models, S961-treated middle-aged mice had 10% β-Gal positive cells [19]. Put together, these results indicate that senolysis would target a range between 5–10% of the β-cell mass. The impact this would have on regulating blood glucose levels can be estimated by reports of the β-cell functional reserve which has been estimated between 20–25% in rats [28] and 50–70% in type 1 [29, 30] and type 2 [7] diabetes. In this study, blood glucose levels were maintained or improved after removal of p16^Ink4a expressing cells in HFD and non-challenged models. However, in S961 treated animals, senolysis significantly increased hyperglycemia suggesting that in the setting of extreme insulin resistance, targeting senescent cells is counterproductive. The fact that there was a different proliferative response between the HFD and the S961 group when both were treated with B/B homodimerizer is intriguing. This is probably due to β-cell proliferation induced by S961 with no further proliferation with senolysis, however HFD by itself did not produce this increase in proliferation and was more amenable to the positive effects of senolytic interventions. Herein, we show that senolysis specifically targeted at p16^Ink4a+ positive cells does not affect the β-cell mass and does not worsen blood glucose levels under physiological conditions.

Lack of proliferative capacity and SASP secretion are two of the hallmarks of senescent cells [31]. Therefore, decreasing the load of senescent cells in a given tissue could hypothetically increase the proliferative capacity of the remaining cells due to the elimination of SASP. Our results show that p16^Ink4a-directed senolysis did not change basal proliferative rates and was able to induce proliferation in response to HFD only in a subset of animals. These results are surprising given previous reports of decreased islet proliferation in p16^Ink4a overexpression mouse models [11]. These contrasting outcomes might be due to cell-autonomous mechanisms in the deletion model that were induced by p16^Ink4a cells and persisted in the remaining neighbouring cells. This concept is further supported by transgenic p16^Ink4a overexpression demonstrating β-cell-autonomous effect of proliferative restriction [32].

Both p21^Cip1 and p16^Ink4a are known markers and effectors of senescence. Analysis of scRNAseq data revealed the presence of different β-cell subpopulations with different levels of p21^Cip and p16^Ink4a. Interestingly, amongst the two, it was the p21^Cip+ subpopulation the one that transcribed most of the known β-cell SASP factors implying that these were not removed in the INK ATTAC transgenic model. These residual p21^Cip+ senescent and SASP producing cells could continue to impair the proliferative capacity of remaining cells and should be targeted before a claim of lack of effects of senolysis upon β-cell proliferation can be conclusively made. Whereas p21^Cip1-null mice reported normal basal β-cell proliferation [33], it would be interesting to subject these animals to a metabolic challenge (HFD, pregnancy) and see whether proliferation is affected in those instances. Interestingly, an increased β-cell replication has been reported in multiple endocrine neoplasia type 1 syndrome associated with loss of p18 and p27 [34]. Additionally, it has been reported that p27 loss is associated with increased β-cell proliferation and mass after HFD or leptin receptor loss [35, 36].

Whereas senolytic therapies are an approach to senescent cell elimination, they can have adverse effects that need to be considered. First, they are not cell-type specific, which can be a challenge for studying their effects on specific tissues, organs and diseases. Second, senolytic therapies can interfere with wound healing and have oncogenic effects [37]. Whereas translation of senolysis would employ senolytic drugs, the INK-ATTAC transgenic model provides a useful and targeted strategy to initially study the potential effects of removing a bona fide senescent subpopulation. An additional transgenic model that can be used to study the effects of senolysis is the p16-3MR [38] However, future studies should employ drugs to obtain a physiological perspective of its effects in vivo.

Further studies will also need to elucidate the effect of p21^Cip1 on β-cell proliferation as well as inhibiting SASP. Senomorphic drugs, which specifically inhibit SASP secretion, would be one option to minimize the effect of senescent cells on healthy, surrounding cells or senolysis specifically directed at this senescent subpopulation.

Methods

Animals

All experiments were conducted at Joslin Diabetes Center with approval of its Animal Care and Use Committee. Mice were kept in a conventional facility in a 12-hour light/dark cycle with water and food ad libitum and a temperature between 22.2–22.7°C. When specified, a high-fat diet (HFD) 60kcal% fat (Research Diets; NJ, USA) was used for the specified amount of time. C57Bl6/J mice (Jackson Laboratory; ME, USA) were used for some experiments as specified. Breeding pairs of INK-ATTAC mice were a gift from Dr. Jan van Deursen [21] and all the animals used came from our colony. Middle-aged animals were ages 6–9 months with 15 animals 6-month-old (n_female = 7 and n_male = 8) and 12 animals 8–9-month-old (n_female = 6 and n_male = 6). 25 one-year-old animals were used for HFD (n_female = 13 and n_male = 12). 9 animals were used for S961 ages 18–19 months (n_female = 2 and n_male = 7). Animals of different ages were used because age is one of the main factors determining both beta-cell proliferation and senescence load. By analyzing animals of different ages, it is possible to evaluate whether the effects of p16^Ink4a removal are dependent on chronological age.

Assessment of glucose homeostasis

Bodyweight and morning fed glucose levels were monitored weekly. Blood glucose values were measured using a glucometer (Contour; NJ, USA) on blood from tail snip. For intraperitoneal glucose tolerance tests, blood samples for glucose levels from mice fasted for 6 hours were collected at 0, 15, 30, 60, 90, and 120 min after intraperitoneal injection of 10% glucose solution (Sigma Aldrich; MO, USA; 0.02 mL/g body weight). Insulin was measured from serum collected time of sacrifice using an insulin ELISA kit (Mercodia; NC, USA). For insulin tolerance tests, mice were fasted for 4 hours, insulin (Humulin R, Eli Lilly; IN, USA; 1 unit/g body weight) was injected intraperitoneally, and blood glucose was measured at 0, 15, 30, and 60 min.

Senolytic treatment

The deletion protocol of p16^Ink4a-expressing cells for INK-ATTAC mice consisted of the administration of 3-day courses of B/B homodimerizer (Takara Bio; CA, USA; 10 mg/kg) repeated every 14 days to activate the caspase-8 moiety. Mice were treated with an intraperitoneal injection of vehicle (Ethanol: polyethylene glycol 400: Tween 2% at 2:5:43) or B/B homodimerizer (in ethanol: polyethylene glycol 400: Tween 2% at 2:5:43).

S961 treatment

S961 was a generous gift from Dr. Lauge Schaffer (Novo Nordisk; Denmark) [39]. Vehicle (PBS) or 20 nmol S961 was loaded into an Alzet osmotic pump and surgically implanted subcutaneously in the back of the anesthetized mice [40] and changed weekly for a total of two weeks.

Pancreas isolation for quantification of mass and proliferation

Mice were injected with 10 mg/mL BrdU (Sigma Aldrich; MO, USA; 10 uL/g body weight) intraperitoneally, 6 hours before sacrifice. Under anesthesia, the pancreas was excised, weighed, and fixed in 4% (para)-formaldehyde (PFA) for 2 hours and embedded in paraffin for sectioning and immunostaining.

Immunostaining and morphometric evaluation

Paraffin sections were deparaffinized with xylene and ethanol gradients, washed with PBS, permeabilized with Triton-X 0.3%, antigen retrieval was completed with heated citric acid for all stains and blocked with normal donkey serum. After washing with PBS and 2% lamb serum, slides were incubated overnight at 4°C with the primary antibody (Table 1). This was followed by one wash with PBS + 2% lamb serum and incubations for 1hr with secondary antibodies (Table 1). The slides were mounted with Fluoroshield + DAPI for nuclear localization (Sigma Aldrich; MO, USA).

Table 1. Antibodies.

Antibody	Species	Manufacturer and item no.	Concentration
Anti-Insulin	Guinea Pig	Abcam; ab195956	1:400
Anti-FLAG	Rabbit	R&D; MAB8529	1:250
Anti-FLAG	Mouse	Sigma Aldrich; F1804	1:500
Anti-BrdU	Mouse	Sigma Aldrich; B8434	1:50
488 anti-Guinea Pig IgG	Donkey	Jackson ImmunoResearch Laboratories; 706-545-148	1:200
594 anti-Rabbit IgG	Donkey	Jackson ImmunoResearch Laboratories; 711-585-152	1:200

For quantification, islet images were captured systematically covering the whole section in confocal mode on a Zeiss LSM 710 microscope. For mass quantification, the entire section was pictured using the tile-scan system and images were quantified using ImageJ (https://imagej.net/ij/index.html) and Adobe Photoshop (https://www.adobe.com/products/photoshop.html). Beta-cell mass was calculated by multiplying the relative area of beta cells by the pancreatic weight. For all other stains, islets were pictured, coded, and read blindly. For BrdU, 589–3295 cells were counted from at least 10 islets per animal. Quantification of FLAG staining was done using ImageJ and selecting the islet area through insulin staining.

TUNEL stain and quantification

TUNEL staining was done using TUNEL Assay Kit-HRP-DAB (ab206386) from Abcam following their protocols. Positive staining, indicating apoptotic nuclei, was determined by a dark brown color (Supplementary Figure 1). Quantification was done using a brightfield microscope and Ocular^® Scientific Image Acquisition Software. For quantification, systematic alternating pictures were taken at 20× magnification of the whole tissue and positive nuclei were manually counted.

Quantitative real-time PCR

RNA was extracted from cells using the RNEasy Plus Mini Kit (QIAGEN; Germany); SuperScript reverse transcriptase (Invitrogen; MA, USA) was used to reverse transcribe RNA and generate cDNA for quantitative PCRs (Table 2). To measure gene expression levels, we used Fast SYBR green (ThermoFisher; MA, USA) and ΔCT values to βActin were calculated.

Table 2. Primer sequences.

Gene	Forward	Reverse
βActin	ACCGTGAAAAGATGACCCAG	GTACGACCAGAGGCATACAG
Cdkn1a	GCAGATCCACAGCGATATCC	CAACTGCTCACTGTCCACGG
Cdkn2a	CCCAACGCCCCGAACT	GCAGAAGAGCTGCTACGTGAA

Single-cell RNA-seq

Single-cell RNA-seq data presented in this article is an analysis of pooled beta-cells from islets from C57BL6/J mice treated with PBS and S961. Beta cells were identified as having high Ins2 expression and raw data was previously published [25] and deposited under accession number GSE149984. Briefly, islets from animals with S961 or PBS pumps (described above) were isolated from mice for scRNA-seq, cultured overnight, and dispersed. Transcriptomic analysis was performed using the 10× Genomics Chromium Single Cell Gene Expression Assay core at Brigham and Women's Hospital. The Illumina NextSeq500 was used for sequencing, and the two libraries were pooled evenly on one lane. Data analysis was performed by the Bioinformatics and Biostatistics Core at Joslin Diabetes Center. Raw sequencing data were demultiplexed, aligned to the mouse genome and UMIcollapsed using CellRanger [41]. The inclusion criteria were: UMI >500, detected genes >1000, and mitochondrial genes < 20%. Data were analyzed using R. Deconvolution of size factors from cell pools, estimation of technical noise, and denoised Principal Component Analysis (PCA) were done using scran [42, 43]. t-Distributed Stochastic Neighbor Embedding (t-SNE) plots were made using scatter [44]. Cells were clustered into putative subpopulations using a shared-nearest-neighbor graph constructed from the PCA coordinates, and the clusters were found using a spin-glass algorithm [45, 46]. High Ins2 expression was used to identify beta cell clusters. Differential gene expression was assessed using linear modeling with limma [47].

Quantification and statistical analysis

Data are shown as mean ± SEM. For statistical analysis, unpaired Student’s t-tests were used to compare two groups and one-way ANOVA followed by post hoc test for more than two groups. A normality test was performed and when not passed non-parametric statistics were run. A p-value ≤ 0.05 was considered significant. Prism 9.0 software by GraphPad (https://www.graphpad.com/scientific-software/prism/) was used for graphs and statistical analysis (significance and distribution). Animals were assigned to either control, intervention, or treatment groups to have equal age and gender distribution among all groups. The B/B homodimerizer treatment was completely composed of animals positive for the INK ATTAC transgene while the intervention and control groups had a mixture of positive and negative animals. No differences have been shown between animals with and without the transgene. Animals were excluded from the analysis if they became sick or developed physical anomalies.

Data availability

scRNASeq data reported in this article can be obtained with the accession number GSE149984.

Supplementary Materials

Supplementary Figures

Author Contributions

EW, LB, CA, SS and NB performed experiments and analyzed data. NB and CAM wrote the manuscript. CAM designed the study.

Acknowledgments

The authors thank Hui Pan and Jonathan Dreyfuss from Joslin’s Bioinformatic core for assisting with data analysis. To Susan Bonner-Weir for insightful discussion and teaching β-cell mass determination. We would like to thank Christopher Cahill of the Joslin Microscopy Core.

Conflicts of Interest

The authors declare no conflicts of interest related to this study.

Ethical Statement

All experiments were conducted at Joslin Diabetes Center with approval of its Animal Care and Use Committee.

Funding

This study was funded by Institutional Startup Funds to CAM, National Institutes of Health grants 1R01DK132535 to CAM, P30 DK036836 to Joslin Diabetes Center (Cores), Thomas J Beatson Jr foundation grant 2020-010 and the Richard and Susan Smith Family Foundation Award to CAM.

References

1. Weir GC, Bonner-Weir S. Islet β cell mass in diabetes and how it relates to function, birth, and death. Ann N Y Acad Sci. 2013; 1281:92–105. https://doi.org/10.1111/nyas.12031 [PubMed]
2. Weir GC, Bonner-Weir S. Five stages of evolving beta-cell dysfunction during progression to diabetes. Diabetes. 2004 (Suppl 3); 53:S16–21. https://doi.org/10.2337/diabetes.53.suppl_3.s16 [PubMed]
3. Porte D Jr. Clinical importance of insulin secretion and its interaction with insulin resistance in the treatment of type 2 diabetes mellitus and its complications. Diabetes Metab Res Rev. 2001; 17:181–8. https://doi.org/10.1002/1520-7560(200105/06)17:3<181::aid-dmrr197>3.0.co;2-1 [PubMed]
4. Prentki M, Joly E, El-Assaad W, Roduit R. Malonyl-CoA signaling, lipid partitioning, and glucolipotoxicity: role in beta-cell adaptation and failure in the etiology of diabetes. Diabetes. 2002 (Suppl 3); 51:S405–13. https://doi.org/10.2337/diabetes.51.2007.s405 [PubMed]
5. Poitout V, Robertson RP. Minireview: Secondary beta-cell failure in type 2 diabetes--a convergence of glucotoxicity and lipotoxicity. Endocrinology. 2002; 143:339–42. https://doi.org/10.1210/endo.143.2.8623 [PubMed]
6. Leahy JL. Pathogenesis of type 2 diabetes mellitus. Arch Med Res. 2005; 36:197–209. https://doi.org/10.1016/j.arcmed.2005.01.003 [PubMed]
7. Butler AE, Janson J, Bonner-Weir S, Ritzel R, Rizza RA, Butler PC. Beta-cell deficit and increased beta-cell apoptosis in humans with type 2 diabetes. Diabetes. 2003; 52:102–10. https://doi.org/10.2337/diabetes.52.1.102 [PubMed]
8. Finegood DT, Scaglia L, Bonner-Weir S. Dynamics of beta-cell mass in the growing rat pancreas. Estimation with a simple mathematical model. Diabetes. 1995; 44:249–56. https://doi.org/10.2337/diab.44.3.249 [PubMed]
9. Scaglia L, Cahill CJ, Finegood DT, Bonner-Weir S. Apoptosis participates in the remodeling of the endocrine pancreas in the neonatal rat. Endocrinology. 1997; 138:1736–41. https://doi.org/10.1210/endo.138.4.5069 [PubMed]
10. Krishnamurthy J, Torrice C, Ramsey MR, Kovalev GI, Al-Regaiey K, Su L, Sharpless NE. Ink4a/Arf expression is a biomarker of aging. J Clin Invest. 2004; 114:1299–307. https://doi.org/10.1172/JCI22475 [PubMed]
11. Krishnamurthy J, Ramsey MR, Ligon KL, Torrice C, Koh A, Bonner-Weir S, Sharpless NE. p16INK4a induces an age-dependent decline in islet regenerative potential. Nature. 2006; 443:453–7. https://doi.org/10.1038/nature05092 [PubMed]
12. Dhawan S, Tschen SI, Bhushan A. Bmi-1 regulates the Ink4a/Arf locus to control pancreatic beta-cell proliferation. Genes Dev. 2009; 23:906–11. https://doi.org/10.1101/gad.1742609 [PubMed]
13. Tschen SI, Dhawan S, Gurlo T, Bhushan A. Age-dependent decline in beta-cell proliferation restricts the capacity of beta-cell regeneration in mice. Diabetes. 2009; 58:1312–20. https://doi.org/10.2337/db08-1651 [PubMed]
14. Köhler CU, Olewinski M, Tannapfel A, Schmidt WE, Fritsch H, Meier JJ. Cell cycle control of β-cell replication in the prenatal and postnatal human pancreas. Am J Physiol Endocrinol Metab. 2011; 300:E221–30. https://doi.org/10.1152/ajpendo.00496.2010 [PubMed]
15. Chen H, Gu X, Su IH, Bottino R, Contreras JL, Tarakhovsky A, Kim SK. Polycomb protein Ezh2 regulates pancreatic beta-cell Ink4a/Arf expression and regeneration in diabetes mellitus. Genes Dev. 2009; 23:975–85. https://doi.org/10.1101/gad.1742509 [PubMed]
16. Campisi J, d'Adda di Fagagna F. Cellular senescence: when bad things happen to good cells. Nat Rev Mol Cell Biol. 2007; 8:729–40. https://doi.org/10.1038/nrm2233 [PubMed]
17. Coppé JP, Desprez PY, Krtolica A, Campisi J. The senescence-associated secretory phenotype: the dark side of tumor suppression. Annu Rev Pathol. 2010; 5:99–118. https://doi.org/10.1146/annurev-pathol-121808-102144 [PubMed]
18. Aguayo-Mazzucato C, van Haaren M, Mruk M, Lee TB Jr, Crawford C, Hollister-Lock J, Sullivan BA, Johnson JW, Ebrahimi A, Dreyfuss JM, Van Deursen J, Weir GC, Bonner-Weir S. β Cell Aging Markers Have Heterogeneous Distribution and Are Induced by Insulin Resistance. Cell Metab. 2017; 25:898–910.e5. https://doi.org/10.1016/j.cmet.2017.03.015 [PubMed]
19. Aguayo-Mazzucato C, Andle J, Lee TB Jr, Midha A, Talemal L, Chipashvili V, Hollister-Lock J, van Deursen J, Weir G, Bonner-Weir S. Acceleration of β Cell Aging Determines Diabetes and Senolysis Improves Disease Outcomes. Cell Metab. 2019; 30:129–42.e4. https://doi.org/10.1016/j.cmet.2019.05.006 [PubMed]
20. Thompson PJ, Shah A, Ntranos V, Van Gool F, Atkinson M, Bhushan A. Targeted Elimination of Senescent Beta Cells Prevents Type 1 Diabetes. Cell Metab. 2019; 29:1045–60.e10. https://doi.org/10.1016/j.cmet.2019.01.021 [PubMed]
21. Baker DJ, Wijshake T, Tchkonia T, LeBrasseur NK, Childs BG, van de Sluis B, Kirkland JL, van Deursen JM. Clearance of p16Ink4a-positive senescent cells delays ageing-associated disorders. Nature. 2011; 479:232–6. https://doi.org/10.1038/nature10600 [PubMed]
22. Baker DJ, Childs BG, Durik M, Wijers ME, Sieben CJ, Zhong J, Saltness RA, Jeganathan KB, Verzosa GC, Pezeshki A, Khazaie K, Miller JD, van Deursen JM. Naturally occurring p16(Ink4a)-positive cells shorten healthy lifespan. Nature. 2016; 530:184–9. https://doi.org/10.1038/nature16932 [PubMed]
23. Wall AA, Phillips AM, Kelly LE. Effective translation of the second cistron in two Drosophila dicistronic transcripts is determined by the absence of in-frame AUG codons in the first cistron. J Biol Chem. 2005; 280:27670–8. https://doi.org/10.1074/jbc.M500255200 [PubMed]
24. Al-Allaf FA, Abduljaleel Z, Athar M, Taher MM, Khan W, Mehmet H, Colakogullari M, Apostolidou S, Bigger B, Waddington S, Coutelle C, Themis M, Al-Ahdal MN, et al. Modifying inter-cistronic sequence significantly enhances IRES dependent second gene expression in bicistronic vector: Construction of optimised cassette for gene therapy of familial hypercholesterolemia. Noncoding RNA Res. 2018; 4:1–14. https://doi.org/10.1016/j.ncrna.2018.11.005 [PubMed]
25. Midha A, Pan H, Abarca C, Andle J, Carapeto P, Bonner-Weir S, Aguayo-Mazzucato C. Unique Human and Mouse β-Cell Senescence-Associated Secretory Phenotype (SASP) Reveal Conserved Signaling Pathways and Heterogeneous Factors. Diabetes. 2021; 70:1098–116. https://doi.org/10.2337/db20-0553 [PubMed]
26. Ebrahimi AG, Hollister-Lock J, Sullivan BA, Tsuchida R, Bonner-Weir S, Weir GC. Beta cell identity changes with mild hyperglycemia: Implications for function, growth, and vulnerability. Mol Metab. 2020; 35:100959. https://doi.org/10.1016/j.molmet.2020.02.002 [PubMed]
27. De Cecco M, Ito T, Petrashen AP, Elias AE, Skvir NJ, Criscione SW, Caligiana A, Brocculi G, Adney EM, Boeke JD, Le O, Beauséjour C, Ambati J, et al. L1 drives IFN in senescent cells and promotes age-associated inflammation. Nature. 2019; 566:73–8. https://doi.org/10.1038/s41586-018-0784-9 [PubMed]
28. Olsson R, Carlsson PO. A low-oxygenated subpopulation of pancreatic islets constitutes a functional reserve of endocrine cells. Diabetes. 2011; 60:2068–75. https://doi.org/10.2337/db09-0877 [PubMed]
29. Damond N, Engler S, Zanotelli VRT, Schapiro D, Wasserfall CH, Kusmartseva I, Nick HS, Thorel F, Herrera PL, Atkinson MA, Bodenmiller B. A Map of Human Type 1 Diabetes Progression by Imaging Mass Cytometry. Cell Metab. 2019; 29:755–68.e5. https://doi.org/10.1016/j.cmet.2018.11.014 [PubMed]
30. Wang YJ, Traum D, Schug J, Gao L, Liu C, Atkinson MA, Powers AC, Feldman MD, Naji A, Chang KM, Kaestner KH, and HPAP Consortium. Multiplexed In Situ Imaging Mass Cytometry Analysis of the Human Endocrine Pancreas and Immune System in Type 1 Diabetes. Cell Metab. 2019; 29:769–83.e4. https://doi.org/10.1016/j.cmet.2019.01.003 [PubMed]
31. Hernandez-Segura A, Nehme J, Demaria M. Hallmarks of Cellular Senescence. Trends Cell Biol. 2018; 28:436–53. https://doi.org/10.1016/j.tcb.2018.02.001 [PubMed]
32. Helman A, Klochendler A, Azazmeh N, Gabai Y, Horwitz E, Anzi S, Swisa A, Condiotti R, Granit RZ, Nevo Y, Fixler Y, Shreibman D, Zamir A, et al. p16(Ink4a)-induced senescence of pancreatic beta cells enhances insulin secretion. Nat Med. 2016; 22:412–20. https://doi.org/10.1038/nm.4054 [PubMed]
33. Cozar-Castellano I, Haught M, Stewart AF. The cell cycle inhibitory protein p21cip is not essential for maintaining beta-cell cycle arrest or beta-cell function in vivo. Diabetes. 2006; 55:3271–8. https://doi.org/10.2337/db06-0627 [PubMed]
34. Karnik SK, Hughes CM, Gu X, Rozenblatt-Rosen O, McLean GW, Xiong Y, Meyerson M, Kim SK. Menin regulates pancreatic islet growth by promoting histone methylation and expression of genes encoding p27Kip1 and p18INK4c. Proc Natl Acad Sci U S A. 2005; 102:14659–64. https://doi.org/10.1073/pnas.0503484102 [PubMed]
35. Uchida T, Nakamura T, Hashimoto N, Matsuda T, Kotani K, Sakaue H, Kido Y, Hayashi Y, Nakayama KI, White MF, Kasuga M. Deletion of Cdkn1b ameliorates hyperglycemia by maintaining compensatory hyperinsulinemia in diabetic mice. Nat Med. 2005; 11:175–82. https://doi.org/10.1038/nm1187 [PubMed]
36. Franklin DS, Godfrey VL, O'Brien DA, Deng C, Xiong Y. Functional collaboration between different cyclin-dependent kinase inhibitors suppresses tumor growth with distinct tissue specificity. Mol Cell Biol. 2000; 20:6147–58. https://doi.org/10.1128/MCB.20.16.6147-6158.2000 [PubMed]
37. Aguayo-Mazzucato C. Functional changes in beta cells during ageing and senescence. Diabetologia. 2020; 63:2022–9. https://doi.org/10.1007/s00125-020-05185-6 [PubMed]
38. Palmer AK, Xu M, Zhu Y, Pirtskhalava T, Weivoda MM, Hachfeld CM, Prata LG, van Dijk TH, Verkade E, Casaclang-Verzosa G, Johnson KO, Cubro H, Doornebal EJ, et al. Targeting senescent cells alleviates obesity-induced metabolic dysfunction. Aging Cell. 2019; 18:e12950. https://doi.org/10.1111/acel.12950 [PubMed]
39. Schäffer L, Brand CL, Hansen BF, Ribel U, Shaw AC, Slaaby R, Sturis J. A novel high-affinity peptide antagonist to the insulin receptor. Biochem Biophys Res Commun. 2008; 376:380–3. https://doi.org/10.1016/j.bbrc.2008.08.151 [PubMed]
40. Dai C, Kayton NS, Shostak A, Poffenberger G, Cyphert HA, Aramandla R, Thompson C, Papagiannis IG, Emfinger C, Shiota M, Stafford JM, Greiner DL, Herrera PL, et al. Stress-impaired transcription factor expression and insulin secretion in transplanted human islets. J Clin Invest. 2016; 126:1857–70. https://doi.org/10.1172/JCI83657 [PubMed]
41. Zheng GX, Terry JM, Belgrader P, Ryvkin P, Bent ZW, Wilson R, Ziraldo SB, Wheeler TD, McDermott GP, Zhu J, Gregory MT, Shuga J, Montesclaros L, et al. Massively parallel digital transcriptional profiling of single cells. Nat Commun. 2017; 8:14049. https://doi.org/10.1038/ncomms14049 [PubMed]
42. Lun AT, Bach K, Marioni JC. Pooling across cells to normalize single-cell RNA sequencing data with many zero counts. Genome Biol. 2016; 17:75. https://doi.org/10.1186/s13059-016-0947-7 [PubMed]
43. Scialdone A, Natarajan KN, Saraiva LR, Proserpio V, Teichmann SA, Stegle O, Marioni JC, Buettner F. Computational assignment of cell-cycle stage from single-cell transcriptome data. Methods. 2015; 85:54–61. https://doi.org/10.1016/j.ymeth.2015.06.021 [PubMed]
44. McCarthy DJ, Campbell KR, Lun AT, Wills QF. Scater: pre-processing, quality control, normalization and visualization of single-cell RNA-seq data in R. Bioinformatics. 2017; 33:1179–86. https://doi.org/10.1093/bioinformatics/btw777 [PubMed]
45. Reichardt J, Bornholdt S. Statistical mechanics of community detection. Phys Rev E Stat Nonlin Soft Matter Phys. 2006; 74:016110. https://doi.org/10.1103/PhysRevE.74.016110 [PubMed]
46. Xu C, Su Z. Identification of cell types from single-cell transcriptomes using a novel clustering method. Bioinformatics. 2015; 31:1974–80. https://doi.org/10.1093/bioinformatics/btv088 [PubMed]
47. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015; 43:e47. https://doi.org/10.1093/nar/gkv007 [PubMed]