Endocrine disrupting chemicals and COVID-19 relationships: A computational systems biology approach.

BACKGROUND: Patients at high risk of severe forms of COVID-19 frequently suffer from chronic diseases, but other risk factors may also play a role. Environmental stressors, such as endocrine disrupting chemicals (EDCs), can contribute to certain chronic diseases and might aggravate the course of COVID-19.
OBJECTIVES: To explore putative links between EDCs and COVID-19 severity, an integrative systems biology approach was constructed and applied.
METHODS: As a first step, relevant data sets were compiled from major data sources. Biological associations of major EDCs to proteins were extracted from the CompTox database. Associations between proteins and diseases known as important COVID-19 comorbidities were obtained from the GeneCards and DisGeNET databases. Based on these data, we developed a tripartite network (EDCs-proteins-diseases) and used it to identify proteins overlapping between the EDCs and the diseases. Signaling pathways for common proteins were then investigated by over-representation analysis.
RESULTS: We found several statistically significant pathways that may be dysregulated by EDCs and that may also be involved in COVID-19 severity. The Th17 and the AGE/RAGE signaling pathways were particularly promising.
CONCLUSIONS: Pathways were identified as possible targets of EDCs and as contributors to COVID-19 severity, thereby highlighting possible links between exposure to environmental chemicals and disease development. This study also documents the application of computational systems biology methods as a relevant approach to increase the understanding of molecular mechanisms linking EDCs and human diseases, thereby contributing to toxicology prediction.

Introduction

The COVID-19 pandemic started in the fall of 2019 and spread to a large part of the world during the winter and spring of 2020. By late September 2020, it had led to more than a million deaths, of which one-fifth in the US and somewhat fewer in the EU ( , ). Despite considerable research activities, there are still many unknowns concerning this infectious disease, especially with regard to the substantial variability of the disease severity. Following an initial infectious phase, a “cytokine storm”, leading to pneumonia is observed in severe cases which may require intensive care. It is still unclear why infections lead to severe cases in some patients and not in others, but both endogenous and exogenous factors can likely influence the outcome of the disease.

In addition to older age and male sex, several comorbidities are associated with severe COVID-19 and increased mortality risk. Disorders such as cardiovascular disease, type II diabetes (T2D), obesity, chronic respiratory disease or hypertension are strongly linked to severe COVID-19 cases (Petrilli et al., 2020, Zhou et al., 2020, Stefan et al., 2020). As has recently been proposed, underlying metabolic and endocrine dysfunctions may be mechanistically linked to the exacerbation of the coronavirus infection (Bornstein et al., 2020), and these observations may inspire new insight into the pathogenesis of this disease, including biological interpretation of the mechanisms involved. Environmental stressors have already been suggested to contribute to the severity of the disease (Bashir et al., 2020, Fattorini and Regoli, 2020, Zhu et al., 2020), but little mechanistic support for this association is available. A relevant approach would be to compare the biological pathways triggered by environmental stressors with those involved in the COVID-19 severity. If similar pathways are found, this would increase the likelihood that such stressors may contribute to critical stages of this disease.

Given the suspected hormonal mode of vulnerability (Drucker, 2020) endocrine disrupting chemicals (EDCs) could represent important triggers of aggravated infection, e.g., in the form of phthalates, bisphenols, organochlorine pesticides, and perfluorinated alkane substances (PFASs) (Trasande et al., 2016, Vandenberg et al., 2016). Exposure to these substances may affect the immune defense, thus potentially increasing the susceptibility to develop COVID-19 (Tsatsakis et al., 2020), as supported by experimental studies (Cipelli et al., 2014, Couleau et al., 2015). For example, epidemiological evidence on children exposed to PFASs show decreased immune responses to routine vaccines (Grandjean et al., 2012) and a greater risk of developing infectious disease (Dalsager et al., 2016, Granum et al., 2013).

As promising tools to gain better insight into the possible risk factors and mechanisms, toxicological and chemical data sources have expanded substantially, thereby enabling network science and computational systems biology methods to become feasible (Audouze et al., 2013, Audouze et al., 2018, Taboureau and Audouze, 2017, Vermeulen et al., 2020, Wu et al., 2020). We have therefore conducted an integrative systems biology exploration to identify overlapping proteins that are both dysregulated by EDCs and involved in comorbidities associated with aggravated COVID-19. Based on this tripartite network, integrating protein-EDC associations and protein-disease annotations, we then performed biological enrichments of pathways to detect the most plausible relationships between EDC exposure and COVID-19 severity.

Material and methods

