<figure><img src="image_1.png" /></figure>

# Distinct transcriptional profiles in rat thyroid glands after developmental exposure to three *in vitro* thyroperoxidase inhibiting chemicals

Anna Kjerstine Rosenmai<sup>a,*</sup>, Terje Svingen<sup>a</sup>, Bertrand Evrard<sup>b</sup>, Khanh Hoang Nguyen<sup>a</sup>, Camilla Nielsen<sup>a</sup>, Marta Axelstad<sup>a</sup>, Frédéric Chalmel<sup>b</sup>, Louise Ramhøj<sup>a</sup>

## ARTICLE INFO

**Keywords:**

- Thyroperoxidase
- Enzyme inhibition
- Transcriptomics
- Thyroid hormone system disruption
- Toxicity profiling
- Thyroxine
- Thyroid stimulating hormone

## ABSTRACT

Thyroperoxidase (TPO) is central in thyroid hormone (TH) synthesis and inhibition can lead to TH deficiency. Many chemicals can inhibit TPO activity *in vitro*, but how this may manifest in the developing thyroid gland at the molecular level is unclear. Here, we characterized the thyroid gland transcriptome of male rats developmentally exposed to the *in vitro* TPO-inhibitors amitrole, 2-mercaptobenzimidazole (MBI), or cyanamide by use of Bulk-RNA-Barcoding (BRB) and sequencing. Amitrole exposure caused TH deficiency and 149 differentially expressed genes in the thyroid gland. The effects indicated an activated and growing thyroid gland. MBI caused intermittent changes to serum TH concentrations in a previous study and this was accompanied by 60 differentially expressed genes in the present study. More than half of these were also affected by amitrole, indicating that they could be early effect biomarkers of developmental TH system disruption due to TPO inhibition. Further work to validate the signature is needed, including assessment of substance independency and applicability domain.

## 1. Introduction

Thyroperoxidase (TPO) is a key enzyme in thyroid hormone (TH) biosynthesis. In the thyroid gland TPO facilitates the oxidation of iodide, iodination of thyroglobulin, and couples iodinated tyrosine residues to form THs [1]. Hence, inhibition of TPO may lead to decreased TH production, low serum TH concentrations, and altered TH signaling in target tissues [2,3].

*In vitro* studies of TPO inhibition suggest that the enzyme is relatively promiscuous. Out of approximately 1000 tested substances, 300 could inhibit enzyme activity [4]. Considering the very high number of chemicals on the market, the number of chemicals that could potentially affect TPO activity could also be very high. This uncertainty, coupled with the fact that we do not fully understand how TPO inhibition manifest *in vivo* at the molecular level, raises serious concerns about the potential impact on human health and hampers development of new approach methodologies for thyroid hormone system disruptors.

Although *in vitro* studies suggest that many chemicals have the potential to inhibit TPO, disruption of the TH system *in vivo* is not a given

[4]. Many factors can contribute to a discrepancy between *in vitro* and *in vivo* effects, for instance toxicokinetics and off-target toxicity in complex organisms. This we illustrated in a recent study where three *in vitro* TPO inhibitors led to markedly different profiles in developmentally exposed male rats [5]. In this study, amitrole caused hypothyroidism and TH-mediated developmental neurotoxicity, 2-mercaptobenzimidazole (MBI) caused intermittent changes to TH levels, whereas cyanamide caused no TH system effects [5].

From animal studies we know that inhibition of TPO in the thyroid gland can cause hypothyroidism. This typically manifest as low serum thyroxine (T4) levels, high serum thyroid stimulating hormone (TSH) levels, and increased thyroid gland weight, along with follicular cell hypertrophy and hyperplasia [6-11]. The effects on the gland are thought to primarily arise from activation of the hypothalamic-pituitary-thyroid (HPT) axis, where a decrease in serum TH levels activate the HPT-axis by a negative feedback loop. In turn, activation of the HPT-axis leads to increased secretion of TSH from the pituitary, activating the TSH receptor (TSHR) on the thyroid follicular cells, and ultimately regulating gene expression [2,3,12]. Changes in gene

<footnote><sup>a</sup> National Food Institute, Technical University of Denmark, 2800 Kgs. Lyngby, Denmark</footnote>

<footnote><sup>b</sup> Univ Rennes, Inserm, EHESP, Irset (Institut de Recherche en Santé, Environnement et Travail) - UMR_S 1085, F-35000 Rennes, France</footnote>

<footnote>* Corresponding author at: Research Group for Molecular and Reproductive Toxicology, National Food Institute Technical University of Denmark, Kemitorvet Building 202, DK-2800 Kgs. Lyngby, Denmark.</footnote>

<footnote>E-mail address: akjro@food.dtu.dk (A.K. Rosenmai).</footnote>
expression downstream of the TSHR are, therefore, potential biomarkers of TPO inhibition.

TSH directly or indirectly controls the transcription of several TH synthesis genes including *Slc5a5*(NIS), *Tg*, *Tpo*, *Tshr*, *Slc26a4* (coding for Pendrin), *Duox* [2,13], as well as the expression of transcription factors such as *Pax8*, *Foxe1*, and *Nkx2.1* [13]. Changes in the expression of several of these genes have been reported after pubertal or developmental exposure to thyroid hormone system disruptors in rats, such as TPO or NIS inhibitors [14–16]. Key genes such as *Slc5a5*(NIS), *Tshr*, *Pax8*, *Slc26a4* (Pendrin), *Tpo*, and *Tg* were differentially expressed in exposed animals, including dams, fetuses, and/or pubertally exposed rats. In these studies, the rats presented with at least one of the aforementioned classical effects: e.g. increased thyroid gland weight, decreased serum TH levels, or increased serum TSH levels [14,15]. This could suggest that the adverse effect profile is linked to the changes in TH synthesis-related gene expression. To follow up on this assertion, we aimed to test whether developmental exposure to the three TPO-inhibitors amitrole, MBI and cyanamide would induce transcriptional profiles in the thyroid gland that support the phenotypic effects, or lack thereof. Ultimately, the aim is to identify common denominators of TPO inhibition that can be used as sensitive predictors and identifiers of effects in the developing thyroid gland.

