Spleen tissues from HS patients with hepatitis B virus (HS-HBV, PH024), portal hypertension (HS-PHT, PH026), hepatitis C virus (HS-HCV, PH027), and hepatocellular carcinoma (HS-HCC, PH037), and one patient with a normal spleen (PN021) were used ( Table 1 ). The average age of patients was 54.4 years (42-67 years), and two were male and three were female. All patients with HS showed splenomegaly on imaging. HS patients demonstrated a decrease in blood cells (single or simultaneous reduction of red blood cells, white blood cells, or platelets). After splenectomy, blood cells returned to a near-normal level.

All samples were digested to a single-cell suspension and studied by 10× Genomics scRNA-seq techniques ( Figure 1A ). After quality control, a total of 43,037 cells were obtained from all samples, of which 24,387 cells (56.7%) were derived from HS patients, about 6,096 cells per HS sample, and 18,650 cells (43.3%) from the normal control sample. The number of cells of control sample was about three times that of the HS sample, and the PH024 and PH037 samples had the least cells. Using t-distributed stochastic neighbor embedding (t-SNE) analysis, all cells were classified in 24 clusters ( Figures S1A–C ), and they were further classified to 9 major cell clusters based on marker genes, including T-cells (61.64%), B-cells (26.72), NK cells (3.68%), hematopoietic stem cells (HSC, 2.47%), neutrophil cells (2.44%), mast cells (0.99%), endothelial cells (0.97%), erythrocyte (0.58%), and dendritic cells (DCs, 0.51%) ( Figures 1B , S1D ). Furthermore, the clusters displayed in t-SNE were labeled according to sample source and disease state ( Figures 1C, D ; Figures S1E, F ). Furthermore, those clusters were assigned to known cell lines through marker genes ( Figure 1E , Figures S1G, H ). Compared with healthy control spleen tissue, there was a deficiency of neutrophils in HS, an increase in dendritic cells, erythrocyte, HSC, mast cells and NK cells ( Figures S1E, F ).

In HS samples, T-cells and B-cells were the main cell subsets (accounting for 86.8%), and significant differences in T-cells and B-cells were found among different types of HS ( Figure 1F ). T-ells accounted for the highest proportion in PH037, and other HS samples had the same proportion of T-cells as the control sample ( Figure S1D ). The proportion of B-cells in PH027 was similar to that of the control sample, followed by PH024 and PH026 samples, while the lowest proportion was found in PH037 at number less than half that of the control. For NK cells, only PH027 was lower than the control sample, and other HS samples contained more than twice that of control sample. The number of HSC in PH037 was three times that of the control, while in PH027 it was only half that of the control. Other cell types with lower proportion also showed great differences among normal spleen and HS, such as dendritic cells (DCs) and mast cells (extremely rare in PH027) ( Figure 1F , Figure S1D ). Overall, T-cells and B-cells are the core cells in the spleen, and the broad changes of transcriptional profiles and cell proportions in HS were characterized, which may be the main cause of immune disorder.

Next, a total of 26,528 T-cells were detected in all samples, and were classified to 29 clusters through t-SNE ( Figures S2A–C ). According to the CD4 and CD8 makers, 12 clusters belonged to CD4⁺ T-cells (44.33%), and 17 clusters belonged to CD8⁺ T-cells (55.67%) ( Figure 2A ). Among all the HS, more CD4⁺ cells than CD8⁺ cells were found in PH307, while in PH027, CD8⁺ cells were more numerous than CD4⁺ cells; Other samples were the same as that in the healthy control ( Figures 2B–D ). In comparison, HBV-infected spleens had more CD4⁺ T-cells; HCV-infected spleens had more CD8⁺ T-cells ( Figures 2B–D ). Then, T-cells were further marked by other genes and clustered to seven groups, including CD4⁺ Naive cells (21.56%, SELL), CD4⁺ T helper cells (20.39%, GZMA), CD4⁺ Tregs cells (2.38%, FOCP3), CD8⁺ T helper cells (31.06%, SELL), CD8⁺ Naïve cells (12.21%, SELL), CD8⁺ Cytotoxic cells (9.56%, GNL1), CD8⁺ MAIT cells (2.84%, SLC4A10) ( Figures 2E, I ; Figure S2D ). Of these seven T-cell subpopulations, CD4⁺ Naïve cell, CD4⁺ T helper cell, and CD8⁺ T helper cell groups were the core cell groups, accounting for 73.01% of all T-cells, and they showed different abundances in different HS samples ( Figures 2F–H ). Only the number of CD4⁺ Naïve cells of PH027 was equivalent to that of the control, and that of other HS samples were about twice as high as the control. The content of CD4⁺ T helper cells in PH037 was higher than that in the control sample, while that in other HS samples was lower than in the control, and it was the lowest in PH026. The content of CD8⁺ T helper cells in PH027 sample was higher than that in control sample, while it in other HS samples was lower than control, and it was the lowest in PH037 ( Figures 2F, G ). T-cells were labeled according to the source of the sample, and HS and normal spleen were roughly separated using t-SNE ( Figure 2H ).

