Long-term transcriptomic and proteomic effects in Sprague Dawley rat thyroid and plasma after internal low dose 131I exposure

Background Radioiodide (131I) is commonly used to treat thyroid cancer and hyperthyroidis.131I released during nuclear accidents, have resulted in increased incidence of thyroid cancer in children. Therefore, a better understanding of underlying cellular mechanisms behind 131I exposure is of great clinical and radiation protection interest. The aim of this work was to study the long-term dose-related effects of 131I exposure in thyroid tissue and plasma in young rats and identify potential biomarkers. Materials and methods Male Sprague Dawley rats (5-week-old) were i.v. injected with 0.5, 5.0, 50 or 500 kBq 131I (Dthyroid ca 1–1000 mGy), and killed after nine months at which time the thyroid and blood samples were collected. Gene expression microarray analysis (thyroid samples) and LC-MS/MS analysis (thyroid and plasma samples) were performed to assess differential gene and protein expression profiles in treated and corresponding untreated control samples. Bioinformatics analyses were performed using the DAVID functional annotation tool and Ingenuity Pathway Analysis (IPA). The gene expression microarray data and LC-MS/MS data were validated using qRT-PCR and ELISA, respectively. Results Nine 131I exposure-related candidate biomarkers (transcripts: Afp and RT1-Bb, and proteins: ARF3, DLD, IKBKB, NONO, RAB6A, RPN2, and SLC25A5) were identified in thyroid tissue. Two dose-related protein candidate biomarkers were identified in thyroid (APRT and LDHA) and two in plasma (DSG4 and TGM3). Candidate biomarkers for thyroid function included the ACADL and SORBS2 (all activities), TPO and TG proteins (low activities). 131I exposure was shown to have a profound effect on metabolism, immune system, apoptosis and cell death. Furthermore, several signalling pathways essential for normal cellular function (actin cytoskeleton signalling, HGF signalling, NRF2-mediated oxidative stress, integrin signalling, calcium signalling) were also significantly regulated. Conclusion Exposure-related and dose-related effects on gene and protein expression generated few expression patterns useful as biomarkers for thyroid function and cancer.


Introduction
The thyroid is an endocrine gland that produces iodine-containing hormones involved in metabolism and brain development. There are several isotopes of iodine, including radioactive 131 I and 123 I. Due to the natural accumulation of iodine in the thyroid, 131 I (half-life 8.04 days [1]) has frequently been used to treat thyroid malignancies, including e.g. hyperthyroidism and thyroid cancer and it may be used in thyroid scintigraphic imaging. Also, 131 I labelled radiopharmaceuticals, e.g. 131 I-MIBG, used for treatment of other diseases contribute to irradiation of healthy thyroid tissue. In the clinic the absorbed dose to thyroid depends on the type of radiopharmaceutical, thyroid function, administered activity and age. According to the Swedish Radiation Safety Authority (2019), the absorbed dose to thyroid in adults was 0.063-7.2 Gy for diagnostics, and 60-2600 Gy for therapy [2,3]. 131 I-MIBG therapy resulted in around 2 Gy absorbed dose to thyroid in young children. Furthermore, 131 I is one of the most commonly released radionuclides in nuclear weapon detonations and nuclear accidents. For medical use 131 I is usually produced by (n, γ) reaction of 130 Te to 131 Te (half-life 25 min), which decays to 131 I. Another way is by separation of fission products from 235 U [4].
After the Chernobyl nuclear accident, individuals living in surrounding areas in Ukraine, Belarus, and Russia were initially exposed to 131 I, and the absorbed dose to thyroid in children was 0-10 Gy, with the majority in the 0-5 Gy range and 1 Gy for pre-school children [5]. This resulted in significantly increased rates of primarily papillary thyroid cancer (PTC) in exposed children and adolescents [6][7][8]. Compared to adults, the biological uptake (concentration) of radioiodine is higher in young children due to the higher proliferation rate of the thyrocytes and the smaller thyroid size, which together with higher milk consumption by children may explain the absorbed dose to the thyroid in children after the Chernobyl accident. However, the increased absorbed dose to thyroid in children cannot explain the difference in thyroid cancer incidence when comparing children and adults [9]. Therefore, we need better knowledge of the effects of ionizing radiation on thyroid tissue, especially in young individuals. It is important to study both the short-and long-term radiobiological effects on thyroid tissue, since tumour initiation occurs early after exposure, while cancer diagnosis occurs years after exposure [10]. Exposure to ionizing radiation results in the activation/inhibition of several different cellular processes, where a number of the activated genes may have altered expression patterns [11,12].
We have previously investigated changes in mRNA expression at early time points after 131 I exposure in thyroid tissue samples from adult mice and rats using microarray techniques [13][14][15][16]. Several potential biomarkers were identified in thyroid tissue, specific for the administered activity level of 131 I and the time after injection, e.g. Klk1, the Klk1b family, Agpat9, Plau, Prf1, and S100a8. Various biological functions were affected, and we showed in adult mice (6-8 mo.) effects related to cell adhesion 24 h after low-medium dose exposure to 131 I [13]. In addition, calcium related gene changes and down regulation of parathyroid hormone and PPARG were found 24 h after 131 I exposure of rat thyroid [14].
Our previous studies were performed on adult animals. It is thus necessary to perform radiobiological studies on young animals to evaluate age-related radiobiological effects. Taken together, transcriptome and proteome profiling may potentially be useful for biomarker discovery, as these techniques detect variations in gene expression over time for different stimuli, e.g. radiation.
The aim of the present work was to investigate biological effects associated with low to medium dose 131 I exposure in young rats by studying changes in gene and protein expression in thyroid tissue and plasma, with special focus on exposure-related and dose-related effects and to identify potential biomarkers related to thyroid function alteration and cancer. Although a large number of regulated genes and proteins were identified, few demonstrated expression patterns (uniformly increased or decreased expression levels for all doses or expression levels increasing or decreasing with dose) useful as biomarkers for thyroid function and cancer.

