This article is mentioned in:

Introduction

Nuclear technology is widely used in various fields, such as medicine, industry and agriculture. Exposure to radiation or nuclear leakage is sometimes unavoidable and potentially catastrophic. Radiation accidents, such as those caused by the Three Mile Island nuclear power plant in the United States (1979), Chernobyl, Russia (1986), and Fukushima, Japan (2011), highlight the potential threats associated with catastrophic nuclear events (1,2). Since the middle of the 20th century, more than 400 radiological accidents have occurred (3), causing thousands of injuries. Ionizing radiation markedly affects individuals, as they are continuously exposed to radiation from various sources, including natural background radiation, diagnostic procedures and nuclear disasters. Radiation therapy, particularly at high doses, can induce inflammation and increase the risk of autoimmune reactions (4). The immune system, which is essential for defending against environmental insults and stress, is profoundly affected by ionizing radiation (5). Experimental studies have indicated that high-dose radiation can suppress immune function, whereas low-dose radiation may have a stimulating effect (6). However, the current evidence regarding the response of the immune system to ionizing radiation remains fragmented and contradictory. Therefore, it is imperative to obtain more comprehensive insights into the effects of ionizing radiation on the immune system.

Clinically, the consequences of radiation exposure are categorized into systemic radiation and localized radiation. Localized irradiation refers to exposure to a relatively limited area of the body, resulting in restricted damage to the directly exposed tissues. The harm experienced by organisms subjected to radiation begins at the skin and progresses towards the tissues in the internal body. The hands of irradiated patients demonstrate a unique and profound relationship with ionizing radiation exposure, as they remain the most affected part of the body. Hand irradiation is typically accompanied by whole-body radiation, which can induce changes in the blood system, as evidenced by incidents such as the Nanjing 192Ir accident (7,8). Clinically, the earliest symptom of localized injury is transient erythema and/or oedema, which is caused by telangiectasia and fluid extravasation within the first week. This may be followed by pruritus, stiffness, stinging, and tenderness (9). The initial appearance of the skin and the affected area does not necessarily correlate with the extent of tissue damage; microcirculatory disturbances often affect much larger areas than initially perceived. The skin contains a network of blood vessels and capillaries, and blood circulates throughout the body via a system composed of arteries, veins, and capillaries. The hematopoietic system has been identified as particularly susceptible to radiation-induced damage, affecting both mature blood cells and hematopoietic stem cells within the bone marrow niche, which play crucial roles in the replenishment of blood cells (10).

One component of healthy tissue that is consistently subjected to irradiation is the circulating peripheral blood mononuclear cells (PBMCs), which demonstrate heightened sensitivity to ionizing radiation (11-14). PBMCs constitute a diverse array of cellular subsets that are traditionally categorized into myeloid and lymphoid lineages (15,16). The myeloid lineage includes monocytes and their progeny, as well as granulocytes, such as neutrophils and eosinophils. Lymphoid cells predominantly consist of T cells, B cells, and natural killer (NK) cells (17-19). RNA is crucial in cellular biological processes, and transcriptomes provide essential information directly linked to cell phenotypes. Single-cell RNA sequencing (scRNA-Seq) is a robust technique for profiling individual cells. The conventional bulk RNA-seq technique measures average gene expression across cells in a sample and identifies differences between sample conditions. By contrast, scRNA-Seq measures the gene expression of individual cells, allowing for the identification of differences between cells in one or more samples. While cells have traditionally been characterized morphologically or by unique molecules, scRNA-Seq facilitates the automatic classification of cells through clustering of transcriptomes, enabling the identification of heterogeneous cell types and molecular states even within groups considered to consist of only one cell type. scRNA-Seq of heterogeneous cell populations has become the standard tool for establishing cellular lineages and tissue compositions in the human body (20,21). This method amplifies minute transcriptome RNA from individual cells, enabling high-throughput sequencing to obtain a comprehensive expression profile of the entire transcriptome at single-cell resolution. This approach reveals the molecular regulatory mechanisms underlying specific biological and disease processes (22-24).

The primary aim of the present study was to characterize the properties of PBMCs from both irradiated patients and their healthy counterparts via scRNA-Seq technology. This approach offered a novel research framework for further investigations into PBMCs in irradiated patients. By analyzing single cells from irradiated patients, the present study identified distinct characteristics of various cell populations and explored the potential biological functions of PBMCs in response to radiation exposure, thereby providing insights into therapeutic possibilities for treating radiation-induced injuries.

Materials and methods

Patient and healthy control blood sample collection

The irradiated patient was a 30-year-old male. The hands of the patient accidentally contacted the radioactive source 192Ir for ~30 sec on August 24, 2019. After a few days, the skin of the hands gradually became red, swollen, painful, ulcerated and sloughed. In the accident, the activity of the radioactive source 192Ir was 59 Ci and ~3 weeks later, the patient was admitted to the hospital for bone marrow aspiration, which indicated hypoplasia of the bone marrow. The physical dose estimation results were as follows: Left hand, from the little finger to the thumb, 2.2, 2.7, 3.4 and 6 Gy; right hand, from the little finger to the thumb, 3.1, 4.8, 6.7, 9.2 and 11.2 Gy, respectively. The blood samples used for scRNA-Seq were taken 113 days and 631 days after irradiation (Fig. 1B). Fresh PBMCs were isolated from the circulating blood of the irradiated patient and 10 healthy controls of matched sex and age (Tables SI, SII and SIII). The patient and the healthy controls provided informed consent for the present study. All patients provided written informed consent for their tissues to be used for scientific research. Ethical approval for the study was obtained from the Second Affiliated Hospital of Chengdu Medical College, Nuclear Industry 416 Hospital [Chengdu, China; approval no.: 2019(25)].