To study T-cell subclusters in different HS types, a functional enrichment analysis was conducted on the enriched genes and the specific marker genes for each T-cell subclusters were screened out. Among these T-cells, some subclusters were specifically present in the sample. For example, one population of CD4⁺ Naïve cells and one population of CD8⁺ Cytotoxic cells were specifically expressed in PH024 ( Figure 2J ). The GSEA enrichment results, which were analyzed based on the enriched genes of specific expressed T-cell subclusters, indicated that the specific T cell in PH037 had stronger red blood cell differentiation signals, and the specific T-cell subclusters in PH024 sample had stronger PD-1 signals (25, 26), Cystic Fibrosis (27), and Polarized epithelium signals (28), and the specific T-cells in PH027 sample had stronger positive regulation signals for natural killer cell-mediated immunity. PH024 and PH027 were simultaneously enriched in cell adhesion mediated by integrin ( Figure 2K ). Among those high-expressed genes in T-cells, we found none were highly expressed in the healthy sample ( Figure 2L ). However, there were 55 highly-expressed genes in HS, and half of them were highly expressed in all HS samples. For instance, RASA3, ADAR, XAF1, and IKZF1 are generally highly expressed in T-cells in HS and are associated with immunosuppression ( Figure 2L ). In particular, ADAR was found to be a target for cancer immunotherapy (29). Furthermore, there were some genes that were only highly expressed in one HS. TNFSF8 was highly expressed in PH037 alone and was also enriched in CD4⁺ Naïve T-cells and CD4⁺ Treg cells ( Figures 2L, M ). Furthermore, previous research have found that TNFSF8 was associated with regulation of immune response (30) and high risk of hepatitis infection (31). In general, our research found that abnormal distribution and proportion of T-cells can be a characteristic of HS; these T-cell enriched functions and genes have certain characteristics, which are mainly manifested through immunosuppression.

A total of 11,500 B-cells were classified into 13 clusters ( Figures S3A–C ), and these cells were further categorized into nine subclusters according to marker genes ( Figures 3A–C, H ). Firstly, HS was found to be significantly deficient in CCL5 ⁺ B-cells, while many SESN1 ⁺ B-cells and SNAI1 ⁺ B-cells are enriched ( Figures 3B, F, G ). In addition, the distribution of B cell subclusters in different HS varies. Among B-cell subclusters, CCL5 ⁺ B-cells were found to be specifically present only in healthy control, while SNAI1 ⁺ B cells in PH027, and SESN1 ⁺ B cells in PH037. And the HLA-DRB5 ⁺ B-cells were highly expressed in PH026 and PH027 ( Figures 3C–E ).

The GSEA enrichment analysis were performed based on the enriched genes in specific expressed B cell subclusters and showed that the specific B-cell subclusters in the healthy sample had stronger cell chemotaxis signals ( Figures 3I, J ). The specific B-cell subclusters in PH024 had stronger PD-1 (25) and polarized epithelium signals (28). The specific B-cell subclusters in PH026 had stronger mRNA splicing regulatory signals, and the specific B-cell subclusters in PH027 had stronger heat shock response signals, and the specific B-cell subclusters in PH037 had stronger cytokine responses and humoral response signals ( Figure 3J ).