## 2. Materials and methods

The animal study, methods and results on reproductive toxicity, serum hormone analysis by radioimmunoassay (RIA) and Luminex and developmental neurotoxicity was reported in Ramhøj et al. [5]. Below is a brief description of the chemical compounds and the animal study along with material and methods specific to the present study.

### 2.1. Test chemicals

Three test substances were used for *in vivo* testing: 2-mercaptobenzimidazole (MBI) (CAS #583-39-1, Sigma-Aldrich, St. Louis, MO, USA, 98 %), 3-Amino-1H-1,2,4-triazole (amitrole) (CAS #61-82-5, Sigma-Aldrich, >97.0 %), and cyanamide (CAS #420-04-2, Sigma-Aldrich, 99 %). They were selected based on ToxCast data, published literature, and REACH registration dossiers with focus on *in vitro* TPO inhibition potency, *in vivo* thyrotoxic potential, presumed lack of carcinogenicity and neurotoxicity through modes of action other than TPO inhibition, and availability of toxicity studies in pregnant animals. The compound selection criteria is further detailed by Ramhøj et al. [5].

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Experimental design and serum thyroid hormones in PD16 male pups perinatally exposed to MBI, amitrole or cyanamide. A) Pregnant rat dams (n = 4) were exposed by oral gavage from GD7 to PD16. Blood and thyroid glands were collected from 16-day old male pups for serum hormone analysis (n = 4) and transcriptomic analysis (n = 7–8). B) Serum thyroxine (T4) and tri-iodothyronine (T3) concentrations in PD16 male pups. Analysis done by UPLC-HRMS. LOQ (T4 = 1.29 nM, T3 = 1.54 nM) and LOD (T4 = 0.05 nM, T3 = 0.06 nM) indicated with dotted line. n = 4. Data shown as mean ± SD. AMI: amitrole, BRB-seq: bulk-RNA-barcoding and sequencing, bw: body weight, CYA: cyanamide, GD: gestation day, LOD: limit of detection, LOQ: limit of quantification, MBI: 2-mercaptobenzimidazole, PD: pup day, UPLC-HRMS: ultra-performance liquid chromatography – high-resolution mass spectrometry. Parts of figure created with **BioRender.com**.</figcaption></figure>

### 2.2. Animals, exposure, and tissue retrieval

Pregnant Sprague-Dawley (Crl:CD(SD)) rats (breeder Charles River Europe, distributor SCANBUR, Karlslunde, Denmark) were delivered on gestation day (GD) 3 of pregnancy to the BioFacility at the Technical University of Denmark. Nulliparous rats with approximate body weights of 240 ± 30 g were time-mated, with day of plug detection denoted GD1. On GD4 dams were weighed and randomly assigned a treatment group (n = 4). Exposure started on GD7 and continued to tissue collection on pup day (PD) 16 with daily doses of MBI (10 and 20 mg/kg body weight (bw)/day (d)), amitrole (25 and 50 mg/kg bw/d) and cyanamide (7 and 15/11.25 mg/kg bw/d) administered by oral gavage (2 mL/kg bw, corn oil (Sigma-Aldrich)) in the morning (Fig. 1A). A vehicle control group was also included. No treatment was given on the day of birth. From GD16 the cyanamide dose was reduced from 15 to 11.25 mg/kg bw/d due to reduced dam weight gain. The high doses were selected for further analysis as part of this study, namely 10 mg/kg bw/d) MBI (high dose group lost due to severe reproductive toxicity), 50 mg/kg bw/d amitrole and 15/11.25 mg/kg bw/d cyanamide. Body weight, organ weight, and TH levels for dams and offspring (males and females) for all dose groups are reported elsewhere [5].

Animals were kept at environmentally controlled conditions of 55 ± 5 % humidity, 21 ± 1 °C, 50 times airchange/h, and reverse light/dark cycles of 12 h. Acidified tap water and a standard diet of Altromin 1314 (soy and alfalfa-free, Altromin GmbH, Lage, Germany) were given *ad libitum*. Pups were weighed, sexed and counted on the day of birth (PD1). At PD16 two male pups/litter were anesthetized by CO₂/O₂ and killed by decapitation. The same procedure was used to kill dams and remaining pups on PD22.

Thyroid glands (n = 7–8, 1–2 male pups per litter) were dissected from the trachea and frozen in liquid nitrogen. Samples were stored at –80 °C until further analysis. Trunk blood was collected and pooled from three male pups/litter (n = 4). Blood was centrifuged for 10 min. at 4 °C and serum was collected and stored at –80 °C until analysis.

Animal experiments were performed with ethical approval from the Danish Animal Experiments Inspectorate (authorization number 2015-15-0201-00553 C3) and overseen by the Animal Welfare Committee at the National Food Institute to assure compliance with relevant regulations and guidelines for animal experimentation. Further, experiments were done in compliance with the ARRIVE guidelines.

### 2.3. Serum hormone analysis

Thawed samples were vortexed and 50 µL transferred to a light

<figure><img src="image_3.png" /></figure>
blocking micro centrifuge tube (LightSafe tubes, Cole Parmer, Vernon Hills, IL, USA). ¹³C-isotope labeled internal standards of T4, T3, and r-T3 were spiked (500 pg) into the sample followed by the addition of 50 µL of 1 N HCL, 250 µL of Milli Q water, and 150 µL of 0.1 % formic acid in Milli Q water:Acetronitrile (50:50). The mixtures were vortexed and then incubated for 2 h at 37 °C. After the incubation, the samples were cooled, then 250 µL of 0.1 % formic acid in Milli Q water was added and mixed vigorously. The samples were then loaded on an SPE cartridge (Evolute CX 25 mg/1 mL, Biotage, Sweden) preconditioned with 400 µL of methanol and 400 µL of 0.1 % formic acid in Milli Q water. The SPE cartridge was washed with 400 µL of 0.1 % formic acid in Milli Q water then 400 µL of methanol. The thyroid hormones were then eluted into an amber glass tube using 500 µL of 5 % ammonium in methanol twice. The eluate was evaporated under a gentle nitrogen gas stream until dryness, reconstituted in 100 µL of 0.1 % formic acid in Milli Q water, and transferred to an amber glass vial.

