Abstract
As the largest salivary gland in oral cavity, the parotid gland plays an important role in initial digesting and lubricating food. The abnormal secretory function of the parotid gland can lead to dental caries and oral mucosal inflammation. In recent years, single-cell RNA sequencing (scRNA-seq) has been used to explore the heterogeneity and diversity of cells in various organs and tissues. However, the transcription profile of the human parotid gland at single-cell resolution has not been reported yet. In this study, we constructed the cell atlas of human parotid gland using the 10× Genomics platform. Characteristic gene analysis identified the biological functions of serous acinar cell populations in secreting digestive enzymes and antibacterial proteins. We revealed the specificity and similarity of the parotid gland compared to other digestive glands through comparative analyses of other published scRNA-seq data sets. We also identified the cell-specific expression of hub genes for Sjögren syndrome in the human parotid gland by integrating the results of genome-wide association studies and bulk RNA-seq, which highlighted the importance of immune cell dysfunction in parotid Sjögren syndrome pathogenesis.
Keywords
Introduction
The parotid gland, as the largest gland of salivary glands, produces 60% to 65% of saliva under physiological conditions in humans (Usuba et al. 2014). Histologically, the parotid gland is mainly composed of serous acini, ductal epithelium, myoepithelium, and fibrous connective tissue. The saliva produced by acini is secreted into the oral cavity through the parotid duct to facilitate mastication, swallowing, and initial digestion (de Almeida Pdel et al. 2008). The impaired function of the parotid gland leads to a decrease in saliva secretion, which brings about uncontrollable dental caries and frequent mucosal inflammation (Mathews et al. 2008).
In the past few years, bulk RNA sequencing has been used to identify the genetic and molecular features of parotid gland tissue at the transcriptional level (Oyelakin et al. 2019). However, the abundant cell composition and cell heterogeneity of the parotid gland prevent further understanding of molecular nature or its biology. Single-cell RNA sequencing (scRNA-seq) provides a powerful research tool to explore the heterogeneity of complex tissues, allowing the transcriptional characteristics of cell populations to be studied at single-cell resolution. Recently, Oyelakin et al. (2019) used scRNA-seq to reveal the degree of molecular and cellular heterogeneity in parotid glands of adult mice for the first time, followed by 2 studies on the developmental trajectory and dynamic changes of the cell population during the development of submandibular glands in mice (Hauser et al. 2020; Sekiguchi et al. 2020). However, the research focusing on the cellular heterogeneity of human parotid gland tissue at a single-cell resolution remains to be reported.
In this study, we constructed a cell atlas of the human parotid gland at a single-cell resolution by scRNA-seq to better understand the molecular and genetic mechanisms of the parotid gland in humans. We identified the biological characteristics of the parotid gland as a serous digestive gland through comparative analysis. In addition, the cell type–specific expression analysis of genes related to Sjögren syndrome susceptibility revealed the key role of immune cells in the parotid gland for the pathogenesis of Sjögren syndrome.
Methods and Materials
For details on materials and methods, see the Appendix.
Results
Characterization of Human Parotid Gland Single-Cell Atlas
The parotid gland tissue was digested into a single-cell suspension and then subjected to sequencing and downstream analysis (Fig. 1A). A total of 16,052 cells were obtained and 48,361 genes were detected. The median value of feature RNA was between 500 and 2,000. After quality control, cells with mitochondrial gene expression below 25% were filtered (Appendix Fig. 1A). In addition, 2,000 hypervariable genes were selected for subsequent dimensionality reduction analysis (Appendix Fig. 1B, C). According to the classic cell gene expression markers, the cells are divided into 10 subgroups: serous acinar cell—lactoperoxidase (LPO) (Huang et al. 2021), ductal epithelial cell—keratin 19 (KRT19) (Sekiguchi et al. 2020), myoepithelium—actin alpha 2 (ACTA2) (Huang et al. 2021), fibroblast cell—collagen type I alpha 1 chain (COL1A1) (Tombor et al. 2021), vascular cell—cadherin 5 (CDH5) (Tombor et al. 2021), T cell—CD3G (Rowe et al. 2018), NK cell—killer cell lectin-like receptor F1 (KLRF1) (Zhang et al. 2020), B cell—CD79A (Zhang et al. 2020), plasma cell—marginal zone B and B1 cell-specific protein (MZB1) (Zhang et al. 2020), and myeloid cell—lysozyme (LYZ) (Devlin et al. 2021) (Fig. 1B–D). The top 5 genes with the highest expression in each category of cells were displayed through the heatmap (Appendix Fig. 1D). Also, the differentially expressed genes (DEGs) represented in each of the 10 cell subgroups are summarized in Appendix Table 1.

