High Glucose Causes Distinct Expression Patterns of Primary Human Skin Cells by RNA Sequencing

Diabetes-related skin problems represent the most common long-term complications in diabetes mellitus patients. These complications, which include diabetic dermopathy, diabetic blisters, necrobiosis lipoidica diabeticorum, and eruptive xanthomatosis, may dramatically impair patients’ quality of life and cause long-lasting disability. However, the cellular and molecular mechanisms linking diabetes-related hyperglycemia and skin complications are still incompletely understood. To assess the role of the various skin-cell types in hyperglycemia-induced skin disorders, we performed RNA sequencing-based transcriptome analysis, measuring gene expression patterns in biological replicates in normal- and high glucose-stimulated skin cells. Three primary human skin-cell types were examined, i.e., epidermal keratinocytes, dermal fibroblasts, and dermal microvascular endothelial cells. For each separate cell type, we identified gene expression. Comparing gene abundances and expression levels revealed that transcription profiles exhibit distinct patterns in the three skin-cell types exposed to normal (i.e., physiological) glucose treatment and high (i.e., supraphysiological) glucose treatment. The obtained data indicate that high glucose induced differential gene expression and distinct activity patterns in signaling pathways in each skin-cell type. We are adding these data to the public database in the hope that they will facilitate future studies to develop novel targeted interventions for diabetic skin complications.

Introduction

Diabetes mellitus (DM) is a highly prevalent metabolic disorder characterized by hyperglycemia and serious complications (1). Diabetes affects ∼9.5% of the world population, and its prevalence was estimated to increase by >50% by 2030 (2). The global increase in diabetes can be expected to pose medical, social, and economic burdens. At present, significant progress is being made in the treatment of diabetes. In addition, extensive attention is being paid to the management of organ complications caused by long-term hyperglycemia, including those in kidneys, eyes, heart, nerves, and skin. Diabetes-related dermal conditions include necrobiosis lipoidica, acanthosis nigricans, digital sclerosis, diabetic dermopathy, skin infections, and ulcers (3). Furthermore, delayed wound healing of diabetes is a major medical problem in the current treatment of diabetic patients and is receiving widespread attention in medical research. Although it is clear that certain skin lesions become more common and are among the first clinically observed signs of diabetes, the exact molecular mechanisms responsible for diabetic skin remain incompletely understood.

The skin is a complex organ, with a structure consisting of two main tissue layers: the overlying epidermis, which primarily contains keratinocytes, and the underlying dermis, a connective tissue that comprises collagens and fibers, along with multiple cell types such as fibroblasts, endothelial cells, macrophages, and mast cells (4, 5). Previous studies demonstrated that these skin cells are responsible for maintaining healthy skin tissue homeostasis and are critically involved in wound healing. Furthermore, more recent data indicate that dermal fibroblasts, microvascular endothelial cells, and epidermal keratinocytes may also influence the outcome of skin diseases in diabetes, in particular, diabetic wound healing (6). However, the exact molecular mechanisms whereby diabetes alters the physiological processes of the skin have not been elucidated. Developmentally, the epidermis originates from the ectoderm, while the dermis emerges from the mesoderm. Although closely related, epidermal keratinocytes, dermal fibroblasts and dermal microvascular endothelial cells have distinct developmental origins, fulfilling unique physiological roles. Thus, there is a need to illuminate how individual cell types function within skin tissue, and to clarify if induced changes in gene expression show similar or distinct patterns across the three cell types in diabetic skin.

Furthermore, patients with DM are also more likely to suffer from various cutaneous manifestations, including skin tightening, thickening, itch, infection, and impaired wound healing (7, 8). Although previous studies demonstrated that certain molecules are abnormally expressed in three cell types of skin under diabetic conditions, other studies indicated that some molecules were expressed differently and inconsistently by the three cell types. One study showed that exposure to high glucose (HG) was followed by increased expression of Notch intracellular domain in fibroblasts, microvascular endothelial cells, and keratinocytes (2). In addition, increased expression of MMP9 under hypoglycemic conditions has been detected not only in human skin fibroblasts and microvascular endothelial cells but also in keratinocytes (9, 10). However, comparisons of normal and diseased skin showed that some genes showed different patterns of differential expression in the three cell types. For instance, following 24 h of incubation with dexamethasone, estrogen receptor-beta was upregulated in cultured epidermal keratinocytes, whereas it was downregulated in dermal fibroblasts (11). Moreover, stromal cell-derived factor-1, which is mainly secreted by dermal fibroblasts and acts as a chemokine, functions as a mitogen in the epidermis, particularly in keratinocyte proliferation (12). Observations such as these suggest that the effects of hyperglycemia on skin fibroblasts, microvascular endothelial cells, and keratinocytes may differ too. This raises the question to which extent shifts in gene expression under HG-induced conditions in these three cell types overlap.

To address this issue, we used RNA sequencing (RNA-Seq) technology to investigate genome-wide changes in gene expression in major skin-cell types during stimulation with HG and normal glucose (NG). For three of these cell types, i.e., dermal fibroblasts, microvascular endothelial cells, and epidermal keratinocytes, we identified several genes that are expressed at physiological glucose levels. Furthermore, we generated comprehensive datasets of genes and pathways whose transcriptional profiles change under HG conditions. Overall, this study was intended to provide an abundant resource that may help improve the understanding of molecular and cellular mechanisms underlying diabetes-related dermatological problems. It may furthermore guide future research of molecular targets for novel therapeutic intervention of diabetic skin complications.

Materials and Methods

Cell Culture and Sample Preparation