T4, T3, and r-T3 were quantified using a liquid chromatography (LC) - high-resolution mass spectrometry (HRMS) system consisting of an ultra-performance liquid chromatography (UPLC) (Vanquish Flex, Thermo Fisher Scientific, Waltham, MA, USA) coupled to an Orbitrap HRMS (Exploris 120, Thermo Fisher Scientific). The analytes were separated using a C18 UPLC column (Accucore Vanquish C18+, 100 mm × 2.1 mm × 1.5 µm, Thermo Fisher Scientific) heated at 50 °C. The mobile phase consists of 0.1 % formic acid in Milli Q water (mobile phase A) and 0.1 % formic acid in acetonitrile (mobile phase B) and run at 0.4 mL/min. The gradient was programmed to start at 5 % B for 2 min, increased to 30 % B at 4.5 min, increased to 38 % B at 5.5 min, slowly increased to 40 % at 6.5 min and kept for 1.5 min, increased to 98 % B at 9 min and remained constant for 3 min, followed by a sharp decrease to 5 % B at 12.1 min and equilibrate for another 1.9 min for a total run time of 14 min. The Orbitrap HRMS was operated in full scan mode (600–800 m/z) at 60000 resolution power, positive ESI mode (3925 V), and vaporization temperature of 350 °C. Accurate mass (with 3 ppm deviation) of the [M + H]⁺ ions of the analytes and their corresponding ¹³C-labeled internal standards were used to quantify them (see supplementary material 1). Calibration curve was constructed within the range of 0.004 ng/mL to 100 ng/mL using at least 5 points. The limit of detection (LOD) and limit of quantification (LOQ) for T3 were estimated at 0.06 and 1.54 nM, respectively, whereas for T4 they were estimated at 0.05 and 1.29 nM, respectively.

Data were evaluated for variance homogeneity and normal distribution by Bartlett's and Shapiro-Wilk test, respectively. Data passed normality test and had homogenous variances. A one-way ANOVA was performed with Dunnett's *post hoc* test to adjust for multiple comparisons between exposed and control groups (significance level set at 0.05).

## 2.4. RNA isolation

RNA was isolated from thyroid glands from the groups exposed to 10 mg/kg bw/d MBI, 50 mg/kg bw/d amitrole, and 15/11.25 mg/kg bw/d cyanamide. Glands were thawed on ice and homogenized on a Tissue Lyser II (Qiagen, Hilden, Germany). Subsequently, total RNA isolation was performed using a RNeasy Micro Kit following manufacturer's protocol with the DNase step included (Qiagen). RNA quantity, quality, and purity was determined by a combination of a Nano-drop-1000 spectrophotometer (ND-1000, Fisher Scientific) and an Agilent 2100 Bioanalyzer System (Agilent Technologies, Santa Clara, CA, USA). Samples with RNA integrity (RIN)-score > 6.8 were used. RNA samples were diluted to 25 ng/µL, and stored at -80 °C.

## 2.5. Bulk RNA barcoding (BRB) library preparation and sequencing

Transcriptome analysis was done by 3Bulk RNA Barcoding and Sequencing (BRB-seq) [17] with library preparation as described by Giacosa et al. [18]. In short, reverse transcription and template switching was done on 4 µL total RNA of a concentration of 2.5 ng/µL,