Construction of cell atlas of the human parotid gland. (
Hematoxylin and eosin (H&E) staining showed the basic structure of human parotid gland tissue from the same sample as scRNA-seq, with Figure 2A displaying negative surgical margin, relatively normal tissue structure, and no cancer cell metastasis. We also collected the parotid gland samples from other 3 patients. The age, gender, and diagnosis of the patients are summarized in Appendix Table 2. Immunohistochemical staining and flow cytometric analysis showed that a large number of CD45+ immune cell infiltrations were observed in these samples from different patients (Appendix Fig. 2). Through immunofluorescence staining, we found that prolactin-inducible protein (PIP) could effectively label the clustered acinar cells, and KRT19 could label the tubular ductal cells. Myofibroblasts could be labeled by ACTA2 as flat cells surrounding acini (Fig. 2B). The immunohistochemical staining results of human parotid gland tissue from the Human Protein Atlas database also verified the specificity of the marker genes of the identified parotid gland parenchymal cells. LPO, amylase alpha 2B (AMY2B), cystatin D (CST5), and aquaporin 5 (AQP5) proteins were distributed in acinar cells; E74-like ETS transcription factor 3 (ELF3), WAP four-disulfide core domain 2 (WFDC2), keratin 8 (KRT8), and keratin 19 (KRT19) proteins were distributed in ductal cells; and ACTA2, myosin heavy chain 11 (MYH11), transgelin (TAGLN), and keratin 14 (KRT14) were distributed in myoepithelium around the acini (Fig. 2C). To better understand the biological functions of the cell populations, we performed Gene Ontology (GO) enrichment analysis on the marker genes of the 3 types of parenchymal cells (serous acinar cells, ductal epithelial cells, and myoepithelium) in the parotid gland. The results suggested that serous acinar cells exhibited amylase and hydrolase activities, ductal epithelial cells were enriched in the epithelial differentiation pathway, and myoepithelium was related to muscle contraction function (Appendix Fig. 3).

Histological staining identifies the spatial expression of different marker genes in parotid gland tissue. (
Comparative Analysis of Parotid Gland Single-Cell Sequencing Data
To better explore the differences in cell heterogeneity among different salivary glands, we combined and analyzed the published human minor salivary gland scRNA-seq data (Huang et al. 2021). Mucous acinar cells were detected in the minor salivary gland, accounting for 14% of the nonimmune cells (Fig. 3A, B). Characteristic gene analysis showed that antimicrobial protein histatin 1 (HTN1), lactotransferrin (LTF), and amylase alpha 1A (AMY1A) were highly expressed in serous acinar cells, while mucin 5B (MUC5B), mucin 16 (MUC16), and protein glycosylase fucosyltransferase 2 (FUT2) were highly expressed in mucous acinar cells (Fig. 3C–E). The expression levels of these 6 genes are highlighted in Appendix Tables 3 and 4. Immunofluorescence staining results showed that HTN1, LTF, and AMY1A were significantly expressed in serous acini of the parotid gland. MUC5B, MUC16, and FUT2 were significantly expressed in mucous acini (Fig. 3F, G). The results of GO enrichment analysis showed that the genes highly expressed in the serous acinar population were related to antibacterial response, while the genes highly expressed in the mucous acinar population were related to protein glycosylation (Appendix Fig. 4). The differentially expressed genes used for GO enrichment analysis are listed in Appendix Tables 3 and 4.

Comparative analysis of human parotid gland and minor salivary gland single-cell RNA sequencing data. (
As a follow-up to the comparative analysis among different salivary glands, we next compared the cell proportion between different species. The comparative analysis of the parotid gland scRNA-seq data between human and adult mice demonstrated that the nonimmune cell proportion of the parotid glands in the 2 species shared high similarity. Serous acinar cells were the major cell type, accounting for about 60% of the nonimmune cells (Appendix Fig. 5A, B). Correlation analysis between human and mouse showed that the serous acinar cells of the 2 species have the highest similarity with a correlation coefficient of 0.52720381 (Appendix Fig. 5C). Similar marker genes are expressed in human and mouse serous acinar cells, such as PIP, LPO, CA6, and EPCAM (Appendix Fig. 5D). Through differential gene analysis, we found that there are some differential genes in serous acinar cells between mice and humans. For example, in humans, SERPINA3 and CHI3L1 are highly expressed, while in mice, MKRN1 and PLA2G7 are highly expressed (Appendix Table 5).
Comparative Analysis of Human Digestive Gland Single-Cell Data
The salivary glands are the most proximal digestive glands in digestive system, while the parotid gland is the largest salivary gland. To compare the similarities and differences between the parotid gland and other major digestive glands at a single-cell level, comparative analysis of the published scRNA-seq data for digestive organs was performed, including salivary glands, esophagus, stomach, small intestine, pancreas, and large intestine (Appendix Fig. 6). We extracted the epithelial cell population and divided the epithelial cells into 5 populations based on the reported epithelial-specific markers of the digestive system: mucous acinar cells (MUC1+), serous acinar cells (PIP+), basal epithelial cells (MKI67+), squamous epithelial cells (KRT5+), and enterocytes (ANPEP+) (Fig. 4A, Appendix Fig. 7). Serous acinar cells predominated in salivary glands (87%), pancreas (69%), and large intestine (45%), while mucous acinar cells were the main acinar cell type in the esophagus (47%), implying the lubricating function of esophageal glands. As for stomach and small intestine, the proportions of serous acinar cells and mucous acinar cells were close, respectively accounting for 43%/40% in the stomach and 20%/16% in the small intestine. In addition, squamous epithelial cells were a unique cell type of the esophagus (22%), while enterocytes were mainly found in the small intestine (62%), indicating absorption and transport functions. The proportion of epithelial cells was consistent with the histological characteristics of different digestive organs (Fig. 4B).

Comparative analysis of single-cell RNA sequencing data of human digestive organs. (
We also tested the expression patterns of digestive enzymes and antimicrobial peptides in different digestive glands. Amylase (AMY) family genes of amylases are expressed in the parotid gland and pancreas, and maltase-glucoamylase (MGAM) is expressed in small intestine. Glucosidase is expressed in a variety of digestive glands. The disaccharide enzymes, including lactase (LCT), sucrase-isomaltase (SI), and trehalase (TREH), are expressed in the small intestine. The pancreatic lipase (PNLIP) and carboxyl ester lipase (CEL) are expressed in the pancreas, and the proteases progastricsin (PGC), serine protease 1 (PRSS1), colipase (CLPS), chymotrypsin-like elastase 3A (CELA3A), carboxypeptidase A1 (CPA1), and carboxypeptidase B1 (CPB1), are expressed in varying degrees among the stomach, small intestine, and pancreas. The ribonuclease A family member 1 (RNASE1) is expressed in most digestive glands. The antimicrobial peptides HNT1 and HNT3 have been shown to be highly expressed in parotid glands, while cathelicidin antimicrobial peptide (CAMP) is highly expressed in minor salivary glands. The expression of defensin beta 1 (DEFB1) can be detected in most digestive glands except stomach (Fig. 4C, Appendix Fig. 8).
Cellular Expression of Susceptibility Genes to Sjögren Syndrome in the Parotid Gland
To better understand the contribution of different cell populations to parotid gland diseases, we explored the cellular expression patterns of reported Sjögren syndrome susceptibility genes. By querying the single-nucleotide polymorphisms (SNPs) related to Sjögren syndrome published in the genome-wide association studies (GWASs) catalog database, 52 SNPs and 47 related protein-coding genes were identified. In addition, in the bulk RNA-seq data set of salivary gland samples from patients with Sjögren syndrome, 2,218 differentially expressed genes were identified. We overlapped the risk genes in GWASs and the identified differentially expressed genes, then verified whether these genes were expressed in at least 1 cell type in human parotid gland scRNA-seq data. Finally, we got 11 hub genes (Fig. 5A), and most of these genes were expressed in immune cells (Fig. 5B, C). By retrieving the Genotype-Tissue Expression (GTEx) database, we found that the leading SNP rs3135394 (A>G) can lead to increased expression of HLA genes (HLA-DMA, HLA-DMB, HLA-DRB1, and HLA-DQB2) (Fig. 5D). For example, across the 33 tissues in the GTEx database, rs3135394 is considered to be positively correlated with the expression of the HLA-DMA gene (Fig. 5E).

Integrated analysis of genome-wide association studies (GWASs) and bulk RNA sequencing (RNA-seq) results for Sjögren syndrome with human parotid gland single-cell RNA-seq (scRNA-seq) data. (
Discussion
Recently, Saitou et al. (2020) analyzed the bulk RNA sequencing data of 3 salivary glands (parotid gland, submandibular gland, and sublingual gland) and proteome data in saliva and identified genes and proteins specific to different gland products. In the current study, we constructed a cell atlas for the human parotid gland and explored the transcriptional characteristics of different cell populations in the parotid gland. Despite a large proportion of lymphocyte infiltration, we detected all parenchymal cell types of the parotid gland, including acinar cells, ductal epithelial cells, and myoepithelium. Although we observed adipocytes in the H&E staining of parotid gland tissue, no adipocytes were detected in the scRNA-seq data, which is consistent with the published scRNA-seq data of salivary glands (Oyelakin et al. 2019; Sekiguchi et al. 2020; Huang et al. 2021). A probable reason is that mature adipocytes are rich in lipid droplets and large in size, which limits the detection of adipocytes by microfluidic-based scRNA-seq technology (Deutsch et al. 2020). For adipocyte scRNA-seq study, the nuclei are usually extracted instead of the whole adipocytes (Rajbhandari et al. 2019). Further studies are needed to determine the biological function of adipocytes in the parotid gland.
The function of saliva has been reported to be related to lubrication, antibacterial function, digestion, and tooth integrity maintenance (Proctor 2016). In different parotid gland cell populations, we identified the most abundantly expressed genes and the encoded secreted proteins. For example, in serous acinar cells, salivary amylases AMY1A and AMY1B, as well as pancreatic amylases AMY2A and AMY2B, are highly expressed. Also, salivary peptides, including statherin (STATH), histatins, cystatins, and proline-rich proteins, which have been reported to attach to the surface of tooth enamel and form a protective pellicle to increase the stability of hydroxyapatite and augment enamel remineralization (Edgar 1990; Yin et al. 2006), are highly expressed. As for the antibacterial function, histatins (Kavanagh and Dowd 2004) and cystatins (Ganeshnarayan et al. 2012), which have antibacterial activity, are highly expressed in serous acinar cells. In addition, we found high expressed immunoglobulins in B cells and plasma cells, which are involved in the immune defense function of saliva (Chang et al. 2021). Combined analysis with human minor salivary gland showed serous acinar cells were the shared predominant populations in the 2 salivary glands, in contrast to the minor salivary glands, in which mucous acinar cells accounted for approximately one-fourth of serous acinar cells. Although very few mucous acinar cells were observed in the parotid gland after integrating the scRNA-seq data of the parotid gland and minor salivary gland, one should be very cautious in concluding that mucous acinar cells exist in the parotid gland due to the possibility of several cells being clustered as other cell types when integrating scRNA-seq data. Based on the previous study regarding minor salivary glands, serous acinar cells are rich in antibacterial proteins, including β-defensins, lactoferrin, and lysozyme (Stoeckelhuber et al. 2016). Similarly, we also observed that multiple genes encoding antibacterial peptides were highly expressed in serous acinar cells. In addition, several mucins and protein glycosylases were expressed in mucous acinar cells, which revealed that mucous acinar cells participated in protein glycosylation, a critical process of mucin production.
The salivary glands are the most proximal digestive glands in the digestive system. As the largest digestive gland in the oral cavity, the parotid gland can secrete a variety of digestive enzymes and lubricate food to facilitate the subsequent digestion process. Subsequently, the food is processed by a series of digestive glands. Comparative analysis proved that both the parotid gland and the pancreas were mainly composed of serous acinar cells. During embryonic development, both salivary glands and pancreas undergo branching morphogenesis involving budding and ductal morphogenesis as stratified epithelium (Shih et al. 2013; Wang et al. 2017). In addition, correlation analysis of bulk RNA-seq also suggested the similarity between the salivary gland and the pancreas at the transcriptome level (Oyelakin et al. 2019). Our study also verified the similarity between the parotid gland and the pancreas from the perspective of cell composition.
We further compared the expression of genes encoding digestive enzymes and antimicrobial peptides in different glands. Different expression patterns were observed among the digestive glands. For example, amylases AMY1A and AMY1B were exclusively expressed in parotid glands, while AMY2A and AMY2B were highly expressed in the pancreas. Lipase was only expressed in the pancreas, while disaccharase was solely expressed in small intestines. As for antimicrobial peptides, we found that salivary glands expressed the most abundant types of antimicrobial peptides among these digestive glands. The HTN family includes antimicrobial peptides mainly secreted by the parotid and submandibular glands, which were involved in promoting endothelial cell migration and angiogenesis (Khurshid, Najeeb et al. 2017; Torres et al. 2017). It is proposed that the C-terminal amino acid sequence SNYLYDN, which is common for HTN-1, -2, and -3, is the minimal region necessary for promoting cell migration and the wound-healing process (Sun et al. 2009). In addition, CAMP (also known as LL-37) was mainly expressed in minor salivary glands, which was proven to regulate the immune response and promote angiogenesis (Khurshid, Naseem et al. 2017; Yu et al. 2018; Shen et al. 2019). Local application of recombinant CAMP can effectively accelerate the wound-healing process (Ramos et al. 2011). Thus, we speculate that the abundant expression of antimicrobial peptides in salivary glands may be related to the role of saliva in accelerating wound healing.
Sjögren syndrome is characterized by immune cell infiltration and acinar atrophy in the salivary gland. The pathogenesis of Sjögren syndrome has not been fully elucidated yet. A large-scale GWAS has identified susceptibility genes for Sjögren syndrome, but the cell-specific contribution of these genes remains to be identified (Imgenberg-Kreuz et al. 2021). By integrating the GWAS and bulk RNA-seq data sets for Sjögren syndrome, as well as the cellular expression of these genes in parotid gland scRNA-seq data, we identified 11 hub genes for Sjögren syndrome pathogenesis, most of which were expressed in immune cells. Genetically, the HLA locus is considered the most susceptible genetic variant to Sjögren syndrome (Nezos and Mavragani 2015; Mingueneau et al. 2016; Rivière et al. 2020). Also, HLA genes have also been shown to be highly expressed in the salivary glands of patients with Sjögren syndrome (Mingueneau et al. 2016; Rivière et al. 2020). The leading SNP rs3135394 exhibits an expression quantitative trait loci (eQTL) effect and is predicted to affect the expression of multiple HLA genes in the GTEx database. The HLA genes are highly expressed in B lymphocytes, and the activation of B lymphocytes is believed to be closely related to Sjögren syndrome (Nocturne and Mariette 2018). In addition to HLA genes, STAT4 (in T cells), PTTG1 (in T cells), and GTF2I (in T and B cells) have also been identified as potential hub genes. The risk gene STAT4, expressed in T cells, has been shown to promote inflammation and increase Th1 cell differentiation (Yang et al. 2020). STAT4 blockade in patients with Sjögren syndrome may reduce the infiltration of inflammatory cells in salivary glands, thereby reducing the destruction of acinar and the progression of the disease. Also, another risk gene, GTF2I, has recently been reported to be upregulated in patients with Sjögren syndrome. The upregulation of GTF2I can activate the NF-κB pathway and promote IL6 expression (Shimoyama et al. 2021). Therefore, inhibition of GTF2I may also become a new therapeutic target for Sjögren syndrome. The identification of hub genes in our study provides potential therapeutic targets for the treatment of Sjögren syndrome.
In summary, our study constructs a single-cell atlas of the human parotid gland through scRNA-seq, explores cell heterogeneity in the parotid gland, and compares its similarities and specificities with the minor salivary gland and other digestive glands. Furthermore, we reveal the important role of the parotid gland in secreting digestive enzymes and antibacterial peptides. In addition, by integrating the GWAS and bulk RNA-seq data sets, we identify the potential hub genes for Sjögren syndrome in the parotid gland, which reveal the important significance of immune cell dysfunction to its etiology. In addition, due to the limitation of clinical sample collection involving older donors undergoing head and neck surgery, a large proportion of immune cell infiltration has been observed. Further comparative study on parotid gland tissue from young donors is still needed.
Author Contributions
M. Chen, contributed to conception, design, data analysis, drafted and critically revised the manuscript; W. Lin, contributed to conception, design, data analysis and interpretation, drafted the manuscript; J. Gan, contributed to data acquisition and analysis, drafted the manuscript; W. Lu, M. Wang, contributed to data analysis, drafted the manuscript; X. Wang, J. Yi, Z. Zhao, contributed to data conception and design, critically revised the manuscript. All authors gave final approval and agree to be accountable for all aspects of the work.
Supplemental Material
sj-docx-1-jdr-10.1177_00220345221076069 – Supplemental material for Transcriptomic Mapping of Human Parotid Gland at Single-Cell Resolution
Supplemental material, sj-docx-1-jdr-10.1177_00220345221076069 for Transcriptomic Mapping of Human Parotid Gland at Single-Cell Resolution by M. Chen, W. Lin, J. Gan, W. Lu, M. Wang, X. Wang, J. Yi and Z. Zhao in Journal of Dental Research
Footnotes
Acknowledgements
The authors thank OE Biotech Company (Shanghai, China) for providing scRNA-seq and Dr. Yao Lu for assistance with bioinformatics analysis.
A supplemental appendix to this article is available online.
Declaration of Conflicting Interests
The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the National Natural Science Foundation of China (grants 81801018 and 81771048), Sichuan Science and Technology Program (grant 2020YFS0170), and Construction of Science and Technology Innovation Base (grant 2021ZYD0104). The authors declare no potential conflicts of interest with respect to the authorship and/or publication of this article.
Data and Materials Availability
The scRNA-seq data set generated in this study has been uploaded to the GEO database (GSE188478). All the scripts and codes for all analyses in this study have been uploaded to the Github website (https://github.com/miao-OvO/PG-scRNA-seq). For publicly available data sets, mouse parotid gland scRNA-seq data can be downloaded from the GEO database with accession number GSE132867 (Oyelakin et al. 2019). The human minor salivary gland scRNA-seq data can be downloaded from the COVID-19 Cell Atlas (https://www.covid19cellatlas.org/) (Huang et al. 2021). The scRNA-seq data of human digestive glands (including large intestine, stomach, esophagus, small intestine, and pancreas) can be downloaded in the GEO database with accession numbers GSE159929 (He et al. 2020) and GSE155698 (Steele et al. 2020). The bulk RNA-seq data of the salivary glands of patients with Sjögren syndrome and healthy controls have been downloaded from the GEO database with accession number GSE159574 (
).
References
Supplementary Material
Please find the following supplemental material available below.
For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.
For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.
