Skip to main content

Six1 promotes skeletal muscle thyroid hormone response through regulation of the MCT10 transporter

Abstract

Background

The Six1 transcription factor is implicated in controlling the development of several tissue types, notably skeletal muscle. Six1 also contributes to muscle metabolism and its activity is associated with the fast-twitch, glycolytic phenotype. Six1 regulates the expression of certain genes of the fast muscle program by directly stimulating their transcription or indirectly acting through a long non-coding RNA. We hypothesized that additional mechanisms of action of Six1 might be at play.

Methods

A combined analysis of gene expression profiling and genome-wide location analysis data was performed. Results were validated using in vivo RNA interference loss-of-function assays followed by measurement of gene expression by RT-PCR and transcriptional reporter assays.

Results

The Slc16a10 gene, encoding the thyroid hormone transmembrane transporter MCT10, was identified as a gene with a transcriptional enhancer directly bound by Six1 and requiring Six1 activity for full expression in adult mouse tibialis anterior, a predominantly fast-twitch muscle. Of the various thyroid hormone transporters, MCT10 mRNA was found to be the most abundant in skeletal muscle, and to have a stronger expression in fast-twitch compared to slow-twitch muscle groups. Loss-of-function of MCT10 in the tibialis anterior recapitulated the effect of Six1 on the expression of fast-twitch muscle genes and led to lower activity of a thyroid hormone receptor-dependent reporter gene.

Conclusions

These results shed light on the molecular mechanisms controlling the tissue expression profile of MCT10 and identify modulation of the thyroid hormone signaling pathway as an additional mechanism by which Six1 influences skeletal muscle metabolism.

Background

The Six family of transcription factors (TFs) are a group of six homeodomain proteins that share a conserved, characteristic Six domain [1, 2]. They were discovered through complementation experiments to the Drosophila sine oculis gene, which is essential for the development of the fruit fly compound eyes [3]. Homologs of the Six family of proteins have been discovered across the animal kingdom [4,5,6] and strong evidence exists that the Six TF-related regulatory network developed evolutionarily prior to the divergence of vertebrates [7]. Six family gene knockout models highlight that these TFs play crucial roles at a global level in the development of several tissue types including skeletal muscle [8], kidney [9, 10], neurones [11], cardiac muscle [12], eyes [13], and gonads [14].

In skeletal muscle, Six1 and Six4 have been shown to be involved in the transcriptional expression of Myogenin [15,16,17], a myogenic regulatory factor (MRF) essential for terminal differentiation of skeletal myofibers [18, 19]. Cooperative interactions between the MRFs and the Six TFs regulate downstream targets of both families of TFs, expressing key genes in the processes of muscle stem cell proliferation, fusion, and differentiation [17, 20,21,22,23]. The role of Six1 in developing muscle is essential, as Six1 knockout mice die at birth due to improperly formed diaphragms [24]. In adult muscle stem cells, Six1 and Six4 are both essential for tissue regeneration following an injury [22, 25]. Comparatively, the role of the Six TFs in adult muscle homeostasis is arguably less well understood. In adult mature skeletal muscle, forced expression of Six1 and its cofactor Eya1 is sufficient to reprogram slow-twitch oxidative muscles toward a fast-twitch glycolytic phenotype [26]. Conversely, myofiber-specific Six1 conditional knockout (Six1-cKO) mice exhibit a switch toward a slow-twitch muscle phenotype, typified by expression of slow-twitch Myosin and Troponin isoforms [27, 28].

The main proposed mechanism for Six1 regulation of adult muscle fiber-type is through direct and positive transcriptional regulation of Linc-MYH, a long non-coding RNA [28]. Linc-MYH knockdown recapitulates a majority of Six1 knockout condition phenotypes, albeit to about 80% of the magnitude of effect, and is incapable of significantly increasing slow-twitch Myosin isoform expression. This partial recapitulation of the Six1 knockout phenotype suggests that additional pathways that regulate fiber-type specification may simultaneously be under Six1 regulatory control. While many regulatory networks are implicated in fiber-type maintenance in adult skeletal muscle, including Calcineurin/NFAT signaling [29], the AMPK axis [30], PGC-1α [31], Sox6 [32], and thyroid hormones (TH) signaling [33], the exact hierarchy and interplay between the different signaling and transcription pathways in the maintenance of adult muscle phenotype is unclear.

Here, we report findings on the implication of Six1 in regulating the TH pathway. Similar to the Six TFs, TH play a crucial role in mammalian health, participating in the development, differentiation, growth, and metabolic homeostasis of various tissues, including the brain, muscles, liver, and pancreas [34]. TH exist in two main forms, T4 and T3 and can affect cellular metabolism either through genomic or nongenomic regulation [35]. While T4 has been shown to have a larger role in nongenomic interactions [36], T3 is the main genomically active form of TH acting as a ligand to activate thyroid receptor TFs (coded by the Thra or Thrb genes) to bind TH response elements (TREs) in the regulatory regions of target genes and modulate their transcription levels [37].

In skeletal muscles, TH activity regulates the processes of myogenesis and muscle regeneration [38, 39], contractile structure [40], energy expenditure [41], mitochondrial thermogenesis, fatty acid oxidation [42], glucose uptake [40], insulin responsiveness [43], and autophagy [44]. In the context of muscle fiber type specification, hyperthyroidism promotes fast-twitch muscle phenotypes, while hypothyroidism is associated with a switch toward slow-twitch phenotype [45, 46]. Mechanistically, TH receptors have been shown to directly control the expression of certain genes associated to the fast-twitch phenotype [47] and to indirectly antagonize the slow-twitch program by upregulating miR-133a1, which in turn downregulates Tead1, a TF driving the slow-twitch program [33, 48]. T3 has been demonstrated to regulate the expression of multiple MRFs, supported by the discovery of functional TREs driving mRNA transcription of Myod1 [49] and Myog [50]. Hyperthyroidism has been shown to push satellite cells toward expressing higher amounts of Myod1 [51], and conversely, hypothyroidism impedes myogenic differentiation [52].

Due to the inability of TH to passively diffuse across lipid bilayers, tissue-specific hormone transporters are required for the cellular efflux and uptake of both T4 and T3 [53]. Several proteins have the ability to transport TH across membranes. They differ in their substrate specificity, tissue expression profile, and whether they operate primarily in hormone influx or efflux [54]. MCT8, MCT10, OATP1C1, LAT1, and LAT2 (encoded by the genes Slc16a2, Slc16a10, Slco1c1, Slc7a5, and Slc7a8, respectively) are the most relevant to TH cellular influx [53]. MCT8 knockout mice present with a regulatory shift away from slow-twitch toward fast-twitch skeletal muscle gene expression [55, 56]. Interestingly, MCT8 is expressed only at low levels in skeletal muscle, and the effect of MCT8 loss-of-function in this tissue are explained by the increased circulating T3 levels in these animals, driving hormone influx and downstream gene activation [55]. MCT8-deficient mice also display impaired muscle regeneration following a myotrauma [57]. MCT8 is detected in muscle satellite cells, the adult stem cells responsible for this process, and specific deletion of MCT8 in this cell type impairs their differentiation, indicating that the regeneration defects in constitutive MCT8 knockout mice are not exclusively due to their above-mentioned serum TH imbalance [57]. MCT10 is a protein structurally related to MCT8 and is expressed at higher levels in skeletal muscle, compared to MCT8 [54]. Yet, less is known about the role MCT10 plays in skeletal muscles. It has been demonstrated recently that in the gastrocnemius muscle, the mRNA levels of MCT10, but not those of MCT8, increase linearly during the aging process of mice, suggesting a potential homeostatic regulatory role for MCT10 [58]. OATP1C1 expression in rodents is mostly restricted to the brain, but its expression has also been detected in muscle satellite cells that are activated by in vitro culture [57]. Finally, LAT1 and LAT2 can also participate in TH uptake [59, 60]. While LAT2 is the most highly expressed of the two in skeletal muscle [53], the mRNA of LAT1 appears more abundant than that of LAT2 in muscle satellite cells [57]. Muscle-specific knockout of LAT1 blunts mTOR-S6K pathway activation by leucine exposure in this tissue [61].

There is paucity of studies on the mechanisms underlying the tissue-specific expression profiles of the various TH transporters [54, 62]. Additionally, to our knowledge, no study has directly investigated links between Six TF and TH functions. The overlapping regulatory roles of both pathways in the regulation of MRF expression [15, 50] and skeletal muscle type specification [28, 33, 63, 64], suggests an interplay between the two regulatory networks worth investigating. With this study, we demonstrate that the thyroid transporter MCT10 is under direct transcriptional control of Six1 in adult skeletal muscle.

Results

Strong correlation between Six1 function and the fast-twitch muscle program

We were initially interested in identifying genes under the direct control of Six1 and which might be associated with the fast-twitch phenotype. For this, we analyzed two gene expression profiling datasets. First, the dataset of Sakakibara et al. provides microarray gene expression profiles in gastrocnemius muscle in wild-type and in muscle-specific Six1 knock-out adult mice (Six1-cKO) [28]. Second, the muscleDB dataset from Terry et al. provides RNA sequencing (RNA-seq) expression profiles in several different murine muscle groups [65]. We started by identifying genes differentially expressed in Six1-cKO gastrocnemius and found 204 genes requiring Six1 for normal expression levels (Fig. 1A; 86 genes less expressed in Six1-cKO and 118 more expressed). To examine in what proportion the genes deregulated in Six1-cKO gastrocnemius are associated with fast- or slow-twitch specific expression programs, we considered their expression in the muscleDB dataset. We limited the analysis to a subset of muscle groups that are predominantly fast-twitch (extensor digitorum longus (EDL), quadriceps, tibialis anterior (TA), gastrocnemius, plantaris, classified in cluster 1 in muscleDB) or slow-twitch (soleus, flexor digitorum brevis (FDB), diaphragm, classified in cluster 2 in muscleDB). We observed a striking association between gene downregulation in the absence of Six1 and genes being normally more highly expressed in fast-twitch muscles, and vice-versa between upregulation in the Six1-cKO and higher gene expression in slow-twitch muscles (Fig. 1B). A Pearson correlation coefficient of −0.72 with two-tailed T test p value smaller than 0.05 was calculated between the log fold-changes (Six1-cKO minus WT and fast minus slow) of these genes (Fig. 1C).

Fig. 1
figure1

Skeletal muscle group distribution and Six1 dependence of gene expression. DNA microarray gene expression profiling from Sakakibara et al. was analyzed to identify differentially expressed genes in the gastrocnemius muscles of wild-type versus Six1-cKO mice. The cut-offs used were abs (log2FC) > 0.58 and Benjamini-Hochberg adjusted p value < 0.05. Genes were annotated based on whether their transcription start sites are within 50 kb of a Six1-binding site in primary myotubes (red bars). Pearson correlation hierarchical clustering was applied to rows (genes) and columns (samples). To further annotate these Six1-dependent genes, their expression from RNA-seq performed in several muscle tissue groups from the muscleDB database (Terry et al.) was also plotted. The rows retained the clustering solution order from the Sakakibara study, and columns were independently clustered by Pearson correlation (clustering performed separately from the Sakakibara data). The log2 fold changes of each gene in the two experiments are shown in the vertical histograms, blue and yellow bars indicating negative and positive values, respectively. A red arrowhead indicates where the Slc16a10 gene appears in these graphs. The gene symbols in this figure are also listed in supplementary Table S7