Study design
Twenty young (5-week-old) male Sprague Dawley rats (Charles River Laboratories International, Inc., Salzfeld, Germany) were randomly divided into five groups containing four individuals each. Four groups were i.v. injected with 0.5, 5, 50 or 500 kBq 131 I (GE Healthcare, Braunschweig, Germany) to obtain an absorbed dose to the thyroid of D thyroid = 1, 10, 100 or 1000 mGy, respectively. One group was not treated and used as control. The animals were under daily supervision and received food (standard rat chow) and water ad libitum. The method for estimation of absorbed dose to the thyroid after i.v. injection of 131 I has been described elsewhere [12]. In brief, the absorbed dose to thyroid was calculated using the Medical Internal Radiation Dose formalism [17] where D thyroid is the absorbed dose to the thyroid, Â is the cumulated 131 I activity in thyroid, E is the mean emitted energy per disintegration, F is the absorbed fraction, m thyroid is the thyroid mass, A inj is the injected activity, a thyroid is the activity in organs at time t, t D is the time to which the absorbed dose contribution is calculated and λ is the radioactivity decay constant. Nine months after 131 I injection, the animals were killed by cardiac puncture under anaesthesia with pentobarbitalnatrium (APL; Kungens kurva, Sweden), followed by thyroid excision and collection of blood samples using heparin syringes. Thyroid samples were flash-frozen in liquid nitrogen and stored at -80oC until further processing. Plasma was collected from individual blood samples via centrifugation and stored at -80oC. The design of this study was approved by the Ethical Committee on Animal Experiments in Gothenburg, Sweden (Permit Number: 146-2015).

Gene expression microarray analysis
Total RNA samples were extracted from flash-frozen thyroid tissue samples with the RNeasy Lipid Tissue Mini Kit (Qiagen; Hilden, Germany). Microarray analysis was performed at the Bioinformatics and Expression Analysis Core Facility at Karolinska Institute (Stockholm, Sweden) using Agilent Sureprint G3 Rat GE 8x60K microarrays (Agilent, Santa Clara, CA, USA). Nexus Expression 3.0 (BioDiscovery; El Segundo, CA, USA) was used to identify differentially expressed transcripts by comparing thyroid tissue from exposed and non-exposed rats. The fold change and FDR-adjusted p-value cut-offs (Benjamini-Hochberg method) were set to > 1.5, <-1.5 and < 0.01, respectively. Functional annotation data based on Gene Ontology (GO) terms with a p-value < 0.05 were compiled with Nexus Expression and associated with higher level cellular functions using an in-house model based on parental GO terms presented in Fig 5, including eight main categories and 31 subcategories (www.geneonology.org) [18]. Thyroidspecific genes were identified using the Human Protein Atlas database (www.proteinatlas.org) [19,20]. The gene expression data discussed in this publication have been deposited in NCBI's Gene Expression Omnibus [21] and are accessible through GEO Series accession number GSE146051 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE146051).

