**3. Discussion**

To our knowledge, no studies to date have assessed the heterogeneity of individual islet cell responses to HFD feeding in mice. The use of sc-RNAseq technologies allows us to employ clustering analysis to segregate populations of cells that exhibit similar gene expression patterns and identify the heterogeneity of cellular responses. In this study, we were interested in the early responses of islet endocrine cells to HFD feeding. Prior studies have shown that one week of 60% kcal from fat a HFD results in mild glucose intolerance in male C57BL/6 mice [10,34], suggesting that β-cell failure in the

setting of HFD-induced insulin resistance may occur early during obesity. Our results from sc-RNAseq analysis show that the minor clusters of β cells (β5, β7, β8, β10, and β11) exhibit an impairment in β cell function and survival, whereas some of the major clusters of β cells (β3 and β4) demonstrate transcriptomic changes consistent with a compensatory response; however, others (β1 and β2) still exhibit no statistically significant change. Our findings, therefore, sugges<sup>t</sup> that the failure of a majority of β cells to compensate for prevailing insulin needs might portend the eventual development of hyperglycemia/T2D. Further studies that include analysis of several time points would provide more insight into whether the defective early response found in the minor β cells clusters triggers T2D or if it is just the result of a compensatory response to an obesogenic diet. Similarly, the robustness of a large proportion of α cells, which exhibit gene expression changes that promote proliferation and survival, might permit the more robust production of glucagon that exacerbates dysglycemia. Our studies support the concept that the balance of di fferent populations of cell types in the islet (those more capable of responding to prevailing stress vs. those less capable) might account for the eventual risk for the development of dysglycemia/T2D in the setting of HFDs.

Several limitations should be noted in our study. First, our study was limited to an early period (one week) following obesogenic diet exposure. Therefore, we cannot know if/how gene expression responses change with time. Do the responsive β-cell clusters remain robustly responsive in their gene expression patterns with more prolonged feeding? Do the less responsive clusters become more responsive with time? Are specific clusters more susceptible to apoptosis or dedi fferentiation over time? Does the islet isolation process influence mRNA expression? In future studies, applying our annotation and clustering analyses should allow us to address these questions using timed feeding cohorts. Second, our study was limited to male C57BL/6J mice, which are known to be the most susceptible to the deleterious e ffects of HFD feeding. The use of mice with more robust islet responsiveness, such as the C57BLKs/J strain [35] or female C57BL/6J mice [36], is likely to provide more insight into the nature and number of β-cell clusters that are responsive to the demand for insulin production. Third, the relevance of our study to human obesity and T2D remains unclear because of species di fferences and the added genetic heterogeneity in humans that compounds the inherent heterogeneity of the islet cell types. Fourth, although our data are compared to controls, because the procedure of islet isolation is a stressful process, we cannot rule out that it could lead to more exaggerated changes in gene expression profile that may impact our findings. Finally, we should note that sc-RNAseq may underestimate the nature of gene expression changes since the depth of sequencing by single-cell technologies is considerably lower than that seen in bulk sequencing technologies. Thus, our findings regarding the minimal gene expression changes seen in the majority of β-cell clusters may not fully report actual changes occurring in these cells. Nevertheless, our studies provide new biologic insight into the early changes in gene expression that occurs in islet cell subsets in the early phase of obesogenic diets.

#### **4. Materials and Methods**

## *4.1. Animal Studies*

Male C57BL/6J mice were purchased from Jackson Laboratories at 7 weeks of age. All mice were kept on a standard 12 h–12 h light–dark cycle with ad libitum access to chow and water. Mice were maintained at the Indiana University vivarium according to protocols approved by the Institutional Animal Care and Use Committee. At 8 weeks of age, animals were fed either a LFD (10% kcal from fat; Research Diets D12450B) or a HFD (60% kcal from fat; Research Diets D12492) for one week.

## *4.2. Islet Isolation*

Mouse islets were isolated from collagenase-perfused pancreata, as previously described [37]. After isolation, mouse islets (approximately 180–200 islets) were handpicked and digested with Accutase (EMD Millipore Corporation, Temecula, CA, USA) containing 2 U/mL of DNAse for 5 min at 37 ◦C sob agitation (1000 rpm). Digested cell islets were washed several times with PBS+2%FBS to

eliminate DNAse and then filtered using a cell strainer (40 μm). Single-cell suspensions and samples with more than 90% viability were used for sc-RNAseq.

#### *4.3. Preparation of Single Cell 3 RNA-seq Library*