cDNA was purified, and double-stranded cDNA was amplified by PCR. We used 50 ng ds cDNA for preparing sequencing libraries built by tagmentation using an Illumina Nextera XT kit (Illumina, #FC-131-1024, San Diego, CA, USA) following manufacturer's protocol. Subsequently, the library was sequenced by an Illumina Hiseq 4000 sequencing system following Illumina's instructions by Plateforme GenomEast (http://genomeast.igbmc.fr/, Illkirch-Graffenstaden, France). Finally, analysis of images and base calling was done using RTA 2.7.7 and bcl2fastq 2.17.1.14, and removal of adapted dimer reads by DimerRemover (https://sourceforge.net/projects/dimerremover/).

## 2.6. BRB-seq preprocessing and analysis

The preprocessing of raw data is described in detail in Giacosa et al. [18]. Briefly, the quality control pipeline requires a first read of 16 bases with a quality score > 10 with the first 6 bp being unique for the sample and the following 10 bp being a unique molecular identifier (UMI). The second reads were aligned with a rat reference transcriptome (UCSC website, release rn6, downloaded in August 2020). When reads were mapped to several positions in the genome these were omitted from further analysis. Finally, a gene count matrix was produced by counting the number of unique UMIs linked with each gene (lines) for each samples (columns) with the matrix undergoing subsequent normalization by rlog transformation package from DeSeq2 [19]. All data, both raw and preprocessed are available at the GEO repository under the accession number GSE242162 [20]. The transcriptomic results were also deposited in the TOXIsgN repository (https://toxsign.genouest.org/) [21].

## 2.7. Differential gene expression and functional analysis

The statistical filtration was done using R v4.0.3. The principal component analysis (PCA) was performed with the FactoMineR package implemented in R [22]. Differentially expressed genes were identified by comparing controls to exposed groups. Two filtrations were done for each comparison omitting all genes that were not expressed above a background cut-off, as well as genes that did not exceed 1.3-fold change compared to the control. A final statistical filtration was done using linear models for Microarray Data (LIMMA) package (Parkville, Australia) [23]. To detect statistically significantly differentially expressed genes (DEGs) a cut-off value of *p* = 0.05 adjusted with Benjamini & Hochberg method [24] was used.

Functional analysis was performed using the Annotation Mapping Expression and Network (AMEN) suite [25] with annotation terms considered enriched when the False Discovery Rate (FDR)-adjusted *p*-value was ≤0.05 and the number of DEGs within each term exceeded five. In total, terms from 27 databases were queried for functional analysis, including Gene Ontology (biological process, molecular function and cellular component), KEGG pathways, WikiPathways and GSEA gene sets.

# 3. Results

Male pups were exposed during gestation and after birth through lactation with dams being dosed by oral gavage from GD7 to PD16 to MBI, amitrole, or cyanamide at 10, 50, and 15/11.25 mg/kg bw/d. In the 16 days old male offspring serum TH concentrations were measured by UPLC-HRMS and thyroid glands were analyzed by BRB-seq (Fig. 1A).

## 3.1. Amitrole reduced serum T4

Previously, we measured T3 and T4 by RIA in serum from male and female pups at PD16 (n = 8, one male and one female sample from each of 4 litters) [5]. To assure that similar effects would be observed in males only and with the smaller sample size (n = 4), serum hormone levels of T4 and T3 were measured by UPLC-HRMS in PD16 male pups after
developmental exposure to MBI, amitrole, or cyanamide.

The absolute levels of T4 were approximately 10 times higher than T3 levels in control animals (Fig. 1B). Serum T3 concentrations were below the limit of quantification, but above the limit of detection. Amitrole caused a significant reduction of approximately 60 % in T4 compared to the control group. There were no statistically significant effects in MBI- and cyanamide-exposed animals. These findings were similar to the levels previously measured by RIA [5].

## 3.2. Amitrole caused marked and unique transcriptional changes

We found 183 DEGs in thyroid glands of male PD16 pups developmentally exposed to MBI, amitrole, or cyanamide (Fig. 2A and Supplementary material Table 2). Amitrole exposure led to the highest number of DEGs, with 149 gene transcripts identified of which about half were upregulated and the other half downregulated. In the MBI group, there were a total of 60 DEGs detected, with 28 upregulated and 32 downregulated DEGs. Cyanamide caused 27 DEGs of which 16 represented upregulated and 11 downregulated gene transcripts.

The DEGs clustered into nine distinct expression patterns (P1-P9) (Fig. 2B). The four largest patterns all contain genes primarily affected by amitrole, either alone or in concert with MBI: namely P2, P3, P5 and P6. The transcripts specific to amitrole-exposed animals belong to P3 (60 upregulated genes) and P6 (35 downregulated genes). P2 and P5 contain genes primarily upregulated and downregulated by MBI and amitrole, respectively. P4 and P8 contain genes not affected by amitrole, but differentially expressed in MBI and/or cyanamide groups: 7 upregulated genes in P4 and 8 downregulated genes in P8.

## 3.3. Amitrole affected genes associated with an active and remodeling thyroid gland

We performed functional analysis to associate DEGs to gene sets from 27 databases including the Gene Ontology or KEGG pathway databases (Table 1 and Supplementary material Table 3). The analysis of all DEGs for each compound showed changes in amitrole-exposed animals with a total of 30 enriched terms in the categories biological process, molecular function, and cellular component. The TH synthesis KEGG pathway was enriched by amitrole (7 genes). There were changes related to tissue structure and remodeling, including vasculature process in the circulatory system, cell cycle, and the extracellular space. Other terms indicated thyroid glands that responded to exogenous and endogenous substances, for instance cellular detoxification, cellular response to chemical stimulus, response to hormone, and response to endogenous stimulus. For MBI, 10 terms were enriched, including the GSEA pathways term “NRF2_pathway”. There were no significantly enriched terms for cyanamide.

## 3.4. MBI and amitrole had the largest intersection of DEGs

Of the 183 DEGs, 103 were specific to amitrole exposed animals, whereas 19 and 14 DEGs were distinct for MBI and cyanamide exposed animals, respectively (Fig. 3A). The greatest intersection was observed for the MBI and amitrole down-regulated and up-regulated genes, with 18 down-regulated DEGs and 15 up-regulated DEGs affected (Fig. 3A). The fold-changes of intersecting genes between amitrole and MBI were generally similar, as were the fold-changes of DEGs only affected by amitrole. For the remaining intersections the number of DEGs were four or below. A full list of intersecting genes between amitrole and MBI is shown in Fig. 3B.

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Work flow for BRB-seq data processing and heatmap of DEGs in thyroid glands of 16-day old male pups developmentally exposed to MBI, amitrole or cyanamide. A) Data analysis flow with filtrations steps based on limit of detection, fold change (FC) cut-off of 1.3 and statistics leading to 183 differentially expressed genes. n = 7-8 (one to two male pups/litter). B) Heat map of all DEGs after exposure to MBI, amitrole and cyanamide. Affected genes clustered in 9 distinct patterns (P1: 5 DEGs, P2: 18 DEGs, P3: 60 DEGs, P4: 7 DEGs, P5: 39 DEGs, P6: 25 DEGs, P7: 9 DEGs, P8: 8 DEGs, P9: 12 DEGs). The left part of the heatmap displays the statistically significant genes (DEGs) for each group. For instance, some genes in P5 are DEGs only in MBI exposed animals, whereas others are specific to amitrole, and some are similarly affected by both compounds. Heatmap legend: scaled log2-fold-change. Blue = significantly downregulated expression, red = significantly upregulated expression, white = no significant change. BRB: bulk-RNA-barcoding and sequencing, FC: fold change, MBI: 2-mercaptobenzimidazole, NDE: not differentially expressed, P: pattern. Parts of figure created with BioRender.com. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<table><caption>**Table 1**
Select enriched terms in male pups developmentally exposed to amitrole.</caption><thead><tr><th>Enriched term</th><th>Category</th><th>No. Genes</th><th>GO-term</th><th>P-values</th></tr></thead><tbody><tr><td>Response to hormone</td><td></td><td>30</td><td>0009725</td><td>0.01685</td></tr><tr><td>Response to endogenous stimulus</td><td></td><td>41</td><td>0009719</td><td>0.02074</td></tr><tr><td>Response to stress</td><td></td><td>62</td><td>0006950</td><td>0.01875</td></tr><tr><td>Cellular detoxification</td><td>Biological process</td><td>8</td><td>1990748</td><td>0.02720</td></tr><tr><td>Vascular process in circulatory system</td><td></td><td>11</td><td>0003018</td><td>0.04329</td></tr><tr><td>Cellular response to chemical stimulus</td><td></td><td>52</td><td>0070887</td><td>0.02074</td></tr><tr><td>Signaling receptor binding</td><td></td><td>31</td><td>0005102</td><td>0.04962</td></tr><tr><td>Protein binding</td><td>Molecular function</td><td>103</td><td>0005515</td><td>0.00001</td></tr><tr><td>Extracellular region</td><td></td><td>41</td><td>0005576</td><td>0.00063</td></tr><tr><td>Collagen trimer</td><td></td><td>5</td><td>0005581</td><td>0.04086</td></tr><tr><td>Extracellular matrix</td><td>Cellular component</td><td>12</td><td>0031012</td><td>0.04086</td></tr><tr><td>Collagen-containing extracellular matrix</td><td></td><td>9</td><td>0062023</td><td>0.04444</td></tr><tr><td>DUTERTRE Estradiol RESPONSE 24HR_UP</td><td></td><td>10</td><td>NA</td><td>0.02990</td></tr><tr><td>ZHOU CELL CYCLE GENES_IN_IR RESPONSE_24HR</td><td>GSEA:</td><td>7</td><td>NA</td><td>0.01703</td></tr><tr><td>WHITFIELD CELL CYCLE G2</td><td>ChemGenPerturbations</td><td>8</td><td>NA</td><td>0.01279</td></tr><tr><td>INTEGRIN3_PATHWAY</td><td></td><td>5</td><td>NA</td><td>0.03290</td></tr><tr><td>NRF2_PATHWAY</td><td>GSEA: PID</td><td>10</td><td>NA</td><td>0.00500</td></tr><tr><td></td><td>GSEA: WIKIPATHWAYS</td><td></td><td></td><td></td></tr></tbody></table>