Liquid chromatography-tandem mass spectrometry (LC-MS/MS) analysis
Protein extraction and TMT (tandem mass tags) protein labelling. The thyroid tissue samples were homogenised in lysis buffer (25 mM Triethylammoinium bicarbonate, TEAB, 2% Sodium dodecyl sulfate (SDS)) using a FastPrep1-24 instrument (MP Biomedicals, OH, USA). Total protein concentration was determined with Pierce™ BCA Protein Assay (Thermo Fisher Scientific, Waltham, MA, USA). The samples were then group wise combined in equal amounts from each animal in the group. Tissue samples, 50 μg of total protein for each pooled group and 50 μg of a reference pool containing equal amounts of all samples, were reduced by DL-Dithiothreitol (DTT) and then trypsin (Pierce Trypsin Protease, MS Grade, Thermo Fisher Scientific, Waltham, MA, USA) digested using the filter-aided sample preparation (FASP) method modified according to Wisniewski et al [22].
Rat plasma samples were similarly combined per group and a total of 30 μl per pooled group was depleted of high-abundance proteins using the Seppro Rat spin column kit (Sigma). In short, samples diluted with 8 M urea were applied on Nanosep 30k Omega filters (Pall Life Sciences, Port Washington, NY, USA) and 8 M urea was used to repeatedly wash away the SDS. Alkylation was performed with methyl methane thiosulfonate (MMTS) diluted in digestion buffer (1% sodium deoxycholate (SDC), 20 mM TEAB) and the filters were repeatedly washed with digestion buffer. Trypsin in a ratio of 1:100 relative to protein amount was added with 25 mM TEAB and the samples were incubated in 37˚C overnight. Another portion of trypsin was added the following morning and the samples were incubated for 4 h at 37˚C. Peptides were subjected to isobaric mass tagging reagent TMT1 according to the manufacturer's instructions (Thermo Fisher Scientific, Waltham, MA, USA). In a set, each sample and a reference was labelled with a unique tag from a TMT 6plex or 10plex isobaric mass tag labelling kit. Following TMT labelling, the samples in a set were pooled and concentrated and acidified to about pH 2 to precipitate SDC.
The combined TMT samples were fractionated by basic-pH reversed-phase chromatography (bRP-LC) on an Ä KTApurifier (Amersham Pharmacia Biotech AB, Uppsala, Sweden) using a reversed-phase XBridge BEH C18 column (3.5 μm, 3.0x150 mm, Waters Corporation), 10 mM ammonium formate buffer at pH 10.0 as solvent A and 90% acetonitrile, 10% 10 mM ammonium formate at pH 10.0 as solvent B and the flow of 400 μl/min. Samples were separated into 27 fractions on a gradient elution from 0% to 7% B for 1 min, from 7% to 50% B for 32 min followed by an increase to 100% B for 1.5 min, and holding 100% B for 2 min. The fractions were evaporated and reconstituted in 15 μl of 3% acetonitrile, 0.2% formic acid for nLC-MS analysis.
Liquid chromatography-mass spectrometry analysis. The fractionated samples were analysed on an Orbitrap Fusion Tribrid mass spectrometer (Thermo Fisher Scientific) interfaced with an Easy-nLC 1000 liquid chromatography system (Thermo Fisher Scientific); solvent A was 0.2% formic acid (FA) in water and solvent B was 0.2% FA in acetonitrile. Peptides were trapped on an in-house packed 3 cm pre-column, and separated on an analytical column (75 μm x 30 cm, both backed with Reprosil-Pur C18 material, particle size 3 μm, Dr. Maisch) using the linear gradient from 5% to 25% B for 45 min, from 25% to 80% B for 5 min, 80% B for 10 min, at 300 nL/min flow rate.
Precursor MS scans were performed at 120 000 resolution, with the 4e5 AGC target in the m/z range 380-1200 with wide quadrupole isolation; the most abundant precursors with charges 2 to 7 were selected for fragmentation over the 3 s cycle time with the dynamic exclusion duration of 30 s. Precursors were isolated with 1.6 window, fragmented by collision induced dissociation (CID) at 30% collision energy with a maximum injection time of 70 ms and AGC target 1e4, and the MS2 spectra were detected in the ion trap followed by the synchronous isolation of the 5 most abundant MS2 fragment ions within m/z range of 400-900, and fragmentation by higher-energy collision dissociation (HCD) at the 55% collision energy; the resulting MS3 spectra were detected in the Orbitrap at 60,000 resolution in the m/z range 100-500 with the maximum injection time 120 ms and AGC target 1e5.

