Gua Sha is an ancient Chinese traditional remedy that has been used for about 2,000 years and is still quite popular in Asia and worldwide among practitioners of traditional medicine to treat colds, fever, flu, heat stroke, respiratory and digestive problems, chronic pain, and other health issues^[1–8]. This safe, inexpensive, and simple technique, also known as scraping, consists of repeated stroking of lubricated skin with a soft-edged instrument to induce the appearance of transient subcutaneous petechiae caused by extravasation of blood, which typically resolves within several days^[4].

In addition to the longstanding support for its effectiveness rooted in the extensive empirical knowledge dating back to ancient times, the effectiveness of Gua Sha has also been demonstrated using randomized controlled trials in numerous studies. For example, several studies from China and the West reported that this technique had beneficial effects on chronic neck and lower back pain with no adverse effects^[9–12]. Furthermore, in a randomized controlled trial of 119 subjects, Gua Sha was more effective than standard procedures in reducing multiple symptoms of diabetic peripheral neuropathy while showing no adverse effects^[13]. A meta-analysis of five randomized controlled trials strongly supported that co-therapy with Gua Sha and modern medicine had a significant positive effect on perimenopausal syndrome, reflected by the changes in the Kupperman Menopausal Index Score, serum levels of follicle-stimulating hormone, and other criteria^[14].

However, despite the very long history and proven effectiveness of this technique, little is known about the molecular mechanisms underlying the effects of Gua Sha on the human body. Several studies on humans or model animals have shown that the healing effects of Gua Sha could be mediated by an increase in immune response and decrease in inflammation, based on measuring changes in the levels of pro-inflammatory and immunosuppressive cytokines^[6,15,16]. However, to our knowledge, no genome-level work has been conducted to study any changes that occur in response to Gua Sha in the transcriptome or other types of molecules. This contrasts with other traditional Chinese medicine techniques, such as acupuncture and moxibustion, which have been found to be effective for various health issues^[17–21], and have also received some attention from modern genomic methods^[22–26]. Therefore, in this study, we performed a proof-of-principle study to investigate whether Gua Sha can induce transcriptome changes that could be directly associated with this technique. Indeed, we could identify very subtle, yet consistent transcriptome changes affecting the expression of three histone genes that could have implications for the immune response and inflammation.

Nine female volunteers between 20 and 26 years of age were recruited. All volunteers were in good health, were not treated during their menstrual period, and were not taking any medications. Each volunteer was given Gua Sha treatment on her back by the same physician, a Gua Sha expert from a traditional Chinese medicine hospital who had been performing this procedure for multiple years (see Authors′ Contributions section), with a smooth-edged Gua Sha instrument on the same day. Gua Sha treatment lasted for approximately 18 min and was terminated after the appearance of red petechiae on the skin.

Peripheral blood samples from each volunteer were drawn by an experienced nurse immediately before and 3 h after the Gua Sha treatment (Supplementary Table S1, available in www.besjournal.com) into Tempus Blood RNA Tubes (Thermo Fisher Scientific) containing reagents that immediately lysed blood cells and inactivated RNase, thus stabilizing the RNA. Blood samples were stored at −80 °C before RNA extraction. All volunteers provided informed consent and the experiments were approved by the ethics review board of the School of Medicine, Huaqiao University.

Total RNA was isolated using the Tempus Spin RNA Isolation Kit (Thermo Fisher Scientific), following the manufacturer′s instructions. The RNA-seq libraries were constructed by Novogene Corporation (Beijing) after removing globin mRNA and rRNA from the total RNA preparations using the Globin-Zero Gold rRNA Removal Kit, followed by the strand-specific lncRNA-seq protocol. As a result, RNA-seq analysis included both polyA+ and polyA− RNA species. Libraries were then sequenced by the Novogene Corporation (Beijing) using the Illumina Hiseq X Ten platform with paired-end 150 bp (PE150) strategy on a 10 gigabase (GB) scale.

The expression levels of genes were estimated based on the RNA-seq data using Salmon software^[27] for the reference human transcriptome (GRCh38) from the Ensembl database^[28] and 2,721 vlincRNA (very long intergenic non-coding RNA) transcripts taken from previous publications^[29,30], as described previously^[31].

