Mouse and Human Single Pancreatic β Cell Transcriptomics Reveal Sexual Dimorphism of Transcriptomes and Identify Sex-dependent Type 2 Diabetes Altered Genes

Type 2 diabetes, characterized by malfunction of pancreatic β cells, is affected by multiple cues including sex differences. Nevertheless, mechanisms of sex differences in type 2 diabetes susceptibility and pathogenesis remain unclear. Using single-cell RNA sequencing (scRNA-seq) technology, we showed that sexual dimorphism of transcriptome exists in both mouse and human β cells, and differs significantly in human under diabetic condition. Our analysis further revealed the existence of sex-dependent type 2 diabetes altered genes in both mice and human beings, suggesting divergences in pathological mechanisms of type 2 diabetes between sexes. Our results indicated that sex should be taken into consideration when treating diabetes, which was further validated by the sex-matched and sex-mismatched islet transplantation in mice. Compared to sex-matched transplants, sex-mismatched transplants showed downregulation of genes involved in the longevity regulating pathway in β cells and led to impaired glucose tolerance in diabetic mice. Taken together, our findings could advance current understanding of type 2 diabetes pathogenesis with sexually dimorphic perspectives and provide new insights to the development of precision medicine. Context and Significance Amounting evidence has shown that sex plays a crucial role for glucose homeostasis and onset risk of diabetes. Based on pancreatic β cell transcriptome profiling at single cell resolution, we focused on identifying sexually dimorphic transcriptomes and sex-dependent type 2 diabetes genes in both mouse and human. Given the importance revealed by our study that different pathological mechanisms of type 2 diabetes exist between sexes, sex differences should be taken into consideration for innovative precision medicine when treating diabetes. Highlights Sex-biased gene expression pattern exists in both mouse and human β cells, and differs dramatically under diabetic condition in human. In both human and mouse β cells, abundant sex-dependent type 2 diabetes altered genes are identified. Sex-specific interspecies conservation of type 2 diabetes altered metabolic pathways exists. Sex-matched islet transplantation in diabetic mice yields better control of glucose homeostasis than sex-mismatched transplantation. Graphic Abstract


INTRODUCTION
The efficacy of current anti-diabetic medication varies significantly among individuals with diabetes, highlighting the importance of personalized treatment for type 2 diabetes (T2D). However, the bottleneck for precision medicine lies in the heterogeneous nature of the disease, which not only hinges on genetic predispositions that have been identified by GWAS (Genome Wide Association Study), but also other cues, including sex, diet, and aging. Among these factors, sex differences should be first considered for personalized therapy since sex is one of the most recognizable traits.
However, in most current studies on metabolism using rodents, female animals are usually neglected, because male animals have the tendency to show better disease phenotypes and are dominantly used (Zucker and Beery, 2010). And this experimental bias on sex hampered novel and comprehensive acknowledgement of metabolic pathological mechanisms. Thus, the National Institutes of Health (NIH) demands sex differences should be emphasized in preclinical studies (Clayton and Collins, 2014;Mauvais-Jarvis et al., 2017), which should be especially stressed in T2D as a global pandemic.
Glucose homeostasis is controlled by pancreatic islet, which is mainly composed of α cells, β cells, δ cells and PP cells. α cells elevate glucose level by secreting glucagon to promote hepatic glucose synthesis, and β cells release insulin to decrease glucose level by stimulating blood glucose uptake by fat, muscle, liver, and intestine cells, etc. δ cells secret somatostatin to downregulate hormones releasing from both α cells and β cells via paracrine signaling (Zhang et al., 2007). Polypeptide from PP cells, the most infrequent islet cell type, has effects on both gastric and pancreatic secretions (Holzer et al., 2012).
Despite the similar cell type components of islet architecture shared by both male and female, there exist profound sex differences in islet physiological function and metabolism, signaling pathways involved in hormone releasing, and diabetes occurrence (Vital et al., 2006). For example, female rats are more susceptible than males to streptozotocin (STZ) induced diabetes, suggesting that female rat β cells are more sensitive to STZ toxicity than male β cells (Vital et al., 2006). Similarly, maternal high fat diet results in insulin resistance and oxidative stress-induced β cell loss specifically in male offsprings, rather than in female offsprings, indicating that female islets might have self-protective machinery against oxidative stress (Yokomizo et al., 2014). In human, type 2 diabetes occurs more frequently in men with younger age and less BMI than women (Kautzky-Willer et al., 2016;Logue et al., 2011;Wandell and Carlsson, 2014). Insulin sexual dimorphism of DNA methylation was also observed in human islets by whole islet genomewide DNA methylation sequencing. It is suggested that sex differences of methylome are associated with differences of islet genes expression and insulin secretion level (Hall et al., 2014). However, sexual dimorphism of pancreatic islet β cell has not been investigated at the single cell level, and these studies could reveal sex differences of gene expression in islet β cells between male and female and provide better treatment plans for diabetic patients from both sex groups.
Sex differences in diabetes susceptibility, development and progression have been previously reported, suggesting the existence of sex-dependent diabetes associated genes. Previous studies showed that androgen receptor specifically expressed in male islet β cells and plays an important role in regulating glucose-stimulated insulin secretion in both mice and humans (Navarro et al., 2016). It has also been reported that KLF14 allele variants show increased female-specific T2D risk, probably via female-specific fat storage and distribution (Small et al., 2018). Recent banding studies reported that single-cell RNA sequencing technology has been applied to identify novel diabetes altered genes (Lawlor et al., 2017;Segerstolpe et al., 2016;Xin et al., 2016). Nevertheless, the diverse sex-dependent diabetes associated genes and molecular pathways have not been comprehensively investigated from these studies.
Here, we systematically analyzed the single cell gene expression profiles of healthy and diabetic β cells from mice and human. We found a considerable number of genes had sex-biased expression in β cells, and showed that this sex dimorphism of human was significantly altered under diabetic condition in terms of the number of sex-biased genes. Furthermore, we identified abundant sex-dependent diabetes altered genes, suggesting that the molecular mechanisms mediating diabetes pathogenesis in males and females have important differences. Based on the recognition of the sex differences in T2D altered genes and pathways in β cells, we concluded that sex as a biological variance should be emphasized in diabetes treatment. And this conclusion was further supported by the sex-matched and sex-mismatched islet transplantation in mice. Compared to sex-matched transplants, we found that genes involved in the longevity regulating pathway tended to be down-regulated in β cells of sex-mismatched transplants, and glucose tolerance notably decreased in diabetic mice transplanted with sex-mismatched islets. Together, our results not only advanced current understanding of T2D pathogenesis, but also provided new insights and targets for developing sex-dependent precision medicine to treat diabetes.