Proteomic data analysis
Peptide and protein identification and quantification were performed using Proteome Discoverer version 1.4 (Thermo Fisher Scientific) with Mascot 2.3 or 2.5.1 (Matrix Science, London, United Kingdom) against the SwissProt database with taxonomy Rattus norvegicus, version 2015/04 (7935 sequences). Trypsin with no missed cleavage was used as a cleavage rule, precursor tolerance was set to 5 ppm and MS2 fragment tolerance was set to 500 mmu, monooxidation on methionine was set as a variable modification, cysteine methylthiolation, TMT-6 label on lysine and peptide N-termini were set as fixed modifications. Percolator was used for the peptide-spectrum match (PSM) validation with the strict false discovery rate (FDR) threshold of 1%.
The TMT reporter ions were identified in the MS3 HCD spectra with a mass tolerance of 3 mmu, and the resulting reporter abundance values for each sample were normalized on protein median in Proteome Discoverer 1.4. In total 3700 and 9747 peptides were found in plasma and thyroid tissue, respectively. Each protein was identified using 1-151 and 1-102 unique peptides for plasma and thyroid tissue, respectively.
The fold change cut-off were set to >1.5 and <-1.5, respectively. GO terms and the DAVID functional annotation tool (https://david.ncifcrf.gov/) were used for functional annotation of the LC-MS/MS data with a p-value cut off of <0.05, using modified Fisher´s exact test (EASE scores) and were categorised in the same manner as the transcripts, using the in-house model [23,24]. Thyroid-specific proteins were identified using the Human Protein Atlas database (https://www.proteinatlas.org) [19,20]. The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE [25] partner repository with the dataset identifier PXD017715.

Ingenuity pathway analysis (IPA)
The Ingenuity Pathway Analysis software (IPA; Ingenuity Systems, Redwood City, USA) was used to analyse affected canonical pathways, biological functions and upstream regulators for both transcripts and proteins. The overlap between the experimental data and the Ingenuity Pathways Knowledge Base was determined using Fisher´s exact test (p<0.05). The z-score was used to determine the activation state of the upstream regulators; z>2 indicates activation, while z<-2 indicates inhibition.

Validation of the microarray and the LC-MS/MS results
The Slc25a5 gene and DSG4 protein were selected for validation of expression microarray and LC-MS/MS results, respectively. In order to validate the microarray results, qRT-PCR was performed on thyroid extracted mRNAs. In brief, the mRNAs were converted to cDNA using SuperScript™ VILO™ cDNA Synthesis Kit (Invitrogen, 11754-050), according to the manufacturer's instructions. Primers corresponding to the Slc25a5, gene were run (Applied Biosystems) on a real-time 7500 HT sequence detection system (Applied Biosystems), with the SYBR green detection system (Applied Biosystems) The beta-actin, Gapdh and HPTR1 genes were used as controls based on their uniform expression in thyroid tissue (Applied Biosystems). The relative expression of the analysed genes was determined in relation to the expression of the control genes using the ΔΔC t method [14]. The intensity value was plotted against the ΔC t value for each animal.
To validate the results from LC-MS/MS analysis, enzyme-linked immunosorbent assay (ELISA) was used for the DSG4 protein (MyBioSource Cat no; MBS9324066). In brief, standards and samples were prepared according to the manufacturer's instructions, and were added into assigned wells. HRP conjugate was added to each well. The plate was incubated at 37 oC for one hour, and then washed according to the instructions. Chromogen solutions (A and B) was added to each well, followed by stop solution. Absorbance was measured for each well at 450 nm. The fold change value for each group was plotted against mean absorbed dose.

Histological analysis
Thyroid tissue samples from each individual rat were fixed in formaldehyde and then imbedded in paraffin. Paraffin slices (4 μm) were stained using haematoxylin-eosin to evaluate thyroid tissue morphology, abnormal tissue structure, and the presence of tumour cells. The slides were evaluated by two certified pathologists.

Regulated transcripts/proteins in thyroid tissue and plasma
In thyroid tissue, the transcriptomic analysis revealed 380 differentially regulated transcripts (corresponding to 347 unique genes) associated with 131 I exposure ( Fig 1A). The majority of the identified transcripts (n = 251) were only found at one activity level (S1A Table), of which eight transcripts (seven genes), 130 transcripts (125 genes), 18 transcripts (17 genes), and 95 transcripts (92 genes) were identified for the 0.5, 5, 50 and 500 kBq group, respectively (S1 Table). Only two transcripts (Afp and Rt1-Bb) were regulated at all activity levels (Fig 2A).
The proteomic analysis identified 792 and 273 differentially expressed proteins associated with 131 I exposure in thyroid tissue ( Fig 1B) and plasma (Fig 1C), respectively. In the thyroid, 247 proteins were regulated at only one activity level, including 54, 112, 47 and 34 proteins identified in the 0.5, 5, 50 and 500 kBq groups, respectively (S1B Table). Furthermore, 133 proteins were regulated at all activity levels (Table 1). Of the 273 proteins differentially regulated in plasma, 164 were only found at one activity level (S1C Table) (47, 19, 50 and 48 proteins were regulated in the 0.5, 5, 50 and 500 kBq groups, respectively) and two proteins (TMSB4X and SPP1) were regulated at all activity levels ( Table 1).

Commonly detected proteins in thyroid tissue and plasma
Interestingly, 63 proteins were commonly regulated in both thyroid and plasma for the same group ( Table 2). The majority of the differentially regulated proteins were found in the 500 kBq and 0.5 kBq groups and were, in general, either over-or underexpressed in both thyroid and plasma samples, irrespective of administered activity. A number of proteins were identified with similar expression patterns in both thyroid and plasma, including 20 proteins in the 0.

Exposure-related candidate biomarkers
Potential biomarkers for 131 I exposure, irrespective of injected 131 I activity, were selected from the transcriptomic and proteomic data as having a uniform increase or decrease in expression levels in all four groups. In thyroid tissue, the Afp and RT1-Bb transcripts were commonly regulated in all four groups, where Afp expression was consistently down-regulated and RT1-Bb up-regulated, but with no clear activity-related response (Fig 2A).
The majority of the 133 differentially regulated proteins identified in thyroid tissue had a uniform increase or decrease in protein expression level in the four groups, and are thus candidates for exposure related biomarkers useful to determinate if an individual is exposed or not, Exposure-related transcripts/proteins identified as candidate biomarkers in thyroid tissue, with their related GO terms. a) Transcript candidate biomarkers associated with response 131 I exposure for all groups. Two transcripts were found in thyroid tissue. b) Seven protein candidate biomarkers for radiation exposure for all groups in thyroid tissue. Associated GO terms for each transcript and protein are displayed.
https://doi.org/10.1371/journal.pone.0244098.g002 Table 1. Proteins in thyroid tissue with differential expression levels above or below 1.5 for all activity groups (0.5, 5, 50 and 500 kBq) (totally 133 proteins). The 20 most abundant proteins for each activity group are marked in grey. Eight of these proteins were present in all four groups (ARF3, DLD, IKBKB, KTR4, NONO, RAB6A, RPN2 and SLC25A5). For the majority of proteins, the expression levels were ether increased or decreased for all activity levels.  Table 2. Expression of statistically significant regulated proteins in both thyroid and plasma. The figure displays all 63 proteins that were statistically significantly regulated in both thyroid tissue and plasma for the same activity level (0.5, 5, 50 or 500 kBq) nine months after 131 I injection. The fold change cut-off was set to ±1.5. A positive or negative fold change is represented by + and-, respectively. If there was a difference between direction of regulation in thyroid and plasma, the thyroid data is presented first and then the plasma data, separated by /. Associated GO terms are displayed for each protein.  (Table 1). Seven proteins found in thyroid tissue (ARF3, DLD, IKBKB, NONO, RAB6A, RPN2, and SLC25A5) had elevated expression level with fold change > 2.0 ( Fig  2B). No protein with homogenous expression levels was obtained in all four groups for plasma. For the three groups with the lowest activities (0.5-50 kBq), six proteins found in thyroid tissue (AGMN, APOD, CAPZA2, GSTZ1, FAM213A, RBM20) and five proteins in plasma (CYP27B1, DSG4, RRAS, TAOK3, TGM3) were identified as potential exposure biomarkers (S2 Table). The corresponding candidate biomarkers for the intermediate activity groups (5- Table). The corresponding candidate biomarkers for the groups with the highest activities (50-500 kBq) included one thyroid transcript (Cd300lg), eleven protein found in thyroid tissue (BLOC1S2, CALCOCO1, COMP, Ester hydrolase C11orf54 homolog, GMFB, HBB, IDH2, LDHB, NDUFS2, PSMC2, YBX1) and eleven proteins in plasma (ATRN, CHGA, CPQ, CTBS, CYCS, ESYT1, F7, LTA4H, PF4, PVALB, SMPX) (S3 Table).