Two-tier principal component analysis (PCA) was performed for all 18 samples using the prcomp function in the R environment^[32]. Genes with very low expression levels (mean of the raw read counts across all 18 samples ≤ 1) were excluded from the analysis. The first tier of PCA (Figure 1A–C) was based on the variance stabilizing transformation^[33] (using the vst function from the DESeq2 package^[34] in the R environment^[32]) of the raw read counts of the 23,768 remaining genes. The second tier of PCA (Figure 1D–F) was performed on the same data after an additional processing step: the variance between individuals was removed as a “batch effect” prior to the analysis using the removeBatchEffect function of the limma package^[35] in the R environment^[32]. Two tiers of uniform manifold approximation and projection (UMAP) analysis^[36] were also performed based on the same data using the same pipeline as the PCA.

Figure 1.  Effects of Gua Sha on the total blood transcriptome. The figure shows the PCA plots of the study participants before and after Gua Sha treatment based on (A–C) the original expression levels of all genes and vlincRNAs or (D–F) after the individual-specific variance was removed.

Differential expression analysis was performed using three packages in the R environment, DESeq2^[34], edgeR^[37], and limma^[35], separately, with the recommended pipelines from the manuals of each package. In the analysis, the individual from whom the sample was derived was included as a covarying factor in the design to eliminate its influence, according to the manual of each package. For example, when using the DESeq2 package, the design formula became “~ individual + treatment”, in which the term “individual” represents the individual-specific transcriptome changes — either related to the Gua Sha treatment or not, but constituting the undesired covarying factor — and the “treatment” term is the factor of interest — effects of the Gua Sha treatment common to all individuals. With such a formula, the packages in the differential expression analysis try to minimize the influence of the former covarying factor and detect the changes due to the factor of interest. In the analysis using DESeq2, genes were filtered to remove those with mean raw read count across all 18 samples ≤ 1, resulting in 23,768 remaining genes; while in the analysis using edgeR and limma, the filterByExpr function from the edgeR package was employed to filter out genes with low expression level as recommended in their manuals, resulting in 16,935 remaining genes. False discovery rate (FDR) calculations were done by the packages using the Benjamini-Hochberg procedure.

Gene co-expression analysis was performed using 18 samples (9 individuals, before and after the treatment) based on the read counts of genes after the variance-stabilizing transformation and the removal of variances between individuals. Spearman′s correlation tests were performed between the expression levels of each of the three histone genes, H1-2, H1-3, and H1-4, and all other genes, using the corr.test function from the psych package^[38] in the R environment. A gene was identified as co-expressed with the histone genes if its expression level significantly correlated with any of the three histone genes under the threshold of Spearman′s ρ > 0.6 or < −0.6 and FDR < 5%. FDR was calculated using the Benjamini-Hochberg procedure, starting with the raw correlation P values and the p.adjust function in the R environment.

The GO^[39] and Reactome pathway^[40] enrichment analyses were performed using the clusterProfiler package^[41] in the R environment^[32]. Significantly enriched terms were identified by the threshold of FDR < 5%, which was calculated by the package using the Benjamini-Hochberg procedure. The absence of results from the analysis in the corresponding figures and tables indicates that no term was significantly enriched in that analysis.

Since the basal human transcriptome, as well as its changes in response to Gua Sha, would likely be influenced by multiple individual-specific factors (e.g., gender, age, nutrition, genetics, and health conditions), we tried to limit the effect of inter-individual variation by selecting participants of the same gender and a narrow age group. Volunteers were asked to fast for 12 h before the Gua Sha to remove as much of the potential transcriptome effects caused by nutritional differences prior to treatment as possible. While Gua Sha is typically performed on people with certain symptoms (e.g., colds, fever, pain), in this study, we focused on apparently healthy participants to remove confounding factors caused by differences in diseases and health conditions (e.g., colds could be caused by different pathogens and pain could have multiple underlying reasons), medications taken to treat them, and so on. Overall, we assumed that if Gua Sha can induce changes in the transcriptomes of healthy people, it can do the same in individuals with a disease. Likewise, after treatment, volunteers were asked to perform the same activity (sit and rest) to avoid as much as possible transcriptome changes caused by post-treatment activities unrelated to Gua Sha treatment.

A previous study showed that Gua Sha induces a number of pro-inflammatory cytokines in serum^[16]; therefore, we assumed that this treatment could also induce certain changes in the whole blood transcriptome. Based on the above considerations, in this study, we analyzed changes in the peripheral blood transcriptome of nine young healthy female volunteers (Supplementary Table S1) caused by Gua Sha treatment (METHODS). Two peripheral blood samples were collected from each individual: one immediately before Gua Sha treatment and one 3 h after treatment. The short time interval was chosen for the following reasons: 1) to ensure that the detected changes in the transcriptome were primary effects of Gua Sha; 2) to be able to limit the study subjects to the same activity (sitting and relaxing) to avoid detecting effects caused by differences in post-treatment activities; and 3) to be able to detect possible transient effects of Gua Sha.

