Abstract Context Despite the well-recognized clinical features resulting from insufficient or excessive thyroid hormone (TH) levels in humans, it is largely unknown which genes are regulated by TH in human tissues. Objective To study the effect of TH on human gene expression profiles in whole blood, mainly consisting of T3 receptor (TR) α-expressing cells. Methods We performed next-generation RNA sequencing on whole blood samples from eight athyroid patients (four females) on and after 4 weeks off levothyroxine replacement. Gene expression changes were analyzed through paired differential expression analysis and confirmed in a validation cohort. Weighted gene coexpression network analysis (WGCNA) was applied to identify thyroid state-related networks. Results We detected 486 differentially expressed genes (fold-change >1.5; multiple testing corrected P value < 0.05), of which 76% were positively and 24% were negatively regulated. Gene ontology (GO) enrichment analysis revealed that three biological processes were significantly overrepresented, of which the process translational elongation showed the highest fold enrichment (7.3-fold, P = 1.8 × 10−6). WGCNA analysis independently identified various gene clusters that correlated with thyroid state. Further GO analysis suggested that thyroid state affects platelet function. Conclusions Changes in thyroid state regulate numerous genes in human whole blood, predominantly TRα-expressing leukocytes. In addition, TH may regulate gene transcripts in platelets. Thyroid hormone (TH) affects virtually all cells and tissues in the human body. The major biologically active TH is T3, and its genomic actions are mediated by binding to nuclear T3 receptors (TRs) that regulate transcription of target genes (1). Several receptor isoforms are encoded by the THRA and THRB genes of which TRα1, TRβ1, and TRβ2 are the truly T3 binding isoforms (2). The TR isoforms have a distinct expression pattern, with a predominance of TRα1 in brain, heart, and bone and TRβ1 in the liver, kidney, and thyroid. TRβ2 is mainly expressed in the hypothalamus and pituitary and is therefore involved in the regulation of the hypothalamus-pituitary-thyroid axis (3). Despite the classical clinical features resulting from insufficient or excessive TH levels that have been recognized for more than 100 years, the underlying molecular mechanisms in humans are not well understood. Knowledge of gene expression modulated by TH is largely derived from animal models or in vitro cellular studies used to explore which genes are regulated by TH (4–7). Expanding knowledge of the effects of TH on gene expression in human tissues will provide more insight in the molecular basis of TH action and may lead to a better understanding of the clinical effects of TH in humans. Progress is limited because most human tissues are not easily accessible. However, blood can be regarded as circulating tissue and contains various cell types including erythrocytes, leukocytes, and platelets (8). RNA in whole blood is largely determined by leukocytes. Peripheral blood mononuclear cells (lymphocytes and monocytes) have been shown to mainly express the TRα isoform (9). Therefore, analysis of gene expression in whole blood may potentially be used as proxy for other TRα-expressing tissues. To study the effects of TH on gene expression in human TRα-expressing cells, we performed next-generation RNA sequencing (RNA-seq) in whole blood cells from athyroid patients off and on levothyroxine (LT4) treatment. Patients and Methods Patients Patients were recruited via the outpatient clinic of the Erasmus Medical Center, which is a tertiary referral center for differentiated thyroid cancer (DTC). Patients with DTC undergoing thyrotropin (TSH)-stimulated 131I therapy after withdrawal of LT4 were asked to participate in the study. Patients were eligible for inclusion if they had no other malignancies or an active inflammatory disease, were not using any drugs known to influence TH metabolism, and were between 18 and 80 years of age. A discovery and a validation cohort were created according to the same protocol. The study protocol was approved by the Medical Ethics Committee of the Erasmus Medical Center (MEC 2012-561). Written informed consent was obtained from all study participants. Sample collection and serum analyses Peripheral blood samples were obtained from all participants when overt biochemical hypothyroidism was achieved by withdrawal of LT4 substitution in thyroidectomized patients and when TSH suppression was achieved after restarting LT4 replacement therapy. Serum free T4 (reference range, 11 to 25 pmol/L), total T4 (reference range, 58 to 128 nmol/L), and total T3 (reference range, 1.4 to 2.5 nmol/L) concentrations were measured by chemoluminescence assays (Vitros ECI Immunodiagnostic System; Ortho-Clinical Diagnostics, Rochester, MI). Serum TSH (reference range, 0.4-4.3 mU/L) was measured by an immunometric assay (Immulite 2000 XPi; Siemens, Den Haag, the Netherlands). Whole blood samples were collected in PAXgene tubes (Hombrechtikon, Switzerland). PAXgene tubes contain a proprietary reagent that immediately stabilizes intracellular messenger RNA (mRNA), thus reducing mRNA degradation and inhibiting gene induction after phlebotomy. The isolated mRNA represents all blood cells, including polymorph nuclear leukocytes, mononuclear cells, platelets and red blood cells (10). RNA was isolated using PAXgene blood miRNA kits (PreAnalytiX, Hombrechtikon). For RNA-seq analysis of the discovery cohort, ribosomal RNA and globin mRNAs were removed from an aliquot of the RNA samples using the Globin-Zero Gold rRNA Removal Kit (Illumina, San Diego, CA) (11). Kapa Stranded RNA library was prepared (Kapa Stranded RNA kit; Kapa Biosystems, Wilmington, MA), followed by sequencing on a HiSeq2500 system, for single end reads, 50 bp in length (Illumina, San Diego, CA). Bioinformatic analysis We analyzed genes with at least five reads in at least six samples. The generated sequencing reads were aligned (stranded alignment) against the GRCh38 version of the human reference genome, with RefSeq gene annotation using Tophat2 (12). Gene counts were generated from the BAM files with HTsEquation (13). Cufflinks was used to compute transcript abundance estimates in fragments per kilobase of transcript per million mapped reads (FPKMs) (14). Differential expression was calculated with the bioinformatics tool DESEquation 2 from Bioconductor, which uses the R statistical programming language (15). Cutoff values for significantly expressed genes were a false discovery rate of 0.05 or less and a fold change of 1.5. To visualize the clustering of the samples, principal component analysis was performed. The normalized data file was transposed and imported into OmniViz, version 6.1.13 (Instem Scientific, Inc., Stone, Staffordshire, UK), for further analysis. The geometric mean of the intensity of all samples was calculated. The level of expression of each gene was determined relative to this geometric mean and 2log transformed. The geometric mean of all samples was used to ascribe equal weight to gene expression levels with similar relative distances to the geometric mean. The data were deposited in National Center for Biotechnology Information/Gene Expression Omnibus reference (https://www.ncbi.nlm.nih.gov/geo/) with accession number GSE103305. We used functional enrichment analysis to explore the biological significance of differentially expressed (DE) genes and of weighted gene coexpression network analysis (WGCNA) module genes (see the following section) using Database for Annotation, Visualization, and Integrated Discovery (Expression Analysis Systemic Explorer) (16), which calculates significantly overrepresentation of gene ontology (GO)–classified biological processes by comparing the number of genes in a gene list for a biological process to the number of genes for that biological process from the RNA-seq analysis (17). Biological processes are shown that were significantly enriched after correction for multiple testing. The independent skeletal muscle dataset for comparative analysis consisted of 10 thyroidectomized patients with DTC 4 weeks after LT4 withdrawal (before TSH-stimulated 131I scintigraphy) and 8 weeks after subsequent LT4 replacement (18). To study the expression of the TR isoforms in different blood cell types, independent datasets on RNA-seq in platelets (National Center for Biotechnology Information Sequence Read Archive SRP028846), lymphocytes (E-MTAB-2319), and neutrophils (Gene Expression Omnibus Series GSE59528) were collected from previous studies (19–21). FPKM values of TRα1 and TRβ1 in the different cell types were compared using a t test. Quantitative polymerase chain reaction To confirm the RNA-seq results by an independent technique, quantitative polymerase chain reaction (qPCR) was used to measure expression of selected genes. qPCR was performed on complementary DNA produced from RNA before applying Globin-Zero Gold rRNA Removal Kits on the discovery cohort and on complementary DNA from the validation cohort. The primer sequences are presented in Supplemental Table 1. RNA levels are expressed relative to the geometric mean of the housekeeping genes GAPDH and ACTB. Statistical analysis Likelihood and significance of the overlap with DE genes in muscle samples off and on LT4 treatment previously reported were tested using the χ2 test (18). Thyroid function tests and leukocyte counts were compared off and on LT4 treatment using the Wilcoxon signed-rank test. FPKM values of the TRα1 gene were compared off and on LT4 using a paired t test. Statistical analysis was performed using SPSS Statistics for Windows, version 22 (IBM Corp., Armonk, NY). WGCNA We performed a WGCNA to discover coexpressed TH-related networks (modules) (22, 23). In short, a weighted adjacency matrix containing pairwise connection strengths was constructed by using the soft-thresholding approach (β = 6) on the matrix of pairwise correlation coefficients. A connectivity measure (k) per gene was calculated by summing the connection strengths with other genes. Modules were defined as branches of a hierarchical clustering tree using a dissimilarity measure (1, topological overlap matrix) (23, 24). Each module was subsequently assigned a color. The gene expression profiles of each module were summarized by the module eigengene (defined as the first principal component of the module expression levels). Each module eigengene was correlated to thyroid status (with age, body mass index, sex, leukocytes, and leukocyte differentiation as covariates) using the WGCNA R package (false discovery rate, 5%). Results Study population We studied whole blood samples from eight athyroid patients (four females and four males) off and on LT4 treatment. Characteristics of the patients are shown in Supplemental Table 2. Table 1 shows the serum TH concentrations and total leukocyte counts, including leukocyte differentiation, off and on LT4 treatment. As expected, the thyroid function tests were significantly different. Neither total leukocyte counts nor leukocyte differentiation showed substantial differences, except for a slight increase in monocyte number during LT4 treatment. Table 1. Thyroid Function Tests and Leukocyte Counts Off LT4 On LT4 P TSH (0.4–4.3) (mU/L) 78.0 (48.5–117.5) 0.07 (0.004–0.11) 0.012 Total T3 (1.4–2.5) (nmol/L) 0.72 (0.63–0.83) 2.11 (1.94–2.28) 0.012 Total T4 (58.0–128.0) (nmol/L) 17.0 (14.0–18.0) 155.0 (135.5–170.5) 0.018 Free T4 (11.0–25.0) (pmol/L) 1.55 (1.12–2.17) 26.2 (22.3–29.5) 0.012 Leukocytes (3.5–10.0) (×109/L) 7.6 (5.5–8.1) 7.4 (5.5–8.7) 0.26 Neutrophils (40.0–80.0) (%) 60.0 (52.0–67.0) 62.8 (58.2–65.7) 0.50 Lymphocytes (15.0–50.0) (%) 32.4 (25.1–35.9) 27.3 (23.5–28.6) 0.13 Monocytes (5.0–14.0) (%) 4.7 (4.6–7.8) 8.1 (6.9–9.6) 0.018 Dosage LT4 (µg/kg) 2.3 (2.0–2.7) Time between tests (wk) 22.5 (17.2–24.0) Off LT4 On LT4 P TSH (0.4–4.3) (mU/L) 78.0 (48.5–117.5) 0.07 (0.004–0.11) 0.012 Total T3 (1.4–2.5) (nmol/L) 0.72 (0.63–0.83) 2.11 (1.94–2.28) 0.012 Total T4 (58.0–128.0) (nmol/L) 17.0 (14.0–18.0) 155.0 (135.5–170.5) 0.018 Free T4 (11.0–25.0) (pmol/L) 1.55 (1.12–2.17) 26.2 (22.3–29.5) 0.012 Leukocytes (3.5–10.0) (×109/L) 7.6 (5.5–8.1) 7.4 (5.5–8.7) 0.26 Neutrophils (40.0–80.0) (%) 60.0 (52.0–67.0) 62.8 (58.2–65.7) 0.50 Lymphocytes (15.0–50.0) (%) 32.4 (25.1–35.9) 27.3 (23.5–28.6) 0.13 Monocytes (5.0–14.0) (%) 4.7 (4.6–7.8) 8.1 (6.9–9.6) 0.018 Dosage LT4 (µg/kg) 2.3 (2.0–2.7) Time between tests (wk) 22.5 (17.2–24.0) View Large Standard gene expression analysis To assess the effects of thyroid state on gene expression in whole blood cells, we performed next-generation RNA-seq on samples drawn from patients off and on LT4 replacement therapy. At least 15 million reads were generated for each sample, and approximately 95% of these reads were aligned. We first quantified both TR isoforms to confirm that TRα is the main receptor isoform expressed in whole blood samples (Fig. 1A). TRα1 expression appeared reduced during LT4 therapy compared with the hypothyroid state (P = 0.008). TRα1 expression quantified by qPCR analysis also showed a trend toward downregulation, although this failed to reach statistical significance (Supplemental Fig. 1). Reanalysis of the TR isoforms in different blood cell types of independent datasets showed that TRα1 expression was much higher than TRβ1 in platelets, lymphocytes, and neutrophils, suggesting that TRα1 is indeed the major TR isoform expressed in whole blood (Supplemental Fig. 2) (19–21). Figure 1. View largeDownload slide (A) Gene expression (FPKMs) of TRα1 and TRβ1 off and on LT4 treatment. *P < 0.05, **P < 0.01, ***P < 0.001. (B) OmniViz Treescape shows the hierarchical clustering of DE genes off and on LT4 treatment. Red, upregulated genes compared with the geometric mean; blue, downregulated genes compared with the geometric mean. The color intensity correlates with the degree of change. Figure 1. View largeDownload slide (A) Gene expression (FPKMs) of TRα1 and TRβ1 off and on LT4 treatment. *P < 0.05, **P < 0.01, ***P < 0.001. (B) OmniViz Treescape shows the hierarchical clustering of DE genes off and on LT4 treatment. Red, upregulated genes compared with the geometric mean; blue, downregulated genes compared with the geometric mean. The color intensity correlates with the degree of change. After filtering for expressed genes, 16,014 genes remained for analysis. Principal component analysis showed a high degree of clustering of patients off and on LT4 treatment. We detected 1227 DE genes (multiple testing corrected P value < 0.05). We selected 486 DE genes with an absolute fold-change >1.5 (Supplemental Table 3), shown as hierarchical clustering in Fig. 1B. This indicates that thyroid state has both positive and negative effects on gene regulation (76% upregulated; 24% downregulated). The 47 genes that differed at least twofold in expression are presented in Table 2. Table 2. Genes That Differed Twofold in Expression (2Log Ratio >1 or <−1) 2Log Ratio Molecules P Adjusted Description Upregulated genes 1.824 CCNJL 2.0 × 10−24 Cyclin J-like 1.377 SHISA4 5.1 × 10−6 Shisa family member 4 1.342 PLVAP 2.6 × 10−6 Plasmalemma vesicle-associated protein 1.257 PDE5A 9.4 × 10−8 Phosphodiesterase 5A 1.247 RIOK3 5.6 × 10−8 RIO kinase 3 1.247 TUBB1 7.8 × 10−8 Tubulin, beta 1 class VI 1.211 ARHGEF12 1.3 × 10−7 Rho guanine nucleotide exchange factor 12 1.193 SH3BGRL2 5.7 × 10−7 SH3 domain binding glutamate-rich protein like 2 1.179 LINC00989 7.3 × 10−6 Long intergenic non-protein coding RNA 989 1.172 DAAM2 0.00012 Dishevelled associated activator of morphogenesis 2 1.169 CA2 6.2 × 10−6 Carbonic anhydrase II 1.168 CTTN 2.1 × 10−6 Cortactin 1.165 RNF182 0.00027 Ring finger protein 182 1.163 ITGA2B 2.7 × 10−6 Integrin, alpha 2b 1.151 EIF2AK1 4.5 × 10−7 Eukaryotic translation initiation factor 2-alpha kinase 1 1.139 RAB27B 2.1 × 10−6 RAB27B, member RAS oncogene family 1.136 PROS1 4.8 × 10−5 Protein S (alpha) 1.123 CXCR2P1 1.4 × 10−5 Chemokine (C-X-C motif) receptor 2 pseudogene 1 1.114 ABCG2 7.8 × 10−5 ATP-binding cassette, sub-family G (WHITE), member 2 1.111 ITGB3 0.00014 Integrin, beta 3 (platelet glycoprotein IIIa, antigen CD61) 1.093 ARL4A 3.0 × 10−6 ADP-ribosylation factor-like 4A 1.082 STON2 7.8 × 10−8 Stonin 2 1.081 CDKN1A 2.1 × 10−6 Cyclin-dependent kinase inhibitor 1A 1.075 CTNNAL1 0.0009 Catenin (cadherin-associated protein), alpha-like 1 1.073 CTDSPL 6.7 × 10−6 CTD (Carboxy-terminal domain, RNA polymerase II, polypeptide A) small phosphatase-like 1.073 MYLK 1.6 × 10−5 Myosin light chain kinase 1.062 ITGB5 1.2 × 10−5 Integrin, beta 5 1.059 TUBB2A 0.0011 Tubulin, beta 2A class IIa 1.041 MFAP3L 4.4 × 10−6 Microfibrillar-associated protein 3-like 1.031 MYOF 4.2 × 10−6 Myoferlin 1.030 GNAZ 6.2 × 10−6 Guanine nucleotide binding protein (G protein), alpha Z polypeptide 1.030 RNF10 4.3 × 10−6 Ring finger protein 10 1.029 MKRN1 0.0008 Makorin ring finger protein 1 1.026 IFIT1B 0.0018 Interferon-induced protein with tetratricopeptide repeats 1B 1.025 LTBP1 6.5 × 10−5 Latent transforming growth factor beta binding protein 1 1.020 MMD 4.2 × 10−6 Monocyte to macrophage differentiation-associated 1.019 GNG11 4.8 × 10−5 Guanine nucleotide binding protein (G Protein), gamma 11 1.011 ELOVL7 2.9 × 10−5 ELOVL fatty acid elongase 7 1.011 HORMAD1 8.4 × 10−5 HORMA domain containing 1 Downregulated genes −1.014 SNORA80B 1.3 × 10−5 Small nucleolar RNA, H/ACA box 80B −1.064 SCARNA21 1.8 × 10−6 Small Cajal body-specific RNA 21 −1.072 COL4A3 0.00017 Collagen, type IV, alpha 3 −1.091 SNORD10 2.4 × 10−8 Small nucleolar RNA, C/D box 10 −1.140 SNORA34 6.7 × 10−6 Small nucleolar RNA, H/ACA box 34 −1.167 SNORA71D 4.0 × 10−6 Small nucleolar RNA, H/ACA box 71D −1.214 RPS18P9 1.1 × 10−6 Ribosomal protein S18 pseudogene 9 −1.397 SNORA49 6.3 × 10−9 Small nucleolar RNA, H/ACA box 49 2Log Ratio Molecules P Adjusted Description Upregulated genes 1.824 CCNJL 2.0 × 10−24 Cyclin J-like 1.377 SHISA4 5.1 × 10−6 Shisa family member 4 1.342 PLVAP 2.6 × 10−6 Plasmalemma vesicle-associated protein 1.257 PDE5A 9.4 × 10−8 Phosphodiesterase 5A 1.247 RIOK3 5.6 × 10−8 RIO kinase 3 1.247 TUBB1 7.8 × 10−8 Tubulin, beta 1 class VI 1.211 ARHGEF12 1.3 × 10−7 Rho guanine nucleotide exchange factor 12 1.193 SH3BGRL2 5.7 × 10−7 SH3 domain binding glutamate-rich protein like 2 1.179 LINC00989 7.3 × 10−6 Long intergenic non-protein coding RNA 989 1.172 DAAM2 0.00012 Dishevelled associated activator of morphogenesis 2 1.169 CA2 6.2 × 10−6 Carbonic anhydrase II 1.168 CTTN 2.1 × 10−6 Cortactin 1.165 RNF182 0.00027 Ring finger protein 182 1.163 ITGA2B 2.7 × 10−6 Integrin, alpha 2b 1.151 EIF2AK1 4.5 × 10−7 Eukaryotic translation initiation factor 2-alpha kinase 1 1.139 RAB27B 2.1 × 10−6 RAB27B, member RAS oncogene family 1.136 PROS1 4.8 × 10−5 Protein S (alpha) 1.123 CXCR2P1 1.4 × 10−5 Chemokine (C-X-C motif) receptor 2 pseudogene 1 1.114 ABCG2 7.8 × 10−5 ATP-binding cassette, sub-family G (WHITE), member 2 1.111 ITGB3 0.00014 Integrin, beta 3 (platelet glycoprotein IIIa, antigen CD61) 1.093 ARL4A 3.0 × 10−6 ADP-ribosylation factor-like 4A 1.082 STON2 7.8 × 10−8 Stonin 2 1.081 CDKN1A 2.1 × 10−6 Cyclin-dependent kinase inhibitor 1A 1.075 CTNNAL1 0.0009 Catenin (cadherin-associated protein), alpha-like 1 1.073 CTDSPL 6.7 × 10−6 CTD (Carboxy-terminal domain, RNA polymerase II, polypeptide A) small phosphatase-like 1.073 MYLK 1.6 × 10−5 Myosin light chain kinase 1.062 ITGB5 1.2 × 10−5 Integrin, beta 5 1.059 TUBB2A 0.0011 Tubulin, beta 2A class IIa 1.041 MFAP3L 4.4 × 10−6 Microfibrillar-associated protein 3-like 1.031 MYOF 4.2 × 10−6 Myoferlin 1.030 GNAZ 6.2 × 10−6 Guanine nucleotide binding protein (G protein), alpha Z polypeptide 1.030 RNF10 4.3 × 10−6 Ring finger protein 10 1.029 MKRN1 0.0008 Makorin ring finger protein 1 1.026 IFIT1B 0.0018 Interferon-induced protein with tetratricopeptide repeats 1B 1.025 LTBP1 6.5 × 10−5 Latent transforming growth factor beta binding protein 1 1.020 MMD 4.2 × 10−6 Monocyte to macrophage differentiation-associated 1.019 GNG11 4.8 × 10−5 Guanine nucleotide binding protein (G Protein), gamma 11 1.011 ELOVL7 2.9 × 10−5 ELOVL fatty acid elongase 7 1.011 HORMAD1 8.4 × 10−5 HORMA domain containing 1 Downregulated genes −1.014 SNORA80B 1.3 × 10−5 Small nucleolar RNA, H/ACA box 80B −1.064 SCARNA21 1.8 × 10−6 Small Cajal body-specific RNA 21 −1.072 COL4A3 0.00017 Collagen, type IV, alpha 3 −1.091 SNORD10 2.4 × 10−8 Small nucleolar RNA, C/D box 10 −1.140 SNORA34 6.7 × 10−6 Small nucleolar RNA, H/ACA box 34 −1.167 SNORA71D 4.0 × 10−6 Small nucleolar RNA, H/ACA box 71D −1.214 RPS18P9 1.1 × 10−6 Ribosomal protein S18 pseudogene 9 −1.397 SNORA49 6.3 × 10−9 Small nucleolar RNA, H/ACA box 49 View Large Validation of selected genes To validate the RNA-seq data from the discovery cohort, qPCR was used as an independent technique to measure the relative expression of selected DE genes. Genes were selected based on robust baseline expression, robust fold change in expression, and low adjusted P value. All upregulated and most downregulated genes showed substantial differential expression using qPCR (Fig. 2). We validated our findings in an independent cohort; the baseline characteristics are shown in Supplemental Table 4. qPCR analysis of the validation cohort confirmed that all genes were similarly regulated by TH as in the discovery cohort, although this was not substantial for some small nucleolar RNAs (snoRNAs) and phosphodiesterase 5A (Fig. 2). Figure 2. View largeDownload slide Verification of RNA-seq results by qPCR. (A) Six upregulated and (B) six downregulated genes were selected from the RNA-seq results. Results are shown as 2log ratio of the fold change in gene expression off and on LT4 treatment. *P < 0.05, **P < 0.01, ***P < 0.001. Figure 2. View largeDownload slide Verification of RNA-seq results by qPCR. (A) Six upregulated and (B) six downregulated genes were selected from the RNA-seq results. Results are shown as 2log ratio of the fold change in gene expression off and on LT4 treatment. *P < 0.05, **P < 0.01, ***P < 0.001. GO analysis Next, we analyzed whether the 486 DE genes were associated with specific biological processes. GO enrichment analysis revealed that seven biological processes were significantly overrepresented after correction for multiple testing (false discovery rate, 0.05) (Supplemental Table 5). Because some GO terms overlapped based on similar groups of genes, the number of biological processes was reduced to three (Fig. 3). Figure 3. View largeDownload slide GO enrichment analyses of DE genes off and on LT4 treatment. The biological processes are shown on the y-axis. Figure 3. View largeDownload slide GO enrichment analyses of DE genes off and on LT4 treatment. The biological processes are shown on the y-axis. Comparative transcriptome analysis To explore to which extent gene expression in whole blood parallels gene expression in other tissues, we investigated the overlap with gene expression in skeletal muscle in different thyroid states (18). This dataset contained 607 DE genes (fold-change >1.5) in muscle samples off and on LT4 treatment as previously reported. Comparative transcriptome analysis of the 486 DE genes (fold-change >1.5) of the current study with the muscle dataset revealed 26 genes that were shared between both tissues, which is a significant 2.3-fold enrichment (P = 4.1 × 10−5) (Supplemental Table 6). Of these 26 shared genes, 6 were regulated in the opposite direction. If removed from the analysis, this resulted in a 1.7-fold enrichment (P = 0.02). WGCNA To further explore TH-dependent transcriptional patterns in the current dataset, we used WGCNA analysis. WGCNA can reveal the underlying organization of the transcriptome based on coexpression relationships. WGCNA complements traditional DE analyses by providing a system level framework for the understanding of transcriptional profiles. This has been shown particularly helpful in transcriptome analysis of samples composed of distinct cell and tissue types (25). Therefore, WGCNA is potentially useful to identify specific cell types responsive to different thyroid states. First, the genes showing most variability between samples (i.e., at least a 2.0-fold change in level of expression from the global mean) were selected. This resulted in 6649 genes, which were determined by their coefficient of variance, rather than any sample characteristics such as disease status. Subsequently, unsupervised hierarchical clustering led to the identification of 17 coexpression modules (Fig. 4A). Modules correspond to branches and are color-coded, ranging in size from 67 genes in the light cyan module to 1737 in the turquoise module (Supplemental Table 7). Next, modules were identified that were significantly associated with clinical parameters, including age, sex, leukocytes, and thyroid state. Therefore, the summary file (eigengene) for each module was correlated with the clinical parameters to select the most relevant associations (Fig. 4B). The strength of WGCNA analysis is well exemplified by the tan module. The unbiased approach identified 120 genes in this module, which were closely linked according to WGCNA. Subsequent inspection learned that this module mainly contained genes expressed from the Y-chromosome. Regression analysis afterward correctly identified (male) sex to this module (r = 0.97, P = 2.0 × 10−10). Using a similar analysis, the blue and the midnight blue modules were significantly associated with thyroid state. The blue module correlated positively with thyroid state (r = 0.65, P = 0.007) and contained numerous genes (213 of the 814 genes in this module) that were also found significantly regulated by TH in the DE analysis. The positive correlation reflects that most blue module genes were upregulated by TH. Figure 4. View largeDownload slide (A) Network construction identifies distinct modules of coexpressed genes. Modules of coexpressed genes were assigned colors corresponding to the branches indicated by the horizontal bar beneath the dendrogram. (B) Heat map plot of the adjacencies in the eigengene network including several traits. Each row and column in the heat map corresponds to one module eigengene (labeled by color) or trait (labeled on the x-axis). In the heat map, green color represents low adjacency (negative correlation), whereas red represents high adjacency (positive correlation). The blue and the midnight blue modules were significantly associated with thyroid state. Figure 4. View largeDownload slide (A) Network construction identifies distinct modules of coexpressed genes. Modules of coexpressed genes were assigned colors corresponding to the branches indicated by the horizontal bar beneath the dendrogram. (B) Heat map plot of the adjacencies in the eigengene network including several traits. Each row and column in the heat map corresponds to one module eigengene (labeled by color) or trait (labeled on the x-axis). In the heat map, green color represents low adjacency (negative correlation), whereas red represents high adjacency (positive correlation). The blue and the midnight blue modules were significantly associated with thyroid state. Next, we performed a GO enrichment analysis for the genes in the blue module (Supplemental Table 5). The biological process hemostasis was significantly enriched in the blue module. Of note, many genes in the blue module appeared transcripts expressed in platelets (e.g., P2RY12 and PF4). The module midnight blue was negatively correlated with thyroid state (r = −0.78, P = 4.0 × 10−4) and contained mostly genes (32 of the 74 genes) that were regulated by TH as well in the DE analysis, of which 26 (75%) were downregulated. Many genes in the module midnight blue belonged to the class of snoRNAs and contained mostly genes without a GO annotation. Thus, WGNCA analysis suggested that gene expression not only in leukocytes but also in platelets is dependent on thyroid state. Discussion In humans, genes that are regulated by thyroid state are largely unknown. Previously, we and others discovered genes that are dependent on thyroid state in human skeletal muscle (18, 26). Ex vivo studies have identified T3-responsive genes in human skin fibroblasts and adipocytes (26–28). The current study identifies numerous genes in human whole blood samples that are regulated by thyroid state. Similar to previous transcriptome analysis studies, TH largely positively regulates gene expression, although a considerable number of genes is downregulated. Because RNA from leukocytes is a major determinant to total RNA in blood (after removal of globin RNA), it is likely that the detected DE genes largely represent the effects of thyroid state on leukocytes. To improve the yield of our genome-wide expression profiling, we performed WGCNA and compared these results with those of a standard analysis based on differential expression. Although a standard analysis typically discloses lists of DE genes, it fails to recognize the different connections between them. Indeed, WGCNA mapped many genes related to thyroid state into two large coexpression modules (modules blue and midnight blue). GO enrichment analyses of the TH-associated module blue showed predominance of the biological process hemostasis of which the involved genes were upregulated by TH. Of note, many genes were platelet specific transcripts (e.g., P2RY12 and PF4). Because platelets are anucleate, their mRNA content derives from nucleate precursors, which is translated into proteins (29). The product of P2RY12 is a purinergic G-protein coupled receptor, which plays a crucial role in thrombus formation (30). P2RY12 inhibitors, such as clopidogrel and ticagrelor, have antithrombotic effects and are widely used in patients with acute coronary syndromes and in the secondary prevention of thrombotic events in vascular diseases (31). PF4 codes for a chemokine (platelet factor 4), which is synthesized in megakaryocytes and stored in platelet alpha granules. When platelets are activated, PF4 is released from the alpha granules facilitating thrombosis (32). The observation that circulating PF4 levels are decreased in patients with subclinical autoimmune hypothyroidism supports our findings (33). F13A1, another upregulated gene by TH, encodes the coagulation factor XIIIA subunit. Coagulation factor XIII is important for stabilization of the fibrin clot. Pietzner et al. recently demonstrated that thyrotoxicosis increased the levels of several coagulation cascade proteins, including coagulation factors IX, XI, and XIII in plasma (34). Together, these data suggest that TH induces transcription of prothrombotic genes. This is in line with the observation that hyperthyroidism increases the risk of thrombosis (35, 36). Even high-normal thyroid function within the reference range is associated with stroke, independent of classical cardiovascular risk factors (37). Previously, it has been shown that the production of several coagulation factors produced in the liver is enhanced by TH (38, 39). Although we were unable to assess thrombocyte counts, other studies have demonstrated that platelet number is independent of thyroid state (40). The current study suggests that increased TH levels also positively regulate prothrombotic factors in platelets. The second WGCNA module that correlated negatively with thyroid state was mainly composed of snoRNAs that were mostly downregulated by TH. snoRNAs are classified into two families (box C/D and box H/ACA) and are required for posttranscriptional modifications of ribosomal RNA (41). Because the relationship between snoRNAs and TH has not been reported before, further studies are needed to understand the consequences of this finding. Although the effects of thyroid state on peripheral blood cells (leukocytes, platelets) are of interest (see the previous section), we explored to which extent those findings were relevant for other tissues. Because whole blood samples mainly contain TRα-expressing cells, our results putatively reflect the effect of TH on gene expression via TRα. Because TRα is abundantly expressed over TRβ in skeletal muscle, we sought to determine overlap between genes dependent on TH in muscle and genes identified in the current study. The overlap with the DE genes from our previous study in muscle samples was small but statistically significant. This finding suggests that a subset of genes is commonly regulated by TRα1 in various tissues. A few of those genes were regulated in the opposite direction. Although counterintuitive, this phenomenon has been described before (e.g., the well-known T3-responsive gene Reelin) (42). However, the lack of a large overlap also indicates that thyroid state affects different genes in different tissues. Together, these findings may limit the use of whole blood samples as a proxy for the effects of TH on other TRα-expressing tissues in humans (18). Several strengths and limitations are worth mentioning to provide context to our findings. A first limitation is that the observed changes in gene expression may not necessarily reflect direct effects of TH. Gene transcription may also be indirectly dependent on TH if it modulates intermediate signaling molecules. Our studies thus reflect the net effect on gene expression upon changes in thyroid state. In this context, it is worth mentioning the somewhat conflicting results on the limited downregulation of TRα1 in the RNA-seq analysis, which is in line with a previous study that observed a downregulation of TRα1 in hematopoietic progenitor cells obtained from the peripheral blood of hyperthyroid patients compared with a euthyroid control group, although qPCR quantification failed to reach statistical significance (43). These observations suggest that the expression of TRα1 in peripheral blood cells is negatively dependent on thyroid state, possibly dampening the effect size of changes in thyroid state. Second, whole blood contains many different cell types. Although WGCNA was relevant in this context, other subtle changes may have gone unrecognized. Third, DE genes in hypothyroid vs mild thyrotoxic state (as in this study) are not necessarily the same DE genes as in hyperthyroid vs mild thyrotoxic state. Finally, because the patients adhered to an iodine deficient diet at the time of hypothyroidism, direct effects of changes in iodine state on gene expression cannot be excluded. Our study has several strengths. First, the study design included paired analyses, which has the advantage to reduce confounders and variability. Second, the results were confirmed in an independent cohort, substantiating the robustness of our findings. Third, we were able to study extreme differences in thyroid state in human subjects without thyroid autoimmunity. In conclusion, we demonstrate that thyroid state regulates numerous genes in human whole blood. Furthermore, we found that thyroid state affects gene transcripts in platelets, which contributes to the understanding of thrombosis in hyperthyroidism. The overlap with previously reported DE genes in muscle samples indicate that a subset of genes is commonly regulated in TRα-expressing tissues in humans. Future studies should explore if specific transcripts in whole blood can be useful biomarkers for tissue-specific thyroid state. Abbreviations: DE differentially expressed DTC differentiated thyroid cancer FPKM fragments per kilobase of transcript per million mapped reads GO gene ontology LT4 levothyroxine mRNA messenger RNA qPCR quantitative polymerase chain reaction snoRNA small nucleolar RNA TH thyroid hormone TR T3 receptor RNA-seq RNA sequencing TSH thyrotropin WGCNA weighted gene coexpression network analysis. Acknowledgments Financial Support: R.P.P. is supported by Zon-MW TOP Grant 91212044 and a grant from the Erasmus MC Medical research advisory committee (Mrace). G.B. was supported by a grant from the Exchange in Endocrinology Expertise (3E) program. Disclosure Summary: The authors have nothing to disclose. References 1. Yen PM. Physiological and molecular basis of thyroid hormone action. Physiol Rev . 2001; 81( 3): 1097– 1142. Google Scholar CrossRef Search ADS PubMed 2. Yen PM, Ando S, Feng X, Liu Y, Maruvada P, Xia X. Thyroid hormone action at the cellular, genomic and target gene levels. Mol Cell Endocrinol . 2006; 246( 1-2): 121– 127. Google Scholar CrossRef Search ADS PubMed 3. Abel ED, Boers ME, Pazos-Moura C, Moura E, Kaulbach H, Zakaria M, Lowell B, Radovick S, Liberman MC, Wondisford F. Divergent roles for thyroid hormone receptor beta isoforms in the endocrine axis and auditory system. J Clin Invest . 1999; 104( 3): 291– 300. Google Scholar CrossRef Search ADS PubMed 4. Gil-Ibanez P, Garcia-Garcia F, Dopazo J, Bernal J, Morte B. Global transcriptome analysis of primary cerebrocortical cells: identification of genes regulated by triiodothyronine in specific cell types. Cereb Cortex . 2017; 27( 1): 706– 717. Google Scholar PubMed 5. Dong H, Yauk CL, Williams A, Lee A, Douglas GR, Wade MG. Hepatic gene expression changes in hypothyroid juvenile mice: characterization of a novel negative thyroid-responsive element. Endocrinology . 2007; 148( 8): 3932– 3940. Google Scholar CrossRef Search ADS PubMed 6. Herwig A, Campbell G, Mayer CD, Boelen A, Anderson RA, Ross AW, Mercer JG, Barrett P. A thyroid hormone challenge in hypothyroid rats identifies T3 regulated genes in the hypothalamus and in models with altered energy balance and glucose homeostasis. Thyroid . 2014; 24( 11): 1575– 1593. Google Scholar CrossRef Search ADS PubMed 7. Feng X, Jiang Y, Meltzer P, Yen PM. Thyroid hormone regulation of hepatic genes in vivo detected by complementary DNA microarray. Mol Endocrinol . 2000; 14( 7): 947– 955. Google Scholar CrossRef Search ADS PubMed 8. Mohr S, Liew CC. The peripheral-blood transcriptome: new insights into disease and risk assessment. Trends Mol Med . 2007; 13( 10): 422– 432. Google Scholar CrossRef Search ADS PubMed 9. Bochukova E, Schoenmakers N, Agostini M, Schoenmakers E, Rajanayagam O, Keogh JM, Henning E, Reinemund J, Gevers E, Sarri M, Downes K, Offiah A, Albanese A, Halsall D, Schwabe JW, Bain M, Lindley K, Muntoni F, Vargha-Khadem F, Dattani M, Farooqi IS, Gurnell M, Chatterjee K. A mutation in the thyroid hormone receptor alpha gene [published correction appears in N Engl J Med. 2012;367(15):1474]. N Engl J Med . 2012; 366( 3): 243– 249. Google Scholar CrossRef Search ADS PubMed 10. Whitney AR, Diehn M, Popper SJ, Alizadeh AA, Boldrick JC, Relman DA, Brown PO. Individuality and variation in gene expression patterns in human blood. Proc Natl Acad Sci USA . 2003; 100( 4): 1896– 1901. Google Scholar CrossRef Search ADS PubMed 11. Mastrokolias A, den Dunnen JT, van Ommen GB, ’t Hoen PA, van Roon-Mom WM. Increased sensitivity of next generation sequencing-based expression profiling after globin reduction in human blood RNA. BMC Genomics . 2012; 13( 1): 28. Google Scholar CrossRef Search ADS PubMed 12. Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol . 2013; 14( 4): R36. Google Scholar CrossRef Search ADS PubMed 13. Anders S, Pyl PT, Huber W. HTSeq--a Python framework to work with high-throughput sequencing data. Bioinformatics . 2015; 31( 2): 166– 169. Google Scholar CrossRef Search ADS PubMed 14. Trapnell C, Hendrickson DG, Sauvageau M, Goff L, Rinn JL, Pachter L. Differential analysis of gene regulation at transcript resolution with RNA-seq. Nat Biotechnol . 2013; 31( 1): 46– 53. Google Scholar CrossRef Search ADS PubMed 15. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol . 2010; 11( 10): R106. Google Scholar CrossRef Search ADS PubMed 16. Huang W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc . 2009; 4( 1): 44– 57. Google Scholar CrossRef Search ADS PubMed 17. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT, Harris MA, Hill DP, Issel-Tarver L, Kasarskis A, Lewis S, Matese JC, Richardson JE, Ringwald M, Rubin GM, Sherlock G; The Gene Ontology Consortium. Gene ontology: tool for the unification of biology. Nat Genet . 2000; 25( 1): 25– 29. Google Scholar CrossRef Search ADS PubMed 18. Visser WE, Heemstra KA, Swagemakers SM, Ozgür Z, Corssmit EP, Burggraaf J, van Ijcken WF, van der Spek PJ, Smit JW, Visser TJ. Physiological thyroid hormone levels regulate numerous skeletal muscle transcripts. J Clin Endocrinol Metab . 2009; 94( 9): 3487– 3496. Google Scholar CrossRef Search ADS PubMed 19. Londin ER, Hatzimichael E, Loher P, Edelstein L, Shaw C, Delgrosso K, Fortina P, Bray PF, McKenzie SE, Rigoutsos I. The human platelet: strong transcriptome correlations among individuals associate weakly with the platelet proteome. Biol Direct . 2014; 9( 1): 3. Google Scholar CrossRef Search ADS PubMed 20. Bonnal RJ, Ranzani V, Arrigoni A, Curti S, Panzeri I, Gruarin P, Abrignani S, Rossetti G, Pagani M. De novo transcriptome profiling of highly purified human lymphocytes primary cells. Sci Data . 2015; 2: 150051. Google Scholar CrossRef Search ADS PubMed 21. Chatterjee A, Stockwell PA, Rodger EJ, Duncan EJ, Parry MF, Weeks RJ, Morison IM. Genome-wide DNA methylation map of human neutrophils reveals widespread inter-individual epigenetic variation. Sci Rep . 2015; 5( 1): 17328. Google Scholar CrossRef Search ADS PubMed 22. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics . 2008; 9( 1): 559. Google Scholar CrossRef Search ADS PubMed 23. Zhang B, Horvath S. A general framework for weighted gene co-expression network analysis. Stat Appl Genet Mol Biol . 2005; 4( 1): Article17. Google Scholar CrossRef Search ADS PubMed 24. Yip AM, Horvath S. Gene network interconnectedness and the generalized topological overlap measure. BMC Bioinformatics . 2007; 8( 1): 22. Google Scholar CrossRef Search ADS PubMed 25. Oldham MC, Konopka G, Iwamoto K, Langfelder P, Kato T, Horvath S, Geschwind DH. Functional organization of the transcriptome in human brain. Nat Neurosci . 2008; 11( 11): 1271– 1282. Google Scholar CrossRef Search ADS PubMed 26. Clément K, Viguerie N, Diehn M, Alizadeh A, Barbe P, Thalamas C, Storey JD, Brown PO, Barsh GS, Langin D. In vivo regulation of human skeletal muscle gene expression by thyroid hormone. Genome Res . 2002; 12( 2): 281– 291. Google Scholar CrossRef Search ADS PubMed 27. Moeller LC, Dumitrescu AM, Walker RL, Meltzer PS, Refetoff S. Thyroid hormone responsive genes in cultured human fibroblasts. J Clin Endocrinol Metab . 2005; 90( 2): 936– 943. Google Scholar CrossRef Search ADS PubMed 28. Visser WE, Swagemakers SM, Ozgur Z, Schot R, Verheijen FW, van Ijcken WF, van der Spek PJ, Visser TJ. Transcriptional profiling of fibroblasts from patients with mutations in MCT8 and comparative analysis with the human brain transcriptome. Hum Mol Genet . 2010; 19( 21): 4189– 4200. Google Scholar CrossRef Search ADS PubMed 29. Rowley JW, Oler AJ, Tolley ND, Hunter BN, Low EN, Nix DA, Yost CC, Zimmerman GA, Weyrich AS. Genome-wide RNA-seq analysis of human and mouse platelet transcriptomes. Blood . 2011; 118( 14): e101– e111. Google Scholar CrossRef Search ADS PubMed 30. Gachet C. ADP receptors of platelets and their inhibition. Thromb Haemost . 2001; 86( 1): 222– 232. Google Scholar PubMed 31. Storey RF, Sanderson HM, White AE, May JA, Cameron KE, Heptinstall S. The central role of the P(2T) receptor in amplification of human platelet activation, aggregation, secretion and procoagulant activity. Br J Haematol . 2000; 110( 4): 925– 934. Google Scholar CrossRef Search ADS PubMed 32. Thachil J. The prothrombotic potential of platelet factor 4. Eur J Intern Med . 2010; 21( 2): 79– 83. Google Scholar CrossRef Search ADS PubMed 33. Görar S, Ademoğlu E, Çarlıoğlu A, Alioğlu B, Bekdemir H, Sağlam B, Candan Z, Üçler R, Culha C, Aral Y. Low levels of circulating platelet factor 4 (PF4, CXCL4) in subclinically hypothyroid autoimmune thyroiditis. J Endocrinol Invest . 2016; 39( 2): 185– 189. Google Scholar CrossRef Search ADS PubMed 34. Pietzner M, Engelmann B, Kacprowski T, Golchert J, Dirk AL, Hammer E, Iwen KA, Nauck M, Wallaschofski H, Führer D, Münte TF, Friedrich N, Völker U, Homuth G, Brabant G. Plasma proteome and metabolome characterization of an experimental human thyrotoxicosis model. BMC Med . 2017; 15( 1): 6. Google Scholar CrossRef Search ADS PubMed 35. Stuijver DJ, van Zaane B, Romualdi E, Brandjes DP, Gerdes VE, Squizzato A. The effect of hyperthyroidism on procoagulant, anticoagulant and fibrinolytic factors: a systematic review and meta-analysis. Thromb Haemost . 2012; 108( 6): 1077– 1088. Google Scholar CrossRef Search ADS PubMed 36. Debeij J, Dekkers OM, Asvold BO, Christiansen SC, Naess IA, Hammerstrom J, Rosendaal FR, Cannegieter SC. Increased levels of free thyroxine and risk of venous thrombosis in a large population-based prospective study. J Thromb Haemost . 2012; 10( 8): 1539– 1546. Google Scholar CrossRef Search ADS PubMed 37. Chaker L, Baumgartner C, den Elzen WP, Collet TH, Ikram MA, Blum MR, Dehghan A, Drechsler C, Luben RN, Portegies ML, Iervasi G, Medici M, Stott DJ, Dullaart RP, Ford I, Bremner A, Newman AB, Wanner C, Sgarbi JA, Dörr M, Longstreth WT Jr, Psaty BM, Ferrucci L, Maciel RM, Westendorp RG, Jukema JW, Ceresini G, Imaizumi M, Hofman A, Bakker SJ, Franklyn JA, Khaw KT, Bauer DC, Walsh JP, Razvi S, Gussekloo J, Völzke H, Franco OH, Cappola AR, Rodondi N, Peeters RP, Thyroid Studies C; Thyroid Studies Collaboration. Thyroid function within the reference range and the risk of stroke: an individual participant data analysis. J Clin Endocrinol Metab . 2016; 101( 11): 4270– 4282. Google Scholar CrossRef Search ADS PubMed 38. Debeij J, Cannegieter SC, VAN Zaane B, Smit JW, Corssmit EP, Rosendaal FR, Romijn JA, Dekkers OM. The effect of changes in thyroxine and thyroid-stimulating hormone levels on the coagulation system. J Thromb Haemost . 2010; 8( 12): 2823– 2826. Google Scholar CrossRef Search ADS PubMed 39. Horacek J, Maly J, Svilias I, Smolej L, Cepkova J, Vizda J, Sadilek P, Fatorova I, Zak P. Prothrombotic changes due to an increase in thyroid hormone levels. Eur J Endocrinol . 2015; 172( 5): 537– 542. Google Scholar CrossRef Search ADS PubMed 40. Yango J, Alexopoulou O, Eeckhoudt S, Hermans C, Daumerie C. Evaluation of the respective influence of thyroid hormones and TSH on blood coagulation parameters after total thyroidectomy. Eur J Endocrinol . 2011; 164( 4): 599– 603. Google Scholar CrossRef Search ADS PubMed 41. McMahon M, Contreras A, Ruggero D. Small RNAs with big implications: new insights into H/ACA snoRNA function and their role in human disease. Wiley Interdiscip Rev RNA . 2015; 6( 2): 173– 189. Google Scholar CrossRef Search ADS PubMed 42. Manzano J, Morte B, Scanlan TS, Bernal J. Differential effects of triiodothyronine and the thyroid hormone receptor beta-specific agonist GC-1 on thyroid hormone target genes in the b ain. Endocrinology . 2003; 144( 12): 5480– 5487. Google Scholar CrossRef Search ADS PubMed 43. Kawa MP, Grymula K, Paczkowska E, Baskiewicz-Masiuk M, Dabkowska E, Koziolek M, Tarnowski M, Kłos P, Dziedziejko V, Kucia M, Syrenicz A, Machalinski B. Clinical relevance of thyroid dysfunction in human haematopoiesis: biochemical and molecular studies. Eur J Endocrinol . 2010; 162( 2): 295– 305. Google Scholar CrossRef Search ADS PubMed Copyright © 2018 Endocrine Society
Journal of Clinical Endocrinology and Metabolism – Oxford University Press
Published: Jan 1, 2018
It’s your single place to instantly
discover and read the research
that matters to you.
Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.
All for just $49/month
Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly
Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.
All the latest content is available, no embargo periods.
“Whoa! It’s like Spotify but for academic articles.”@Phil_Robichaud