Figure 1 Diverse cell types in PBMCs delineated by single-cell transcriptomics. (A) Results of dicr frequencies in irradiated patient blood. (B) Schematic of the workflow showing the overall strategy of scRNA-Seq to create a peripheral blood mononuclear cell atlas and subjected to droplet-based 10X Genomics. The hands of the patient accidentally contacted the radioactive source 192Ir. (C) The t-SNE plot displays main cell types in PBMCs. Cell types are summarized in the right panel. Each dot represents only one cell. The number of single cells in each sample captured by single-cell RNA-Seq 10X Genomics. (D) Bubble plot shows the expression levels of key linage-defining genes among all clusters. The size of bubble indicates the percentage of cells expressing a specific gene and the color of bubble indicates the average level of gene expression (E) The heatmap displays the top 10 SDE genes in each cell type. Each column represents a cell, and expression values are scaled by log2fold change. (F) t-SNE plot embedding of all PBMCs from scRNA-Seq data, colored by each group. (G) The violin plot shows the gene number of each sample. (H) Venn diagrams of all expressed genes for each cell type. PBMCs, peripheral blood mononuclear cells; dicr, dose-response of induced chromosomal aberrations; scRNA-Seq, single-cell RNA sequencing; t-SNE, t-distributed stochastic neighbor embedding; SDE, differentially expressed; Con, control.

Chromosome aberration analysis

Venous blood (2 ml) was taken from the subject and placed into a heparin tube. Next, the blood was transferred to a culture container, gently mixed and placed in an incubator with RPMI 1640 medium (Gibco; Thermo Fisher Scientific, Inc.) at 37˚C for 24 h. Subsequently, 20 µl of 10 µg/ml colchicine (cat. no. ST1173, Beyotime) was added to the culture medium (resulting in a final concentration of 0.04 µg/ml) and the cells cultured for an additional 24 h before harvesting. After terminating the culture, the supernatant was gently removed from the culture bottle using a pipette. KCl (0.075 mol/l; 8 ml), pre-warmed to 37˚C, was added to each bottle. The cell clumps were pipetted and transferred to a 10-ml tube. The tube was incubated in a 37˚C constant temperature water bath for 30 min. The hypotonic centrifuge tubes were removed and add 5-10 drops of freshly prepared fixative to each tube (with a volume ratio of methanol to glacial acetic acid of 3:1). They were mixed by gently blowing with a pipette and centrifugation at room temperature 200 x g for 10 min. The clear liquid was aspirated with a pipette, 8 ml of fixative (with a volume ratio of methanol to glacial acetic acid of 3:1) added and the cell clumps mixed thoroughly by pipetting. Fixation was at room temperature for 30 min, followed by another centrifugation at room temperature 200 x g for 10 min twice. After the final centrifugation, the supernatant was discarded and 3-6 drops of fixative added to adjust the cell concentration. The cell suspension was dropped from a height of 10-30 cm onto a clean glass slide that had been pre-cooled at 4˚C and allowed to dry naturally in room temperature. The sample was stained with 10% Giemsa stain at room temperature for 8 min. The stain was gently rinsed with distilled water and the slide placed on a slide rack to dry naturally at room temperature (25). Finally, the chromosomal aberrations were analyzed under an optical microscope. Under the x10 objective, each stained chromosome slide is scanned from right to left, column by column or row by row, to search for analyzable metaphase cells. After locating the target, the cells are counted and analyzed under the oil immersion lens (magnification, x100). At least 100 metaphase cells are analyzed for each examination.

Preparation of single-cell suspension of human PMBC samples

Human peripheral blood (5 ml) was placed in an EDTA anticoagulant tube and diluted with an equal volume of 1X PBS. An equal volume of Lymphocyte Separation Solution (Ficoll; cat. no. P4350; Beijing Solarbio Science & Technology Co., Ltd.) was added to a 50-ml centrifuge tube, and the diluted blood carefully layered onto the Lymphocyte Separation Solution. It was then centrifuged at 20˚C and 400 x g using a horizontal rotor for 20 min, with the brake set to 0. After centrifugation, the middle buffy coat layer containing the PBMC cells was carefully pipetted into a new 15 ml centrifuge tube. The buffy coat cells were washed with 10 ml of 1X PBS and centrifuged at room temperature 300 x g for 10 min. The supernatant was discarded and the cells resuspended in 5 ml of 1X PBS. Following another centrifugation at room temperature 300 x g for 10 min and two additional washes, the supernatant was discarded and the cells resuspended in 1 ml of RPMI-1640 medium (Corning, cat. no. 10-040-CVR) supplemented with 0.04% BSA. Cell concentration and viability was assessed using trypan blue staining (room temperature, 3 min) to obtain a single-cell suspension.

Single-cell capture, library building and sequencing

The freshly prepared single-cell suspension was adjusted to 700-1,200 cell/µl according to the 10X Genomics Chromium Next GEM Single Cell 3'Reagent Kit v3.1 (10X Genomics, cat. no. 1000268) operation manual for computer and library construction. The constructed library was sequenced using Illumina Nova 6000 PE150 platform (Illumina, Inc.) (26).

ScRNA-Seq bioinformation analysis process

The library construction, sequencing and data analysis were conducted by Shanghai OE Biotech Co., Ltd. The raw data generated from high-throughput sequencing were in FASTQ format. The official 10X Genomics software, CellRanger (v5.0.0; 10X Genomics), was used to perform data quality statistics on the raw data and to compare it with the reference genome (https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-GRCh38-2020-A.tar.gz). The software identified the barcode markers that distinguished cells in the sequence, and each unique molecular identifier (UMI) marker corresponding to different intracellular mRNA molecules is employed to quantify high-throughput single-cell transcriptome data, yielding high-quality statistical information such as cell count, gene median, and sequencing saturation.

Following the preliminary quality control performed by CellRanger, the Seurat (v3.1.2) (27) software package was employed for additional quality control processing of the data. Theoretically, the number of genes expressed by most cells, the number of UMIs, and the expression proportion of mitochondrial transcripts will cluster within a specific range. Consequently, low-quality cells were filtered based on the distribution of three indicators: nUMI, nGene, and percent.mito. The criteria established were to retain cells with gene and UMI counts within the median ±2 times the absolute median absolute deviations (MAD) and with a mitochondrial transcript ratio of <20%, which were classified as high-quality cells. Additionally, DoubletFinder (v2.0.2) (28) software was utilized to remove doublets.

Top variable genes across single cells were identified using the method described in Macosko et al (29). The most variable genes were selected using FindVariableGenes function (mean.function=FastExpMean, dispersion.function=FastLogVMR) in Seurat. To remove the batch effects in single-cell RNA-sequencing data, the mutual nearest neighbors presented by Haghverdi et al (30) was performed with the R package batchelor (v1.3.4). Graph-based clustering was performed to cluster cells according to their gene expression profile using the FindClusters function in Seurat. Cells were visualized using a 2-dimensional t-distributed stochastic neighbor embedding (t-SNE) algorithm with the RunTSNE function in Seurat.

