Toxicogenomic analysis of publicly available transcriptomic data can predict food, drugs, and chemical-induced asthma

Abstract

Background

: With the increasing incidence of asthma, more attention is focused on the diverse and complex nutritional and environmental triggers of asthma exacerbations. Currently, there are no established risk assessment tools to evaluate asthma triggering potentials of most of the nutritional and environmental triggers encountered by asthmatic patients.

Purpose

 The objective of this study is to devise a reliable workflow, capable of estimating the toxicogenomic effect of such factors on key player genes in asthma pathogenesis.

Methods

Gene expression extracted from publicly available datasets of asthmatic bronchial epithelium were subjected to a comprehensive analysis of differential gene expression to identify significant genes involved in asthma development and progression. The identified genes were subjected to Gene Set Enrichment Analysis using a total of 31,826 gene sets related to chemical, toxins, and drugs to identify common agents that share similar asthma-related targets genes and signaling pathways.

Results

Our analysis identified 225 differentially expressed genes between severe asthmatic and healthy bronchial epithelium. Gene Set Enrichment Analysis of the identified genes showed that they are involved in response to toxic substances and organic cyclic compounds and are targeted by 41 specific diets, plants products, and plants related toxins (eg adenine, arachidonic acid, baicalein, caffeic acid, corilagin, curcumin, ellagic acid, luteolin, microcystin-RR, phytoestrogens, protoporphyrin IX, purpurogallin, rottlerin, and salazinic acid). Moreover, the identified chemicals share interesting inflammation-related pathways like NF-κB.

Conclusion

Our analysis was able to explain and predict the toxicity in terms of stimulating the differentially expressed genes between severe asthmatic and healthy epithelium. Such an approach can pave the way to generate a cost-effective and reliable source for asthma-specific toxigenic reports thus allowing the asthmatic patients, physicians, and medical researchers to be aware of the potential triggering factors with fatal consequences.

Keywords:

• toxicogenomic
• transcriptomic
• GSEA
• chemical-induced asthma

Introduction

It is widely accepted that asthma has a multifactorial etiology, where genetic predisposition plays an essential role in disease susceptibility,¹ while environmental factors play a critical role in disease development and progression.² Due to the rise in the incidence of asthma, there is a growing concern over the environmental exposures that may trigger asthma exacerbations.³ Although many theories were suggested regarding how exposure to drugs, toxins, chemicals, and infections can participate in asthma development and/or exacerbation, the exact mechanism is still not fully understood.⁴

Asthma was linked to food allergy as children with food allergies have a higher risk of developing food-induced episodes of asthma that can end up with anaphylaxis; nevertheless, this link is not fully understood yet.⁵ Changes in dietary habits were suggested as a possible cause of increased asthma prevalence⁶ in developed and developing countries. Besides food, air pollution can adversely influence lung function in asthmatic individuals,⁷ but which particles in the air can precisely trigger such an effect is still a matter of debate between researchers except for few well-studied examples.⁸ Exposure to chemicals at work is a significant risk factor for occupational asthma and should be brought to the attention and awareness of every asthmatic patient.⁹ Occupational asthma should be distinguished from the non-immunologic asthma-like syndrome¹⁰ called Reactive Airways Dysfunction Syndrome (RADS), which develops after a single high-level exposure to a pulmonary irritant.¹¹ Many substances used in consumer products are associated with occupational asthma or asthma-like syndromes.¹² Besides occupation-induced asthma, common household chemicals can be another uncountable trigger for asthma in adults.¹³ Drug-induced asthma, especially aspirin-induced asthma, is well-defined, relatively common, and often an underdiagnosed asthma phenotype.^14,15

Currently, there are no ideal asthma risk assessment tools for food, drugs, occupational, and household chemicals. Moreover, there is no means of prediction of potential respiratory sensitization for all possible food or environmental items that we encounter in our daily life.¹⁶ Only a few of these tools are available in the clinical setting, with a limited list of items.¹⁷ Recently, toxicogenomic investigation of different toxic agents’ interaction with the cellular genome improved our understanding of the effect of different chemicals, hazardous agents, drugs, and environmental stressors on different cellular and biological systems. Through multi-omics analysis, the response of all genes to chemical exposure can be examined in order to gain a more comprehensive insight into the potential hazards of that toxicant.¹⁸ Although toxicogenomics was proposed to be a useful tool in health risk assessment,¹⁸ this approach has not been tried yet for asthma triggers’ assessment. Since the bronchial epithelium is the key player in asthma initiation and progression that orchestrates airway inflammation and remodeling, toxicogenomic analysis of bronchial epithelium in asthma is mandated.¹⁹

In this study, we used an in-house bioinformatics pipeline that has shown a remarkable performance in clustering complex diseases previously using publicly available omics data.²⁰ We aimed at identifying the effect of dietary, environmental, and occupational influences on genes that are differentially expressed between healthy and severe asthmatic bronchial epithelium. Therefore, this appraoch can facilitate the development of a comprehensive toxicogenomic database that can link and predict asthma susceptibility or progression in response to a given chemical.

Materials and methods

Bioinformatics approach: microarray analysis

To identify differentially expressed genes in asthmatic patients’ bronchial epithelium (in both small and large airways) compared to healthy controls, publicly available transcriptomic datasets from Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) were extracted. We decided to use dataset (GSE64913) due to its appropriate design, a complete characterization and proper categorization of patients as well as being a representative of the two extremes of the disease (healthy versus severe asthma). The study was done using Affymetrix Human Genome U133 Plus 2.0 Array, which has the advantages of complete coverage of over 53,000 transcripts for analysis. Additionally, this dataset shows the effect of sampling of the bronchial tree as it included central and peripheral airway samples from each participant. Accordingly, we hypothesized that genes that are differentially expressed between severe asthmatic and healthy bronchial epithelium in both central and peripheral airways must have a role in the initiation or progression of the disease.