Among genes typically associated with the TH synthesis system, only *Slc5a5* and *Tg* were affected significantly, and only by amitrole exposure. These effects clustered into pattern 3 (Fig. 2B), where amitrole generally led to an increased expression, while MBI and cyanamide caused few effects. No effects on classical TH synthesis genes were observed in the MBI- and cyanamide-exposed groups (Supplementary material Table 2).

## 4. Discussion

We found distinct transcriptional changes in thyroid glands of amitrole-exposed male rat pups, with 149 genes affected out of a total of 183 DEGs. In comparison, there were 60 and 27 DEGs in MBI- and cyanamide-exposed pups, respectively. Amitrole also caused a different phenotype than MBI and cyanamide, with a ~ 50 % increase in thyroid gland weight, 7-fold increase in serum TSH levels, and decrease in serum TH levels [5], accompanied by pronounced histological changes [6]. Hence, it is perhaps not surprising that amitrole also caused unique changes at the transcriptional level. The associated functional analysis pointed towards an active and remodeling thyroid gland: e.g., response to stress, response to toxic substance, extracellular matrix, and vasculature, which may be linked to the increased tissue weight and morphological changes in amitrole-exposed animals. Our data suggests that amitrole exposure induced phenotypic changes associated with TH system disruption and that these were accompanied by a markedly affected transcriptome compared to in MBI- and cyanamide-exposed animals.

Effects on the transcriptome have previously been examined with amitrole exposure in pubertal rats [14]. Pubertal rats were exposed for 7 days to high doses of amitrole, as well as three other TPO inhibitors and five indirect TH system disruptors. Though the total number of genes affected by amitrole in the present study far exceeds that of MBI and cyanamide, the total number of DEGs affected was much smaller than that reported in the Ohara et al. study [14], where approximately 8900 genes were affected, of which we identified approximately 50 [14]. Thus, very different results have been obtained in the two studies. We cannot determine the reasons for the discrepancies, but suspect that some may be due to differences in study designs. Very different technologies (microarray versus BRB-seq), life stages (pubertal versus developmental), doses (300 versus 50 mg/kg bw/d) were used, also noting potential differences in absorption, distribution, metabolism and excretion between direct exposure for 7 days and exposure of dams throughout gestation and lactation.

