The olfactory receptor OR51E2 regulates prostate cancer aggressiveness and modulates STAT3 in prostate cancer cells and in xenograft tumors

Abstract

Background

Despite advancements in the detection and treatment of prostate cancer, the molecular mechanisms underlying its progression remain unclear. This study aimed to investigate the role of the receptor OR51E2, which is commonly upregulated in prostate cancer, in the progression of this disease.

Methods

We investigated the physiological effects of OR51E2 through CRISPR-Cas9-induced monoclonal OR51E2 knockout. We assessed in vitro and in vivo tumorigenicity and conducted transcriptomic and proteomic analyses of xenograft tumors derived from these knockout cells. Furthermore, we analyzed the effects of differences in OR51E2-expression levels in patients from a TCGA cohort.

Results

OR51E2-knockout cells exhibited increased proliferation, migration, adhesion, anchorage-independent colony formation, and tumor growth rates, resulting in a more aggressive cancer phenotype. Omics analyses revealed several potential pathways associated with significant molecular changes, notably an aberration in the STAT3 pathway linked to IL-6 signaling, highlighting a connection to inflammatory pathways. TCGA cohort analysis revealed that prostate cancer patients with low tumor OR51E2 expression had a worse prognosis and a higher average Gleason grade than those with higher expression levels. Additionally, this analysis supported the putative OR51E2-related modulation of the STAT3 pathway.

Conclusions

OR51E2 is regulated throughout prostate cancer progression and actively influences cancer cell physiology affecting cancer aggressiveness. Reduced OR51E2 expression may adversely affect patient outcomes, potentially through alterations in the STAT3 pathway that impact cellular responses to inflammatory signaling.

Background

Prostate cancer is a growing health challenge worldwide and with 1 in 8 men predicted to receive the diagnosis throughout their lifetime, prostate cancer is the second most common noncutaneous cancer in men [1]. In recent decades, both the detection and treatment of prostate cancer have improved substantially [2], however, to further improve patient outcomes, more knowledge on the interplay between different factors and dysregulations that trigger prostate cancer progression is needed. The heterogeneity of prostate cancer and the diverse clinical phenotypes observed among patients may clutter such investigations. Still, by investigating common mutations and regulations, the findings are likely to be more applicable in a clinical setting.

The ectopic olfactory receptor OR51E2 is abundant in the prostate (and was originally termed the prostate-specific G protein-coupled receptor; PSGR) and is further upregulated in prostate cancer in approximately 2 out of 3 patients [3,4,5]. It is a member of the G protein-coupled receptor family (GPCR), a large family of signaling receptors that affect the progression of several cancer types [6,7,8]. In addition to being a candidate gene for prostate cancer biomarker panels [9,10,11], this upregulation indicates a putative role for OR51E2 in prostate cancer development and progression. Indeed, previous investigations utilizing murine models revealed that overexpression of OR51E2 induces chronic inflammation in the prostate and prostatic intraepithelial neoplasia (PIN) and stimulation with β-ionone, a ligand for OR51E2, results in increased tumor growth rates and metastasis [5, 12, 13]. Furthermore, cell studies have indicated that β-ionone affects cellular proliferation, migration, and invasion through OR51E2 signaling and have suggested a role in prostate cancer progression toward a more aggressive, metastatic phenotype [14,15,16,17].

Despite accumulating evidence implicating OR51E2 in prostate cancer pathogenesis, our understanding of its signaling pathways and mechanistic involvement in tumor progression remains incomplete. Among the identified downstream targets of OR51E2 are members of the mitogen-activated protein kinase (MAPK)-family and nuclear factor κB (NF-κB) [18], which are involved in inflammatory responses and are commonly dysregulated in the tumor environment [19]. How these potential targets are linked to cancer physiology in response to OR51E2 activation and how OR51E2 signaling interacts with other molecular players within the complex dysregulated signaling network in cancer have yet to be determined. As a heavily upregulated G protein-coupled receptor in prostate cancer, OR51E2 could provide a key link in understanding the mechanisms underlying prostate cancer progression, which may aid in uncovering novel therapeutic targets associated with this progression. Additionally, much of the current information on OR51E2 signaling is obtained experimentally by β-ionone exposure, which may exert effects not directly linked to OR51E2 binding [18], necessitating studies discriminating the effects of β-ionone and OR51E2 in prostate cancer.

Here, we investigated the role of OR51E2 in prostate cancer progression by constructing an OR51E2 knockout prostate cancer cell line and examined physiological and molecular alterations as well as tumorigenicity, both in vitro and in vivo. Additionally, we analyzed data from a TCGA cohort to investigate the effects of low and high OR51E2 expression cancer phenotypes to assess the clinical significance of OR51E2. By identifying the changes directly associated with OR51E2 regulation and by analyzing the proteomic and transcriptomic profiles associated with the knockout, we aimed to identify specific pathways affected and assess the clinical potential of these OR51E2-associated pathways.

Methods

Cell culture and media

LNCaP human prostate carcinoma cells were obtained from Sigma-Aldrich (FGC clone, RRID: CVCL_1379, St. Louis, MO, US) and were cultured in RPMI 1640 (ATCC modification) medium (Gibco, Gaithersburg, MD, USA) supplemented with 10% heat-inactivated fetal bovine serum (FBS) (Thermo Fisher Scientific, Waltham, MA, USA) and 1% penicillin-streptomycin-amphotericin B (Sigma-Aldrich), hereafter referred to as standard growth medium. The cell cultures were incubated at 37 °C in 5% CO₂. The cells were subcultured at 60–90% confluence and were kept at < 20 passages. Similarly, knockout cell lines were kept at < 20 passages after the knockout was established (details below). The cells were visually inspected for contamination and were screened for mycoplasma contamination by PCR analysis regularly.

Generation of knockout cell line

Transfection of SgRNA and the Cas9 enzyme