Dose-related candidate biomarkers
To identify dose-related biomarkers, we analysed common transcripts or proteins with a monotonic dose-response relationship (either increasing or decreasing expression levels for all groups) in all four or in three of the groups (0.5-50 kBq and 5-500 kBq). No differentially regulated transcripts with dose-response features were found in the transcriptome analysis of rat thyroid tissue. In contrast, the LDHA and APRT proteins were significantly overexpressed for all four groups in a dose-dependent manner in the proteomic analysis of thyroid tissue (Fig 3). Furthermore, seven dose-related proteins were identified in three of the groups (0.5-50 kBq or 5-500 kBq). For the three highest activities (5-500 kBq), two proteins with underexpression (NOL3 and YBX3) and two proteins with overexpression, (SCL3A2 and DSTN) were identified. For the three lowest activities (0.5-50 kBq) AGRN was underexpressed, while two proteins were overexpressed (MYBPC3 and HADHA). In plasma, the TGM3 and DSG4 proteins were overexpressed in a dose-dependent manner for the three lowest activities (0.5-50 kBq).

Validation of methods and histological evaluation of thyroid samples
The validation of microarray data using RT-qPCR on extracted thyroid mRNA for the Slc25a5 gene revealed a correlation (R2 = 0.52). (Fig 4). The validation of the LC-MS/MS analysis using ELISA of the DSG4 protein demonstrated a significant correlation (R2 = 0.71) (Fig 4). Three of the sixteen thyroid tissue samples from exposed animals showed in some parts abnormal tissue structures (neoplastic tissue).

Gene ontology annotations
Enriched biological processes for thyroid transcripts (BioDiscovery Nexus Expression 3.0 analysis) and proteins detected in thyroid and plasma tissues (DAVID functional annotation tool) were categorised into GO terms associated with cellular function and the extent of regulation calculated for each GO term (Fig 5). A clear association with specific GO terms was obtained for the thyroid transcripts and the highest number of associated terms was seen for proteins found in thyroid tissue. Some of the groups differed from the rest, e.g. proteins found in thyroid tissue were predominantly associated with apoptosis, suggesting cellular damage. However, the main cellular functions were similarly affected for all tissues, with metabolism (nucleic acid-related response) displaying the most pronounced response to 131 I irradiation. Less pronounced response to irradiation was shown for DNA integrity.
Gene Ontology association analysis was also performed for transcripts and proteins identified as candidate biomarkers. The analysis showed that the Afp and RT1-Bb exposure-related Table 3. Genes/proteins associated with thyroid function and thyroid cancer for all groups in thyroid tissue and plasma in young rats injected with 0.5-500 kBq 131 I. The fold change for each transcript/protein is presented in parenthesis.

Fig 4. Validation of microarray and LC-MS/MS methods using the Slc25a5 gene (left) and DSG4 protein (right).
The intensity from the microarray of the Slc25a5 gene for each animal is displayed against the ΔCT value from the RT-qPCR analysis. The fold change value for the DSG4 protein was plotted against the absorbance from the ELISA analysis. The R2 value was calculated using linear regression.
https://doi.org/10.1371/journal.pone.0244098.g004 transcripts were associated with signal transduction, systemic regulation, ontogenesis, lipid metabolism, cellular integrity and reproduction, and cellular integrity and immune response, respectively (Fig 2A). The GO terms associated with the seven exposure-related protein candidates were mainly associated with cellular integrity (general) and supramolecular maintenance ( Fig 2B). The associated GO terms for the eleven proteins in thyroid and plasma related to dose were mainly involved in apoptotic cell death, cell death, cellular integrity, metabolism, nucleic acid related metabolism, ontogenesis, supramolecular maintenance and systemic regulation (Fig 3). Furthermore, the majority of the proteins obtained in both thyroid tissue and plasma for the same activity level were involved in cellular integrity and metabolism (Table 2). IPA was used to analyse canonical pathways, functions and diseases and upstream regulators for transcripts and proteins in thyroid and proteins in plasma. In general, significantly affected canonical signalling pathways, functions and diseases and predicted upstream regulators were predominantly associated with proteins found in thyroid tissue, with pathway activation (z-score �2) being most prevalent (Tables 4-6). Pathway activation/inhibition was common for proteins in plasma in the group with the lowest activity (0.5 kBq). The majority of the signalling pathways were only significantly affected in one of the groups (Table 4). Several other signalling pathways were common in two groups, including actin cytoskeleton signalling, B cell receptor signalling and HGF signalling (0.5 and 50 kBq groups), integrin signalling (5 and 50 kBq groups) and calcium signalling (50 and 500 kBq groups). The NRF2-mediated oxidative stress response was common in three groups (5, 50 and 500 kBq). For the biological functions and diseases analyses, the majority of activated functions were seen in the 5 and 50 kBq groups and the total number of associated functions was similar for two groups (Table 5). Few up-stream regulators were significantly regulated, with only one or  Table 4. Predicted canonical signalling pathways affected in rat thyroid and plasma after exposure to low doses of 131 I. All analyses were performed for all twelve samples and the only obtained canonical signalling pathway for transcripts was "Dopamine Receptor Signalling" for the 0.5 kBq group. The direction of regulation of each target molecule is marked by bold for increased expression and underscore for decreased expression.   two for each group. Only the PPARG was identified as an up-stream regulator in all four groups (Table 6).