Identification of Male and Female Pancreatic β Cell Transcriptomes in Mouse and Human
To obtain expression profiles of mouse pancreatic β cells, we employed flow cytometry to isolate single islet cells from dissociated mouse islets with live dye staining ( Figure 1A). Totally, we collected and sequenced 5,472 islet cells (3,264 male; 2,208 female) from both male and female mice, including 1,056 islet cells (1,056 male; 768 female) of 8-week-old healthy mice, 2,208 islet cells (1,152 male; 1,056 female) of 9-month-old healthy and diabetic mice, 672 male transplanted islet cells (9 months post-transplant) in kidney capsules of both male and female recipient mice, and 768 endogenous pancreatic islet cells (384 male; 384 female) of recipient mice (detailed sample information in Table 1). Single-cell RNA sequencing library was constructed by a modified Smart-seq2 protocol. The average library size is 370K reads per cell, and the average number of genes detected in these cells is about 1500 per cell. We retained 4,662 cells that have at least 500 genes detected ( Figure S1A). The retained cells had high total counts of unique molecular identifiers (UMIs) mapped to gene exons (Figure S1B), and low fractions of UMI counts of mitochondrial genes ( Figure S1C), suggesting we obtained a set of qualified scRNA-seq profiles of mouse islet for the main purpose of this study.
From these 4,662 cells, we went on to identify β cells for downstream analysis. Firstly, we applied an adjusted CPM method (adjCPM) to normalize our scRNA-seq data by excluding the union of the top 2 genes with the highest expression from all cells while calculating the normalization factors ( Figure S2A and S2B; see Methods). Secondly, hierarchal clustering of scRNA-seq profiles based on the union of the top 10 highly expressed genes from all cells was used to discriminate the identity of the cells ( Figure S3A). After performing principal component analysis with the identified hyper variable genes (HVGs) and subsequent visualization by t-distributed stochastic neighbor embedding (t-SNE) ( Figure S3B; see Methods), the retained islet cells of healthy (8-week-old; 9-month-old) and diabetic mice were found to be aggregated into two clusters exhibiting differential expression of β cell marker gene Ins2 ( Figure 1B and 1C). Then, 3,912 β cells (2,197 male; 1,715 female) from mice in different conditions were identified through high expression of Ins2 ( Figure 1C and S4A). The sex of mouse β cells was labeled on the t-SNE map (male in blue; female in red; Figure 1D), and further confirmed by X and Y chromosome genes ( Figure  S4B).
To obtain expression profiles of both healthy and diabetic β cells in men and women, we used a previously published human islet scRNA-seq data set (Xin et al., 2016). The single cells in this data set were from 7 healthy male, 5 healthy female and 3 T2D male, 3 T2D female donors, and we reanalyzed the data by adding the factor of sex. The result of t-SNE showed that all the cells were segregated into 4 distinct clusters, including 460 β cells ( Figure 1E and 1F) which were used for further analysis. Similar to the result of mouse β cells ( Figure 1D), the human male and female β cells were also not completely overlapped on the t-SNE plot ( Figure 1G and S4C). To verify the existence of sexually dimorphic gene expression, we applied Kolmogorov-Smirnov test (K-S test) to analyze the differential expression of Ins2 between male and female β cells, and found its expression was significantly different between two sexes both in mouse and human ( Figure 1D and 1G), indicating sexual dimorphism of β cell transcriptomes exist in both species and could be important to β cell functions.