Single-cell 3 RNA-seq experiments were conducted using the Chromium single-cell system (10 x Genomics, Inc., Pleasanton, CA, USA) and Illumina sequencers at the Center for Medical Genetics (CMG) of the Indiana University School of Medicine. Each cell suspension, a collection of 180–200 islets isolated from an individual mouse, was first inspected under a microscope for cell number, cell viability, and cell size. For the initial cell suspension with low viability, the single-cell preparation was further processed, which included centrifugation, resuspension, and filtration to remove cell debris, dead cells, and cell aggregates. Single-cell capture and library preparation were performed according to the Chromium Single Cell 3 Reagent Kits V3 User Guide (10x Genomics PN-1000075, PN-1000073, PN-120262). An appropriate number of cells were loaded on a multiple-channel microfluidics chip of the Chromium Single Cell Instrument (10x Genomics) with a targeted cell recovery of 8000 to 10,000. Single-cell gel beads in an emulsion containing barcoded oligonucleotides and reverse-transcriptase reagents were generated with the V3 single-cell reagen<sup>t</sup> kit (10x Genomics). Following cell capture and cell lysis, cDNA was synthesized and amplified. Illumina sequencing libraries were then prepared with the amplified cDNA. The resulting libraries were assessed with Agilent TapeStation. The final libraries were sequenced using a custom program on Illumina NovaSeq 6000. About 550 million read pairs were generated for each sample, with 28 bp of cell barcode and unique molecular identifier reads and 91 bp RNA reads.

#### *4.4. Analysis of sc-RNAseq Sequence Data*

CellRanger 3.0.2 (http://support.10xgenomics.com/) was utilized to process the raw sequence data generated. Briefly, CellRanger used bcl2fastq (https://support.illumina.com/) to demultiplex raw base sequence calls generated from the sequencer into sample-specific FASTQ files. The FASTQ files were then aligned to the mouse reference genome mm10 with RNA-seq aligner STAR. The aligned reads were traced back to individual cells, and the gene expression level of individual genes was quantified based on the number of UMIs detected in each cell.

The filtered gene–cell barcode matrices generated from CellRanger were used for further analysis with the R package Seurat (Seurat development version 3.0.0.9200) [38,39] with Rstudio version 1.1.453 and R version 3.5.1. Quality control (QC) of the data was implemented as the first step in our analysis. We filtered out genes that were detected in less than five cells and cells with less than 200 genes. To further exclude low-quality cells in downstream analysis, we used the Outlier function from the R package scatter [40] together with a visual inspection of the distributions of the number of genes, UMIs, and mitochondrial gene content. Cells with an extremely high or low number of detected genes/UMIs were excluded. In addition, cells with a high percentage of mitochondrial reads were also filtered out. After removing likely multiplets and low-quality cells, the gene expression levels for each cell were normalized with the NormalizeData function in Seurat. Highly variable genes were subsequently identified.

To integrate the single-cell data of the HFD and LFD samples, functions FindIntegrationAnchors and IntegrateData from Seurat were applied. The integrated data were scaled, and PCA was performed. Clusters were identified with the Seurat functions FindNeighbors and FindClusters. The FindConservedMarkers function was subsequently used to identify cell cluster marker genes. Cell cluster identities were manually defined with the cluster-specific marker genes or known marker genes. The cell clusters were visualized using the t-Distributed Stochastic Neighbor Embedding (t-SNE) plots and Uniform Manifold Approximation and Projection (UMAP) plots.

To investigate cell cluster/type-specific responses triggered by di fferent diet treatments, a pseudobulk count approach was applied. The counts of each gene of all cells within a cell cluster for each biological replicate were aggregated. In this way, the gene–cell count matrix was

transformed into a gene–sample level count matrix. The R package muscat (version 1.0.1, [41]) was employed for generating pseudobulk counts and differential gene expression analysis between conditions with the same cluster. In more details, the Seurat object generated from the integrative analysis was first converted into a SingleCellExperiment object using the as.SingleCellExperiment function in Seurat. The SingleCellExperiment object was then used as input for muscat. Gene–cell count data were aggregated with the function aggregateData, and differential gene expression analysis of the pseudocounts was performed with edgeR [42]. GO gene set enrichment analysis was performed using the R package clusterProfiler [43]. The datasets generated and analyzed in the present study are available in the GEO repository (accession number: GSE162512).

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2218-1989/10/12/513/s1, Table S1: Percentage of cells per cluster identified by sc-RNAseq, Figure S1: Genes used for identification of cell clusters, Figure S2: Identification of differentially expressed genes of the major β cell clusters.

**Author Contributions:** Conceptualization, A.R.P., W.W., S.A.T., and R.G.M.; methodology, A.R.P. and S.A.T.; formal analysis, H.G. and Y.L.; writing—original draft preparation, A.R.P.; writing—review and editing, A.R.P., S.A.T., R.G.M.; funding acquisition, R.G.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Institutes of Health, gran<sup>t</sup> number R01 DK60581, R01 DK 105588, P30 DK097512, and P30 DK020595.

**Acknowledgments:** The authors wish to acknowledge K. Orr and X. Xuei at Indiana University for technical assistance in the studies.

**Conflicts of Interest:** The authors declare no conflict of interest.