Discussion
In this study, we investigated long-term (9 months) effects of low to medium dose 131 I exposure on expression patterns in thyroid tissue and plasma collected from young rats. The gene and protein expressions were compared with corresponding data from non-exposed controls of the same batch and age. Then the expression patterns were compared with absorbed dose, and potential relation to thyroid function and cancer were examined for the identified genes and proteins. However, biomarker discovery is not easy as ideal biomarkers must fulfil several criteria. Biomarkers for exposure should give unidirectional differential expression levels, Table 5. Diseases or functions annotation of thyroid and plasma using IPA. Four, six, six and one diseases or functions were seen for 0.5, 5, 50 and 500 kBq groups, respectively. All diseases or functions were activated except for "fatty acid metabolism" and "apoptosis". All but three diseases or functions, "fatty acid metabolism", "neuronal cell death" and "apoptosis" were obtained in the thyroid proteomic analysis. The direction of regulation of each target molecule is marked by bold for increased expression and underscore for decreased expression. while biomarkers for dose relationships should have monotonic differential expression levels that correlate to absorbed dose. For thyroid function, a biomarker needs to be thyroid-specific and dependent on cellular response. Lastly, tumour biomarkers must be a) detectable at an early stage in cancer development, b) absent in healthy individuals, and c) tumour marker abundance should correlate directly with tumour size [26]. In addition, biomarkers should also ideally be detectable in body fluids, which would simplify sample collection and be less invasive for patients. In total, we found nine potential exposure-related biomarkers in thyroid (transcripts: Afp, RT1-Bb, and proteins: ARF3, DLD, IKBKB, NONO, RAB6A, RPN2 and SLC25A5), and eleven differentially expressed dose-related proteins (in thyroid; AGRN, APRT, DSTN, HADHA, LDHA, MYBPC3, NOL3, SLC3A2 and YBX3; in plasma; DSG4 and TGM3). The Afp gene had lower expression after 131 I exposure in thyroid. Afp is highly expressed in foetal liver cells, foetal yolk sac and the gastrointestinal tract [27]. The Afp gene is also involved in the regulation of cellular growth and its mRNA has been detected in both human PTC and in normal thyroid [28]. Elevated expression of AFP has been found in radiation-induced hepatocellular carcinomas [29]. The RT1-Bb gene had higher expression after 131 I exposure in thyroid, and belongs to the MHC II class of proteins that are involved in the immune system [30]. All seven exposure-related proteins in thyroid tissue were overexpressed after irradiation. The ARF3 protein is GTP binding and involved with protein trafficking and vesicle binding, the Arf3 has also been suggested as a marker in radiation transformed breast cancer [31,32]. The DLD protein is a part of several multi-enzyme complexes involved in energy metabolism, and it has also been proposed as a biomarker for UV radiation induced skin aging [32,33]. IKBKB has an important role in the NF-kappa beta signalling pathway, which is triggered by radiation and can cause resistance for irradiated cancer cells [34]. NONO is an RNA binding protein that is part of transcription and RNA splicing, and has been associated with UV radiation induced cellular damages [35,36]. RAB6A is located in the Golgi apparatus and is involved in the targeting and fusion of transport carriers [32,37]. RPN2 is a membrane protein that mediates Table 6. Upstream regulators for thyroid and plasma using IPA. PPARG was activated in all four groups (0.5-500 kBq) found in thyroid tissue protein analysis. TSC2 was activated for the 0.5 kBq group found in protein analysis of plasma. In the protein analysis of thyroid tissue TNF and NFKB (complex) were inhibited and activated for the 5 and 50 kBq group, respectively. The direction of regulation of each target molecule is marked by bold for increased expression and underscore for decreased expression. protein translocation in the endoplasmic reticulum and the SLC25A5 protein transports ADP from the cytoplasm to the mitochondria and ATP in the opposite direction [32]. All eleven dose-related proteins, nine in thyroid tissue and two in plasma, were overexpressed after 131 I exposure. Among those, the most promising response relationship was found for LDHA and APRT in thyroid tissue (0.5-500 kBq) and TGM3 and DSG4 in plasma (0.5-50 kBq). LDHA catalyses the last step in aerobic glycolysis and is primarily located in muscle tissue, and APRT is an enzyme involved in the adenine nucleotide metabolism [32]. Furthermore, LDHA overexpression has been found in BRAF V600E-mutated thyroid cancer and inactivation of this protein can make cells more sensitive to radiation [38,39]. TGM3 is involved in the transformation of the cell envelope in epidermis and hair follicles, but it is also a tumour suppressor and has been proposed as a progression marker for radio-chemotherapy of head and neck cancers [32,40]. DSG4 is involved in cell adhesion, primarily in epithelial cells [32].

