Transcriptomic Profile and Expression of Immune Checkpoints in Peripheral Blood Mononuclear Cells of Women with Breast Fibroadenoma

Introduction

Fibroadenoma of the breast (FA) is the most frequently diagnosed type of benign (non-cancerous) fibroepithelial tumors characterized by a biphasic structure demonstrating proliferative changes in both epithelial and stromal (connective tissue) components in varying proportions. The main histological variants of FA are well known, and in most cases this diagnosis implies conservative treatment or, depending on specific medical indications, surgical removal of the tumor [1]. However, upon closer examination, fibroepithelial tumors represent a heterogeneous group of diseases with highly variable and largely overlapping morphological features, which can complicate differential diagnosis and treatment decisions, and may increase the likelihood of overestimation or, conversely, underestimation of disease severity and prognosis, as well as the risk of developing a malignancy [2, 3]. Besides FA, a variety of histopathological changes have been described, affecting not only stromal but also epithelial compartments, sometimes including features of atypical hyperplasia and carcinoma in situ [4].

Clearly, molecular genetic causes must underlie this morphological and clinical variability rooted in dysregulated signaling mechanisms. However, the mechanisms underlying the pathogenesis of FA remain largely unknown [5] due to a rather paradoxical situation: FA is rarely chosen as the primary and independent object of molecular genetic research due to its favorable clinical prognosis in most cases; women with this diagnosis usually constitute a comparison group (or control group) used to study various aspects of the molecular pathogenesis of malignant breast tumors, and therefore FA is considered a normal variant. Changes in immune status in patients with FA are also poorly understood due to the prevailing opinion that fibroepithelial diseases have a (dys)hormonal etiology. Very little is known about possible changes in the immune microenvironment or the profile of circulating immune-related biomarkers in FA. However, based on current knowledge of bidirectional immunoendocrine interactions (including the mutual influence of the estrogen and progesterone, and also between lymphocytes, monocytes and granulocytes), it is possible to hypothesize the involvement of the immune system in the development of FA.

Some experimental data directly or indirectly indicate that immunoregulatory mechanisms may be involved in the pathogenesis of FA. First, despite its benign nature, the development of FA is probably associated with mutational processes. Large-scale genomic screening of fibroepithelial tumors (fibroadenomas and phyllodes tumors) has recently revealed a specific landscape of recurrent genetic alterations in key proto-oncogenes and tumor suppressor genes [6], with MED12, KMT2D, and RARA being the most frequently mutated in ordinary FAs. Next-generation sequencing detected mutations in genes encoding various components of the DNA repair system [7, 8], implying the emergence of neoantigens and, consequently, the existence of immune surveillance of these antigens, including the so-called immune checkpoints [9, 10]. Importantly, this mutational landscape was linked with the tumor growth pattern and some other morphological characteristics of FAs [1]. With regard to the genome-wide methylation profile, fibroepithelial tumors (FAs and phyllodes tumors) form a distinct cluster that differs from both normal breast tissue and cancer tissue. As for the pattern of copy number variations (CNVs), abnormalities commonly observed in malignant pathology cannot be ruled out [11]. Clearly, the structural genomic and epigenomic aberrations observed in FA may manifest themselves in alterations of the activation of various signaling pathways, including those mediating relationships with immune mechanisms [7]. Various changes in the secreted proteome were identified in FA samples: specifically, the content of FA-derived extracellular vesicles differed much from that of normal tissue, but was largely consistent with breast cancer in terms of signaling pathway profile [12], contributing to the formation of a specific local immune environment in FA. Indeed, very few published transcriptomic studies focusing on the issue of benign tumors have revealed the involvement of immune-related signaling pathways [13, 14]. Some studies reported increased expression of some specific markers of epithelial-mesenchymal transition (EMT) and extracellular matrix remodeling known for their immunomodulatory activity [15]. Apart from the molecular genetics data, there is some direct evidence demonstrating local changes in the distribution of different populations of immune cells in FA, which most likely reflects dysregulations in the immunosurveillance. For instance, some studies described the changes in the levels of CD4/CD8 Т lymphocytes and antigen-presenting cells infiltrating FA vs. the normal breast tissue [16, 17]. Terzeman et al. showed that the connective tissue compartment in FA was infiltrated by mast cells producing vasoactive mediators and proinflammatory immune factors [18].

Local molecular genetic alterations can extend to the systemic level and can be detected in the peripheral circulation as abnormalities in the phenotypic characteristics of peripheral leukocytes or immunological parameters of blood plasma. For example, proteomic analysis of plasma samples from patients with FA or other benign breast lesions revealed systemic alterations associated with innate immune signaling pathways [19]. In another study, whole blood transcriptome analysis allowed the characterization of specific gene expression profiles encompassing adaptive immune mechanisms in a number of high-risk breast pathologies, including FA [20]. There are also data demonstrating altered levels of some circulating microRNAs known to be involved in immune regulation and inflammation mechanisms [21]. Regarding the variations in immunophenotypes of different lymphocyte populations, there is increasing attention paid to the expression of immune checkpoint markers and the possibilities of their clinical application in various breast neoplasms, since these markers are believed to reflect not only the state of cellular immunity activation but also its exhaustion (immunosuppression), which may affect the prognosis of the disease and the effectiveness of a particular therapeutic approach. In FA, a marked increase in the expression of the PD-1 checkpoint marker in the circulating regulatory T cell (Treg) population and a trend towards increased PD-L1 expression in regulatory B cells without changes in cell number was demonstrated [22]; however, the expression/co-expression patterns of other immune checkpoint molecules have not yet been investigated. Similarly, there are no published studies examining transcriptomic changes in various peripheral blood cell subsets, primarily lymphocytes/monocytes, in patients with FA.

