Comparative Transcriptome Analysis Reveals the Potential Cardiovascular Protective Targets of the Thyroid Hormone Metabolite 3-Iodothyronamine (3-T1AM)

Background The thyroid hormone metabolite 3-iodothyronamine (3-T1AM) is rapidly emerging as a promising compound in decreasing the heart rate and lowering the cardiac output. The aim of our study was to fully understand the molecular mechanism of 3-T1AM on cardiomyocytes and its potential targets in cardiovascular diseases. Materials and Methods In our study, we utilized RNA-Seq to characterize the gene expression in H9C2 cells after 3-T1AM treatment. Comparative transcriptome analysis, including gene ontology, signaling pathways, disease connectivity analysis, and protein-protein interaction networks (PPI), was presented to find the critical gene function, hub genes, and related pathways. Results A total of 1494 differently expressed genes (DEGs) were identified (192 upregulated and 1302 downregulated genes) in H9C2 cells for 3-T1AM treatment. Of these, 90 genes were associated with cardiovascular diseases. The PPI analysis indicated that 5 hub genes might be the targets of 3-T1AM. Subsequently, eight DEGs characterized using RNA-Seq were confirmed by RT-qPCR assays. Conclusions Our study provides a comprehensive analysis of 3-T1AM on H9C2 cells and delineates a new insight into the therapeutic intervention of 3-T1AM for the cardiovascular diseases.


Introduction
3-Iodothyronamine (3-T1AM) is a metabolite product of the thyroid hormones (TH) which exerts an opposite physiological function of classical TH such as a decrease of body temperature and metabolic rate in rodents [1]. Several studies have shown that 3-T1AM is presumed to act on various organs and different receptors rapidly [2]. Moreover, it has therapeutic potential for treatment of metabolic diseases [3], myocardial vascular diseases [4], and neurological diseases [5].
Thyroid hormones and their metabolite products play a critical role in the development of cardiovascular dis-eases [6][7][8]. Of note, 3-T1AM, a decarboxylated and deiodinated thyroid hormone derivative, has displayed a potency to decrease the heart rate and lower the cardiac output [9], but the underlying mechanism between 3-T1AM and myocardial vascular diseases has remained unknown. It would therefore be of interest to perform RNA-Seq assay to predict the potential transcriptional mechanism used by 3-T1AM to protect from the myocardial injuries. In our present study, we aim to identify differently expressed genes in H9C2 cells and cells treated with 3-T1AM and investigate the underlying molecular mechanism of its effect on myocardial protection.
2.2. Cell Viability Assessments. The cell viability was assessed by CCK-8 assay. Briefly, cells were seeded into a 96-well culture plate at a density of 6 × 10 3 /well. Cells were treated with various doses of 3-T1AM for 24 h. H9C2 cells cultured in 3-T1AM-free growth medium were used as controls. Then, cells in each well were incubated with CCK-8 at 37°C for 4 hours. The absorbance was determined by a plate reader.
2.3. Total RNA Extraction. Total RNA in both groups (control and 3-T1AM treatment) was extracted with RNAiso Plus (Takara) according to the manufacturer's instructions. The cDNA library was constructed based on the instructions of Illumina. RNA-Seq of each samples was performed by Sangon Biotech Co. Ltd. (Shanghai). Methods of RNA-Seq analysis were described in previous studies [10].

Construction of RNA-Seq Library for Illumina
Sequencing. Sequencing library preparations were constructed using the manufacturer's protocol (VAHTSTM mRNA-seq V2 Library Prep Kit for Illumina®). The digested RNA samples were used for firstand second-strand cDNA synthesis with random hexamer primers. After the firstand second-strand cDNA synthesis, the double-stranded products were end repaired, a single "A" base was added, and adaptors were ligated onto the cDNA products. cDNA libraries were constructed. Then, the paired-end sequencing was used on an Illumina HiSeq™ 2500 [11].