We had expected that the transcriptomic changes in amitrole-exposed pups would reveal genes that are typically associated with the TH synthesis apparatus. TSH levels were increased by approximately 7-fold [5] and TSH is known to control either directly or indirectly the expression of genes such as *Slc5a5* (NIS), *Tg*, *Tpo*, *Tshr*, *Slc26a4* (Pendrin), *Duox* [2,26], as well as the transcription factors *Pax8*, *Foxe1*, and *Nkx2.1* [26]. However, only *Tg* and *Slc5a5* (NIS) were upregulated after amitrole exposure. A total of 7 DEGs belonged to the “Thyroid hormone synthesis”-KEGG pathway (https://www.genome.jp/kegg-bin/shortest_pathway?rno04918), including *gs*, *grp94*, *asgpr*, *gsh-px*, and *atpase* in addition to *Slc5a5* (NIS) and *Tg*. Amitrole exposure increased *Tg* expression in FRTL-5 cells [27], whereas the expression was down-regulated in pubertally exposed rats [14]. The expression of *Slc5a5* (NIS) was not affected in these studies. In the Ohara et al. study [14], 15 DEGs belonged to the “Thyroid hormone synthesis”-KEGG pathway, but of these only *Tg* was also affected in our study. A less pronounced phenotype was observed in MBI-exposed animals, with no effects observed for genes related to TH synthesis. Taken together, our data suggest that TH synthesis genes may not be the most suitable indicators of disruption of the thyroid hormone system and HPT-axis activation in developing animals.

To identify potential global biomarkers of both mild and severe TH system disruption by developmental TPO inhibition, we turned to the MBI- and amitrole-exposed animals. Amitrole exposure caused obvious TH system disruption while the effects of MBI were less pronounced, and transient, in the pups [5]. Early in gestation the high dose of 20 mg/kg bw/d MBI decreased serum T3 and increased serum TSH levels in the pregnant dams. However, the 20 mg/kg bw/d dose group was terminated due to severe reproductive toxicity. Thus, only pups from the 10 mg/kg bw/d dose group were available for postnatal endpoints and thyroid gland transcriptomics. At this dose, there were no effects of MBI in the pregnant dams but there were increased T4 in the 6-day old pups. Taken together, our data indicate that developmental exposure to 10 mg/kg bw/d MBI cause a very mild TH system disruption. This is in line with previous reports from *in vivo* studies in adult animals [28–30]. Thus, both amitrole and MBI caused some degree of TH disruption and displayed the greatest overlap of DEGs. Of the intersecting genes only few were among the genes identified as a TPO inhibitor signature in pubertal rats [14]. Collectively, these findings indicate that either MBI exposure did not affect the TH system sufficiently to induce transcriptional changes characteristic of TPO inhibition or that the differences were due to exposure during different life stages, since the overlap for amitrole affected genes between this study and Ohara et al. [14] were also modest.

For the 33 intersecting genes observed between MBI- and amitrole-exposed animals, there are some interesting observations. Several genes involved in calcium signaling were affected. Of these, parathyroid hormone (*Pth*) and calcium sensing receptor (*Casr*) are particularly interesting. Both were downregulated with similar fold-changes in MBI-
<figure><img src="image_5.png" /></figure>

<table><caption>**Fig. 3.** Intersection plot and intersecting gene transcripts. A) Intersection plot for up- and down-regulated DEGs in cyanamide-, MBI-, and amitrole-exposed PD16 male pups. B) Intersecting genes, both up and down regulated, common between MBI and amitrole exposed animals. AMI: amitrole, DEGs: differentially expressed genes, FC: fold change, MBI: 2-mercaptobenzimidazole.</caption><thead><tr><th colspan="3">UP</th><th colspan="2">FC</th><th colspan="3">DOWN</th><th colspan="2">FC</th></tr><tr><th>Symbol</th><th>Full name</th><th>MBI</th><th>AMI</th><th>Symbol</th><th>Full name</th><th>MBI</th><th>AMI</th><th>Symbol</th><th>Full name</th><th>MBI</th><th>AMI</th></tr></thead><tbody><tr><td>Crbn</td><td>cereblon</td><td>1.41</td><td>1.32</td><td>Adipoq</td><td>adiponectin, C1Q and collagen domain containing</td><td>0.76</td><td>0.73</td><td>Aif1l</td><td>allograft inflammatory factor 1-like</td><td>0.64</td><td>0.66</td></tr><tr><td>Cryab</td><td>crystallin, alpha B</td><td>1.35</td><td>1.43</td><td>Casr</td><td>calcium-sensing receptor</td><td>0.75</td><td>0.73</td><td>Casr</td><td>calcium-sensing receptor</td><td>0.75</td><td>0.73</td></tr><tr><td>Fth1</td><td>ferritin heavy chain 1</td><td>1.67</td><td>1.74</td><td>Chga</td><td>chromogranin A</td><td>0.51</td><td>0.46</td><td>Col15a1</td><td>collagen type XV alpha 1 chain</td><td>0.73</td><td>0.64</td></tr><tr><td>Gstp1</td><td>glutathione S-transferase pi 1</td><td>1.38</td><td>1.91</td><td>Creld2</td><td>cysteine-rich with EGF-like domains 2</td><td>0.73</td><td>0.72</td><td>Col15a1</td><td>collagen type XV alpha 1 chain</td><td>0.73</td><td>0.64</td></tr><tr><td>Higd1b</td><td>HIG1 hypoxia inducible domain family, member 1B</td><td>1.38</td><td>2.08</td><td>Cryba2</td><td>crystallin, beta A2</td><td>0.65</td><td>0.56</td><td>Dgat2</td><td>diacylglycerol O-acyltransferase 2</td><td>0.68</td><td>0.69</td></tr><tr><td>Inmt</td><td>indolethylamine N-methyltransferase</td><td>1.32</td><td>1.64</td><td>Eno2</td><td>enolase 2</td><td>0.76</td><td>0.67</td><td>Eno2</td><td>enolase 2</td><td>0.76</td><td>0.67</td></tr><tr><td>Mfsd13a</td><td>major facilitator superfamily domain containing 13A</td><td>1.37</td><td>1.39</td><td>Herpud1</td><td>homocysteine inducible ER protein with ubiquitin like domain 1</td><td>0.71</td><td>0.75</td><td>Herpud1</td><td>homocysteine inducible ER protein with ubiquitin like domain 1</td><td>0.71</td><td>0.75</td></tr><tr><td>Muc20</td><td>mucin 20, cell surface associated</td><td>1.39</td><td>1.41</td><td>Mia</td><td>MIA SH3 domain containing</td><td>0.72</td><td>0.71</td><td>Mia</td><td>MIA SH3 domain containing</td><td>0.72</td><td>0.71</td></tr><tr><td>Osgin1</td><td>oxidative stress induced growth inhibitor 1</td><td>1.32</td><td>1.34</td><td>Mvk</td><td>mevalonate kinase</td><td>0.76</td><td>0.76</td><td>Mvk</td><td>mevalonate kinase</td><td>0.76</td><td>0.76</td></tr><tr><td>Pgd</td><td>phosphogluconate dehydrogenase</td><td>1.33</td><td>1.36</td><td>Nusap1</td><td>nucleolar and spindle associated protein 1</td><td>0.77</td><td>0.71</td><td>Nusap1</td><td>nucleolar and spindle associated protein 1</td><td>0.77</td><td>0.71</td></tr><tr><td>Ptgr1</td><td>prostaglandin reductase 1</td><td>1.42</td><td>1.58</td><td>Pth</td><td>parathyroid hormone</td><td>0.59</td><td>0.58</td><td>Pth</td><td>parathyroid hormone</td><td>0.59</td><td>0.58</td></tr><tr><td>Rhbg</td><td>Rh family B glycoprotein</td><td>1.37</td><td>1.53</td><td>Rhcg</td><td>Rh family, C glycoprotein</td><td>0.66</td><td>0.58</td><td>Rhcg</td><td>Rh family, C glycoprotein</td><td>0.66</td><td>0.58</td></tr><tr><td>Sat1</td><td>spermidine/spermine N1-acetyl transferase 1</td><td>1.35</td><td>1.36</td><td>Scgb1c1</td><td>secretoglobin, family 1C, member 1</td><td>0.61</td><td>0.75</td><td>Scgb1c1</td><td>secretoglobin, family 1C, member 1</td><td>0.61</td><td>0.75</td></tr><tr><td>Smim22</td><td>small integral membrane protein 22</td><td>1.33</td><td>1.37</td><td>Smoc2</td><td>SPARC related modular calcium binding 2</td><td>0.58</td><td>0.64</td><td>Smoc2</td><td>SPARC related modular calcium binding 2</td><td>0.58</td><td>0.64</td></tr><tr><td>Vnn1</td><td>vanin 1</td><td>1.42</td><td>1.53</td><td>Tspan6</td><td>tetraspanin 6</td><td>0.67</td><td>0.73</td><td>Tspan6</td><td>tetraspanin 6</td><td>0.67</td><td>0.73</td></tr></tbody></table>