Sex-biased Gene Expressions in Mouse β Cells Under Healthy and T2D Conditions
To identify genes that display sex-biased expression pattern in mouse pancreatic β cells, we first compared the transcriptional profiles of male β cells with female β cells from 8-week-old C57BL/6J mice ( Figure 2A). Differentially expressed genes were sorted out by MAST (Finak et al., 2015), with cutoffs of FDR<=0.05 and |log2(fold change)|>0.2 (the same cutoffs were used for all the differential analysis using MAST in our study). In total, we obtained 162 differentially expressed genes (DEGs), comprising 37 genes expressed higher in male and 125 genes expressed higher in female ( Figure 2A). Among them, only four genes were on the sex chromosomes, including X chromosome (Xist) and Y chromosome (Eif2s3y, Ddx3y, Uty) ( Figure 2B), suggesting that considerable level of sexual dimorphism exist in healthy β cell transcriptomes.
To investigate whether such sexually dimorphic β cell transcriptomes also exist in diabetic conditions, β cells from high-fat-diet (HFD) induced diabetic mice were collected for scRNA-seq. Firstly, we used HFD to induce T2D model with β cell failure in C57BL/6J mice (fed with HFD from 8-week-old) as previously described (Sone and Kagawa, 2005). Consequently, characteristic T2D phenotypes were observed in those diabetic mice of both male and female, including high body weight, increased fasting insulin level, and impaired glucose tolerance. Notably, sex differences of T2D were also observed: the fasting serum insulin level of healthy male mice is significantly higher than healthy female mice ( Figure 2C). In addition, the area under curve (AUC) of glucose tolerance test showed that the impairment of glucose tolerance in female T2D mice was more significant than that in male T2D mice, as p value of female T2D <0.001 and male T2D<0.01 ( Figure 2D). Then, to explore the genetic basis for such differences, we obtained single β cell transcriptomes from 9-month-old diabetic mice (HFD feeding up to 7 months) and age-matched healthy mice (normal diet, ND) of both sexes. The comparison of single β cell transcriptional profiles of male and female from either healthy or T2D mice was carried out by MAST with the same cutoffs as described above. As shown in the results, 394 DEGs between male and female in 9-month-old healthy mice were identified, including 200 genes expressed higher in male and 194 genes expressed higher in female. In parallel, 81 male highly expressed genes, and 52 female highly expressed genes were identified in β cells from T2D mice ( Figure 2E).
To further elucidate sex-biased pathways based on these comparisons, gene set enrichment analysis (GSEA) was performed with the cutoff set as FDR<=0.25 (the same cutoff was used for all GSEA in our study) (Subramanian et al., 2005). In β cells from both 8-week-old and 9-month-old healthy animals, the longevity regulating pathway was enriched in males, and related genes (Ins1, Ins2, Irs2, Hspa1a, Hspa8) were expressed significantly higher in males ( Figure 2F; S5A and S5B). Intriguingly, the ribosome pathway was consistently enriched in females in β cells from both healthy (8-week-old; 9-month-old) and diabetic mice ( Figure S5A; S5B and S5C). Importantly, in β cells from diabetic mice, both the N-Glycan biosynthesis pathway and Notch signaling pathway were enriched in males. Conversely, the JAK-STAT signaling pathway, ferroptosis pathway, spliceosome pathway, carbohydrate digestion and absorption pathway were enriched in female β cells from diabetic mice ( Figure S5C). Expression patterns of representative sex-biased DEGs involved in the GSEA results were shown in the heatmap ( Figure 2F). Above all, our results validated the existence of sex-biased gene expression in β cells of both healthy and diabetic mice, indicating that the pathological mechanism of T2D might differ between males and females.

Abundant Sex-dependent T2D Altered Genes Were Found in Mouse β Cells
Given the sex-biased gene expressions detected in T2D β cells, we hypothesized T2D development may differ in males and females. Previous studies compared the single-cell transcriptional profiles of T2D islets with that of healthy islets to identify T2D altered genes (Lawlor et al., 2017;Segerstolpe et al., 2016;Xin et al., 2016) without considering the factor of sex. To dissect that in a sex-dependent manner, we compared the single-cell transcriptome of β cells from T2D mice with that from age-matched healthy mice. Firstly, we did a sex-independent comparison, in which β cells from both sexes were used for differential analysis, and we defined 98 sex-independent genes. Then, we compared the diabetic β cells with healthy β cells of the same sex to obtain DEGs in either male or female specific manner, and consequently compared the DEGs of each sex with the 98 sex-independent genes. The non-overlapping parts were defined as sexdependent T2D altered genes, including 57 female-dependent and 66 male-dependent genes ( Figure 3A and 3B), among which the Spc5 was upregulated in female diabetic β cells but down-regulated in male diabetic β cells ( Figure S6A).
To note, among the female-dependent T2D altered genes, the top 5 up-regulated genes were Fxyd2, Hsp8, Cpe, G6pc2 and Serp1, while the top 5 down-regulated genes were Sdf2l1, Ssr4, Shisal2b, Chac1, Ewsr1( Figure S6A). Notably, according to previous reports, the Fxyd2 knock-out mice showed significantly improved glucose tolerance and enlarged islet size, suggesting its negative regulation of islet growth (Arystarkhova et al., 2013). Hsp8/HSC70 encodes the heat-shock cognate protein 70, which was recognized as an endogenous danger signal in β cell to trigger inflammation in diabetes (Alam et al., 2009). Cpe encodes Carboxypeptidase E, which plays a role in insulin synthesis by regulating proinsulin to insulin conversion for insulin synthesis (Davidson and Hutton, 1987). The upregulation of Cpe was in line with the increased fasting serum insulin level observed in the T2D mice ( Figure 2C). The G6pc2 has been reported to be a negative regulator for glucose stimulated insulin secretion in mice, and its deletion results in female specific body fat reduction (Pound et al., 2013), which supports our concept of sex-dependent T2D altered gene.
In parallel, among the male-dependent T2D altered genes, the top 5 genes up-regulated in T2D were Iapp, P4hb, Rps29, Trpm5, and Gabbr2, while the top 5 genes down-regulated in T2D were Cox4i1, Malat1, Ubb, Nisch, Ndufb7 ( Figure S6A). Iapp encodes islet amyloid polypeptide, whose aggregation results in β cell failure and T2D onset (Cooper et al., 1987). P4hb encodes protein disulfide isomerase A1, whose expression is essential for maturation of proinsulin and survival of β cell in high-fat diet fed mice (Jang et al., 2019). Trpm5 encodes transient receptor potential channel type M5, whose deletion can reduce body weight gain and improve glucose tolerance in mice fed with high calorie diet (Larsson et al., 2015). Collectively, these results demonstrated that abundant genes' expression was altered under T2D diabetic condition in a sex-dependent manner.
To further systematically investigate the sex differences in T2D altered pathways, GSEA was performed based on comparison of β cells from T2D and healthy mice of each sex. In females, several pathways enriched in T2D-β cells such as pancreatic secretion pathway and longevity regulating pathway (multiple species) have direct links to β cell function and onset of T2D. Conversely, oxidative phosphorylation pathway and ribosome pathway were enriched in healthy β cells. At the same time, in males, T2D-β cells had ribosome pathway enriched and healthy β cells had type 2 diabetes mellitus pathway enriched ( Figure 3C and 3D). The expression patterns of genes involved in the pathways mentioned above were shown in the heatmap ( Figure 3E) and violin plots ( Figure S6B). Altogether, both the distinctive patterns of sex-dependent T2D altered genes and different T2D associated pathways suggested divergence of T2D pathogenesis between sexes.