To detect changes induced by Gua Sha, blood RNAs were subjected to RNA-seq analysis that preserved both polyadenylated and non-polyadenylated RNA species (METHODS). We then calculated the expression levels of all annotated human genes and the widespread class of vlincRNAs discovered by our group^[29,42]. The reason for including the latter transcripts was that their expression in whole blood reflects changes in the physiological status of an organism, such as chronological age or the presence of non-blood cancers^[31]. As shown in Figure 1, the differences between the blood transcriptome profiles were dominated by the differences between the individuals, as revealed by the PCA based on the expression levels of genes and vlincRNAs (Figure 1A–C). This observation is consistent with the large amount of inter-individual variation in transcriptome profiles that has been extensively documented in the human population^[43] and could also represent individual-specific differences in response to Gua Sha.

The consistent effect of the treatment on the blood transcriptome became obvious only after the removal of the individual-specific variance (Figure 1D–F; METHODS). This was particularly apparent in the PC2 and PC3 dimensions, as evident from the obvious separation of the samples before and after Gua Sha treatment (Figure 1E), which was absent before the removal of the variance caused by the individual-specific effects (Figure 1B). UMAP analysis revealed essentially the same results as PCA (the effects of the Gua Sha treatment were observed in UMAP3 and UMAP4 after the individual-specific variance removal; Supplementary Figure S1, available in www.besjournal.com). These results indicated that although Gua Sha treatment shifted the peripheral blood transcriptome profile, the common changes induced in all individuals were rather subtle.

Figure S1.  Effects of Gua Sha on the total blood transcriptome. The figure shows the UMAP analysis plots of the study participants before and after Gua Sha based on (A–C)the original expression levels of all genes and vlincRNAs or (D–F) after the individual-specific variance was removed.

To further characterize the blood transcriptome changes induced by Gua Sha treatment, differential expression analyses were independently performed with DESeq2, edgeR, and limma (METHODS), the 3 packages representing some of the most popular tools for differential expression analysis^[44–46]. DESeq2 and edgeR were designed for RNA-seq data, whereas limma was initially designed for microarrays and then expanded for RNA-seq data^[45–47]. The main difference among them is that DESeq2 and edgeR rely on negative binomial models to derive differentially expressed genes (DEGs) while limma uses a linear model^[45–47]. Since only few DEGs were identified with a conservative FDR cutoff of < 5% (7, 6 and 2 by DESeq2, edgeR and limma, respectively), we relaxed it to a moderate threshold of < 10%^[48] which is still in the range of FDR used in differential expression analysis^[49]. As the result, the DESeq2, edgeR, and limma packages identified 12, 22, and 5 DEGs, respectively (Figure 2A; Supplementary Table S2, available in www.besjournal.com). Three histone genes, H1-2, H1-3, and H1-4, were shared by the three DEG sets and were upregulated after Gua Sha treatment (Figure 2A and B). The expression fold changes (after vs. before the treatment) of the three histone genes were small (the log₂ fold changes of the three histone genes were ~0.30–0.33, see Supplementary Table S2), which is in line with the above finding of the subtlety of the global effects of Gua Sha on the transcriptome. However, the trends of increased expression in response to Gua Sha were statistically significant (Figure 2B). Specifically, histone H1 genes H1-2, H1-3, and H1-4 were upregulated in response to Gua Sha treatment with P values of 2.4 × 10^-7, 1.4 × 10^-9, and 2.7 × 10^-6, respectively, as calculated using the two-tailed paired Wilcoxon test (Figure 2B).

Figure 2.  Identification of the histone genes induced by Gua Sha treatment and their potential regulatory networks. (A) Venn diagram of the DEG sets identified by DESeq2, edgeR and limma methods. (B) Expression levels of the 3 histone genes, H1-2, H1-3, and H1-4, before and after the Gua Sha treatment. The P values were determined using a two-tailed paired Wilcoxon test are shown in the figure. (C and D) Venn diagrams of genes (C) positively and (D) negatively co-expressed with the 3 histone genes.