To increase the likelihood that the genes identified by this analysis are direct transcriptional targets of Six1, the genes were annotated for proximity to a Six1-binding site. Six1-bound genomic loci were identified by chromatin immunoprecipitation followed by high-throughput sequencing (ChIP-seq). Primary myoblasts were isolated and kept in growth phase (MB) or induced to differentiate for 48 h, a duration yielding large and fully differentiated primary myotubes (MT). ChIP-seq was performed with an anti-Six1 antibody or control IgG in MB and MT samples. We identified 11,664 high confidence peaks in MB and 7624 peaks in MT; 3985 regions were bound by Six1 in both MB and MT (Fig. 2A and Tables S1 and S2). Approximately 20% of binding sites were found to be close to the transcription start sites of genes (within 5 kb of the start site); the majority of binding sites however were further from the start sites, suggesting that a sizeable fraction of Six1 target sites represent distal transcriptional enhancers (Fig. 2B and C). In general, the results were highly comparable to results of ChIP-on-chip experiments performed in C2C12 myoblasts [17]; examples of the ChIP-seq read density coverage are provided (Fig. S1A). De novo DNA sequence motif analysis identified the MEF3 sequence element as the most enriched among Six1-binding sites (motif M1, Fig. 2D), a result also found in our earlier study from C2C12 cells [66]. Among 7679 regions bound only in MB, 4311 have at least one occurrence of the MEF3 motif (56%). In sequences bound only in myotubes, the proportion is similar, at 58%. A larger proportion (72%) of regions bound by Six1 in MB and MT has at least one MEF3 site. Other motifs recognized by TFs known to be implicated in myogenesis are among the most enriched, including AP-1 (M3 [67,68,69,70]), the MRFs (E-box element, M2 [71, 72]), Runx1 (M5 [73,74,75]), and Mef2 factors (M6 [76,77,78]). The MEF3 and E-box elements are most enriched among the Six1-binding sites found in both MB and MT while the AP-1 motif is most abundant in Six1-binding sites in MB only or MB and MT. This is corroborated by the observation that anti-c-Jun ChIP-seq read density [73] is highest among the loci bound by Six1 in MB only or MB and MT. Likewise, the signal of anti-MyoD ChIP-seq is highest among the sites bound by Six1 in MB and MT or MT only (Fig. S1B). This proximity on chromatin between Six1 and MyoD highlights the fact that the two factors cooperate in regulating gene expression [17]. The condition-specific Six1-binding sites were mirrored by the profile of H3K4 mono-methylation (H3K4me1) in MB and MT (Fig. S1B): sites bound in MB only are characterized by the highest H3K4me1 ChIP-seq signal in MB, while those bound in MT only have higher H3K4me1 signal in the corresponding sample type. A comparable situation was witnessed with histone H3 lysine 27 acetylation (H3K27ac), where Six1 binding in MB is associated with highest H3K27ac signal in MB, and the highest H3K27ac signal is seen at Six1-binding sites in MT. The Six1-binding sites in cultured MT are relevant to mature skeletal muscle tissue, as they display noticeable H3K4me2 and H3K27ac signal, and open chromatin feature, in gastrocnemius and soleus muscles. The presence of H3K4me1 and H3K4me2 signal at Six1-bound loci reinforces the notion that this TF binds functional distal enhancers. Finally, gene ontology (GO) term enrichment analysis was used to examine the biological processes that the Six1-bound genes participate to (Fig. 2E). As expected, terms related to muscle development and function show a significant enrichment. GO terms associated with several signaling pathways are significantly represented among Six1-bound genes, including the Wnt, BMP, and Hippo pathways and the MAPK cascade. Interestingly, and as noted in our earlier study performed in C2C12 cells, GO terms associated with the development of other tissues known to be affected in Six1 knockout mice, such as the ear, kidney, and thymus, are also significantly enriched [8, 9, 79].

Fig. 2
figure2

Discovery and analysis of Six1-binding sites in primary myoblasts and myotubes. A Venn diagram showing the overlap between Six1-binding peaks in MB and MT. B Proportion of Six1-binding sites, in MB or MT, that overlap various gene features. In this analysis, promoters are defined as the 2 kb region centered at the transcription start site of genes, and the downstream region is defined as the 300 base pairs immediately following the transcription end sites. C Distribution of Six1-binding sites around transcription start sites of genes, in MB and MT. D Results from de novo motif finding using MEME. Redundant motifs and those with lower MEME scores were removed. For each new motif, the most similar motif from the JASPAR database is shown. The enrichment of each motif, in the indicated set of binding site sequences, was calculated. %TP and %FP represent the true positive and false positive percentages, respectively. The Bonferroni-adjusted p value is given under “adj_p-value.” E Biological process GO term enrichment among genes bound by Six1 (all genes within 50 kb of a site bound in MB or bound in MT). Selected terms are shown. N, number of Six1-bound genes with a given annotation; FDR, false discovery rate. The total number of genes assigned to the Six1 peaks is indicated at the top of the respective columns

We annotated the gene expression profiling data with information on whether the genes have a nearby Six1-binding site (Fig. 1A, red bars in the row annotations). Nearly half of the differentially expressed genes in Six1-cKO muscle (100 out of 204) have a Six1-binding site in MT in their vicinity. Interestingly, genes upregulated in Six1-cKO were as likely to be bound by Six1 in MT (57 of 118 genes have a Six1-binding site) as the genes downregulated in Six1-cKO (43 of 86). This analysis reveals that genes likely to be directly activated by Six1 in the gastrocnemius muscle are strongly associated with the fast-twitch phenotype.

Six1 directly regulates expression of the MCT10 transporter in skeletal muscle

We focused our attention on identifying regulators of the fast-twitch program that might be under direct control of Six1. Our attention was drawn to Slc16a10 gene (marked with red arrowhead in Fig. 1B and C), which became the focal point of the present study. Slc16a10 encodes the MCT10 transmembrane transporter for TH, which are linked to fast-twitch skeletal muscle formation [34]. Our analyses found that the expression of MCT10 is downregulated in Six1-cKO gastrocnemius muscle (Fig. 1). Among the five major genes encoding TH influx transporters, MCT10 is expressed in skeletal muscle to considerably higher levels than the other four tested [54]. Further, analysis of muscleDB RNA-seq data showed that MCT10 expression is noticeably higher in the fast-twitch muscle groups compared to the slow ones (Figs. 1 and 3A). The other TH transporters (MCT8, Slco1c1, Slc7a5, and Slc7a8) are expressed at noticeably lower levels and do not show an expression profile aligned with fiber type compositions (Fig. 3B-E). The fiber type specific expression profile of MCT10 was validated using quantitative reverse-transcription PCR (qRT-PCR) on RNA isolated from gastrocnemius, soleus, and TA of mice (Fig. 3F, along with comparison of the profiles of Myh1, Myh2, and Myh7 in Fig. 3G-I as markers of fast type IIX, fast type IIA, and slow type I fibers, respectively). These results, along with prior knowledge on the role of TH in muscle metabolism, suggest that MCT10 may be a mediator of Six1’s effect in establishing the fast-twitch phenotype.

Fig. 3
figure3

Expression of thyroid hormone transporter-coding genes in select murine skeletal muscle groups. RNA-seq expression data from muscleDB were quantitated and normalized by fragments per kilobase of transcript per million mapped reads (FPKM). The expression of the five transporters is shown: A Slc16a10 (MCT10), B Slc16a2 (MCT8), C Slco1c1 (Oatp1c1), D Slc7a5 (Lat1), and E Slc7a8 (Lat2). Each of six replicates is represented by a dot and the statistical distribution of FPKM values is represented by the boxplot. F qRT-PCR validation of Slc16a10 expression in soleus (Sol), gastrocnemius (Gas), and tibialis anterior (TA) muscles. Expression was normalized over the geometric mean of control genes Actnb and 18S rRNA. Expression in soleus and tibialis anterior is lower than in gastrocnemius, by paired, two-tailed t test (p value < 0.05 shown in red). G-I The expression profiles of Myh1, Myh2, and Myh7 are shown for comparison

By examining our ChIP-seq data, we noted the presence of a Six1-binding site located 48 kilobases upstream of the transcription start site of the MCT10 gene (Fig. 4A-B). To obtain additional evidence that this Six1-binding site represents a regulatory element, we analyzed histone mark ChIP-seq and DNA accessibility (assay for transposase-accessible chromatin (ATAC-seq)) data acquired in fast-twitch quadriceps and in slow-twitch soleus [80]. The Six1-binding site is located within a likely gene enhancer, as it is characterized by high DNA accessibility, the presence of H3K3me2 and H3K27ac, and these three marks are all more abundant in the quadriceps samples than the soleus samples (Fig. 4A-B). The Slc16a2/MCT8 locus did not show any Six1-binding site (Fig. 4C). By comparison, Myh1 (fast, type IIX) and Myh4 (fast, type IIB) display higher abundance of these histone marks in the quadriceps, while Myh2 (fast, type IIA) and Myh7 (slow type I) appear more abundant in the soleus (Fig. 4D and E). Binding of Six1 at the Myog locus was also detected, and both histone marks are more abundant at that locus in the quadriceps compared to the soleus (Fig. 4F).

Fig. 4
figure4

Six1-binding profile and epigenetic marks at the genes coding the two main thyroid hormone transporters. At each locus, the gene name, genomic position (mm9 coordinates), and length of the interval shown are given. The Six1 ChIP-seq signal represents data from primary proliferating myoblasts or differentiated myotubes, quantitated in sliding windows across the genome in reads per million sequenced, and subtracted from the signal in input chromatin control sample. The ATAC-seq data in quadriceps and soleus is shown for inferred nucleosome-free regions (insert sites between 10 and 130 base pairs) and quantitated in bins of 10 base pairs in counts per millions sequenced (CPM). The H3K4me2, indicative of enhancers and promoters, and H3K27ac, indicative of active promoters and enhancers, are quantitated in bins of 10 and CPM. A Signal at the Slc16a10 (MCT10) locus. The pink box upstream of the Slc16a10 start site shows the location of a putative Six1-bound enhancer, active predominantly in fast-twitch quadriceps, compared to slow-twitch soleus. B Magnification of the −48 kb enhancer at the Slc16a10 locus. C Signal at the Slc16a2 locus. D Genomic signal at the locus comprising the Myh2 to Myh8 genes. E Genomic signal at the Myh7 locus. F Signal at the Myogenin gene. The blue box represents the area amplified by PCR, in Fig. 5

To validate the Six1-binding site identified by ChIP-seq, we performed a ChIP assay with anti-Six1 and control rabbit IgG on chromatin isolated from pooled mouse hindlimb skeletal muscle tissue. This was followed by quantitative PCR on the ChIP eluates (qChIP) using primers specific for this −48 kb enhancer, the Myog proximal promoter as positive control (blue box in Fig. 4F, [17]), and a negative control locus (Hoxd10 locus, heterochromatic in skeletal muscle [81, 82]). This showed a significant enrichment of Six1 binding (Fig. 5). Our results confirm that the −48kb putative enhancer at the MCT10 locus is a bona fide Six1-binding site.

Fig. 5
figure5

Confirmation of Six1 binding at the Slc16a10 enhancer. Chromatin from hindlimb muscles of four mice was harvested, was fixed, and ChIP was carried out using either anti-Six1 antibody or normal rabbit IgG. Results represent qPCR quantification of IP enrichment (expressed as percentage of input DNA) at the 48 kb binding site upstream of the Slc16a10/MCT10 locus (region corresponding to the pink box in Fig. 4A). Enrichment at the Myog proximal promoter is shown as positive control (region corresponding to the blue box in Fig. 4F) while lack of enrichment at the Hoxd10 proximal promoter is shown as specificity control. Results from each biological replicate are represented by different colors. Two-tailed paired T test were conducted to compare enrichment with IgG and with anti-Six1. One-tailed paired T tests also indicate that the degree of enrichment with anti-Six1 is significantly higher (p < 0.05) at the Myog and Slc16a10 loci compared to the Hoxd10 locus

The microarray transcriptome analysis by Sakakibara et al. (Fig. 1A) were acquired in mice where Six1 was ablated around embryonic day 9, when the HSA-Cre transgene is first detected in the myotome [28, 83]. To explore the effect of Six1 loss-of-function specifically in adult muscles, we employed intra-muscular short interfering RNA (siRNA) duplex electroporation to knockdown Six1 expression in the TA muscle (another predominantly fast-twitch muscle from muscleDB cluster 1). The TA was chosen for its small size, compared to the gastrocnemius, making electroporation more efficient. Total proteins and RNA from electroporated tissues were harvested. A western blot was performed on siRNA-treated muscle extracts demonstrating a corresponding 88% reduction of Six1 protein levels (Fig. 6A and B). The effect on target gene expression was evaluated using qRT-PCR. Six1 mRNA expression levels decreased by approximately 57% under siSix1 knockdown condition; no compensatory effect was seen with Six4, a related family member (Fig. 6C). A significant decrease in MCT10 expression was observed in the Six1-knockdown samples, validating that the expression of this TH transporter depends on Six1 (Fig. 6C). The effects are specific to MCT10, as MCT8 expression was invariant with Six1 knockdown. Together, these results indicate that TH transporter MCT10 is a direct target gene of Six1 and depends on this TF to reach its normal expression level in fast-twitch skeletal muscle.

