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, T₄ and T₃ and can affect cellular metabolism either through genomic or nongenomic regulation [35]. While T₄ has been shown to have a larger role in nongenomic interactions [36], T₃ 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]. T₃ 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

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 (log₂FC) > 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 log₂ 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

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

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

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

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

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

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

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. Log₂ 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 MgCl₂, 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.