Among those highly expressed genes screened in B cells, only POLR2L was highly expressed in the healthy sample ( Figure 3K ), which was reported as a part of the core element for RNA polymerases (32). Among the B cells, POLR2L was mainly expressed in B cell-1, CCL5 ⁺ B cells and plasma cells ( Figures 3K, I ). There were 35 genes that were only highly expressed in HS samples, and most of them were highly expressed in more than one HS samples, but some genes were only highly expressed in one HS sample. BTK, SYK, and HLA-DRB5 were only highly expressed in PH024, and BTK is a small-molecule inhibitor which exhibited impressive anti-tumor activity in clinical studies in patients with various B-cell malignancies (33). NMT1 and RERE were only highly expressed in PH026, and NMT1 was studied widely for breast cancer and showed suppression on initiation, proliferation and invasion of breast cancer cells (34). CD70 was only highly expressed in PH027, and it was implicated in tumor cell and regulatory T-cell survival through interaction with its ligand (35). PIK3IP1, RNASE6, and SESN1 were only highly expressed in PH037, and those genes also highly expressed in SESN1 ⁺ B-cells ( Figures 3K, L ). While SESN1 has been studied and found to be associated with p53 mutations (36). Collectively, the abnormality of B cells could be more clearly observed in HS; furthermore, we found the lack of CCL5 ⁺ B cells and the abnormal presence of HLA-DRB5 ⁺, SESN1 ⁺ and SNAI1 ⁺ B cells in HS were related to the etiology of HS.

To detect the evolutionary trajectory of splenic cells, Monocle 2 was used to perform a pseudo-time analysis for all cells, T-cells, and B-cells, respectively. The results of cell trajectory analysis on all cells showed that no difference existed among different HS samples ( Figure 4A ). HSC and endothelial cells were seen to be precursors to other cells ( Figure S4A ). The pseudo-time analysis based on T-cells and B-cells demonstrated different clear evolutionary trajectories. In the evolution trajectory of T-cells, these cells were distributed on the trajectory chart composed of three branches (progressing from Branch 1 to Branch 3, Figure 4B ). The developmental trajectories of PH024 and PH037 are relatively consistent. Compared with the control PN021, PH026 T-cells were overly naïve, and PH027 T-cells were senescent ( Figure 4B ). The composition of the corresponding cell subsets developed from CD4⁺ Naïve and CD8⁺ Naïve to CD4⁺ T helper and CD8⁺ T helper, and finally to CD8⁺ MAIT and CD8⁺ cytotoxic cells ( Figure 4C ). All B-cells were distributed on a trajectory diagram composed of four branch points, starting from Branch 2 to Branch 4. B-cells from healthy spleen were evenly distributed on the developmental trajectory, while the trajectory of B cells in HS showed a lack of development ( Figure 4D ). Analysis of the developmental trajectory of B-cell subgroups shows that Bcell-2, IGKV1-39 ⁺ B-cells, and SESN1 ⁺ B-cells were in the early developmental trajectory, while CCL5 ⁺_B-cells and Plasma_B-cell were in a mature functional state ( Figure 4E ). In summary, the developmental trajectory of all cells in HS was not different from that of the control spleen; but the T-cells and B-cells in HS were in a state of incomplete development.

Next, CellPhoneDB (37) was used to determine the ligand-receptor interactions among HS cells, showing that there were significant differences in ligand-receptor interactions between HS and normal control. The cell group having strongest interaction with the other cell groups involved neutrophil cells in normal spleen, but the HSC and endothelial cells were the strongest interaction cell groups in HS, especially in PH037 ( Figure 5A ). This finding was consistent with the phenomenon observed by Wang and Katztong such that endothelial cells directly contact HSCs (38). For the ligand-receptor analysis, there were differences among HS and normal spleen ( Figure 5B ). T cells were mainly regulated by B-cells, DCs, and neutrophil cells through CD74-MIF. Compared with the normal spleen, NK cells exerted a negative regulatory effect on T-cells through CD94-HLA-E in HS samples ( Figure 5B ). B-cells in the spleen play a major regulatory role, whether in HS or control. B-cells mainly interacted with other cells through CD74. In addition, B cells in HS exert regulatory effects on a variety of cells through CD74-COPA. Importantly, it was found that both T cells and B cells in PH were negatively regulated by NK cells through CD94-HLA-E (39) ( Figure 5B ). In summary, abnormal enhancement of the interaction between epithelial cells and HSC in HS, and T/B-cells in HS was found to play a regulatory role through CD94-HLA-E and NK cells.