Sex Bias of Gene Expression Expanded in Human β Cells Under T2D Conditions
Since we observed sex-biased gene expression in mouse β cells, we further examined whether sex-biased gene expression also exists in human β cells by comparing transcriptomes of male and female β cells from healthy or T2D donors, respectively. In β cells from healthy people, we identified 12 genes with sex-biased expression, of which 10 genes were highly expressed in males, such as CASC4, DDX3Y, FKBP9, RAB21, RIN2, RPL7, RPS26, RPS28, SCD, SMIM8. Except for the Y chromosome gene DDX3Y, the others are located on autosomes. By contrast, only 2 genes were identified as being expressed higher in females, including one X-linked gene XIST and one autosomal gene KNG1 ( Figure 4A). In β cells from T2D patients, we found 143 and 155 genes were highly expressed in males and females, respectively ( Figure 4A and 4B). Notably, compared with β cells from healthy people, the results showed that the sex bias of gene expression in β cells from T2D patients expanded dramatically ( Figure 4B).
Next, in β cells from T2D patients, we applied GSEA to investigate the sex-biased pathways. The significantly enriched pathways in β cells from T2D men included ribosome, neurodegenerative disease, oxidative phosphorylation, protein export and proteasome. Interestingly, pathways enriched in β cells from T2D women, were more directly related to β cell function and T2D, including estrogen signaling, autophagy, insulin signaling, insulin secretion, AMPK signaling, type 2 diabetes mellitus, mTOR signaling, insulin resistance and longevity regulating pathway ( Figure 4C). The expression of representative sexbiased genes involved in the above pathways was shown in the heatmap ( Figure 4D). Collectively, in both male and female T2D patients, the escalated sex bias of gene expression and pathway regulation in β cells indicated that there might also exist sex-dependent alteration of gene expression caused by T2D occurrence in human.

Sex-dependent T2D Altered Genes Were Identified in Human β Cells
We used the same method to investigate sex-dependent T2D altered genes in human β cells as we did with mouse data. Firstly, we identified 282 sex-independent genes by sex-independent comparison between β cell transcriptomes of T2D and healthy donors. Secondly, we obtained DEGs by comparing the β cell transcriptomes of diabetic and healthy donors in males and females, respectively. And the nonoverlapping parts of these DEGs and sex-independent genes were defined as sex-dependent T2D altered genes, including 26 female-dependent, and 142 male-dependent ( Figure 5A and 5B). Intriguingly, one of the sex-dependent T2D altered genes, SCD (encoding an enzyme for oleic acid synthesis), was altered oppositely. It was down-regulated in female diabetic β cells but up-regulated in male diabetic β cells ( Figure S7A).
For the female-dependent T2D altered genes, compared with the β cells of healthy women, the top 5 genes that significantly up-regulated in T2D female β cells included CNN3, CDC42SE1, CTNNB1, LSM12, PPP1CB, while the top 5 genes significantly down-regulated were DNAJC15, SCNM1, HLA-L, PPWD1, PSMG3 ( Figure S7A). Among them, it has been reported that MCJ/DnaJC15 (DnaJ heat shock protein family (Hsp40) member C15) negatively affects respiratory chain, and its depression could inhibit lipid accumulation in mouse liver (Hatle et al., 2013). In parallel, as for the male-dependent T2D altered genes, compared with the β cells of healthy men, the top 5 significantly up-regulated genes in male T2D β cells were MT-TS2, RELL1, PCDH7, CDH10 and FAM229B, while the top 5 significantly down-regulated genes were TM4SF1, MCOLN3, LOC100996447, GPX4 and COA3 ( Figure S7A). According to previous reports, the gene variants of GPX4 (glutathione peroxidase 4) determine the onset of diabetic nephropathy in type 1 diabetes, and lead to higher risk in men (Monteiro et al., 2013), suggesting GPX4 may also play a crucial role for the onset of type 2 diabetes in a male-dependent manner. Moreover, we found a male-dependent T2D up-regulated gene GSTT1 (log2(fold change)=1.340; Table S2-2) may be a potential candidate gene for diabetic precision medicine, as its null genotype combined with GSTM1 null genotype contributes to highly susceptible to troglitazone triggered hepatotoxicity in type 2 diabetes patients of Japan (Watanabe, 2003).
Furthermore, GSEA analysis was performed to explore the sex differences in T2D altered pathways in human β cells. The results showed that several diabetes-related pathways were significantly enriched, among which ribosome, type 1 diabetes, endocrine and other factor-related, calcium reabsorption, and insulin resistance pathways were enriched in the β cells of T2D women, and thermogenesis and oxidative phosphorylation pathways were enriched in the β cells of healthy women. Concomitantly, the ribosome pathway was the most significantly enriched in β cells of T2D men, which was consistent with the result of T2D women β cells. Moreover, pathways including peroxisome, insulin resistance, mTOR signaling pathway, insulin secretion, type 2 diabetes mellitus, longevity regulating pathway-multiple species, fatty acid degradation, glycolysis/ gluconeogenesis, and pancreatic secretion were significantly enriched in β cells of healthy men ( Figure 5C), which were different from the results of women cells, indicating there is sex divergence in the pathogenesis of T2D in human.

