The human hepatocellular carcinoma (HCC) cell lines Huh7, HepG2, and SMMC-7721, and human fetal hepatocyte line LO2 were cultured in DMEM (HCC lines) or RPMI-1640 (LO2). The HCC cell lines Huh7, HepG2, and SMMC-7721 were chosen because they originate from the liver cancer tissues of patients of different ages and with different etiologies, therefore potentially representing distinct histological subtypes of HCC. This diversity allows for a broader investigation of lipid metabolic abnormalities across varying HCC contexts and enhances the generalizability of the findings.

We conducted preliminary experiments with varying the concentration of FICZ (50, 100, 200, and 400 nM) and the treatment duration (12, 24, and 48 h). We found that the cells treated with 200 nM FICZ for 24 h showed the strongest expression of the AHR target genes TIPARP and CYP1A1 based on the qPCR results, justifying the choice of this condition for further study. Cells were grown to 70%–80% confluency and treated with various concentrations of FICZ (50, 100, 200, 400 nM; n=3) for 12 or 24 h. Based on the IC50 value, treatment with 200 nM FICZ for 24 h was chosen for the further experiments. Meanwhile the control groups were treated with DMSO.

HCC or LO2 cells were cultured in 24-well dishes containing poly-L-lysine-treated coverslips. The cells were treated with either FICZ or DMSO for 24 h and then fixed with 4% paraformaldehyde (in phosphate-buffered saline, PBS) at room temperature for 15 min. After fixation, the cells were washed three times with PBS, for 5 minutes each wash. The coverslips were then permeabilized with 0.1% Triton X-100 in PBS for 15 min and subsequently blocked with 5% fetal bovine serum (FBS) to reduce non-specific binding. The cells were then incubated with a rabbit anti-AHR primary antibody (Abcam, 1:300 dilution) at 4°C overnight. Following the primary antibody incubation, the coverslips were washed three times with PBS containing 0.5% Tween 20 (PBST), for 5 min each wash. Next, cross-absorbed Alexa Fluor 488-conjugated goat anti-rabbit secondary antibody (ThermoFisher, 4 µg/mL) was then applied and the samples were incubated at room temperature for 1 h in the dark. After this secondary antibody incubation, the samples were washed three times with PBST, for 5 min each wash. Finally, the coverslips were mounted onto microscope slides using 20 μl ProLong Gold Antifade Mountant with DAPI (ThermoFisher) to stain the cell nuclei. After air drying, the samples were visualized under a fluorescence microscope.

Total RNA was extracted from the treated cells and control cells using Trizol. After lysing the cells with Trizol, chloroform was added, and the aqueous phase was collected and precipitated with isopropanol. The obtained RNA was washed with ethanol, dissolved in DEPC water, and then stored at -80°C. The RNA purity and concentration were assessed using a NanoPhotometer^® system and by gel electrophoresis.

Total RNA samples were processed for preparation of the cDNA library, with sequencing performed by BGI Group using the Illumina HiSeq platform. Eight samples (three replicates each) were sequenced, generating an average of 10.17 Gb of data per sample with an average alignment rate of 88.77%. A total of 124,841 transcripts were identified, including 12,597 novel lncRNAs, 7,264 novel mRNAs, 69,206 known lncRNAs, and 35,774 known mRNAs.

Clean reads were aligned to the human reference genome (hg38, GRCh38) using HISAT2 (14). Transcripts were assembled with StringTie (15), and their coding potential was evaluated using CPC (16), txCdsPredict, and CNCI (17). Transcripts with coding potential and alignments with protein through Pfam-scan (18) were regarded as mRNAs; while the other transcripts were regarded as lncRNAs. The gene expression levels were quantified with RSEM. Functional annotation of the genes was performed using the GEO, Ensembl, NONCODE, and UCSC databases.