Primary HDF (Cat # 2320), HDMEC (Cat # 2020), and HEK (Cat # 2110) were obtained from ScienCell and cultured in their special medium provided by ScienCell. Primary HDFs were cultured in a fibroblast medium (8 mM glucose), supplemented with 2% fetal bovine serum (FBS) and 1% fibroblast growth supplement. Primary HDMECs were cultured in an endothelial cell medium (5.55 mM glucose), supplemented with 5% FBS and 1% endothelial cell growth supplement. Primary HEKs were cultured in a fibroblast medium (6 mM glucose) supplemented with 1% keratinocyte growth supplement. All media were supplemented with 100 IU/ml penicillin. All types of the cells were maintained at 37°C in a 5% CO₂ humidified incubator, and cells from the fourth to ninth passages were used for experiments. In the experiments, cells were initially plated on 3-cm plates in their special medium. We uniformly define the glucose concentration of their special medium as normal glucose (NG). At 80% confluence, the cells were switched to HG (30 Mm) medium for 24 h. For each condition, three biological replicates were generated.

Total RNA Extraction and RNA Sequencing Procedures

RNAs from cultivated skin cells were extracted using RNAiso Plus (Cat # 9108, TaKaRa). Extracted RNA was measured and checked for concentration, purity, and integrity of RNA on a Fragment Analyzer Bioanalyzer using a standard sensitivity RNA analysis Kit. RNA-Seq was then performed on a commercially available service (service ID # F19FTSCCWLJ6421, BGI, Wuhan, China). Briefly, after poly (A)-based mRNA enrichment with oligo (dT) magnetic beads using 18 total RNA samples as the input material, the mRNA was fragmented. After the first and second strands, cDNA synthesis, end-repair, 3′-adenylation, and ligation of the fork-tailed adapters were performed. Following PCR amplification, library-quality was verified on a Fragment Analyzer bioanalyzer for a fragment size of 150 bp and for concentration. Paired-end libraries were sequenced using the BGISEQ-2000 (2 × 150 nucleotide read length) with a mean sequencing coverage depth of up to 300 ×, and coverage of the target region close to 99.9%. Meanwhile, the FASTQ data has been deposited in the CNGB Nucleotide Sequence Archive (CNSA) repository (CNP0000999).

Quality Control and Alignment of Sequencing Datasets

Sequencing of cDNA libraries generated over 41 to 53 million paired reads per sample. Low-quality reads of FASTQ files were filtered using SOAPnuke software v1.5.2 to obtain clean reads. After filtering, the high-quality reads were aligned against the human reference genome (GRCh38.p12) using an RNA-Seq-spliced read mapper HISAT2 (version 2.0.4). FPKM values for gene expression levels were calculated for annotated RefSeq genes using RSEM (version 1.2.8). To estimate the heterogeneity of our samples, we performed Pearson correlation coefficient analysis and PCA, both using R software, to map clusters within the gene expression profile from all 18 samples. Normalization and differential expression analyses were estimated from count data using DESeq2. The heat map analysis was performed with the “pheatmap” R package with default parameters selecting the DEGs after FDR-corrected p values (q values ≤ 0.05) and fold changes ≥1.2.

Functional Enrichment Analysis

GO and KEGG enrichment analysis of DEGs was performed using the Dr. TOM approach, which is an online database of BGI. The GO and KEGG functional enrichment apply the hypergeometric test to find GO term and Pathway that are significantly enriched in candidate genes’ map to the entire genome background. Subsequently, q value was conducted. We performed bubble diagram to show the significantly enriched KEGG signaling pathway and GO. X-axis represents the Rich Ratio, which means the ratio of selected gene number annotated to a particular item to the total gene number annotated in this item. The calculating formula is Rich Ratio = Term Candidate Gene Num/Term Gene Num. Y-axis represents GO term or KEGG Pathway. The GO terms and KEGG pathways of which q value ≤ 0.05 were defined as significantly enriched.

Protein–Protein Interaction Network Construction and Analysis

To estimate the interactive relationship of DEGs, we constructed PPI networks using the online tool STRING (13, 14).

Statistics

Results are presented as fold change. GO terms or KEGG pathways with q values (corrected p values) ≤0.05 were considered to be significantly enriched by candidate genes. All the figures were performed using GraphPad Prism 7.0 and Adobe Illustrator.

Results

Exprimental Design of RNA Sequencing

Human dermal fibroblasts (HDFs), human dermal microvascular endothelial cells (HDMECs), and human epidermal keratinocytes (HEKs) are the main contributors to diabetic dermopathy; therefore, we used these three cell types as models to assess the impact of hyperglycemia on the condition of the skin. To investigate transcriptome abundance, RNA-Seq was performed in primary HDF, HDMEC, and HEK exposed to either NG or HG (30 mM) medium for a duration of 24 h (Supplemental Figure 1A). To ensure the quality of the RNA obtained for sequencing, we evaluated its purity, concentration, and integrity (Supplemental Table 1). Next, RNA library construction, sequencing, and enrichment analysis were performed (Supplemental Figure 1B). The final dataset contained 41 to 53 million raw reads per sample. Following the filtering of the raw reads, sequencing showed that the average sample RNA-Seq library resulted in 42 million clean reads with an average mapping rate of 92% when aligned to human genomes (Supplemental Table 2).

Feature Analysis of the Full Sequencing Data

For each cell type, the biological duplicates yielded reproducible normalized gene expression results showing Pearson correlation coefficients (R) in the following range: 0.90<R<0.99 (Supplemental Figure 2A), indicating excellent experimental specificity and technical reproducibility. To explore the ontogeny and relationships among the three cell types, we performed dimension reduction and clustering using Principal Component Analysis (PCA) (Supplemental Figure 2B). Reassuringly, the biological replicates from each cell type were found to cluster. Next, we used the Fragments per Kilobase of exon per Million mapped reads (FPKM) method to normalize the dataset of RNA-Seq expression levels. To identify the distribution of gene expression levels in each sample, we created Boxplots showing the degree of dispersion of the data distribution (Supplemental Figure 2C). To better delineate patterns of gene abundance, we created density maps of FPKM, thereby visualizing more clearly the range of gene expression levels (Supplemental Figure 2D). To visualize the number of genes within the different FPKM ranges, we examined three ranges, namely FPKM ≤1, FPKM 1 ~10, and FPKM ≥10 (Supplemental Figure 2E). As evident from the figure, the numbers in the FPKM ≤1 range account for a large proportion of gene expression variability. Overall, these outcomes confirmed the reliability of our sequencing data.

Bioinformatic Analysis of Ubiquitously Expressed Genes in the Three Skin-Cell Types

To investigate the differences and similarities among HDF, HDMEC, and HEK, we performed transcriptomic profiling of these cell types under physiological conditions. Genes were classified as ubiquitously expressed or specifically expressed via entropy filtering. We found that a larger proportion of the identified genes were globally rather than specifically expressed, based on a total of 14,411 identified genes that were expressed in all cell types taken together (Figure 1A). To gain insight into the biological roles of the globally expressed genes, gene ontology (GO)-term and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were conducted. At the level of the GO cellular component (GO_CC), the top three enriched items included cytoplasm, nucleoplasm, and cytoskeleton (Figure 1B), whereas protein binding, transferase activity, and nucleic acid binding were most enriched in the GO molecular functions component (GO_MF) (Figure 1C). At the level of GO biological process (GO_BP), phosphorylation, cell cycle, and regulation of transcription by RNA polymerase II were highly enriched (Figure 1D). KEGG analysis suggested that the most significantly altered pathways comprised metabolic pathways, pathways in cancer, and thermogenesis (Figure 1E).

FIGURE 1

Figure 1 Bioinformatics analysis of ubiquitously expressed genes among three kinds of human skin cells. (A) Venn diagrams depicting the numbers of genes identified in HDF, HDMEC, and HEK under normal glucose stimulation. (B–D) Genes identified in all type cells are classified as ubiquitously expressed, whereas those identified in single-cell types were classified as cell type-specific. GO term analyses were performed for ubiquitously expressed genes. And GO analysis containing three sub-ontologies. CC for cellular component (B), MF for molecular function (C), and BP for biological process (D). The color and size of the bubble represent the enriched q value and the number of genes, respectively. (E) Bubble chart of top 20 KEGG enrichments involving ubiquitously expressed genes.

Bioinformatic Analysis of Cell Type-Specific Expression Genes

In addition to the analysis of the genes globally expressed across the three cell types, we included genes that were specifically expressed, reflecting heterogeneity among the three cell types. In the GO_CC category, HDF specifically expressed several genes involved in the extracellular region, extracellular space, and integral components of the plasma membrane. In the GO_MF category, HDF cells highly expressed genes involved in chemokine activity, CCR chemokine receptor binding, and cytokine activity, whereas in the GO_BP category, genes involved in G protein-coupled receptors, multicellular organism development, and regulation of signaling receptor activity were represented strongly (Figure 2A). Cell-specific gene representation in HDMEC was significantly enriched in cell type-specific genes related to integral components of the cell membrane, extracellular region, and extracellular space for the GO_CC category, whereas, cytokine activity, CCR chemokine receptor binding, and carbohydrate binding were mainly enriched for the GO_MF category. Within GO_BP, highly enriched terms included regulation of signaling receptors, G protein-coupled receptor, and lymphocyte chemotaxis (Figure 2B). It is worth noting that genes specifically expressed in HEK were significantly enriched in the cornified envelope, intermediate filament, and keratin filament in the GO_CC category. Meanwhile, GO_BP-related genes in HEK were specifically enriched in structural molecule activity and G protein-coupled receptor activity; moreover, serine-type peptidase activity was most enriched among molecular functions and for the GO_BP category, the most enriched terms were cornification, keratinization, and epidermis development (Figure 2C). We further compared the KEGG pathway of specifically expressed genes among the different cell types. The results showed that the cell type-specific expression pattern of HDF was significantly enriched in cytokine–cytokine receptor interaction, neuroactive ligand–receptor interaction, and chemokine signaling pathway, which is similar to HDMEC. In contrast, the cell type specificity of HEK was significantly enriched in neuroactive ligand-receptor interaction, retinol metabolism, and chemical carcinogenesis (Figure 2D).

FIGURE 2

Figure 2 Bioinformatics analysis of cell type specificity of HDF, HDMEC, and HEK genes. (A–C) GO term analysis of genes specifically expressed in HDF (A), HDMEC, (B) and HEK (C) cells. Enriched GO terms for cellular component, molecular function, and biological process are shown. The color and size of the bubble represent the enriched q value and the genes number, respectively. (D) Bubble chart of top twenty KEGG enrichments involving specifically expressed genes. Left, HDF; middle, HDMEC; right, HEK.

Identification of Overlapping, Differentially Expressed Genes in Human Dermal Fibroblast, Human Dermal Microvascular Endothelial Cells, and Human Epidermal Keratinocytes

Given the different gene expression patterns in the three skin-cell types under physiological conditions, we wondered whether HG-induced shifts in gene expression are also similar in these cell types. Therefore, we compared the effects of HG stimulation on DEGs across cell types. When subjected to filtering with statistical tests (q value ≤ 0.05) and fold change thresholds (≥1.2), the HDF was found to have a total of 2,122 DEGs, whereas the HDMEC and HEK cells yielded 393 and 860 of such genes, respectively. Overall, the majority of these DEGs were observed exclusively in single-cell types, while only 16 genes were found equally in all cell types (Figure 3A). To elucidate the biological processes regulated by HG levels, we performed a network analysis on the 16 overlapping genes whose expression was changed in HDF, HDMEC, and HEK. We gathered protein–protein interaction data among the 16 genes by means of the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (15). We determined densely connected genes underlying protein–protein interactions within the network to identify nodes predicted to be involved in specific regulatory processes using functional enrichment analysis. This identified cathepsin D (CTSD) protein as the most highly connected nodes among 16 overlapping genes in three cell types (Figure 3B). Subsequently, we provided histograms of the overlapping genes expression, which were up- and downregulated in response to the HG treatment, separately displayed per cell line and direction of change (Figure 3C). Remarkably, Keratin 17 (KRT17) is the only gene whose expression was upregulated consistently in all cell types, while the others were regulated in diverse directions. To understand better the functions of overlapping genes, we performed GO enrichment analysis and KEGG pathway classifications. For the GO_CC category, the top three enriched items were lysosomal lumen, endoplasmic reticulum, and lysosome (Figure 3D), whereas serine-type endopeptidase activity, homoaconitate hydratase activity, and MHC class II receptor activity were the most enriched in the GO_MF category (Figure 3E). We furthermore observed high values for GO_BP, such as MHC class I and class Ib protein complex assembly and chondrocyte development (Figure 3F). Moreover, these genes were mainly enriched in antigen processing and presentation and lysosome signaling pathway (Figure 3G).

FIGURE 3

Figure 3 Bioinformatics analysis of overlapping DEGs among three cell types. (A) Venn diagram mapping intergroup comparisons (HDF_NG vs. HDF_HG, HDMEC_NG vs. HDMEC_HG, HEK_NG vs. HEK_HG) showing the number of differentially expressed genes with a q value ≤0.05 and fold change ≥1.2) (B) Protein interaction network of 16 overlapping genes showing mutual interactions. Each node indicates the individual gene, and edges indicate the protein–protein interactions among the overlapping genes. (C) Histogram depicting the logbase two-fold change for overlapping genes in the three cell types, showing up- and downregulation in separate columns. The orange columns represent HDF; red columns represent HDMEC; green columns represent HEK. (D–F) GO enrichment analysis of overlapping genes in the three cell types. Enriched GO terms for cellular component (D), molecular function (E), and biological process (F) are shown. (G) KEGG analysis of overlapping genes across the three cell types.

A strong relationship between DEGs and cell type-specific expression genes has been reported in previously published studies (16, 17). This urges us to study the cell type-specific gene expression changes in our sequencing data. However, the results suggested that a handful of DEGs was cell type-specific in each skin-cell type following HG stimulation. And the proportion of overlapping genes in cell type-specific expression genes was 3.91% in HDF, 0.83% in HDMEC, and 2.56% in HEK respectively (Supplemental Figure 3). Thus, the subsequent bioinformatic analyses of HG-induced DEGs in each cell type were performed to elucidate the potential mechanisms underlying diabetic skin complications.

Identification of High Glucose-Induced Changes in Human Dermal Fibroblast

To gain comprehensive, more precise insights into the effects of HG levels on HDF, we analyzed the transcription profiles of dermal fibroblasts in NG and HG conditions and screened for DEGs. Of the resulting genes, we selected those whose expression increased or decreased by more than 1.2-fold in HG-induced HDF, with a q value of ≤0.05; this yielded 2,122 DEGs (Figure 4A), 1,077 of which were upregulated and 1,045 downregulated. Next, a hierarchical cluster analysis of biological replicates was performed to visualize the experimental reproducibility of the selected genes (Figure 4B). Genes are color-coded according to their normalized read counts, which is a relative measure for gene expression, and the highlighted red areas and blue areas represent the high and low expression, respectively. Next, the identified genes were subjected to GO-annotation and KEGG pathway analysis. Cytoplasm, cytosol, and cytoskeleton were the most enriched terms in the GO_CC category (Figure 4C). The most highly represented terms of the GO_MF category were protein binding, transferase activity, and hydrolase activity (Figure 4D). In the GO_BP category, the cell cycle, DNA replication, and immune system process were the most highly enriched items (Figure 4E). KEGG pathways, including ribosome, DNA replication, and cell cycle signaling pathways, were among the most enriched (Figure 4F).

FIGURE 4

Figure 4 Bioinformatics analysis of DEGs in HG-induced HDF. (A) Volcano plots representing up- and downregulated DEGs of HDF in response to high glucose stimulation conditions. The size of each data point has been plotted against the negative log₁₀ of the q value × log₂ of the fold change, q value ≤0.05, and fold change ≥1.2. Red dots represent upregulated genes; blue dots represent downregulated genes. (B) Heat map of normalized expression values across all samples for genes differentially expressed between HDF_HG (n = 3) and HDF_NG (n = 3). Red and blue represent high and low expression, respectively. (C–E) GO enrichment analysis of differentially expressed genes in normal and high glucose-induced HDF. Enriched GO terms are shown for cellular component (C), molecular function (D) and biological process (E). (F) KEGG enrichment of differentially expressed genes in HDF.

Identification of High Glucose-Induced Changes in Human Dermal Microvascular Endothelial Cells

Next, we performed analyses on HDMEC, which exhibited the smallest number of DEGs. A total of 393 genes were filtered to be differentially expressed between HDMEC_NG and HDMEC_HG, among which 173 were upregulated and 220 were downregulated (Figure 5A). Furthermore, hierarchical cluster analysis was performed based on normalized read counts reflecting gene expression, following log z-transformation in the row direction (Figure 5B). To identify the biological functions in which these DEGs partake, we performed GO enrichment analysis. In the GO_CC category, the most abundant genes were associated with ribosome, extracellular region, and extracellular space (Figure 5C). Within the GO_MF category, the genes were mostly associated with the structural constituent of ribosome, protein binding, and integrin binding (Figure 5D). Moreover, SRP-dependent co-translational protein targeting the membrane, translational initiation, and nuclear-transcribed mRNA catabolic process were the main groups within the GO_BP category (Figure 5E). A KEGG pathway enrichment analysis was also carried out, and its results showed that most of the DEGs were enriched in the ribosome, PI3K–Akt signaling, and focal adhesion signaling pathways (Figure 5F).

FIGURE 5

Figure 5 Bioinformatics analysis of DEGs in HDMEC following HG stimulation. (A) Venn plots depicting upregulated (red) and downregulated (blue) DEGs of HDMEC across two glucose-stimulation conditions. (B) Heat map of 393 genes that were differentially expressed between HDMEC_HG (n = 3) and HDMEC_NG (n = 3). Red and blue represent high and low expression, respectively. (C–E) GO enrichment analysis of differentially expressed genes in normal and high glucose-induced HDMEC. Enriched GO terms are shown for cellular component (C), molecular function (D), and biological process (E). (F) KEGG enrichment of differentially expressed genes in HDMEC.

Identification of High Glucose-Induced Changes in Human Epidermal Keratinocytes

Subsequently, we examined DEGs involved in HEK cells. We finally identified 860 DEGs, 820 of which were upregulated, and 40 were downregulated (Figure 6A). To visualize the experimental reproducibility, we performed hierarchical clustering (Figure 6B). To evaluate their biological roles, GO enrichment analysis and KEGG pathway classifications were performed. The GO_CC category was most enriched for the integral component of the membrane, the membrane, and the endoplasmic reticulum (Figure 6C). For the GO_MF, the protein binding, transferase activity, and glycotransferase activity were highly represented (Figure 6D). Additionally, these DEGs were similarly enriched in GO_BP categories, such as the lipid metabolic process, neutrophil degranulation, and transmembrane transport (Figure 6E). The KEGG pathway enrichment analysis demonstrated that most of the genes were enriched in the lysosome, protein processing, and sphingolipid metabolism (Figure 6F).

FIGURE 6

Figure 6 Bioinformatics analysis of DEGs affected by HG in HEK. (A) Volcano plots representing significant DEGs (q value ≤0.05, fold change ≥1.2) in HEK_HG versus HEK_HG. Red dots represent upregulated genes; blue dots the downregulated ones. (B) Heat map of 860 genes differentially expressed between HEK_HG (n = 3) and HEK_NG (n = 3). Red represents high expression, and blue represents low expression. (C–E) GO enrichment analysis of differentially expressed genes in normal and high glucose-induced HEK. Enriched GO terms are shown for cellular component (C), molecular function (D) and biological process (E). (F) KEGG enrichment of differentially expressed genes in HEK.

Discussion

In this study, we generated abundant transcriptome datasets that can be interrogated for a detailed understanding of the molecular and cellular bases of differences between normal and diabetic skin. We compared the gene expression profiles of separated cohorts of HDF, HDMEC, and HEK, which revealed distinct differential gene expression patterns. Specifically, for physiological conditions, we found that a larger proportion of the identified genes were expressed uniformly instead of specifically by three cell types. In addition, we globally compared transcriptome profiles in skin cells upon stimulation with two different levels of glucose. We observed many DEGs in the different cell populations. Importantly, it is plausible that certain genes and pathways are differentially regulated across cell types because of their distinct roles in dermal homeostasis. Moreover, these data may provide further insights into the transcriptional organization of the skin cells as well as skin biology and -development. In addition, predictive biomarker panels could be identified that could improve tissue analyses of normal and diabetic skins.

When transcriptome profiles were compared across the different skin-cell types under physiological conditions, no significant differences were found. KEGG analysis revealed global patterns of gene expression associated with metabolism and thermogenesis, which presumably enable the skin to perform normal physiological functions. Clearly, genes specifically expressed in a certain cell type are likely to be closely related to its unique physiological functions. Importantly, in the examination of the list of KEGG pathways, we noted that fibroblasts and endothelial cells were associated with immune system functioning (namely, cytokine–cytokine receptor interaction and chemokine signaling pathway). Consistent with previously published studies, the abnormal expression of cytokines, such as IL-20 and IL24, was also found in the skin of patients with dermopathy (18, 19). In contrast, the clearest trend in keratinocytes was a significant enrichment of specific genes in pathways related to metabolism. Our observations in the skin cells correspond to the anatomical and physiological diversity of the skin cells across dermal and epidermal tissue layers. Moreover, while the membrane and extracellular region were the major cellular components of fibroblasts and endothelial cells, the cornified envelope and intermediate filament were found to be a critical cellular component for keratinocytes, consistent with previous reports (20).

We then compared cells stimulated by two glucose levels (HG vs. NG) and observed a distinct transcriptional variation within all cell types. Unexpectedly, the three cell types had only 16 DEGs in common, most prominently in the lysosome signaling pathway. Interestingly, KRT17 is the only gene that was upregulated consistently in all three examined cell types, while the other identified genes were regulated in diverse directions in the different cell types. The lysosome was initially considered the cell’s “suicide bag” that degrades and processes cellular waste (21). Lysosomes are now recognized as signaling hubs that are crucially involved in energy metabolism and are key regulators of cell homeostasis (22, 23). However, the potential impact of lysosomal signaling in diabetic skin requires further verification. Surprisingly, the proportion of overlapping genes to DEGs was low, i.e., 0.75%. In view of this observation, we analyzed the DEGs of all cell types separately. A clear pattern observed in fibroblasts was the significant enrichment of DEGs in the cell cycle and DNA replication signaling pathways. During the proliferative phases of wound healing, fibroblasts migrate into the wound site and differentiate into myofibroblasts, which synthesize and remodel the extracellular matrix to promote wound contraction (24). The observed effects of hyperglycemia on the cell cycle and DNA repair may be one of the mechanisms responsible for delayed wound healing in diabetes. The clearest trend found in fibroblasts was the significant enrichment of DEGs in the inflammation and infectious disease pathways. The family of 2′, 5′ oligoadenylate synthetases (OAS; 2′, 5′ AS), including OAS1, OAS2, OAS3, and OASL, are generally upregulated in HDF (Supplemental Figure 4). Comparing with previous study, the expression levels of OAS3 and OASL were also significantly increased in whole blood of diabetic patients (25). Meanwhile, this finding is in line with results obtained by Wu et al. (26), and our results further indicate that this upregulation mainly occurs in the fibroblasts. Moreover, consistent with data from the literature (10), we detected an elevated MMP9 expression in HDF after stimulation by HG levels (Supplemental Figure 5). Previous studies have observed increased MMP9 in chronic wounds (10, 27), which is thought to produce prolonged inflammation in diabetic wound healing models. The clearest trend in microvascular endothelial cells was the significant enrichment of DEGs in the immune and endocrine system pathways, such as the phosphatidylinositol 3-hydroxy kinase/protein kinase B (PI3K/AKT) signaling pathway. PI3K/AKT signaling is significant in the normal wound healing process and plays a vital role in the proliferation, regeneration, remodeling, and reepithelialization of injured tissue (28–30). In addition, much evidence points to a critical role of the PI3K/AKT pathway in wound healing in diverse pathological states, including hyperglycemia (31, 32). We found that the DEGs in HEK were significantly enriched in the metabolic pathway, especially lipid metabolism and glycan metabolism, and that the formyl peptide receptor 1 (FPR1) was significantly upregulated (Supplemental Figure 6). FPR1 is a member of the family of G protein-coupled receptors and plays a role following tissue damage (33). The FPR1 signal induces the production of reactive oxygen species (34), which is also a mediator of necroptosis (35). In a more recent study, Saito et al. reported that necroptosis in keratinocytes can be triggered via the RIP1/RIP3 complex by the interaction of Annexin A1 and FPR1 (36). Another study indicates that the level of FPR1 in diabetic mice is elevated and directly impaired glucose homeostasis (37). Annexin A1/FPR1 are candidate markers of disease occurrence and may be promising therapeutic targets.

Our study identified KRT17 as the only gene whose expression was upregulated in all three skin-cell types. KRT17 is a type I intermediate filament, which is expressed in the basal cells of epithelia. As a multifaceted cytoskeletal protein, KRT17 regulates multiple biological processes, including cell proliferation, migration, differentiation, inflammatory, and immune responses. Abnormal expression of KRT17 has been found in various dermatoses, including psoriasis, alopecia, and pachyonychia congenita (PC) (38). The genetic mutation in KRT17 is related to PC, which mainly manifests itself as hyperkeratosis of the nails, palm, and pelma (39). Meanwhile, microarray analysis of non-healing chronic ulcers identified that the KRT17 gene was significantly downregulated relative to expression levels in healing chronic ulcers (40), suggesting that dysregulation of KRT17 is involved in defective healing. Another research found KRT17 was significantly enriched in diabetic foot ulcers when compared to oral mucosa (41). In addition, keratin supplementation has been shown to speed wound healing and epithelialization in mice (42, 43). In general, abnormal expression of KRT17 was found to play a key role in a variety of skin diseases. This time, we first identified the overexpression of KRT17 in diabetic skin. In conclusion, KRT17 holds great potential as a future target for the treatment of diabetic dermopathy.

Here, we conducted transcriptomic profiling on various skin-cell types and found diverse patterns of gene expression both in NG and HG stimulation. However, the generalizability of these results is subject to certain limitations. Firstly, our experiments selected the major cellular components of the skin, including keratinocytes, fibroblasts and endothelial cells. However, it has been reported that besides these cells, melanocytes also play a prominent role in dermopathy (44, 45). Secondly, the experiments in this study were only performed in cell-based in vitro models. The differences between in vitro and in vivo need to be considered, further work needs to be done to establish whether the significant genes and pathways of our findings were involved in real-life conditions. Thirdly, the three cell types were cultured in different base media, which may cause certain influence on the results. Therefore, these findings need to be interpreted with caution in skin problems associated with diabetes mellitus. Moreover, even if 16 overlapping genes were differentially expressed in all the three examined skin-cell types, most of them exhibited inconsistent directions of gene expression changes, which is far beyond our expectations. In addition, previous studies have often focused on entire tissues, containing complex mixtures of cells. Whether such approaches provide accurate information about transcription patterns in different cells is debatable. Because of technological advances, more and more scholars have reverted to single-cell sequencing, through which, researchers have become aware of the heterogeneity of the cells and started to pay attention to the behavior of single cells and cell types rather than entire tissues. Some challenges remain in single-cell RNA-Seq, in which measurements are highly sensitive to the efficiency of required enzymatic steps. Additionally, the required amplification from single cells in RNA-Seq may introduce considerable errors (46). To conclude, we performed transcriptome analysis on three skin-cell types and compared their gene expression patterns. Our analyses may support a series of hypotheses about the genes and pathways that cause diabetic skin lesions. Some genes within these pathways have already been validated as regulators of diabetic dermopathy, whereas other predictions remain to be tested. We anticipate that panels of genes derived from studies of this type will serve as useful biomarkers for clinical diagnosis.

Conclusion

Normal- and high glucose levels induced cell type-specific shifts in gene expression in three primary skin-cell types, i.e., epidermal keratinocytes, dermal fibroblasts, and dermal microvascular endothelial cells. This study uncovers a molecular mechanism of diabetes-related dermatological problems.

Data Availabilty Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics Statement

The study was reviewed and approved by the Animal Ethics Committee of Huazhong University of Science and Technology.

Author Contributions

YL and QL conceptualized the study. SZ, ZK, and QL developed the methodology. SZ, ZK, CY, PZ, and HJ were in charge of the investigation. SZ, ZK, PZ, LC, and QL wrote the manuscript. QL acquired the funding. CY, YL, and QL provided the resources. CY, YL, and QL supervised the study. All authors contributed to the article and approved the submitted version.

Funding

This work was supported financially by grants from the National Natural Science Foundation of China (No.81900432), Natural Science Foundation of Hubei Province (2019CFB499), and Science Foundation of Union Hospital (02032018-244).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Shuai Tian (BGI group, China) for the helpful advice during the conceiving of the experiment and supporting.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fendo.2021.603645/full#supplementary-material

References

1. Quondamatteo F. Skin and diabetes mellitus: what do we know? Cell Tissue Res (2014) 355:1–21. doi: 10.1007/s00441-013-1751-2

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Zheng X, Narayanan S, Sunkari VG, Eliasson S, Botusan IR, Grunler J, et al. Triggering of a Dll4-Notch1 loop impairs wound healing in diabetes. Proc Natl Acad Sci USA (2019) 116:6985–94. doi: 10.1073/pnas.1900351116

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Lima AL, Illing T, Schliemann S, Elsner P. Cutaneous Manifestations of Diabetes Mellitus: A Review. Am J Clin Dermatol (2017) 18:541–53. doi: 10.1007/s40257-017-0275-z

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Sorrell JM, Caplan AI. Fibroblast heterogeneity: more than skin deep. J Cell Sci (2004) 117:667–75. doi: 10.1242/jcs.01005

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Jevtić M, Löwa A, Nováčková A, Kováčik A, Kaessmeyer S, Erdmann G, et al. Impact of intercellular crosstalk between epidermal keratinocytes and dermal fibroblasts on skin homeostasis. Biochim Biophys Acta (BBA) Mol Cell Res (2020) 1867:118722. doi: 10.1016/j.bbamcr.2020.118722

CrossRef Full Text | Google Scholar

6. Xu F, Zhang C, Graves DT. Abnormal cell responses and role of TNF-alpha in impaired diabetic wound healing. BioMed Res Int (2013) 2013:754802. doi: 10.1155/2013/754802

PubMed Abstract | CrossRef Full Text | Google Scholar

7. de Macedo GM, Nunes S, Barreto T. Skin disorders in diabetes mellitus: an epidemiology and physiopathology review. Diabetol Metab Syndr (2016) 8:63. doi: 10.1186/s13098-016-0176-y

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Demirseren DD, Emre S, Akoglu G, Arpaci D, Arman A, Metin A, et al. Relationship between skin diseases and extracutaneous complications of diabetes mellitus: clinical analysis of 750 patients. Am J Clin Dermatol (2014) 15:65–70. doi: 10.1007/s40257-013-0048-2

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Liang Y, Yang C, Lin Y, Parviz Y, Sun K, Wang W, et al. Matrix metalloproteinase 9 induces keratinocyte apoptosis through FasL/Fas pathway in diabetic wound. Apoptosis (2019) 24:542–51. doi: 10.1007/s10495-019-01536-w

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Lerman OZ, Galiano RD, Armour M, Levine JP, Gurtner GC. Cellular dysfunction in the diabetic fibroblast: impairment in migration, vascular endothelial growth factor production, and response to hypoxia. Am J Pathol (2003) 162:303–12. doi: 10.1016/S0002-9440(10)63821-7

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Pomari E, Dalla Valle L, Pertile P, Colombo L, Thornton MJ. Intracrine sex steroid synthesis and signaling in human epidermal keratinocytes and dermal fibroblasts. FASEB J (2015) 29:508–24. doi: 10.1096/fj.14-251363

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Quan C, Cho MK, Shao Y, Mianecki LE, Liao E, Perry D, et al. Dermal fibroblast expression of stromal cell-derived factor-1 (SDF-1) promotes epidermal keratinocyte proliferation in normal and diseased skin. Protein Cell (2015) 6:890–903. doi: 10.1007/s13238-015-0198-5

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, et al. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res (2019) 47:D607–13. doi: 10.1093/nar/gky1131

PubMed Abstract | CrossRef Full Text | Google Scholar

14. He F, Qi G, Zhang Q, Cai H, Li T, Li M, et al. Quantitative Phosphoproteomic Analysis in Alpha-Synuclein Transgenic Mice Reveals the Involvement of Aberrant p25/Cdk5 Signaling in Early-stage Parkinson’s Disease. Cell Mol Neurobiol (2020) 40:897–909. doi: 10.1007/s10571-019-00780-7

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Szklarczyk D, Franceschini A, Wyder S, Forslund K, Heller D, Huerta-Cepas J, et al. STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res (2015) 43:D447–52. doi: 10.1093/nar/gku1003

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Velmeshev D, Schirmer L, Jung D, Haeussler M, Perez Y, Mayer S, et al. Single-cell genomics identifies cell type-specific molecular changes in autism. Science (2019) 364:685–9. doi: 10.1126/science.aav8130

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Hu Y, Ranganathan M, Shu C, Liang X, Ganesh S, Osafo-Addo A, et al. Single-cell Transcriptome Mapping Identifies Common and Cell-type Specific Genes Affected by Acute Delta9-tetrahydrocannabinol in Humans. Sci Rep (2020) 10:3450. doi: 10.1038/s41598-020-59827-1

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Keermann M, Koks S, Reimann E, Abram K, Erm T, Silm H, et al. Expression of IL-36 family cytokines and IL-37 but not IL-38 is altered in psoriatic skin. J Dermatol Sci (2015) 80:150–2. doi: 10.1016/j.jdermsci.2015.08.002

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Ratsep R, Kingo K, Karelson M, Reimann E, Raud K, Silm H, et al. Gene expression study of IL10 family genes in vitiligo skin biopsies, peripheral blood mononuclear cells and sera. Br J Dermatol (2008) 159:1275–81. doi: 10.1111/j.1365-2133.2008.08785.x

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Marshall D, Hardman MJ, Nield KM, Byrne C. Differentially expressed late constituents of the epidermal cornified envelope. Proc Natl Acad Sci USA (2001) 98:13031–6. doi: 10.1073/pnas.231489198

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Ballabio A. The awesome lysosome. EMBO Mol Med (2016) 8:73–6. doi: 10.15252/emmm.201505966

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Settembre C, Zoncu R, Medina DL, Vetrini F, Erdin S, Erdin S, et al. A lysosome-to-nucleus signalling mechanism senses and regulates the lysosome via mTOR and TFEB. EMBO J (2012) 31:1095–108. doi: 10.1038/emboj.2012.32

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Tapia D, Jimenez T, Zamora C, Espinoza J, Rizzo R, Gonzalez-Cardenas A, et al. KDEL receptor regulates secretion by lysosome relocation- and autophagy-dependent modulation of lipid-droplet turnover. Nat Commun (2019) 10:735. doi: 10.1038/s41467-019-08501-w

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Willer MK, Carroll CW. Substrate stiffness-dependent regulation of the SRF-Mkl1 co-activator complex requires the inner nuclear membrane protein Emerin. J Cell Sci (2017) 130:2111–8. doi: 10.1242/jcs.197517

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Cabrera SM, Chen YG, Hagopian WA, Hessner MJ. Blood-based signatures in type 1 diabetes. Diabetologia (2016) 59:414–25. doi: 10.1007/s00125-015-3843-x

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Wu C, Chen X, Shu J, Lee CT. Whole-genome expression analyses of type 2 diabetes in human skin reveal altered immune function and burden of infection. Oncotarget (2017) 8:34601–9. doi: 10.18632/oncotarget.16118

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Mirastschijski U, Impola U, Jahkola T, Karlsmark T, A.G MS, Saarialho-Kere U. Ectopic localization of matrix metalloproteinase-9 in chronic cutaneous wounds. Hum Pathol (2002) 33:355–64. doi: 10.1053/hupa.2002.32221

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Peiris TH, Ramirez D, Barghouth PG, Oviedo NJ. The Akt signaling pathway is required for tissue maintenance and regeneration in planarians. BMC Dev Biol (2016) 16:1–15. doi: 10.1186/s12861-016-0107-z

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Chen G, He L, Zhang P, Zhang J, Mei X, Wang D, et al. Encapsulation of green tea polyphenol nanospheres in PVA/alginate hydrogel for promoting wound healing of diabetic rats by regulating PI3K/AKT pathway. Mater Sci Eng C Mater Biol Appl (2020) 110:110686. doi: 10.1016/j.msec.2020.110686

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Xiao W, Tang H, Wu M, Liao Y, Li K, Li L, et al. Ozone oil promotes wound healing by increasing the migration of fibroblasts via PI3K/Akt/mTOR signaling pathway. Biosci Rep (2017) 37:1–11. doi: 10.1042/BSR20170658

CrossRef Full Text | Google Scholar

31. Liu C, Ke X, Wang Y, Feng X, Li Q, Zhang Y, et al. The level of netrin-1 is decreased in newly diagnosed type 2 diabetes mellitus patients. BMC Endocrine Disord (2016) 16:1–5. doi: 10.1186/s12902-016-0112-z

CrossRef Full Text | Google Scholar

32. Zhu KX, Nie SP, Li C, Gong D, Xie MY. Ganoderma atrum polysaccharide improves aortic relaxation in diabetic rats via PI3K/Akt pathway. Carbohydr Polymers (2014) 103:520–7. doi: 10.1016/j.carbpol.2013.12.080

CrossRef Full Text | Google Scholar

33. Kim SD, Kim JM, Jo SH, Lee HY, Lee SY, Shim JW, et al. Functional Expression of Formyl Peptide Receptor Family in Human NK Cells. J Immunol (2009) 183:5511. doi: 10.4049/jimmunol.0802986

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Wentworth CC, Alam A, Jones RM, Nusrat A, Neish AS. Enteric commensal bacteria induce extracellular signal-regulated kinase pathway signaling via formyl peptide receptor-dependent redox modulation of dual specific phosphatase 3. J Biol Chem (2011) 286:38448–55. doi: 10.1074/jbc.M111.268938

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Wang Z, Jiang H, Chen S, Du F, Wang X. The Mitochondrial Phosphatase PGAM5 Functions at the Convergence Point of Multiple Necrotic Death Pathways. Cell (2012) 148:228–43. doi: 10.1016/j.cell.2011.11.030

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Saito N, Qiao H, Yanagi T, Shinkuma S, Nishimura K, Suto A, et al. An annexin A1–FPR1 interaction contributes to necroptosis of keratinocytes in severe cutaneous adverse drug reactions. Sci Trans Med (2014) 6:245ra95. doi: 10.1126/scitranslmed.3008227

CrossRef Full Text | Google Scholar

37. Wollam J, Riopel M, Xu Y-J, Johnson AMF, Ofrecio JM, Ying W, et al. Microbiota-Produced &lt;em<N&lt;/em<-Formyl Peptide fMLF Promotes Obesity-Induced Glucose Intolerance. Diabetes (2019) 68:1415. doi: 10.2337/db18-1307

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Yang L, Zhang S, Wang G. Keratin 17 in disease pathogenesis: from cancer to dermatoses. J Pathol (2019) 247:158–65. doi: 10.1002/path.5178

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Samuelov L, Smith FJD, Hansen CD, Sprecher E. Revisiting pachyonychia congenita: a case-cohort study of 815 patients. Br J Dermatol (2020) 182:738–46. doi: 10.1111/bjd.18794

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Charles CA, Tomic-Canic M, Vincek V, Nassiri M, Stojadinovic O, Eaglstein WH, et al. A gene signature of nonhealing venous ulcers: potential diagnostic markers. J Am Acad Dermatol (2008) 59:758–71. doi: 10.1016/j.jaad.2008.07.018

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Sawaya AP, Stone RC, Brooks SR, Pastar I, Jozic I, Hasneen K, et al. Deregulated immune cell recruitment orchestrated by FOXM1 impairs human diabetic wound healing. Nat Commun (2020) 11:4678. doi: 10.1038/s41467-020-18276-0

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Konop M, Sulejczak D, Czuwara J, Kosson P, Misicka A, Lipkowski AW, et al. The role of allogenic keratin-derived dressing in wound healing in a mouse model. Wound Repair Regener (2017) 25:62–74. doi: 10.1111/wrr.12500

CrossRef Full Text | Google Scholar

43. Pechter PM, Gil J, Valdes J, Tomic-Canic M, Pastar I, Stojadinovic O, et al. Keratin dressings speed epithelialization of deep partial-thickness wounds. Wound Repair Regener (2012) 20:236–42. doi: 10.1111/j.1524-475X.2012.00768.x

CrossRef Full Text | Google Scholar

44. Kingo K, Aunin E, Karelson M, Ratsep R, Silm H, Vasar E, et al. Expressional changes in the intracellular melanogenesis pathways and their possible role in the pathogenesis of vitiligo. J Dermatol Sci (2008) 52:39–46. doi: 10.1016/j.jdermsci.2008.03.013

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Reemann P, Reimann E, Ilmjarv S, Porosaar O, Silm H, Jaks V, et al. Melanocytes in the skin–comparative whole transcriptome analysis of main skin cell types. PloS One (2014) 9:e115717. doi: 10.1371/journal.pone.0115717

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Klein AM, Mazutis L, Akartuna I, Tallapragada N, Veres A, Li V, et al. Droplet barcoding for single-cell transcriptomics applied to embryonic stem cells. Cell (2015) 161:1187–201. doi: 10.1016/j.cell.2015.04.044

PubMed Abstract | CrossRef Full Text | Google Scholar