We used a novel in-house R Bioconductor based pipeline as described previously by Hamoudi et al.²⁰ The pipeline is composed of 5 steps: (1) preprocessing and QC assessment of the downloaded raw microarray image files, (2) normalization to remove background noise and (3) filtration of nonvariant probes between severe asthmatics and healthy controls to (4) precisely identify differentially expressed genes (DEG). Finally, the DEG between the two groups will be used for (5) Gene Set Enrichment Analysis (GSEA) to identify top pathways where the identified genes are enriched. Such an approach will give us a clear list of genes that may participate in the pathogenesis of severe asthma. Figures 1 and 2 outline the pipeline steps used in this study.

Figure 1 Flowchart outlining the steps of the bioinformatics approach to identify differentially expressed genes in severe asthmatic bronchial epithelium compared to healthy controls.

Abbreviations: GEO omnibus, Gene Expression Omnibus; RMA, Robust Multiarray Averaging; GC-RMA, GeneChip RMA; MAS5, Affymetrix Microarray Suite 5.

Figure 2 The flowchart of the bioinformatics approach to identify gene sets related to chemical, toxins, and drugs.

Raw microarray image processing and normalization

Raw CEL files (n=70) that stores the results of the intensity calculations on the pixel values were extracted, then the dataset underwent pre-processing and normalization separately. Affy, Robust Multiarray Averaging (RMA), GeneChip RMA (gcRMA), Affymetrix Microarray Suite 5 (MAS5) packages of R Bioconductor statistical software version 3.0.2 were applied to normalize and remove the background noise. gcRMA and MAS5 expression values were used for the next non-specific filtering based on the coefficient of variation (CV). The CV was calculated as the mean/standard deviation of each probe across all cases.

Non-specific filtration

