Linking Coregulated Gene Modules with Polycyclic Aromatic Hydrocarbon-Related Cancer Risk in the 3D Human Bronchial Epithelium

Abstract

Exposure to polycyclic aromatic hydrocarbons (PAHs) often occurs as complex chemical mixtures, which are linked to numerous adverse health outcomes in humans, with cancer as the greatest concern. The cancer risk associated with PAH exposures is commonly evaluated using the relative potency factor (RPF) approach, which estimates PAH mixture carcinogenic potential based on the sum of relative potency estimates of individual PAHs, compared to benzo[a]pyrene (BAP), a reference carcinogen. The present study evaluates molecular mechanisms related to PAH cancer risk through integration of transcriptomic and bioinformatic approaches in a 3D human bronchial epithelial cell model. Genes with significant differential expression from human bronchial epithelium exposed to PAHs were analyzed using a weighted gene coexpression network analysis (WGCNA) two-tiered approach: first to identify gene sets comodulated to RPF and second to link genes to a more comprehensive list of regulatory values, including inhalation-specific risk values. Over 3000 genes associated with processes of cell cycle regulation, inflammation, DNA damage, and cell adhesion processes were found to be comodulated with increasing RPF with pathways for cell cycle S phase and cytoskeleton actin identified as the most significantly enriched biological networks correlated to RPF. In addition, comodulated genes were linked to additional cancer-relevant risk values, including inhalation unit risks, oral cancer slope factors, and cancer hazard classifications from the World Health Organization's International Agency for Research on Cancer (IARC). These gene sets represent potential biomarkers that could be used to evaluate cancer risk associated with PAH mixtures. Among the values tested, RPF values and IARC categorizations shared the most similar responses in positively and negatively correlated gene modules. Together, we demonstrated a novel manner of integrating gene sets with chemical toxicity equivalence estimates through WGCNA to understand potential mechanisms.

Figure 1.. Module-trait relationships and genes of interest (A) Module-trait relationships matrix. WGCNA identification of gene modules correlated to PAH cancer risk, calculated through RPF.

Seven module eigengenes (MEs) were highly correlated with RPF. Degree of correlation is colored, with red as positive and blue as negative, and p-values in parentheses. (B) MEturquoise genes of Interest. The most positively correlated module eigengene (“MEturquoise”) was prioritized for evaluation; genes of interest that had expression levels correlated to MEturquoise are shown. Expression levels are Z-score normalized. (C) MEred genes of interest. The most significantly negatively correlated module eigengene (“MEred”) was priotized for evaluation; genes of interest shown, ranked by RPF, with expression levels Z-score normalized.

Figure 2.. Network category table of gene modules MEturquoise and MEred that are most significantly correlated to PAH RPF.

Pathway enrichment analysis was conducted, and a broader network category-level significant enrichment score was calculated. Numbers indicate number of pathways enriched under each network category, and ave FDR (average false discovery rate) shown in parentheses. Distinct biological networks were enriched by MEturquoise and MEred, as well as several network categories commonly enriched, notably Inflammation.

Figure 3.. Gene networks showing most significantly enriched networks associated with prioritized gene modules significantly correlated to PAH RPF. (A) MEturquoise significant enrichment of the Cell cycle S phase gene network, and (B) MEred significant enrichment of the Cytoskeleton actin filament gene network.

Node sizes denote betweennenss centrality measures calculated in each network. A larger node, or “bottleneck gene” thus represents a gene or biomolecule that stands in between the shortest paths connecting the maximum number of neighboring nodes to each other.

Figure 4.. Module-trait relationship correlation matrix linking module eigengenes with risk values.

WGCNA identification of gene modules significantly correlated to risk values and cancer categorizations (RfC, RfD, inhalation unit risk, oral slope factor, RPF, and IRIS and IARC carcinogen classes). Correlation matrix shows that several module eigengenes (MEs) were highly correlated with PAH cancer risk values (RPF), as well as correlated in common between RPF, OSF, and IARC classes. Unsupervised hierarchical clustering was performed on the correlation values and the dendrogram identifies RPF and IARC as having the most similar module correlation patterns. Degree of correlation is colored, with red as positive and blue as negative, and p-values in parentheses.