Cross-species Comparison of T2D Altered Pathways in Both Sexes
With an interspecific perspective, we also compared sex differences in the pathogenesis of T2D between mouse and human. In both species, significant enrichment of the thermogenesis and oxidative phosphorylation pathways in healthy β cells was found only in females ( Figure 5D and 5E), while the significant enrichment of type 2 diabetes pathway in healthy β cells was found only in males ( Figure 5F), suggesting interspecific conservation of sex dimorphism in the pathogenesis of T2D. Nevertheless, there were also interspecific discrepancies. For example, the ribosome pathway was significantly enriched in healthy β cells of female mice, but it was also enriched in T2D β cells of male mice. Whereas in humans, the pathway was significantly enriched in T2D β cells of both men and women ( Figure 5G).

Mice with Sex-matched β Cell Transplant Exhibited Better Control of Glucose Homeostasis
Altogether, the existence of sexual dimorphism in mouse and human β cell transcriptome informed us that sex as an important factor in β cell function and pathogenesis of T2D should be emphasized when treating diabetes. In order to validate this critical concept, islet transplantations (an experimental treatment for insulin insufficient diabetes mellitus) with sex-matched and sex-mismatched islets were performed in mice. Male islets were isolated from 6-8 weeks old mice and transplanted into kidney capsules of age-matched male and female mice, respectively. Transplanted islets (TX-islet) and endogenous islets (Endo-islet) from recipient mice were all collected and dissociated for scRNA-seq profiling, 9 months post-transplant ( Figure 6A). Firstly, the β cells were identified by high expression of Ins2 as previously described ( Figure 6B), and the transcriptomes of sex-matched and sex-mismatched transplant-β cells were compared with that of male and female endogenous β cells, respectively ( Figure  6C). The correlation analysis revealed that the overall sex-matched transplant-β cell transcriptomes were closer to endogenous islet β cells from both male and female recipients ( Figure 6C). Secondly, transplantβ cell transcriptomes were directly compared between sex-matched and sex-mismatched groups for DEGs. Sex-matched transplant-β cells had 3 genes (Kap, Sfrp5, and Akr1c21) significantly up-regulated and 2 genes (Ovol2, Matn2) significantly down-regulated ( Figure S8A and S8B). Notably, previous study reported that overexpression of Sfrp5 (down-regulated in obesity and T2D) can ameliorate impaired glucose tolerance in mice (Ouchi et al., 2010). In humans, high level of serum SFRP5 was correlated with lower risk of T2D onset (Carstensen-Kirberg et al., 2017). Moreover, the results of GSEA showed that "longevity regulating pathway-multiple species" was significantly enriched in β cells of sex-matched transplant-islets and among the 16 leading-edge genes, Ins2, Sod1, Sod2 and Foxa2 are directly linked to β cell secretion ( Figure 6D). Collectively, these results suggested that sex-matched islet transplantation could be more beneficial to the function of β cells of the transplants, thus better controlling glucose homeostasis.
To further validate the advantage of sex-matched transplantation, we transplanted islets from male and female mice into streptozotocin induced female diabetic mice, respectively ( Figure 6E). Although the hyperglycemia (>350 mg/dl) of the diabetic recipient with sex-matched or sex-mismatched transplantation restored to normal levels 1 month post-transplant (<144 mg/dl; Figure S8C), the glucose tolerance of the diabetic mice with sex-matched transplantation was significantly better as evidenced by the oral glucose tolerance test ( Figure 6F). Above all, our scRNA-seq analysis of β cells and experimental validation concluded that sex should be taken into consideration in diabetes treatment.