We employed a computational systems biology approach to explore putative linkages between EDCs and COVID-19 as presented in Fig. 1 . First, a tripartite network was created based on known associations between proteins and either COVID-19 comorbidities or EDCs, as compiled from existing databases (CompTox, DisGeNET, GeneCards) (A). Then, biological enrichment was performed with the jointly identified proteins (i.e., those retrieved in both association studies) (B) by over-representation analysis (ORA) to identify the pathways that were the highly linked to both the diseases and the EDCs (C). As a final step, the biological pathways were explored with available knowledge regarding COVID-19 mechanisms (from the literature and the AOP-Wiki database), thereby allowing consideration of hypothetical linkages between EDCs and COVID-19 (D).

Fig. 1

Overview of the integrative systems toxicology approach. A: Human proteins known to be dysregulated by endocrine-disrupting chemicals (EDCs) were extracted from the CompTox database; human proteins linked to obesity or to comorbidities or metabolic dysfunction known to be associated with obesity were compiled using DisGeNET and GeneCards. These compiled data were used to develop a tripartite network. B: A set of proteins was identified that was common to both association studies (proteins targeted by the EDCs and also involved in comorbidities). C: Biological enrichment was performed for pathways for each of the four databases, by over-representation analysis (ORA) to identify potential mechanisms of action related to these proteins, where the biological pathways were ranked by their statistical significance. D: The most relevant of the potential pathways were compared to known COVID-19 dysregulated pathways from the literature and the AOP-Wiki database.

Endocrine-disrupting chemical dataset