The aim of this study was to analyze putative changes in cellular immunity in women diagnosed with FA in two aspects: (a) changes affecting the transcriptomic profile of the peripheral blood mononuclear cell (PBMC) fraction and the spectrum of activated/dysregulated signaling pathways that form a specific immune and stromal environment of the tumor; (b) changes in (co-)expression of immune checkpoint markers – CD279/PD-1, CD274/PD-L1, CD366/TIM3 and CD223/LAG3 – in the CD4 and CD8 subsets of circulating T lymphocytes. Data were collected using RNA sequencing technology, followed by bioinformatic analysis and flow cytometry.

Material and Methods

Patients and samples

Peripheral blood samples were collected from women diagnosed with breast FA (23-54 years of age, mean age 35.8 years, n=16) undergoing surgery at the private medical institution, RZhD-Medicine Clinical Hospital (Petrozavodsk, Russia), as well as from healthy women of the control group (22-52 years of age, mean age 38.6 years; n=16; no breast pathology and no acute infectious, inflammatory or chronic diseases). Venous blood was collected immediately before surgical removal of the breast tumor; samples were transported to the laboratory at a temperature of +4 ℃ and processed no later than 2 h after blood collection. In all cases, the primary diagnosis was confirmed by histopathological examination. Patient exclusion criteria were as follows: signs of a malignant neoplasm of the breast, a history of oncological diseases, pregnancy or lactation. All women were informed about the purposes of blood collection and provided written voluntary consent. The study was approved by the Medical Ethics Committee of the Institute of Medicine, Petrozavodsk State University and the Republic of Karelia Ministry of Healthcare (Protocol No. 3, approval date: October 21, 2024) and was conducted in accordance with the Declaration of Helsinki and good clinical practice guidelines.

RNA sequencing (RNA isolation and cDNA library creation).

PBMCs were isolated by density gradient centrifugation in Ficoll-Paque PLUS density gradient medium (density: 1.077 g/mL), GE Healthcare, USA. The total RNA fraction was purified from PBMCs by phenol-chloroform extraction using TriZOL reagent (Invitrogen, USA) according to the manufacturer’s instructions. RNA samples were then treated with DNase I (New England Biolabs, USA) to remove contaminating genomic DNA, followed by the spin column-based purification of the DNase reaction products with the CleanRNA Standard kit (Eurogen, Russia). Total RNA concentration was determined fluorometrically using a selective RNA dye and the Equalbit RNA BR Assay kit (Vazyme, China). After DNase treatment and purification, the integrity of the total RNA samples was monitored in a Fragment Analyzer automated capillary electrophoresis system (Advanced Analytical Technologies, Agilent, USA) using the HS RNA reagent kit (Agilent Technologies, USA). RNA quality assessment was based on the 28S:18S rRNA ratio. The presence of residual contaminants in RNA samples was determined based on the absorbance ratio at 260/280 nm and 260/230 nm using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, USA). Solely RNA samples that met quality requirements were selected for subsequent cDNA library preparation and sequencing.

cDNA libraries were created using the TruSeq Stranded Total RNA with RiboZero single index reagent kit (Illumina, USA) supplemented with SuperScript III reverse transcriptase (Invitrogen, USA) and Agencourt RNAClean XP or Agencourt AMPure XP magnetic beads for purification of RNA fragments or cDNA and PCR products, respectively (Beckman Coulter, USA). The amount of RNA used for cDNA library preparation was 500-1,000 ng. Quantification of the synthesized cDNA libraries was performed by fluorometric analysis based on an intercalating dye using the QuDye BR reagent kit (Lumiprobe, Russia). Correspondence of the obtained DNA fragments to the theoretically expected size of 260 bp was checked using the Fragment Analyzer system and the HS NGS reagent kit (Agilent Technologies, USA). After normalization, denaturation, and dilution steps, the amplified cDNA libraries were loaded into a flow cell; paired-end reads of 75 bp in length were obtained using the MiSeq v3 reagent kit (150 cycles) on the MiSeq sequencing platform (Illumina, USA).

Bioinformatic analysis