Potential lncRNA targets were predicted by constructing a co-expression network of lncRNAs and mRNAs. By calculating the Pearson correlation coefficient between known annotated lncRNAs and mRNAs, lncRNA–mRNA pairs were selected based on an absolute value of the correlation ≥ 0.9, and a significance threshold of p < 0.05. These selected pairs were used to construct the lncRNA–mRNA co-expression network. Data visualization was performed using Cytoscape.

DEGs were identified using the limma package in R (19). Genes with significant expression changes (|fold change| ≥ 2, FDR ≤ 0.05) were selected. The DEG analysis included comparisons between the DMSO- and FICZ-treated groups across different cell lines. We applied a threshold fold change ≥ 2 and FDR ≤ 0.05 to maintain strong, reliable signals, as we found a lower fold-change cutoff would result in an overwhelming number of DEGs.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using the R package clusterProfiler (20).

For the reverse transcription, the RNA and primers were added to a PCR tube with a total volume of 10 µL. The mixture was incubated at 70°C for 10 min, and then quickly cooled on ice for 2 min. Subsequently, to the reaction mixture containing the RNA/primer denaturation solution, a 10 mM dNTP mixture and other reagents were added to obtain a total volume of 20 µL. The mixture was incubated at 42°C for 60 min, followed by 15 min at 72°C. The resulting cDNA was stored at -20°C for later use.

Next, the cDNA was diluted fivefold and mixed with forward and reverse primers, SYBR^® Premix Ex Taq™ (Tli RNaseH Plus) (2×), in a total volume of 20 µL. The PCR reaction conditions were as follows: 95°C for 5 min; 95°C for 10 s, followed by 60°C for 34 s (when the fluorescence signal was collected), for 39 cycles. Melt curve analysis was performed by heating the mixture from 60°C to 95°C, with fluorescence measurements taken every minute.

The reaction mixture was loaded into a 96-well PCR plate and centrifuged to ensure the contents were well-settled at the bottom. The reactions were conducted using a Bio-Rad CFX96 real-time PCR system under the following conditions: 95°C for 2 min; 43 cycles at 95°C for 15 s; 58°C for 5 s; and 72°C for 20 s. Fluorescence data were collected at the 72°C step for each cycle.

The primers were designed using PRIMER5 software or the NCBI online tool and are listed in Table 1 .

Survival analysis was conducted using data from The Cancer Genome Atlas (TCGA). The expression levels of the selected lncRNAs were correlated with the patient survival outcomes. Statistical analysis was used to determine the prognostic significance of these lncRNAs, for identifying potential biomarkers for liver cancer prognosis.

All the data were statistically analyzed using SPSS software (version 13.0, SPSS Inc., Chicago, IL, USA), GraphPad Prism 8.0, or R. One-way analysis of variance (ANOVA) was performed, followed by unpaired t-tests where appropriate. A p-value < 0.05 was considered statistically significant. Using the survival package in R, we performed a proportional hazards assumption test and Cox regression analysis on Overall Survival (OS) using liver cancer data from the TCGA database (GDC Portal).

To study the role of AHR in liver cancer, we induced AHR activation in HCC cells using 200 nM FICZ. We selected three hepatocellular carcinoma cell lines (HepG2, Huh7, and SMMC-7721) that have enriched plasma unactivated AHR, as well as the human fetal hepatocyte line LO2. We treated these cells with FICZ, creating four FICZ-treatment groups, as well as four DMSO control groups for comparison.

We assessed AHR activation and nuclear translocation in the liver cancer cell lines after FICZ treatment using immunofluorescence assays. The results showed a significant increase in AHR fluorescence within the nuclei of the FICZ-treated cells. In contrast, the DMSO-treated controls displayed no such nuclear fluorescence, confirming that FICZ treatment successfully activated AHR, as indicated by the strong fluorescence burst seen in both the nucleus and cytoplasm, while the control DMSO group showed only weak AHR fluorescence, primarily localized in the cytoplasm ( Figure 1 ).

