Front. Endocrinol.Frontiers in EndocrinologyFront. Endocrinol.1664-2392Frontiers Media S.A.10.3389/fendo.2023.1200211EndocrinologyOriginal ResearchInvestigation of the effects of phthalates on in vitro thyroid models with RNA-Seq and ATAC-SeqNazzariMarta1RomittiMírian2HauserDuncan1CarvalhoDaniel J.3GiselbrechtStefan3MoroniLorenzo4CostagliolaSabine2CaimentFlorian1*1Department of Toxicogenomics, GROW School for Oncology and Reproduction, Maastricht University, Maastricht, Netherlands2Institute of Interdisciplinary Research in Molecular Human Biology (IRIBHM), Université Libre de Bruxelles, Brussels, Belgium3Department of Instructive Biomaterials Engineering, MERLN Institute for Technology-Inspired Regenerative Medicine, Maastricht University, Maastricht, Netherlands4Department of Complex Tissue Regeneration, MERLN Institute for Technology-Inspired Regenerative Medicine, Maastricht University, Maastricht, Netherlands
Edited by: Terry Francis Davies, Icahn School of Medicine at Mount Sinai, United States
Reviewed by: Taisen Iguchi, Graduate University for Advanced Studies (Sokendai), Japan; Rauf Latif, Icahn School of Medicine at Mount Sinai, United States
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
Introduction
Phthalates are a class of endocrine-disrupting chemicals that have been shown to negatively correlate with thyroid hormone serum levels in humans and to cause a state of hyperactivity in the thyroid. However, their mechanism of action is not well described at the molecular level.
Methods
We analyzed the response of mouse thyroid organoids to the exposure to a biologically relevant dose range of the phthalates bis(2-ethylhexyl) phthalate (DEHP), di-iso-decylphthalate (DIDP), di-iso-nonylphthalate (DINP), and di-n-octylphthalate (DnOP) for 24 h and simultaneously analyzed mRNA and miRNA expression via RNA sequencing. Using the expression data, we performed differential expression and gene set enrichment analysis. We also exposed the human thyroid follicular epithelial cell line Nthy-ori 3-1 to 1 µM of DEHP or DINP for 5 days and analyzed changes in chromatin accessibility via ATAC-Seq.
Results
Dose-series analysis showed how the expression of several genes increased or decreased at the highest dose tested. As expected with the low dosing scheme, the compounds induced a modest response on the transcriptome, as we identified changes in only mmu-miR-143-3p after DINP treatment and very few differentially expressed genes. No effect was observed on thyroid markers. Ing5, a component of histones H3 and H4 acetylation complexes, was consistently upregulated in three out of four conditions compared to control, and we observed a partial overlap among the genes differentially expressed by the treatments. Gene set enrichment analysis showed enrichment in the treatment samples of the fatty acid metabolism pathway and in the control of pathways related to the receptor signalling and extracellular matrix organization. ATAC-Seq analysis showed a general increase in accessibility compared to the control, however we did not identify significant changes in accessibility in the identified regions.
Discussion
With this work, we showed that despite having only a few differentially expressed genes, diverse analysis methods could be applied to retrieve relevant information on phthalates, showing the potential of in vitro thyroid-relevant systems for the analysis of endocrine disruptors.
Phthalates are a class of manmade compounds used in the manufacturing industry as solvents or added as plasticizers, mainly polyvinyl chloride (PVC) or other polymers, to confer flexibility and softness (1, 2). Phthalates are alkyl or dialkyl esters of phthalic acid, and their functional groups can be linear, branched, or circular (3). Depending on their size, phthalates are classified into low and high molecular weights (MW) (4, 5). Low-MW phthalates include benzyl butyl phthalate (BBP), diethyl phthalate (DEP), di-iso-butyl phthalate (DiBP), dimethyl phthalate (DMP), and di-n-butyl phthalate (DnBP), while high-MW ones comprise bis(2-ethylhexyl) phthalate (DEHP), di(2-propylheptyl) phthalate (DPHP), di-iso-decylphthalate (DIDP), di-iso-nonylphthalate (DINP), and di-n-octylphthalate (DnOP) (6). They are found in common household items, medical devices, construction materials, and consumer products (6). Since they are not covalently bound to the plastic matrix they are contained in, phthalates can leach or gas out and contaminate either the environment or be ingested via contaminated food (5). Indeed, food constitutes one of the biggest sources of human exposure to phthalates (7). Human biomonitoring studies conducted on the general population in Asia, Europe, and North America show a widespread exposure of the general population to phthalates (6, 8–11). Despite their broad use and pervasive environmental presence, they have been recognized as toxic substances both in humans and other organisms (12). Once ingested, they are rapidly metabolized in the digestive tract to their monoester form, which is the species responsible for the phthalates’ toxicity. Low-MW phthalate metabolites are then excreted through the urine, while high MW metabolites are excreted both via the urine and feces (1). While they do not bioaccumulate, their persistent exposure of the population is a matter for concern. Short- and medium-chain phthalates have been associated with higher toxicity than long-chain ones, which has led to their banning or restriction in children’s toys or teething products (13, 14). Some of the examined compounds have been reported to cause chronic or subchronic toxicity in several organs and systems, namely the liver, kidney, immune system, testes, uterus, ovary, central nervous system, and thyroid in vivo (2, 15–17). Phthalates can also negatively interfere with the endocrine system and are thus considered endocrine-disrupting chemicals (EDCs). They have been shown to interfere with prenatal and postnatal development in animal models (18), with the female and male reproductive systems (19–21), as well as being possibly linked to obesity and type 2 diabetes (22–24).
The thyroid is an endocrine gland positioned in the lower part of the anterior neck and is responsible for the production of the thyroid hormone (TH), whose receptors are expressed throughout the body (25). The TH is essential for normal growth and development and metabolism regulation (26). Its production is mainly regulated by the thyroid-stimulating hormone (TSH), which is secreted by the adenohypophysis. In turn, TSH production is regulated both by circulating TH levels and the thyrotropin-releasing hormone (TRH), synthesized in the hypothalamus. The main cell type of the thyroid is constituted by thyrocytes, which organize in small hollow spheres called follicles and are responsible for synthesizing the TH, which is stored in the center of the follicle (the lumen) in a dense matrix termed colloid. The synthesis of the TH starts with the active transport of iodine inside the thyrocyte via the sodium iodide symporter (NIS in humans). In the follicle lumen, it is covalently bound via oxidation to the tyrosyl (Tyr) residues of the protein thyroglobulin (TG) via the action of the membrane-bound enzyme thyroid peroxidase (TPO) (27, 28). Following TSH stimulation, TG is degraded in the lysosomes, freeing TH, which can be transported outside the thyrocyte.
In the thyroid, phthalate treatment has been shown to have an effect in vitro and in vivo, causing histological changes such as reduced follicle size and colloid density, hypertrophy of the Golgi apparatus, an increase in the number and size of lysosomes, and alteration of the TH levels (2, 15, 16, 29, 30). DEHP has been shown to downregulate Tshr (Tsh receptor) expression and interfere with the Tsh/Tshr signaling pathway in vivo (31, 32). In humans, the presence of phthalate metabolites in urine has been observed in association with alterations in TH and TSH serum levels (22, 29). In addition, there is evidence for phthalates altering the methylation status in sperm cells (33) and adrenal glands (34) of the offspring of exposed rats, as well as the expression or activity of histone deacetylases and histone methyltransferases (35, 36).
Over the years, great effort has been made to develop thyroid organoids using both embryonic and induced pluripotent stem cells (37) that can be used for developing thyroid disease models (38) and performing cancer (39) and toxicological and drug screening (40). In the context of toxicology, in vitro models can offer high throughput capability, and mechanistic insight into endocrine disruption and reduce the use of animal testing, in line with the 3Rs principles for animal welfare (replacement, reduction, and refinement) (41).
In this work, we analyzed our ability to identify alterations induced by phthalate treatment by using two in vitro models of the thyroid. To this end, we exposed mouse embryonic stem cell (mESC)-derived thyroid follicles (42) to the high-MW phthalates DEHP, DIDP, DINP, and DnOP for 24 h and analyzed the transcriptome via RNA-Sequencing (RNA-Seq) using the Combo-Seq library prep kit for simultaneous analysis of mRNA and miRNA expression. Data analysis revealed the upregulation of the growth protein 5 (Ing5) gene in three out of four tested compounds (DEHP, DINP, DnOP) compared to the control. ING5 is a component of the histones H3 and H4 acetyltransferase complexes HBO1-JADE, HBO1-BRPF1, and MOZ/MORF (43, 44). To investigate the potential effect of phthalate treatment on the chromatin status of thyroid models in vitro, we exposed the human thyroid follicular epithelial cell line Nthy-ori 3-1 to 1 µM of DEHP or DINP for 5 days and analyzed the genome accessibility with Assay for Transposase-Accessible Chromatin (ATAC)-Seq. We used maSigPro to analyze gene expression across the dose series and performed gene set enrichment analysis (GSEA) to identify enriched pathways.
Materials and methodsChemicals information
The following phthalates were used for the experiments described in this paper: DEHP (CAS 117-81-7; purity: 99.8% ± 0.4%) (67261, Sigma-Aldrich, St. Louis, MO, USA), DINP (CAS 28553-12-0; ester content: ≥ 99% mixture of C9 isomers) (376663, Sigma-Aldrich), DIDP (CAS 26761-40-0; purity: ≥ 99.0%) (80135, Supelco, St. Louis, MO, USA), and DnOP (CAS 117-84-0; purity ≥ 99.5%) (D201154, Sigma-Aldrich).
Organoids differentiation
Thyroid organoids were differentiated from the A2Lox.Cre_TRE-Nkx2-1/Pax8_Tg-EGFP mouse ESC and enriched as previously described (42, 45). For more information, see Supplementary Methods.
Exposure to phthalates and RNA-Seq library preparationExposure to phthalates
In low-adhesion 48-well cell culture plates, 1,000 follicles per well were seeded in triplicate using supplemented differentiation medium (Supplementary Methods) and 1–10–100 nM to 1–10 μM of DEHP, DINP, or DIDP or 2–20–200 nM to 2–20 μM of DnOP dissolved in DMSO (1029521000, Merck Millipore, Burlington, MA, USA) (final DMSO concentration: 0.5%). Of note, the slightly different dose range for DnOP was caused by an unwanted dilution error. We decided to still consider DnOP not differently than the other three phthalates in our following data analysis, considering the dose range still maintains the same scaling between each dose and its order of magnitude is comparable to the others. As a control, 1,000 follicles per well were seeded in a supplemented differentiation medium and 0.5% DMSO (n = 5) or medium alone (n = 3). The plated follicles were incubated at 37°C, 5% CO2, and > 95% humidity for 24 h.
RNA isolation
After 24 h, the follicles were collected, washed once with PBS, and lysed in QIAzol Lysis Reagent (79306, Qiagen, Venlo, The Netherlands). Total RNA was extracted using the miRNAeasy Micro Kit (217084, Qiagen). All samples had a RNA integrity number (RIN) of 8 or higher.
RNA-Seq library preparation
The NEXTFLEX® Combo-Seq™ mRNA/miRNA Kit (NOVA-5139-53, PerkinElmer, Waltham, MA, USA) was used to prepare RNA-Seq libraries using 20 ng of total RNA. To deplete tRNA and Y RNA fragments, the NEXTFLEX® tRNA/YRNA Blocker was used during the library preparation following the manufacturer’s instructions. In total, 14 cycles of PCR were performed during the protocol. For some samples, the final library concentration was below the pooling concentration used for sequencing (1.6 nM). In these cases, the library was prepared again to perform 16 cycles. For three samples (10 μM of DEHP: replicate 3; 1 nM of DIDP: replicate 3; untreated control: replicate 1), there was not enough RNA to repeat the library preparation and could thus not be sequenced. The prepared libraries were sequenced on an S4 Illumina flowcell 35 cycles (v1.5) (Illumina) in single-end mode.
Exposures to DEHP or DINP and ATAC-Seq library preparationExposure to DEHP or DINP
The human thyroid follicular epithelial cell line Nthy-ori 3-1 was plated at a density of 10,000 cells/cm2 on six-well plates and cultured in RPMI 1640 Medium with GlutaMAX™ Supplement (61870036, Gibco, Waltham, MA, USA), 10% FBS, and penicillin–streptomycin (15140122, Gibco) and incubated at 37°C, 5% CO2, and > 95% humidity. Cells were left 1 day to adhere, and the following day, the medium was changed to culture medium with DEHP (n = 6) or DINP (n = 6) at 1 μM in 0.5% DMSO. As solvent control, the culture medium was added with just 0.5% DMSO (n = 6). Cells were incubated for 5 days, refreshing the media with the compound or DMSO only at day 3. At the end of the incubation period, cells were collected and counted manually.
ATAC-Seq libraries preparation
To prepare ATAC-Seq libraries, 50,000 cells per sample were used. Libraries were prepared following the Omni-Atac protocol of Corces et al. (46) with the replacement of NP40 from the original protocol with IGEPAL (I8896: 50 ml, Merck Millipore). The tegmental kit used was the Illumina Tagment DNA Enzyme and Buffer Small Kit (20034197, Illumina, San Diego, CA, USA) and the indexes IDT® for Illumina® DNA/RNA UD Indexes Set A, Tagmentation (96 indexes, 96 samples) (20027213, Illumina). Seven PCR cycles were used for all samples. The prepared libraries were sequenced on an SP Illumina Flowcell v1.5 (100 cycles) (Illumina) in paired-end mode.
All RNA samples and sequencing library concentrations were measured with the Qubit 2.0 Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA), and quality control was performed on a BioAnalyzer 2100 expert (Agilent, Santa Clara, CA, USA) or a 2200 TapeStation System (Agilent).
Data analysis
All the scripts used for RNA-Seq and ATAC-Seq data analysis have been collected in a markdown file available at https://github.com/marta-nazzari/phthalates_rnaseq_atacseq.
RNA-Seq data processing
The.fastq files were processed according to our previously published CODA pipeline (47). Briefly, reads were trimmed from the 5′ 4N and 3′ 8A adapters using Cutadapt (v3.7) (48), as recommended by the manufacturer (49). To obtain gene read counts, trimmed reads were aligned to the mouse transcriptome (GRCm39 v27) and quantified using RSEM (v1.3.3) with the “–STAR” parameter (v2.7.10a), following the ENCODE3’s STAR-RSEM pipeline (50, 51). To analyze miRNAs, the trimmed files were used as input for miRge3.0 (v0.0.9) (52) using miRBase mouse annotations (v22).
ATAC-Seq data processing
The.fastq files were preprocessed using the PEPATAC pipeline (v0.10.4) (53) using bowtie2 (v2.4.2) (54) as mapper, samtools (v1.4) (55) as deduplicator, and the included Python tool “pyadapt” as trimmer. The human genome GRCh38 v38 build was used for alignment.
RNA-Seq sample biotype mapping and outlier identification
Quantified RNA species were mapped to their respective biotypes using the R (56) package biomaRt (57). We calculated the percentage of mapped reads per biotype and retained only those constituting at least 1% in at least one sample. Outliers for each biotype were identified per treatment group (DEHP, DIDP, DINP, DnOP, DMSO, and untreated) and calculated as being 1.5 times the interquartile range (IQR) below the 25th percentile or above the 75th percentile:
Biotype x in sample y<25th percentile (biotype xin group z)−1.5*IQR (biotypex in group z)
or
Biotype xin sample y>75th percentile (biotype x in group z)+1.5*IQR(biotype x in group z)MaSigPro analysis
Normalized gene counts were used for maSigPro (58) analysis according to the maSigPro user’s guide for next-generation sequencing data (59) for a single-series course experiment. We set the “tetha” (θ) value to 10 (default), the FDR to 0.05 (default), and the “degree” parameter to 3 (this corresponds to a cubic polynomial regression model). The variable “time” with values 0, 1, 2, 3, 4, and 5 was used in the model to represent the “dose” values of 0 (DMSO control), 1 nM, 10 nM, 100 nM, 1 μM, and 10 μM (or 2 nM, 20 nM, 200 nM, 2 μM, and 20 μM for DnOP).
Differential expression analysis
Differential gene and miRNA expression analysis was performed comparing the phthalate-treated samples to the DMSO solvent control using the R package DESeq2 (60), following a slightly modified version of the guidelines of the Omics Data Analysis Framework for regulatory application (R-ODAF) pipeline developed by our group (61, 62). Briefly, a first filtering step (“relevance threshold”) was applied to select the expressed genes/miRNAs by retaining only those whose normalized expression is ≥ 1 count per million (CPM) in at least 75% of the samples in either group (i.e., treatment versus control). To increase statistical power, all doses of a single compound were grouped together and compared to the DMSO control. The RUVg function from the RUVSeq package (63) (k = 2) was used on the genes/miRNAs passing the relevance threshold filter to remove unwanted variation. Differential expression analysis on the expressed genes/miRNAs was then performed setting the FDR to 0.01. The resulting differentially expressed (DE) genes/miRNAs/snoRNAs were subjected to an additional filtering step (“spurious spikes”) to identify those cases in which a very high expression value in only one replicate in a group is responsible for a certain gene/miRNA/snoRNA to be differentially expressed. To this end, the following formula was applied to every DE gene/miRNA/snoRNA for both treatment and control groups: read count gene/miRNAitotal read count gene/miRNAi in groupj<1.4×(number of replicates in groupj)−0.66, where i refers to any gene/miRNA/snoRNA, and j refers to either the treatment or control group. The expression of such genes/miRNAs/snoRNAs was manually checked in all replicates to determine whether a gene that failed this spurious spike filter was indeed a technical artifact or could instead be biologically relevant.
Gene set enrichment analysis
GSEA was performed using the R package ReactomePA (v1.40.0) (64) and Reactome as a database (65) using the DESeq2 “stat” value for gene ranking. For significance, we set a q-value threshold of 0.05.
Differential accessibility analysis
The alignment files (.bam) output by the PEPATAC pipeline were shifted with the deepTools (v3.5.1) (66) utility alignmentSieve to account for the Tn5 transposase duplication at the cut site. To identify differentially accessible (DA) regions, we used a sliding window approach with the R package csaw (v1.32.0) (67) and a modified version of the script made available by Sheikh and Blais on bioRxiv (68). For quantification, we used the five-prime reads, a sliding window of 50 bp without overlap, and a minimum number of 50 counts for a window to be retained. To calculate the background, we binned the genome into 10 kb bins. To distinguish the signal from the background, we compared each region against the global background and set a fold change compared to the background to 3. The differential accessibility analysis was performed with the R package edgeR (v3.4.0) (69), and we performed batch correction using the RUVs function (k = 5) from the RUV-Seq package (v1.32.0) (63). As multiple testing corrections should be performed on regions and not windows (70), we merged the regions identified as “signal” that are at most 500 bp apart, reaching up to a maximum merged region width of 5 kbp, and performed multiple-testing correction using the Benjamini–Hochberg method (FDR = 0.01).
Identified differentially accessible regions were annotated with HOMER (v3.13) (71) (genome version hg38). Regulatory region annotations were retrieved from the ENCODE Candidate Cis-Regulatory Elements (cCREs) registry (72). Coverage tracks were normalized using BeCorrect (v1.1.0) (73) and visuals extracted from the Integrative Genomics Viewer (IGV) (v2.13.2) (74).
ResultsRNA-Seq results
To investigate the effect of phthalates on the transcriptome, we generated RNA-Seq data from mESC-derived thyroid organoids exposed to four phthalates using in vivo relevant concentrations (1–10–100 nM to 1–10 μM DEHP, DIDP, DINP; 2–20–200 nM to 2–20 μM DnOP) for 24 h. A schematic representation of the exposure regimen is shown in Supplementary Figure S1. In the following paragraphs, we provide some dataset quality control (QC) metrics followed by the results of gene expression analysis.
RNA-Seq QC and outlier identification
Combo-Seq libraries had a median of 51.8 million (M) reads per sample (min = 17.7, max = 92.4 M) (Supplementary Figure S2A), with a median of 97% of sequenced reads with a quality score of 30 or more (min = 96.4%, max = 97.3%) (Supplementary Figure S2B). The median number of reads mapped to mRNAs was 41.8 M (min = 14.3, max = 78.8 M) and to miRNAs was 1.30 M (min = 0.16, max = 2.2 M) (Supplementary Figure S2C). As explained in the Methods section, we performed 16 cycles for some RNA-Seq libraries to reach the required concentration for sequencing (1.6 nM) (Supplementary Table S1). In consequence, this increased the percentage of snoRNA-mapping reads (Supplementary Figure S3A). As the read count of the protein-coding genes would be underestimated during DESeq2 normalization, we removed the snoRNA-mapping reads from the main dataset and performed the analysis of snoRNA genes separately. Boxplot of mapped read distributions per gene biotype after snoRNA removal revealed that one DMSO control replicate (DMSO_ctrl_1) was a clear outlier in multiple biotypes (Supplementary Figure S3B). As the DMSO samples would be used as a control for all comparisons, this outlier would have had a major impact in all downstream analyses and, importantly, in the most important biotypes (“protein coding” and “miRNA” in particular”). Although other samples were flagged as outliers in other biotypes (“processed pseudogene” or “rRNA”), this was less consistent and did not warrant further samples removal.
MaSigPro analysis
The MaSigPro R package, initially developed to identify changes in gene expression along a time series, can also be used to analyze the evolution of the gene expression level across a dose range exposure. Next, we investigated using MaSigPro whether some genes would show dose regulation across our six doses (untreated plus five doses). We then allowed the significant genes to be grouped into nine clusters, which include the genes that have a similar trend in change in expression over the dose series. For every compound, we observed some clusters with a nonmonotonic dose–response curve (DEHP: clusters 1 and 3–5; DIDP: clusters 1 and 3–9; DINP: clusters 1, 3–5, and 9; DnOP: clusters 1, 3, 6, and 8) (Figure 1). The genes belonging to the various clusters are reported in Supplementary Table S2. In those indicated clusters, the highest dose (10 μM or 20 μM) was consistently shown to be different from the other four. In this cellular system, this dose could be used for phthalates to derive a point of departure (PoD) metric, which in the toxicology field represents a dose at which a biological response is first observed and can be used to make extrapolations for risk assessment (75).
MaSigPro analysis of the gene expression over the dose series employed in the study (0–1–10–100 nM to 1–10 μM for DEHP, DIDP, and DINP and 0–2–20–200 nM to 2–20 μM for DnOP). The curves for each compound are color-coded, and the compound name is reported in every plot. Each gene was analyzed to fit up to a cubic polynomial regression model. Fits that passed multiple testing corrections (FDR = 0.05) were selected and clustered in nine groups using hierarchical clustering. If a cluster comprises only one gene, the gene name is indicated on top of the plot; otherwise, the cluster name is reported. The x-axis reports the dose range used (nM, nanomoles per liter; μM, micromoles per liter). The dots represent the expression values of each replicate (or the average of each gene if a cluster comprises more genes), and the dotted line shows the fit. The genes belonging to the various clusters are reported in Supplementary Table S2.
Differentially expressed genes and miRNAs
Considering the divergent nature of the highest dose compared to the other four, we decided to exclude it from the differential expression analysis. Given that for each dose we had triplicates or duplicates, by consolidating all the remaining doses together and comparing them to the solvent control, we aimed at increasing the statistical power and detecting gene, miRNA, and snoRNA expression alterations specifically attributable to phthalate treatment. By doing so, we could focus on identifying changes at the compound level while accounting for the different responses observed with the highest dose.
Differential expression analysis revealed how all the treatments had moderate effects on the cells in terms of the number of differentially expressed genes (DEGs), miRNAs, and snoRNAs (Figure 2): the number of DEGs compared to the control was 5, 5, 10, and 49 for DEHP, DIDP, DINP, and DnOP, respectively (FDR< 0.01) (Table 1). Only DIDP treatment influenced miRNA expression, with mmu-miR-143-3p being downregulated. No effect was observed on snoRNA expression or on thyroid markers (Supplementary Figure S4).
Volcano plots of the expressed genes (A–D), miRNAs (E–H), and snoRNAs (I–L) in each phthalate vs. DMSO control. Every dot represents a gene. Elements not differentially expressed (not significant – n.s.) are in black. The false discovery rate (FDR) threshold of 0.01 is indicated as a dotted line on the y-axis. Upregulated elements (“FDR < 0.01, log2(fold change) ≶ 0”) are indicated in red, and downregulated ones (“FDR < 0.01, log2(fold change)< 0”) in cyan. The genes above this line colored in black are the ones that fail to pass the “spurious spike” filter as described in the Methods section.
List of differentially expressed genes (DEGs) and miRNAs (DE miRNAs) in every phthalate vs. DMSO comparison (FDR< 0.01).
DEHP
DIDP
DINP
DnOP
DEGs
Gpd1
Plekha3
Acaa2
Acaa2
Hsd17b10
Plekha3
Tmem80
Ing5
Cxcl14
Cops5
App
Hspa1b
Pole2
Trmt61b
Myh14
Zgpat
Idh3g
Arhgef10l
Idh3g
Ppp1r7
Tspan1
Acaa2
Rpl38-ps2
Ing5
Ccnd2
Ifitm3
Ptp4a1
Ube2g2
Gm15516
Gm10323
Mid1
Ccni
Ing5
Rab5a
Uqcrc2
Rab5a
Cops5
Mid1
Rars
Vps25
Srrm1
Dnajc8
Mtpn
Rnf128
Wasf2
Tnrc6b
Exoc6
Mzt1
Rpl19-ps11
Zc3h11a
Ttc32
Fbp2
Npm1
Shroom3
Zfand1
Zfp960
Gm21293
Nrip1
Sox4
Zfp330
Gm5641
Nucb2
Tjp1
Zfp960
Hbb-bs
Phlda1
Tmed7
Zswim6
Herpud2
DE miRNAs
mmu-miR-143-3p
Downregulated genes are reported in blue text, and upregulated ones in red. DEGs that appear in more than one comparison are in bold.
Interestingly, despite the weak effects on gene expression, the inhibitor of the Ing5 gene was consistently upregulated in three out of four treatments (FDR: DEHP vs. DMSO = 4.44 * 10−3, DIDP vs. DMSO = 0.14, DINP vs. DMSO = 1.37 * 10−3, DnOP vs. DMSO = 1.23 * 10−3) (Figure 3).
Ing5 normalized expression in each phthalate and DMSO control sample. The different doses are reported in the legend. The darkest shade corresponds to the highest concentration (10 μM or 20 μM), while the lightest shades corresponds to the lowest (1 nM or 2 nM). The DMSO control samples are reported in black.
Other genes differentially expressed in more than one condition were identified: Acaa2 (DEHP, DINP, and DnOP vs. DMSO), Plekha3 (DIDP and DnOP vs. DMSO), and five genes (Cops2, Idh3g, Mid1, Rab5a, Zpf960) dysregulated in DINP and DnOP vs. DMSO.
Gene set enrichment analysis
GSEA was performed using the Reactome database with a q-value threshold of 0.05. We identified 123 enriched pathways in the DEHP vs. DMSO comparison, 79 in DIDP vs. DMSO, 173 in DINP vs. DMSO, and 311 in DnOP vs. DMSO (Supplementary File 1; Supplementary Figure S5). In both DEHP and DIDP vs. DMSO comparisons, most pathways were enriched in the control (normalized enrichment score NES< 0). Conversely, we observed a balance between pathways enriched in the treatment (NES > 0) and in the control in the DINP and DnOP vs. DMSO comparisons. To identify common effects across the treatments, we focused on the pathways that appeared in all comparisons (Figure 4), thus retrieving 23 terms, one enriched in the treatment and 22 enriched in the DMSO control.
Results of gene set enrichment analysis (GSEA) on Reactome pathways. The reported pathways have a q-value of< 0.05 and appear in all four phthalates vs. DMSO comparisons.
Among the selected pathways with NES< 0, we identified several terms related to signal transduction and extracellular matrix (ECM) organization. The only term with NES > 0 is “fatty acid metabolism”. In Supplementary Figure S6, the terms are reported with their respective position in the Reactome term hierarchy for a better understanding of their relationships.
ATAC-Seq results
As explained in the previous section, Ing5 upregulation was observed in three phthalate exposures out of four. ING5 is a component of the histone acetyltransferase complexes HBO1-JADE, which mediates histone H4 acetylation in vivo, and HBO1-BRPF1 and MOZ/MORF, which mediate histone H3 acetylation (43, 44). For this reason, we investigated whether phthalate treatment could have an impact on the chromatin status with ATAC-Seq. To this end, we selected two of the four phthalates, DEHP and DINP, and the highest dose included in the differential expression analysis (1 μM). The exposure was increased to 5 days to allow time for any chromatin rearrangements, if any, to take place, accounting for any delay between gene upregulation of Ing5 and an actual observable effect on the epigenome. For ATAC-Seq library preparation, a viability of at least 90% was required. Unfortunately, we were not able to recover enough cells from our thyroid follicles model with this viability. For this reason, we selected the human epithelial thyroid cell line Nthy-ori 3-1 (Supplementary Figure S7). In the next sections, some quality control metrics of the ATAC-Seq libraries and the results of the differential accessibility analysis are reported.
ATAC-Seq QC
ATAC-Seq libraries had a median of 82.2 M reads per sample (min = 26.5, max = 237.9 M) (Supplementary Figure S8A), with a median of 77.36% (min = 75.46%, max = 79.08%) of sequenced reads being successfully aligned to the GRCh38 nuclear genome (Supplementary Figure S8B). The transcription start site (TSS) enrichment score had a distribution between 10.7 and 19.5 (median: 18.1) (Supplementary Figure S8C). The distribution of nucleosome-free regions (NFR), mono-, di-, tri-, or poly-nucleosome regions, was consistent across samples (Supplementary Figure S8D), as were the read length distribution profiles typical of ATAC-Seq libraries (Supplementary Figure S8E). The library complexity metrics were within the accepted values recommended by the Encode Project (Supplementary Figure S8F) (76).
Differential accessibility analysis by ATAC-Seq
We identified 111,133 genomic regions when comparing DEHP-treated and DMSO samples, and 118,855 regions in the DINP vs. DMSO comparison were tested for differential accessibility. In both treatments, we observed a general increase in accessibility compared to the control, but none of the regions passed multiple-testing correction (Figure 5).
Number of differentially accessible regions with p-value of< 0.001 or 0.01 and |log(fold change)|>1 in the DEHP- or DINP-treated samples vs. DMSO control. After multiple-testing correction, none of the regions had an FDR of< 0.01.
We annotated the regions with a p-value of< 0.01 with HOMER to identify the closest gene to each region and looked for overlap between DEHP vs. DMSO and DINP vs. DMSO. We identified four regions with decreased accessibility and 17 with increased accessibility that overlap with regulatory regions (Supplementary Table S3).
We looked at which of these regions fall within the TSS or transcription termination site (TTS). We identified four regions, all with increased accessibility: two located at the TSS of the cell division cycle associated 2 (CDCA2) and complexin 1 (CPLX1) genes, and two located at the TTS of cyclin I (CCNI) and dystrobrevin-binding protein 1 (DTNBV1). However, when inspecting the normalized read coverage on the IGV, there did not seem to be a clear increase in accessibility compared to the control (Supplementary Table S4; Supplementary File 2).
Discussion
In this work, we analyzed the response of mouse embryonic stem cell-derived thyroid follicles after exposure to the phthalates DEHP, DIDP, DINP, and DnOP in a range of concentrations from 1 nM to 10 μM (2 nM to 20 μM for DnOP) for 24 h. The low dose range was selected to reflect the low daily intake of phthalates measured in the general population (77) and the relatively short exposure time to detect the initial reaction to phthalate exposure by identifying the early changes in the transcriptome. In this way, we aimed to test whether our model would be able to capture the molecular initiating event (MIE) of these phthalates, which would then be followed by the key events (KEs) to ultimately lead to an adverse outcome (integrated in the concept of an adverse outcome pathway (AOP)) (78). We performed RNA-Seq analysis and simultaneously analyzed both mRNA and small RNAs from the same samples. The dose series analysis showed how most of the identified genes either increased or decreased sharply in expression at the highest dose, setting it apart from the others and possibly indicating it as a dose to determine a PoD for those genes, which is used in toxicology to establish a threshold dose for risk assessment (75, 79).
The compounds showed a modest effect on the cells at the time and doses of exposure in terms of the number of differentially expressed genes and miRNAs, while no effect was observed on snoRNA expression. DIDP was the only compound where a microRNA (mmu-miR-143-3p) was downregulated. This microRNA, together with mmu-miR-143-5p, has been observed to be downregulated in several cancers and is thought to have tumor-suppressing activity and be a negative regulator of cell proliferation (80–82). Despite the low number of DEGs, we observed a partial overlap across treatments (Acaa2 and Plekha3 in three treatments, Cops2, Idh3g, Mid1, Rab5a, and Zpf960 in two treatments). It is possible that the higher number of DEGs in the DnOP vs. DMSO comparison could be explained by the doses used being twice as high as the other phthalates, though still within the same order of magnitude. Acetyl-CoA acyltransferase 2 (Acaa2) is one of the enzymes that catalyzes the last step of the mitochondrial beta-oxidation pathway. Pleckstrin homology domain containing A3 (Plekha3) is involved in the regulation of vesicular cargo transport from the trans-Golgi network to the plasma membrane and is predicted to be involved in ceramide transport and intermembrane lipid transfer (83, 84). Cops2 is a member of the COP9 signalosome complex (CSN), which is involved in decreasing the ubiquitin ligase activity of the SCF-type E3 ligase complexes. Idh3g is an enzyme that takes part in the Krebs cycle and performs the decarboxylation of isocitrate into alpha-ketoglutarate. Midline 1 (Mid1) is likely involved in the formation of multiprotein structures acting as anchor points to microtubules. It has also E3 ubiquitin ligase activity toward the protein Igbp1, promoting its degradation. Rab5a is a member of the RAS oncogene family and is a small GTPase that, in its active form, recruits proteins responsible for vesicle formation, movement, tethering, and fusion (83).
Via GSEA, we looked for enriched pathways shared by the four treatments to try and identify effects that could be attributed to the “phthalate” EDC class. Only the pathway fatty acid metabolism was enriched in all treatments. Interestingly, phthalates have been shown to increase the metabolism of fatty acids not only in the liver (85, 86) but also in nonliver tissue such as cardiomyocytes, where increased use of fatty acids for energy production was suggested (87). To our knowledge, our analysis is the first observation of the conservation of these mechanisms in an in vitro thyroid model. Additionally, among the pathways downregulated in the treatment groups, we found many related to cell–ECM organization and receptor signaling. It is also interesting to note that despite the low number of DE genes due to the low doses used, were still able to detect relevant enriched pathways using GSEA.
Taken together, the results of differential gene expression gene analysis and GSEA seem to point to an effect of phthalates on energy production, with genes involved in cellular respiration being dysregulated and lipid metabolism increasing.
Ing5 was upregulated in three treatments (DEHP, DINP, DnOP) compared to the control. The ING family comprises five genes (ING1 to ING5), which have a role in cell cycle regulation and cell proliferation by interacting with several partners, such as p53, p300, and histone acetylation complexes (43, 88). ING5 is a tumor-suppressor gene that is downregulated in several types of cancer, including thyroid (89), colorectal (90), breast (91), and lung (92). Its protein is a component of the histone acetyltransferase HBO1-JADE, which acetylates histone H4 at Lys residues 5, 8, and 12 (H4K5ac, H4K8ac, H4K12ac), MOZ/MORPH, which performs histone H3 acetylation, and HBO1-BRPF (H3K14ac) (44, 93).
As a consequence of Ing5 overexpression, we hypothesized that phthalate treatment could have an impact on chromatin status. For this reason, we exposed the human thyroid follicular epithelial cell line Nthy-ori 3-1 to 1 μM of DEHP or DINP for 5 days and analyzed the genomic accessibility by ATAC-Seq. We reasoned that, since Ing5 is not a thyroid-specific gene and its expression is not limited to the thyroid, we would be able to observe changes in a different cell model since we would be investigating a general phthalate mechanism rather than a model-specific response. Differential accessibility analysis resulted in a general increase in accessibility in the treatment group, but none of the identified regions passed multiple-testing correction. Among the regions with a p-value of< 0.01, we identified four common ones with increased accessibility in the DEHP vs. DMSO and DINP vs. DMSO comparison localized on regulatory regions in the TSS or TTS. However, the signal did not seem to reflect a real change in accessibility.
In this work, we showed that even with a stem cell-derived in vitro thyroid model exposed to a range of low, biologically relevant concentrations of four phthalates, we were able to detect some of the effects that have been previously reported in vivo. Our analysis demonstrates that it is not necessary to use cytotoxic doses in toxicological experiments to obtain observable results and that low-dose exposure can be analyzed without lowering the statistical stringency. We are convinced that 3D in vitro systems, such as organoids, can be a valid alternative to animal studies even for EDCs, provided that enough datasets are generated to allow regulators to infer risk assessment thresholds.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ebi.ac.uk/arrayexpress/, E-MTAB-12830.
Author contributions
Conceptualization: MN and FC. Methodology: MN, MR, DH, and DC. Formal analysis and investigation: MN. Writing—original draft preparation: MN and FC. Writing—review and editing: all coauthors. Funding acquisition: FC, SC, SG, and LM. Resources: FC, SC, SG, and LM. Supervision: FC.
Funding
The SCREENED project has received funding from the European Union’s Horizon 2020 research and innovation program under Grant Agreement No. 825745.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher's note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fendo.2023.1200211/full#supplementary-material
ReferencesBabichMA. Overview of phthalates toxicity. Bethesda, MD 20814, US: US Consumer Product Safety Commission, Bethesda, MD 20814. (2010) Available at: https://www.cpsc.gov/s3fs-public/phthalover.pdf.CarlsonKR. Toxicity review of Di(2-ethylhexyl) Phthalate (DEHP). Bethesda, MD 20814, US: US Consumer Product Safety Commission, Bethesda. (2010) Available at: https://www.cpsc.gov/s3fs-public/ToxicityReviewOfDEHP.pdf.U.S EPA. Phthalates action plan2012 . Available at: https://www.epa.gov/sites/default/files/2015-09/documents/phthalates_actionplan_revised_2012-03-14.pdf.HollandM. Socio-economic assessment of phthalates. Paris: OECD Publishing (2018). doi: 10.1787/a38a0e34-en.Phthalates NRCUCotHRo. 2 - phthalate exposure assessment in humans. In: Phthalates and cumulative risk assessment: The tasks ahead. Washington, DC, US: Washington (DC): National Academies Press (US (2008). Available at: https://www.ncbi.nlm.nih.gov/books/NBK215044/.WittassekMKochHMAngererJBruningT. Assessing exposure to phthalates - the human biomonitoring approach. Mol Nutr Food Res (2011) 55(1):7–31. doi: 10.1002/mnfr.201000121Versar Inc. Review of exposure data and assessments for select dialkyl ortho-phthalatesFebruary 2010 . Available at: https://www.cpsc.gov/s3fs-public/pthalexp.pdf.(CDC) CfDCaP. Fourth national report on human exposure to environmental chemicals. Washington, DC, US: U.S. Department of Health and Human Services. (2019). doi: 10.15620/cdc75822.HainesDASaravanabhavanGWerryKKhouryC. An overview of human biomonitoring of environmental chemicals in the Canadian Health Measures Survey: 2007-2019. Int J Hyg Environ Health (2017) 220(2 Pt A):13–28. doi: 10.1016/j.ijheh.2016.08.002WittassekMAngererJ. Phthalates: metabolism and exposure. Int J Androl. (2008) 31(2):131–8. doi: 10.1111/j.1365-2605.2007.00837.xKochHMAngererJ. Chapter 3A. Phthalates: biomarkers and human biomonitoring. In: Biomarkers and human biomonitoring. Cambridge: Royal Society of Chemistry. (2011) 179–233. doi: 10.1039/9781849733373-00179NetSSempereRDelmontAPaluselliAOuddaneB. Occurrence, fate, behavior and ecotoxicological state of phthalates in different environmental matrices. Environ Sci Technol (2015) 49(7):4019–35. doi: 10.1021/es505233bConsumer product safety improvement act of 2008, title I—Children's product safety, sec. 108, public law (2008). Available at: https://www.cpsc.gov/s3fs-public/pdfs/blk_pdf_cpsia.pdf.Official Journal of the European Union. Document L:2005:344:TOC, 27.12.2005, p. 43. Available at: https://eur-lex.europa.eu/legal-content/EN/ALL/?uri=OJ:L:2005:344:TOCBabichMAOsterhoutCA. Toxicity review of Diisononyl Phthalate (DINP). Bethesda, MD 20814, US: US Consumer Product Safety Commission, Bethesda, MD (2010). p. 20814. Available at: https://www.cpsc.gov/s3fs-public/ToxicityReviewOfDINP.pdf.CarlsonKR. Toxicity review of Di-n-Octyl Phthalate (DnOP). Bethesda, MD 20814, US: US Consumer Product Safety Commission, Bethesda, MD (2010). p. 20814. Available at: https://www.cpsc.gov/s3fs-public/pdfs/blk_media_toxicityDNOP.pdf.OsterhoutCA. Toxicity review of di(isodecyl) phthalate. Bethesda, MD 20814, US: US Consumer Product Safety Commission, Bethesda, MD (2010). p. 20814. Available at: https://www.cpsc.gov/s3fs-public/toxicityDIDP.pdf.Phthalates NRCUCotHRo. 3 - toxicity assessment. Phthalates and cumulative risk assessment: The tasks ahead. Washington, DC, US: National Academies Press (US (2008). Available at: https://www.ncbi.nlm.nih.gov/books/NBK215030/.FletcherEJSantacruz-MárquezRMourikesVENeffAMLawsMJFlawsJA. Effects of phthalate mixtures on ovarian folliculogenesis and steroidogenesis. Toxics. (2022) 10(5):251. doi: 10.3390/toxics10050251OECD. Detailed review paper on the state of the science on novel in vitro and in vivo screening and testing methods and endpoints for evaluating endocrine disruptors. Paris: OECD Series on Testing and Assessment, No 178, OECD Publishing, Paris (2012). doi: 10.1787/9789264221352-enRadkeEGBraunJMMeekerJDCooperGS. Phthalate exposure and male reproductive outcomes: A systematic review of the human epidemiological evidence. Environ Int (2018) 121(Pt 1):764–93. doi: 10.1016/j.envint.2018.07.029GoreACChappellVAFentonSEFlawsJANadalAPrinsGS. EDC-2: The endocrine society's second scientific statement on endocrine-disrupting chemicals. Endocr Rev (2015) 36(6):E1–E150. doi: 10.1210/er.2015-1010LindPMLindL. Endocrine-disrupting chemicals and risk of diabetes: an evidence-based review. Diabetologia. (2018) 61(7):1495–502. doi: 10.1007/s00125-018-4621-3TuculinaMJPerleaPGheorghitaMCumpataCNDascaluITTurcuA. Diabetes mellitus: Plasticizers and nanomaterials acting as endocrine-disrupting chemicals (Review). Exp Ther Med (2022) 23(4):288. doi: 10.3892/etm.2022.11217ChengSYLeonardJLDavisPJ. Molecular aspects of thyroid hormone actions. Endocr Rev (2010) 31(2):139–70. doi: 10.1210/er.2009-0007Hiller-SturmhöfelSBartkeA. The endocrine system: an overview. Alcohol Health Res World. (1998) 22(3):153–64.BenvengaSTuccariGIeniAVitaR. Thyroid gland: anatomy and physiology. Encyclopedia Endocrine Diseases. (2018) 4:382–90. doi: 10.1016/B978-0-12-801238-3.96022-7CarvalhoDPDupuyC. Thyroid hormone biosynthesis and release. Mol Cell Endocrinol (2017) 458:6–15. doi: 10.1016/j.mce.2017.01.038BereketogluCPradhanA. Plasticizers: negative impacts on the thyroid hormone system. Environ Sci pollut Res Int (2022) 29(26):38912–27. doi: 10.1007/s11356-022-19594-0GhisariMBonefeld-JorgensenEC. Effects of plasticizers and their mixtures on estrogen receptor and thyroid hormone functions. Toxicol Lett (2009) 189(1):67–77. doi: 10.1016/j.toxlet.2009.05.004KimMJKimHHSongYSKimOHChoiKKimS. DEHP down-regulates tshr gene expression in rat thyroid tissues and FRTL-5 rat thyrocytes: A potential mechanism of thyroid disruption. Endocrinol Metab (Seoul). (2021) 36(2):447–54. doi: 10.3803/EnM.2020.920WuHZhangWZhangYKangZMiaoXNaX. Novel insights into di−(2−ethylhexyl)phthalate activation: Implications for the hypothalamus−pituitary−thyroid axis. Mol Med Rep (2021) 23(4):290. doi: 10.3892/mmr.2021.11930ManikkamMTraceyRGuerrero-BosagnaCSkinnerMK. Plastics derived endocrine disruptors (BPA, DEHP and DBP) induce epigenetic transgenerational inheritance of obesity, reproductive disease and sperm epimutations. PLoS One (2013) 8(1):e55387. doi: 10.1371/journal.pone.0055387Martinez-ArguellesDBPapadopoulosV. Prenatal phthalate exposure: epigenetic changes leading to lifelong impact on steroid formation. Andrology. (2016) 4(4):573–84. doi: 10.1111/andr.12175GuidaNLaudatiGGalganiMSantopaoloMMontuoriPTriassiM. Histone deacetylase 4 promotes ubiquitin-dependent proteasomal degradation of Sp3 in SH-SY5Y cells treated with di(2-ethylhexyl)phthalate (DEHP), determining neuronal death. Toxicol Appl Pharmacol (2014) 280(1):190–8. doi: 10.1016/j.taap.2014.07.014JacobsMNMarczyloELGuerrero-BosagnaCRüeggJ. Marked for life: Epigenetic effects of endocrine disrupting chemicals. Annu Rev Environ Resources. (2017) 42(1):105–60. doi: 10.1146/annurev-environ-102016-061111OgundipeVMLPlukkerJTMLinksTPCoppesRP. Thyroid gland organoids: Current models and insights for application in tissue engineering. Tissue Eng Part A. (2022) 28(11-12):500–10. doi: 10.1089/ten.tea.2021.0221LiYTangPCaiSPengJHuaG. Organoid based personalized medicine: from bench to bedside. Cell Regen. (2020) 9(1):21. doi: 10.1186/s13619-020-00059-zSondorpLHJOgundipeVMLGroenAHKelderWKemperALinksTP. Patient-derived papillary thyroid cancer organoids for radioactive iodine refractory screening. Cancers (Basel). (2020) 12(11):3212. doi: 10.3390/cancers12113212Caipa GarciaALArltVMPhillipsDH. Organoids for toxicology and genetic toxicology: applications with drugs and prospects for environmental carcinogenesis. Mutagenesis. (2022) 37(2):143–54. doi: 10.1093/mutage/geab023RussellWMSBurchRL. The principles of humane experimental technique. Med J Australia. (1960) 1(13):500–. doi: 10.5694/j.1326-5377.1960.tb73127.xAntonicaFKasprzykDFOpitzRIacovinoMLiaoXHDumitrescuAM. Generation of functional thyroid from embryonic stem cells. Nature. (2012) 491(7422):66–71. doi: 10.1038/nature11525DoyonYCayrouCUllahMLandryAJCoteVSelleckW. ING tumor suppressor proteins are critical regulators of chromatin acetylation required for genome expression and perpetuation. Mol Cell (2006) 21(1):51–64. doi: 10.1016/j.molcel.2005.12.007LalondeMEAvvakumovNGlassKCJoncasFHSaksoukNHollidayM. Exchange of associated factors directs a switch in HBO1 acetyltransferase histone tail specificity. Genes Dev (2013) 27(18):2009–24. doi: 10.1101/gad.223396.113RomittiMEskiSEFonsecaBFGillotayPSinghSPCostagliolaS. Single-cell trajectory inference guided enhancement of thyroid maturation in vitro using TGF-beta inhibition. Front Endocrinology. (2021) 12:2021.657195. doi: 10.3389/fendo.2021.657195CorcesMRTrevinoAEHamiltonEGGreensidePGSinnott-ArmstrongNAVesunaS. An improved ATAC-seq protocol reduces background and enables interrogation of frozen tissues. Nat Methods (2017) 14(10):959–62. doi: 10.1038/nmeth.4396NazzariMHauserDvan HerwijnenMRomittiMCarvalhoDJKipAM. CODA: a combo-Seq data analysis workflow. Brief Bioinform (2022) 24(1):1–14. doi: 10.1093/bib/bbac582MartinM. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnetjournal. (2011) 17(1):10–12. doi: 10.14806/ej.17.1.200PerkinElmer Inc. NEXTFLEX® Combo-seq analysis guidelines 2020 . Available at: https://perkinelmer-appliedgenomics.com/wp-content/uploads/2020/06/NOVA-5139-AG_v01_NEXTFLEX-Combo-seq-Analysis-Guideline.pdf.LiBDeweyCN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinf (2011) 12:323. doi: 10.1186/1471-2105-12-323LiB. rsem-prepare-reference documentation page . Available at: https://deweylab.github.io/RSEM/rsem-prepare-reference.html.PatilAHHalushkaMK. miRge3.0: a comprehensive microRNA and tRF sequencing analysis pipeline. NAR Genom Bioinform (2021) 3(3):lqab068. doi: 10.1093/nargab/lqab068SmithJPCorcesMRXuJReuterVPChangHYSheffieldNC. PEPATAC: an optimized pipeline for ATAC-seq data analysis with serial alignments. NAR Genom Bioinform (2021) 3(4):lqab101. doi: 10.1093/nargab/lqab101LangmeadBSalzbergSL. Fast gapped-read alignment with Bowtie 2. Nat Methods (2012) 9(4):357–9. doi: 10.1038/nmeth.1923DanecekPBonfieldJKLiddleJMarshallJOhanVPollardMO. Twelve years of SAMtools and BCFtools. Gigascience. (2021) 10(2):1–4. doi: 10.1093/gigascience/giab008Core TeamR. R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing (2021).DurinckSSpellmanPTBirneyEHuberW. Mapping identifiers for the integration of genomic datasets with the R/Bioconductor package biomaRt. Nat Protoc (2009) 4(8):1184–91. doi: 10.1038/nprot.2009.97ConesaANuedaMJ. maSigPro: significant gene expression profile differences in time course gene expression data. R package version. (2020)ConesaANuedaMJ. maSigPro user’s guide4 september (2017). Available at: https://www.bioconductor.org/packages/release/bioc/vignettes/maSigPro/inst/doc/maSigProUsersGuide.pdf.LoveMIHuberWAndersS. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol (2014) 15(12):550. doi: 10.1186/s13059-014-0550-8CEFIC C4 team. Omics Data Analysis Framework for Regulatory application (R-ODAF) (2021). Available at: https://github.com/R-ODAF/Main.VerheijenMCMeierMJAsensioJOGantTWTongWYaukCL. R-ODAF: Omics data analysis framework for regulatory application. Regul Toxicol Pharmacol (2022) 131:105143. doi: 10.1016/j.yrtph.2022.105143RissoDNgaiJSpeedTPDudoitS. NorMalization of RNA-seq data using factor analysis of control genes or samples. Nat Biotechnol (2014) 32(9):896–902. doi: 10.1038/nbt.2931YuGHeQY. ReactomePA: an R/Bioconductor package for reactome pathway analysis and visualization. Mol Biosyst (2016) 12(2):477–9. doi: 10.1039/C5MB00663EGillespieMJassalBStephanRMilacicMRothfelsKSenff-RibeiroA. The reactome pathway knowledgebase 2022. Nucleic Acids Res (2022) 50(D1):D687–D92. doi: 10.1093/nar/gkab1028RamirezFRyanDPGruningBBhardwajVKilpertFRichterAS. deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Res (2016) 44(W1):W160–5. doi: 10.1093/nar/gkw257LunATSmythGK. csaw: a Bioconductor package for differential binding analysis of ChIP-seq data using sliding windows. Nucleic Acids Res (2016) 44(5):e45. doi: 10.1093/nar/gkv1191SheikhAABlaisA. Improved sensitivity and resolution of ATAC-seq differential DNA accessibility analysis. bioRxiv. (2022) 2022:3. doi: 10.1101/2022.03.16.484118RobinsonMDMcCarthyDJSmythGK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. (2010) 26(1):139–40. doi: 10.1093/bioinformatics/btp616LunATLSmythGK. De novo detection of differentially bound regions for ChIP-seq data using peaks and windows: controlling error rates correctly. Nucleic Acids Res (2014) 42(11):e95. doi: 10.1093/nar/gku351HeinzSBennerCSpannNBertolinoELinYCLasloP. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol Cell (2010) 38(4):576–89. doi: 10.1016/j.molcel.2010.05.004AbascalFAcostaRAddlemanNJAdrianJAfzalVAiR. Expanded encyclopaedias of DNA elements in the human and mouse genomes. Nature. (2020) 583(7818):699–710. doi: 10.1038/s41586-020-2493-4GontarzPFuSXingXLiuSMiaoBBazylianskaV. Comparison of differential accessibility analysis strategies for ATAC-seq data. Sci Rep (2020) 10(1):10150. doi: 10.1038/s41598-020-66998-4RobinsonJTThorvaldsdóttirHWincklerWGuttmanMLanderESGetzG. Integrative genomics viewer. Nat Biotechnol (2011) 29(1):24–6. doi: 10.1038/nbt.1754JohnsonGESoeteman-HernandezLGGollapudiBBBodgerOGDearfieldKLHeflichRH. Derivation of point of departure (PoD) estimates in genetic toxicology studies and their potential applications in risk assessment. Environ Mol Mutagen (2014) 55(8):609–23. doi: 10.1002/em.21870Terms and definitions (2023). Available at: https://www.encodeproject.org/data-standards/terms/.FrommeHGruberLSchlummerMWolzGBohmerSAngererJ. Intake of phthalates and di(2-ethylhexyl)adipate: results of the Integrated Exposure Assessment Survey based on duplicate diet samples and biomonitoring data. Environ Int (2007) 33(8):1012–20. doi: 10.1016/j.envint.2007.05.006OECD. ENV/JM/MONO. 6 revised guidance document on developing and assessing adverse outcome pathways 27-jul-2017 (2013). Available at: https://one.oecd.org/document/env/jm/mono(2013)6/en/pdf.SturlaSJ. Point of departure. Chem Res Toxicol (2018) 31(1):2–3. doi: 10.1021/acs.chemrestox.7b00341WadaMGotoYTanakaTOkadaRMoriyaSIdichiT. RNA sequencing-based microRNA expression signature in esophageal squamous cell carcinoma: oncogenic targets by antitumor miR-143-5p and miR-143-3p regulation. J Hum Genet (2020) 65(11):1019–34. doi: 10.1038/s10038-020-0795-xSanadaHSekiNMizunoKMisonoSUchidaAYamadaY. Involvement of Dual Strands of miR-143 (miR-143-5p and miR-143-3p) and Their Target Oncogenes in the Molecular Pathogenesis of Lung Adenocarcinoma. Int J Mol Sci (2019) 20(18):4482. doi: 10.3390/ijms20184482LiDHuJSongHXuHWuCZhaoB. miR-143-3p targeting LIM domain kinase 1 suppresses the progression of triple-negative breast cancer cells. Am J Transl Res (2017) 9(5):2276–85.StelzerGRosenNPlaschkesIZimmermanSTwikMFishilevichS. The geneCards suite: From gene data mining to disease genome sequence analyses. Curr Protoc Bioinf (2016) 54:1 30 1–1 3. doi: 10.1002/cpbi.5GodiADi CampliAKonstantakopoulosADi TullioGAlessiDRKularGS. FAPPs control Golgi-to-cell-surface membrane traffic by binding to ARF and PtdIns(4)P. Nat Cell Biol (2004) 6(5):393–404. doi: 10.1038/ncb1119HintonRHMitchellFEMannAChescoeDPriceSCNunnA. Effects of phthalic acid esters on the liver and thyroid. Environ Health Perspect (1986) 70:195–210. doi: 10.1289/ehp.8670195HowarthJAPriceSCDobrotaMKentishPAHintonRH. Effects on male rats of di-(2-ethylhexyl) phthalate and di-n-hexylphthalate administered alone or in combination. Toxicol Lett (2001) 121(1):35–43. doi: 10.1016/S0378-4274(01)00313-7PosnackNGSwiftLMKayMWLeeNHSarvazyanN. Phthalate exposure changes the metabolic profile of cardiac muscle cells. Environ Health Perspect (2012) 120(9):1243–51. doi: 10.1289/ehp.1205056ShisekiMNagashimaMPedeuxRMKitahama-ShisekiMMiuraKOkamuraS. p29ING4 and p28ING5 bind to p53 and p300, and enhance p53 activity. Cancer Res (2003) 63(10):2373–8.GaoWHanJ. Overexpression of ING5 inhibits HGF-induced proliferation, invasion and EMT in thyroid cancer cells via regulation of the c-Met/PI3K/Akt signaling pathway. BioMed Pharmacother. (2018) 98:265–70. doi: 10.1016/j.biopha.2017.12.045XinHWangCChiYLiuZ. MicroRNA-196b-5p promotes Malignant progression of colorectal cancer by targeting ING5. Cancer Cell Int (2020) 20(1):119. doi: 10.1186/s12935-020-01200-3XuJZhaoJJiangMYangLSunMWangH. MiR-193 promotes cell proliferation and invasion by ING5/PI3K/AKT pathway of triple-negative breast cancer. Eur Rev Med Pharmacol Sci (2020) 24(6):3122–9. doi: 10.26355/eurrev_202003_20679LiuXLMengJZhangXTLiangXHZhangFZhaoGR. ING5 inhibits lung cancer invasion and epithelial-mesenchymal transition by inhibiting the WNT/beta-catenin pathway. Thorac Cancer. (2019) 10(4):848–55. doi: 10.1111/1759-7714.13013HanJLachanceCRickettsMDMcCulloughCEGeraceMBlackBE. The scaffolding protein JADE1 physically links the acetyltransferase subunit HBO1 with its histone H3-H4 substrate. J Biol Chem (2018) 293(12):4498–509. doi: 10.1074/jbc.RA117.000677