Conclusions
According to previously published work, it has been well recognized that sexual dimorphism exists in many organs or systems, such as heart, kidney and immune responses (Klein and Flanagan, 2016;Ransick et al., 2019;Skelly et al., 2018). What's more, a recent study in humans has shown that sexual dimorphism not only exhibits in neuron cells, but also is associated with susceptibility of mental diseases (Lobentanzer et al., 2019). Herein, T2D is a complex metabolic disorder characterized with islet β cell failure and impacted by sex. Sex differences in islet β cell physiological function and diabetes prevalence have been recognized, but still need to be better understood. To decipher the sexually dimorphic T2D pathogenesis, we used scRNA-seq to comprehensively measure the transcriptomes of healthy and diabetic β cells from mice and humans of both sexes. We found a considerable number of genes have significantly sex-biased expressions in β cells in both mouse and human, and further revealed that in human, this sex bias of gene expressions was significantly enlarged in diabetics.
Furthermore, we identified 122 sex-dependent T2D altered genes in mouse and 167 sex-dependent T2D altered genes in human, suggesting important differences in the molecular mechanisms of diabetes pathogenesis between males and females in both mouse and human. Moreover, consistent results of species conservation were uncovered, including 3 sex-specific conserved T2D altered pathways. Specifically, oxidative phosphorylation and thermogenesis pathways were significantly enriched in the healthy female β cells, and this was consistent with the result that expression levels of several genes involved in oxidative phosphorylation cooperatively decreased in the islets of T2D patients (Olsson et al., 2011). Meanwhile, the T2D pathway was significantly enriched in the healthy male β cells. Collectively, these results provided innovative sex-specific targets for future studies on the precision treatment of T2D.
In current clinical trials of islet transplantation, sex matching between donor and recipient has not been stressed. Based on the recognition of sex differences in the transcriptome of β cells, we concluded that sex as a crucial biological variance should be emphasized in the treatment of diabetes. And this conclusion was further supported by the sex-matched and sex-mismatched islet transplantation in mice. Compared with the sex-mismatched group, the β cells of transplants in the sex-matched group showed significant enrichment of the longevity regulating pathway and glucose tolerance notably improved. For better curative effect, our result raised the necessity for sex-matched islet transplantation, and even for sexmatched stem cell-based cell replacement therapy for diabetes treatment.
Beyond the islet β cells focused in this study, sex differences of T2D susceptibility were also associated with sex steroid hormones. It has been found that endogenous estrogens are protective against T2D in females and the risk of T2D increases following menopause (Mauvais-Jarvis et al., 2013;Tramunt et al., 2020). Estrogen not only improves islet β cell function and survival (Tiano and Mauvais-Jarvis, 2012), but also stimulates the secretion of GLP-1 by both islet α cells and intestine L cells to maintain glucose homeostasis (Handgraaf et al., 2018). Along with that, deficiency of the male sex hormone, testosterone, increases the T2D risk in male (Basaria et al., 2006;Mauvais-Jarvis, 2011). Furthermore, a recent report revealed that testosterone improves insulin secretion through androgen receptor on male islet β cells in both mouse and human (Navarro et al., 2016). Taken together, sex differences in both islet β cells and sex steroid hormones need to be highlighted in the development of sex-specific precision medicine.