Fig. 6
figure6

Six1 is needed for MCT10 expression and activity of the TH pathway. A Western blots assaying Six1 protein expression levels under siRNA knockdown condition. siNS represents a control, non-silencing RNA duplex. Note that animals 1 to 4 and 5 to 7 were processed on different days and the western blots were done on different gels and membranes. B Densitometric quantitation of Six1 abundance shown in panel A. Error bars indicate the SEM. The difference in Six1 protein abundance is statistically significant by a one-tailed paired T test relative to siNS (*p value < 0.05). The experiment was performed with seven mice treated identically. C mRNA expression levels of Six1 and MCT10 in electroporated mouse TA muscle after Six1 knockdown, quantified by qRT-PCR. Six4 and MCT8 are shown as members of the same gene families with invariant expression in this experiment. mRNA levels were normalized to the geometric mean of those of 18S rRNA and Actnb. Data represent an average of 7 biological replicates (animals). The normalized expression data is given, with each individual replicate (mouse) shown in a distinct color. Error bars indicate the mean ± SEM of the replicates. P values in red indicate statistically significant decreases in levels of expression as determined by a one-tailed paired Student’s T test relative to siNS (p value < 0.05). D mRNA expression level of a panel of four TH pathway target genes in the control and Six1 knockdown samples from panel C. Two-way ANOVA (testing mRNA levels as a function of which gene was tested and which siRNA was received) indicates significant reduction of the TH signaling gene panel with Six1 knockdown (p value = 0.001)

Thyroid pathway regulation is impacted negatively in Six1 loss-of-function conditions

To determine the extent of the effects of Six1 knockdown and its regulation of MCT10 expression on the thyroid regulatory pathway, a panel of four TH receptor target genes was selected for qRT-PCR analysis: Atp2a1, Slc2a4 (Glut4), Me1, and Ucp3 were chosen as they are among the best characterized direct targets of TH signaling and TH receptors in skeletal muscles [40, 84,85,86,87,88,89]. The downregulation of MCT10 expression with Six1 knockdown was paralleled by that of the four TH pathway genes (Fig. 6D); the gene panel as a whole was significantly less expressed with Six1 loss-of-function by two-way ANOVA. The panel’s reduced expression suggests a decrease in overall thyroid genomic regulatory pathway activity after Six1 knockdown.

We further examined the link between Six1 function and TH signaling in skeletal muscle by performing gene set enrichment analysis (GSEA) on the Six1-cKO gene expression profiling data of Sakakibara et al. For this purpose, and in order to work with gene sets that are most relevant to TH action in skeletal muscle, we generated two gene sets from the study of Nicolaisen et al. [41], where the impact of T3 treatment and muscle-specific knockout of the thyroid hormone receptor alpha gene (Thra-cKO) was examined by RNA-seq. Specifically, we generated a first gene set representing genes whose expression is significantly lower in T3-treated muscle of Thra-cKO compared to wild-types, and a second one with the genes showing the opposite behavior, with expression higher in the Thra-cKO tissue. We found that genes whose expression is lower in Six1-cKO are enriched for genes whose expression is lower in Thra-cKO (Fig. 7A and B). Among the genes in this group, we note Sox6, a repressor of slow-twitch muscle gene transcription [90], and Myoz1 (Calsarcin 2), a protein blocking the slow-twitch program promoting function of the NFAT TF [29, 91, 92]. Four genes involved in cellular respiration were also among this group: Idh3a, Ndufs6, Ndufa5, and Uqcrq. Conversely, we observed that genes upregulated in the Six1-cKO tend to overlap with genes that are upregulated by conditional Thra loss in muscle (Fig. 7C and D). Among these, we noted the calcium homeostasis proteins Atp2a2, Ryr3, and Casq2, the regulators of glucose homeostasis Pdk1 and Pdk3, the slow troponins Tnnc1, Tnnt1, and Tnni1, and genes involved in fatty acid metabolism Abhd1, Acot1, Acsf2, and Tecrl. The results of this analysis show that Six1 activity is directly correlated with that of Thra. Together, these observations indicate that Six1 is functionally and positively associated with TH pathway activity in skeletal muscle.

Fig. 7
figure7

Six1 function is correlated with that of the TH receptor alpha. Gene set enrichment analysis performed on the gene expression profiling of wild-type or Six1-cKO skeletal muscle (Sakakibara et al.), using two custom sets representing genes that are significantly less or more expressed in Thra knock-out skeletal muscle treated with T3 compared to wild-type T3-treated muscle (from Nicolaisen et al.). A Enrichment score graph showing that when genes are ranked from the most downregulated in Six1-cKO to the most upregulated, the beginning of the list is enriched in genes that are more expressed in T3-treated wild-type muscle compared to Thra-cKO muscle. B Heat map of the genes that contribute the most to the enrichment shown in panel A (“core enrichment” genes). The gene order from top to bottom in the heatmap follows the order from left to right in panel A. C and D Similar analyses as in A and B but for the set representing genes that are more expressed in Thra-cKO compared to WT. The gene symbols shown in panels B and D are listed in supplementary files Code_S3 and Code_S4

MCT10 is required for activity of the TH signaling pathway in skeletal muscle

To examine the contribution of MCT10 down-regulation to the TH signaling phenotype observed after Six1 loss-of-function, we used siRNA to knock down the MCT10 transporter directly. Similar to the approach taken with Six1 knockdown, siRNA duplexes targeting the MCT10 transcript (siMCT10) were electroporated into adult TA muscle and mRNA expression levels were assessed using qRT-PCR. siMCT10 electroporation reduced MCT10 mRNA expression levels significantly and had no notable effect on the related MCT8 gene (Fig. 8A). Despite modest downregulation of MCT10 expression after siRNA knockdown (29% reduction) when compared with the efficiency of Six1 knockdown (Fig. 6C), MCT10 knockdown was accompanied by a reduced expression of the TH signaling pathway gene panel (Fig. 8B), with trends similar to those observed under Six1 knockdown conditions (statistical significance as a group by two-way analysis of variance (ANOVA), Fig. 6D). Finally, to further confirm that MCT10 knockdown is associated with a decrease in transcriptional response to TH, we assayed the activity of a TH-responsive luciferase reporter gene electroporated in TA muscles at the same time as the siRNA duplexes. MCT10 knockdown caused a significant reduction in the activity of the TH-dependent reporter gene (Fig. 8C). These results indicate that MCT10 is required for normal TH signaling in skeletal muscle.

Fig. 8
figure8

Recapitulation of TH pathway effects with MCT10 knock-down. A Gene expression in TA muscles electroporated with siRNA against MCT10 or a non-silencing RNA duplex. Expression levels of MCT10, targeted by the siRNA employed, and those of the related gene MCT8 and Six1 are shown. mRNA levels are normalized to the geometric mean of 18S and Actnb. Data were obtained from 4 biological replicates (animals). Error bars indicate the mean ± SEM of the replicates. P values in red indicate statistically significant changes in levels of expression as determined by one-tailed paired Student’s T tests relative to siNS (p value < 0.05). B mRNA expression level of a panel of four TH pathway target genes in the control and MCT10 knockdown samples from panel A. Two-way ANOVA (testing mRNA levels as a function of which gene was tested and which siRNA was received) indicates significant reduction of the TH signaling gene panel with MCT10 knockdown (p value = 0.005). C MCT10 knock-down leads to a decrease in TH-dependent gene transcription. Relative TRE-dependent luciferase activity in TA muscle electroporated with siRNA targeting MCT10. N = 5 biological replicates. The decrease in reporter gene expression is significant as determined by paired one-tailed T test, compared to non-silencing control (siNS)

Conclusion

By analyzing in combination of three genomic studies on gene regulation in skeletal muscle (Six1 ChIP-seq, expression profiling in Six1-cKO and the muscleDB database), we have highlighted the strong association between Six1 function and the fast-twitch phenotype. The results of this analysis not only confirm that Six1 is implicated in establishing the proper gene expression profile of several genes in gastrocnemius muscle [28], but also reaffirm Six1’s function as both activator and repressor of gene expression. Previous reports have shown that Six1 and other Six family TFs may repress transcription in cooperation with co-regulators Sobp [93] and members of the Dach or Groucho/TLE families [15, 94], and activate transcription in association with co-factors of the Eya family or the SWI/SNF complex [15, 79, 95,96,97]. This association between Six1 binding and both gene up- and downregulation after loss-of-function was previously observed in cultured myoblasts and myotubes [17] but the precise mechanisms underlying transcriptional repression by Six1 in skeletal muscle, and how they may be implicated in determining the slow- or fast-twitch expression programs, remain to be uncovered. A limitation of our study is that our Six1 ChIP-seq analysis was performed in cultured primary cells differentiated in vitro, representing a heterogeneous mixture of muscle precursors coming from different muscle groups. It is possible that our Six1 genomic binding data reflects not only maintenance of a differentiated myofiber phenotype such as the fast-twitch program but also the implication of Six1 in upregulating these genes during the myogenic differentiation process. This could in part explain why a binding site for Six1 is found near the slow twitch specific Myh7 gene (Fig. 4E). A genome-wide study of Six1 binding performed in myofibers isolated from specific muscle groups will be required to further assess the function and mode of action of this TF in establishing and maintaining the slow- versus fast-twitch gene expression programs.

Our genomics analysis has revealed that in fast-twitch skeletal muscle, Six1 directly regulates the expression of the TH transporter MCT10. It will be interesting to determine if other tissue types also depend on Six1 for MCT10 expression. In that context, it is worth noting that thyroid gland development is defective in Eya1 knockout mice [98] and that the Six1 and MCT10 mRNAs are both detected in the human adult thyroid [99]. We note that Atp2a1 is a known direct target of Six1 [17, 63], and the analysis of ChIP-seq and gene expression profiling data in Fig. 1 suggests that Me1 represents another direct target of this TF. Whether the downregulation of these genes after Six1 loss-of-function is strictly due to decreased Six1 activity or a combination of direct and indirect effects (e.g., through deregulation of TH signaling) remains to be formally established, but these observations raise the possibility of combinatorial regulation of gene expression by Six1 and TH receptors. Our results also show that continued Six1 expression is required to maintain proper gene expression programs in adult muscle.

We have drawn a parallel between Six1 deficiency and a hypothyroidic condition by showing that Six1 loss-of-function and TH signaling pathway impairment (by Thra gene ablation) share similarities in terms of effects on gene expression. Our findings are consistent with the idea that Six1 deficiency causes a hypothyroidic state by lowering TH entry in skeletal muscle. Conversely, hyperthyroidism and Six1 gain-of-function are both associated with acquisition of a fast-twitch program [26, 45, 46]. However, it remains to be seen to what extent the transcriptomic changes elicited by each condition overlap and if MCT10 induction is part of such a mechanism. The possibility of a crosstalk between TH signaling and Six1-regulated gene program also deserves attention in the future.