To further explore the properties of genes responding to Gua Sha treatment, we identified GO terms and Reactome pathways enriched in the 29 DEGs found by at least one package (Figure 2A) relative to the background of all genes (METHODS). The only terms and pathways identified were associated with the 3 histone genes and were related to multiple biological processes and pathways such as histone methylation, regulation of gene silencing, apoptosis, cellular senescence, and programmed cell death, in which products of these genes have been previously implicated^[50–59] (Supplementary Figure S2, available in www.besjournal.com). The diversity of processes and pathways is reflective of the fundamental biological functions of the products encoded by these three genes, which represent members of the H1 family of histone proteins that bind to linker DNA between nucleosomes to form chromatin fibers, and are thus necessary for the condensation of nucleosome chains into highly ordered structures^[51], and regulation of gene expression through epigenetic modifications, nucleosome spacing, and chromatin remodeling^{[51,52,60,61]}. Some H1 functions are not related to chromatin organization; for example, in response to DNA damage, the histone protein H1.2 can translocate from the nucleus to the cytosol where it activates the pro-apoptotic protein Bak, leading to apoptosis^[53,62–64]. The basic biological functions of the H1 family, combined with multiple reports that associate members of this family with disease^[56,65–69], prompted us to further investigate the potential mechanisms of the function of the 3 H1 histones in the Gua Sha response.

Figure S2.  Enriched GO terms and Reactome pathways of the 3 histone genes. (A) Enriched GO terms in the Biological Process (BP) category. (B) Enriched Reactomepathways. (C&D) Diagrams of the (C) “Programmed Cell Death” and (D) “Cell Senescence” Reactome pathways with the enriched pathways highlighted. These two panels are derivatives offigures from Reactome website which are licensed under Creative Comments Attribution 4.0 International (CC BY 4.0) License.

Since the products of the 3 histone genes can affect expression of other genes, to further investigate how the 3 H1 family members might mediate the effect of Gua Sha treatment, we attempted to identify potential members of their networks using a co-expression analysis (METHODS). We identified 763 and 1,020 genes that were positively and negatively co-expressed with the three histone genes, respectively (Figure 2C and D; Supplementary Table S3, available in www.besjournal.com). Interestingly, genes positively co-expressed with the histone genes showed enrichment in the cellular component GO term “T cell receptor complex” (FDR = 0.074%; Supplementary Table S4, available in www.besjournal.com). In addition, they showed enrichment in such Reactome pathway terms, as “phosphorylation of CD3 and TCR zeta chains”, “translocation of ZAP-70 to immunological synapse” and “generation of second messenger molecules” — all under the parent term of T cell receptor signaling pathway — and “co-stimulation by the CD28 family” with its child term “PD-1 signaling” (Figure 3A; Supplementary Table S5, available in www.besjournal.com). The above terms represent important signaling pathways related to the activation of the adaptive immune system. These results indicate that the induction of histone genes by Gua Sha treatment could result in the activation of T lymphocyte-mediated adaptive immunity (Figure 3B), which is consistent with the observed increase in active immune cells, including activated T lymphocytes, found in mice treated with Gua Sha^[16]. Negatively co-expressed genes were enriched in multiple GO and Reactome pathway terms related to “platelet activation, signaling and aggregation”, as well as in several GO terms of signal transduction in response to DNA damage (Figure 4; Supplementary Tables S4 and S5). In general, platelet activation is a complex phenomenon associated with hemostasis and inflammation (reviewed in^[70]) that could be relevant to the effects of Gua Sha treatment (see DISCUSSION). Additionally, 83 and 79 vlincRNAs were identified as positively and negatively co-expressed with at least one of the three histone genes, respectively (Supplementary Table S3), indicating that non-coding transcripts might also be involved in mediating the effect of Gua Sha treatment.

Figure 3.  Enriched GO terms and Reactome pathways of genes positively co-expressed with the 3 histone genes. (A) Enriched Reactome pathways. (B) Diagram of the “Adaptive Immune System” Reactome pathways with the enriched pathways highlighted. This panel is a derivative from a figure from the Reactome website which is licensed under Creative Commons Attribution 4.0 International (CC BY 4.0) License.

Figure 4.  Enriched GO terms and Reactome pathways of genes negatively co-expressed with the 3 histone genes. (A–C) Enriched GO terms corresponding to (A) BP, (B) MF and (C) CC categories. The X-axes represents the number of genes in each category. (D) Enriched Reactome pathways. (E) Diagram of the “Platelet activation, signaling and aggregation” Reactome pathways with the enriched pathways highlighted. This panel is a derivative from a figure from the Reactome website which is licensed under Creative Commons Attribution 4.0 International (CC BY 4.0) License. GO, gene ontology; BP, biological process; MF, molecular function; CC, cellular component.