and amitrole-exposed pups. In the parathyroid, CaSR sense small changes in serum calcium levels and induce parathyroid hormone secretion to maintain concentrations [31]. The calcium-signaling pathway was affected at the gene level in thyroid glands from rats exposed to acrylamide, though this did not specifically include *casr* and *pth* [32]. Whether the calcium signaling pathway is interesting as a biomarker for TPO inhibition should be investigated further, including consideration of whether the changes can be ascribed to changes in ratio of different cell populations of the parathyroid and thyroid glands. The patterns in intersecting genes may represent some of the first observable changes in thyroid glands of developing animals exposed to TPO inhibitors and may serve as sensitive and compound-independent

indicators of early TH system disruption and TPO inhibition. This hypothesis should be investigated further by using low doses of a panel of TPO inhibitors.

Phenotypical effects of cyanamide exposure were different from those observed for MBI and amitrole despite all compounds being *in vitro* TPO inhibitors. There were no changes in any thyroid-related parameter in cyanamide exposed animals [5]. Accordingly, there were few unique DEGs for cyanamide and little overlap with amitrole and MBI groups. Thus, our study did not identify any sign of TH system disruption by cyanamide. Noteworthy, there are indications that cyanamide can cause thyroid hormone system disruption *in vivo*; however, these findings were in studies using higher doses and/or longer duration of exposure in adult
animals [33,34].

We have identified a set of genes that may be early effect biomarkers of developmental TH system disruption in animals exposed to *in vitro* TPO inhibitors. Upon sufficient validation this signature could be useful for molecular identification of TPO inhibiting properties in toxicity studies. Further studies are necessary, including clarification of whether the major effects are dictated by changes in tissue development and organization, if the transcriptional signature is substance independent, and what the applicability domain is. Furthermore, it is of interest whether this signature can also be used as an *in vitro* indicator of TH disruption and TPO inhibition, for instance in advanced 3D thyroid gland models or in simpler cell line assays.

### CRediT authorship contribution statement

**Anna Kjerstine Rosenmai:** Writing – original draft, Data curation.
**Terje Svingen:** Writing – review & editing, Funding acquisition.
**Bertrand Evrard:** Writing – review & editing, Formal analysis.
**Khanh Hoang Nguyen:** Writing – review & editing, Methodology, Formal analysis.
**Camilla Nielsen:** Writing – review & editing, Formal analysis.
**Marta Axelstad:** Writing – review & editing, Resources.
**Frédéric Chalmel:** Writing – review & editing, Visualization, Methodology, Data curation.
**Louise Ramhøj:** Writing – review & editing, Project administration, Conceptualization.

### 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.

### Data availability

BRB-seq data is available at the GEO repository under the accession number GSE242162. Other raw data presented in this manuscript will be made available upon request.

### Acknowledgements

This work was supported by the EU Horizon 2020 project ATHENA: Assays for the identification of Thyroid Hormone axis-disrupting chemicals: Elaborating Novel assessment strategies (grant number 825161; Kortenkamp et al. 2020 [35]) and the Danish Environmental Protection Agency as a project under The Danish Centre for Endocrine Disrupters (CeHoS).

TH analysis was performed through accessing research infrastructure at The National Food Institute, Technical University of Denmark, including FOODHAY (Food and Health Open Innovation Laboratory, Danish Roadmap for Research Infrastructure). Some of the figures were generated using Biorender. We would like to thank Mette Voigt Jessen, Heidi Broksø Letting, Lillian Sztuk and Dorte Lykkegaard Korsbech for technical assistance. We thank GenomEast for performing the sequencing as well as the BioFacility at the Technical University of Denmark for animal care and husbandry.

### Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.ygeno.2024.110938.

### References

[1] D.P. Carvalho, C. Dupuy, Thyroid hormone biosynthesis and release, Mol. Cell. Endocrinol. 458 (2017) 6–15.

[2] T.M. Ortiga-Carvalho, M.I. Chiamolera, C.C. Pazos-Moura, F.E. Wondisford, Hypothalamus-pituitary-thyroid axis, Compr. Physiol. 6 (2016) 1387–1428.

[3] R.T. Zoeller, S.W. Tan, R.W. Tyl, General background on the hypothalamic-pituitary-thyroid (HPT) axis, Crit. Rev. Toxicol. 37 (2007) 11–53.

[4] K.P. Friedman, et al., Tiered high-throughput screening approach to identify thyroperoxidase inhibitors within the toxcast phase I and II chemical libraries, Toxicol. Sci. 151 (2016) 160–180.

[5] L. Ramhøj, et al., Testing for heterotopia formation in rats after developmental exposure to selected in vitro inhibitors of thyroperoxidase, Environ. Pollut. 283 (2021).

[6] L. Ramhøj, et al., Perinatal exposure to the thyroperoxidase inhibitors methimazole and amitrole perturbs thyroid hormone system signaling and alters motor activity in rat offspring, Toxicol. Lett. 354 (2022) 44–55.

