1 Introduction
Thousands of years ago, the famous script of traditional Chinese medicine (TCM), “Huangdi Neijing” stated that “The superior doctor prevents illness; the mediocre doctor attends to impending sickness; the inferior doctor treats actual illness.” Interestingly, the concept of “preventative medicine” was acknowledged very early in human history. In the present time, to enable preventative medicinal practice, one must capture physiologic indices including genetic predispositions, environmental, and/or lifestyle influences prior to disease onset. The ancient Chinese used pulse touching, verbal inquiries, and extensive observations to capture biometric information such as those related to the smoothness of meridian channels about individuals, thereby defining their patho-physiologic states. Apparently, the present and ancient methods and concepts have to converge because they reflect the same principles about human health. Immune cells in the circulation provide surveillance in the body to monitor disturbances in organs and elicit proper immune responses to clear out pathogens and/or damaged or mutated cells. Individuals’ immunological characteristics are closely associated with their abilities to fight off all kinds of pathological insults. Interestingly, recent studies on the neuro-anatomical bases for classic acupunctures in the realm of TCM have also revealed that interactions between the autonomic nervous system and immune system are the major therapeutic targets for acupunctures [
1,
2]. Therefore, the ability to profile people’s immunological state may help unveil thousands of years’ worth of myths associated with TCM and may greatly help early diagnosis aimed at curing diseases in the early stages and/or even prevent disease onset.
One TCM theory divides people into nine different basic/major constitutions [
3], presumably reflecting their different genetic predispositions [
4], as well as environmental influences and lifestyles. People with certain constitutions, except the normal kind, are prone to the development of certain illnesses. For example, Tanshi (phlegm dampness), as one of the eight basic subhealth-constitution types, is relatively easy to identify [
5]. Individuals with the Tanshi constitution are usually overweight, prone to suffering from diabetes and other metabolic diseases [
6]. Sleep-apnea syndrome is frequently detected in Tanshi individuals [
7], who also often feel heavy and lethargic. Previous investigations have reported peripheral-blood gene expression, single nucleotide polymorphisms (SNPs) [
4], serum proteomic features [
8], and DNA methylation-related epigenetic changes in Tanshi individuals [
5]. These studies have aimed to identify biomarkers, assist in precise diagnosis, and understanding the genetic and epigenetic bases for Tanshi constitution. However, the peripheral blood comprises various heterogeneous cell populations, and gene-regulatory programs vary among different blood-cell subtypes. Immune-cell subtype compositions in human populations are also highly variable. All of these confounding factors confer extreme difficulty in identifying biomarker identification based on whole-blood samples extremely challenging.
With the development of high-throughput large scale single-cell transcriptome profiling technology, the precise profiling of the gene expression of each mono-nucleated immune-cell subtype with single-cell resolution is possible. Single-cell RNA sequencing (scRNA-seq) of peripheral blood mononuclear cells (PBMCs) not only reveals the composition of immune cell subtypes but also provides insights into subtype-specific gene expression. These findings largely reflect an individual’s immune functional state, commonly referred to as their “immune state” [
9]. In the present study, we used scRNA-seq to analyze the “immune states” of individuals with Tanshi constitution and those who did not belong to the Tanshi population. Data revealed significantly reduced peripheral-blood mucosal-associated invariable T (MAIT) cell contents in Tanshi individuals and increased CD14
+ classic monocyte population. We also observed increased expression of genes associated with classic inflammatory responses including those responding to tumor necrosis factor (TNFα) and IL1β, as well as interferon responses. These findings indicated increased systemic and perhaps chronic inflammation levels in Tanshi individuals, which appeared to be quite characteristic of the Tanshi constitution. Gene-expression programs related to cell apoptosis and hypoxic responses also appeared to be characteristic to Tanshi individuals, which may result from frequent sleep apnea associated with Tanshi constitution. By providing high-resolution data on “immune state” specific to people with different constitutions, this study began to pave a feasible path to unveil the mysterious “constitution theory” in TCM, which has been in practice for thousands of years in China and is essentially an ancient Chinese way toward preventative medicine.
2 Materials and methods
2.1 Participants
Adult volunteers aged 18–60 years were recruited to the program. All participants were categorized into Tanshi and non-Tanshi groups based on TCM diagnosis and a scaled questionnaire. The Tanshi group comprised four individuals (2 females and 2 males) with a Tanshi constitution, whereas the non-Tanshi group comprised nine individuals (6 females and 3 males).
2.2 scRNA-seq of PBMCs
The procedures of PBMC isolation and scRNA-seq were the same as previously described [
11]. PBMCs were isolated from heparinized venous blood of Tanshi (phlegm-dampness constitution) or non-Tanshi donors using Ficoll–Paque
TM (GE Healthcare) density gradient centrifugation, subsequently frozen in freezing media (70% RPMI-1640, 20% FBS, and 10% DMSO), and stored in liquid nitrogen until use. Single-cell capture and library construction were performed using the Chromium Single Cell 3′ Reagent Kits v2 (10x Genomics) according to the manufacturer’s specifications. Libraries were sequenced on an Illumina Novaseq 6000 platform.
2.3 scRNA-seq data analysis and statistics
scRNA-seq data were aligned and quantified using the Cellranger pipeline (10x Genomics) against the GRCh38 human reference genome downloaded from the 10x Genomics official website. Preliminary counts were then subjected to downstream analyses, and the procedures were same as described in a previous study [
11]. In a typical procedure, cells with less than 200 genes were filtered out, and the logarithmic normalization of raw counts and top 3000 highly variable genes (HVGs) selection were performed using Scanpy [
20]. Specific genes from HVGs including mitochondrial genes, immunoglobulin genes, and genes linked to poorly supported transcriptional models (annotated with the prefix “Rp-”) were excluded. The remaining HVGs were then subjected to principal component analysis (PCA). Harmony algorithm was used to remove batch effects [
21]. To identify clusters and cluster specific features, we used PARC approach [
22] and “FeatureSelectionByEnrichment” function from cytograph2 algorithm [
23], respectively. Subsequently, another round of PCA, Harmony, and PARC were performed, followed by k-nearest neighbors’ determination in a KNN graph, uniform manifold approximation and projection (UMAP) via Pegasuspy [
24], and clustering by PARC. Additionally, we applied Scrublet [
25] to determine potential doublets.
2.4 Comparing immune-cell proportion
We calculated immune-cell proportions for each major cell type in Tanshi and non-Tanshi PBMCs as ratios of number of cells in certain cell type divided by the total number of cells. To reveal changes in cell compositions between samples in different groups, we performed Wilcoxon test on the proportions of major cell types across different groups.
2.5 Differential expression analysis, gene set enrichment analysis (GSEA), and score signature modules
To investigate immunological-feature alterations, we identified differentially expressed genes (DEGs) using muscat algorithm [
26] with default parameters. In a typical procedure, we first collapsed the data, summing UMIs across cells for each volunteer, to produce a bulk RNA-seq style UMI profile for each sample. Afterwards, the aggregated counts were loaded onto pbDS function to identify DEGs, which were then visualized through heatmaps produced by pbHeatmap function. GSEA of DEGs (log
2FC > 0.5 and adjusted
P-value < 0.05) were performed using one-sided Fisher’s exact test (as implemented in the ‘gsfisher’ R package) with the “HALLMARK,” “KEGG,” and “REACTOME” gene sets derived from the MsigDB database. Gene sets with
P value < 0.05 were considered to be significant.
2.6 Antibody staining and fluorescent activated cell sorting (FACS) analysis
The antibodies used in this study were PerCP-Cy™5.5 Mouse Anti-Human CD3 (BD Pharmingen, 55282, clone SP34-2, 1:250), FITC Mouse Anti-Human TCR Vδ2 (BioLegend, 331406, clone B6, 1:20), APC Mouse Anti-Human TCR γ/δ (Miltenyi, 130-113-500, clone 11F2, 1:50), PE/Cyanine7 Mouse Anti-Human TCR Vα7.2 (BioLegend, 351712, clone 3C10, 1:20), PE Mouse Anti-Human CD161 (BioLegend, 339904, clone HP-3G10, 1:20), PerCp Mouse Anti-Human CD45 (BioLegend, 304025, clone HI30, 1:20), APC Mouse Anti-Human CD14 (BioLegend, 301808, clone M5E2, 1:20), PE Mouse Anti-Human CD16 (BioLegend, 302056, clone 3G8, 1:20), FITC Mouse Anti-Human CD4 (BDIS, 340133, clone SK3, 1:5), PE Mouse Anti-Human CD127 (BD Pharmingen, 557938, clone HIL-7R-M21, 1:50), and APC Mouse Anti-Human CD25 (BDIS, 662525, clone 2A3, 1:20). Aliquots of 1 × 106 PBMC or erythrocyte-lysed whole blood cells were stained with antibodies for 30 min in darkness at room temperature. Samples were then washed with 1 mL of PBS and centrifuged at 1500 rpm for 3 min. The washing step was repeated with another 1 mL of PBS. Finally, 50 μL of PBS was added to resuspend the cells for flow cytometry. Immunophenotyping was conducted using a BD FACSVerse™ Cell Analyzer, with 10 000 events recorded per run. Flow-cytometry data were gated using FlowJo (version v10.0.7) software. T cells were classified by CD3 expression, and two groups of γ/δ T cells were selected based on their expression of Vd2. Some Vd2-positive cells showed negative expression for TCR γ/δ likely due to steric hindrance; thus, they were still considered Vd2 cells. MAIT cells were identified by the co-expression of two markers, Vα7.2 and CD161. Classical and non-classical monocytes were distinguished by different expression levels of CD14 and CD16. CD4+ T cells were isolated from lymphocytes based on CD4 expression, and those with a CD25+/highCD127–/low profile were defined as CD4+ regulatory T (Treg) cells.
3 Results
3.1 scRNA-seq revealed the specific “immune state” of Tanshi individuals
The experimental design of the study can be divided into two parts, namely, the “discovery” part and the “validation” part (Fig.1). Volunteers went through the classic four TCM diagnosis methods, “Wang (looking), Wen (listening), Wen (questioning), Qie (sensing the pulse),” and a scaled questionnaire (Fig.1) for calculating “Tanshi (phlegm–dampness)” scores. Four individuals (2 females and 2 males) with “Tanshi” constitution and 9 non-Tanshi individuals (6 females and 3 males) were subjected to PBMC scRNA-seq analyses.
A total of 51 196 PBMCs from 13 individuals were sequenced with good-quality data (Fig.2). UMAP plot demonstrated the clustering of 15 different cell subtypes, namely, B cells, plasma-blasts, two types of dendritic cells (DCs; myeloid DC/mDC and plasmacytoid DC/pDC), two types of monocytes (classic CD14
+ monocytes and non-classic CD16
+ monocytes), NK cells, CD4 T helper cells, Naïve CD4 T cells, CD4 regulatory T cells (Tregs), CD8 T cells, and Naïve CD8 T cells, as well as three types of innate T cells including MAIT cells, Vd2 γδ T cells, and non-Vd2 γδ T cells. Cell cluster subtypes were defined based on cell-type specific gene expressions (Fig.2 and S1) and referenced from a large body of literature [
10–
12]. When non-Tanshi samples (colored in orange) were overlaid on Tanshi samples (colored in blue), the blue dots that did not overlap represented cells in the Tanshi samples with transcriptomes substantially different from the controls. These were referred to as Tanshi-unique cells (Fig.2). The proportion of each of the 15 immune-cell subtypes was color coded, and averaged results from Tanshi and non-Tanshi groups are presented in Fig.2. Statistical analysis of the content of each of the 15 subtypes of cells in Tanshi and non-Tanshi groups revealed statistically significant decreases in Tregs, MAIT cells in samples with Tanshi constitution (Fig.2). Regarding monocyte population, we found statistically significant increases in classic (CD14
+) monocyte content and a trend of decreases in non-classic (CD16
+) monocytes.
3.2 Validation of scRNA-seq results by FACS analysis
scRNA-seq measures gene expression only at mRNA levels, not protein levels and thus cannot be directly linked to function. However, immune-cell subtype specification based on scRNA-seq is quite reliable because the definition does not rely on single genes but a group of genes in the form of transcriptomes. Given that the sample size for the original scRNA-seq was relatively small, we decided to validate the scRNA-seq findings using FACS analysis based on cell type-specific markers for MAIT, CD4, and CD25 double positive Tregs, and CD14
+/CD16
+ PBMCs in additional Tanshi (5) and non-Tanshi (5) individuals. MAIT cells were mucosal-associated invariable T cells often residing in organs (such as intestines, lung epithelia, etc.), as well as in peripheral blood. They are a group of T cells primarily involved in innate immunity, unlike most T cells that are involved in adaptive immunity. During the COVID-19 pandemic, MAIT cells in peripheral blood reportedly dropped dramatically as they entered the lung for anti-viral action. When patients entered the convalescent period, MAIT contents in the periphery reverted [
13].
The FACS-based detection of MAIT cells is not as straightforward [
10] and involves multiple steps. After enrichment for single live cells, PBMC specimens were first selected for CD3 positive T cells (Fig.3–3C). The total CD3 positive T cell population were then excluded for γδ T cells, which were the other T cells primarily involved in innate immunity. Lastly, after excluding γδ T cells, CD161 and Vα7.2 double-positive MAIT cells were defined (Fig.3 and 3E). When the percentage of MAIT cells over total PBMC were calculated, the content of MAIT cells in non-Tanshi and Tanshi groups demonstrated statistically significant differences between groups, and Tanshi individuals had much more reduced MAIT cell contents in peripheral blood (Fig.3).
For CD14+ and CD16+ PBMC populations, the FACS analysis strategy using CD14 and CD16 immunolabeling was straightforward. Results demonstrated statistically significant differences in CD16+ PBMCs between the Tanshi and non-Tanshi groups. This result did not fully agree with the scRNA-seq finding. However, when CD14+ and CD16+ PBMCs were taken into an account, the two different types of analyses appeared to be consistent. In Tanshi individuals compared with non-Tanshi ones, decreases in the ratios of CD16+/CD14+ cells were significant (Fig.4–4G). Interestingly, the content of CD25 and CD4 double-positive Tregs in Tanshi and non-Tanshi individuals were not statistically significant (Fig. S2), suggesting that a decrease in Treg content is not a robust biomarker for Tanshi constitution.
3.3 Increased IL-1β and TNFα-based inflammatory immune states are associated with Tanshi constitution
scRNA-seq provides information on immune-cell subtype composition and cell-type-specific gene-expression changes between Tanshi and non-Tanshi individuals. We first identified DEGs by muscat algorithm with default parameters. By summing UMIs across cells for each volunteer, the data was collapsed into a pseudo-bulk RNA-seq style UMI profile for each sample. Afterwards, the aggregated counts were loaded onto pbDS function to identify DEGs, and heatmaps were plotted by the pbHeatmap function. Volcano plot demonstrated DEGs between Tanshi and non-Tanshi individuals (Fig.5). Gene ontology (GO) analyses based on biological processes demonstrated that samples from Tanshi individuals were enriched in “regulation of lipid metabolism,” “regulation of ERK1/2 cascade,” “cellular responses to TNF,” and “regulation of cell adhesion/cell growth/cytokine production,” whereas samples from non-Tanshi individuals were enriched in “regulation in defense responses,” “mitochondrial transportation,” “regulation of endopeptidase activity,” etc. These finding indicated that these processes were downregulated in Tanshi individuals (Fig.5).
scRNA-seq contains cell-type specific gene-expression information, and MAIT and monocytes are the two immune-cell subtypes mostly affected in Tanshi individuals. Accordingly, we performed DEG analyses specifically on MAIT and monocytes (Fig.5). The GSEA of DEGs (log2FC > 0.5 and adjusted P-value < 0.05) was performed using one-sided Fisher’s exact test (as implemented in the ‘gsfisher’ R package) with the “HALLMARK,” “KEGG,” and “REACTOME” gene sets derived from MSigDB. The “HALLMARK” enrichment of Tanshi-specific MAIT cells highlighted “TNFα signaling via NFκB,” “Apoptosis,” and “Hypoxia,” whereas enrichment for Tanshi monocytes were “TNFα signaling via NFκB,” “Interferon-α responses,” and “Interferon-γ responses” (Fig.5). We selected six representative genes including TNF, a main factor in “TNFα signaling via NFκB,” IRF1, IFI44L, and IFITM3, which are major interferon responsive factors, as well as SOCS3 and IL-1β, major factors in “JAK-STAT” signaling. We then plotted their expression in UMAP and found cell-type-specific changes (Fig.5). Altogether, these analyses consistently suggested that the immune-state for Tanshi individuals was heightened with inflammation indicated by elevated “TNFα-NFκB” signaling, “IL1-JAK/STAT” signaling, and interferon-mediated immune responses. Moreover, increased “Hypoxia” may result from frequent episodes of sleep apnea, and “Apoptosis” may be associated with hypoxia-related insult.
4 Discussion
Research efforts have been exerted toward delineating the biological bases underlying Tanshi (phlegm-dampness) constitution based on peripheral liquid biopsies, including genetic (SNPs), epigenetic (DNA methylation), bulk-sample based transcriptome analyses, and serum proteomic analyses [
5,
8,
14]. Some candidate biomarkers have been postulated through these studies, but confirming the putative biomarkers or revealing their associated biological/pathological processes has encountered some challenges. This study is perhaps the first one to use PBMC scRNA-seq and discover, with single-cell resolution, the immune-cell composition and molecular characteristics associated with Tanshi constitution. Multiple data analysis methods and biological confirmations by using FACS analyses consistently pointed out that Tanshi individuals compared with non-Tanshi people had reduced content of MAIT cells, as well as increased ratio of CD14
+ classic over CD16
+ non-classic monocytes. Moreover, TNFα-NFκB signaling, JAK-STAT signaling, and IFN signaling appeared to be heightened in Tanshi individuals, indicating heightened chronic inflammatory immune state associated with Tanshi constitution. Moreover, “hypoxia” and “apoptosis” appeared to be linked to Tanshi, and as aforementioned, may result from the frequent sleep-apnea episodes that Tanshi people are known to experience.
The high-throughput scRNA-seq of PBMCs has opened a new way to holistically describe a person’s immune-system function. Even a single blood sample from an individual can facilitate scRNA-seq analysis of PBMCs, yielding data from thousands of immune cells and capturing the expression profiles of thousands of genes per cell. It is almost like a “molecular microscope,” allowing for detailed scrutiny of a person’s immune functional state. With this approach, and from “big-data/deep-data” processing, reliable and interesting patterns can often be identified with relatively small-sized patient samples [
15,
16], as revealed by many COVID-19-related studies describing specific immune features linked to severe versus mild or asymptomatic individuals infected with SARS-CoV-2 [
13,
17]. Based on some reports, low MAIT cell content is associated with the development of severe symptom upon SARS-CoV-2 infection. Based on our study, Tanshi individuals had low MAIT cell content and should be prone to develop severe symptoms upon SARS-CoV-2 infection, and indeed, this was the case [
18], confirming the validity of this study.
TCM has been in practice for thousands of years. It has very ancient theoretical/philosophical roots, deep and difficult to grasp. TCM is also quite empirical, but often lacks the long-term accumulation of treatment records and statistical analyses. Apparently, at an early stage, ancient Chinese had already recognized that the human body is a heavily interconnected system, and problems occurring in one organ can be rooted in not easily captured problems in another organ. TCM places significant emphasis on modulating the immune system. Today’s medical knowledge tells us that the circulating immune system and the nervous system, particularly the autonomic nervous system, strongly connect the body to the mind and vice versa. For example, stimulation of the vagal-adreno axis through acupuncture has profound effects on the immune system [
1,
2]. Therefore, the ability to scrutinize the immune system using a “molecular microscope” with single-cell resolution ought to yield many new insights regarding the “mind–body” connection and the mechanisms underlying how TCM works and the biological correlates for TCM terms such as “Yin-Yang,” “Qi (vital energy),” “Xue (blood),” “Huo (fire),” “Shi (dampness),” etc. This study was the pilot research on the power of scRNA-seq in revealing immunological features associated with Tanshi (phlegm-dampness) constitution. The same method can certainly be used to study other constitutions, as well as dynamic changes in one’s health, such as “Shang-Huo (with elevated fire),” which can often be affiliated with allergies or acute inflammation [
19]. Moreover, using the scRNA-seq of PBMCs before and after treatment with TCM herbs can reveal what a particular TCM prescription does to the body or how it influences the immune system.
We foresee that in the near future, the usage of scRNA-seq of PBMC, together with other big-data-based omics such as metabolome, gut microbiota, and sympathetic–parasympathetic activities, will unveil the myth behind the thousand-year-old wisdom of TCM, including those associated with preventative medicine, so that we can all be “superior doctors.”