Bioinformatic analysis was performed on the Galaxy web portal (https://usegalaxy.eu) [23]. The quality of the raw sequence data was assessed using the FastQC quality control tool. The reads were then aligned to the human reference genome (assembly version: GRCh38/p13), and BAM files were generated using RNA STAR mapper and the GTF genome annotation file downloaded from the GENCODE server (https://www.gencodegenes.org/human). A summary of the read alignment statistics is presented in Table 1. In addition, the distribution of mapped reads across the genome and its main structural elements, viz., exons (including 5'- and 3'-untranslated regions), intronic and intergenic regions, was examined using the Read Distribution tool from the RSeQC package based on the resulting BAM files. As expected, the majority of sequenced reads (%) was mapped to exonic regions, while the proportion of reads aligned to non-coding intergenic regions indicated the amount of residual contamination of genomic DNA.

Table 1. Summary statistics of RNA STAR read alignment results

Quantification of relative transcript abundance (gene counts) was performed using the HTSeq package. DESeq2 was then applied to the gene count matrix to identify genes differentially expressed between PBMCs from FA patients and healthy controls, and to obtain normalized expression values ​​for the entire panel of samples. Genes that exhibited |log₂FC| values ​​>1.0, where FC is the fold change in gene expression and the adjusted p-value (adj p-val) < the threshold, were considered significant differentially expressed genes (DEGs); the choice of significance threshold (0.05 or 0.25) is specified in the Results section. One of the DESeq2 outputs (a PCA scatterplot illustrating the results of a bivariate principal component analysis) was used to assess the degree of similarity/dissimilarity in global gene expression profiles between samples, as well as to identify datasets for comparison and identify DEGs or differentially activated signaling pathways between them. To better visualize differences in expression patterns between the compared sample groups, heatmaps were constructed using unsupervised hierarchical clustering based on Z-score normalized gene expression values, as well as Volcano plots using Galaxy’s built-in tools (heatmap2). Ensembl gene IDs were similarly converted to corresponding gene symbols using Galaxy’s built-in function for annotating DESeq2 output tables.

Enrichment analysis of DEGs by the Gene Ontology (GO) terms categorized into three categories – Biological Process (BP), Molecular Function (MF), and Cellular Component (CC) – as well as by KEGG terms was performed in the GOSeq package using the Benjamini-Hochberg correction method to control for the false discovery rate (FDR) in multiple comparisons. Terms with an FDR-corrected p-value <0.05 or <0.25 were considered significantly enriched for various comparison groups, as described in the Results section. To study the differences in the spectrum of activated signaling pathways for DEGs, gene set enrichment analysis (GSEA) was applied using the built-in FGSEA (fast preranked GSEA) package (with a minimum gene set size of 20, a maximum gene set size of 500, and a number of permutations of 1,000). The following gene sets were downloaded from the curated collections of MSigDB (Molecular Signature Database, https://www.gsea-msigdb.org/gsea/msigdb): “C2 subcollection CGP – chemical and genetic perturbations”, “C4 cancer-oriented genes”, “C6 oncogenic signature gene sets”, “C7 immunologic signature gene sets (subcollection ImmuneSigDB)”, and “H collection (hallmark gene sets)”. Before GSEA analysis, all genes in the DESeq2 output table were ranked by the sign of their FC and statistical significance values ​​according to the formula: (−log10(adj p-val)) × sign(FC). GSEA scores were considered statistically significant at FDR<0.05.

Flow cytometry

Venous blood samples were collected in special vacuum tubes containing 3.8% sodium citrate immediately before surgery. Erythrocytes were removed using osmotic shock by treating whole blood with ammonium chloride buffer containing 0.15 M NH₄Cl, 0.01 M KHCO₃, and 0.0001 M Na₂EDTA (RBC lysis buffer). Leukocytes were then washed and resuspended in phosphate-buffered saline. To analyze the expression of CD markers on the cell surface, the latter were treated with fluorochrome-conjugated monoclonal antibodies, the characteristics of which are presented in Table 2. The optimal working dilution was determined empirically for each antibody clone. The incubation time was 30 min at ambient temperature in the dark according to the manufacturer’s specifications. Measurements were performed using a MACSQuant flow cytometer (Miltenyi Biotec, Germany). The data were further analyzed using MACSQuantify v.2.11software (Miltenyi Biotec, USA). Statistical analysis of the flow cytometry results, and construction of violin plots were performed using the RStudio environment (version R-4.3.2). The lymphocyte population was gated based on forward and side scatter (FSC/SSC) parameters. To more accurately determine the gating boundaries for positively/negatively stained cell populations, ‘fluorescence minus one’ controls were prepared for each of the fluorophores used. In the case of PD-L1 marker expression, in addition to the percentage of positively stained cells, the level of PDL1 antibody binding was calculated as the median increase in fluorescence intensity (DMFI) relative to the autofluorescence level of unlabeled cells in the corresponding detector channel. At least 100,000 events were measured for each probe. The statistical significance of the differences between the FA group and the control group was determined by nonparametric Wilcoxon-Mann-Whitney test (U test; p<0.05 was considered statistically significant).

Table 2. Characteristics of fluorescently labeled monoclonal antibodies used for phenotyping peripheral blood lymphocytes

Results

Transcriptome analysis

As a first step, complete transcriptome profiles of PBMCs isolated from seven FA patient samples (PBMC-FA) and seven control samples were analyzed using principal component analysis (PCA) to assess the distance between sequenced datasets and identify groups of samples for comparison. The similarity/dissimilarity of global gene expression profiles is illustrated by the PCA plot in Figure 1. As can be seen in the plot, the FA datasets form a cluster separated from the two control datasets (Control_1) along the first principal component axis (PC1) and from the three control datasets (Control_2) along the second principal component axis (PC2). Pairwise comparisons of the resulting clusters were then performed, and genes differentially expressed in the FA patient PBMCs were identified using DESeq2 analysis.

Figure 1. Principal component analysis (PCA) of global RNA expression profiles of peripheral blood mononuclear cells (PBMCs) obtained from women with fibroadenoma (FA) and healthy controls (c). The two-dimensional PCA scatterplot shows the distribution of sample transcriptomic profiles along two principal components (PC1 and PC2), which together explain 45% of the data variability; each dot represents one sample, the color of which corresponds to one of the two study groups. As can be seen in the plot, the FA samples form a cluster (indicated by the blue dashed line) that is variably separated from the transcriptomic datasets of the two control clusters (indicated by the red lines and labeled Control_1 and Control_2; the two samples are grouped together with the FA samples).

Comparative transcriptome analysis of PBMC-FA vs. Control_1 identified 286 DEGs meeting the significance criteria of adj p-val<0.05 and |log2FC|>1, most of which (266) were upregulated in PBMC-FA group compared to Control group, while only 20 DEGs were significantly downregulated in FA datasets (Figure 2A, C). Comparison of PBMC-FA datasets with Control_2 resulted in the identification of 48 differentially expressed genes (with adj p-val<0.25 and |FC|>1.5), most of which (39 DEGs) were downregulated in PBMC-FA, while 9 DEGs were found upregulated at the given significance level (Figure 2B, C). Further hierarchical clustering analysis of DEGs confirmed a clear division of the datasets into subgroups of FA patients and controls (Figure 2D, E). The resulting heatmaps indicate, on the one hand, a certain degree of heterogeneity in the PBMC transcriptomic profiles across all compared groups, and on the other hand, the existence of distinct expression patterns of up- and downregulated genes.

Figure 2. Analysis of differentially expressed genes (DEGs) between the transcriptomes of peripheral blood mononuclear cells (PBMCs) from women with fibroadenoma (FA) and the control group. (A) and (B) Volcano plot showing the distribution of differentially expressed genes (DEGs) according to statistical significance and magnitude of change. The y-axis represents the -log10 transformed adjusted p-values ​​ (adj p-val); the x-axis represents the log₂FC values, where FC denotes the fold change of expression obtained by comparing the FA datasets to the Control_1 subset (A) or to Control_2 (B). Individual genes are depicted as dots colored according to the significance threshold and FC of expression threshold: red dots indicate DEGs that showed significantly upregulated expression in PBMCs of the FA patient group (Up); blue dots indicate genes with downregulated expression, and gray dots represent non-DEGs. The plot also shows the gene names for the 25 (A) and 10 (B) DEGs that showed the largest expression change in the FA datasets compared to Control_1 and Control_2, respectively. (C) Total numbers of up- and downregulated DEGs identified via comparing the PBMC sample groups: FA vs. Control_1 and FA vs. Control_2; the full list of genes is available upon request. (D) and (E) Results of hierarchical cluster analysis of DEGs identified between the FA and control datasets; heatmap was generated based on the standard score (or Z-score) values; red color corresponds to upregulated genes, while blue color indicates downregulated genes.

We also noted that the list of upregulated genes in the FA datasets included, on the one hand, a large number of regulators of inflammatory responses and modulators of innate immune response mechanisms, including those mediated by natural killer (NK) cells, macrophages, and dendritic cells. On the other hand, some of the upregulated genes indicate possible dysfunction in a T cell subset (e.g., activation of apoptosis, response to glucocorticoids and oxidative stress, differentiation of Treg lymphocytes). Meanwhile, among the genes showing reduced expression in the FA datasets, we identified B cell markers and various stimulators of T cell function (such as activation, expansion/proliferation, differentiation, or migration), as well as negative modulators of inflammation. These data suggest that abnormalities in the cellular immune compartment, presumably accompanying FA development, may arise from various molecular mechanisms underlying its pathogenesis. To elucidate which biological functions may be affected by the pathological process, we additionally performed GO and KEGG functional enrichment analysis for DEGs, followed by comparative signaling pathway activation pattern analysis via GSEA method.

According to the GOSeq analysis of genes differentially expressed between the PBMC-FA and Control_1 datasets, 486 BP GO terms (with FDR<0.01), 24 MF GO terms (FDR<0.05), 68 CC GO terms (FDR<0.01), and 16 KEGG signaling pathways (FDR<0.25) were overrepresented in the PBMC-FA datasets. Given the hierarchical modular architecture of the GO dictionary and its redundancy (overlapping synonymous terms), the ShinyGO v.0.81 graphical tool (https://bioinformatics.sdstate.edu/go [24]) was used to reduce the dimensionality of the GO-BP enrichment data, visualize it, and highlight the most representative terms. Figure 3 shows a large cluster of terms covering processes regulating T cell activation and proliferation (primarily Th1), monocyte maturation, and inflammatory responses. A closer inspection revealed that these processes were enriched primarily in IFN-γ and IL-12 signaling pathways, while enrichment in IL-1, IL-6, IL-8, and TNFα-dependent pathways indicates activation of inflammatory responses. Another cluster of BP terms represents the processes of mononuclear cell adhesion and directed cell migration. Notably, downregulated genes (representing the minority of DEGs) were associated with negative regulation of mononuclear cell migration, as well as with processes of T cell selection, activation, and differentiation, including IL-4 and IL-10-dependent mechanisms, apoptosis, and Treg functions. Another group of overrepresented terms was associated with various functions of antigen-presenting cell populations (macrophages, dendritic cells), including phagocytosis, cytokine production, reactive oxygen/nitrogen generation, and pattern recognition receptor signaling. Key enrichment results for the MF and CC categories are presented in Table 3. The protein expression products of DEGs involved in the above-mentioned processes are primarily components of the major histocompatibility complex (MHC) and members of signaling cascades triggered by innate immune receptors, antigen recognition receptors, and cytokine receptors, including signaling kinases/phosphatases and adapter proteins/signal transducers. Their subcellular localization relates primarily to cellular structures involved in secretion or, conversely, receptor-mediated endocytosis, or to membrane-bound receptor complexes, as well as to the membrane structures responsible for cell adhesion/movement (Table 3).

Figure 3. Network diagram illustrating the results of Gene Ontology (GO) functional enrichment analysis of differentially expressed genes between the PBMC-FA and Control_1 datasets; significantly enriched terms from the Biological Process (BP) category were visualized using the ShinyGO online tool [24]. To create the diagram, GO-BP terms meeting the FDR<0.01 criterion were sorted by fold enrichment, and the 40 most represented terms were defined (the full list of terms is available upon request). Connecting lines indicate the degree of relatedness between terms (>20% overlapping genes), and their thickness indicates the proportion of shared genes. The brightness of the nodes indicates the level of statistical significance; the size of the node indicates the number of genes annotated with a specific term.

Table 3. Results of Gene Ontology (GO) functional enrichment analysis for differentially expressed genes (DEGs) between the PBMC-FA and Control_1 sample datasets; showing the top 10 significantly enriched terms (ranked by enrichment score) in the Molecular Function and Cellular Component categories (FDR: false discovery rate)

Comparison of the PBMC-FA and Control_2 datasets using GOSeq revealed functional enrichment for 81 GO terms of the BP category (FDR<0.01), 10 GO terms of the MF category (FDR<0.05), and 8 GO terms of the CC category (FDR<0.01), as well as for 12 KEGG signaling pathways (FDR<0.25). The most statistically significantly enriched GO BP categories were humoral regulatory mechanisms mediated by immunoglobulins or some other soluble factors of innate immunity (Figure 4A, B). In terms of fold enrichment score, the Th2-like inflammatory profile, pyrogenic/thermogenic responses, and the associated mechanisms of NO generation, chemotaxis, and phagocytosis were the most over-represented. Another group of terms reflects the general cellular response to stress, including apoptotic processes and antioxidant defense responses, which is interesting since among DEGs, we found the glucocorticoid-inducible gene encoding the zinc-binding protein metallothionein MT1E, as well as the SOD1 gene, the role of which in immune responses is becoming the subject of active research [25]. The enrichment results for the MF category suggest that the DEGs involved in the above-mentioned processes are components of cytokine/chemokine or immunoglobulin receptor complexes (Table 4). KEGG pathway analysis reveals that different functional activities of antigen-presenting cells (macrophages), cytotoxic T/NK cell responses to auto-/alloantigens, pattern recognition receptor-mediated inflammatory responses, and a spectrum of specific cytokines/chemokines may also contribute to the transcriptional landscape of PBMC-FA (Table 5).

Figure 4. Gene Ontology (GO) functional enrichment analysis of differentially expressed genes between the PBMC-FA and Control_2 datasets, with terms belonging to the Biological Process (BP) sub-ontology. (A) The bubble chart shows the 10 terms with the minimum adj p-val (FDR). (B) The hierarchical clustering dendrogram demonstrates the relationships between statistically significant terms: those sharing a higher number of genes in common are grouped together; the size of the blue dots at the branch ends is determined by the significance level (adj p-val) indicated alongside. To construct the dendrogram, GO-BP terms satisfying the FDR<0.01 criterion were sorted by fold enrichment, and a set of the 40 most significantly enriched terms is shown (the full list of terms is available upon request). The ShinyGO online tool [24] was employed to visualize the GO analysis results.

Figure 5. Distribution of differentially expressed genes (DEGs) across chromosomal loci. DEG positions are marked with red dots. Purple segments indicate genomic regions that are statistically significantly enriched in consistently overexpressed DEGs compared to the average gene density in the human genome. Genome screening was performed using the sliding window technique; the window size was set to 4 Mb, and the number of steps in each window was 2. The ShinyGO online tool [24] was employed to visualize the scanning results.

Table 4. The set of most significantly enriched GO terms from the Molecular Function and Cellular Component functional categories in differentially expressed genes (DEGs) between the PBMC-FA and Control_2 datasets; (FDR: false discovery rate)

Table 5. The KEGG signaling pathway enrichment analysis of the differentially expressed genes (DEGs) in peripheral blood mononuclear cells (FA: fibroadenoma; FDR: false discovery rate)

We also investigated whether the identified DEGs could be associated with specific chromosomal regions (or genomic loci). As shown in Figure 5, highly enriched regions were located on chromosomes 2, 5, 11, and 19. This is expected given that chromosome 19, for example, contains a large collection of immunoglobulin superfamily genes, which are expressed in leukocytes and encode receptors that recognize different classes of antigenic molecules or molecular patterns, as well as genes encoding cytochrome and serine protease subfamilies involved in inflammatory responses [26]. Chromosome 5 contains the interleukin gene cluster encoding several hematopoietic cytokines [27], and chromosome 2 includes the IL-1 gene family, which are inflammatory mediators [28]. Thus, it can be concluded that the distribution of genes differentially expressed in the PBMCs of patients with FA is non-random.

Since GO analysis can only extract information about the functions of overrepresented genes but does not allow us to infer the contribution of entire signaling pathways to the observed differences, we performed GSEA of the identified PBMC clusters using gene sets curated by the molecular signature database (MSigDB) [29] to better understand the biological processes that may be involved in the development of FA at the level of cellular immunity. Of primary interest was, of course, the ImmuneSigDB collection of immune-associated gene signatures (subset C7) [30]. Among these gene sets, the following pathways were most enriched in PBMC-FA samples vs. Control_1 (adj p-val/FDR=0.006, NES>1.8, where NES is the normalized enrichment score): (1) innate proinflammatory pathways activated in monocytes/macrophages and dendritic cells in response to stimulation of pattern recognition receptors specific for exo/endogenous ligands (such as NOD2, TLR2, TLR1, TLR4/CD14, TLR7/TLR8, TLR9, TREM1 and lectins) and mediated by MYD88, DUSP1/p38, WPB7/MLL4, CSF2RA and other regulators; (2) pathways responsible for the suppression of Treg differentiation (including LEF1/TCF10-mediated pathways); (3) signaling networks controlling IL-6/STAT3-dependent Th17 cell differentiation, SAP (SH2D1A)-dependent Tfh cell differentiation, and CD4 effector and memory T cell differentiation (Teff_mem, T_cent_mem) in general; (4) pathways associated with inhibition of CD8 T cell activation (downstream of type I and type II interferons, CD3/CD28 costimulatory receptors, IL-2, IL-21, STAT4, and STAT5), CD8 T cell anergy, and CD8 memory T cell formation. The signaling pathways identified by comparing PBMC-FA with Control_2 (FDR<0.2, |NES|>1.4) are likely alternative pathways that can cause inflammatory immune dysregulation. These include: (1) PPARG- and STAT6-dependent mechanisms that trigger the polarization of M2 macrophages (exhibiting profibrotic and proliferative activity) [31] and activate the type 2 immune response in dendritic cells (a recognized factor in persistent inflammation) [32]; (2) mechanisms that control B cell activation and the induction of T cell anergy (mediated by IRF4 and NFATC1); and (3) mechanisms that induce the differentiation of or Tregs (including adipose tissue-resident Tregs) and are controlled by PPARG, DNMT1, FOXP3, and STAT5 [33].

Of particular interest are the differences found in CD8 subsets characterized by the PD1^low/PD1^high phenotype. PD1 expression status (low/high) was shown to correlate with the state of T cell activation/exhaustion [34]. In line with this, GSEA data indicate a general upregulation of PD-1-dependent processes and a shift in the PD-1 status of CD8 T cells during FA development. Further evidence is provided by analysis of the Hallmark (H) collection of key biological processes and physiological states [29]: among the most enriched gene sets, inflammation and T cell activation processes involving interferons, TNFα, and IL6, as well as induction of apoptosis, showed the highest enrichment score (with NES>1.7 and FDR=0.033) (Figure 6).

Figure 6. Plotted results of gene set enrichment analysis (GSEA) comparing the transcriptomic profiles of PMBC-FA vs. the control group (Control_1) using gene sets obtained from the Hallmark MSigBD collection. Nine upregulated gene sets are shown, ranked by the highest enrichment score, normalized enrichment score (NES), and the lowest false discovery rate (FDR). The black lines at the bottom of each plot correspond to the position of the given gene set in the ranked list (ranking was performed as described in the Materials and Methods section); the green lines represent the enrichment score (y-axis) for each position in the ranked list.

Given that the immune microenvironment is an integral part of malignant tumors, a similar role for the immune compartment can be hypothesized in benign tumors (such as FA). Although our transcriptome analysis was limited to peripheral blood and did not include the tumor compartment, we decided to investigate whether the transcriptomic differences identified in PBMCs could be linked to any cancer-related metaprograms or gene networks using the C4 collections from MSigDB (3CA: Curated Cancer Cell Atlas and CGN: cancer gene neighborhoods). When comparing the transcriptional profiles of the PBMC-FA and Control_1 datasets using tumor-specific metaprograms, the highest enrichment score (|NES|>1.5, FDR<0.05) was obtained for inflammatory signaling pathways initiated by macrophages through interferon secretion and upregulation of proteasomal degradation and antigen presentation mechanisms (Figure 7A), which is consistent with the results of the C7 ImmuneSig and Hallmark analyses. It is noteworthy that metaprograms specific to epithelial cells and fibroblasts (in particular, myofibroblasts) were also identified, as well as an EMT metaprogram, which is most likely directly related to the pathology under study, viz., breast fibroepithelial tumors. It is also interesting that FGSEA revealed an increase in interferon-associated antigen presentation activity in epithelial and stromal cells. Furthermore, FGSEA revealed a positive gene signature of PI16, which is recognized as a marker of a proinflammatory fibroblast population capable of attracting immune cells [35] and a negative modulator of Treg differentiation [36]. Enrichment of the signature of cancer-associated fibroblasts is also noteworthy. Only naïve CD4 and CD8 T cell-specific gene sets were negatively enriched (Figure 7A). The enrichment profile of the 3CA gene pathways obtained by comparing the PBMC-FA transcriptome with Control_2 largely overlapped with that when comparing PBMC-FA with Control_1. However, some specific and interesting mechanisms were identified that deserve mentioning: (1) mechanisms underlying the impaired cellular response to stress; (2) mechanisms of cell cycle regulation, including those specific to fibroblasts and epithelial cells; (3) signaling pathways directing the activation of the cellular response to metal ions mediated by class I metallothioneins (MT1) with Zn/Mg binding activity. The involvement of the latter gene family seems quite interesting given the fact that, in addition to their metabolic and antioxidant roles, MT1 are currently considered important regulators of immune responses (including the formation of an immunosuppressive environment) that can occur in response to a wide range of proinflammatory factors [25], as well as members of regulatory networks of carcinogenesis relevant to the biology of both benign and malignant breast tumors [37] (Figure 7B). Regarding gene clusters coregulated by any of the 380 tumor-associated genes in the CGN collection (MSigDB's C4 gene pool), FGSEA identified gene networks organized around key mediators of monocyte/macrophage-initiated inflammatory responses (CASP1, CARD15/NOD2, HCK, TNFRSF1B, S100A4, CD14, MYD88, TNFSF10/TRAIL, MCL1/BCL2L3, SPI1), some of which also act as regulators of apoptosis in lymphocytes/monocytes. 

Figure 7. Gene set enrichment analysis (GSEA) results for gene sets obtained from the MSigBD C4 collection (3CA: Curated Cancer Cell Atlas and CGN: Cancer gene neighborhoods). The PMBC-FA transcriptome was compared with the Control_1 (A) or Control_2 (B) transcriptome, false discovery rate (FDR)<0.05; the horizontal axis shows the normalized enrichment score (NES).

Given that, according to GO and GSEA, certain signs of T lymphocyte exhaustion were observed at the transcriptome level (probably caused by their activation), we decided to use flow cytometry for evaluating possible changes in the expression patterns of four markers of immune exhaustion: immune checkpoint molecules PD-1, PD-L1, TIM-3 and LAG-3 in CD4 and CD8 T cell subsets.

Expression of immune checkpoints in circulating lymphocytes (based on flow cytometry)

Under conditions of chronic stimulation (and, in particular, chronic inflammation), the simultaneous (co-)expression of several immune checkpoints is more important and representative than the expression of any of the immune checkpoint markers individually, since it more deeply characterizes the degree of lymphocyte exhaustion [38]. That is why we decided to test whether it was possible to detect changes in the number of T cells simultaneously positive for three key checkpoint markers (PD1⁺PD-L1⁺TIM3⁺ or PD1⁺PD-L1⁺LAG3⁺) in patients with FA. As can be seen in Figure 8, in the group of patients with FA, there was an average increase in the fraction of cells with the specified phenotype, while no significant differences between the CD4 and CD8 subsets were detected. In the patient group, a relatively wide range of variations in cell counts was also noted, whereas in the control group, cells with the PD1⁺PD-L1⁺TIM3⁺ or PD1⁺PD-L1⁺LAG3⁺ phenotypes were virtually undetectable.

Figure 8. Co-expression levels of three immune checkpoint markers in peripheral blood T lymphocytes obtained from women with fibroadenoma (FA) and healthy controls, measured by flow cytometry. Hereinafter, the group median and interquartile range are shown in box and whisker plots overlaid onto violin plots.

When analyzing the mono-expression of each of the selected markers in the total T lymphocyte population or in its CD4/CD8 subsets separately, a trend towards an increased share of PD-1⁺ T cells in the FA sample set was observed (Figure 9), although it did not reach statistical significance. Taking into account the results of the signaling pathway analysis described above, as well as the published data reporting heterogeneity of the PD-1⁺ T cell population in terms of PD-1 expression level (high, low/intermediate) [34, 39, 40], we decided to investigate whether a shift in the intensity of PD-1 expression could be observed in FA group vs. the control. CD4⁺PD-1⁺ or CD8⁺PD-1⁺ T cell populations were divided into PD-1^high and PD-1^low/int expression groups (high or low/intermediate expression, respectively), and the ratio of %PD-1^low/int to %PD-1^high cells was measured (Figure 9). No statistically significant differences were found between the FA group and the control group, even though there was a clear trend towards an increased relative proportion of CD8⁺PD-1^high T cells. The latter were reported to exhibit more pronounced immune exhaustion [34, 39, 40]. At the same time, the results for the CD4 subset indicated a potential increase in the relative frequency of the PD-1^low/int phenotype.

Figure 9. Expression levels of the PD-1 marker in peripheral blood T lymphocytes from women with fibroadenoma (FA) vs. the control (based on flow cytometry data).

When PD-L1 mono-expression was assessed, a significant increase in the percentage of PD-L1⁺ T cells was detected, as opposed to PD-1 expression: in fact, most T cells in the FA group were PD-L1-positive (Figure 10). However, it is important to point out that the fluorescent signal intensity was very weak, thereby implying a low overall cell surface expression density of this checkpoint marker in both study groups. TIM3 and LAG3 checkpoints also exhibited marked increases in expression levels, with a similar pattern of changes observed in the CD4 and CD8 subsets.

Figure 10. Changes in the expression levels of PD-L1, LAG3, and TIM3 markers in peripheral blood T lymphocytes from women with fibroadenoma (FA) vs. the control (based on flow cytometry data).

Hence, in our study, the results of transcriptome analysis and flow cytometry suggest that changes in the number of T cells (co-)expressing immune checkpoint markers that accompany the development of FA potentially reflect the disruptions in signaling pathways and gene expression patterns described above.

Discussion

Transcriptome analysis of the PBMC fraction is currently considered a promising approach to studying the pathogenetic mechanisms underlying the formation of benign and malignant breast lesions, helping to decipher the heterogeneity of this group of diseases and, consequently, expand diagnostic capabilities and treatment monitoring [41]. Transcriptomic studies of circulating blood cell populations can provide a deeper understanding of the mechanisms of the body’s immune response to the onset and progression of malignant neoplastic processes at the systemic level. However, as already mentioned, although such studies are actively conducted in relation to breast cancer, information on benign breast pathologies is virtually nonexistent (experimental research data presented in the literature concern mainly malignant phyllodes tumors). In this study, we attempted to characterize the changes in the transcriptomic profile of PBMCs obtained from patients with FA as a proof-of-concept study, and for the above reasons, we can compare our results exclusively with the data obtained for invasive breast cancer, of course, taking into account the fundamental differences between malignant neoplasms and benign proliferative lesions of the breast. For example, by sequencing the PBMC transcriptome, Ming et al. demonstrated the existence of two independent immune subtypes of breast cancer that do not overlap with the clinical classification of breast cancer, which is based on the immunohistochemical expression of certain receptors by tumor cells [42]. In some cases, the PBMC transcriptome in breast cancer shows a strong signature of inflammatory gene expression [41]. In our study, we did not aim to identify distinct FA immune subtypes and define their specific gene signatures; however, we were able to demonstrate major intrinsic heterogeneity in the PBMC-FA transcriptome compared to the transcriptomic profile of control PBMCs, which likely reflects different modes of interaction between FA tumors and the immune system. The consistency of these changes is supported by the fact that, following bioinformatic analysis, we were able to identify not only differentially expressed genes but also entire signaling pathways. In this regard, Ming et al. reported an interesting finding that the major transcriptome differences they found between the two immune subtypes of breast cancer PBMCs were primarily associated with cytokine-dependent activation mechanisms, in particular, TNFα-dependent monocyte/macrophage-mediated inflammatory pathways [42]. Several transcriptomic studies of peripheral blood mononuclear cells (PBMCs) in breast cancer revealed an important role of innate immune response regulators, such as Toll-like receptors (TLRs) [41], and reported the possibility of detecting an inhibitory T cell phenotype that can be defined by altered immune checkpoint-dependent signaling. Consistent with this, we also noted a significant impact of inflammation-related mechanisms and innate immune responses on the PBMC-FA transcriptome, along with evidence of inhibition of T cell functions.

Changes in the transcriptomic profile of PBMCs reflect changes in the phenotypic composition of circulating lymphocytes/monocytes. The expression of immune checkpoint molecules in different PBMC populations is currently of great interest, not only because of their potential as targets for cancer immunotherapy, but also because the analysis of immune checkpoint molecules allows assessing the balance between the states of lymphocyte activation and exhaustion as a natural mechanism for controlling the immune response. However, even for breast cancer, changes in immune checkpoint expression in PBMCs are not yet well understood [43], let alone for benign breast diseases. We found an increased frequency of circulating PD-L1⁺ T cells in a set of FA samples, which have a tolerogenic effect on tumor cells and an immunosuppressive effect on effector T cells [44]. Recent studies also showed that increased expression of the PD-L1 molecule on CD4 T cells may be a consequence of their activation upon stimulation by proinflammatory factors and may interfere with the cytotoxic T cell response [45]. In experimental (organoid) models of breast cancer, LAG3 was shown to have an immunosuppressive effect as well, in particular on CD8 T cells [46]. It can be speculated that this effect may be extrapolated to the case of FA, given the fact that we found not only an increase in LAG3 expression but also an increase in its co-expression with PD-1/L1, which may be an indicator of an exhausted T cell phenotype (as in the case of triple-negative breast cancer [47]). Similar observations were made for TIM3 expression levels. It is widely believed that while the expression of a single immune checkpoint molecule is a marker of transient T cell activation, co-expression of multiple checkpoints is a sign of T cell exhaustion, which progresses in the context of chronic inflammation [48]. Indeed, according to our flow cytometry results, peripheral blood immune cells from patients with FA exhibited an increased frequency of simultaneous expression of three checkpoint molecules on CD4 and CD8 T cells and also displayed features of chronic inflammation at the PBMC transcriptome level.

Unlike malignant tumors, benign breast lesions do not pose a serious clinical problem and are not directly associated with immune dysfunction. Nevertheless, studying the immune profile (including the transcriptome landscape) of women diagnosed with FA may provide insight into impairment of putative immunoendocrine pathways which, in turn, may lead to various disruptions in the mechanisms that control epithelial and stromal cell proliferation, tissue architecture, and the normal processes of development and involution of glandular breast tissue. As follows from our results, transcriptomic and phenotypic changes in peripheral blood lymphocytes/monocytes in women with FA coincide to a certain extent with oncogenesis-associated features and metaprograms of cancer. Therefore, further research on this issue will expand our understanding of oncogenic risk factors for early cancer detection and help in the development of differentiated approaches to the clinical management of FA per se.

Limitations of the study

The authors fully acknowledge the limited sample size analyzed at this stage of the study and the need for further expansion, given the particularly high interpatient variability of the peripheral blood immune profile, which is typically influenced by numerous intrinsic and extrinsic (environmental) factors. Another important point is the need for long-term follow-up of patients in the postoperative period to identify correlations between molecular profiling, immunophenotyping, and clinicopathological parameters. In the flow cytometry portion of our study, we did not analyze the expression of immune checkpoints in Treg, which is a goal for future research.

Conclusion

To our knowledge, this is the first study to examine the transcriptomic features of circulating lymphocytes and monocytes in women diagnosed with breast FA and the first to analyze the co-expression pattern of immune checkpoint markers (PD-1, PD-L1, TIM-3, and LAG-3). Our results convincingly demonstrate the involvement of specific immune molecular pathways in the pathogenesis of FA, including inflammatory mechanisms, regulation of innate immunity, humoral response, and antioxidant defense. All of these pathways converge on the immune checkpoint regulatory network that controls the activation/suppression of lymphocytes. For the first time, we demonstrated signs of T cell exhaustion determined by changes in the expression of these immune checkpoint molecules.

Conflict of interest

The authors declare that they have no conflicts of interest.

Funding

The study received no external funding.

Acknowledgments

The authors state that they did not use AI or AI-assisted technologies in writing the article.

Ethical approval

All procedures involving human participants were conducted in accordance with the ethical standards and with the 1964 Declaration of Helsinki and its subsequent amendments. The study was approved by the Committee on Medical Ethics at the Institute of Medicine, Petrozavodsk State University, Republic of Karelia Ministry of Healthcare (protocol No.3, approval date: 21 Oct 2024), and was completed in accordance with the good clinical practice guidelines.