The SingleR (v0.2.2) (31) package was used for cell type identification. Based on the single-cell reference expression quantitative public data set, the correlation between the cell expression profiles and the reference data were calculated. This process assigned the cell type with the highest correlation in the reference data set to the cell to be identified, minimizing human subjectivity and interference. The identification principle calculated the Spearman correlation between the expression profile of each cell in the sample and the annotated cells in the reference data set. The cell type with the greatest expression correlation was chosen as the final cell type to be identified.

The FindMarkers function in the Seurat package was used to perform differential gene screening, screen out significantly differential genes based on the conditions of P<0.05 and difference fold >1.5 and Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of significantly differential genes conducted through the hypergeometric distribution test. The probability model of the hypergeometric distribution is used to assess the enrichment of gene sets in specific biological processes or pathways.

SCENIC analysis

SCENIC analysis was conducted using the RcisTarget motifs database and GRNboost [SCENIC (32) v1.1.2.2, RcisTarget v1.2.1, and AUCell v1.4.1; Bioconductor)] with default parameters. The process involves the following steps: i) Identification of potential target genes of each transcription factor based on co-expression, ii) Using the RcisTarget package to ascertain the actual transcription factors and their corresponding target genes through motif analysis and iii) Evaluation the activity of each regulator in each cell using the AUCell package. To assess the cell type specificity of each regulon, the regulon specificity score (RSS) (33) was calculated based on the Jensen-Shannon divergence (JSD) using the scFunctions package (https://github.com/FloWuenne/scFunctions/), along with the relatedness index (CSI) of all regulators.

Pseudo-time analysis

The Monocle2 (v2.9.0) (34) package was used to infer cell differentiation trajectories. The specific steps were as follows: First, the importCDS function of the Monocle2 package was used to convert the Seurat object into a CellDataSet object, and the differentialGeneTest function used to filter out the genes used to sort cells (ordering genes qval<0.01), then the reduceDimension function was used to perform dimensionality reduction clustering and, finally, the orderCells function used to infer the differentiation trajectory.

Gene set variation analysis (GSVA) enrichment analysis

Gene set files were initially obtained and structured from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database via the gene set enrichment analysis base package (v1.44.0; Kyoto University Bioinformatics Center; https://www.kegg.jp/). Subsequently, the GSVA package (35) (v1.30.0) was employed for individual cells to assess pathway activity. The scoring of pathway activity values was performed, followed by the utilization of the LIMMA software package (v3.38.3; http://www.R-project.org/) (36) to determine the variance in signaling pathway activity across distinct groups.

Results

Information on the irradiated patient

Human PBMCs were collected from an irradiated patient at different time points (113 and 631 days after irradiation). The control group included a combination of peripheral blood samples obtained from 10 healthy individuals. The irradiated patient had been diagnosed with mild bone marrow acute radiation sickness. The present study collected peripheral blood samples from the patient for examination and experimentation. Chromosomal aberrations in peripheral blood lymphocytes were detected via the peripheral blood lymphocyte culture method. The overall rate of chromosome abnormalities in patients exposed to radiation was 0.3%, which was not significantly different from that in the control group (Fig. 1A).

Analysis of the heterogeneity of PBMC populations by scRNA-Seq

Subsequently, scRNA-Seq was employed (Fig. 1B) to examine the cellular transcriptome profile of the PBMCs. Following quality control, 6,447 and 7,892 cells were obtained from the PBMCs of the irradiated patient on the 113th (Stage 1 ) and 631st (Stage 2) days post radiation, respectively. Overall, 9.101 cells were acquired from 10 healthy controls. Visualization analysis via t-distributed stochastic neighbor embedding (t-SNE) clustering revealed that 23,440 cells from the control, Stage 1, and Stage 2 groups were categorized into 17 clusters. The results of the subgroup correlation analysis revealed that the 17 clusters were roughly separated into 5 major cell types (Fig. 1C), including T cells, NK cells, B cells, neutrophils and monocytes. The markers for each cluster and cell type are shown in Fig. 1D (37). A comparison of each cluster was conducted to identify distinct gene signatures, and the heatmap shown in Fig. 1E displays the top 10 genes that were significantly differentially expressed (SDE) within each cluster. The cells were all sorted according to the quantity of genes detected, with an average of ~2,000 genes recognized in the three samples. In Stage 2, more genes were detected relative to the remaining samples, as illustrated in Fig. 1F. Separate analysis of the two groups revealed the presence of all clusters and comparable cluster marker gene expression in each group, with significant differences observed in cell numbers. It was hypothesized that T cells, B cells and NK cells were the main cells involved in radiation injury. A statistical analysis of the percentages of the five types of PBMCs was subsequently conducted (Fig. 1G) and the results revealed that the percentages of T cells, B cells and NK cells significantly changed after exposure; thus, these cells may be the main cells involved in the radiation process.

Ionizing radiation is recognized for its immunomodulatory properties, which influence both localized and systemic immune responses. The subsequent investigation focused on the immune cell distribution among patients in Stage 1 and Stage 2, in addition to healthy control subjects, as demonstrated in Fig. 1G. Notably, certain immune subpopulations, including B cells and NK cells, were found to be more prevalent in Stage 2 patients and healthy controls than in Stage 1 patients. This observation suggested a gradual restoration of the typical circulating levels of these immune cell types in irradiated patients towards a distribution akin to that of healthy individuals. T cells were present in the control group but few in the Stage 1 and Stage 2 patient. Some cell subsets, such as monocytes and neutrophils, were very few in number between the control and Stage 1 groups (Fig. 1G). Examinations revealed the existence of two distinct cellular subsets that were specifically identified in the patient subjected to irradiation, comprising a B-cell subfraction and an NK cell contingent. Given the aim of the present study of radiation-sensitive cell populations of irradiated patients, the preponderance of the subsequent examinations was directed towards a detailed evaluation of these particular cell groups.

A Venn diagram revealed the numbers of all expressed genes for each cell type. Omitting the regions of overlap, a total of 599 distinct genes were identified in NK cells, 908 genes in B cells, 835 genes in monocytes, 2,462 genes in neutrophils and 722 genes in T cells (Fig. 1H). As expected, the KEGG pathway analysis of the variably expressed genes among disparate groups revealed that the upregulated genes, as ascertained by scRNA-Seq, were associated with notably increased pathways, indicative of comprehensive immune system stimulation. These pathways include pathways related to cell adhesion molecules, the sphingolipid signaling pathway, platelet activation, the cytotoxicity mediated by natural killer cells, the signaling cascade of the T-cell receptor, and the Rap1 signaling pathway (Fig. S1A and B). The downregulated pathways included ribosome enrichment. Exposure to high levels of ionizing radiation has been demonstrated to initiate a profound reorganization of overall protein synthesis, enabling cellular energy conservation through the inhibition of unnecessary protein production (38). The process of protein synthesis, which is inherently energetically costly, results in heightened sensitivity to external stressors, including exposure to ionizing radiation (39). Using gene-based clustering techniques, the analysis identified and stratified the five cellular subsets within the samples onto a two-dimensional plane, concurrently illustrating the relative abundance of each subset across the samples (Fig. S1C). Overall, the quantitative findings indicated a scarcity of B cells in Stage 1 and a dearth of NK cells in Stage 2. Among the PBMCs of Stage 1 - and Stage 2-irradiated patients, the proportion of T cells was notably lower than that of normal controls (Fig. S1C).

Analysis of NK cell heterogeneity

NK cells, recognized alternatively as large granular lymphocytes, constitute a subset of cytotoxic lymphocytes endowed with the innate capacity to discern and eliminate deleterious cellular entities in the absence of major histocompatibility complex molecule participation and antibody mediation, thereby playing an essential role in innate immune defense mechanisms. The present study further analyzed NK cell clusters by comparing transcriptome expression patterns between radiation-exposed patients and control patients. By t-SNE analysis, the participating NK cells were divided into 9 clusters: NK1 to NK9 (Fig. 2A). The percentage distribution of each cell subpopulation in each tissue is shown in Fig. 2B. The NK1, NK2, NK6 and NK7 subsets were present mainly in Stage 1, while the NK5 subsets were present mainly in Stage 2. To further clarify the characteristics of each subgroup, the marker genes of each cluster was analyzed and the top 10 genes used to construct a heatmap (Fig. 2C).

Figure 2 scRNA-Seq analysis reveals NK cell transcriptome in irradiated patient blood. (A) t-SNE plot visualization of NK cell subclusters colored by each cluster, and between-group differences. (B) Bar plots showing the proportion of each cluster contributed by each group. Percentage of cells in Control, Stage 1 and Stage 2 for nine different subsets of NK cells. (C) The heatmap shows the top 10 SDE genes in each subset of NK cells. scRNA-Seq, single-cell RNA sequencing; NK, natural killer; SDE, differentially expressed; Con, control.

The results revealed that the NK2, NK5, NK7, NK8 and NK9 genes highly expressed IL-7R, CD3D and CD3G, which are involved in lymphocyte development, T-cell development and signal transduction. The increased expression of these genes suggested that these five cell subpopulations may be in an activated state. KLRC2 and KLRC3 are highly expressed in NK3 and NK6 cells, which are involved in immune activation (36).

Functional enrichment analysis of NK cell subsets

NK cells were most abundant in the Stage 1 subgroup, mainly in the NK1, NK2, NK6 and NK7 subpopulations. To further clarify the function and molecular markers of each subpopulation, the DEGs in each subpopulation were analyzed and further functional enrichment analysis performed via GSVA (Figs. S2 and 1A). The results revealed that thiamine metabolism and retinoid metabolism in the NK1 subgroup were markedly activated. Of intracellular thiamin ~80% is phosphorylated and most of it is bound to proteins. Most of the thiamine in serum is bound to proteins, mainly albumin. Thiamin is primarily a transporter form of vitamins (40), suggesting that the NK1 subgroup is crucial for energy metabolism. The hematopoietic cell lineage pathway of the NK2 subpopulation was markedly activated. Hematopoiesis is initiated by hematopoietic stem cells (HSCs), which have the potential for self-renewal or progression into progenitor cells with multilineage potential, such as a common lymphoid progenitor (CLP) or a common myeloid progenitor (41). A CLP gives rise to the lymphoid lineage of white blood cells or NK cells and T and B lymphocytes (42). Staphylococcus aureus infection and DNA replication pathways were significantly activated in the NK6 subpopulation. Within the NK7 subpopulation, a multitude of biological pathways exhibited significant activation. These include the ribosome pathway, pathways associated with infections by human T-cell leukemia virus-1, RNA transport mechanisms, thermogenesis, infection by human immunodeficiency virus-1, the T-cell receptor signaling pathway, oxidative phosphorylation and infection by Epstein-Barr virus, among others. The significant activation of viral infection-related signaling pathways in human blood after radiation exposure may indicate a possible relationship between human exposure to irradiation and exposure to infection.

Analysis of the developmental trajectories of NK cell subsets

To delve deeper into the developmental trajectory of cellular subsets, a pseudotime sequence analysis was conducted utilizing the scRNA-Seq dataset with the Monocle 2 algorithm. Fig. 3A shows the results of a trajectory analysis of nine relevant clusters and indicates the pseudotime direction. The trajectory of pseudotime encompasses three distinct stages, with various cellular clusters distinctly positioned along different points of the pseudotemporal continuum (Fig. 3B). The NK1, NK6 and NK7 subsets were in the late stage. The NK2, NK3, NK4, and NK5 subsets were identified at divergent stages along the pseudotemporal trajectory. The NK8 and NK9 subsets were in the early stage. The control cohort predominantly occupied the initial phase of the pseudotime trajectory, which was predominantly localized within State 1. The Stage 1 cohort was predominantly situated at the terminal stage of the pseudotime trajectory, with the majority of its distribution concentrated in State 3. The Stage 2 group was in the whole state of the pseudotime path. NK cells recovered in the Stage 2 group. For the purpose of examining the pseudotime patterns of definitive marker genes, a selection comprising the signature genes characteristic of each respective cluster was made. As shown in Fig. 3C, the heatmap of pseudotime can be classified into three modules for the marker genes. Module 1 is related mainly to apoptosis and transport and increases gradually with pseudotime. By contrast, module 3 decreased at the end of the pseudotime axis and is related mainly to ribosome processing, antigen presentation and translation. The expression of all the marker genes increased and then decreased as the pseudotime increased in module 2. This is related mainly to the inflammatory response and immunity. The data of the present study revealed dynamic gene expression profiles during the different stages of NK cells.

Figure 3 Analysis of developmental pseudo-time of NK cell subsets. (A) Pseudo-time state visualization of NK cell subclusters, colored by each cluster. Each dot represents a single cell. The black arrow indicates the start and direction of the trajectory. (B) Trajectories of NK cells by state, pseudotime and groups. (C) A heatmap displays the SDE genes during progression of NK cells. Color key from blue to red indicates relative expression levels from low to high. The GO and KEGG analyses reveal enriched functions and pathways of the SDE genes. (D) Heat map of RAS activity of regulons in each cell population. Rows indicate different regulons and columns indicate different cell populations. The color changes from blue to red indicating a low to high RAS activity score, with higher RAS scores indicating stronger regulon activity in that cell population. NK, natural killer; SDE, differentially expressed; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; RAS, regulon activity score; Con, control.

Transcription factors (TFs) and their regulated genes form intricate networks that dictate cell identity. The present study applied single-cell regulatory network inference and clustering (SCENIC) to deduce the regulon activities of the NK cell clusters (Fig. 3D), identifying key regulators and their targets. Fig. 2A shows that the number of cells in the NK1, NK2, and NK6 cell subpopulations significantly increased in Stage 1. The STAT1 signaling pathway was identified as the most actively regulated pathway in NK1, NK2 and NK6 cells (Fig. 3D). STAT1 plays a critical role in the development and functionality of NK cells, influencing various aspects such as maturation, survival, and proliferation. Activation of STAT1 can enhance the cytotoxicity and cytokine production capacity of NK cells (43). These results suggested that STAT1 may act as an essential driver of irradiation in NK cells.

Analysis of B cell heterogeneity

The participating B cells were divided into three clusters, B1-B3 (Fig. 4A). Statistical evaluation of B cell cluster proportions revealed a post irradiation decline followed by recovery to baseline two years postexposure in the irradiated cohort. The percentage distribution of each cell subpopulation in each tissue is shown in Fig. 4B. The B1 subset was predominantly observed in Stage 2. The B2 subset was present mainly in the Stage 2 and the control groups. The B3 subset was present mainly in the control group, and a lower amount was also present in Stage 1 and Stage 2. To further clarify the characteristics of each subgroup, the marker genes of each cluster were analyzed and the top 10 genes used to construct a heatmap (Fig. 4C). The findings indicated that the B2 and B3 subgroups presented increased expression of genes associated with inflammatory and immune response pathways, including IGHA1, IGHG1, TUBA1A and MYADM. Moreover, the B1 subgroup presented significant upregulation of genes such as AREG, IRS2, LY9, NR4A3 and DUSP5. These genes associated with proliferation and differentiation suggest that this B1 subpopulation is in an activated state. Compared with those in the control and Stage 1 groups, the number of B cells was greater in the Stage 2 group, whereas 99% of the B1-cell subpopulation was present in Stage 2. To further clarify the function of each subgroup, functional enrichment analysis was conducted via GSVA with the Gene Ontology (GO) database as a reference (Figs. 4D and S3A). The results revealed that metabolism, Th1 and Th2 cell differentiation, ribosome biogenesis and T-cell receptor signaling pathways were significantly activated in B1-cell subpopulations. Oxidative phosphorylation, thermogenesis and metabolic pathways were significantly activated in B2 cell subpopulations. Significant activation of metabolism-related pathways was observed in the B3 cell subpopulation.

Figure 4 scRNA-Seq analysis reveals B cell transcriptome in irradiated patient blood. (A) t-SNE plot visualization of B cell subclusters colored by each cluster, and between-group differences. (B) Bar plots showing the proportion of each cluster contributed by each group. Percentage of cells in Control, Stage 1 and Stage 2 for three different subsets of B cells. (C) Heatmap showing the top 10 SDE genes in each subset of B cells. (D) The function of B cell subsets was analyzed by GSVA enrichment. scRNA-Seq, single-cell RNA sequencing; t-SNE, t-distributed stochastic neighbor embedding; SDE, differentially expressed; GSVA, gene set variation analysis.

Analysis of the developmental trajectories of B cell subsets

To further analyze the evolution of the cell subsets, pseudotime sequence analysis was performed via Monocle 2 software. Fig. 5A shows the results of trajectory analysis of three relevant clusters, and Fig. 5B shows the direction of pseudotime; the black arrow signifies the initiation and progression direction of the trajectory. The pseudotime trajectory encompasses three distinct stages, with various cell clusters distinctly aligned at the corresponding points along the path (Fig. 5B). The B1 subsets were in the late stage. The B2 subsets were mostly in the early stage. The B3 subsets were mostly in the middle stage. The control group was mainly in the early stage of the pseudotime path and was mainly distributed in States 1 and 2. The Stage 1 group was mainly in the last stage of the pseudotime path and was mainly distributed in State 3. The Stage 2 group was mainly in the last stage of the pseudotime path and was mainly distributed in State 3. For the analysis of pseudotime trends among marker genes, key genes were identified within each cluster. Fig. 5C presents a heatmap indicating that these gene patterns can be divided into three modules. In module 1, which is associated with apoptosis and RNA processing, marker gene expression decreases towards the terminus of the pseudotime axis. Conversely, module 3 displays a gradual increase in marker gene expression with increasing pseudotime. Module 3 is related mainly to B-cell activation, B-cell receptor signaling activation and immune regulation. In the second module, the expression of all marker genes increased and then decreased as the pseudotime increased. Module 2 is related mainly to mRNA catabolism, the interferon-γ-mediated signaling pathway and the tumor necrosis factor-mediated signaling pathway. The data revealed dynamic gene expression profiles during different stages of B cell development. The GO and KEGG analyses revealed enriched functions and pathways of the SDE genes (Fig. S3B). These findings showed the dynamic profiling of gene expression at 631 days following radiation exposure and imply that the deregulation of multiple pathways is related to radiation exposure.

Figure 5 Analysis of developmental pseudo-time of B cell subsets. (A) Pseudo-time state visualization of B cell subclusters, colored by each cluster. Each dot represents a single cell. The black arrow indicates the start and direction of the trajectory. (B) trajectories of B cells by state, pseudotime and groups. (C) Heatmap displays the SDE genes during progression of B cells. Color key from blue to red indicates relative expression levels from low to high. The GO and KEGG analyses reveal enriched functions and pathways of the SDE genes. (D) Heat map of RAS activity of regulons in each cell population. Rows indicate different regulons and columns indicate different cell populations. The color changes from blue to red indicating a low to high RAS activity score, with higher RAS scores indicating stronger regulon activity in that cell population. Regulon specificity ranking graph. The horizontal coordinates indicate the ranking and the vertical coordinates indicate the RSS score. Higher RSS may be associated with regulon specificity in this cell population. SDE, differentially expressed; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; RAS, regulon activity score; RSS, regulon specificity score; Con, control.

Transcription factors, along with the genes they regulate, form a sophisticated network that dictates cell identity. The present study conducted SCENIC analysis to deduce the activity levels of regulons within B cell clusters, as depicted in Fig. 5D. The numbers of different B cell subpopulations changed after irradiation, but the regulon-specific sorting results shown in Fig. 5D revealed that HMGB2 was consistently the top protein. HMGB2 is a member of the high mobility group protein family. HMGB2 is associated with the activation, differentiation and function of B cells through multiple mechanisms. Although there are few direct studies, HMGB2 may play an important role in regulating B cell-mediated immune responses (44). These results suggested that HMGB2 may act as an essential driver of B cells during radiation exposure.

Analysis of T cell heterogeneity

By reanalyzing T cell subsets, the present study found that T cells could be divided into 10 clusters (T1-T10; Fig. 6A). The T cells were divided into CD4+ T cells and CD8+ T cells, as shown in Fig. 6B. CD4+ T cells are also part of the adaptive immune system, where they assist in coordinating the immune response by stimulating other immune cells, such as macrophages, B cells, and CD8+ T cells, to fight infection. Statistical analysis was conducted on the proportion of each cluster of T cells, and the results revealed that the number of cells in the Stage 1 and Stage 2 groups decreased after irradiation but did not return to normal levels. In the control group, the T1, T2, T3 and T9 subsets accounted for a large proportion. In the Stage 1 group, the T4, T6, T8 and T10 subsets accounted for a large proportion. The T5 and T7 subsets accounted for a large proportion of the Stage 2 subgroup (Fig. 6C). To further clarify the characteristics of each subgroup, the present study analyzed the marker genes of each cluster and used the top 10 genes to construct a heatmap (Fig. 6D). The results revealed that CD8B was highly expressed in the T2, T3 and T4 subsets. The T2 cell subpopulation also highly expressed GZMB. The expression of the CEBPD and NCR3 genes was high in the T7 and T8 subpopulations. MAL and LEF1 were highly expressed at T1, T3, T5 and T6. LMNA and GZMK were highly expressed in the T10 subpopulation (Fig. 6D). These preferentially expressed genes may be associated with cell differentiation.

Figure 6 scRNA-Seq analysis reveals T cell transcriptome in irradiated patient blood. (A) t-SNE plot visualization of B cell subclusters colored by each cluster, and between-group differences. (B) T cells were divided into CD4T cells and CD8T cells. (C) Bar plots showing the proportion of each cluster contributed by each group. Percentage of cells in Control, Stage 1 and Stage 2 for nine different subsets of T cells. (D) The heatmap shows the top 10 SDE genes in each subset of T cells. scRNA-Seq, single-cell RNA sequencing; t-SNE, t-distributed stochastic neighbor embedding; SDE, differentially expressed; GSVA, gene set variation analysis.

Functional enrichment analysis of each subpopulation of T cells

To further clarify the function of each subgroup, functional enrichment analysis was conducted via GSVA with the GO database as a reference (Fig. S4). The T4 subset significantly activated the hematopoietic cell lineage pathway, and the T6 subpopulation significantly activated the RNA transport and HIF-1 signaling pathways. The T8 and T10 cell subpopulations markedly activated the IL-17, MAPK and TNF pathways. The T5 subpopulation showed significant activation of DNA replication, ribosome biogenesis, metabolism and mismatch repair pathways. This subpopulation may be related to cellular repair. The T7 subpopulation of cells significantly activated cytotoxicity, the cell cycle and the T-cell receptor signaling pathway.

Analysis of the developmental trajectories of T cell subsets

For deeper insight into cell subset evolution, a pseudotime analysis was conducted with the scRNA-Seq dataset using Monocle 2. Fig. 7A-B shows the developmental stage of each subpopulation of CD4+ and CD8+ T cells. In CD4+ T cells, the pseudotime trajectory is segmented into three stages, with distinct cell clusters evident at each stage. The control group was mainly in the early state of the pseudotime path. The predominant distribution of cells in the control group at the putative time point was State 1. The T1 subset was the predominant subgroup. The Stage 1 group was mainly in the last state of the pseudotime path. The predominant distribution of cells in the Stage 1 group at the initial time point was State 3. The T6 and T10 subsets were the predominant subgroups. The Stage 2 group was mainly in the last State of the pseudotime path. The predominant distribution of cells in the Stage 2 group at the initial time point was State 3. The number of T6 and T10 subset cells increased after irradiation. The data revealed the dynamic gene expression profiles of CD4+ T cells after radiation exposure (Fig. 7C). GO and KEGG analyses revealed that post radiation exposure, immune response, apoptosis, and mRNA processing genes were significantly perturbed.

Figure 7 Analysis of developmental pseudo-time of T cell subsets. (A) Pseudo-time state visualization of CD4 T cell subclusters, colored by each cluster. Each dot represents a single cell. The black arrow indicates the start and direction of the trajectory. Trajectories of CD4 T cells by state, pseudotime and groups. (B) Pseudo-time state visualization of CD8 T cell subclusters, colored by each cluster. Each dot represents a single cell. The black arrow indicates the start and direction of the trajectory. Trajectories of CD8 T cells by state, pseudotime and groups. (C) Heatmap displays the SDE genes during progression of CD4 T cells. Color key from blue to red indicates relative expression levels from low to high. The GO and KEGG analyses reveal enriched functions and pathways of the SDE genes. (D and F) Heat map of RAS activity of regulons in each cell population. Rows indicate different regulons and columns indicate different cell populations. The color changes from blue to red indicating a low to high RAS activity score, with higher RAS scores indicating stronger regulon activity in that cell population. Regulon specificity ranking graph. The horizontal coordinates indicate the ranking and the vertical coordinates indicate the RSS score. Higher RSS may be associated with regulon specificity in this cell population. (E) Heatmap displays the SDE genes during progression of CD8 T cells. Color key from blue to red indicates relative expression levels from low to high. The GO and KEGG analyses reveal enriched functions and pathways of the SDE genes. SDE, differentially expressed; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; RAS, regulon activity score; RSS, regulon specificity score; Con, control.

The intricate gene regulatory network, which is composed of transcription factors and their regulated genes, is vital for preserving cell identity (45). SCENIC analysis was applied to ascertain the regulon activities within CD4+ T cell clusters, as illustrated in Fig. 7D. CD4+ T cells were increased in the T6 and T10 subsets after irradiation, and the regulon-specific sorting results in Fig. 7D revealed that IRF1 and RUNX2 were strongly active in both subsets. IRF1 enhances STAT1 DNA binding through the promotion of STAT1 phosphorylation, indicating a potential involvement of IRF1 in a positive feedback mechanism within JAK-STAT signaling. The RUNX transcription factor family plays a significant role in T cell development by both positively and negatively influencing the progression of T cells from double negative to double positive to single positive stages. IRF1 and RUNX2 are both crucial in T cell development and functionality, affecting T cell differentiation and immune responses through the modulation of gene expression and involvement in signaling pathways (46). These results suggested that IRF1 and RUNX2 may act as essential drivers of CD4+ T cells following irradiation.

The CD8+ T cell pseudotime trajectory is characterized by three phases, with discrete clusters delineated along its course (Fig. 7B). The predominant distributions of cells in the control group at the putative time points were State 2 and State 3. The predominant distribution of cells in the Stage 1 group at the putative time was State 1. The number of T4 and T8 subset cells increased after irradiation. The T4 and T8 subsets were the predominant subgroups. The predominant distribution of cells in the Stage 2 group at the initial time point was State 3. The data of the present study revealed dynamic gene expression profiles of CD8+ T cells following radiation exposure (Fig. 7E). GO and KEGG analyses highlighted significant disruptions in the immune response, apoptosis, and cytokine signaling pathways among genes with dynamic expression following radiation exposure.

A complex gene regulatory network, which is critical for cell identity, is formed by TFs and their downstream genes (47). SCENIC analysis revealed regulon activities across CD8+ T cell clusters, as shown in Fig. 7F. CD8+ T cells were increased in the T4 and T8 subsets after irradiation, and the regulon-specific sorting results in Fig. 7F revealed that ATF3 and JUN were strongly active in both subsets. ATF3 and JUN act as transcription factors that govern the expression of various genes through heterodimer formation (48). Their influence on AML cell response to endoplasmic reticulum stress, via UPR signaling modulation, affects cell viability. These findings suggested a potential parallel impact of ATF3 and JUN on stress response and survival in T cells. These two transcription factors are probably capable of modulating the expression of the identified genes associated with inflammation. These results suggested that ATF3 and JUN may act as essential drivers of CD8+ T cells after irradiation.

Discussion

Over the past few years, the revolution in scRNA-Seq has enabled unbiased quantification of gene expression in thousands of individual cells, which provides a more efficient tool to decipher the progression of human diseases (49). The present study collected PBMCs and analyzed their profiles via the scRNA-Seq platform (22). It subsequently isolated NK, T, and B cells from the transcriptomic data and conducted subcluster analyses to explore their properties. scRNA-Seq revealed significant heterogeneity in the profiles of these cell subsets between the irradiated cohort and healthy controls. Following radiation exposure, the host immune cells exhibited an imbalance, with reductions observed in certain T cells, neutrophils, and monocytes. As shown in Fig. S1C, the numbers of monocytes and neutrophils in the control and Stage 1 groups were notably low, whereas there was a sharp increase in the numbers of these cells in Stage 2. Conventionally, the blood half-life of human neutrophils is considered brief, ranging from 4-8 h (50). These cells are highly susceptible to apoptosis with any manipulation outside the body. Given the short lifespan of neutrophils in the bloodstream, the present study suggested that neutrophil precursors nearing maturation may be more sensitive to radiation, whereas those further along the lineage demonstrate greater radioresistance. It has been reported that monocyte chemotactic protein-1 increases in the serum of individuals exposed to radiation. Additionally, the reduction in T cell numbers after irradiation is consistent with the findings of previous studies (39).

NK cells serve as frontline defenders against infections and malignancies, initiating inflammation through the release of cytokines and chemokines (50). A notable increase in NK cell populations is typically observed following irradiation; however, the dynamics of circulating NK cells after exposure remain poorly understood. The present study revealed an increase in NK cell populations (Stage 1) and a decrease in B cells (Stage 2) among irradiated patients, as confirmed by both transcriptomic profiling and cell surface marker analysis. Using scRNA-Seq technology, a transcriptome analysis of PBMCs was conducted and revealed that T cells, B cells and NK cells may markedly contribute to the immune response following radiation exposure. Through tSNE analysis, differential gene enrichment analysis, and other comprehensive evaluations of PBMCs, it was determined that B cells and NK cells are the key subsets involved in the radiation response.

Although radiation induces a dose-dependent reduction in hematopoietic and peripheral immune cells, the specific phenotypic changes in the cellular and molecular pathways of PBMCs that persist following radiation injury are less well understood. The TRGC1 constant region of the T cell receptor (TR) γ chain plays a role in antigen recognition (51). Γ-δ TRs are capable of recognizing a diverse array of self- and foreign nonpeptide antigens, which are frequently expressed at the epithelial interfaces between the host and the external environment. This antigen recognition triggers rapid, innate-like immune responses that are crucial for pathogen clearance and tissue repair (52,53). AREG is an autocrine growth factor as well as a mitogen for astrocytes, Schwann cells and fibroblasts. It is related to epidermal growth factor (EGF) and transforming growth factor alpha (TGF-α). The protein interacts with the EGF/TGF-α receptor to promote the growth of normal epithelial cells, and it inhibits the growth of certain aggressive carcinoma cell lines (54-56). A previous study indicated that blockade of AREG signaling is a promising therapeutic strategy to mitigate radiation-induced kidney fibrosis (57).

Compared with that in nonirradiated individuals, the number of B cells in PBMCs decreased markedly after exposure to radiation. For example, NR4A3 plays a role in the regulation of the proliferation, survival and differentiation of a number of different cell types and in metabolism and inflammation. By binding to an NBRE site, SKP2 mediates the proliferation of vascular smooth muscle, myeloid progenitor cells and type B pancreatic cells and promotes mitogen-induced vascular smooth muscle cell proliferation through transactivation of the SKP2 promoter (58,59). NR4A3 induces apoptosis and inhibits cell proliferation by increasing intracellular ROS levels. NR4A3 is critical for increasing the efficacy of radiotherapy in combination with hyperthermia (60).

Delving into PBMCs provides a clearer understanding of a patient's immune status, with both diagnostic and therapeutic value. In PBMCs from the irradiated patient, a significant decrease in T cell and B cell populations was observed. It was found that the numbers of CD4+ T cells and CD8+ T cells were lower in the patient who received radiation. Lower levels of T cells suggest a role for dysregulated immune responses in radiation pathogenesis. T cells, especially CD4+ T cells and CD8+ T cells, play important roles in the immune response following radiation (7). Dose-dependent decreases in CD4+ and CD8+ T cells have been recorded (58). Ionizing radiation is known to have an immune modulatory effect on the skin by inducing a T cell response. Low-dose radiation-induced quantitative and functional alterations in immune parameters were reported by Ilienko et al (61), who reported decreased CD4+/CD8+ ratios in those exposed to low-dose radiation (62). At present, the enrichment of GO and KEGG genes differentially expressed in B cells, T cells and NK cells was analyzed in the present study, and the results indicated a number of similarities between radiation exposure processes and viral infections and that the underlying mechanism still needs to be further studied.

Human blood lymphocytes are crucial components of the immune system and are markedly affected by radiation exposure. Among cell groups, lymphocytes are particularly sensitive to radiation and serve as important markers for the early detection of radiation injury and assessment of exposure levels. Radiation can alter the radiosensitivity of different lymphocyte subsets, influencing immune response dynamics. In general, B cells are more susceptible to radiation than are T cells, whereas NK cells demonstrate greater resistance (63). These findings are consistent with the present study. Individuals such as nuclear accident victims, those exposed to radiation over an extended period, or patients receiving radiotherapy for tumors often experience severe immune deficiencies, with changes in lymphocytes playing a key role in radiation-induced immune damage.

The present study had several limitations, including a small sample size for scRNA-Seq analysis, which may have affected the statistical power for differential abundance and expression analyses. The present study focused on NK, T and B cells without detailing the transcriptome features of other members (for example, bone marrow mononuclear cells, neutrophils and monocytes), as well as changes in inflammatory cytokines. Owing to the specificity of the samples, the present study lacked flow cytometry results to verify the clustering of peripheral blood cells together with single-cell sequencing results. Furthermore, the present study focused on PBMCs and may not fully represent localized immune responses. For improved understanding of the relationship between immune responses of different cell types and radiation, further research involving a larger and more diverse patient population with various clinical presentations is necessary.

In conclusion, the present study demonstrated the differences in several cell subsets between an irradiated patient and healthy controls via scRNA-Seq technology. Differential cell subpopulations in PBMCs and preferentially expressed genes in NK cells and B cells of the irradiated patient were both identified, which may offer new insights into disease diagnosis and therapeutic intervention. Limitations exist because longitudinal investigations to further determine the specific roles of the obtained genes in the development of irradiation are lacking, which will be addressed in future clinical and basic research.

Supplementary Material

Verification of the differences in subsets of PBMCs in irradiated patients. (A) Heat map of differentially expressed genes in samples from irradiated patient and normal controls. (B) KEGG analysis of pathways involved in upregulated genes and downregulated genes. (C) The t.SNE map shows the different distribution of 6 clusters in Con, Stage 1, and Stage 2. The bars on the right show the number of cells in Con, Stage 1 and Stage 2 for the five cell types. PBMCs, peripheral blood mononuclear cells; KEGG, Kyoto Encyclopedia of Genes and Genomes; t.SNE, t.distributed stochastic neighbor embedding; Con, control.

The potential function of NK cell subsets. The abscissa represents different cluster information, and the ordinate is the pathway name. The color from blue to red indicates the higher enrichment of the pathway in this cluster. The cell subsets were analyzed by GSVA enrichment. NK, natural killer; GSVA, gene set variation analysis.

KEGG analysis of differential gene expression. (A) KEGG analysis of differential gene expression in NK cells. (B) KEGG analysis of differential gene expression in B cells. KEGG, Kyoto Encyclopedia of Genes and Genomes; NK, natural killer.

The potential function of T cell subsets. The abscissa represents different cluster information, and the ordinate is the pathway name. The color from blue to red indicates the higher enrichment of the pathway in this cluster. The cell subsets were analyzed by GSVA enrichment. GSVA, gene set variation analysis.

Characteristics of the patient and healthy controls.

Blood routine examination at Stage 1.

Blood routine examination in Stage 2.

Acknowledgements

Not applicable.

Funding

Funding: The present study was supported by the National Natural Science Foundation of China (grant nos. 82473574 and 82373523), Science and Technology Project of Sichuan Province (grant nos. 20YYJC035 2023NSFSC0648 and 2024ZYD0126), Chengdu Innovation Project (grant no. 2021-YF05-01603-SN), NHC Key Laboratory of Nuclear Technology Medical Transformation (Mianyang Central Hospital; grant no. 2023HYX022) and Natural Science Project of Chengdu Medical College (grant no. CYZYB23-02).

Availability of data and materials

The data generated in the present study may be found in the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) and are accessible through GEO Series under accession numbers GSE166902 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE166902) and GSE190439 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE190439).

Authors' contributions

DY, HZ and SZ designed and supervised the present study. ZJ, WT and KF collected the clinical samples and clinical measurement information. ZJ and WH analyzed the clinical data. TY, WT and XX conducted scRNA-Seq analyses. WH, JC, ZJ, SZ and DY interpreted the data. TY, WT and ZJ wrote the manuscript. DY, SZ, HZ, JC, XX and WH edited the manuscript. SZ, WT and DY confirm the authenticity of all the raw data. All authors read and approved the final manuscript.

Ethics approval and consent to participate

The present study was approved by the Ethics Committee of the Second Affiliated Hospital of Chengdu Medical College, Nuclear Industry 416 Hospital [approval number: 2019(25)]. Written informed consent was obtained from all participants.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References