Next, we performed RNA-seq to analyze the transcriptomic changes induced by AHR activation, with the workflow customized to analyze both mRNAs and lncRNAs. Given our interest in both mRNAs and long non-coding RNAs (lncRNAs), we used StringTie for de novo transcript assembly, achieving an average genome alignment rate of 88.77%. We classified the transcripts as mRNAs if they were predicted to have coding potential by at least three of the following four methods: txCdsPredict, CNCI, Pfam-scan, and CPC ( Figure 2A ). Conversely, transcripts that were identified as non-coding by at least three of these methods were classified as lncRNAs ( Figure 2B ).

We then cross-referenced these transcripts with annotated databases to distinguish between annotated and novel mRNAs and lncRNAs. In total, we identified 124,841 transcripts, comprising 69,206 annotated lncRNAs, 12,597 novel lncRNAs, 35,774 annotated mRNAs, and 7,264 novel mRNAs. The distribution of these transcript categories was consistent across the different cell lines ( Figure 2C ).

To further investigate the common AHR target genes in the tested HCC cells, we conducted differential gene expression (DEG) analysis by comparing tumor cell lines with or without AHR activation. Additionally, we included a control group based on a human fetal hepatocyte line to filter out genes that were altered upon AHR activation in non-cancerous conditions. This approach allowed us to focus specifically on AHR-related changes in cancer cells. Our analysis identified 427 significantly differentially expressed lncRNA transcripts with known gene annotations, of which 167 (39.1%) were upregulated and 260 (60.9%) were downregulated ( Table 2 , Figures 3A, B ). We also found 413 significantly differentially expressed mRNA transcripts, with 321 (77.7%) upregulated and 92 (22.3%) downregulated ( Table 3 , Figures 4A, B ). This analysis revealed that AHR exerts broad regulatory effects on both non-coding and coding genes, with a common gene expression signature associated with liver cancer.

Next, we validated the RNA-seq findings by performing qRT-PCR on selected differentially expressed genes. We also included two well-characterized AHR target genes (TIPARP and CYP1A1) as positive controls. We selected the top 10 differentially expressed mRNAs and lncRNAs based on their absolute log fold change values and proceeded with those for which high-fidelity primers could be reliably designed. Consequently, 8 lncRNAs (ASAP1-IT1, NONHSAT221345.1, RMDN2-AS1, lnc-TMEM232-4, lnc-TP53TG5-6, lnc-FGA-2, lnc-DGKK-1, and lnc-FAM237B-2), and 9 mRNAs (BBC3, ANPEP, DIPK1B, ASAP1, RhoBTB1, CPA4, ANKRD1, IGFBP3, and DEFB1) were selected for validation.

We conducted qRT-PCR on Huh7 cells treated with either FICZ or DMSO. The qRT-PCR results demonstrated that, except for lnc-TP53TG5-6, the expression trends for both the lncRNAs and mRNAs in the FICZ-treated cells were consistent with those identified in the RNA-seq analysis ( Figure 5 ). This confirmation supported the reliability of our high-throughput sequencing data and verified the accuracy of the observed gene expression changes.

Next, we analyzed the functional pathways of the differentially expressed RNAs following AHR activation in HCC cells. Gene Ontology (GO) analysis identified significant enrichment in 19 cellular component (CC) terms, 14 molecular function (MF) terms, and 30 biological process (BP) terms ( Figure 6A ). The analysis highlighted the prominence of genes involved in membrane components, developmental and immune processes, signal transduction, biological regulation, and metabolic pathways ( Figure 6B ).

The further classification of 72,686 transcripts using the EuKaryotic Orthologous Groups (KOG) database revealed 25 functional groups. The major categories included signaling mechanisms (12,637 transcripts) and general function prediction (11,545 transcripts), with additional representation in transcription (4,565 transcripts), post-translational modification, protein turnover, and molecular chaperones (4,338 transcripts), cytoskeleton (4,132 transcripts), and intracellular transport/secretion (2,933 transcripts). Less represented categories included cell motility (203 transcripts) and coenzyme transport/metabolism (300 transcripts) ( Figure 6C ).