To generate CRISPR-Cas9-mediated knockout cell lines, we incubated specific sgRNAs (Supplementary Table S1) with SpCas9 nuclease (Synthego, Redwood City, CA, USA) in Opti-MEM (Gibco) for 15 min to form ribonucleoproteins (RNPs) complexes before transfection. The sgRNAs were designed using Synthego’s CRISPR design tool (https://www.synthego.com/products/bioinformatics/crispr-design-tool) and the guides were chosen with an emphasis on avoiding off-target effects. The cells were transfected with RNPs via nucleofection (4D-Nucleofector System, Lonza, Cologne, Germany) using program P3 CM-138. Each nucleofection was performed in a total volume of 30 µL containing 1.5 × 10⁵ cells, 20 pmol SpCas9, and 180 pmol sgRNA in opti-MEM. To optimize and verify the transfection conditions, we used a pMAXGFP-vector (Lonza), and as a positive control, we used TRAC multi-guide sgRNAs as well as a mock transfection and a negative (Cas9 only) control.

Production of stable monoclonal knockout cell lines

To obtain monoclonal knockout cells following nucleofection, the cells were seeded into 96-well plates at an average density of 0.5 cells per well in standard growth medium. The cells were observed by microscopy twice a week for the first three weeks to determine which wells contained proliferating cells originating from a single cell, and they were then gradually moved to larger containers as they became nearly confluent.

Knockout verification

When the cells had grown sufficiently for further analysis, they were harvested for genomic DNA extraction using QuickExtract DNA extraction solution (Lucigen, Middleton, WI, USA) following the supplier’s protocol, and the target loci were PCR amplified (AmpliTaq Gold 360 Master Mix, Thermo Fisher). The PCR products were purified (Invitrogen PureLink PCR Purification Kit, Thermo Fisher) and Sanger sequenced (Eurofins Genomics, Ebersberg, Germany) to detect putative knockout cell lines. The sequences were analyzed using the ICE CRISPR analysis tool available at https://ice.synthego.com/. All primer sequences are available in Supplementary Table S2.

To verify the OR51E2 knockout at the protein level, we performed a sandwich ELISA with OR51E2-specific antibodies (Nordic BioSite, Täby, Sweden) to measure the OR51E2-protein concentration in the cell lysate from the putative knockout cell lines determined via Sanger sequencing. To prepare the samples for the ELISA, we washed the cells with ice-cold PBS and harvested them with a cooled cell scraper in ice-cold RIPA buffer (Thermo Fisher) supplemented with 1x HALT protease inhibitor cocktail (Thermo Fisher) and 2 mM activated sodium orthovanadate (Sigma-Aldrich). The cell lysates were then transferred to microcentrifuge tubes, and kept on ice for 60 min, while briefly vortexed every 10–15 min. The tubes were then spun down at 16,000 x g for 20 min at 4 °C and the supernatants were collected. The total protein concentration was determined using a BCA protein assay kit (Thermo Fisher) and 15 µg of total protein (diluted with distilled water) was used for the sandwich ELISA. The ELISA was carried out according to the manufacturer’s protocol, including the recommended negative controls.

Potential off-target effects (all genomic sites with < 6 mismatches) were examined using mRNA-seq data from the knockout cells to investigate the sequences of all potential off-target sites within genes that were expressed in the cells. The details regarding this mRNA-seq dataset can be found below.

Cell proliferation and viability assay

We performed an MTT (3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide) assay to assess the putative effects of OR51E2 on proliferation using an MTT assay kit (Abcam, Cambridge, UK). We seeded 5000 cells in wells of 96-well plates at time 0 and incubated both wildtype and knockout cells (and blank media as a control) for three hours with MTT substrate after 0, 24, 48, 72, and 96 h following the manufacturer’s recommendations. For each cell line, we prepared a standard curve with known cell numbers, which was used to calculate the number of proliferating cells at the specified time points. An MTT assay was also used to examine androgen-induced proliferation using 10 nM dihydrotestosterone (DHT) for 24 h following the same procedure, to examine if the lack of OR51E2 altered this response, as earlier studies have suggested an androgenic component of OR51E2-signalling [20]. All the assays were replicated three times with two technical replicates.

Wound healing migration assay

To examine the effects of OR51E2-signalling on the migration capacity of the cells, we performed a wound healing (scratch) assay, using inserts to form the scratch area (Abcam). Cells were seeded densely (5 × 10⁵ cells per well in a 24-well plate) and were allowed 24 h to attach and form a monolayer. Then, the inserts were removed, the cells were carefully washed with PBS, and the medium was replaced with 1% FBS growth medium to limit proliferation during the assay. The wells were imaged with two images per well and 3–4 wells per treatment (technical replicates) after 0, 24, 48, 72, and 96 h. The images were analyzed using ImageJ (v. 1.53) to automatically delineate the edges of the migrating cells and calculate the migration area. The experiment was replicated three times.

Xenograft tumor generation

Male CIEA NOG mice (Taconic Biosciences A/S, Ejby, Denmark) were used to establish xenograft tumors, to examine in vivo tumorigenicity and for molecular analyses. Briefly, 4 × 10⁶ cells of either wildtype LNCaP or LNCaP OR51E2^−/− were suspended in 1:1 growth medium: Geltrex (Gibco) and injected subcutaneously into the flank region of the mice in a volume of 0.3 mL. Tumor growth was monitored by measuring the length (L) and width (W) of the tumor externally with a caliper twice a week after the first appearance of the tumor. Tumor volumes (V) were calculated as V = π/6 x L x W², which assumes an ellipsoid shape of the tumors. These tumor volumes were used to follow the growth and to calculate the specific growth rate (SGR) of the xenograft tumors using the formula SGR = (V2-V1)/Δt, where V1 and V2 are the volumes of the xenograft tumors at the first and last measurement, respectively and Δt is the days between the two measurements. Mice were euthanized by cervical dislocation when the tumor volume reached an estimated size of 1000 mm³. Tumors were harvested, snap-frozen in isopentane cooled to -40 °C with dry ice, and stored at -80 °C.

All mouse studies were conducted according to the animal welfare policy of Aarhus University (http://dyrefaciliteter.au.dk) and with the Danish Animal Experiments Inspectorate’s approval (License number: 2021-15-0201-00828).

Soft-agar colony formation assay

To determine the anchorage-independent growth capacity of the cells, as a proxy for in vitro tumorigenesis to compare with the in vivo tumorigenicity, we performed a soft-agar colony formation assay [21]. The bottom layer of agar (Chemsolute, Roskilde, Denmark) was prepared with 0.6% agar in growth medium and the top agar with 0.3%. The cells were seeded in the top layer with 1000 cells per well in a 48-well plate and fresh medium was supplied every second day for 14 days, after which the colonies were manually counted under a microscope. The experiment was replicated three times with 12 wells per condition each time.

Cell-to-matrix adhesion assay

To assess potentially altered adhesive capacities following OR51E2 knockout, based on observations from the xenograft tumor generation (see Results), we examined the cell-to-matrix adhesion of the cells, by coating a 48-well plate with a 1:30 dilution of Geltrex basement membrane matrix (LDEV-Free, hESC-Qualified, Reduced Growth Factor, Gibco) in RPMI 1640 for one hour at room temperature and washed gently with PBS. Thirty thousand cells were added to each well in 300 µL of growth medium and allowed to attach for 30 min in the incubator. The cells were carefully washed with PBS, fixed with freezer-cold methanol, washed, stained with crystal violet, and washed again. The cell-bound crystal violet in the adhered cells was recovered with methanol (15 min at room temperature with constant agitation) and the absorbance of this methanol/crystal violet solution at 590 nm was measured, which provided an indirect measurement of the number of adhered cells. We visually inspected the wells under a microscope between each step to verify that the adhered cells were not lost during the protocol.

3D geltrex drop invasion assay

The potential changes in the cells’ invasive capacity induced by the lack of OR51E2-signalling were assessed using a 3D Geltrex drop invasion assay following the protocol of the 3D Matrigel drop invasion assay described in detail previously [22]. Briefly, the cells were mixed well with undiluted Geltrex, and 10 µL droplets with 5 × 10⁴ cells were quickly placed centrally in wells of a 24-well plate. The droplets were allowed 15–20 min to solidify and were then carefully covered with growth medium. The droplets were imaged every 24 h with 4 images per well and the cells’ invasive capacity and migration/expansion were analyzed by delineating the droplets’ edges and the leading edges of the migrating cells using ImageJ’s freehand tool, which were used to calculate the area covered by cells. We used 3–4 wells as technical replicates and repeated the experiment three times.

Western blotting

Proteins were extracted from cell cultures, and concentrations were analyzed as in the ELISA protocol. For STAT3-phosphorylation analysis, cells were pretreated with 10 ng/mL IL-6 for 30 min before extraction. Proteins were then mixed with SDS sample buffer, denatured at 95 °C for 5 min, and 15 µg was loaded into a 4–15% TGX Stain-free gel (Bio-Rad Laboratories, Hercules, CA, USA) for electrophoresis. Proteins were transferred to PVDF membranes, blocked with 5% skimmed milk in TBST for 2 h, and probed with primary antibodies overnight at 4 °C (Supplementary Table S3). After washing, membranes were incubated with HRP-conjugated secondary antibodies for 2 h at room temperature. Protein bands were visualized using an ECL detection system, with GAPDH as a loading control. Densitometric analysis was performed using Image Lab software v. 5.0 (Bio-Rad).

Quantitative real-time PCR

To verify mRNA levels of target genes, we extracted total RNA from the xenograft tumors or cell lines using an RNeasy spin-column kit (Qiagen, Hilden, Germany) and removed genomic DNA using the gDNA Eliminator spin column from the same kit following the manufacturer’s instructions. We measured the RNA content and purity using a NanoPhotometer system (Implen GmbH, Munich, Germany) and performed cDNA-synthesis of the mRNA using SuperScript IV First-Strand Synthesis System (Invitrogen) with oligo d(T)₂₀-primers. We designed two sets of primers for each gene using the Primer3 online utility (https://primer3.ut.ee/) and subsequently performed a primer-BLAST using the online NCBI BLAST tool [23] to further verify the uniqueness of the primers (Supplementary Table S2). The RT-PCR was performed using a PowerTrack SYBR Green Master Mix (Applied Biosystems) on a Stratagene Mx3000P system (Agilent Technologies, Santa Clara, CA, US) with a total volume of 20 µL per well. The program used was: enzyme activation (2 min at 95 °C), two-step amplification x 40 (5 s at 95 °C, 30 s at 60 °C), followed by a melting curve acquisition (15 s at 95 °C, 1 min at 60 °C, 15 s at 95 °C with a ramp rate of 0.075 °C/s). ROX was used as a reference dye and the primer efficiency of each primer pair was calculated using the LinRegPCR software v. 2017.0. All samples and genes were run in triplicate (technical replicates) with the inclusion of a housekeeping gene (GAPDH) in each run for normalization as well as negative controls without either reverse transcriptase or without template added. The full protocol was replicated three times. The relative expression of the target mRNA to GAPDH was calculated and the Log₂ fold changes relative to the control (wildtype cell line in standard growth medium) were calculated using the formula:

$$\:{Log}_{2}\:fold\:change={Log}_{2}\left(\frac{{{E}_{target}}^{{\varDelta\:Ct}_{target}}}{{{E}_{GAPDH}}^{{\varDelta\:Ct}_{GAPDH}}}\right)\:$$

,

where E is the primer efficiency and ΔCt is the difference in Ct between the average Ct of the control replicates of the gene and the sample Ct.

RNA sequencing

Sample Preparation and sequencing

To identify the OR51E2-associated changes in the xenograft tumors, we extracted total RNA from the tumors using an RNeasy spin-column kit (Qiagen). We performed quality control of the extracted RNA using a Bioanalyzer 2100 system (Agilent) and shipped the samples on dry ice to Novogene (Sacramento, CA, USA), who prepared the mRNA library, sequenced to ~ 50 M read-pairs of 150 bp paired-end reads using an Illumina NovaSeq system, and performed additional quality control of both the input material and output data.

RNA sequencing data handling

After receiving the raw sequencing data, we built a pipeline using the public server at https://usegalaxy.org for the initial analyses of the data [24]. Briefly, we used Trimmomatic, FastQC, and HISAT2 to prepare the data, for quality control, and to align the reads to a reference genome, respectively. Since the tumors contained both human and mouse cells, we aligned the reads to both human and mouse genomes (hg38 and mm10), imported these aligned files into R, and used the R package ‘XenofilteR’ v. 1.6 to remove all mouse reads from the dataset. The output files were re-uploaded to https://usegalaxy.org, where featureCounts was used to quantify gene expression and DESeq2 was used to perform a differential expression analysis and to normalize the gene counts. This dataset was downloaded and low reads (< 10 average counts in either of the groups) were filtered out. This dataset was used for further analyses as described below.

Proteomic analysis

Sample preparation

To investigate the molecular changes associated with OR51E2 at the protein level, we performed shotgun proteomics on tumor samples from both WT and KO tumors. We lysed the tissue samples in 1.5% sodium dodecyl sulfate (SDS) with 1x protease inhibitors (Sigma-Aldrich), homogenized the ice-cooled tissue with a homogenizer probe for 3 × 15 s, and sonicated the homogenate with a Branson probe sonicator (Thermo Fisher) with an amplitude of 40% for 3 × 15 s. The extracted protein was quantified using a BCA kit (Thermo Fisher). For each sample, an aliquot of 25 µg total protein was used. We further prepared the samples for proteomics analysis by reducing the disulfide bonds with 10 mM Tris (2-carboxyethyl) phosphine (TCEP) (Sigma-Aldrich) at 65 °C for 1 h. Then the samples were alkylated by 12.5 µM iodoacetamide (Acros Organics, Belgium) at room temperature in the dark for 45 min. Proteins were then precipitated with freezer-cold acetone and stored overnight at -20 °C. The precipitated proteins were pelleted by centrifugation at 14,000 x g at 4 °C for 10 min, were allowed to dry, and were then digested with 1 µg of sequencing grade modified trypsin enzyme (Thermo Fisher) in 50 µL 50 µM ammonium bicarbonate at 37 °C overnight. The samples were then dried down using a speed vacuum and were reconstituted in 50 µL 0.1% formic acid in water for LC/MS analysis.

LC/MS analysis

Two µL of reconstituted tryptic peptides were loaded onto a C18 trap column (Thermo Fisher) coupled to a Dionex Ultimate Rapid Separation Liquid Chromatography system (Thermo Fisher) at a rate of 5 µL/min for 10 min. The peptides were separated by reverse-phase chromatography on a 25 cm C18 analytical column (New Objective Inc, Littleton, MA, US) packed in-house with BEH C18, 130 Å, 1.7 μm particle size (Waters Corp, Milford, MA, US). A column heater (PST Phoenix S&T) was used to heat the column to 60 °C. The chromatographic gradient was established by changing the mixture composition of mobile phase A (0.1% formic acid in water) and mobile phase B (0.1% formic acid in acetonitrile) at a constant flow rate of 0.3 µL/min. The gradient conditions were as follows: mobile phase B kept at 2% for the first 6 min then slowly ramped up to 35% over 124 min followed by a rapid increase to 85% over the next 5 min with a 5 min hold time. The analytical column was re-equilibrated for 10 min before the next sample injection. Each sample was analyzed in triplicates. Eluted peptides were ionized by applying 2.2 kV using a nano-spray source (Thermo Fisher Scientific) coupled to an Orbitrap Eclipse Tribrid mass spectrometer (Thermo Fisher Scientific). For the MS1 scan settings, the orbitrap mass resolution was set to 240,000, the scan mass range between 375 and 1800 m/z, and the maximum injection time was set to 35 ms with an intensity threshold of 5000. The precursor ions were subject to higher collisional disassociation (HCD) with a fixed collisional energy of 28%. MS2 scans were acquired using Top-Speed with a cycle time of 1 ms. Dynamic exclusion was enabled for 30 s using a mass tolerance of ± 10 ppm. The fragmented ions were detected in a linear ion-trap using turbo scan rate.

Proteomics data handling

For each LC/MS run, the resulting raw data file was processed through two separate searches using Byonic 2.11.0 (Protein Metrics, San Carlos, CA): one search was conducted against a Swiss-Prot database featuring the reference human proteome (2022; 20,645 entries) and the other against a Swiss-Prot database of the reference mouse proteome (2022; 17,380 entries). The search parameters were set to include trypsin digestion allowing up to two missed cleavages, a precursor mass tolerance of 0.5 Da, and a fragment mass tolerance of 10 ppm. We applied fixed modifications for cysteine carbamidomethylation and variable modifications for methionine oxidation and asparagine deamination. Peptide identifications were refined by excluding any peptides with a false discovery rate (FDR) higher than 1%. Additionally, peptides identified in both human and mouse database searches were excluded to ensure a stringent analysis of uniquely human proteins with non-homologous mouse peptides using an in-house R script. Spectral counts of all identified peptides were extracted and analyzed. These counts were normalized using z-score transformation and differential protein abundances across samples were assessed using Student’s t-tests. The P values were adjusted using a Benjamini-Hochberg P value adjustment. In our analysis, only proteins detected in at least 40% of the samples were considered, and proteins with the lowest counts (maximum count < 2) were excluded.

Omics analyses

For both the RNA-seq data and the proteomics data, we identified the differentially expressed genes and proteins as described above. From the normalized data, we prepared complete linkage Euclidean-clustered heatmaps of the top differentially expressed genes and proteins (by P values) using the R package ‘ComplexHeatmap’ v. 2.12.1, volcano plots using the R package ‘EnhancedVolcano’ v. 1.14.0, and MA plot using the ‘ggmaplot’ function of the ‘ggpubr’ R package v. 0.6.0. For both the RNA seq data and the proteomics data, we investigated the significantly up- and downregulated terms through the use of a pre-ranked Gene Set Enrichment Analysis (GSEA) where the genes/proteins were ranked by their fold change. The GSEA was run against all gene sets containing between 10 and 500 genes in the mSigDB Hallmark gene sets collection using the R package ‘fgsea’ v. 1.22.0. To further investigate the dataset, we performed enrichment analyses using ‘genekitr’ v. 1.2.5 on the GO: BP databases on all terms smaller than 1000 to identify both up- and downregulated terms. For these analyses, we used a cutoff of FDR < 0.05 for both proteins and mRNA and log₂ FC greater than ± 1 for the proteomics dataset and ± 1.5 for the transcriptomics dataset. The Venn diagrams for visually comparing the two datasets were made using the ‘VennDiagram’ R package v. 1.7.3. Selected genes from the RNA-seq data were experimentally verified by qPCR (Supplementary Fig. S1).

Database analyses

We accessed the prostate adenocarcinoma dataset ‘Firehose Legacy’ generated by the TCGA Research Network (https://www.cancer.gov/tcga) from the cBioportal for Cancer Genomics (http://cbioportal.org) and obtained the counts data from the same dataset from the Broad Institute at https://gdac.broadinstitute.org/. This dataset contains clinical and analytical data from 499 prostate cancer patients including pathological descriptions (i.e. Gleason scores), disease-free survival data, mRNA expression data, and more. We analyzed the mRNA expression data of this data set by grouping the patients based on their expression level of OR51E2. This grouping was done by ranking the patients by their OR51E2 mRNA expression level and picking the lower and upper quartiles. We examined these patients’ disease progression, tumor Gleason grades, and OR51E2 correlations with other genes (for the full dataset) as well as performed a GSEA as described above for the RNA-seq data. We used the R packages ‘survival’ v. 3.5-8 and ‘survminer’ v. 0.4.9 for the survival analysis and plots.

Statistics and figures

Figures were prepared in R v. 4.3.3 (mainly using ‘ggplot2’ v. 3.5.0 or specific packages already mentioned above) and GraphPad Prism v. 10, where statistics were also performed. The normality of residuals was assessed by qqplots prior to the downstream statistical analyses. In addition to the statistics mentioned specifically for the handling of the omics data as described above, we used a two-tailed unpaired Student’s t-test when comparing two groups with single measurements and an ANOVA with a Tukey post hoc analysis when comparing more than two groups. For measurements with several time points from the same cells, we used a mixed-model ANOVA with replicate as a random effect and time and treatment as fixed effects, and we performed a pairwise comparison for each time point with a Tukey P value correction as a post hoc test. We used the R package ‘nlme’ v. 3.1–164 and ‘emmeans’ v. 1.10.1 for these mixed model ANOVAs. We used a Benjamini-Hochberg P value adjustment to control the false discovery rate for multiple comparisons in all other cases where multiple comparisons were performed. We analyzed the difference in the tumor growth rates as survival data, counting a size of > 1000 mm³ (where the tumors were harvested) as an event, and we used a Cox proportional hazard model with a likelihood ratio test to investigate the effect; the same test was used for the TCGA survival data. All results in the text and figures are reported as means ± standard deviations unless otherwise noted. Differences were considered significant when P (or FDR) < 0.05.

Results

OR51E2 regulation affects prostate cancer cell physiology and morphology

The knockout of OR51E2 in LNCaP cells (Supplementary Fig. S2) resulted in apparent differences from the wildtype cell line. Consistent with previous findings, OR51E2 had a significant role in regulating cell proliferation, as evidenced by the 1.3x increased proliferation rate in LNCaP OR51E2^−/− cells relative to wildtype LNCaP cells (Fig. 1a, Mixed model ANOVA, P_{cell type} = 0.006, P_interaction < 0.001). This increased proliferation was also apparent in their anchorage-independent growth and colony formation (Fig. 1b and c, t-test, P = 0.01). LNCaP OR51E2^−/− cells exhibited distinct, more “mesenchymal-like” morphology, with notably less clumping of cells, despite an increased proliferation, likely due to decreased cell-cell adhesion capacity (Fig. 1d). However, the knockout resulted in an increased cell-to-matrix adhesive capacity (Fig. 1e, t-test, P = 0.03). These discrepancies in the knockout’s effect on adhesive capabilities suggest that OR51E2 influences cellular adhesion differently depending on the specific conditions.

The invasive capacity following OR51E2 knockout was unchanged (Fig. 2a, c, Mixed-model ANOVA, P_{cell type} = 0.55), as determined by the 3D Geltrex drop assay. The wound-healing assay revealed substantially increased migration capacity in OR51E2^−/− cells (Fig. 2b, d, Mixed-model ANOVA, P_{cell type} < 0.001, P_interaction < 0.001). Overall, these in vitro results indicate a central regulatory role of OR51E2 in prostate cancer cells, at least in cells with high OR51E2 expression levels, and further suggest that the elimination of OR51E2 regulation may lead to a more aggressive cancer phenotype.

Fig. 1

OR51E2 -knockout elicits apparent morphological and physiological changes in LNCaP cells. (a) Proliferation assay results measured with an MTT assay at the indicated time points relative to the seeded cell number at time 0. Asterisks indicate a significant difference between the cell types of P < 0.05 at the indicated time point. (b) Soft agar colony formation results. The count is the total count from the 12 wells used in each replicate. (c) Representative images of the colony formations. The scale is the same for both images. The scale bar represents 500 μm. (d) Representative images of the cell morphology of the two cell types, the scale is the same for the two images and the scale bar represents 200 μm. (e) Cell-to-matrix adhesion assay results measured on plates coated with Geltrex. For all three assays, the black symbols and bars and the solid line indicate wildtype LNCaP measurements, whereas the white symbols and bars and dashed line indicate LNCaP OR51E2^−/− measurements. n = 3 for all assays

Fig. 2

OR51E2-knockout affects the migration but not invasion in LNCaP cells. (a) Representative images of the Geltrex drop assay. The edge of the Geltrex drop is indicated with a yellow dashed line and the edge of the invading cells is marked with a blue dotted line. Scale bar = 500 μm. (b) Representative images of the wound healing assay. The assay was carried out in 1% FBS growth medium with 10 ng/ml IL-6 or vehicle. Scale bar = 500 μm. (c) Results of the Geltrex drop assay. Black symbols and solid lines represent the wildtype LNCaP measurements whereas white symbols and dashed lines represent the LNCaP OR51E2^−/− measurements. No statistical differences in invasive capacity were identified. (d) Results of the wound healing assay. The colors of the symbols and lines represent the same as in c. The circles represent control treatment (DMSO), and the squares represent cells treated with 10 ng/mL IL-6. n = 3 for both c and d

Xenograft tumor development is slower without OR51E2 despite an increased growth rate

Counterintuitively, xenograft tumors developed markedly slower when mice were inoculated with LNCaP OR51E2^−/− cells compared to wildtype LNCaP xenografts, and the median time from inoculation to full tumor development, defined as a size > 1000 mm³, was 90.5 days and 65 days, respectively (Fig. 3a, b, Cox proportional hazard model, likelihood ratio, P < 0.001). During xenograft tumor development, biweekly volume measurements revealed that despite a slower overall tumor development, the growth rate of OR51E2^−/− tumors, once established, was increased by 2.55 times in the knockout cell line (Fig. 3c 3d, t-test, P = 0.007). Thus, the proliferation-reducing effect of OR51E2 appears to be similar in cell culture and xenografts and the slower overall tumor development stems from a regulation of mechanism(s) in the establishment of the tumor seed – potentially mechanisms similar to those involved in metastasizing cancers in patients, such as epithelial-mesenchymal transition (EMT).

Fig. 3

Xenograft tumor growth is affected by OR51E2-knockout in LNCaP. (a) The growth curves of all individual tumors measured externally in the flank with a caliper. The tumors were harvested shortly after reaching 1000 mm³. (b) Same tumor growth data as panel a, shown as mean ± S.D. (c) The time between the first detection of a tumor in the mice and the harvest date, i.e. approximately the time the tumor grows from ~ 200 mm³ to > 1000 mm³. (d) The specific growth rate of the xenograft tumors in the time between the first and last measurement of the tumors. Black symbols and bars represent wildtype LNCaP xenograft measurements and white symbols and bars represent LNCaP OR51E2^−/− xenograft measurements. n = 9 for wildtype and 8 for knockout

OR51E2 knockout impacts molecular pathways associated with cancer development and progression in prostate cancer cells and xenografts

To investigate the molecular changes linked to OR51E2 regulation, we examined the proteome and the mRNA transcriptome of the xenograft tumors. A total of 6472 of 16,023 mRNA transcripts and 1245 of 3336 proteins were differentially expressed after filtering out low reads. Of these, 1416 genes and 548 proteins, had a Log₂ fold change greater than the threshold of ± 1.5 for mRNA data and ± 1.0 for proteomics data (Figs. 4a and 5a). These substantial molecular differences underscore OR51E2’s regulatory role in the cancer cells, reinforcing the findings from the physiological cell studies. Enrichment analyses and gene set enrichment analyses (GSEA) revealed different aspects of the OR51E2-linked cellular regulation at the mRNA and protein levels. At the protein level, the top downregulated processes in the OR51E2-absent tumors were associated with various categories of metabolic processes as well as pathways including mTOR signaling, G2-M checkpoint, and E2F targets. The top upregulated processes included amino acid metabolism-related pathways, cell adhesion, extracellular matrix (ECM) organization, and actin-related filament and cytoskeletal organization; the top pathway upregulated in the GSEA analysis was EMT (Fig. 4bd and Supplementary Fig. S3). Notably, OR51E2 has previously been associated with metabolomic changes in cancer cells (as our data support) following treatment with an OR51E2 agonist, but we have not investigated this part of the dataset further [16].

Fig. 4

Proteomics analysis of the xenograft tumors. (a) Volcano plot of the dataset with the top-regulated proteins based on both FC and P values labeled. (b) Top up- and downregulated Hallmark gene sets determined by a Gene Set Enrichment Analysis. (c) Heatmap of the top 150 regulated proteins by P value. (d) Results of the enrichment analysis on the GO: BP gene sets. For all figures, the direction of the comparison is OR51E2 knockout vs. wildtype, i.e. the direction indicated is the change in the cells harboring the OR51E2 knockout

For the mRNA transcriptomic dataset, the metabolic processes that dominated the proteomic profile were largely absent among the top downregulated genes. Instead, retinoic acid metabolism and steroid hormonal processes were upregulated in addition to adhesion, ECM, and structural pathways as in the proteomics dataset. Processes involved in cell adhesion were also found among the significantly downregulated terms, supporting the in vitro finding that the adhesive properties of the cells were regulated in both directions depending on the circumstances. Among the other downregulated pathways at the mRNA level, most are associated with intercellular signaling, usually involved in synaptic pathways. As in the proteomics profile, mTOR-signaling also appeared slightly downregulated at the mRNA level, and additionally, we identified a strong downregulation of genes associated with hypoxia (Fig. 5c e and Supplementary Fig. S4).

While the differences between the two omics-datasets may seem large judging by the difference in the top regulated datasets and genes/proteins, overall, they are similar, especially when comparing the significantly different genes/proteins, where only 62 of 16,151 genes/proteins were significantly different in opposing directions in the two datasets (Fig. 5f). These differences could be a result of true biological differences between mRNA and proteins or outliers in the analysis. These genes/proteins did not have a substantial effect on the overall analysis.

Fig. 5

mRNA transcriptomics analysis of the xenograft tumors. (a) Volcano plot of the dataset with the top-regulated genes based on both FC and P values labeled. (b) MA plot of the dataset with the same genes as in A labeled. (c) Top downregulated Hallmark gene sets determined by a Gene Set Enrichment Analysis. No gene sets were statistically upregulated in this analysis. (d) Heatmap of the top 150 regulated genes by P value. (e) Results of the enrichment analysis on the GO: BP gene sets. (f) Venn diagrams comparing the mRNA transcriptomics analysis to the proteomics analysis. The first Venn diagram shows the total list of input genes/proteins after filtering out low measurements, whereas the second diagram shows the distribution of significant regulations within the overlapping genes/proteins. For all figures, the direction of the comparison is OR51E2 knockout vs. wildtype, i.e. the direction indicated is the change in the cells harboring the OR51E2 knockout

The omics results suggested that the knockout of OR51E2 resulted in a more mesenchymal-like phenotype, and we examined this effect by stimulating EMT in vitro. We used 10 nM DHT for androgen stimulation instead of the commonly used TGF-β, as LNCaP is TGF-β receptor-negative. The results were not conclusive regarding OR51E2’s involvement in EMT, as the panel of EMT markers we examined elicited differing results (Supplementary Fig. S5). If OR51E2 has a role in regulating EMT in prostate cancer, as the omics data suggest, it is likely affected by other currently unknown factors in concert with the effects of OR51E2. Without additional support, we cannot confidently determine this putative role, and determining this will require substantial additional experimental evidence warranting a study of its own.

The mTOR pathway was downregulated at both the mRNA and protein levels following OR51E2 knockout. However, the involvement of this pathway was previously examined pharmacologically and was determined not to be directly affected by OR51E2 stimulation, so the regulatory effects we observed are likely secondary effects [15]. Thus, we investigated another potential pathway underlying these results: the STAT3 pathway. STAT3 is an essential transcription factor involved in several cellular pathways and has been linked to changes in migration, proliferation, adhesion, and apoptosis, among other functions [25,26,27]. In prostate cancer, STAT3 is involved in promoting metastasis and EMT and pharmacological inhibition of STAT3 reduces tumor growth [28, 29]. In this dataset, the STAT3 pathway was slightly downregulated at a transcriptional level (P = 0.0267, FDR = 0.154, GSEA). We detected an almost threefold increase in the amount of STAT3 protein in the LNCaP OR51E2^−/− cell line relative to the wildtype but found a similar level of phospho-STAT3(Tyr705) following 10 ng/mL IL-6 exposure for 30 min in both cell lines (Fig. 6). While we have not examined the cause of this difference further, modulation of the STAT3 pathway may underlie several of the physiological effects associated with OR51E2.

Interestingly, the knockout of OR51E2 also led to a significant decrease in key genes and proteins associated with prostate cancer, such as KLK3 (PSA), ETV1, and AMACR (Figs. 4 and 5). These findings suggest that OR51E2 may directly regulate these genes or influence a pathway that affects their levels, further implicating OR51E2 in prostate cancer progression.

Fig. 6

STAT3 regulation is altered in OR51E2 -knockout LNCaP cells. (a) Western blot showing the difference in STAT3 and phospho-STAT3(Tyr705) between the wildtype LNCaP and LNCaP OR51E2^−/− cells with and without 10 ng/mL IL-6 treatment for 30 min. The uncropped blots are shown in Supplementary Fig. S6. (b) and (c) Summarized data of the densitometric analysis of band intensity from western blots of STAT3 and pSTAT3, respectively, relative to the intensity of GAPDH. Black symbols and bars represent wildtype cells and white symbols and bars represent knockout cells

Database analysis revealed a worse prognosis for prostate cancer patients with low OR51E2 levels

To examine potential clinical aspects of OR51E2, we used the TCGA Firehose prostate adenocarcinoma cohort to investigate the correlations between OR51E2 mRNA expression levels and patients’ clinical attributes. The Gleason scores of the patient cohort were correlated with OR51E2 mRNA expression, with lower Gleason-score patients having higher OR51E2 levels than higher-grade patients (Fig. 7a, one-way ANOVA, P < 0.001). The variation in the mRNA expression levels in a patient cohort data is large, and the OR51E2 expression level is not a good predictor of Gleason grade in isolation. Still, the data highlight a change in OR51E2 expression levels as the prostate cancer progresses – and as seen in the in vitro prostate cancer cell line experiments, a lower level of OR51E2 invokes a more aggressive cancer phenotype. We do not propose that OR51E2 expression alone determines prostate cancer aggressiveness, but OR51E2 regulation is likely involved in modulating cancer molecular physiology in patients as opposed to being a passively regulated bystander in the cancer cell. The same pattern applies to the disease-free survival data, which indicates a more severe outcome for cancer patients with low OR51E2 levels (Fig. 7b, Cox proportional hazard model, likelihood ratio, P = 0.02). The size of the survival dataset is, however, limited due to the usual slow progression of prostate cancer. Several genes associated with prostate cancer were among the genes mostly correlated to OR51E2 expression on the mRNA level, which could suggest common upstream (or downstream) regulation (Supplementary Fig. S7 and S8).

The GSEA results of the TCGA data suggested that low OR51E2 leads to an upregulation of genes associated with EMT (Fig. 7c). Among the other results, genes associated with STAT3 regulation, inflammation, NF-κB signaling, apoptosis, and the p53 pathway were significantly upregulated. The androgen response gene set was the most downregulated gene set in this analysis. Collectively, data from the TCGA cohort strongly support the role of OR51E2 in regulating several processes associated with cancer aggressiveness.

Fig. 7

Effects of differences in OR51E2 expression in the TCGA Firehose dataset analysis. The patients were divided into two groups based on their expression level of OR51E2, with the “High” group representing the upper quartile and the “Low” group the lower quartile. (a) Boxplot showing the distribution of patients by their OR51E2 expression level (RSEM values) and tumor Gleason score. The boxes are colored by the proportion of high/low OR51E2 expression within the Gleason-grade group. The green symbols indicate patients with high OR51E2 expression, and the blue symbols indicate patients with low OR51E2 expression. The P value is the result of a one-way ANOVA. (b) Survival curve analysis of the groups with the 95% confidence interval indicated as well as the number at risk and censoring plots shown. n = 244. The P value is the result of the likelihood ratio from a Cox proportional hazard model. (c) Gene Set Enrichment Analysis results showing all significantly up- or downregulated Hallmark gene sets based on a differential gene expression analysis of the “Low” vs. “High” patients’ mRNA expressions (i.e. the effects of a low OR51E2 level)

Discussion

Identifying drivers of prostate cancer progression is essential for advancing therapy. In this study, we focused on the ectopic olfactory receptor OR51E2, which is significantly upregulated in prostate cancer. Through physiological measurements on an OR51E2 knockout prostate cancer cell line, xenograft tumors, and transcriptomic and proteomic analyses alongside data from a prostate cancer patient cohort, we found that OR51E2 serves as a central regulator of prostate cancer cell physiology by affecting a wide range of pathways and processes. Loss of OR51E2 regulation in LNCaP cells led to a distinct, more aggressive cancer phenotype, and impacted the STAT3 pathway, potentially explaining novel parts of the mechanisms underlying the altered physiology in prostate cancer.

Stimulation of OR51E2 in prostate cancer cell lines by the agonist β-ionone inhibits proliferation [14, 15, 20], suppresses migration [30], and enhances invasion [13, 15]. Our use of OR51E2^−/− LNCaP cells provided robust validation of these findings; the reversed responses on proliferation and migration observed in these cells relative to wildtype cells support that these effects are indeed mediated by OR51E2. We did, however, not detect differences in the invasive capacity, whereas the tendency was similar to previous findings [13, 15]. Methodological differences or the relatively small sample number could explain this divergence in our study. The counterintuitive finding of reduced initial growth of OR51E2^−/− xenograft tumors complements previous observations of faster growth rates in OR51E2 overexpressing xenograft tumors [5]. This reduced growth rate appeared to result from slow initial growth or establishment rather than reduced proliferation once the tumor was established, consistent with our in vitro findings. While the exact mechanisms involved in this change are unknown, we noticed a substantial reduction in the clumping of the cells in vitro and found a significant effect on several genes and proteins involved in cell adhesion in both omics datasets. The in vitro tumorigenesis appeared to reflect this in vivo difference poorly, as indicated by the soft agar colony formation assay, thus the differences in xenograft tumor development following OR51E2 knockout could be associated with an altered response to an undetermined stimulus in vivo.

The observed differences in the protein levels of STAT3, a transcription factor associated with increased metastasis, proliferation, and resistance to apoptosis [25], could represent a key link in explaining how aberrations in OR51E2 expression impact these parameters in vitro. However, the comparable effects of IL-6 stimulation on STAT3 activation in both cell types suggest that the IL-6/JAK/STAT3 pathway, a primary pathway associated with STAT3 activation, remains functional despite regulatory disparities. A potential explanation for the altered STAT3 levels could be linked to mitogen-activated protein kinases (MAPK). OR51E2 impacts several MAPK-pathways, including p38, p44/42, SAPK/JNK, and ERK1/2 [14, 31,32,33] and MAPK and STAT3 signaling may be either directly or indirectly connected, thus the changes in STAT3 could be a result of an affected shared upstream regulation, a coupling of pathways, or a compensatory mechanism [19]. The precise role of OR51E2 in regulating STAT3 thus remains unclear and elucidating their connection requires further experimental investigation.

OR51E2 is involved in regulating the inflammatory pathway through the activation of NF-κB [5]. In a transgenic mouse model, overexpression of human OR51E2 induced a state of chronic prostatic inflammation, a recognized hallmark of cancer [34]. Subsequently, this overexpression led to PIN and, in conjunction with experimentally induced loss of phosphatase and tensin homolog (PTEN), prostate cancer [35]. NF-κB regulates several EMT-related genes and plays a pivotal role in this process. Thus, NF-κB activation may account for at least part of the differential expression of EMT-associated proteins observed between wildtype and OR51E2^−/− xenograft tumors in this study. Moreover, the incomplete effects on all EMT genes and proteins observed in this and other studies may stem from the fact that NF-κB constitutes only one branch of the broader EMT pathway. Follow-up studies are needed to establish which role OR51E2-regulation plays in this complex transition in cancer cells.

While most studies have focused on stimulation or overexpression of OR51E2, a previous study examining decreased OR51E2 levels using a knockdown approach in a prostate cancer cell line found that migration and proliferation were suppressed, which contrasts with our knockout approach and the typical effects of OR51E2 stimulation [35]. This discrepancy could be due to differences in experimental methods, residual OR51E2 expression, or off-target effects. However, our full gene knockout approach and examination of potential off-target sites reduce these risks.

OR51E2 upregulation in prostate cancer tissue was one of the initial findings once the receptor was identified [36]. Despite discrepancies in the exact OR51E2 levels measured in prostate cancer patients, studies indicate a tendency for OR51E2 levels to be highest in low-grade prostate cancer and PIN, compared to high-grade prostate cancer [15, 20], as supported by the TCGA dataset analyzed in this study. While caution is warranted in drawing firm conclusions from these results, the lower OR51E2 regulation in higher Gleason-grade cancers likely influences the cancer cells’ physiology and contributes to prostate cancer progression and aggressiveness based on observed in vitro and in vivo effects in this study and other studies. Future studies should further investigate the expression of OR51E2 in primary prostate samples to establish a detailed expression profile across different Gleason grades. Although inter-patient variation appears substantial ( [15, 20] and Fig. 7a), OR51E2 has the potential to enhance the current panel of prostate cancer biomarkers. This is similar to findings using OR51E2 mRNA in tissue or urine samples [9,10,11], and particularly promising for early detection, where the levels are generally the highest.

Using a single cell line limits our study’s generalizability, however, due to its very high endogenous OR51E2 expression relative to other well-studied prostate cancer cell lines [13], which was necessary for examining gene function, we chose to use LNCaP as our cell line model. While other cell lines with lower OR51E2 levels were considered, they were deemed less informative for this purpose (refs [5, 15, 37]. and Supplementary Fig. S9). Future studies expanding on the findings from LNCaP should target exploring multiple cell lines with varying OR51E2 levels, including overexpression models and rescuing OR51E2 after knockout. Previous research indicates that further overexpression in LNCaP yields modest effects, suggesting it already has near-maximal OR51E2 regulation [35].

Our findings underscore OR51E2’s significant regulatory role using a prostate cancer cell line and bring novel insights into potential mechanisms to investigate further in future studies. The interplay between OR51E2, MAPK family members like ERKs, and STAT3 signaling, which is critical in cancer-associated gene regulation, remains underexplored and could reveal novel therapeutic targets with future investigations [25, 33]. Additionally, OR51E2’s upstream regulation by factors like IL-6 and potential androgen sensitivity warrants further exploration [37, 38], as well as investigating the impact of prostate cancer-associated mutations on OR51E2 expression. The discrepancy between the in vivo and in vitro tumorigenicity of the cells is peculiar, and further investigation, either by stimulation with a drug panel or by investigating the different properties of the xenograft tumors and colonies, could reveal interesting mechanisms underlying tumorigenicity. Increasing our understanding of these mechanisms could offer deeper insights into the heterogeneity and progression of prostate cancer.

Conclusions

OR51E2 is among the most regulated genes in prostate cancer; influencing multiple cancer-related pathways, including the STAT3 pathway, as demonstrated in this study. Our findings suggest that the role of OR51E2 in prostate cancer progression is likely bigger than previously understood. Beyond its established functions, we identified novel regulatory roles, particularly concerning the STAT3 pathway, which may contribute to the cancer’s aggressiveness. This expanded understanding of OR51E2 offers new insights into its potential as a therapeutic target once the generalizability and mechanisms are uncovered. Given the heterogeneity of prostate cancer and the trend towards personalized medicine, further research into OR51E2’s upstream regulation and its impact on different stages of cancer progression is essential. This study also raises new hypotheses regarding tumor formation and cellular mechanisms for further testing in prostate cancer cells with varying levels of OR51E2 expression that could provide valuable insights into the findings’ generalizability and the clinical potential of this receptor. Building on findings like ours, such studies could lead to innovative therapeutic approaches that improve patient outcomes.