In this study, we showed that transcriptome changes induced by Gua Sha were subtle. This contrasts somewhat with the longstanding empirical knowledge proving the effectiveness of this technique. However, there could be objective reasons to explain this apparent contradiction. First, Gua Sha was conducted on healthy volunteers with no apparent health issues, such as cold and fever, while Gua Sha is typically performed on people who have health problems that might exhibit stronger transcriptome changes. Second, we focused on immediate changes that could be directly attributed to Gua Sha by keeping the participants in the same controlled environment post-treatment to minimize the risk of detecting changes unrelated to Gua Sha that could occur if the volunteers were allowed to continue with their normal routines. However, it is conceivable that longer times are needed to allow transcriptome changes related to the treatment to become more apparent. Third, many of the responses to Gua Sha could be individual-specific, while we focused on common transcriptome changes across individuals. Indeed, as the analysis has shown, the latter concern was valid because most of the transcriptomic differences, even after a short interval post-treatment, had individual-specific signatures. Fourth, in this study, we performed a bulk transcriptome analysis; however, the actual changes could affect only a small fraction of specific immune cells in the blood. Therefore, single-cell transcriptome analysis of the blood might be a good alternative approach to further study the molecular changes induced by Gua Sha.

Despite these caveats, we identified three genes whose expression changes could be associated with treatment. These genes encode members of the H1 histone family of proteins that have been widely implicated in the control of chromatin state and gene expression. However, in addition to their well-known roles in maintaining chromatin structure, H1 histones can drive an inflammatory response in microglia^[71], and the expression of H1 subtypes plays a role in neutrophil differentiation^[59]. Furthermore, using co-expression analysis, we found that the three histone genes may affect the expression of genes associated with the T cell receptor complex and platelet activation involved in the immune response and inflammation. The association between Gua Sha and inflammation and the immune system found in this study using a genomics approach was also consistent with previous studies relying on traditional biochemical and cell-based assays. For example, Chen et al. found an increase in active immune cells, levels of the pro-inflammatory cytokines tumor necrosis factor-alpha (TNF-α), IL-6, and IL-1β, and a decrease in the immunosuppressive cytokine IL-10 in healthy mice following Gua Sha treatment^[16]. In contrast, a negative effect of Gua Sha on proinflammatory cytokines has been reported in other studies. For example, Yuen et al. found a reduction in TNF-α following treatment in elderly individuals with chronic lower back pain^[15]. The inhibitory effect of Gua Sha on the expression of TNF-α and other pro-inflammatory cytokines, such as IL-1β and IL-6, was also observed in rats with lumbar disc herniation induced by autologous nucleus pulposus^[6]. The differences among studies might be due to differences in the species used and the health status of the subjects. Nonetheless, these studies strongly suggest that Gua Sha exerts its effects, at least in part, by affecting the immune system, which is consistent with our findings. In addition, the enrichment of functions related to platelet activation among the genes negatively co-expressed with the three histones is consistent with a previous study in which Gua Sha was shown to improve blood flow via blood vessel expansion in mice^[16].

Overall, this study shows that changes in the transcriptome profile can occur in response to Gua Sha, and they likely represent the physiological mechanisms behind the healing effect of this procedure. However, given the subtlety of the effects and/or individual-specific variation, more extensive studies on the effects of Gua Sha on transcriptome that involve 100′s or even more participants of both genders, different age groups, and health conditions and conducted at different time points after the treatment are required to comprehensively understand the effects of Gua Sha. We hope that this work will be an important stepping stone leading to such endeavors.

KAPRANOV Philipp conceived the project and supervised the analytical and wet laboratory parts of the project. QI Fei performed all bioinformatic analyses with the help of CHEN Jun Jie and XIA Qiu. CAI Ye and HAN Xue Er organized the participation in the study and performed the molecular biology part of the project. CHEN Chun Li performed the Gua Sha treatments. KAPRANOV Philipp and QI Fei wrote the manuscript with help from CAI Ye.

The data presented in this study were deposited in the GSA-Human repository (accession number HRA002107). Custom R scripts for analysis were deposited at Github (https://github.com/qifei9/guasha_code) and are also available at Zenodo (https://doi.org/10.5281/zenodo.7046437).

All participants provided informed consent and the experiments were approved by the ethics review board of the School of Medicine, Huaqiao University.