In parallel, KEGG pathway analysis further elucidated the involvement of the AHR target genes in various cancer-related pathways. Specifically, these genes were prominently linked to PI3K-Akt, VEGF, Notch, and PPAR signaling pathways, as well as cancer-related microRNAs ( Figures 6D–F ). Additionally, AHR activation was associated with pathways related to fatty acid synthesis and metabolism, immune responses, and hormonal signaling, including estrogen, thyroid hormone, and oxytocin pathways. These results indicate that AHR activation significantly influences metabolic and signaling pathways pertinent to cancer progression, although the detailed mechanisms and downstream targets remain to be elucidated. This comprehensive classification underscores the broad impact of AHR activation on a diverse array of cellular functions and processes.

To explore the regulatory roles of the candidate lncRNAs, we constructed a functional network centered around the differentially expressed (DE) lncRNAs and mRNAs. We screened the regulatory relationships between the DE lncRNAs and mRNAs in the TCGA and TargetScan databases and then visualized their co-expression network with Cytoscape ( Figure 7 ).

We identified several significant co-expression pairs. For instance, RMDN2-AS1 was found to be co-expressed with death-associated protein kinase 3 (DAPK3), FAM99B with ferritin heavy chain (FTH1), and ASAP1-IT1 with plexin A1 (PLXNA1), a plasma membrane protein regulated by low-density lipoprotein receptor-related protein 1 (LRP1). Other notable associations included TP53TG5-6 with transmembrane protein adipocyte-related 1 (TPRA1), lnc-BMP6-106 with LDL receptor-related protein 5 (LRP5), and DGKK-1 with insulin-like growth factor binding protein 3 (IGFBP3). Notably, these genes, and consequently their interacting lncRNAs, are involved in the regulation of lipid metabolism, and hence may potentially contribute to the lipid metabolic abnormalities observed in HCC cells. These findings aligned with our KEGG pathway analysis, which highlighted the impact of AHR activation on metabolic pathways, particularly those related to glucose and lipid metabolism.

Overall, our co-expression network analysis suggested that AHR activation modulates glucose and lipid metabolism in HCC at the transcriptional level. This network provides insights into the potential mechanisms underlying hepatocarcinogenesis and suggests directions for future functional studies aimed at understanding the role of lncRNAs in cancer metabolism.

To investigate the prognostic significance of the AHR-dysregulated lncRNAs in HCC, we conducted survival analysis using data from The Cancer Genome Atlas (TCGA). We analyzed the expression profiles of all the differentially expressed lncRNAs following AHR activation in 424 TCGA HCC samples and integrated the associated clinical data. Through this analysis, we identified 10 lncRNAs with significant survival implications and evaluated their potential impact on patient outcomes.

Our analysis revealed that several of these lncRNAs were significantly associated with the prognosis in HCC patients ( Figure 8 ). Specifically, higher expression levels of ASAP1-IT1, RMDN2-AS1, RNF208, and TP53TG5-6 were correlated with a poorer overall survival, while elevated FAM99B expression was linked to improved outcomes. In contrast, among the downregulated lncRNAs, CDC6-1, DGKK-1, NIFK-AS1, and ASH1L-AS1 were associated with an unfavorable prognosis, whereas a higher expression of ADORA2A-AS1 was related to better patient survival.

Notably, ASAP1-IT1, RMDN2-AS1, and TP53TG5-6 showed a particularly strong correlation with reduced overall survival in primary liver cancer patients, emphasizing their potential role as biomarkers for prognosis. These findings complement our earlier results, which implicated these lncRNAs in AHR-mediated regulatory networks, highlighting their dual importance in both disease progression and prognosis.

In summary, these AHR-dysregulated lncRNAs offer valuable insights into the molecular mechanisms of hepatocellular carcinoma and present potential targets for future therapeutic strategies.