As a Six1 target gene, MCT10 was identified as a mediator of the influence of Six1 on the TH signaling pathway in skeletal muscle. The implication of other transporters in TH uptake in skeletal muscle fibers cannot be ruled out. However, MCT10 appears to play a significant role in the process and other transporters are unable to fully compensate the effect of MCT10 knockdown. The situation in mature muscle appears different from that prevailing during regeneration, where MCT8 and OATP1C1 play a more important role [57]. A prior report suggested that MCT10 can provide a compensatory role in the absence of MCT8/Oatp1c1 [57], a condition under which satellite cell-mediated muscle regeneration is temporarily delayed. MCT8/Oatp1c1 knockout mice similarly display impaired neural development, and the two transporters are much more highly expressed in neural stem cells when compared to differentiated lineages [100]. In this context, MCT10 mRNA is much more highly expressed in differentiated lineages than either of the other transporters. This is supported by data at the protein level demonstrating increased MCT10 expression in terminally differentiated neuronal tissue [101]. These data suggest a larger role for MCT10 in adult tissues, with the potential for a compensatory role during early phases of development and differentiation. It has also been reported that the mRNA levels of MCT10, but not those of MCT8, significantly increase in gastrocnemius muscle during the aging process of mice [58]. Thus, the available evidence suggests that MCT10 is involved more closely with adult tissue homeostatic regulation of the TH genetic regulation program. So far, MCT10 knockout has been achieved in mice only in a constitutive fashion; the animals display higher concentration of aromatic amino acids in certain peripheral organs including skeletal muscle [102] but normal circulating TH axis values [103]; the phenotype of MCT10-null skeletal muscle was not thoroughly examined. It will be important to characterize the phenotype of mice with a conditional, muscle-specific deletion of MCT10 to determine its role in the development of this tissue and also in the establishment and maintenance of its metabolic phenotype in the adult.

In addition to TH transporters, other classes of proteins are well-known to participate in TH signaling. The transcriptional effects of TH are mediated by their nuclear receptors, the TFs Thra and Thrb. The intracellular deiodinase enzymes Dio2 and Dio3 convert T4 to the more potent T3, and T3 to the lower-potency T4, RT3, or T2, respectively. Mu-Crystallin, encoded by the Crym gene, is a high-affinity intracellular TH binding protein which may regulate the access of TH to their nuclear receptors and Crym-null mice show hypertrophy of fast type IIb myofibers [104]. In querying the muscleDB mRNA expression dataset, we found that all these proteins are expressed at varying levels in the different muscle groups examined, without any obvious expression bias toward fast- or slow-twitch muscle groups; in fact, Slc16a10/MCT10 remains the gene with the clearest expression bias (Fig. S2). Of these genes, MCT10 was also the only one to be affected significantly by Six1 conditional knockout muscle (Fig. 1).

Considering the importance of TH signaling for muscle function and metabolism, the findings presented here therefore suggest that TH transporter expression must be tightly controlled, and we have shown that Six1 plays an important part in the regulatory mechanism by upregulating MCT10 expression.

Methods

Animals

Female C57BL/6J mice 6–8 weeks old were used for all work and were purchased from Charles River. Mice were sacrificed using cervical dislocation.

ChIP-sequencing

ChIP-seq was performed for Six1 in primary mouse myoblasts as previously described [23]. Sequencing data have been deposited on the NCBI Gene Expression Omnibus (GEO) under accession GSE175999. Primary myoblasts were isolated from 8-week-old female mice and were harvested sub-confluent or grown to confluence and cultured for 48 h in differentiation medium (DMEM supplemented with 10% horse serum) to generate primary myotubes. Preparations were nearly pure, as 95-100% of cells expressed MyoD, as detected by immunofluorescence in growth phase (not shown). Differentiation for 48 h was sufficient to obtain very large myotubes where more than 80% of nuclei were found (comparable to 120 h of differentiation for C2C12 cells, data not shown). To identify Six1-binding sites, MACS2 [105] was used with default parameters in BED mode and using fragment size estimated using SSP [106]. Peaks were refined by splitting large, compound peaks into individual ones using PeakSplitter with a cut-off of 70 reads and a valley of 0.8 [107]. Six1 peaks were further refined by eliminating those that overlap peaks discovered by MACS2 using the input sample, peaks discovered using the IgG against the input control, and mm9 genome blacklisted regions. BED files with peak coordinates are provided in Tables S1 (myoblasts) and S2 (myotubes). Six1 ChIP-seq peaks were annotated to mouse genes using the packages TxDb.Mmusculus.UCSC.mm9.knownGene and ChIPseeker [108], setting a distance cut-off of 50 kb. Gene annotations are provided in Tables S3 and S4. Histone marks ChIP-seq and ATAC-seq in quadriceps and soleus were from the study by Barish et al. [80] and were obtained from the NCBI Short Reads Archive (accessions SRP199043 and SRP173476). Histone marks in C2C12 cells were from Asp et al. [81] and Blum et al. [69] (SRP006743 and SRP012465). TF binding in primary myoblasts was from Umansky et al. [73] (SRP040422). To simplify the analysis, all available biological and technical replicates for each condition were combined into a single file. Reads were trimmed of adapters and low-quality sections using fastp [109], they were aligned to the mouse mm9 genome using STAR [110] in local mode with maximum intron size set to 1 base pair and using a genome index built without gene model GTF file. Only reads aligning to a single genomic location were retained. Duplicates were marked using Picard [111]. Genome-wide sequencing coverage was calculated using deepTools bamCoverage, excluding reads marked as duplicates [112]. In the specific case of ATAC-seq, bamCoverage was restricted to paired-end fragments between 10 and 130 base pairs, assumed to represent nucleosome-free regions. For histone mark ChIP-seq, available as single-end data, bamCoverage was performed using a read extension length calculated using SSP as the most prevalent chromatin fragment size in the respective library. Genomic coverage plots were prepared using R/Bioconductor [113] and the karyoploteR package [114]. ChIP-seq read density (coverage) heatmaps were prepared with deepTools computeMatrix and plotHeatmap.

De novo motif finding was performed within R/Bioconductor using the “memes” package [115] and the STREME program [116]. The sequences under the 2500 tallest peaks (in number of ChIP-seq reads) for the MB-only, MT-only, and the MB and MT sets were combined. The control sequences were the same set, shuffled to retain di-nucleotide frequencies using the “universalmotif” package with the “euler” method, k = 2 and a random seed of 100. Similarity of novel motifs to known motifs from the JASPAR database [117] was identified using Tomtom [118]. Enrichment of the motifs was calculated using Ame [119], using as test set the DNA sequences under the various groups of Six1-binding peaks. The control set of sequences was the union of sequences under all Six1-binding peaks, shuffled with the Euler method, k = 2 and a random seed of 200, to make the control sets for motif discovery and enrichment different.

GO term enrichment analysis was performed in R/Bioconductor using the “topGO” package [120], using the entire set of mouse genes as the background. The test sets were the genes whose transcription start site is closest to each peak, tolerating a maximum distance of 50 kb. Only GO terms with at least 10 mouse genes were included. P values were corrected to the false discovery rate using the Benjamini-Hochberg algorithm [121].

Gene expression profiling data analysis

Gene expression data from the Sakakibara study [28] was retrieved from the NCBI GEO (accession GSE46151). Only gastrocnemius muscle samples were included in the analysis. Data were analyzed using R/Bioconductor and the oligo and limma packages [122, 123], using RMA normalization and a differential expression testing design based on genotype. Log2 fold changes were calculated, and P values were adjusted using the Benjamini-Hochberg algorithm; for significance thresholding, cut-offs of 0.58 and 0.05 were used, respectively. The list of differentially expressed genes is provided in supplementary Table S5. Gene expression data from the Terry et al. muscleDB study [65] was retrieved from the NCBI GEO (accession GSE100505) and SRA accession SRP110541. Reads were trimmed with fastp and aligned to the mm10 mouse genome using STAR in local mode and a genome index built using a gene model GTF file (Mus_musculus.GRCm38.102.gtf obtained from ENSEMBL). Only reads aligning to a single genomic location were retained. Duplicates were marked, but retained for all downstream analyses, using Picard. Read summarization to genes was performed using Subread featureCounts [124]. R/Bioconductor and the DESeq2 package were used to quantitate and normalize gene expression and obtain variance stabilization expression estimates [125]. Data clustering and heatmap generation were accomplished using the pheatmap package [126]. The R script (including package versions and session information) used to analyze data and generate the figures is provided as supplementary files Code_S1. Gene set enrichment analysis was performed using the Sakakibara data and two gene sets from the Nicolaisen et al. study [41]. Briefly, RNA-seq data (read counts over genes) were downloaded from NCBI GEO accession GSE146336 and analyzed in R/Bioconductor using the edgeR package [127]. Data were normalized using TMM and batch effects were modeled using RUVseq with the RUVr algorithm and k = 2 [128]. Genes significantly up- or downregulated in T3-treated Thra-cKO muscle compared to T3-treated wild-type muscle were identified (FDR < 0.05 with the glmQLFTest function) and saved as two separate gene sets. GSEA was performed in R/Bioconductor with the clusterProfiler package [129], using the fgsea algorithm with default parameters, except for the maximum gene set size limit which was adjusted to 1000. The two gene sets (with Entrez Gene ID as identifiers) are provided in supplementary file Code_S2.

siRNA electroporation

Electroporation experiments were performed in the tibialis anterior (TA) muscle on 6–8-week-old female C57BL/6 mice. For knockdown experiments, on day 0 of the experiment, the inter-connective tissue of the TA was digested with an intramuscular injection of 25 μL hyaluronidase (Worthington) at 0.4 unit/μL concentration through the skin of the hindlimb just above the ankle tendon, in order to ensure genetic material can surround and successfully transfect a majority of muscle fibers once injected [130]. One hour later, an intramuscular injection of 22.5 μL siRNA at 20 μM was administered through the skin of the hindlimb into the belly of the TA. Electroporation of nucleic acids was performed with Tweezertrodes (2-paddle electrode assemblies, BTX-Harvard Apparatus 45-0165), used with a 7-mm gap between the electrodes. After siRNA injection, ultrasound gel was applied to the skin, and two electroporation rounds were applied superficially, with the paddles in the first electroporation oriented transverse and sagittal to the hindlimb, and the second electroporation oriented transverse and contra-sagittal to the hindlimb. Electrical stimulation was carried out using a BTX ECM 630 Electroporation system programmed at a setting of 50 volts, 6 pulses, 50 ms/pulse, and 200 ms interval between pulses. Tweezertrodes electroporations were repeated on day 3 of the experiment. Non-silencing control and Six1- or Mct10-targeting siRNA duplexes were obtained from Invitrogen (Stealth siRNA technology) and reconstituted at 20 μM. The sequences are siSix1: GCGAGGAGACCAGCUACUGCUUUAA; siMct10: GCGUCUUCACAAUCCUGCUCCCUUU. The siNS negative control was Stealth RNAi siRNA Negative Control, Med GC sequence 1 (catalog number 12935300). To reduce the effect of biological variability, in all experiments, one leg received the siNS duplex while the other leg of the same animal received the Six1- or Mct10-targeting siRNA. All mice were sacrificed on day 7 of the experiment.

In vivo luciferase assays