To filter out non-variant genes, only probes with a MAS5 value of 50 or more and CV value of 10–100% in the gcRMA across all cases, were passed and intersected to obtain a common set of variant probes. Out of the 54,675 probes present in the chip, only 9682 probes passed the filtration process. These filtered probes were annotated, collapsed to their corresponding genes using GSEA software (http://software.broadinstitute.org/gsea/downloads.jsp) by choosing probes with the maximum expression for each gene.²¹ The housekeeping probes, along with those that are not assigned to a gene, were excluded. Hence the resultant filtered probes were the only variant probes as per the GSEA manual.

Limma package to identify DEG

R Bioconductor Limma package was used to identify DEG between severe asthma and healthy controls. Out of the 6014 filtered genes, 225 genes with an adjusted p-value less than 0.05 were identified to be differentially expressed between severe asthma and healthy controls. To visualize top pathways and biological processes shared by the DEG gene list, a simplified and customizable web portal (http://www.metascape.org) was used.²² The gene list enrichment analysis was carried out with the following ontology sources: KEGG Pathway, GO Biological Processes, Reactome Gene Sets.

GSEA

The resultant 6014 filtered gene list was used as input for the GSEA to identify the significantly enriched pathways among gene sets related to chemical, toxins, and drugs, as shown in Figure 2. 31,826 gene sets, were downloaded from two major resources : DSigDB (http://tanlab.ucdenver.edu/DSigDB/DSigDBv1.0/) and DrugMatrix (ftp://anonftp.niehs.nih.gov/ntp-cebs/datatype/Drug_Matrix/) databases. DSigDB organizes drugs and small molecules-related gene sets into four collections based on quantitative inhibition, and drug-induced gene expression changes data.²³ The DrugMatrix database is one of the world’s most massive toxicogenomic reference resources. Table 1 shows the details of the gene sets, the gene coverage, and the number of sets included in each set. The results of GSEA were ranked according to the nominal p-value (<0.05) and false discovery rate (≤0.25) as described previously.²⁴

Table 1 Details of datasets extracted from DSigDB and DrugMatrix and used for GSEA

Collection	Description	Unique Number of Genes	Number of Gene Sets
D1: FDA Approved	FDA Approved Drug Gene Sets.	1,288	1,202
D2: Kinase Inhibitors	Kinase Inhibitors Gene Sets based on in vitro kinase profiling assays.	407	1,220
FDA	FDA Approved Kinase Inhibitors.	341	28
HMS LINCS	Kinase inhibition assays extracted from HMS LINCS database.	381	90
MRC	Kinase inhibition assays extracted from MRC Kinome Inhibition database.	137	157
GSK	GSK Published Kinase Inhibitor Set (PKIS), kinase inhibitors used as chemical probes.	116	204
Roche	Kinase Inhibitors profiled by Roche.	153	570
RBC	Kinase Inhibitors profiled by Reaction Biology Corporation.	246	99
KinomeScan	Kinase Inhibitors profiled by DiscoveryRx using KinomeScan technology.	374	72
D3: Perturbagen Signatures	7,064 gene expression profiles from three cancer cell lines perturbed by 1,309 compounds from CMap (build 02).	11,137	1,998
CMAP	7,064 gene expression profiles from three cancer cell lines perturbed by 1,309 compounds from CMap (build 02).	11,137	1,998
D4: Computational Drug Signatures	Drug signatures extracted from literature using a mixture of manual curation and by automatic computational approaches.	18,854	18,107
BOSS	Text mining approaches of drug-gene targets using Biomedical Object Search System (BOSS).	3,354	2,114
CTD	Curation of targets from Comparative Toxicogenomics Database (CTD).	18,700	5,163
TTD	Manual curation of targets from the Therapeutics Targets Database (TTD).	1,389	10,830
DrugMatrix database	The DrugMatrix database is one of the world’s largest toxicogenomic reference resources	5209	7876

Cell of origin

ARCHS4 is a web resource that makes the majority of published RNA-sequencing data from human and mouse available at the gene and transcript levels. This resource was used to determine which cell type or tissue can express the genes that are differentially expressed between severe asthmatic and healthy bronchial epithelium and are enriched in a given gene set.

Finding a common pathway between identified chemicals

In order to identify common pathways targeted by most of the identified chemicals in the GSEA step, we used the Comparative Toxicogenomics Database (CTD) batch query webtool (http://ctdbase.org/tools/batchQuery).²⁵ All the earlier identified drugs and chemicals were uploaded to the query tool to search for genes and pathways that were reported to be affected by the queried chemicals. The tool will generate a list of pathways where the given chemical affects genes related to that pathway significantly (adjusted p-value <0.05). Only pathways that are shared by at least 50 percent and above of the identified chemicals are selected. As illustrated in Figure 3, a schematic flowchart of this step is outlined.

Figure 3 The flowchart outline using the Comparative Toxicogenomics Database (CTD) batch query tool (http://ctdbase.org/tools/batchQuery) to identify common pathways targeted by most of the GSEA-identified chemicals. (A) All the earlier identified drugs and chemicals were uploaded to the query tool to search for genes and (B) pathways that were documented to be affected by queried chemicals. The tool will generate (C) a list of pathways where the given chemical affects genes related to that pathway significantly (adjusted p-value <0.05). (D) Only pathways that are shared by at least 50 percent and above of the identified chemicals are selected.

Results & discussion

Transcriptomic analysis reveals significant enrichment of genes related to cell division between asthmatic and healthy bronchial epithelial cells

Our analysis identified 225 differentially expressed genes between severe asthmatic bronchial epithelium and healthy bronchial epithelium, as shown in Figure 4A and B. Furthermore, the identified genes shared common pathways related to epithelial cell differentiation, response to growth factors, extracellular stimulus, mechanical stimulus, and wounding (Figure 4C). Interestingly, pathways related to the response to toxic substances and organic cyclic compounds were among the top enriched pathways. These findings indicate that genes altered by environmental substances might play a significant role in asthma development and/or progression to severe asthma.

Figure 4 Gene Set Enrichment Analysis (GSEA) of the differentially expressed genes between severe asthmatic bronchial epithelium (n=22) and healthy bronchial epithelium (n=37) in GSE64913. (A) Distribution of the identified genes ranked according to their position (B) Heatmap image generated from the 2952 DEG between severe asthma and healthy controls which were later filtered into 225 genes (C) the top enriched pathways whether upregulated or downregulated in severe asthma compared to healthy controls using metascape (http://metascape.org): a gene annotation and analysis online resource that generates a graphical representation.

Genes that are differentially expressed in the asthmatic bronchial epithelium are targeted by specific diets, plants products, and plants related toxins

Our further analysis revealed that the significant differentially expressed genes in asthmatic epithelium compared to healthy controls are targets for many substances that have not been previously associated or documented to trigger asthma, as shown in Table 2. Additionally as shown in Table 3, these substances can be categorized into three subgroups: (1) Occupational hazards, (2) Drugs, (3) Dietary factors: plant, plant toxins and food. This is substantial as most of the asthmatic individuals are not explicitly aware that such factors might have a potential effect on their disease status.

Table 2 List of the significantly enriched pathways related to chemicals, toxins, and drugs for the genes that showed significant differential expression in severe asthmatic bronchial epithelium compared to healthy controls

#	Gene Set Name	Size	Enrichment score	Normalized Enrichment score	Nominal p-value	False Discovery Rate q-value	Familywise-error rate p-value	Rank at Max	Leading Edge
1	Purpurogallin	28	0.630629	1.827431	0.002024292	0.0122542	0.033	1538	tags=64%, list=26%, signal=86%
2	Tyrphostin AG 538	16	0.631126	1.846969	0.003937008	0.0136133	0.028	797	tags=38%, list=13%, signal=43%
3	7-amino-4-hydroxy-2-naphthalenesulfonic acid	15	0.640253	1.799482	0.002016129	0.0146326	0.049	176	tags=27%, list=3%, signal=27%
4	Corilagin	17	0.676499	1.8472	0.003898636	0.02042	0.028	964	tags=53%, list=16%, signal=63%
5	4-hydroxytamoxifen	17	0.588731	1.724107	0.012219959	0.0285716	0.096	699	tags=35%, list=12%, signal=40%
6	Salazinic Acid	15	0.646698	1.672565	0.011857707	0.0403463	0.147	1700	tags=73%, list=28%, signal=102%
7	Ellagic Acid	27	0.559782	1.639618	0.004008016	0.0463795	0.183	1067	tags=48%, list=18%, signal=58%
8	Baicalein	15	0.607701	1.556171	0.032719836	0.076446	0.273	1067	tags=47%, list=18%, signal=57%
9	Curcumin	25	0.542949	1.509353	0.024439918	0.0802608	0.333	978	tags=48%, list=16%, signal=57%
10	SB 202190	22	0.489577	1.518684	0.032258064	0.082731	0.323	699	tags=36%, list=12%, signal=41%
11	Acrylamide	15	0.586047	1.661663	0.009451796	0.1950874	0.477	125	tags=20%, list=2%, signal=20%
12	Myricetin	40	0.547651	1.654645	0.01713062	0.1982278	0.487	1174	tags=48%, list=20%, signal=59%
13	Cupric Oxide	131	0.413138	1.663593	0.012526096	0.2041864	0.471	936	tags=29%, list=16%, signal=34%
14	Microcystin RR	20	0.547977	1.66728	0.00990099	0.2108996	0.462	462	tags=25%, list=8%, signal=27%
15	Arachidonic Acid	78	0.524402	1.64341	0.004115226	0.212851	0.518	998	tags=38%, list=17%, signal=46%
16	CLOFOP [ISO] (2-[4-(4-chlorophenoxy)phenoxy]propanoic acid)	16	0.667372	1.550415	0.02586207	0.2191059	0.681	1328	tags=69%, list=22%, signal=88%
17	Luteolin	44	0.54766	1.634669	0.024844721	0.2199198	0.532	1012	tags=43%, list=17%, signal=52%
18	Ammonium Hexachloroplatinate (IV)	16	0.600401	1.550776	0.040733196	0.2241297	0.681	648	tags=38%, list=11%, signal=42%
19	Bisindolylmaleimide I	23	0.587615	1.552711	0.030991735	0.2265287	0.679	964	tags=39%, list=16%, signal=46%
20	Phytoestrogens	16	0.768827	1.667617	0.022	0.2282088	0.462	856	tags=63%, list=14%, signal=73%
21	Thapsigargin (67526-95-8)	327	0.333608	1.533748	0.002061856	0.2255021	0.711	1303	tags=34%, list=22%, signal=41%
22	1-(5-deoxypentofuranosyl)-5-fluoropyrimidine-2,4 (1 h, 3 h)-dione	23	0.565554	1.552936	0.0385439	0.2324028	0.679	1330	tags=52%, list=22%, signal=67%
23	Nelfinavir	18	0.513173	1.529705	0.029661017	0.2336185	0.716	1169	tags=56%, list=19%, signal=69%
24	Vinblastine	75	0.419995	1.523846	0.036885247	0.234146	0.724	995	tags=31%, list=17%, signal=36%
25	Thimerosal	53	0.525926	1.601744	0.018867925	0.2357184	0.592	1071	tags=47%, list=18%, signal=57%
26	MG-132 (N-Benzyloxycarbonyl-L-leucyl-L-leucyl-L-leucinal)	67	0.455003	1.554736	0.017391304	0.2358129	0.675	1012	tags=39%, list=17%, signal=46%
27	Thalidomide	61	0.512241	1.525544	0.025641026	0.2359595	0.719	1101	tags=48%, list=18%, signal=58%
28	Protoporphyrin IX	22	0.540459	1.536768	0.048625793	0.2360099	0.704	1355	tags=50%, list=23%, signal=64%
29	Phosphine	37	0.573179	1.534214	0.038793102	0.2361448	0.71	1479	tags=51%, list=25%, signal=68%
30	Adenine	25	0.575044	1.752054	0.007968128	0.2367139	0.286	650	tags=40%, list=11%, signal=45%
31	Caffeic Acid	24	0.516087	1.539089	0.032388665	0.2367951	0.7	1803	tags=67%, list=30%, signal=95%
32	Lucanthone	71	0.547325	1.609604	0.046653144	0.2390304	0.571	742	tags=44%, list=12%, signal=49%
33	Dronabinol	134	0.414731	1.603898	0.008230452	0.241325	0.586	1067	tags=34%, list=18%, signal=40%
34	Antimony Potassium Tartrate	37	0.50915	1.613491	0.012711864	0.24176	0.561	1181	tags=49%, list=20%, signal=60%
35	Rottlerin	35	0.511257	1.554947	0.020408163	0.2423282	0.674	1221	tags=40%, list=20%, signal=50%
36	DMNQ (2,3-Dimethoxy-1,4-naphthoquinone)	62	0.493426	1.592946	0.022916667	0.2429922	0.604	874	tags=44%, list=15%, signal=50%
37	Rapamycin	89	0.435861	1.58848	0.012765957	0.2433278	0.61	1346	tags=45%, list=22%, signal=57%
38	Gefitinib	33	0.622531	1.769392	0.006276151	0.2467885	0.242	753	tags=42%, list=13%, signal=48%
39	Atenolol	17	0.598956	1.667728	0.006012024	0.2487788	0.462	120	tags=24%, list=2%, signal=24%
40	Antimony	35	0.540104	1.615625	0.02096436	0.2491837	0.556	1181	tags=51%, list=20%, signal=64%
41	Acrolein	43	0.538176	1.555005	0.0392562	0.2494786	0.674	1067	tags=47%, list=18%, signal=56%

Notes:  Column Headings as per GSEA website (https://software.broadinstitute.org/gsea/): Size; Number of genes in the gene set. Enrichment score for the gene set; the degree to which this gene set is overrepresented at the top or bottom of the ranked list of genes in the expression dataset. Normalized enrichment score; the enrichment score for the gene set after it has been normalized across analyzed gene sets. Nominal p-value; the statistical significance of the enrichment score. The nominal p-value is not adjusted for gene set size or multiple hypothesis testing; therefore, it is of limited use in comparing gene sets. False discovery rate; the estimated probability that the normalized enrichment score represents a false-positive finding. Familywise-error rate; a more conservatively estimated probability that the normalized enrichment score represents a false-positive finding. Rank at Max; the position in the ranked list at which the maximum enrichment score occurred. Three statistics are used to define the leading edge subset. Tags; the percentage of gene hits before (for positive ES) or after (for negative ES) the peak in the running enrichment score. This indicates the percentage of genes contributing to the enrichment score. List; the percentage of genes in the ranked gene list before (for positive ES) or after (for negative ES) the peak in the running enrichment score. This indicates where in the list, the enrichment score is attained. Signal, the enrichment signal strength that combines the two previous statistics.

Table 3 The top enriched chemicals by GSEA categorized into different subgroups

Occupational	Drugs	Plant/Plant toxins/Food
1. Ammonium Hexachloroplatinate (IV)	1. Nelfinavir	1. Adenine
2. Phosphine	2. Thalidomide	2. Arachidonic Acid
3. Acrylamide	3. Antimony Potassium Tartrate	3. Baicalein
Pesticide/herbicide	4. 4-Hydroxytamoxifen	4. Caffeic Acid
5. Acrolein	5. SB 202190 (4-[4-(4-fluorophenyl)-5-(4-pyridinyl)-1H-imidazol-2-yl]-phenol)	6. Corilagin
7. CLOFOP [ISO] (2-[4-(4-chlorophenoxy)phenoxy]propanoic acid)	6. Myricetin	8. Curcumin
Chemical compounds	7. Lucanthone	9. Ellagic Acid
10. Bisindolylmaleimide I	8. Dronabinol	11. Luteolin
12. Thapsigargin (67526-95-8)	9. Rapamycin	13. Microcystin RR
14. MG-132 (N-Benzyloxycarbonyl-L-leucyl-L-leucyl-L-leucinal)	10. Atenolol	15. Phytoestrogens
16. DMNQ (2,3-Dimethoxy-1,4-naphthoquinone)	Chemotherapy	17. Protoporphyrin IX
18. Tyrphostin AG 538	19. Vinblastine	20. Purpurogallin
21. 1-(5-Deoxypentofuranosyl)-5-Fluoropyrimidine-2,4 (1 h, 3 h)-Dione	22. 7-Amino-4-Hydroxy-2-Naphthalenesulfonic Acid	23. Rottlerin
24. Antimony	25. Gefitinib	26. Salazinic Acid

The identified chemicals share exciting immune/inflammation-related pathways

In order to examine which pathways are associated with the largest number of the identified 41 chemicals, we used the Comparative Toxicogenomics Database (CTD) batch query webtool (http://ctdbase.org/tools/batchQuery). The tool can generate a report listing the pathways that show significant association with the given chemical, thus having a potentially significant effect on a proportion of the genes of that pathway. More than 70% of the 41 identified chemicals are associated with common pathways mainly involved in the immune response, as shown in Table 4. Those pathways are: Immune system, Cytokine signaling in immune system, IL-17 signaling pathway, Pathways in cancer, Signaling by interleukins, Innate immune system, Apoptosis, TNF signaling pathway, Cellular responses to stress, Toll-like receptor signaling pathway, Influenza A, Adaptive immune system, Downstream signaling events of B Cell Receptor (BCR), Senescence-Associated Secretory Phenotype (SASP), Fc epsilon receptor (FCERI) signaling, Signaling by EGFR, Th17 cell differentiation, Toll-Like receptors cascades, Cellular senescence, Interleukin-10 signaling, Activated TLR4 signaling.

Table 4 List of pathways significantly associated with the largest number of the identified 41 chemicals using Comparative Toxicogenomics Database (CTD) batch query webtool (http://ctdbase.org/tools/batchQuery). Only pathways that are shared by more than 50% of the identified chemicals were listed

Significant Chemicals Associated Pathways	Shared by how many chemicals (Total=41)	Percentage
Immune System	33	80%
Cytokine Signaling in the Immune system	31	76%
IL-17 signaling pathway	31	76%
Pathways in cancer	31	76%
Signaling by Interleukins	31	76%
Signal Transduction	31	76%
HTLV-I infection	30	73%
Interleukin-4 and 13 signaling	30	73%
Chagas disease (American trypanosomiasis)	29	71%
Fluid shear stress and atherosclerosis	29	71%
Hepatitis B	29	71%
Innate Immune System	29	71%
Metabolism	29	71%
Apoptosis	28	68%
Pertussis	28	68%
TNF signaling pathway	28	68%
Toxoplasmosis	28	68%
Tuberculosis	28	68%
AGE-RAGE signaling pathway in diabetic complications	27	66%
Endocrine resistance	27	66%
Gene Expression	27	66%
MAPK signaling pathway	27	66%
PI3K-Akt signaling pathway	27	66%
Signaling by NGF	27	66%
Viral carcinogenesis	27	66%
Cellular responses to stress	27	66%
Hemostasis	27	66%
Herpes simplex infection	27	66%
Non-alcoholic fatty liver disease (NAFLD)	27	66%
Platinum drug resistance	27	66%
MicroRNAs in cancer	26	63%
NOD-like receptor signaling pathway	26	63%
Proteoglycans in cancer	26	63%
Toll-like receptor signaling pathway	26	63%
Amoebiasis	26	63%
HIF-1 signaling pathway	26	63%
Influenza A	26	63%
Legionellosis	26	63%
Progesterone-mediated oocyte maturation	26	63%
Prostate cancer	26	63%
Amyotrophic lateral sclerosis (ALS)	25	61%
Adaptive Immune System	25	61%
Bladder cancer	25	61%
Cell cycle	25	61%
Colorectal cancer	25	61%
Downstream signaling events of B Cell Receptor (BCR)	25	61%
Generic Transcription Pathway	25	61%
Leishmaniasis	25	61%
Senescence-Associated Secretory Phenotype (SASP)	25	61%
FoxO signaling pathway	25	61%
Measles	25	61%
p53 signaling pathway	25	61%
Rheumatoid arthritis	25	61%
Developmental Biology	24	59%
Downstream signal transduction	24	59%
Epstein-Barr virus infection	24	59%
Estrogen signaling pathway	24	59%
Fc epsilon receptor (FCERI) signaling	24	59%
Metabolism of lipids and lipoproteins	24	59%
NGF signaling via TRKA from the plasma membrane	24	59%
Osteoclast differentiation	24	59%
Signaling by EGFR	24	59%
Signaling by PDGF	24	59%
Th17 cell differentiation	24	59%
Toll-Like Receptors Cascades	24	59%
Cellular Senescence	24	59%
Insulin resistance	24	59%
Interleukin-10 signaling	24	59%
Transcriptional misregulation in cancer	24	59%
Transcriptional Regulation by TP53	24	59%
Cell Cycle, Mitotic	24	59%
Activated TLR4 signalling	23	56%
Breast cancer	23	56%
Central carbon metabolism in cancer	23	56%
DAP12 interactions	23	56%
DAP12 signaling	23	56%
MyD88 cascade initiated on the plasma membrane	23	56%
MyD88 dependent cascade initiated on endosome	23	56%
MyD88-independent TLR3/TLR4 cascade	23	56%
MyD88:Mal cascade initiated on plasma membrane	23	56%
Neurotrophin signaling pathway	23	56%
Prolactin signaling pathway	23	56%
Salmonella infection	23	56%
Signaling by SCF-KIT	23	56%
Th1 and Th2 cell differentiation	23	56%
Toll Like Receptor 10 (TLR10) Cascade	23	56%
Toll-Like Receptor 2 (TLR2) Cascade	23	56%
Toll-Like Receptor 3 (TLR3) Cascade	23	56%
Toll-Like Receptor 5 (TLR5) Cascade	23	56%
Toll-Like Receptor 7/8 (TLR7/8) Cascade	23	56%
Toll-Like Receptor 9 (TLR9) Cascade	23	56%
Toll-Like Receptor TLR1: TLR2 Cascade	23	56%
Toll-Like Receptor TLR6: TLR2 Cascade	23	56%
TRAF6 mediated induction of NFkB and MAP kinases upon TLR7/8 or 9 activation	23	56%
TRIF-mediated TLR3/TLR4 signaling	23	56%
Chronic myeloid leukemia	23	56%
Cytokine-cytokine receptor interaction	23	56%
Disease	23	56%
Hepatitis C	23	56%
Inflammatory bowel disease (IBD)	23	56%
Pancreatic cancer	23	56%
Signaling by the B Cell Receptor (BCR)	23	56%
Alzheimer’s disease	23	56%
Cell Cycle	23	56%
ErbB signaling pathway	22	54%
FCERI mediated MAPK activation	22	54%
GAB1 signalosome	22	54%
Metabolism of proteins	22	54%
PI3K/AKT activation	22	54%
PIP3 activates AKT signaling	22	54%
Shigellosis	22	54%
Signaling by VEGF	22	54%
T cell receptor signaling pathway	22	54%
Toll-Like Receptor 4 (TLR4) Cascade	22	54%
VEGF signaling pathway	22	54%
Apoptosis	22	54%
Diseases of signal transduction	22	54%
Epithelial cell signaling in Helicobacter pylori infection	22	54%
Intrinsic Pathway for Apoptosis	22	54%
Malaria	22	54%
Mitotic G1-G1/S phases	22	54%
NF-kappa B signaling pathway	22	54%
Programmed Cell Death	22	54%
Ras signaling pathway	22	54%
Small cell lung cancer	22	54%
Prion diseases	22	54%
Activation of the AP-1 family of transcription factors	21	51%
Autophagy - animal	21	51%
B cell receptor signaling pathway	21	51%
cAMP signaling pathway	21	51%
Endometrial cancer	21	51%
MAPK family signaling cascades	21	51%
Thyroid hormone signaling pathway	21	51%
VEGFA-VEGFR2 Pathway	21	51%
Acute myeloid leukemia	21	51%
Adipocytokine signaling pathway	21	51%
Apoptosis - multiple species	21	51%
Chemokine signaling pathway	21	51%
EGFR tyrosine kinase inhibitor resistance	21	51%
Glioma	21	51%
Jak-STAT signaling pathway	21	51%
Longevity regulating pathway	21	51%
Oxidative Stress-Induced Senescence	21	51%
Renal cell carcinoma	21	51%
Role of LAT2/NTAL/LAB on calcium mobilization	21	51%
Sphingolipid signaling pathway	21	51%
Platelet activation, signaling, and aggregation	21	51%

22 of the identified chemicals are enriched for NF-kB pathways

Notably, 22 out of the total 41 identified chemicals showed significant enrichment for the NF-κB pathway and those were ellagic acid, baicalein, curcumin, SB 202190, acrylamide, myricetin, arachidonic acid, luteolin, ammonium hexachloroplatinate(iv), bisindolylmaleimide I, 67526-95-8, vinblastine, thimerosal, MG-132, thalidomide, adenine, caffeic acid, dronabinol, rottlerin, rapamycin, gefitinib and acrolein.

The transcription factor NF-κB regulates innate and adaptive immune functions through upregulation of pro-inflammatory genes, and when deregulated, it can contribute to the pathogenic processes of various inflammatory diseases.²⁵ NF-κB was shown to be activated predominantly in the epithelial cells of the conducting airways, which have been reported to be the main source of NF-κB-dependent mediators that play a role in asthma.²⁶ Any inhalational stimuli can activate bronchial epithelial NF-κB pathway sufficiently to promote allergic sensitization to innocuous inhaled antigens.²⁷ The mainstay of therapy for asthma is the anti-inflammatory glucocorticoids that act mainly by inhibiting NF-κB induced gene transcription.²⁸ These reports indicate the central role of the NF-κB pathway in asthma pathogenesis and hence propose it as an important therapeutic target.²⁹

Plants, plant toxins and food-related asthma triggers

In this study, we focused on the identified plants, plant toxins, and food-related chemicals. This is due to the fact that there is no specific risk assessment for their potential effect on asthma development or exacerbation, even though asthmatic patients are in contact with one or more of these triggers in their close environment. Drugs, chemotherapy, chemical compounds, and occupational hazards are usually associated with a specific warning and awareness regarding asthma, although the exact underlying mechanisms are not fully understood.

Food

DNA, present in food, can survive harsh processing and be absorbed to circulate through the blood to other tissues of human and animals.³⁰ Dietary purines like adenine are found in virtually all foods.³¹ Adenine and guanine comprise more than 60% of total purine-rich foods (such as cereals, beans, soybean products, and seaweeds),³¹ with greater bioavailability of adenine than guanine.³² There are many circumstantial pieces of evidence that purine and its metabolites might have a role in asthma with no conclusive findings. Allergic asthmatic plasma metabolomics showed aberrant purine metabolism that may change the consequence of having a more purine-rich diet in such patients.³³ It has been previously reported that allergy to purine-rich wheat flour is the leading cause of serious occupational asthma among bakery workers called baker’s asthma.³⁴ Another possible indirect link between purine-rich diet and asthma is through gout, sleep apnea, and circadian rhythm. A purine-rich diet is associated with a high risk of gout,³⁵ which in turn is linked to sleep apnea.³⁶ Unrecognized obstructive sleep apnea (OSA) can potentiate poor asthma control despite optimal therapy.³⁷ The endogenous circadian system prolongs respiratory events across the night and can modulate sleep apnea.³⁸ Circadian regulation of de novo purine synthesis is an important mechanism conferring circadian rhythmicity on the cell cycle.³⁹ Intriguingly, our results showed that genes affected by dietary adenine were differentially expressed in severe asthmatic bronchial epithelium. Most of these genes are lung epithelial cell tissue-specific genes (BLM; CHAF1A; GDF15; KIF20A; CDC25A). Of interest, one of these genes (PRKAB1) is involved in circadian rhythmicity.⁴⁰ Furthermore, one of the characteristics of asthma is worsening of symptoms overnight that has been linked to circadian variations controlled by clock genes.⁴¹ Therefore, activation of circadian rhythm genes by purine-rich food can be the link between gout, sleep apnea, and asthma.³⁸

Another interesting food-related chemical identified by our method is Arachidonic acid(ARA), an omega-6 polyunsaturated fatty acid found in the phospholipids of the cell membranes and is abundant in the brain, muscles, and liver.⁴² Arachidonic acid occurs in the animal source diet such as eggs, poultry, and meat.⁴³ This could explain and suggest a possible contributing factor to the increased incidence of asthma in western societies due to their consumption of such a pro-inflammatory diet and thus promoting the release of pro-inflammatory arachidonic acid metabolites (leukotrienes and prostanoids).⁴⁴ Prostaglandins and leukotrienes are arachidonic acid-derived lipid mediators converted via cyclooxygenase and lipoxygenase, respectively and play a major role in asthma.⁴⁵ However, there are insufficient studies to draw any firm conclusions about the relationship between ARA and asthma risk.⁴⁶ Surprisingly, our results showed that genes targeted by arachidonic acid are specific to alveolar macrophages (ABCA1; JUN; SERPINB2; HPGD; IL1B; PLA2G4A; PPARG; FOS; PTGS2; PHLDA1; PLA2G7; ATF3). Alveolar macrophages serve as the first line of defense against foreign invaders to the lung tissue and have a critical role in asthma.⁴⁷ Unlike blood-borne monocytes, resident alveolar macrophages have a suppressive role to inflammation but could gain pathogenic functions after repeated exposures.⁴⁸

Strawberries are considered as functional food and nutraceutical source, mainly because of their high concentration of ellagic acid (EA) and its precursors.⁴⁹ EA is derived from ellagitannins (ETs) and is found in some nuts, seeds, and fruits, especially berries and fruit juices.⁵⁰ Dietary ETs are partially hydrolyzed in the gut to EA then to urolithin A (UA) by colonic microflora to enter the circulation.⁵¹ ETs are natural polyphenolic compounds that show potent anti-inflammatory properties in various diseases such as that observed in OVA-induced asthma mouse model, possibly through inhibition of NF-κB activation.⁵² Furthermore, EA has an anti-eosinophilic activity in a murine model of asthma⁵³ and was suggested as a potential therapeutic agent for accelerating the resolution of allergic airways inflammation.⁵⁴

Another identified food component is the flavone luteolin, which is found in several plant products, including broccoli, pepper, thyme, and celery. Due to its anti-inflammatory and neuroprotective function, plants rich in luteolin have been used in Chinese traditional medicine for treating various diseases such as hypertension, inflammatory disorders, and cancer.⁵⁵ Through intrinsic and extrinsic signaling pathways, luteolin as an active compound showed anti-oxidant, anti-tumor, anti-inflammatory, and anti-apoptotic activities.⁵⁶ Our results showed that targeted genes by luteolin in asthmatic epithelium are related to inflammation pathways like TNF signaling pathway (NFKBIA; JUN; IL1B; FOS; PTGS2; JUNB), Th17 cell differentiation and IL-17 signaling pathway (NFKBIA; JUN; IL1B; FOS, FOSB; PTGS2), Arachidonic acid metabolism (PTGS2; CBR3), Toll-like receptor signaling pathway (NFKBIA;JUN; IL1B; FOS) and NF-κB signaling pathway (NFKBIA; GADD45B; IL1B; PTGS2).

Caffeic acid is an active anti-oxidative component, that has been shown to have beneficial effects on several respiratory disorders, such as chronic obstructive pulmonary disease and lung cancer.⁵⁷ Caffeic acid has powerful antimicrobial, antioxidant activities, and can influence collagen production and block premature aging.⁵⁸ It was shown previously that caffeic acid has immunoregulatory effects by inhibition of cytokine and chemokine production as well as enhancement of transforming growth factor-beta 1 production in asthmatics.^59,60 Of interest, another enriched pathway between severe asthmatic and healthy epithelium and related to coffee is Pyrogallol which is converted under alkaline conditions into purpurogallin,⁶¹ generating reactive oxygen species. Another source for purpurogallin is Galls, the abnormal growth in plants. Galls are induced by viruses, bacteria, fungi, nematodes, arthropods, or even other plants, which are similar to cancers in fauna and used in folkloric medicine.⁶² Purpurogallin was shown to exert antiplatelet, antithrombotic⁶³ and anti‑inflammatory effects by inhibiting NF‑κB and MAPK signaling pathways in lipopolysaccharide‑stimulated cells.⁶⁴ Due to these anti-inflammatory activities, it was suggested to be a therapeutic target for various systemic inflammatory diseases.⁶⁵ Our results showed that genes affected by purpurogallin and upregulated in severe asthmatic epithelium (TNNI3, HPGD, UHRF1, TNNC1, SELL, BLM, GFER, PLA2G7, TDP1, CDK5, SENP8, MCL1, GAPDH, RNASEH1, GSK3A, RUNX1, PABPC1, BCL2L1) are related to positive regulation of apoptotic process, leukocyte cell-cell adhesion and DNA repair.

Protoporphyrin IX (PPIX) is a heterocyclic organic compound, which consists of four pyrrole rings, and is the final intermediate in the heme biosynthetic pathway. It is ubiquitously present in all living cells in small amounts.⁶⁶ PPIX is a naturally occurring pigment in meat products that is increased by higher pH conditions in the context of nitrite reduction.⁶⁷ Also, PPIX is the main pigment resulting in the brown coloration of eggshells.⁶⁸ PPIX can induce heme oxygenase which was shown to inhibit Th17 cell-mediated immune response and prevent ovalbumin-induced neutrophilic airway inflammation.⁶⁹

Plants

Baicalein (5,6,7 trihydroxyflavone) is a famous phenolic flavonoid present in the dry roots of Scutellaria baicalensis plant and is a component of the traditional herbal remedy known as Chinese skullcap (or Huang Qin).⁷⁰ It was shown to attenuate inflammatory responses by suppressing TLR4 mediated NF-κB and MAPK signaling pathways.⁷¹ Furthermore, baicalein protects cells from hydrogen peroxide by inhibiting 12-lipoxygenase thus blocking the increase in ROS levels.⁷² Our results have shown that genes (HPGD, SELL, BLM, CYP2D6, PLA2G7, TDP1, GAPDH) that are affected by baicalein are significantly enriched in asthmatic bronchial epithelium. HPGD (15-Hydroxyprostaglandin Dehydrogenase) gene contributes to the regulation of events that are under the control of prostaglandin levels, and its expression is affected by aspirin.⁷³ This can be part of the Aspirin-Exacerbated Respiratory Disease (AERD) which is a syndrome that includes asthma, recurrent nasal polyps, and pathognomonic reactions to aspirin and other nonselective cyclooxygenase inhibitors.⁷⁴ On the other side, SELL (selectin S) gene was previously shown to be upregulated in different lung inflammatory diseases.⁷⁵

Corilagin is one of the major active components of many ethnopharmacological plants isolated from Caesalpinia and was reported to exhibit anti-tumor and anti-inflammatory activities.⁷⁶ Corilagin was shown to inhibit the release of cytokines such as TNF-α, IL-1β, and IL-6 as well as the production of nitric oxide.⁷⁷ More specifically, corilagin possess anti-anaphylactic and anti-allergic activities by inhibiting the release of mediators from mast cells and by decreasing the serum concentration of immunoglobulin E (IgE).⁷⁸ Furthermore, the potent inhibition of the Corilagin on the phagocytic activity of neutrophils makes it an interesting herbal asthma remedy.⁷⁹ Our results showed that the upregulated genes in the asthmatic epithelium and are part of Corilagin targets include TNNI3, HPGD, TNNC1, BLM, GFER, PLA2G7, SQLE, MCL1, PPARG. Four of them are alveolar macrophage-specific genes (HPGD; PPARG; PLA2G7; MCL1).

Phytoestrogens are plant-derived compounds found in a wide variety of foods and plants, being most abundant in soy,⁸⁰ which is known to induce allergy, affecting approximately 0.4% of children.⁸¹ This could be the reason why it is discouraged to use soya protein in children in the first six months of life to avoid sensitization and exposure to phytoestrogens.⁸² Soy sauce is a traditional fermented seasoning of Japan, that is made from soybeans and wheat, both of which are established food allergens.⁸⁰ On the other hand, phytoestrogens were reported to have a protective role against heart disease, breast cancer, and menopausal symptoms of osteoporosis.⁸³ Phytoestrogens have structural similarities to estrogen and hence can bind to its receptor causing (anti)oestrogenic effects⁸⁴ and could cause potential adverse health effects as well.⁸³ With regards to asthma, it was found that increased consumption of phytoestrogens may help prevent or treat asthma and allergic disease.⁸⁵ Furthermore, phytoestrogens can reduce antigen-induced eosinophilia in the lung.⁸⁶ It was not surprising that the genes affected by phytoestrogens and upregulated in the severe asthmatic epithelium in our analysis (TOP2A, ANLN, MKI67, UHRF1, BIRC5, TFF1, MND1, RRM2, VWF, NUSAP1), are related to cell nuclear division, mitotic nuclear division, and cell cycle. TOP2A, ANLN, RRM2, UHRF1, NUSAP1, BIRC5, MND1, and MKI67 are enriched specifically in bronchial epithelial cells.

Rottlerin, also called mallotoxin, is the principal phloroglucinol constituent of the Mallotus Philippinensis (known as Kamala Tree). Previous studies have shown that rottlerin induces apoptosis, autophagy, and suppresses NF-κB and PKCδ in cancer cells, such as lung cancer.^87,88 Rottlerin was shown to inhibit microvascular endothelial cells tube formation, block cell senescence, and intracellular ROS generation in psoriasis.⁸⁹ In the lung, rottlerin was shown to be anti-inflammatory, airways smooth muscles relaxant⁹⁰ and suppressant of airway hyperreactivity in mouse models of experimental asthma.⁹¹ Additionally, rottlerin induces apoptosis of human blood eosinophils , hence, can attenuate allergic reactions.⁹² On the other hand, rottlerin increases barrier dysfunction in pulmonary endothelial cell monolayers and causes pulmonary edema in rats.⁹³

Salazinic acid can be isolated from Xanthoparmelia camtschadalis, Rimelia cetrata, and Parmelia caperata. It can be used as an antioxidant agent which plays an important role in macrophage killing of bacteria and tumors.⁹⁴ Beside the anti-oxidative effect, salazinic acid has immunostimulatory, antimicrobial and antiproliferative potentials.^94,95 Lichens contain large amounts of salazinic acid and have been used since ancient times as a therapeutic agent for the treatment of bronchitis, asthma, and inflammation.⁹⁶ In our analysis, genes targeted by salazinic acid and upregulated in severe asthmatic bronchial epithelium (TNNI3, TNNC1, GFER, PLA2G7, TDP1, MCL1, PIN1, RUNX1, PABPC1, BCL2L1, RGS12) are related to the regulation of neuron apoptotic process and cardiac muscle tissue development.

Plant toxins

Microcystins (MCs) are hepatotoxins, produced by various species of cyanobacteria, whose occurrence is increasing worldwide owing to climate change and anthropogenic activities.⁹⁷ Various edible aquatic organisms, plants, and food supplements based on algae can bioaccumulate these toxins.⁹⁸ This occurs at times when blooms form and accumulate as scum on the water surface after which the death and decay of cells release large amounts of cyanotoxins which become toxic to eukaryotic organisms, including humans.⁹⁹ Contact dermatitis, asthma-like symptoms, and symptoms resembling hay fever have been attributed to microcystins chemical sensitivity.¹⁰⁰ Hence, water-based recreational activities can expose people to very low concentrations of aerosol-borne microcystins¹⁰¹ or even aerosolized cyanotoxins, making inhalation a potential route of exposure.¹⁰² Exposure to such aerosolized toxins in asthmatic subjects can have adverse effects.³ Our results showed that targeted genes by Microcystin in the asthmatic epithelium are specific to the trachea (ZNF57, LMNA, SERPINB5) and that the gene GSTT1 showed significant association with asthma risk.¹⁰³

Conclusion

Our analysis using the publicly available gene expression data and linking it to toxicological omics’ data was able to explain and predict the toxicity in terms of affwcting the differentially expressed genes between severe asthmatic and normal epithelium. Many of the identified chemicals using this approach have no special warnings or precautions to avoid them by asthma patients. Even if some of the identified genes were reported earlier and linked to asthma, the exact mechanism is still poorly understood. The enriched pathways shared by most of the chemicals identified were related to significant players in the signaling pathways that are associated with triggering or exacerbation of asthma development. Such an approach can pave the way to generate a cost-effective and reliable source for asthma-specific toxigenic reports thus allowing the asthmatic patients, physicians, and medical researchers to be aware of the potential triggering factors with fatal consequences.

Acknowledgment

R.H. is funded by the Sharjah Research Academy (Grant code: MED001), University of Sharjah (Grant code: 1901090254) and Al-Jalila Foundation (Grant code: AJF201741).

Disclosure

The authors report no conflicts of interest in this work.