Identification of Differentially Expressed Genes (DEGs).
After sequencing, raw data were obtained in the fastq format. FastQC was used for validating the quality of the data. Trimming of sequences was performed using Trimmomatic. Then, quality check using FastQC was performed again on the trimmed sequences. Reads per kilobase of exon mode per million mapped reads (RPKM) were employed to quantify transcript levels [12]. Differentially expressed genes in four samples were sought, and conditions for differences in genes were filtered. The adjusted P value < 0.05 and Log2ð fold changeÞ > 1 were set as the cutoff criterion.
2.6. Gene Ontology and KEGG Pathway Enrichment Analysis. The Database for Annotation, Visualization and Integrated Discovery (DAVID 6.8) gene annotation tool (http://david .ncifcrf.gov/tools.jsp) was used to perform GO and KEGG analysis [13]. Gene ontology term analyses of the DEGs, including biological process (BP), cellular component (CC), and molecular function (MF), were identified using the gene ontology (GO) project (http://www.geneontology.org). KEGG pathway enrichment analysis was performed based on the Kyoto Encyclopedia of Genes and Genomes database (http://www.kegg.jp) [14]. In addition, GO terms and pathway with P values less than 0.05 were considered to be statistically significant.

Construction of Protein-Protein Interaction (PPI)
Network. The PPI network was generated based on STRING 10.0 (http://String-db.org) [15]. A confidence score of 0.9 was set as the cutoff criterion. All of the networks were visualized and analyzed with Cytoscape 3.4.0 (http://cytoscape.org). Moreover, the significant network modules with degree cutoff = 2, node score cutoff = 0:2, and k − score = 2 were constructed using the Molecular Complex Detection (MCODE) plugin for Cytoscape.
2.9. Quantitative Real-Time PCR (RT-qPCR). In order to confirm the DEGs of RNA-Seq, the expression levels of 8 genes, such as Camk2d, Ppp2Ca, Nf1, Yes1, Psma6, Psmb6, Jak2, and Sirt4, were chosen and analyzed by RT-qPCR. Total RNA in both groups (control and 3-T1AM treatment) was extracted with RNAiso Plus (Takara) according to the manufacturer's instructions. RT-qPCR was conducted by using the Vazyme HiScript II One Step qRT-PCR SYBR Green Kits (Vazyme, Nanjing) on a 7500 real-time PCR system (ThermoFisher Scientific). GAPDH was used as a control for normalization of RT-qPCR results. The sequences of primers designed by NCBI Primer Blast are shown in Supplementary Table S1. Three independent replicates were conducted for this experiment. The fold change was calculated by 2−ΔΔCt.

Illumina Sequencing and Preprocessing of Raw Data.
RNA sequencing was performed on poly(A)-enriched RNA extracted from H9C2 cells and H9C2 cells treated with 3-T1AM. The raw data of sequencing was evaluated and is shown in Table 1. The Q30 data, the probability of an error base call 1 in 1000 times, were over 91%. The proportion of Q20, the probability of an error base call in 1 in 100 times, was more than 96%, which suggests that the quality of sequencing data was appropriated for the following data analysis.

Identification of Differentially Expressed Genes
Associated with 3-T1AM Treatment. The transcriptome profiles of H9C2 cells treated with 3-T1AM were compared with those of H9C2 cells. Differently expressed genes (DEGs) were identified according to Log2 fold changes and P value. Volcano plots were used to analyze DEGs, in which red spots and green spots represent upregulated and downregulated genes separately (Figure 2  KEGG enrichment analysis was performed to identify relative pathway in which the 1494 DEGs were involved. A scatter plot was used to depict significantly enriched pathway terms. Pathways such as RNA transport, P53 signaling pathway, and ribosome were most significantly enriched terms in H9C2 cells treated with 3-T1AM compared with the control group( Figure 2(c)).

Disease Connectivity Analysis.
We conducted an additional analysis of the DEGs using the Comparative Toxicogenomics Database (CTD), which permits the development of novel hypotheses about the relationships between 3-T1AM and cardiovascular diseases. To reduce interference of unrelated genes and identify the cardiovascular diseaseassociated genes, 90 genes were determined as hub genes ( Table 2). In addition, we also verified that 68 genes were associated with heart diseases, 7 genes were associated with MELAS syndrome, 50 genes were associated with vascular diseases, and 22 genes were associated with cardiomyopathies.
3.6. PPI Network Construction. In order to further analyze a protein-protein interaction network of these 90 DEGs associated with cardiovascular diseases, a PPI network was generated based on the STRING database and Cytoscape 3.4.0 software. Our results indicated that a module that consisted     BioMed Research International of 16 nodes was identified from 90 genes associated with cardiovascular diseases (Figure 3(a)). By using the MCODE analysis, the module that consisted of 5 nodes and 10 edges (MCODE score = 4) was screened out (Figure 3(b)).

Verification of Differentially Expressed Genes by RT-qPCR.
To verify the reliability of the RNA-Seq results, eight cardiovascular-related candidates were chosen in both groups to verify differential expression. The results strengthened that variation tendencies between the RNA-Seq and RT-qPCR data were identical, verifying the reliability of the RNA-Seq data (Figure 4).

Discussion
Previous studies have indicated that 3-T1AM is an effective metabolite product of the TH that can decrease the heart rate and lower the cardiac output [16][17][18]. We speculated that cardiomyocytes would be an essential target site for 3-T1AM. RNA-Seq profile analysis is an efficacious assay for uncovering the potential molecular mechanisms of 3-T1AM and has been extensively accepted as a method to find the therapeutic targets of drugs [19]. In our investigation, we explored the DEGs and mechanisms of cultured H9C2 cell lines after 3-T1AM treatment using RNA-Seq and appropriate transcriptome analysis.
Our research substantiated that 3-T1AM can reduce the cell viability of H9C2 cells in a time-and dose-dependent manner using a CCK8 assay. Studies showed that doses of 3-T1AM used in in vitro studies were different in various cell lines [3,[20][21][22]. Then, in our study, 20 μM 3-T1AM that could not affect cell viability after treatment for 24 h was chosen for RNA-Seq. A limitation of our study is that time-and concentration-related responses cannot be monitored according to the results of RNA-Seq.
A total of 1494 DEGs were identified in cultured H9C2 cells exposed to 3-T1AM for 24 h based on the screening criteria of DEGs in our study. Among them, only 90 DEGs associated with cardiovascular diseases were identified. To confirm the reliability of the RNA-Seq results, a number of DEGs, such as Camk2d, Ppp2ca, Nf1, Yes1, Psmb7, Sirt4, and Jak2, were selected and tested by RT-qPCR, which indicated that these two methods were in good agreement.
The comparative analysis of DEGs provides the therapeutic intervention for the cardiovascular diseases and delineates novel insight into other diseases. Several differentially expressed genes identified in this study, including HMGB1, ROCK2, and Jak2, are amenable in principle to diagnostic biomarkers and therapeutic targets. Previous studies indicated that upregulation of HMGB1 aggravated mechanical stress-induced cardiomyocyte hypertrophy via the RAGE/ERK1/2 signaling pathway [23] [24]. Moreover, upregulation of Jak2 activated the Jak2-STAT3 pathway, thereby alleviating myocardial ischemiareperfusion injury, cardiomyocyte apoptosis, and ROS production in vitro [25].

BioMed Research International
It is widely known that 3-T1AM exerts a physiological function of lowering body temperature in rodents. Studies suggest that 3-T1AM acts through sirtuin-mediated pathways to metabolically reprogram fatty acid and glucose metabolism possibly through small-molecule signaling [26]. To understand biological functions of the 3-T1AM-induced DEGs in H9C2 cells, several GO terms which are related to metabolism were explored including the biological process category "metabolic function" and the molecular function category "catalytic activity". In addition, glutathione metabo-lism, pyrimidine metabolism, and purine metabolism were found by KEGG analysis. These finding suggested that 3-T1AM could affect cardiomyocytes by metabolic function.
Previous studies found that hormone-binding domain was constituted by two binding sites, one of which interacts with thyroid receptors, resulting in a decreased transcriptional activity of numerous proteins via phosphorylate p53 [27]. Moreover, the p53 signaling pathway was also found by KEGG analysis of 3-T1AM-induced DEG. These results indicated that 3-T1AM is involved in the p53 signaling pathway, but whether it is involved in the p53 pathway needs to be further studied.

Conclusion
In conclusion, we identified 1494 DEGs in H9C2 cells for 3-T1AM treatment in our research. Of these, 90 genes had associated with cardiovascular diseases. Interestingly, we further employed comparative analysis of biological functions, related signaling pathways, and PPI network of these identified DEGs. The 5 hub genes in cardiovascular diseases were identified according to the PPI network analysis. In addition, our results provide a new insight into the therapeutic intervention of 3-T1AM for the cardiovascular diseases.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no competing interests.