A list of commonly used substances known or suspected to act as EDCs was established, based on knowledge from three data sources: the endocrine disruptor assessment list from ECHA (https://echa.europa.eu/fr/ed-assessment, as of April 24, 2020), the list from NIEHS (https://www.niehs.nih.gov/health/topics/agents/endocrine/index.cfm as of April 28, 2020), and the TEDX list (https://endocrinedisruption.org/interactive-tools/tedx-list-of-potential-endocrine-disruptors/search-the-tedx-list, as of April 24, 2020).

To explore as much as possible the chemical diversities, 34 EDCs chosen for this study were selected to represent different chemical classes (Table 1 ). The CAS numbers were used for data integration. The substances selected are known to activate different receptors, and their various impacts have been shown in epidemiological studies. Still, these EDCs do not necessarily represent the most abundant chemicals in human matrices, and bioaccumulation and biotransformation may affect the degree of human exposures and their effects.

Table 1

List of the 34 major substances known or suspected to be endocrine-disrupting chemicals.

CAS	Chemical name	Abbreviation	Chemical structure
35065-27-1	2,2′,4,4′,5,5′-hexachlorobiphenyl	PCB 153
1746-01-6	2,3,7,8- tetrachlorodibenzodioxin	TCDD
1912-24-9	atrazine	-
131-56-6	benzophenone-1	-
117-81-7	bis (2-ethylhexyl)phthalate	DEHP
620-92-8	bisphenol F	BPF
80-05-7	bisphenol A	BPA
80-09-1	bisphenol S	BPS
94-26-8	butyl-paraben	BUPA
57-74-9	chlordane	–
2921-88-2	chlorpyrifos	–
210880-92-5	clothianidin	–
52315-07-8	cypermethrin	–
486-66-8	daidzein	–
84-74-2	dibutyl phthalate	DBP
72-55-9	Dichlorodiphenyldichloroethylene	DDE
50-29-3	dichlorodiphenyltrichloroethane	DDT
446-72-0	genistein	–
3194-55-6	hexabromocyclododecane	HBCD
118-74-1	hexachlorobenzene	HCB
138261-41-3	imidacloprid	–
625-45-6	methoxyacetic acid	MAA
99-76-3	methyl-paraben	MEPA
68412-53-3	nonylphenol ethoxylate	NPEO	–
103-90-2	paracetamol	–
68631-49-2	PBDE-153	PBDE-153
5436-43-1	PBDE-47	PBDE-47
14797-73-0	perchlorate	–
1763-23-1	perfluorooctane sulfonic acid	PFOS
335-67-1	perfluorooctanoic acid	PFOA
67747-09-5	prochloraz	–
153719-23-4	thiamethoxam	–
688-73-3	tributyltin	TBT
3380-34-5	triclosan	–

Disease dataset

Comorbidities known to be associated with obesity or otherwise leading to severe COVID-19 were extracted from a recent study (Stefan et al., 2020), and resulted in a total of 13 disorders for exploration in the integrative systems toxicology (Table 2 ).

Table 2

List of the 13 diseases.

Obesity
Impaired respiratory mechanisms	respiratory dysfunction
Increased airway resistance	respiratory dysfunction
Impaired gas exchange	respiratory dysfunction
Low lung volume	respiratory dysfunction
Low muscle strength	respiratory dysfunction
Cardiovascular disease	comorbidities
Diabetes mellitus	comorbidities
Kidney disease	comorbidities
Hypertension	metabolic risk
Prediabetes	metabolic risk
Insulin resistance	metabolic risk
Dyslipidemia	metabolic risk

Endocrine-disrupting chemical-protein associations

Human proteins known to be associated with each of the 34 EDCs were extracted from the U.S. Environmental Protection Agency web-based CompTox Chemistry dashboard (https://comptox.epa.gov/dashboard), which contains a wide range of data related to chemical toxicity, physico-chemical properties, human exposure, and in vitro bioassay data (agonist, antagonist, up- and down- regulation) for over 87,500 chemicals (as of April 30, 2020) (Williams et al., 2017).

Each linked protein was matched to a gene symbol and classified using the Panther (protein analysis through evolutionary relationships) classification system (http://www.pantherdb.org/about.jsp) (version 15, released February 14, 2020) (Mi et al., 2013). The Panther database is a biological database of gene/protein families, and their functionally related subfamilies that can be used to classify and identify the function of gene products. This database results from a human curation and advanced bioinformatics algorithms, and the current version contains 15,702 protein families, divided into 1,239,989 functionally distinct protein subfamilies.

Disease-protein associations

From two human protein-disease databases, proteins known to be linked to the 13 studied diseases were listed (as of April 29, 2020 for both data sources). The DisGeNet database is a discovery platform containing one of the largest publicly available collections of genes and variants associated with human diseases (https://www.disgenet.org/) (Pinero et al., 2015). DisGeNet integrates data from expert curated repositories, GWAS catalogues, animal models and the scientific literature. The current version contains 1,134,942 gene-disease associations between 21,671 genes and 30,170 diseases. The GeneCards database contains manually curated information for substances and their associations to genes and proteins, that are scored (https://www.genecards.org/) (Safran et al., 2010). GeneCards integrates gene-centric data from more than 150 web sources, including genomics, transcriptomics, proteomics, clinical and functional annotations. The current version contains information for 270,168 genes, and among them 18,871 are known to be linked to human diseases. For the present study, only associations were kept only for those between human diseases and proteins categorized as coding proteins, and all non-human information, including gene clusters, genetic locus, pseudogenes, RNA genes and those uncategorized were disregarded. All listed proteins were matched to their gene symbol to facilitate further analysis. Each identified protein from both databases, was categorized into the protein class using the Panther classification (version 15).

Pathways enrichment analysis

To decipher biological pathways potentially linked to the selected EDCs and explore if they might overlap with the ones known for COVID-19, an ORA was done. Four major sources of protein-pathway information were independently integrated, i.e., using the Kyoto Encyclopedia of Genes and Genomes (KEGG), the Reactome, the Wiki-pathways and the Panther databases (Fabregat et al., 2018, Kanehisa et al., 2019, Mi et al., 2013, Slenter et al., 2018). The KEGG database is a collection of manually drawn pathways maps representing existing knowledge of the molecular interaction, reaction and relation networks for several levels of the biological systems ( ). The current version contains 537 pathway maps that include 724,592 references. The Reactome database is a manually curated and peer-reviewed source of pathway information (https://reactome.org/), which provides data related to several species. Regarding human data, 2423 pathways are currently described involving 10,923 proteins. The wiki-pathways is an open and collaborative platform dedicated to the curation of biological pathways. It has information related to 2887 pathways for several species (1200 are for human). Finally, the panther pathways database ( ) contains information for various species, among them 6976 human genes are annotated to 2608 pathways. To assess the statistical significance of the protein-pathway relationships, a hypergeometric test was used for each of the four sources, followed by a multiple testing correction of the p-values with the Benjamini-Hochberg method. The ORA was performed on the common proteins identified to identify the most strongly linked proteins that are affected by the EDCs and also associated with at least one the 13 comorbidities. As a last step, manual curation allowed us to consider relevant outcomes for interpretation. The four data sources provided complementary information, with some overlapping findings.

COVID-19 and biological mechanism of action

Linkage between COVID-19 and potential biological targets and affected pathways were extracted from the literature using the PubMed database (as of May 22, 2020) and the AOP-Wiki database (as of May 22, 2020). The PubMed database comprises more than 30 million references to biomedical literature from MEDLINE, life science journals, and online books (https://pubmed.ncbi.nlm.nih.gov. The AOP-Wiki database is part of the larger OECD AOP knowledge base that represents the central repository for all AOPs developed ( ). The current version contains 306 AOPs involving 582 stressors (chemicals).

Results

From the CompTox database, information on the links between chemicals and human proteins were compiled. Data for 30 of the 34 chemicals could be retrieved, and a total of 208 unique human proteins were involved via 1632 associations. No information was retrieved for hexachlorobenzene, nonylphenol ethoxylate, perchlorate and tributyltin. Perfluorooctane sulfonic acid (PFOS) targeted the highest number of proteins (1 1 3), and 2,2′,4,4′,5,5′-hexachlorobiphenyl (PCB 153) was associated with only one biological target (the progesterone receptor). The most frequently affected proteins included the androgen receptor (AR) and the estrogen receptor-alpha (ESR1), which were each linked to 23 EDCs, whereas 61 individual proteins were associated with only one EDC.

To identify the biological targets that are most often affected by EDCs, proteins were grouped in clusters according to their families, as based on the Panther classification system (Fig. 2 ). The majority of the 208 proteins related to EDCs belonged to 12 classes among the 23 present in Panther, while the remaining proteins were classified as ‘uncategorized protein class’. Each protein was assigned to only one category, although only one of them, HLA-DRA, (HLA class II histocompatibility antigen, DR alpha chain) belonged to the defense/immunity group. Other immunity-related proteins. such as interleukin 6 (IL-6) or interleukin 1 alpha (IL-1A), were not associated with any class in the Panther classification. We therefore manually added all immune system-related proteins to the “uncategorized class”. Given that Bisphenol A (BPA) increases the release of these proteins(Ben-Jonathan et al., 2009), and because antibodies to the IL-6 receptor (such as tocilizumab) or to the IL-1 receptor (such as Anakinra) are currently tested for the treatment of COVID-19 patients (Zhou et al., 2020), we also explored if the proteins selected could be mapped to defense and/or immunity biological categories. For this purpose, we used the Gene Ontology (GO) classification (as of May 26, 2020), and among the 208 proteins dysregulated by EDCs, 58 were associated with inflammatory response, 75 with defense response, and 66 with regulation of immune system process.

Fig. 2

Tripartite network representation of endocrine-disrupting chemicals-proteins-diseases relationships. First, a bipartite network of the 208 human proteins known to be dysregulated by the 30 endocrine-disrupting chemicals (EDCs) was created as extracted from the CompTox database. Each yellow diamond node represents an EDC, and edges are the interactions between EDCs and proteins. Then, a second bipartite network was generated for the 4184 human proteins known to be linked to the 13 predisposing diseases, as extracted from the DisGeNET (3262 links) and GeneCards (2079 links) databases. Each red square node represents a disease, and edges are the interactions between diseases and proteins. A total of 1156 proteins were overlapping. All proteins were grouped using the Panther classification system (version 15) and are represented by circles (colors are according to their Panther family classes). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Regarding diseases associated with human proteins, two databases were screened. From the DisGeNET database, we were able to retrieve information for 8 of the 13 diseases, which were connected to 3262 unique proteins via 7195 links (as of April 29, 2020). The proteins were categorized in 22 protein classes using the Panther classification (version 15) (Figure S1). Proteins that did not belong in any class were again grouped into the uncategorized class. Obesity and diabetes were linked to proteins belonging to each of the 22 categories, whereas insulin resistance and dyslipidemia were linked to only half of the categories.

From the GeneCards database, all 13 predisposing diseases were retrieved (as of 29 April 2020), and a total of 115,289 associations were identified between the diseases and 29,094 unique human proteins were extracted. Among them, only protein-coding information according to HGNC, Ensembl or Entrez Gene were kept (proteins data related to biological regions, gene clusters, genetic loci, pseudogenes, non-coding RNA genes and uncategorized elements were not considered), thereby reducing the total number of unique protein to 18,931, representing 97,855 disease-protein links. As a next step, grouping of the proteins using the Panther classification system allowed identification of 23 clusters correspond to the 23 different protein classes (Figure S2). Each protein was assigned to only one category, except for ameloblastin (AMBN), which was associated with both ‘extracellular matrix protein’ and ‘structural protein’. Proteins not associated with Panther classes, were again grouped into the uncategorized class. Excluding the viral or transposable element protein class, all diseases (except dyslipidemia) were associated with all the other Panther classes.

In order to keep the most relevant protein-disease associations obtained from the GeneCards database, data were filtered based on their scores. The GeneCards scores are calculated based on publications mentioning a protein and a disease, using a Boolean model. The higher the score, the more relevant the protein-disease association is. Among the 97,855 links between the 13 diseases and 18,931 proteins, the score values ranged between 0.13 (representing very low association) to 228 (very high evidence for a protein-disease connection). After evaluation of the extracted data (number of proteins by GeneCards scores), we selected associations with a score ≥ 20 (see Figure S3). Within this threshold, a total of 5732 associations were retained that link the 12 diseases with 2079 unique human proteins (no information was retained for ‘dyslipidemias’ from the GeneCards database).

Generating a tripartite network of protein-EDC-disease associations

A human bipartite associative network of proteins and the 13 diseases was created. Among the 3262 unique proteins from the DisGeNET, and the 2079 proteins from the GeneCards databases, 1157 were overlapping proteins and only 922 and 2105 proteins were uniquely associated with GeneCards or DisGeNET, respectively. All 4184 unique proteins were again grouped into 23 clusters using the Panther classification (the class ‘viral or transposable element protein’ was not kept after the cleaning step. Among the groupings, we retrieved a cluster of proteins linked to the ‘defense/immunity’ category. These results were merged with the bipartite protein-EDCs network to develop a tripartite network (Fig. 2).

Translation into pathways

To identify biological pathways that may be involved in the predisposing diseases while also being dysregulated by the EDCs, we first analyzed the overlaps between the two sets of proteins. Among the proteins identified from the three data sources, 98 were common (Figure S4), and all of them were mapped to unique Entrez GeneID, and could therefore be used for biological enrichment analyses, which were performed independently using four data sources (KEGG, Reactome, Wiki-pathways and Panther). The ORA analysis revealed several statistically significant pathways linked to interleukins/cytokines signaling, intracellular signaling pathways and, regulation of metabolic pathways (Table 3 ). Interestingly, the different data sources showed very significant associations with common pathways, such as interleukins (IL) related pathways: IL-4 and IL-13, IL-10 signaling pathways (padj < E-16, and padj of 2.85E-09 respectively, Reactome), IL-17 signaling pathway (padj of 1.05E-10, KEGG), IL-3, IL-5 and IL-18 signaling pathways (padj of 1.09E-09, 2.52E-09, 1.49E-08 respectively, Wiki-pathways), the IL-signaling pathways (padj of 1.10E-05, Panther); or the Toll-like receptor signaling pathway (padj 3.91E-09 for KEGG, padj of 0.99 for Panther and padj 1.93E-08 for Wiki-pathways).

Table 3

Pathway enrichment for the set of proteins that are linked to both the predisposing diseases and to the EDCs. The pathways were extracted from the KEGG, Panther, Reactome and the Wikipathways database.

Data sources	Name of pathways	Proteins*	P-value	FDR**
KEGG	AGE-RAGE signaling pathway in diabetic complications	22	< E−16	< E−16
KEGG	Fluid shear stress and atherosclerosis	20	2.22E−16	1.81E−14
KEGG	TNF signaling pathway	18	8.8E−16	5.79E−14
KEGG	Insulin resistance	17	1.04E−14	5.67E−13
KEGG	Endocrine resistance	15	7.81E−13	2.40E−11
KEGG	MAPK signaling pathway	23	8.84E−13	2.40E−11
KEGG	HIF−1 signaling pathway	15	1.06E−12	2.66E−11
KEGG	Non-alcoholic fatty liver disease (NAFLD)	17	2.83E−12	5.77E−11
KEGG	FoxO signaling pathway	16	5.17E−12	9.36E−11
KEGG	IL-17 signaling pathway	14	6.11E−12	1.05E−10
KEGG	EGFR tyrosine kinase inhibitor resistance	13	1.14E−11	1.85E−10
KEGG	PI3K-Akt signaling pathway	23	3.88E−11	5.51E−10
KEGG	Prolactin signaling pathway	12	4.50E−11	6.12E−10
KEGG	Ras signaling pathway	19	4.91E−11	6.40E−10
KEGG	Thyroid hormone signaling pathway	14	1.32E−10	1.39E−09
KEGG	Toll-like receptor signaling pathway	13	4.08E−10	3.91E−09
KEGG	Insulin signaling pathway	14	1.25E−09	1.07E−08
KEGG	Human T-cell leukemia virus 1 infection	18	1.93E−09	1.57E−08
KEGG	Chronic myeloid leukemia	11	2.02E−09	1.58E−08
KEGG	B cell receptor signaling pathway	10	1.50E−08	1.09E−07
KEGG	T cell receptor signaling pathway	11	4.34E−08	3.01E−07
KEGG	C-type lectin receptor signaling pathway	11	5.90E−08	4.01E−07
Panther	Interleukin signaling pathway	12	1.94E−07	1.10E−05
Panther	Insulin/IGF pathway-protein kinase B signaling cascade	7	1.51E−05	3.42E−04
Panther	Ras Pathway	9	3.61E−05	6.80E−04
Panther	T cell activation	9	6.34E−05	8.95E−04
Panther	PI3 kinase pathway	7	1.12E−04	0.0012
Panther	Insulin/IGF pathway-mitogen activated protein kinase kinase/MAP kinase cascade	5	5.81E−04	0.0050
Panther	Inflammation mediated by chemokine and cytokine signaling pathway	13	8.29E−04	0.0062
Panther	B cell activation	6	0.0026	0.0154
Panther	FGF signaling pathway	8	0.0034	0.0183
Panther	EGF receptor signaling pathway	8	0.0063	0.0324
Panther	Interferon-gamma signaling pathway	2	0.1423	0.5544
Panther	JAK/STAT signaling pathway	1	0.3137	0.9401
Panther	Toll receptor signaling pathway	2	0.3521	0.9947
Reactome	Signaling by Interleukins	31	< E−16	< E−16
Reactome	Interleukin-4 and Interleukin-13 signaling	21	< E−16	< E−16
Reactome	Cytokine Signaling in Immune system	33	5.55E−16	3.20E−13
Reactome	Interleukin-10 signaling	10	1.05E−11	2.85E−09
Reactome	Negative regulation of the PI3K/AKT network	13	1.16E−11	2.85E−09
Reactome	PIP3 activates AKT signaling	18	1.83E−11	3.51E−09
Reactome	PI3K/AKT Signaling in Cancer	12	6.79E−11	1.17E−08
Reactome	Cytochrome P450 - arranged by substrate type	10	3.70E−10	5.81E−08
Reactome	PI5P, PP2A and IER3 Regulate PI3K/AKT Signaling	11	1.43E−09	2.07E−07
Reactome	Signaling by Receptor Tyrosine Kinases	19	2.43E−08	2.70E−06
Reactome	Insulin receptor signalling cascade	8	2.79E−08	2.84E−06
Reactome	Signaling by Insulin receptor	8	5.25E−07	3.27E−05
Reactome	MAPK family signaling cascades	13	2.57E−06	1.14E−04
Reactome	Constitutive Signaling by Aberrant PI3K in Cancer	7	3.66E−06	1.58E−04
Reactome	Immune System	37	7.14E−06	2.94E−04
Wiki-pathway	Netrin-UNC5B signaling Pathway	15	2.22E−16	1.18E−13
Wiki-pathway	Nonalcoholic fatty liver disease	20	3.06E−14	2.71E−12
Wiki-pathway	Aryl Hydrocarbon Receptor Netpath	12	1.32E−12	6.36E−11
Wiki-pathway	AGE/RAGE pathway	13	3.93E−12	1.61E−10
Wiki-pathway	Insulin Signaling	18	5.35E−12	2.03E−10
Wiki-pathway	RAC1/PAK1/p38/MMP2 Pathway	13	7.18E−12	2.38E−10
Wiki-pathway	Relationship between inflammation, COX-2 and EGFR	9	2.42E−11	7.32E−10
Wiki-pathway	IL-3 Signaling Pathway	11	4.31E−11	1.09E−09
Wiki-pathway	Ras Signaling	18	5.70E−11	1.38E−09
Wiki-pathway	IL-5 Signaling Pathway	10	1.09E−10	2.52E−09
Wiki-pathway	PI3K-Akt Signaling Pathway	23	2.10E−10	4.13E−09
Wiki-pathway	Aryl Hydrocarbon Receptor Pathway	10	7.67E−10	1.13E−08
Wiki-pathway	IL-18 signaling pathway	20	1.09E−09	1.49E−08
Wiki-pathway	Cells and Molecules involved in local acute inflammatory response	7	1.48E−09	1.91E−08
Wiki-pathway	Toll-like Receptor Signaling Pathway	13	1.52E−09	1.93E−08

Number of proteins from the studied set that is involved in a pathway.

False discovery rate.

Among the most significant pathways, several were retrieved from each of the data sources with relation to the AGE/RAGE pathway (i.e. Advanced Glycation End products and its receptor), which is known to cause cellular stress and inflammation. The AGE are formed non-enzymatically, by Maillard reaction products (carbohydrates with proteins and/or lipids) and bind to the RAGE. Formation of AGE has been associated with chronic diseases such as type 2 diabetes (Cai et al., 2012, Menini et al., 2018). Similarly, the stress or inflammatory pathways (e.g. shear stress, defined as the tangential force exerted by the blood flow on the vascular endothelium, TNF-alpha) are highlighted by our analysis; the shear stress activates the AhR signaling pathway, which is also involved in the regulation of IL-17 production by the Th17 lymphocytes; interleukin 17 has been suspected to be involved in the pathogenesis of COVID19 (Guti19nvolved in the regulation of IL-17 production by Pacha et al., 2020). Interestingly, inflammation is suspected to influence insulin resistance.

We went one step further and explored at the chemical level the linkage between the main identified pathways via the proteins that chemicals are known to target. Many chemicals target proteins that are common to several pathways, and thus most of the EDCs examined overlap with several pathways, as illustrated by four main KEGG pathways identified (Fig. 3 ). Regarding the pathways AGE/RAGE (AR in Fig. 3) and IL-17 (IL in Fig. 3), most of the linked chemicals were common, probably due to the high number of overlapping proteins in these pathways (Fig. 3). Indeed, among the 22 proteins targeted by EDCs and involved in the AGE/RAGE pathways, and the 14 proteins for the IL-17 pathways, eight were common to both pathways (JUN, CCL2, CXCL8, MAPK3, NFKB1, IL6, MAPK1, TNF). However, our studies revealed that bisphenol compounds and PFAS substances may be linked to both pathways.

Fig. 3

Venn diagram Illustrating the distribution of chemical links to different KEGG pathways identified after biological enrichment. The following four pathways were considered as an example: IR (Insulin resistance), ER (endocrine resistance), IL (IL-17) and AR (AGE/RAGE). Number represent the number of EDC chemicals linked to each pathway. The Venn diagram shows that the majority of chemicals are linked to the four pathways and few are specific. This may be related to the overlapping proteins in these pathways (see text).

Exploration of EDCs linkage to COVID-19

To explore putative links between COVID-19 and exposure to EDCs, we first screened the AOP-Wiki database, and then further examined the pathways identified using literature references.

In the AOP-Wiki database, only one AOP was related to COVID-19, and it involves several key events, such as ‘increased pro-inflammatory mediators’ (KE 1496), ‘increased inflammatory immune responses’ (KE 1750), which leads to the adverse outcome ‘increased mortality’ (AO 351). Such knowledge-based linear chain of events highlights the importance of the link between COVID-19 and inflammatory processes.

Discussion

In order to investigate possible links between exposure to EDCs and the severity of COVID-19, we explored a computational systems biology approach. The tripartite network model first linked EDCs to targeted proteins and then proteins related to diseases that predispose to more serious COVID-19 development, thereby allowing us to identify common signaling pathways. The identification of such joint pathways and their role as possible targets of EDCs highlights the potential links between exposure to environmental chemicals and COVID-19 severity.

This integrative approach can be easily applied as a new approach methodology (NAM) (Bopp et al., 2019), which may offer support to methods alternative to animal testing or to identify biological pathways that require more focused laboratory study. Previous studies have demonstrated that systems chemical toxicology models combined with computational network biology may help in understanding chemical toxicity in humans (Hartung et al., 2017, Nie et al., 2015, Taboureau and Audouze, 2017). Our tripartite network supports the notion that exposure to EDCs may contribute to aggravation of COVID-19. Although major links were identified at extremely low p values, the approach relies on existing information available in within the very substantive data sources, but some causal associations may have been overlooked or disregarded because of missing or incomplete information. For example TCDD was connected to the endocrine resistance pathways via the ESR1 protein, and bisphenol compounds and PFASs were linked to the IL-17 signaling pathways (via four proteins for BPs, ten proteins for PFAS) and to the AGE/RAGE signaling pathways (via five proteins for BPs, 17 proteins for PFAS) (Fig. 3).

To assess the validity of our approach, a more focused expert analysis was attempted, where we selected the IL-17 and the AGE/RAGE signaling pathways because of their pathophysiological relevance in the context of COVID-19. The interleukin-17 (IL-17) signaling pathway plays several important roles, and IL-17 is produced by a pro-inflammatory subtype of T helper lymphocytes named Th17 cells, located at mucosal barriers where they contribute to pathogen clearance. The IL-17 produced stimulates the synthesis of cytokines (IL1ß, TNF-alpha…) and chemokines (MCP-1…) by other cell types, thereby favoring the recruitment of monocytes and neutrophils at inflammatory sites. However, an over-activation of Th17 cells can lead to a hyper-inflammatory state which is deleterious (Pacha et al., 2020).

The highly variable symptomatology associated with the infection by SARS-CoV-2 depends on the levels of IL-17 and of other cytokines including IL-1ß, IL-6, IL-15, TNF-alpha and IFNγ. The most deleterious effect of SARS-CoV-2 in humans is an acute lung injury leading to a severe acute respiratory syndrome (SARS) that is partly due to IL-17-related excessive recruitment of pro-inflammatory cells and production of pro-inflammatory cytokines. Therefore, an increased basal level of IL-17 (in the absence of infection, for example due to obesity or to induction by a chemical) might represent a lung injury risk associated with SARS-CoV-2 infection. Our finding of EDC linkage to this pathway is therefore of high pathogenetic relevance.

Obesity promotes a high basal level of inflammation which contributes to insulin resistance and type 2 diabetes (Goldberg, 2009). This phenomenon is due to an infiltration of the adipose tissue (AT) by macrophages and T cells and their production of various pro-inflammatory cytokines, including IL-1ß, TNF-alpha, IL-17 and IL-6. Several EDCs are suspected to be obesogenic (and are subsequently named obesogens). This has been demonstrated for several substances (e.g. tributyltin) and linked to the stimulation of pro-adipogenic signaling pathway (e.g. through PPARγ) (Egusquiza and Blumberg, 2020). Similarly, the aryl hydrocarbon receptor (AhR) is highly expressed in Th17 cells and is an essential contributor to the production of IL-17(Veldhoen et al., 2008). The AhR, known as the receptor of dioxins and dioxin-like PCBs, is also activated by shear stress (SS), another pathway highlighted in our computational analysis. Indeed, several studies have shown using various endothelial models that laminar SS leads to the activation of two target genes of the AhR, namely CYP1A1 and CYP1B1 (Conway et al., 2009). Two recent studies suggest an indirect link between SARS-CoV-2 and SS by showing that the expression of ACE2 (angiotensin-converting enzyme 2), the receptor of the virus, is increased by SS (Song et al., 2020).

These observations support a dual impact of EDCs on IL-17 production and inflammatory state; this impact could be indirect due to the effect of these chemicals on obesity or through a direct stimulation of several signaling pathways, such as AhR or PPARγ, leading to an overproduction of systemic IL-17; the shear stress pathway represents an additional link between AhR activation and the EDC/disease connection. The implication of shear stress also suggests a possible contribution of increased expression of ACE2, the receptor of the SARS-CoV-2. While the role of these pathways at the nexus between exposure to EDCs and COVID-19 severity appears to be relevant, their actual contribution remains to be demonstrated and their putative role as therapeutic targets remains to be further substantiated.

Our integrative systems biology study also indicates a strong statistical association between the AGE/RAGE signaling pathway, chronic diseases and EDC effects. This is likely due to the well-known links between this pathway and type 2 diabetes (Ravichandran et al., 2019). Indeed, hyperglycemia leads to increased amounts of glycation products and their metabolites which results in the activation of the RAGE receptors. The latter are highly expressed in endothelial cells, and their activation leads to increased oxidative stress and inflammation and ultimately to endothelial damage, thrombotic disorders and vascular diseases (Egaña-Gorroño et al., 2020). Other endogenous ligands can also activate RAGE, among them HMGB1 (high-mobility group box 1), an extra-cellular protein also linked to a variety of inflammatory responses (Andersson et al., 2020). Interestingly, the AGE/RAGE signaling pathway is highly expressed in the lung vasculature and has been implicated in several pulmonary diseases (Oczypok et al., 2017). All these observations support the implication of the AGE/RAGE signaling pathway in vascular, thrombotic and lung diseases which are the hallmarks of COVID-19 severity. Interestingly, there are also complex connections between HMGB1 and ACE2 which is the receptor for SARS-Cov2 and other coronaviruses (Luft, 2016). These results are in accordance with recent proposals in published commentaries of environmental chemical impacts on COVID-19 progress (Andersson et al., 2020, Rojas et al., 2020).

The three-way approach did not attempt to identify direct immunotoxic effects due to environmental chemicals otherwise considered to be EDCs. However, some of the EDCs selected, i.e., PCB-153, PFOA and PFOS, are known to have immunotoxic properties (Heilmann et al.), and the same is true for some common air pollutants (Tsatsakis et al., 2020). Accordingly, the impact of environmental chemicals on COVID-19 severity demands attention.

Inflammation appears to be a critical mechanism for both EDCs and COVID-19, with many possible implications that could be foreseen. First, non-EDCs pollutants also trigger inflammation, and the question of a mixture effects between these pollutants and EDCs is highly relevant. From a risk assessment perspective, it would be useful to test whether inflammatory markers constitute relevant effect markers of such contaminants and thus provide a link to the pathogenesis of relevant diseases. The present study therefore opens new perspectives for research on the interaction between chemicals, chronic and infectious diseases. As an illustration, the immune system may serve as a link between chemical and biological stressors. While evidence already demonstrates the role of the AhR system in regard to the regulation of TH17 cells and other immune targets, our study suggests a need to examine the role of the other signaling pathways triggered by EDCs and on the interaction of these pathways with, e.g., the AhR pathway.

The increased availability of high throughput data, will continue to improve the precision and robustness of computational modeling. For example, in a previous pilot study, the authors used three levels of evidence concerning chemical-disease links to perform a human environmental disease network model (Taboureau and Audouze, 2017). Other possible extensions of such models would be to integrate more data, quantitative information such as dose levels, and the biological complexity and organization at several levels (cells, tissues and organs) (Taboureau et al. 2020) in order to evaluate more closely the potential toxicity of EDCs on human health. In addition, novel in vitro approaches using human organoids or 3D cell systems are also critical to improve our knowledge on the effects of chemicals in agreement with the 3R principle. These studies will be complementary to the ones presented here, thus highlighting the potential of applying and further developing in silico models to identify potential harmful effects from chemical exposure.

Conclusion

The results of this computational study appear as a promising initial step toward systematically linking a major group of environmental chemicals to the severity of COVID-19, although the findings need to be further supported by high-throughput screening tests, clinical and experimental data. Nevertheless, these observations bridge environmental stressors and infectious diseases and support an integrated exposome approach. Preliminary focus on the AGE/RAGE and IL-17 pathways illustrates the potential connection between exposure to EDCs and diseases predisposing to COVID-19 severity.

Funding

This work was supported by OBERON (https://oberon-4eu.com/), a project funded by the European Union’s Horizon 2020 research and innovation program under grant agreement no. 825712; and by the National Institute of Environmental Health Sciences, NIH (P42ES027706).

CRediT authorship contribution statement

Qier Wu: Resources, Methodology, Software, Visualization. Xavier Coumoul: Investigation, Writing - review & editing. Philippe Grandjean: Writing - review & editing, Funding acquisition. Robert Barouki: Investigation, Writing - review & editing. Karine Audouze: Conceptualization, Methodology, Visualization, Writing - original draft, Supervision, Funding acquisition.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.