Except for T-cells and B-cells, other cells accounted for 13% of the total cell count, but there were also specific marker genes enriched in these clusters. A total of 184 specific marker genes were screened from these cells. RPS4Y1 was highly-expressed in a variety of cells in the normal spleen ( Figure 6 ). Through database analysis, we found RPS4Y1, a member of the S4E family of ribosomal proteins (40), was poorly expressed in a variety of tumors including LIHC ( Figures S5A, B ). The high expression of RPS4Y1 was associated with high fibrosis and high immune infiltration in the liver ( Figure S5C, D ), and the expression of RPS4Y1 was positively correlated with Tcm, Tem and pDC ( Figure S5F ). High expression of PRS4Y1 indicated a better prognosis ( Figure S5E ). Furthermore, FKBP5 was found to be specifically expressed in a variety of cells in PH037 ( Figure 6 ), and it was also enriched in Treg cells in HS ( Figure 2M ). Previous research indicated that FKBP5 was related to immune response, and the high expression of FKBP5 was linked to a pro-inflammatory profile and altered NF-κB-related gene networks (41). In our study, FKBP5 is highly expressed in a variety of tumor tissues ( Figure 7A ). In liver cancer, FKBP5 was related to tissue type and immune infiltration, but not to the degree of fibrosis ( Figures 7B–D ). FKBP5 is negatively correlated with B cells, DCs and T cells in the microenvironment ( Figure 7E ). In the analysis of immune cells and immune factors, FKBP5 specifically and negatively correlated with a variety of immune effector cells and immune effector factors in the microenvironment of liver cancer ( Figures 7F–I ). In addition, HLA-DRB5 was enriched in HSC and neutrophils of PH024 and PH026; XIST was enriched and expressed in PH026 and PH037. Collectively, only a small number of genes in HS were common in different HS samples, while the expression of RPS4Y1 is generally absent in HS, and FKBP5 is enriched and expressed in HCC-HS.

In this study, four patients with HS (group: PH) and a patient with normal spleen (group: PN) were selected. The splenomegaly of PH group patients was caused by different reasons, including hepatitis B virus infection (HS-HBV, PH024), portal hypertension (HS-PHT, PH026), hepatitis C virus infection (HS-HCV, PH027), and hepatocellular carcinoma (HS-HCC, PH037). The patient with a normal spleen was used as control (PN027), whose spleen was ruptured due to trauma.

This study was approved by the Ethics Committee of Shanxi Bethune Hospital, and complied with all relevant ethical regulations (SBQLL-2020-038). Each participants signed an informed consent proforma.

Samples were freshly collected from patients during surgical resection. Firstly, the spleen in the abdominal cavity was removed from its blood supply and quickly excised by surgical instruments. To maintain cell viability, each specimen was immediately immersed in culture medium and stored in separate 50 ml centrifuge tubes with 10% fetal bovine serum (FBS). Then, about 0.1-1 g spleen tissue was transferred into a 100-mm Petri dish and was cut into 1-2 mm tissue blocks using sterilized scissors and tweezers. The tissue blocks were transferred to a 10 ml centrifuge tube and were washed two times by adding PBS solution, centrifuging at 300g for 5 minutes, and removing the supernatant. Thereafter, the tissue blocks were transferred to cryopreservation tube and 1 ml of cryopreservation fluid was added (CELLBANKERTM2, ZENOAQ). After standing at room temperature for ten minutes, the samples were stored at -80 °C for further investigation. All laboratory work in this study was conducted in a biological safety cabinet to prevent contamination.

To maintain cell structure, spleen tissue was dissociated into single-cell suspensions by enzymatic degradation of the extracellular matrix. Firstly, the tissue was washed by using ice-cold PBS buffer and soaked in 30 mM EDTA solution for 20 mins on ice to inhibit nuclease activity and reduce the stability of the cell membrane. The tissue was then washed again using PBS buffer, and 2 mg/ml collagenase IV (Worthington Biologicals) was used to dissociate tissue into single cells at 37°C for 20 minutes. To maintain cell viability, the single-cell suspension was put at ice and 10% FBS solution added. After dissociation, the suspension was filtered through a 40-um cell strainer (BD Falcon, BD Biosciences, San Jose, CA, USA) to remove large particles.

Single-cell barcoding and scRNA-seq library preparation were performed based on 10X Genomics single-cell RNA sequencing platform (10×Genomics, Pleasanton, CA, USA). Before barcoding, the cell concentrations of single-cell suspensions were counted using a hemocytometer (TC20, Bio-Rad, Hercules, CA, USA) and adjusted to 1,000 cells/μl. Single-cell suspensions were loaded on the Chromium Single Cell Controller Instrument (10× Genomics, Pleasanton, CA, USA) to generate single cell gel beads in emulsions (GEMs), and 10×Genomics Chromium barcoding system was adopted to construct 10×barcoded cDNA library according to the manufacturer’s instructions. All libraries were sequenced on the Illumina HiSeq X Ten sequencing platform (Illumina, San Diego, CA, USA) with 150 bp pair-end module.