Limitations of Study
Plenty of sex-biased genes and pathways were identified in β cells of mouse and human either under healthy or type 2 diabetes conditions. Meanwhile, a series of novel sex-dependent T2D altered genes were also uncovered. However, the physiological significance of sex-biased genes and sex-dependent T2D altered genes have not been validated in this stud
See also Table S1.      In 86.7% of the mouse islet cells, the UMI count fraction of the top 2 most expressed genes are more than 0.5. The union of the top 2 most expressed genes from all these cells includes 9 genes (Gcg,Ins1,Ins2,Malat1,Ppy,Pyy,Rn45s,Ttr,, which are excluded while calculating total count for normalization using adjCPM method.

(B) t-SNE map of retained mouse cells colored by expression of sex chromosome linked genes.
The reds for X-chromosome linked gene Xist, and the blues for Y-chromosome linked genes Eif2s3y, Ddx3y, Uty.

Animals and High Fat Diet Induced Diabetic Mice
All experiments were performed in accordance with the University of Health Guide for the Care and Use of Laboratory Animals and were approved by the Biological Research Ethics Committee of Tongji University. We housed all the mice under the specific pathogen free grade environment of animal facility at Tongji University, Shanghai, China. Adult male and female mice (6-8 weeks old) were purchased from Shanghai Slac Laboratory Animal. To induce mouse T2D model, male and female C57BL/6J mice (n=5 per sex) were fed with high-fat diet from 8-week-old to 9-month-old as a high-fat-diet (HFD) group, and equal number of male and female mice were fed up with normal diet (ND) in an age-matched way to a control potential confounding factor, age.

Islet Transplantation
Islets were isolated from 6-8 weeks old male ICR mice (n=20), and the detailed process of islet isolation was as previous description (Zmuda et al., 2011). Similar size islets were handpicked after purification under a stereomicroscope, and 3 age-matched healthy male and 3 female ICR mice were selected as a recipient with about 300-400 islets were transplanted under the kidney capsule. About 9 months after, transplanted islets were dissected and scraped from the kidney capsule, and the pellet was collected and dissociated into a single cell for scRNA-seq. At the same time, the endogenous islets of recipient mice were also isolated and dissociated into a single cell for scRNA-seq. Diabetic mice and islet transplantation were both performed as previous description (Zmuda et al., 2011). Adult male and female C57BL/6J mice were selected as islet donor, and age-matched female C57BL/6J mice were chosen as recipient. C57BL/6J mice were injected with STZ at the dose of 170 mg/kg after 6 hours of fast. Mice that exhibited non-fasting hyperglycemia (>350 mg/kg) with 3 consecutive detection were regarded as diabetic mice for islet transplantation. Each diabetic mouse was transplanted with about 350 islets, and mice that non-fasting blood glucose (<144 mg/kg) recover normal were selected for physiological experiments 1 month post-transplant.

Preparation of Single Islet Cells
Mouse islets were pooled from either 6-8 weeks old male or female C57BL/6J and ICR mice. For islets isolation, 0.5 mg/mL collagenase P (Roche) was poured into pancreas by perfusion of the common bile duct. After digestion, islets were purified through Histopaque (Sigma) gradient centrifugation. The Histopaque buffer was made with mixing Histopaque-10771 and Histopaque-11191 together as the ratio of 5:6. Purified islets were dissociated into single cells as follows: washed the islets with cold PBS at least two times, spun at 1000 rpm for 2 minutes, then replace PBS with 1mL TrypLE Express (Gibco) and incubated at 37 • C for 10-15 minutes with pipetting cell pellet by using P1000 pipette gently. Stop the reaction with low glucose DMEM (containing 10% FBS, 1% HEPES, 1% PenStrep), and centrifuged cells at 4 • C, 1000 rpm for 2 min. Then washed the cell pellet with cold PBS for 1 time, and the resuspended cells in PBS (containing 0.5%BSA) was filtered with a 40-micrometer strainer to get single cell suspension.
To obtain live single cell, Calcein Blue, AM (Invitrogen) was used to measure the viability of cell and high viability single cell was sorted by using BD FACS Aria II flow cytometry. Single islet cell was sorted into 96-well plates containing lysis buffer.

Library Construction and NGS Sequencing
Single-cell RNA-seq libraries were constructed according to Smart-seq2 protocol except that oligod (T) primers comprising 16bp cell barcode sequence and 9bp molecular barcode sequence were used to allow sample pooling and molecular counting. In addition, part of Truseq read2 sequence was used to replace ISPCR sequence in the oligod (T) primer thus to be compatible with Illumina sequence platform. Cells in lysis buffer were denatured, reverse transcribed in the presence of template switching oligos and preamplified by adding both ISPCR and ISPCR-read2 primers with 24 PCR cycles. cDNA of cells with different barcodes were then pooled together and were purified using 0.8x AMPure XP beads. Homemade Tn5 enzyme was used to tagment cDNA. Final amplification was processed using P7-index primers (Truseq) and P5-index primers (Nextera). Sequencing libraries were purified twice using 0.6x AMPure XP beads and once with 1x AMPure XP beads, and were sequenced on Illumina HiSeq X10 platform with default parameters. Primer and adaptor sequences were all listed in supplementary Table S3 (information of barcode).

Physiology Experiment
For intraperitoneal injection glucose tolerance test, normal diet and high-fat diet feeding mice fasted for 16 hours. 1 g/Kg body weight of glucose was intraperitoneally injected, and blood glucose was measured at time point of 0 min, 15 min, 30 min, 60 min, 120 min after glucose injection using glucometer (ACCU-CHEK). Blood samples were collected from the mice that fasted for 6 hours, to measure the insulin level by using Mouse Ultrasensitive Insulin ELISA Kit (ALPOC).
For oral glucose tolerance test, mice were fasted for 6 hours, and glucose was gavaged at a dose of 2 g/Kg. Tail blood was collected for blood glucose detection, at time point of 0 min, 15 min, 30 min, 60 min, 120 min, using glucometer (ACCU-CHEK).

Preprocessing Before Normalization
Reads were stored in paired-end fastq format. Reads of one end of fragment contain cell barcode and UMI information which was subsequently extracted and added to the name of corresponding reads of the other end in fastq file containing molecular sequence. That barcode/UMI information was in fastq-R2 files. Next, the generated single-end fastq files were cleaned by Trim Galore! (Krueger, 2012;Martin, 2011) with the parameter --length 30.
Then, the FastQC (Andrews, 2012) was applied to check reads quality before subsequent alignment on mm9 genome using STAR (Dobin et al., 2013). The parameters of STAR software were --AlignEndsType EndToEnd --out FileterMismatchNoverReadLmax 0.04 --outSAMattrIHstart 0 --outSAMmultNmax 1 --outFilterMultimapNmax 1. GTF file of mm9 reference genome was derived from the RefSeq gene annotation (Pruitt et al., 2005) file downloaded from UCSC genome browser database (Karolchik et al., 2003), http://genome.ucsc.edu/cgi-bin/hgGateway?db=mm9. After alignment, SAM files underwent demultiplexing by Catadapt (Martin, 2011) with --overlap 16 -no-indels --match-read-wildcards parameters; reads were removed if they were assigned to more than one cells. Then, we removed PCR duplicates by UMI-tools (Smith et al., 2017). Next, featureCounts (Liao et al., 2014) was utilized to quantify gene expression levels according to the above mentioned GTF file. Once the expression profiles were generated, cells with fewer than 500 expressed genes (UMI count > 1) were considered "low-quality" and were removed ( Figure S1). The fraction of transcripts from mitochondrial genes in each single cell was investigated as it was used as an indicator for general quality control. Part of our ICR cells had high transcripts fraction of mitochondrial genes, which can be accounted for that kidney cells intrinsically express high-level mitochondrial genes (Park et al., 2018). So we kept those cells for analysis because our transplantation cells were under renal capsule and either kidney cells or affected transplant cells may well be introduced in our data. Next, only genes expressed (UMI count > 1) in 5 or more cells were used for further analysis. Eventually, there are 14,152 genes and 4,662 cells retained in mouse scRNA-seq data.
The sample information of the Illumina high-throughput sequencing data is listed as following table.

Normalization
Here we used an adjusted count per million (CPM) normalization method, named adjCPM.The CPM method assumes that the total molecule reads are equal among cells. Similarly, our adjCPM method has an assumption except excluding a few genes in calculating the total number of molecular reads of each cell. It is based on the observation that sometimes a few top expressed genes can possess more than 50% total UMI count ( Figure S2). These genes were selected as union set of the top 2 expressed genes of single cells in which the sum of the corresponding top 2 genes' UMI counts is beyond 50% of the total count of the cell. Specifically, we obtained 9 genes (Gcg,Ins1,Ins2,Malat1,Ppy,Pyy,Rn45s,Ttr, and they were excluded when calculating the total count of each cell. It is well known that the scRNA-seq experiments are affected by drop-out events especially for UMI method. Here, we also observed that many genes have zero expression value because of insufficient detection power or intrinsically under-expressed.
To address this issue, SAVER  was applied to recover gene expression in the entire matrix.

Cell Type Identification
Since cell type identification is critical for the rest analysis, two steps were carried out to robustly define the cell type. Firstly, we selected the top 10 expressed genes of each cell, resulting in a total of 118 genes to perform the subsequent hieratical clustering analysis. Accordingly, we primarily assigned 4 clusters of single cells after removing 111 cells which had ulterior distance in clustering dendrogram and annotated each cluster based on marker gene expression ( Figure S3). Secondly, by a lowess regression between the log2 mean and log2 coefficient of variance of each gene's normalized UMI count across all single cells, we selected 3,000 genes according to the residue of lowess regression as hypervariable genes (HVGs). Next, principal component analysis (PCA) was applied to the log2 transformed normalized expression of these HVGs using sklearn (Pedregosa et al., 2011) after centering and scaling, and top 25 PCs were selected based on the statistical significance of the fraction of total variance explained by each of them, which was estimated using jackstraw (Chung and Storey, 2015). We used the top 25 PCs as input for subsequent t-Distributed Stochastic Neighbor Embedding (t-SNE) analysis. To have a stable t-SNE visualization, we tried a lot of combinations of parameters (perplexing [10,15,20,25,30], early exaggeration [12,15,20], learning rate [200,300,400,500,600,800] and 50 random seeds) using sklearn (Pedregosa et al., 2011), and found they generally gave very similar results. We finally chose perplexing:15, early exaggeration:12, learning rate:500 to generate the t-SNE plots shown in the figures. Next, Density-based spatial clustering algorithm, DBSCAN (Ester et al., 1996) with parameters eps=5 and MinPts=5 was applied on the resulting two-dimensional t-SNE map to group the single cells into cliques ( Figure S3B). After removing singletons, we examined each single cell clique in order of clique size with following criteria: 1. only cells with consistent cluster/clique labels generated by hierarchical clustering using the top expressed genes and by DBSCAN were retained; 2. differential expression analysis of single cells in clique 3 (C3) compared the single cells in other cliques demonstrated that the single cells in C3 should be excluded because they highly expressed cell cycle related genes ( Figure S3D). Finally, we kept 4,467 high-quality mouse single cells for downstream analysis, and the cell type labels were inferred from the marker genes highly expressed in each single cell clique ( Figure S3C).

Human Data Acquisition and Processing
The raw reads of full-length scRNA-seq data of human pancreatic islets were downloaded from https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE81608 (Xin et al., 2016). As afore mentioned steps for mouse data processing, we firstly generated raw read count matrix using hg19 reference genome and the GTF file was derived from the RefSeq gene annotation file acquired from http://genome.ucsc.edu/cgi-bin/hgGateway?db=hg19. And we calculated Transcripts Per Million (TPM) as normalized value of gene expression.
Then, we performed cell type identification on 1,600 human islet cells in a process similar to that applied to mouse data. Firstly, we used the union of the 10 top expressed genes in each single cell across all cells (in total 102 genes) to perform hierarchical clustering, through which 62 cells were removed. Afterwards, we performed HVG selection, PCA and t-SNE visualization and applied DBSCAN on the resulted t-SNE map with parameters eps=5 and MinPts=10. After removing 13 singleton cells, there were 6 cliques of cells grouped by DBSCAN. The cell types of single cells in 4 of the 6 cliques were labeled according to highly expressed hormone genes (GCG known as a marker gene of α cells, INS known as a marker gene of β cells, SST known as a marker gene of δ cells and PPY known as a marker gene of PP cells), while cells in the other 2 cliques were considered as contaminated ones because of the high expression of non-islet tissue markers like SPARC, COL1A1, COL1A2, COL1A3, FABP5, LGALS1 (Xin et al., 2016). After this step, 1,479 cells with consistent cell-type labels given by hierarchical clustering and by DBSCAN were retained. In addition, due to Xin et al. had defined the cell type for these cells, we added an additional filtering step, that is, cells were filtered if their cell type labels given by us were inconsistent with that by Xin et al. Finally, 1,477 human islet cells were retained for downstream analysis.

Differential Expression Analysis and Gene Set Enrichment Analysis
MAST (Finak et al., 2015) was used to perform differential expression analysis between male and female or between T2D and healthy scRNA-seq profiles of β cells. Cellular detection rate (CDR, which means the number of genes detected in each single cell) was controlled as a covariate while estimating treatment effects. Significantly regulated genes were identified by using |log2(fold change)| >0.2 and false discovery rate <=0.05 as cutoffs. For identification of functional pathways enriched in the detected differentially expressed genes between different sexes or between T2D and healthy conditions in our scRNA-seq data, we performed gene set enrichment analysis (GSEA) (Subramanian et al., 2005)on genes ranked based on their z-statistics, which were derived by mapping their MAST p-values to the standard normal distribution and using the sign of the log2(fold change) of their expression values to represent the direction of regulation. Technically, the z-statistic of each gene was calculated as: qnorm(p) × sign(lfc) where p and lfc are p-value and log2(fold change) of this gene obtained from MAST output, qnorm is the standard normal quantile function and sign is the signum function. KEGG pathways (Kanehisa and Goto, 2000) were used as input gene sets for GSEA, and fdr cutoff of 0.25 was used to select statistically significant pathways.