Diseases or Functions Annotation
As far as we know, data on short-or long-term effects on rats after 131 I exposure have not been published. When comparing our present long-term data in young rats with publicly available short-term data in adult mice, few common differentially expressed transcripts and proteins were identified [13][14][15][16]. Interestingly, the Klk1 (0.5 kBq thyroid), CPA3 (0.5 kBq thyroid), and S100A8 (50 kBq plasma) markers were commonly identified in both the present study and previous studies on the short-term effects of 131 I exposure on thyroid. Klk1, CPA3 and S100A8 are important for thyroid function in both mice and rats [19,20]. Genome homology between humans and rodents is high (almost 85% at the expression level) and they share several cellular mechanisms [41]. We believe that comparing results generated between species despite differences in species, age and time is valuable for the transfer of the generated results to humans. Furthermore, identified markers between species with similar function are conserved during evolution reflect important function. Therefore, the commonly identified markers in rats and previously in mice adult and young as well as persist with time after exposure initiation, might be potential exposure-related biomarkers also for humans.
The ACADL and SORBS2 protein, were regulated at all dose levels studied and are the most promising candidate biomarkers, since they have been previously associated with thyroid function [19,20] Among the thyroid-specific proteins, elevated levels of thyroid peroxidase (TPO) were found in the three groups with lowest administered activities (0.5-50 kBq). This finding is in contrast to results from a study using 60 Co gamma irradiated cultured human thyroid epithelial cells, where TPO expression levels were reduced, possibly indicating short-term effects in vitro [42]. Since the TPO and TG thyroid associated proteins were differentially expressed in some of our exposed groups (primarily at low doses), these proteins might be potential biomarkers for exposure at low dose-levels since no monotonic dose-response relationships were found.
In total, five transcripts (Bhlhe41, Calca, Calcb, Dio2 and Krt7) related to thyroid function were found in some of the groups. Bhlhe41 is related to circadian rhythm [32]. The Calca gene codes for calcitonin production, which is involved in the calcium regulation and phosphor metabolism [32]. Calcb is a paralog to Calca and is associated with medullary thyroid carcinoma [32]. The Dio2 (Deiodinase 2) gene is highly associated with thyroid function and encodes for a protein that catalyses the deiodination of T4 (thyroxine) to the active thyroid hormone T3 (triiodothyronine). The DIO2 protein is also found in several tissue types and is thought to be a local producer of T3 [32]. Overexpression of the Dio2 gene has been seen in follicular adenoma and Grave´s disease [32]. The Krt7 gene belongs to the keratin family and stimulates DNA synthesis in cells by blocking interferon-dependent interphase [32].
In total, eight proteins associated with thyroid function were identified in thyroid tissue and/or plasma (ACADL, BCL2, CALCA, CYP27B1, PLCB4, SORBS2, TG and TPO). The ACADL and SORBS2 proteins were overexpressed in all groups in thyroid tissue. ACADL is an enzyme involved in the metabolism of fatty acids and amino acids, while SORBS2 belongs to a family of non-receptor protein-tyrosine kinases [32]. Thyroglobulin (TG) and thyroid peroxidase (TPO) are thyroid-specific proteins. TG is a protein primarily produced by the thyroid gland, involved in the production of the iodine-containing hormones T3 and T4, and the storage of T4 [32]. TG is clinically relevant, especially to predict thyroid residues, since TG blood levels are related to thyroid function, rather than thyroid malignancy [43]. Elevated TG expression levels can therefore indicate that thyroid function is impaired after exposure. TPO is an enzyme essential for thyroid gland function involved in T4 and T3 production, and mutations in this gene can lead to hyperthyroidism [32]. Of the observed proteins, the strongest association with thyroid function was seen for TG and TPO. These proteins were not differentially regulated in all groups, but upregulated at the two lowest activities (0.5 and 5 kBq) and at the highest activity (500 kBq). This is interesting since it suggests that thyroid function after irradiation is not only affected by dose, but other biological effects as well, as previously suggested iodine deficiency, age at exposure, gender and the issue with increased screening [44].
For clinical use, biomarkers in plasma would be preferable, and only four proteins were differentially expressed regulated among those known to be related to thyroid function (TG, CYP27B1, PLCB4 were overexpressed and CLIP2 was underexpressed). TG is interesting, but is a protein that might be influenced by general physiological status. CYP27B1 is involved in transforming vitamin D to its active form, and PLCB4 is an enzyme, which transduces signals over the plasma membrane by using second message molecules diacylglycerol (DAG) and inositol 1,4,5-trisphosphate (IP3), both of which have functions that are important for normal thyroid and body function [32]. CLIP2 has been proposed as a biomarker for thyroid cancer in several studies on children in the Chernobyl cohort [45][46][47]. In the present study, we identified decreased CLIP2 expression in exposed thyroid tissue and plasma (all doses), while Selmansberger et al. showed an increase in both protein and mRNA expression levels in thyroid tumour tissue in children [48].
The higher association of transcripts than proteins with GO terms is most probably explained by investigational bias, since transcriptomic profiling is hitherto far more common than proteomic analyses, with less data for protein enrichment. A pronounced association with metabolism, immune system, apoptosis and cell death was established for all samples. Increased metabolism is caused by increased thyroid hormone production, and might suggest abnormal thyroid function (hyperthyroidism). The strong association with apoptosis and cell death may be due to thyroid cell impairment and function. The immune system is highly involved in cancer and other diseases. All these effects indicate that thyroid cells are affected by irradiation (observed at as low as 1 mGy), and that the thyroid function differs from that of non-irradiated thyroid.
The significant canonical signalling pathways that were commonly identified in more than one group are all necessary for functions related to cellular survival. NRF2-oxidative stress leads to the binding of antioxidant responsive elements and can thus neutralise reactive molecules [32]. The upstream regulator PPARG is involved in differentiation, cell growth and nuclear regulation. It is primarily found in adipose tissue, is associated with diabetes 2 and colon cancer but is also related to follicular thyroid cancer and abnormal thyroid function [32,49]. It is also interesting to note that activation of PPARG was not influenced by dose-level. We also previously showed downregulation of PPARG expression in our short-term studies on rats [14].
Other studies on 131 I-induced effects in vitro and in vivo have focused on specific effects on non-coding RNAs and selected proteins. Long non-coding RNA MEG 3 has been a suggested as a biomarker for thyroid cancer and 131 I radio-resistance [50]. Proteomic studies on PTC tissue from humans have suggested NTRK1, metalloproteinases (MMP-1, MMP-9 and MMP-13) and Cathepsins (-W and -X), NF-KB and other apoptosis associated proteins as radiationinduced PTC biomarkers, together with SPANXA as an important protein for thyroid cancer development [51][52][53][54][55][56]. However, none of these were found in the present study. In vitro studies on thyroid cells have shown that DNA damage after 131 I exposure increases in the presence of high levels of TSH [57]. Thyroid function has also been evaluated, by measurements of the IL-6 TNF-αand NO proteins in serum samples using a rat model using four week old rats irradiated with 5.5 MBq 131 I [58]. Some of these proteins were found also in the present study, but the expression levels were not in the same direction (increased or decreased), probably due to the longer time after administration and the much lower amount of 131 I administered.
In the present study, several regulated genes and proteins for exposure and dose-response are presented. However, no single biomarker was found that seems to fulfil the necessary criteria for being an optimal biomarker related to dose response, with clearly differentiated expression level. These findings are in line with our previous studies on the short-term effects of exposure to ionizing radiation in adult mice and rats, where no clear dose-response relationships were found for any single transcript [14,59]. We and others then suggested the use of a panel of several genes to evaluate radiation response [11,14,16]. This situation seems to also be the case for long-term effects in young rats.