The CellRanger ™ pipeline (v3.0.0) was used to pre-process the raw sequencing data, including demultiplexing of cellular barcodes, evaluation of the total number of cells, and applying map reads to the genome. Briefly, after demultiplex cellular barcodes, the total number of cells captured in each sample were evaluated based on CellRanger ‘cell ranger count’ embedded function. The generated data were then mapped to the GRCh38.p5 human reference genome using the STAR aligner (https://github.com/alexdobin/STAR). Finally, the gene expression matrices of each cell were obtained for downstream visualization. After use of the CellRanger™ pipeline, a filter criterion was used to obtain high-quality cells. To remove multiple captured cells, the unique molecular identifier (UMI) count matrix (64) was calculated using R package Seurat (v2.3.4) (https://satijalab.org/seurat/), and the cells with UMI from 500 to 10,000 were selected. The double T-cells generated from two or more cells, were also removed using ‘DoubletFinder’ in Seurat. For each sample, we filtered the genes detected in fewer than three cells and cells with the total number of detected genes (nGenes) less than 500 and more than 3000. Following visual inspection of the distribution of cells by the fraction of mitochondrial genes expressed, low-quality cells where more than 25% of the counts belonged to mitochondrial genes, were further discarded.

The remaining cells and genes were used for cell clustering, and the normalization of gene expression was performed on the filtered matrix using by ‘NormalizeData’ function in Seurat. The total expression of each cell was 10,000. To compare transcriptome profiles along five different samples, we merged five different datasets using ‘merge’ of Seurat (65). The FindVariableFeatures tool was used to identify variable genes whose expression was between 0.125-3 times the average level of gene expression. All variable genes were ranked in descending order and the top 2,000 variable genes selected. The Principal component analysis (PCA) on cells was performed on cells on the basis of the expression of variable genes using ‘RunPCA’ in Seurat. Based on the PCA results, the top 30 principal components were selected for clustering analysis based on t-distributed stochastic neighbor embedding analysis, and the FindCluster tool was used to cluster the cells.

To identify cell subtypes, the cell clusters were annotated with the previously reported maker genes. The specifically expressed genes of each cell subtype were studied using the ‘FindMarkers’ function in Seurat, and the differential expressed genes were screened with the criterion (qvalue: 0.01, N0: 20, N1: 1.5, N2: 30%, N3: 3) and divided into three groups: significantly enriched genes (qvalue < 0.01), strictly enriched genes (qvalue < 0.01), and specific molecular marker genes of all cell groups. To sub-group cell clusters of interest for a more in-depth analysis, the ‘SubsetData’ function in Seurat was used to extract each cluster of interest. The extracted subclusters were also annotated with the previously reported maker genes to provide higher resolution for dissecting cellular heterogeneity among particular types of cells. Based on the significantly enriched genes for all samples, gene enrichment was conducted on the specific subclusters of each sample, and the GO and KEGG Pathways for the enriched gene were enriched using the ‘clusterProfiler’ software package.

Based on the significantly enriched genes for each subcluster, cell differentiation trajectories were analyzed. Monocle (v 2.10.0) software (66) was used to order single cells along pseudo-time axis. The pseudo-time analysis tool Monocle (v 2.10.0) can simulate the dynamic changes of time development using machine learning methods (67). First, the significantly enriched genes were selected and performed spatial dimensionality reduction based on their expression. Then, the minimum spanning tree was constructed to determine the differentiation trajectory of the cell. Finally, the continuous evolution process of cells was revealed, and the position and relationship of different cell types in the evolutionary process were also distributed.

Cell-cell interactions between different cell types were predicted based on known ligand–receptor pairs by Cellphone DB v2.1.0 (37). Permutation number for calculating the null distribution of average ligand-receptor pair expression in randomized cell identities was set to 1000. Individual ligand or receptor expression was subject to a threshold based on a cut-off based on the average logarithm of the gene expression distribution for all genes across each cell type. Predicted interaction pairs with a p-value below 0.05 and of average log expression exceeding 0.1 were deemed significant. The interaction between cells was visualized using Cytoscape (68).