[7] A. Hood, Y.P. Liu, V.H. Gattone, C.D. Klaassen, Sensitivity of thyroid gland growth to thyroid stimulating hormone (TSH) in rats treated with Antithyroid drugs, Toxicol. Sci. 49 (1999) 263–271.

[8] I. Hassan, et al., Neurodevelopment and thyroid hormone synthesis inhibition in the rat: quantitative understanding within the adverse outcome pathway framework, Toxicol. Sci. 160 (2017) 57–73.

[9] I. Hassan, et al., Extrapolating in vitro screening assay data for Thyroperoxidase inhibition to predict serum thyroid hormones in the rat, Toxicol. Sci. 173 (2020) 280–292.

[10] M. Axelstad, et al., Developmental neurotoxicity of Propylthiouracil (PTU) in rats: relationship between transient hypothyroxinemia during development and long-lasting behavioural and functional changes, Toxicol. Appl. Pharmacol. 232 (2008) 1–13.

[11] S. Spaulding, Biological Actions of Thyrotropin. In the Thyroid: A Fundamental and Clinical Text 227–233, Lippincott Williams and Wilkins, Philadelphia, 2000.

[12] L.P. Fernández, A. López-Márquez, P. Santisteban, Thyroid transcription factors in development, differentiation and disease, Nat. Rev. Endocrinol. 11 (2015) 29–42.

[13] A. López-Márquez, C. Carrasco-López, C. Fernández-Méndez, P. Santisteban, Unraveling the complex interplay between transcription factors and signaling molecules in thyroid differentiation and function, From Embryos to Adults, Front. Endocrinol. 12 (2021).

[14] A. Ohara, F. Yamada, T. Fukuda, N. Suzuki, K. Sumida, Specific alteration of gene expression profile in rats by treatment with thyroid toxicants that inhibit thyroid hormone synthesis, J. Appl. Toxicol. 38 (2018) 1529–1537.

[15] M.E. Gilbert, et al., Gestational exposure to perchlorate in the rat: thyroid hormones in fetal thyroid gland, serum, and brain, Toxicol. Sci. 188 (2022) 117–130.

[16] J. Dong, et al., Effects of perinatal di (2-ethylhexyl) phthalate exposure on thyroid function in rat offspring, Environ. Toxicol. Pharmacol. 67 (2019) 53–60.

[17] D. Alpern, et al., BRB-seq: ultra-affordable high-throughput transcriptomics enabled by bulk RNA barcoding and sequencing, Genome Biol. 20 (2019) 1–15.

[18] S. Giacosa, et al., Cooperative blockade of CK2 and ATM kinases drives apoptosis in VHL-deficient renal carcinoma cells through ROS overproduction, Cancers (Basel) 13 (2021) 1–27.

[19] M.I. Love, W. Huber, S. Anders, Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2, Genome Biol. 15 (2014) 1–21.

[20] R. Edgar, M. Domrachev, A.E. Lash, Gene expression omnibus: NCBI gene expression and hybridization array data repository, Nucleic Acids Res. 30 (2002) 207–210.

[21] T.A. Darde, et al., TOXISgN: a cross-species repository for toxicogenomic signatures, Bioinformatics 34 (2018) 2116–2122.

[22] S. Lè, J. Josse, A. Rennes, F. Husson, FactoMineR: an R package for multivariate analysis, J. Stat. Softw. 25 (2008) 1–18.

[23] G.K. Smyth, E. Hall, Linear models and empirical Bayes methods for assessing differential expression in microarray experiments linear models and empirical Bayes methods for assessing differential expression in microarray experiments, Stat. Appl. Genet. Mol. Biol. 3 (2004) 1–25.

[24] Y. Benjamini, Y. Hochberg, Controlling the false discovery rate: a practical and powerful approach to multiple, J. R. Stat. Soc. Ser. B Methodol. 57 (1995) 289–300.

[25] F. Chalmel, M. Primig, The annotation, mapping, expression and network (AMEN) suite of tools for molecular systems biology, BMC Bioinformat. 9 (2008).

[26] A. López-Márquez, C. Carrasco-López, C. Fernández-Méndez, P. Santisteban, Unraveling the complex interplay between transcription factors and signaling molecules in thyroid differentiation and function, from embryos to adults, Front. Endocrinol. 12 (2021) 654569, https://doi.org/10.3389/fendo.2021.654569.

[27] Pan Hongmei, Sun Yulan, Zhang Lishi, The effects of amitrole on thyroglobulin and iodide uptake in FRTL-5 cells, Toxicol. Ind. Health 27 (2011) 187–192.

[28] ECHA, Benzimidazole-2-Thiol, Registration Dossier, 2021 https://echa.europa.eu/da/registration-dossier/-/registered-dossier/6188/1.

[29] T. Yamano, T. Noda, M. Shimizu, S. Morita, The adverse effect of oral MBI expo in pregnant rats and their fetuses, Fundam. Appl. Toxicol. 25 (1995) 218–223.

[30] Y. Kawasaki, et al., Toxicity study of a rubber antioxidant, 2-mercaptobenzimidazole, by repeated oral administration to rats, J. Toxicol. Sci. 23 (1998) 53–68.

[31] G. Díaz-Soto, A. Rocher, C. García-Rodríguez, L. Núnez, C. Villalobos, The calcium-sensing receptor in health and disease, in: International Review of Cell and Molecular Biology 327, Elsevier Inc., 2016, pp. 321–369.

[32] N.L. Chepelev, et al., Transcriptional profiling of male F344 rats suggests the involvement of calcium signaling in the mode of action of acrylamide-induced thyroid cancer, Food Chem. Toxicol. 107 (2017) 186–200.
[33] ECHA, Cyanamide, Registration Dossier, 2023 https://echa.europa.eu/da/ registration-dossier/-/registered-dossier/15823/7/9/1.

[34] National Cancer Institute, Bioassay of calcium Cyanamide for possible carcinogenicity, in: Technical Report Series No. 163, 1979.

[35] A. Kortenkamp, et al., Removing critical gaps in chemical test methods by developing new assays for the identification of thyroid hormone system-disrupting chemicals—the athena project, Int. J. Mol. Sci. 21 (2020) 1–13.