For in vivo electroporation experiments that included a transcription reporter plasmid, the siRNA injection on experimental day 0 was performed with a mixture of 22.5 μL siRNA combined with 7.5 μL of the pdV-L1 plasmid at 3 μg/μL in sterile half-saline solution (0.45% NaCl). This plasmid contains two T3 response elements and the SERCA1 basal promoter upstream of the firefly luciferase gene as well as the control renilla luciferase gene preceded by the SERCA1 promoter (kindly provided by Dr. W.S. Simonides (VU Medical Center, Amsterdam, NL) [131, 132]. Experimental day 3 injections contained only siRNA as outlined above. Treated muscles were collected 7 days later, pulverized in liquid nitrogen using a chilled mortar and pestle (Plattner’s) and re-suspended in 100 μL of passive lysis buffer (Promega) per 100 mg of tissue mass. Luciferase results were obtained by reading 5 μL of lysate with 50 μL luciferase assay reagents (Promega) in a Glomax Luminometer according to manufacturer’s specifications. Firefly Luciferase signal readings were normalized to Renilla Luciferase readings to control for variability of electroporation efficiency.

RNA extraction

All samples were extracted in 1 mL of TRIzol Reagent (Invitrogen). Dissected and pulverized TA samples were placed in TRIzol in Lysing Matrix D Tubes (MP Biomedical) and were homogenized in a MagNa Lyser machine (Roche) programed to 7000 RPM for a total of three 20-s bursts separated by 10-s cool-downs on ice. The samples were then spun down for 5 min at 12,000 g at 4 °C to remove fat content. RNA extraction proceeded as recommended by the TRIzol manufacturer. RNA was further purified by treatment with DNase I (RNase-free, New England Biolabs) and heat inactivation of the enzyme as recommended by the manufacturer.

Reverse transcription

Reverse transcription of RNA into cDNA was performed using 500 ng of DNAse-treated total RNA with the SuperScript III First-Strand Synthesis System (Thermo Fisher) using the random hexamers method and post-reaction RNAse H treatment, according to the manufacturer’s recommendation. cDNA samples were diluted to 50 μL with 10 mM Tris pH 8.0 prior to quantitative PCR.

Chromatin immunoprecipitation

For ChIP experiments, each replicate was performed with all hindlimb muscles from two mice. Tissue was dissected and homogenized on ice in hypotonic buffer (25 mM pH 7.8, 1.5 mM MgCl2, 10 mM KCl, 0.1% (v/v) NP-40 and protease inhibitors), using a Teflon Potter-Elvehjem homogenizer mounted on a benchtop drill. Formaldehyde was added from an 11X solution (11% formaldehyde, 0.1 M NaCl, 1.0 mM EDTA, 0.5 M EGTA, 50 mM HEPES pH 8.0) for 10 min at room temperature. Quenching was achieved using 0.125 M glycine for 5 min at room temperature. Samples were spun down at 1000 g for 5 min at 4 °C, resuspended in fresh hypotonic buffer, filtered through 70 μm cell strainer, spun down again at 1000 g for 5 min at 4 °C and re-suspended in 650 μL sonication buffer (10 mM EDTA, 50 mM Tris-HCl pH 8, 0.1% (w/v) SDS, and protease inhibitors). The nuclear pellet was collected by centrifugation and lysed in nuclear lysis buffer (200 mM NaCl, 1 mM EDTA, 0.5 mM EGTA, 10 mM Tris-HCl pH 8, and protease inhibitors). Nuclei were pelleted and resuspended in 650 μL sonication buffer (1 mM EDTA, 0.5 mM EGTA, 0.5% (w/v) sodium sarkosyl, 10 mM Tris-HCl pH 8.0, and protease inhibitors). All sonications were carried out using a probe sonicator at an amplitude of 40% power, with 60 cycles of 1 s on to shear DNA and 4 s off to cool down. The rest of the assay was performed as previously published [17], using 25 μg of chromatin per sample, pre-clearing with Protein A sepharose beads, using 2 μg of antibody per sample with overnight incubation, seven washes of the beads, elution in SDS combined with cross-link reversal and DNA clean-up by phenol extraction followed by ethanol precipitation. Primary antibodies used in ChIP were anti-Six1 (rabbit polyclonal HPA001893, Sigma) and normal rabbit IgG (Jackson ImmunoResearch). The final pellet was re-suspended in 50 μL 10 mM Tris-HCl pH 8.0 for downstream qPCR.

qPCR

qPCR on cDNA and ChIP samples was carried out in reactions containing SYBR Green I HotStarTaq enzyme (Qiagen). Relative quantitation was performed using standard curves made by pooling aliquots of each cDNA sample being examined. Since absolute mRNA abundance of each gene varies in the pool, the expression of a given gene can be compared across samples, but the expression levels of different genes cannot be directly compared. Instead, standard curves for qChIP were made with input genomic DNA; assuming every gene is present in the same allele number in the genome, enrichment levels of different genes can be compared directly. All samples were run as technical triplicates and their geometric mean was taken. For qRT-PCR, relative expression was obtained by dividing by the geometric mean of 18S ribosomal RNA and Actnb, which were found invariable under experimental conditions [133]. For each experiment, the graphs show individual biological replicates (experimental mice) as separate points, along with their mean and standard errors. Oligonucleotide sequences are provided in Table S6.

Protein extraction, quantitation, and immunodetection

Muscle proteins were extracted from the TRIzol organic phase following the procedure recommended by the manufacturer. Protein pellets were resolubilized in Urea-SDS Sample Buffer (6 M Urea, 1% SDS, 20 mM Tris pH 6.8) at a ratio of about 3 μL buffer per 1 mg of tissue and were allowed to solubilize with gentle rocking overnight at 4 °C. Sample protein concentrations were determined using the BCA Protein Assay Kit (Thermo Scientific). Primary antibodies were anti-Six1 (rabbit polyclonal HPA001893, Sigma) and anti-beta-Tubulin (mouse monoclonal hybridoma clone E7, DSHB). Secondary antibodies conjugated with HRP were used. Signal was acquired on X-ray films, which were subsequently digitized and analyzed in the FIJI software for densitometric quantitation.

Statistical analyses

Required number of replicates was determined empirically for sufficient statistical power. No samples were excluded specifically from analysis. R and the rstatix package were used to carry out statistical analyses [134]. Statistical significance for qPCR, densitometry, and luciferase assays were determined by paired T tests, always pairing the experimental muscle electroporated with knockdown siRNA, with its contralateral muscle from the same animal electroporated with siNS. In the presence of prior data suggesting the direction of the effects, one-tailed tests were used. P values ≤ 0.05 were considered to represent means that were significantly different. Where applicable, multiple hypothesis adjustment of p values was performed using the Benjamini-Hochberg method [121]. Unless otherwise stated, no fewer than 3 biological or technical replicates were obtained per experiment. Statistically analyzed control sample groups and experimental sample groups always utilized the same number of technical and biological replicates. In the case of mRNA expression of the TH signaling gene panel (composed of Atp2a1, Me1, Slc2a4, and Ucp3) after Six1 or MCT10 knockdown, a two-way ANOVA test of the effect of independent categorical variables “Gene tested” and “siRNA duplex used” on the dependent continuous variable “Normalized expression” was performed. Normal distribution and variance homogeneity assumptions were confirmed by Shapiro-Wilk’s and Levene’s tests, respectively. P values for the effects of “siRNA duplex used” were significant and are reported.

Availability of data and materials

The ChIP-seq data for the Six1 genomic binding profile has been deposited on the NCBI Gene Expression Omnibus (GEO) under accession GSE175999.

Abbreviations

ANOVA:

Analysis of variance

ATAC-seq:

Assay of transposase-accessible chromatin

ChIP:

Chromatin immunoprecipitation

ChIP-seq:

Chromatin immunoprecipitation followed by high-throughput sequencing

EDL:

Extensor digitorum longus muscle

FDB:

Flexor digitorum brevis muscle

GO:

Gene ontology

GSEA:

Gene set enrichment analysis

H3K4me1:

Monomethylation of histone H3 on lysine 4

H3K4me2:

Dimethylation of histone H3 on lysine 4

H3K27ac:

Acetylation of histone H3 on lysine 27

MB:

Myoblasts

MRF:

Myogenic regulatory factor

MT:

Myotubes

PCR:

Polymerase chain reaction

qChIP:

Chromatin immunoprecipitation followed by quantitative polymerase chain reaction

qRT-PCR:

Reverse-transcription followed by quantitative polymerase chain reaction

RNA-seq:

RNA sequencing

RT-PCR:

Reverse-transcription followed by polymerase chain reaction

siMCT10:

Short-interfering RNA targeting MCT10

siNS:

Non-silencing short-interfering RNA

siRNA:

Short-interfering RNA

siSix1:

Short-interfering RNA targeting Six1

Six1-cKO:

Six1 conditional knockout

TA:

Tibialis anterior muscle

TF:

Transcription factor

TH:

Thyroid hormone

Thra-cKO:

Thyroid hormone receptor alpha gene conditional knockout

TRE:

Thyroid hormone response element

References

  1. 1.

    Kawakami K, Ohto H, Takizawa T, Saito T. Identification and expression of six family genes in mouse retina. FEBS Lett. 1996;393(2–3):259–63.

    CAS  PubMed  Google Scholar 

  2. 2.

    Kawakami K, Sato S, Ozaki H, Ikeda K. Six family genes–structure and function as transcription factors and their roles in development. BioEssays News Rev Mol Cell Dev Biol. 2000;22(7):616–26.

    CAS  Google Scholar 

  3. 3.

    Serikaku MA, O’Tousa JE. sine oculis is a homeobox gene required for Drosophila visual system development. Genetics. 1994;138(4):1137–50.

    CAS  PubMed  PubMed Central  Google Scholar 

  4. 4.

    Pineda D, Gonzalez J, Callaerts P, Ikeo K, Gehring WJ, Salo E. Searching for the prototypic eye genetic network: sine oculis is essential for eye regeneration in planarians. Proc Natl Acad Sci. 2000;97(9):4525–9.

    CAS  PubMed  PubMed Central  Google Scholar 

  5. 5.

    Bebenek IG, Gates RD, Morris J, Hartenstein V, Jacobs DK. sine oculis in basal Metazoa. Dev Genes Evol. 2004;214(7):342–51.

    CAS  PubMed  Google Scholar 

  6. 6.

    Dozier C, Kagoshima H, Niklaus G, Cassata G, Bürglin TR. The Caenorhabditis elegans Six/sine oculis Class Homeobox Gene ceh-32 is required for head morphogenesis. Dev Biol. 2001;236(2):289–303.

    CAS  PubMed  Google Scholar 

  7. 7.

    Abitua PB, Gainous TB, Kaczmarczyk AN, Winchell CJ, Hudson C, Kamata K, et al. The pre-vertebrate origins of neurogenic placodes. Nature. 2015;524(7566):462–5.

    CAS  PubMed  PubMed Central  Google Scholar 

  8. 8.

    Laclef C, Souil E, Demignon J, Maire P. Thymus, kidney and craniofacial abnormalities in Six1 deficient mice. Mech Dev. 2003;120(6):669–79.

    CAS  PubMed  Google Scholar 

  9. 9.

    Xu P-X, Zheng W, Huang L, Maire P, Laclef C, Silvius D. Six1 is required for the early organogenesis of mammalian kidney. Development. 2003;130(14):3085–94.

    CAS  PubMed  Google Scholar 

  10. 10.

    Self M, Lagutin OV, Bowling B, Hendrix J, Cai Y, Dressler GR, et al. Six2 is required for suppression of nephrogenesis and progenitor renewal in the developing kidney. EMBO J. 2006;25(21):5214–28.

    CAS  PubMed  PubMed Central  Google Scholar 

  11. 11.

    Lagutin OV, Zhu CC, Kobayashi D, Topczewski J, Shimamura K, Puelles L, et al. Six3 repression of Wnt signaling in the anterior neuroectoderm is essential for vertebrate forebrain development. Genes Dev. 2003;17(3):368–79.

    CAS  PubMed  PubMed Central  Google Scholar 

  12. 12.

    Boucher CA, King SK, Carey N, Krahe R, Winchester CL, Rahman S, et al. A novel homeodomain-encoding gene is associated with a large CpG island interrupted by the myotonic dystrophy unstable (CTG) n repeat. Hum Mol Genet. 1995;4(10):1919–25.

    CAS  PubMed  Google Scholar 

  13. 13.

    Oliver G, Mailhos A, Wehr R, Copeland NG, Jenkins NA, Gruss P. Six3, a murine homologue of the sine oculis gene, demarcates the most anterior border of the developing neural plate and is expressed during eye development. Development. 1995;121(12):4045–55.

    CAS  PubMed  Google Scholar 

  14. 14.

    Larder R, Clark DD, Miller NLG, Mellon PL. Hypothalamic dysregulation and infertility in mice lacking the homeodomain protein Six6. J Neurosci. 2011;31(2):426–38.

    CAS  PubMed  PubMed Central  Google Scholar 

  15. 15.

    Li X, Oghi KA, Zhang J, Krones A, Bush KT, Glass CK, et al. Eya protein phosphatase activity regulates Six1-Dach-Eya transcriptional effects in mammalian organogenesis. Nature. 2003;426(6964):247–54.

    CAS  PubMed  Google Scholar 

  16. 16.

    Spitz F, Demignon J, Porteu A, Kahn A, Concordet J-P, Daegelen D, et al. Expression of myogenin during embryogenesis is controlled by Six/sine oculis homeoproteins through a conserved MEF3 binding site. Proc Natl Acad Sci. 1998;95(24):14220–5.

    CAS  PubMed  PubMed Central  Google Scholar 

  17. 17.

    Liu Y, Chu A, Chakroun I, Islam U, Blais A. Cooperation between myogenic regulatory factors and SIX family transcription factors is important for myoblast differentiation. Nucleic Acids Res. 2010;38(20):6857–71.

    CAS  PubMed  PubMed Central  Google Scholar 

  18. 18.

    Cheng TC, Tseng BS, Merlie JP, Klein WH, Olson EN. Activation of the myogenin promoter during mouse embryogenesis in the absence of positive autoregulation. Proc Natl Acad Sci U S A. 1995;92(2):561–5.

    CAS  PubMed  PubMed Central  Google Scholar 

  19. 19.

    Molkentin JD, Olson EN. Combinatorial control of muscle development by basic helix-loop-helix and MADS-box transcription factors. Proc Natl Acad Sci. 1996;93(18):9366–73.

    CAS  PubMed  PubMed Central  Google Scholar 

  20. 20.

    Bryson-Richardson RJ, Currie PD. The genetics of vertebrate myogenesis. Nat Rev Genet. 2008;9(8):632–46.

    CAS  PubMed  Google Scholar 

  21. 21.

    Buckingham M, Rigby PW. Gene regulatory networks and transcriptional mechanisms that control myogenesis. Dev Cell. 2014;28(3):225–38.

    CAS  PubMed  Google Scholar 

  22. 22.

    Chakroun I, Yang D, Girgis J, Gunasekharan A, Phenix H, Kærn M, et al. Genome-wide association between Six4, MyoD, and the histone demethylase Utx during myogenesis. FASEB J Off Publ Fed Am Soc Exp Biol. 2015;29(11):4738–55.

    CAS  Google Scholar 

  23. 23.

    Liu Y, Chakroun I, Yang D, Horner E, Liang J, Aziz A, et al. Six1 regulates MyoD expression in adult muscle progenitor cells. PloS One. 2013;8(6):e67762.

    CAS  PubMed  PubMed Central  Google Scholar 

  24. 24.

    Laclef C, Hamard G, Demignon J, Souil E, Houbron C, Maire P. Altered myogenesis in Six1-deficient mice. Dev Camb Engl. 2003;130(10):2239–52.

    CAS  Google Scholar 

  25. 25.

    Le Grand F, Grifone R, Mourikis P, Houbron C, Gigaud C, Pujol J, et al. Six1 regulates stem cell repair potential and self-renewal during skeletal muscle regeneration. J Cell Biol. 2012;198(5):815–32.

    PubMed  PubMed Central  Google Scholar 

  26. 26.

    Grifone R, Laclef C, Spitz F, Lopez S, Demignon J, Guidotti J-E, et al. Six1 and Eya1 expression can reprogram adult muscle from the slow-twitch phenotype into the fast-twitch phenotype. Mol Cell Biol. 2004;24(14):6253–67.

    CAS  PubMed  PubMed Central  Google Scholar 

  27. 27.

    Hetzler KL, Collins BC, Shanely RA, Sue H, Kostek MC. The homoeobox gene SIX1 alters myosin heavy chain isoform expression in mouse skeletal muscle. Acta Physiol Oxf Engl. 2014;210(2):415–28.

    CAS  Google Scholar 

  28. 28.

    Sakakibara I, Santolini M, Ferry A, Hakim V, Maire P. Six homeoproteins and a Iinc-RNA at the fast MYH locus lock fast myofiber terminal phenotype. PLoS Genet. 2014;10(5):e1004386.

    PubMed  PubMed Central  Google Scholar 

  29. 29.

    Chin ER, Olson EN, Richardson JA, Yang Q, Humphries C, Shelton JM, et al. A calcineurin-dependent transcriptional pathway controls skeletal muscle fiber type. Genes Dev. 1998;12(16):2499–509.

    CAS  PubMed  PubMed Central  Google Scholar 

  30. 30.

    Röckl KS, Hirshman MF, Brandauer J, Fujii N, Witters LA, Goodyear LJ. Skeletal muscle adaptation to exercise training. Diabetes. 2007;56(8):2062–9.

    PubMed  Google Scholar 

  31. 31.

    Lin J, Wu H, Tarr PT, Zhang C-Y, Wu Z, Boss O, et al. Transcriptional co-activator PGC-1α drives the formation of slow-twitch muscle fibres. Nature. 2002;418(6899):797–801.

    CAS  PubMed  Google Scholar 

  32. 32.

    Jackson HE, Ono Y, Wang X, Elworthy S, Cunliffe VT, Ingham PW. The role of Sox6 in zebrafish muscle fiber type specification. Skelet Muscle. 2015;5(1):2.

    PubMed  PubMed Central  Google Scholar 

  33. 33.

    Zhang D, Wang X, Li Y, Zhao L, Lu M, Yao X, et al. Thyroid hormone regulates muscle fiber type conversion via miR-133a1. J Cell Biol. 2014;207:jcb-201406068.

    Google Scholar 

  34. 34.

    Zhang J, Lazar MA. The mechanism of action of thyroid hormones. Annu Rev Physiol. 2000;62:439–66.

    CAS  PubMed  Google Scholar 

  35. 35.

    Cheng S-Y, Leonard JL, Davis PJ. Molecular aspects of thyroid hormone actions. Endocr Rev. 2010;31(2):139–70.

    CAS  PubMed  PubMed Central  Google Scholar 

  36. 36.

    Yonkers MA, Ribera AB. Molecular components underlying nongenomic thyroid hormone signaling in embryonic zebrafish neurons. Neural Develop. 2009;4:20.

    Google Scholar 

  37. 37.

    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–55.

    CAS  PubMed  Google Scholar 

  38. 38.

    Dentice M, Marsili A, Ambrosio R, Guardiola O, Sibilio A, Paik J-H, et al. The FoxO3/type 2 deiodinase pathway is required for normal mouse myogenesis and muscle regeneration. J Clin Invest. 2010;120(11):4021.

    CAS  PubMed  PubMed Central  Google Scholar 

  39. 39.

    Dentice M, Ambrosio R, Damiano V, Sibilio A, Luongo C, Guardiola O, et al. Intracellular inactivation of thyroid hormone is a survival mechanism for muscle stem cell proliferation and lineage progression. Cell Metab. 2014;20(6):1038–48.

    CAS  PubMed  PubMed Central  Google Scholar 

  40. 40.

    Simonides WS, van Hardeveld C. Thyroid hormone as a determinant of metabolic and contractile phenotype of skeletal muscle. Thyroid. 2008;18(2):205–16.

    CAS  PubMed  Google Scholar 

  41. 41.

    Nicolaisen TS, Klein AB, Dmytriyeva O, Lund J, Ingerslev LR, Fritzen AM, et al. Thyroid hormone receptor α in skeletal muscle is essential for T3-mediated increase in energy expenditure. FASEB J. 2020;34(11):15480–91.

    CAS  PubMed  Google Scholar 

  42. 42.

    Lombardi A, de Lange P, Silvestri E, Busiello RA, Lanni A, Goglia F, et al. 3,5-Diiodo-l-thyronine rapidly enhances mitochondrial fatty acid oxidation rate and thermogenesis in rat skeletal muscle: AMP-activated protein kinase involvement. Am J Physiol-Endocrinol Metab. 2009;296(3):E497–502.

    CAS  PubMed  Google Scholar 

  43. 43.

    Grozovsky R, Ribich S, Rosene ML, Mulcahey MA, Huang SA, Patti ME, et al. Type 2 deiodinase expression is induced by peroxisomal proliferator-activated receptor-gamma agonists in skeletal myocytes. Endocrinology. 2009;150(4):1976–83.

    CAS  PubMed  Google Scholar 

  44. 44.

    Zhou J, Gauthier K, Ho JP, Lim A, Zhu X-G, Han CR, et al. Thyroid Hormone Receptor α Regulates Autophagy, Mitochondrial Biogenesis, and Fatty Acid Use in Skeletal Muscle. Endocrinology. 2021;162(8):bqab112.

  45. 45.

    Kirschbaum BJ, Kucher HB, Termin A, Kelly AM, Pette D. Antagonistic effects of chronic low frequency stimulation and thyroid hormone on myosin expression in rat fast-twitch muscle. J Biol Chem. 1990;265(23):13974–80.

    CAS  PubMed  Google Scholar 

  46. 46.

    Ianuzzo CD, Hamilton N, Li B. Competitive control of myosin expression: hypertrophy vs. hyperthyroidism. J Appl Physiol. 1991;70(5):2328–30.

    CAS  PubMed  Google Scholar 

  47. 47.

    Bloise FF, Cordeiro A, Ortiga-Carvalho TM. Role of thyroid hormone in skeletal muscle physiology. J Endocrinol. 2018;236(1):R57–68.

    PubMed  Google Scholar 

  48. 48.

    Tsika RW, Schramm C, Simmer G, Fitzsimons DP, Moss RL, Ji J. Overexpression of TEAD-1 in transgenic mouse striated muscles produces a slower skeletal muscle contractile phenotype. J Biol Chem. 2008;283(52):36154–67.

    CAS  PubMed  PubMed Central  Google Scholar 

  49. 49.

    Carnac G, Albagli-Curiel O, Vandromme M, Pinset C, Montarras D, Laudet V, et al. 3, 5, 3’-Triiodothyronine positively regulates both MyoD1 gene transcription and terminal differentiation in C2 myoblasts. Mol Endocrinol. 1992;6(8):1185–94.

    CAS  PubMed  Google Scholar 

  50. 50.

    Downes M, Griggs R, Atkins A, Olson E, Muscat G. Identification of a thyroid hormone response element in the mouse myogenin gene: characterization of the thyroid hormone and retinoid X receptor heterodimeric binding site. Cell Growth Differ. 1993;4:901.

    CAS  PubMed  Google Scholar 

  51. 51.

    Anderson JE, McIntosh LM, Moor AN, Yablonka-Reuveni Z. Levels of MyoD protein expression following injury of mdx and normal limb muscle are modified by thyroid hormone. J Histochem Cytochem Off J Histochem Soc. 1998;46(1):59–67.

    CAS  Google Scholar 

  52. 52.

    Jacobs S, Bär P, Bootsma A. Effect of hypothyroidism on satellite cells and postnatal fiber development in the soleus muscle of rat. Cell Tissue Res. 1996;286(1):137–44.

    CAS  PubMed  Google Scholar 

  53. 53.

    Hennemann G, Docter R, Friesema EC, de Jong M, Krenning EP, Visser TJ. Plasma membrane transport of thyroid hormones and its role in thyroid hormone metabolism and bioavailability. Endocr Rev. 2001;22(4):451–76.

    CAS  PubMed  Google Scholar 

  54. 54.

    Groeneweg S, van Geest FS, Peeters RP, Heuer H, Visser WE. Thyroid hormone transporters. Endocr Rev. 2020;41(2):bnz008.

    PubMed  Google Scholar 

  55. 55.

    Di Cosmo C, Liao X-H, Ye H, Ferrara AM, Weiss RE, Refetoff S, et al. Mct8-deficient mice have increased energy expenditure and reduced fat mass that is abrogated by normalization of serum T3 levels. Endocrinology. 2013;154(12):4885–95.

    PubMed  PubMed Central  Google Scholar 

  56. 56.

    Soukup T, Smerdu V. Effect of altered innervation and thyroid hormones on myosin heavy chain expression and fiber type transitions: a mini-review. Histochem Cell Biol. 2015;143(2):123–30.

    CAS  PubMed  Google Scholar 

  57. 57.

    Mayerl S, Schmidt M, Doycheva D, Darras VM, Hüttner SS, Boelen A, et al. Thyroid hormone transporters MCT8 and OATP1C1 control skeletal muscle regeneration. Stem Cell Rep. 2018;10(6):1959–74.

    CAS  Google Scholar 

  58. 58.

    Wang L, Sheng Y, Xu W, Sun M, Lv S, Yu J, et al. Mechanism of thyroid hormone signaling in skeletal muscle of aging mice. Endocrine. 2021;72(1):132–9.

    CAS  PubMed  Google Scholar 

  59. 59.

    Friesema EC, Docter R, Moerings EP, Verrey F, Krenning EP, Hennemann G, et al. Thyroid hormone transport by the heterodimeric human system L amino acid transporter. Endocrinology. 2001;142(10):4339–48.

    CAS  PubMed  Google Scholar 

  60. 60.

    Ritchie JW, Peter GJ, Shi YB, Taylor PM. Thyroid hormone transport by 4F2hc-IU12 heterodimers expressed in Xenopus oocytes. J Endocrinol. 1999;163(2):R5–9.

    CAS  PubMed  Google Scholar 

  61. 61.

    Poncet N, Mitchell FE, Ibrahim AFM, McGuire VA, English G, Arthur JSC, et al. The catalytic subunit of the system L1 amino acid transporter (slc7a5) facilitates nutrient signalling in mouse skeletal muscle. PloS One. 2014;9(2):e89547.

    PubMed  PubMed Central  Google Scholar 

  62. 62.

    Felmlee MA, Jones RS, Rodriguez-Cruz V, Follman KE, Morris ME. Monocarboxylate transporters (SLC16): function, regulation, and role in health and disease. Pharmacol Rev. 2020;72(2):466–85.

    CAS  PubMed  PubMed Central  Google Scholar 

  63. 63.

    Niro C, Demignon J, Vincent S, Liu Y, Giordani J, Sgarioto N, et al. Six1 and Six4 gene expression is necessary to activate the fast-type muscle gene program in the mouse primary myotome. Dev Biol. 2010;338(2):168–82.

    CAS  PubMed  Google Scholar 

  64. 64.

    Sakakibara I, Wurmser M, Dos Santos M, Santolini M, Ducommun S, Davaze R, et al. Six1 homeoprotein drives myofiber type IIA specialization in soleus muscle. Skelet Muscle. 2016;6(1):30.

    PubMed  PubMed Central  Google Scholar 

  65. 65.

    Terry EE, Zhang X, Hoffmann C, Hughes LD, Lewis SA, Li J, et al. Transcriptional profiling reveals extraordinary diversity among skeletal muscle tissues. eLife. 2018;7:e34613.

  66. 66.

    Liu Y, Nandi S, Martel A, Antoun A, Ioshikhes I, Blais A. Discovery, optimization and validation of an optimal DNA-binding sequence for the Six1 homeodomain transcription factor. Nucleic Acids Res. 2012;40(17):8227–39.

    PubMed  PubMed Central  Google Scholar 

  67. 67.

    Andreucci JJ, Grant D, Cox DM, Tomc LK, Prywes R, Goldhamer DJ, et al. Composition and function of AP-1 transcription complexes during muscle cell differentiation. J Biol Chem. 2002;277(19):16426–32.

    CAS  PubMed  Google Scholar 

  68. 68.

    Tobin SW, Yang D, Girgis J, Farahzad A, Blais A, McDermott JC. Regulation of Hspb7 by MEF2 and AP-1: implications for Hspb7 in muscle atrophy. J Cell Sci. 2016;129(21):4076–90.

    CAS  PubMed  Google Scholar 

  69. 69.

    Blum R, Vethantham V, Bowman C, Rudnicki M, Dynlacht BD. Genome-wide identification of enhancers in skeletal muscle: the role of MyoD1. Genes Dev. 2012;26(24):2763–79.

    CAS  PubMed  PubMed Central  Google Scholar 

  70. 70.

    Bengal E, Ransone L, Scharfmann R, Dwarki VJ, Tapscott SJ, Weintraub H, et al. Functional antagonism between c-Jun and MyoD proteins: a direct physical association. Cell. 1992;68(3):507–19.

    CAS  PubMed  Google Scholar 

  71. 71.

    Blais A, Tsikitis M, Acosta-Alvear D, Sharan R, Kluger Y, Dynlacht BD. An initial blueprint for myogenic differentiation. Genes Dev. 2005;19(5):553–69.

    CAS  PubMed  PubMed Central  Google Scholar 

  72. 72.

    Cao Y, Yao Z, Sarkar D, Lawrence M, Sanchez GJ, Parker MH, et al. Genome-wide MyoD binding in skeletal muscle cells: a potential for broad cellular reprogramming. Dev Cell. 2010;18(4):662–74.

    CAS  PubMed  PubMed Central  Google Scholar 

  73. 73.

    Umansky KB, Feldmesser E, Groner Y. Genomic-wide transcriptional profiling in primary myoblasts reveals Runx1-regulated genes in muscle regeneration. Genomics Data. 2015;6:120–2.

    PubMed  PubMed Central  Google Scholar 

  74. 74.

    Wang X, Blagden C, Fan J, Nowak SJ, Taniuchi I, Littman DR, et al. Runx1 prevents wasting, myofibrillar disorganization, and autophagy of skeletal muscle. Genes Dev. 2005;19(14):1715–22.

    CAS  PubMed  PubMed Central  Google Scholar 

  75. 75.

    Umansky KB, Gruenbaum-Cohen Y, Tsoory M, Feldmesser E, Goldenberg D, Brenner O, et al. Runx1 transcription factor is required for myoblasts proliferation during muscle regeneration. PLoS Genet. 2015;11(8):e1005457.

    PubMed  PubMed Central  Google Scholar 

  76. 76.

    Molkentin JD, Black BL, Martin JF, Olson EN. Cooperative activation of muscle gene expression by MEF2 and myogenic bHLH proteins. Cell. 1995;83(7):1125–36.

    CAS  PubMed  Google Scholar 

  77. 77.

    Sebastian S, Faralli H, Yao Z, Rakopoulos P, Palii C, Cao Y, et al. Tissue-specific splicing of a ubiquitously expressed transcription factor is essential for muscle differentiation. Genes Dev. 2013;27(11):1247–59.

    CAS  PubMed  PubMed Central  Google Scholar 

  78. 78.

    Wales S, Hashemi S, Blais A, McDermott JC. Global MEF2 target gene analysis in cardiac and skeletal muscle reveals novel regulation of DUSP6 by p38MAPK-MEF2 signaling. Nucleic Acids Res. 2014;42(18):11349–62.

    CAS  PubMed  PubMed Central  Google Scholar 

  79. 79.

    Ahmed M, Xu J, Xu P-X. EYA1 and SIX1 drive the neuronal developmental program in cooperation with the SWI/SNF chromatin-remodeling complex and SOX2 in the mammalian inner ear. Dev Camb Engl. 2012;139(11):1965–77.

    CAS  Google Scholar 

  80. 80.

    Ramachandran K, Senagolage MD, Sommars MA, Futtner CR, Omura Y, Allred AL, et al. Dynamic enhancers control skeletal muscle identity and reprogramming. PLoS Biol. 2019;17(10):e3000467.

    CAS  PubMed  PubMed Central  Google Scholar 

  81. 81.

    Asp P, Blum R, Vethantham V, Parisi F, Micsinai M, Cheng J, et al. Genome-wide remodeling of the epigenetic landscape during myogenic differentiation. Proc Natl Acad Sci U S A. 2011;108(22):E149–58.

    PubMed  PubMed Central  Google Scholar 

  82. 82.

    Blais A, van Oevelen CJC, Margueron R, Acosta-Alvear D, Dynlacht BD. Retinoblastoma tumor suppressor protein-dependent methylation of histone H3 lysine 27 is associated with irreversible cell cycle exit. J Cell Biol. 2007;179(7):1399–412.

    CAS  PubMed  PubMed Central  Google Scholar 

  83. 83.

    Miniou P, Tiziano D, Frugier T, Roblot N, Le Meur M, Melki J. Gene targeting restricted to mouse striated muscle lineage. Nucleic Acids Res. 1999;27(19):e27.

    CAS  PubMed  PubMed Central  Google Scholar 

  84. 84.

    Simonides WS, Thelen MH, van der Linden CG, Muller A, van Hardeveld C. Mechanism of thyroid-hormone regulated expression of the SERCA genes in skeletal muscle: implications for thermogenesis. Biosci Rep. 2001;21(2):139–54.

    CAS  PubMed  Google Scholar 

  85. 85.

    Zorzano A, Palacin M, Guma A. Mechanisms regulating GLUT4 glucose transporter expression and glucose transport in skeletal muscle. Acta Physiol. 2005;183(1):43–58.

    CAS  Google Scholar 

  86. 86.

    Dozin B, Magnuson MA, Nikodem VM. Thyroid hormone regulation of malic enzyme synthesis. Dual tissue-specific control. J Biol Chem. 1986;261(22):10290–2.

    CAS  PubMed  Google Scholar 

  87. 87.

    Reitman ML, He Y, Gong DW. Thyroid hormone and other regulators of uncoupling proteins. Int J Obes Relat Metab Disord J Int Assoc Study Obes. 1999;23(Suppl 6):S56–9.

    CAS  Google Scholar 

  88. 88.

    Desvergne B, Petty KJ, Nikodem VM. Functional characterization and receptor binding studies of the malic enzyme thyroid hormone response element. J Biol Chem. 1991;266(2):1008–13.

    CAS  PubMed  Google Scholar 

  89. 89.

    Solanes G, Pedraza N, Calvo V, Vidal-Puig A, Lowell BB, Villarroya F. Thyroid hormones directly activate the expression of the human and mouse uncoupling protein-3 genes through a thyroid response element in the proximal promoter region. Biochem J. 2005;386(3):505–13.

    CAS  PubMed  PubMed Central  Google Scholar 

  90. 90.

    Hagiwara N, Yeh M, Liu A. Sox6 is required for normal fiber type differentiation of fetal skeletal muscle in mice. Dev Dyn Off Publ Am Assoc Anat. 2007;236(8):2062–76.

    CAS  Google Scholar 

  91. 91.

    Frey N, Frank D, Lippl S, Kuhn C, Kögler H, Barrientos T, et al. Calsarcin-2 deficiency increases exercise capacity in mice through calcineurin/NFAT activation. J Clin Invest. 2008;118(11):3598.

    CAS  PubMed  PubMed Central  Google Scholar 

  92. 92.

    Frey N, Richardson JA, Olson EN. Calsarcins, a novel family of sarcomeric calcineurin-binding proteins. Proc Natl Acad Sci. 2000;97(26):14632–7.

    CAS  PubMed  PubMed Central  Google Scholar 

  93. 93.

    Tavares ALP, Jourdeuil K, Neilson KM, Majumdar HD, Moody SA. Sobp modulates the transcriptional activation of Six1 target genes and is required during craniofacial development. Dev Camb Engl. 2021;148(17):dev199684.

  94. 94.

    López-Ríos J, Tessmar K, Loosli F, Wittbrodt J, Bovolenta P. Six3 and Six6 activity is modulated by members of the groucho family. Development. 2003;130(1):185–95.

    PubMed  Google Scholar 

  95. 95.

    Xu J, Li J, Zhang T, Jiang H, Ramakrishnan A, Fritzsch B, et al. Chromatin remodelers and lineage-specific factors interact to target enhancers to establish proneurosensory fate within otic ectoderm. Proc Natl Acad Sci U S A. 2021;118(12):e2025196118.

    CAS  PubMed  PubMed Central  Google Scholar 

  96. 96.

    Ohto H, Kamada S, Tago K, Tominaga S-I, Ozaki H, Sato S, et al. Cooperation of six and eya in activation of their target genes through nuclear translocation of Eya. Mol Cell Biol. 1999;19(10):6815–24.

    CAS  PubMed  PubMed Central  Google Scholar 

  97. 97.

    Ikeda K, Watanabe Y, Ohto H, Kawakami K. Molecular interaction and synergistic activation of a promoter by Six, Eya, and Dach proteins mediated through CREB binding protein. Mol Cell Biol. 2002;22(19):6759–66.

    CAS  PubMed  PubMed Central  Google Scholar 

  98. 98.

    Xu P-X, Zheng W, Laclef C, Maire P, Maas RL, Peters H, et al. Eya1 is required for the morphogenesis of mammalian thymus, parathyroid and thyroid. Development. 2002;129(13):3033–44.

    CAS  PubMed  Google Scholar 

  99. 99.

    GTEx Consortium. The Genotype-Tissue Expression (GTEx) project. Nat Genet. 2013;45(6):580–5.

  100. 100.

    Luongo C, Butruille L, Sébillot A, Le Blay K, Schwaninger M, Heuer H, et al. Absence of both thyroid hormone transporters MCT8 and OATP1C1 impairs neural stem cell fate in the adult mouse subventricular zone. Stem Cell Rep. 2021;16(2):337–53.

    CAS  Google Scholar 

  101. 101.

    Mayerl S, Heuer H, Ffrench-Constant C. Hippocampal neurogenesis requires cell-autonomous thyroid hormone signaling. Stem Cell Rep. 2020;14(5):845–60.

    CAS  Google Scholar 

  102. 102.

    Mariotta L, Ramadan T, Singer D, Guetg A, Herzog B, Stoeger C, et al. T-type amino acid transporter TAT1 (Slc16a10) is essential for extracellular aromatic amino acid homeostasis control. J Physiol. 2012;590(24):6413–24.

    CAS  PubMed  PubMed Central  Google Scholar 

  103. 103.

    Müller J, Mayerl S, Visser TJ, Darras VM, Boelen A, Frappart L, et al. Tissue-specific alterations in thyroid hormone homeostasis in combined Mct10 and Mct8 deficiency. Endocrinology. 2014;155(1):315–25.

    PubMed  Google Scholar 

  104. 104.

    Seko D, Ogawa S, Li T-S, Taimura A, Ono Y. μ-Crystallin controls muscle function through thyroid hormone action. FASEB J Off Publ Fed Am Soc Exp Biol. 2016;30(5):1733–40.

    CAS  Google Scholar 

  105. 105.

    Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 2008;9(9):R137.

    PubMed  PubMed Central  Google Scholar 

  106. 106.

    Nakato R, Shirahige K. Sensitive and robust assessment of ChIP-seq read distribution using a strand-shift profile. Bioinforma Oxf Engl. 2018;34(14):2356–63.

    CAS  Google Scholar 

  107. 107.

    Salmon-Divon M, Dvinge H, Tammoja K, Bertone P. PeakAnalyzer: genome-wide annotation of chromatin binding and modification loci. BMC Bioinformatics. 2010;11:415.

    PubMed  PubMed Central  Google Scholar 

  108. 108.

    Yu G, Wang L-G, He Q-Y. ChIPseeker: an R/Bioconductor package for ChIP peak annotation, comparison and visualization. Bioinforma Oxf Engl. 2015;31(14):2382–3.

    CAS  Google Scholar 

  109. 109.

    Chen S, Zhou Y, Chen Y, Gu J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinforma Oxf Engl. 2018;34(17):i884–90.

    Google Scholar 

  110. 110.

    Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinforma Oxf Engl. 2013;29(1):15–21.

    CAS  Google Scholar 

  111. 111.

    Picard Tools - By Broad Institute [Internet]. [cited 2020 Feb 24]. Available from: http://broadinstitute.github.io/picard/index.html

  112. 112.

    Ramírez F, Ryan DP, Grüning B, Bhardwaj V, Kilpert F, Richter AS, et al. deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Res. 2016;44(W1):W160–5.

    PubMed  PubMed Central  Google Scholar 

  113. 113.

    Huber W, Carey VJ, Gentleman R, Anders S, Carlson M, Carvalho BS, et al. Orchestrating high-throughput genomic analysis with Bioconductor. Nat Methods. 2015;12(2):115–21.

    CAS  PubMed  PubMed Central  Google Scholar 

  114. 114.

    Gel B, Serra E. karyoploteR: an R/Bioconductor package to plot customizable genomes displaying arbitrary data. Bioinforma Oxf Engl. 2017;33(19):3088–90.

    CAS  Google Scholar 

  115. 115.

    Nystrom S. memes: motif matching, comparison, and de novo discovery using the MEME Suite [Internet]. Bioconductor version: Release (3.13); 2021 [cited 2021 Oct 10]. Available from: https://bioconductor.org/packages/memes/

  116. 116.

    Bailey TL. STREME: Accurate and versatile sequence motif discovery. Bioinforma Oxf Engl. 2021;37(18):2834–40.

    Google Scholar 

  117. 117.

    Fornes O, Castro-Mondragon JA, Khan A, van der Lee R, Zhang X, Richmond PA, et al. JASPAR 2020: update of the open-access database of transcription factor binding profiles. Nucleic Acids Res. 2020;48(D1):D87–92.

    CAS  PubMed  Google Scholar 

  118. 118.

    Gupta S, Stamatoyannopoulos JA, Bailey TL, Noble WS. Quantifying similarity between motifs. Genome Biol. 2007;8(2):R24.

    PubMed  PubMed Central  Google Scholar 

  119. 119.

    McLeay RC, Bailey TL. Motif enrichment analysis: a unified framework and an evaluation on ChIP data. BMC Bioinformatics. 2010;11:165.

    PubMed  PubMed Central  Google Scholar 

  120. 120.

    Alexa A, Rahnenfuhrer J. topGO: enrichment analysis for gene ontology [Internet]. Bioconductor version: Release (3.13); 2021 [cited 2021 Oct 10]. Available from: https://bioconductor.org/packages/topGO/

  121. 121.

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

    Google Scholar 

  122. 122.

    Carvalho BS, Irizarry RA. A framework for oligonucleotide microarray preprocessing. Bioinforma Oxf Engl. 2010;26(19):2363–7.

    CAS  Google Scholar 

  123. 123.

    Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47.

    PubMed  PubMed Central  Google Scholar 

  124. 124.

    Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinforma Oxf Engl. 2014;30(7):923–30.

    CAS  Google Scholar 

  125. 125.

    Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

    PubMed  PubMed Central  Google Scholar 

  126. 126.

    CRAN - Package pheatmap [Internet]. [cited 2021 Aug 19]. Available from: https://cran.r-project.org/web/packages/pheatmap/index.html

  127. 127.

    Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinforma Oxf Engl. 2010;26(1):139–40.

    CAS  Google Scholar 

  128. 128.

    Risso D, Ngai J, Speed TP, Dudoit S. Normalization of RNA-seq data using factor analysis of control genes or samples. Nat Biotechnol. 2014;32(9):896–902.

    CAS  PubMed  PubMed Central  Google Scholar 

  129. 129.

    Yu G, Wang L-G, Han Y, He Q-Y. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics J Integr Biol. 2012;16(5):284–7.

    CAS  Google Scholar 

  130. 130.

    McMahon JM, Signori E, Wells KE, Fazio VM, Wells DJ. Optimisation of electrotransfer of plasmid into skeletal muscle by pretreatment with hyaluronidase – increased expression with reduced muscle damage. Gene Ther. 2001;8(16):1264–70.

    CAS  PubMed  Google Scholar 

  131. 131.

    van Mullem AA, van Gucht ALM, Visser WE, Meima ME, Peeters RP, Visser TJ. Effects of thyroid hormone transporters MCT8 and MCT10 on nuclear activity of T3. Mol Cell Endocrinol. 2016;437:252–60.

    PubMed  Google Scholar 

  132. 132.

    Pol CJ, Muller A, Zuidwijk MJ, van Deel ED, Kaptein E, Saba A, et al. Left-ventricular remodeling after myocardial infarction is associated with a cardiomyocyte-specific hypothyroid condition. Endocrinology. 2011;152(2):669–79.

    CAS  PubMed  Google Scholar 

  133. 133.

    Vandesompele J, De Preter K, Pattyn F, Poppe B, Van Roy N, De Paepe A, et al. Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 2002;3(7):RESEARCH0034.

    PubMed  PubMed Central  Google Scholar 

  134. 134.

    Kassambara A. rstatix: pipe-friendly framework for basic statistical tests [Internet]. 2021 [cited 2021 Aug 25]. Available from: https://CRAN.R-project.org/package=rstatix

Download references

Acknowledgements

The authors would like to thank the following individuals or groups for their contributions: Sarah Hemens and Jack Guthrie for experimental assistance; the staff of the animal care and veterinary services of the University of Ottawa for technical assistance; the McGill University and Génome Québec Innovation Centre for ChIP-seq library preparation and high-throughput sequencing; Compute Canada and technical staff for high-performance computing access and support; and Warner Simonides (Amsterdam, The Netherlands) for the TRE reporter plasmid. The E7 hybridoma against beta-tubulin, developed by M. Klymkowsky was obtained from the Developmental Studies Hybridoma Bank, created by the NICHD of the NIH, and maintained at The University of Iowa, Department of Biology, Iowa City, IA 52242. The authors acknowledge the financial support of the Canadian Institutes of Health Research.

Funding

This work was supported by a Canadian Institutes of Health Research operating grant (MOP 119458) to A. B. The funder had no role in the design of the study, data collection, analysis, interpretation of data, or in writing the manuscript.

Author information

Affiliations

Authors

Contributions

AB designed the study, secured funding, performed genomics experiments, analyzed data, and wrote the manuscript. JG performed experiments and analyzed data, DY performed experiments, YL and IC performed experiments and analyzed data. All authors read, edited, and approved the final manuscript.

Corresponding author

Correspondence to Alexandre Blais.

Ethics declarations

Ethics approval and consent to participate

Surgical procedures were performed in compliance with the University of Ottawa Animal Care and Use Committee, the Guidelines of the Canadian Council on Animal Care, and the Animals for Research Act.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Figure S1

. Six1 ChIP-seq results analysis. A) Example of Six1 target loci, showing the Six1 ChIP-seq read coverage in MB and MT. At the Myod1 locus, the core enhancer region (CER) is highlighted in blue and the distal regulatory region (DRR) is highlighted in green. Coverage is shown normalized to counts per million in each library and represented as the difference between the ChIP and input samples, calculated in bins of 10 base pairs. The coverage track for the negative control rabbit IgG ChIP-seq sample is also shown. For each genomic locus, all three tracks are shown at the same scale. B) ChIP-seq coverage at Six1 binding sites, separated in three groups: bound by Six1 in MB only (green), in MB and MT (beige) or MT only (orange). Coverage is shown in a 5 kb window centered on the center of the Six1 binding site. For the same groups of regions, coverage is also given for additional genomics datasets from the indicated studies: Umansky et al. ChIP-seq for c-Jun, MyoD and Runx1 in primary myoblasts (73), Asp et al. and Blum et al. histone mark ChIP-seq in C2C12 MB and MT (69,81), and Barish et al. histone marks and ATAC-seq signal in skeletal muscle groups (gastrocnemius and soleus) (80). For each antibody used (each TF or mark), the heatmap color and average plot scales are maintained so they can be directly compared.

Additional file 2: Figure S2

. Expression of TH pathway genes in different muscle groups. The expression of the indicated genes was estimated from the muscleDB dataset and is represented as RPKM. Data analysis was performed as for Fig. 3.

Additional file 3: Table_S1

.Six1_1MB_peaks(mm9).bed.

Additional file 4: Table_S2

.Six1_1MT_peaks(mm9).bed.

Additional file 5: Table_S3

.Six1_1MB_peaks(mm9).annotated.

Additional file 6: Table_S4

.Six1_1MT_peaks(mm9).annotated.

Additional file 7: Table_S5

.Expression_profiling_cluster_genes.

Additional file 8: Table_S6

.Primers.

Additional file 9: Table_S7

.Expression_profiling_numeric_data.

Additional file 10.

Code_S1.

Additional file 11.

Code_S2.

Additional file 12.

Code_S3.

Additional file 13.

Code_S4.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Girgis, J., Yang, D., Chakroun, I. et al. Six1 promotes skeletal muscle thyroid hormone response through regulation of the MCT10 transporter. Skeletal Muscle 11, 26 (2021). https://doi.org/10.1186/s13395-021-00281-6

Download citation