Conclusions
In the present study, the gene and protein expression profiles for rat thyroid tissue and plasma were investigated in young rats nine months after exposure to low amounts of 131 I. The results showed that only a few genes or proteins altered their expression in a clear dose-related manner. Expression patterns mostly demonstrated altered (up and down) regulation at certain dose levels in a non-monotonic manner. The most promising dose-related candidate biomarkers were the LDHA and APRT proteins (in thyroid), and TGM3 and DSG4 proteins (in plasma, only for low doses). The most promising exposure-related candidate biomarkers were the Afp and RT1-Bb transcripts in thyroid, both showing uniform differential expression patterns, irrespective of dose level.
The GO term association analysis revealed that most of the regulated genes/proteins were associated with immune response, which is influenced in many diseases, including cancer. For the transcripts, a pronounced association with cellular processes related to cell cycle regulation, metabolism and cancer was found.
Thyroid function-related proteins overexpressed in thyroid and/or plasma included SORBS2, ACADL, TG, and TPO (low doses). CLIP2, which has previously been found to be overexpressed in children with PTC in the Chernobyl cohort, was underexpressed in thyroid and plasma in the present study on young rats.
Preferably, biomarker examinations should be performed in blood samples for less invasive tissue sampling. However, in the present study the majority of the candidate biomarkers were found only in thyroid tissue, and thyroid needle biopsies are routinely collected clinically and thus possible to collect.
In conclusion, this work resulted in increased knowledge of late radiation-induced cellular effects after low to medium dose exposure from 131 I. Novel candidate biomarkers related to exposure or absorbed dose were identified and previously proposed candidate biomarkers were not detected in this work. The identified biomarkers need further validation using human tissue samples.
Supporting information S1 Table. Uniquely regulated a) transcripts in thyroid tissue b) proteins in thyroid tissue and c) proteins in plasma for each group. (DOCX) S2 Table. Transcripts and proteins in thyroid and plasma tissue with statistically significant differential expression levels for three of four groups. (DOCX) S3 Table. Transcripts and proteins in thyroid and plasma tissue with statistically significant differential expression levels for two out of four groups. (DOCX)