Increased biological relevance of transcriptome analyses in human skeletal muscle using a model-specific pipeline

Human skeletal muscle responds to weight-bearing exercise with significant inter-individual differences. Investigation of transcriptome responses could improve our understanding of this variation. However, this requires bioinformatic pipelines to be established and evaluated in study-specific contexts. Skeletal muscle subjected to mechanical stress, such as through resistance training (RT), accumulates RNA due to increased ribosomal biogenesis. When a fixed amount of total-RNA is used for RNA-seq library preparations, mRNA counts are thus assessed in different amounts of tissue, potentially invalidating subsequent conclusions. The purpose of this study was to establish a bioinformatic pipeline specific for analysis of RNA-seq data from skeletal muscles, to explore the effects of different normalization strategies and to identify genes responding to RT in a volume-dependent manner (moderate vs. low volume). To this end, we analyzed RNA-seq data derived from a twelve-week RT intervention, wherein 25 participants performed both low- and moderate-volume leg RT, allocated to the two legs in a randomized manner. Bilateral muscle biopsies were sampled from m. vastus lateralis before and after the intervention, as well as before and after the fifth training session (Week 2). Bioinformatic tools were selected based on read quality, observed gene counts, methodological variation between paired observations, and correlations between mRNA abundance and protein expression of myosin heavy chain family proteins. Different normalization strategies were compared to account for global changes in RNA to tissue ratio. After accounting for the amounts of muscle tissue used in library preparation, global mRNA expression increased by 43–53%. At Week 2, this was accompanied by dose-dependent increases for 21 genes in rested-state muscle, most of which were related to the extracellular matrix. In contrast, at Week 12, no readily explainable dose-dependencies were observed. Instead, traditional normalization and non-normalized models resulted in counterintuitive reverse dose-dependency for many genes. Overall, training led to robust transcriptome changes, with the number of differentially expressed genes ranging from 603 to 5110, varying with time point and normalization strategy. Optimized selection of bioinformatic tools increases the biological relevance of transcriptome analyses from resistance-trained skeletal muscle. Moreover, normalization procedures need to account for global changes in rRNA and mRNA abundance.

For transcriptome data to provide adequate biological information about a given experimental set-up, numerous bioinformatic steps need to be adopted in a customized manner [24,27]. Of these steps, data normalization is particularly decisive [26], as it aims to transform naïve transcript counts into biologically meaningful results. This essentially means expressing them as per-cell abundances [28]. For most experimental models, this is equivalent to providing transcript-to-total RNA ratios, given the fulfillment of the assumption that total RNA levels remain stable between conditions on a per-unit-cell or per-unit-tissue basis [28]. In cell models that exhibit high degrees of plasticity, gene expression events result in increased amounts of total RNA and mRNA transcripts per cell [29], specifically violating the assumption that most genes are not differentially expressed [25,28]. We are not aware of any study that has addressed the need to account for such perspectives during transcriptome analyses of skeletal muscle subjected to mechanical stress, such as resistance training. Indeed, this assumption can be expected to be violated, as total RNA content increases markedly on a per-unit-weight basis [3], with potential global changes also occurring for the mRNA pool, though this remains unknown. The extent to which total RNA, and therefore ribosomal RNA, increases, coincides with the increase in muscle mass [3,7], underlining its importance for cellular growth but also its inevitable presence as a potential confounding factor in RNA sequencing experiments.
In this study, we aimed to (1) establish a bioinformatic pipeline specific for analysis of RNA-seq data from skeletal muscles, (2) explore the effects of using different normalization strategies for analyzing skeletal muscle tissue subjected to resistance training, and (3) identify genes responding to moderate versus low resistance exercise volume, while simultaneously identifying genes whose expression changes with time. To achieve these aims, we utilized RNA-seq data generated from a within-participant study, comparing the effects of low and moderate resistance training volume, as previously described [3]. Also, myosin heavy chain protein expression, quantified using immunohistochemistry, was used to validate RNA quantification tools.

Results
For the RNA-seq analyses presented here, a subset of participants was selected based on RNA quality measurements from a previously reported study comparing the functional and biological efficacy of low-and moderate-volume resistance training [3] (Fig. 1a). Twenty-five participants (out of 34) had a full set of RNA-samples with RNA quality indicator (RQI) scores ≥ 7 , which were subjected to bioinformatic data analysis (Fig. 1b). RQI scores were not associated with RNA yield (Fig. 1c). In these participants, twelve weeks of training with moderate volume led to greater increases in limb lean-mass than low volume (3.5% vs. 2.0%, pre-training MOD mean (SD) 8.9 (2.2), to post-training 9.2 (2.3) kg; pre-training LOW 8.9 (2.2), to post-training 9.0 (2.2) kg, Fig. 1d), corresponding well with MRI-based muscle cross-sectional area data from the full data set [3]. Similar responses were seen in the excluded participants (Fig. 1d). This coincided with greater strength gains (~ 25% vs. ~ 19%, Fig. 1f ), which also agrees with observations made in the full cohort, accompanied by greater changes in muscle fiber proportions (type IIX fibers ↓ ) [3]. Biopsies from participants who completed > 85% of prescribed sessions were used for RNA extraction (n = 34; A). RNA quality was assessed (b), and participants with RNA quality indicator (RQI) scores > 7 were submitted for RNA-seq (n = 25). RNA quality was not associated with muscle tissue weight (c), and participants included in RNA-seq experiments did not differ from excluded in terms of limb lean-mass gains (d). Higher training volume led to greater gains in limb lean mass (e) and strength (f) in the lower extremities (n = 25). RNA-seq data were quality filtered using trimgalore and trimmomatic and reads were compared to unfiltered reads (g). Read alignment was performed using five tools of which RSEM, kallisto, and Salmon showed greater fractions of genes with robust expression after removing low-abundance genes (expression filtering; H) compared to HISAT2 and STAR. RSEM, kallisto and Salmon also showed less Log2-differences between biological replicates in a subset of genes with known robust expression (see text for details, i)

Bioinformatic pipeline for analysis of RNA-seq data from skeletal muscles
To select the most appropriate tools for bioinformatic analyses, we first compared quality filtering using Trimmomatic and Trim Galore, both of which are commonly used [30,31]. Quality scores were generally better with Trimmomatic (Fig. 1g). Trim Galore did not improve scores over non-filtered data (Fig. 1g). Filtered reads were then aligned to the human genome using three alignment-based methods (including HISAT2, STAR, RSEM, all used together with Bowtie 2) and two non-alignmentbased methods (kallisto and Salmon). RSEM, Salmon and kallisto all showed similar characteristics in terms of gene counts, resulting in the expected bimodal distribution of counts and a larger subset of detected genes after expression filtering compared to STAR and HISAT2 (Fig. 1h). For a selection of genes with known robust expression across tissues [32], Salmon, kallisto, and RSEM resulted in higher proportions of genes with high count numbers (Fig. 1i). RSEM was found to be associated with lower technical variation than Salmon, kallisto, HISAT2, and STAR expressed as log-differences in the expression of these genes between bilateral biopsies sampled prior to the intervention (Fig. 1j). For HISAT2 and STAR, this distorted correlations between RNA-seq based myosin heavy chain family mRNA and myosin heavy chain protein profiles ( Fig. 2a and b), with the latter identified in Hammarström et al. [3], which are established hallmarks of skeletal muscle biology [33][34][35]. Overall, RSEM, kallisto, and Salmon thus displayed superior technical performance than HISAT2 and STAR, resulting in data with lower degrees of technical variation and higher degrees of biological validity. RSEM displayed slightly lower average variation between paired samples and was thus chosen for downstream analyses. Week 0 Week  Fig. 3b), but after normalization to tissue weight, the two training modalities resulted in similar increases in effective library size (low volume, 43%, and 53% at Week 2 and 12, respectively; moderate volume, 43%, and 49%; Fig. 3c). Overall, this suggests marked increases in global mRNA expression in response to both low-and moderate-volume resistance training.

Identification of genes responding to moderate, compared to low exercise volume
During subsequent identification of differentially expressed (DE) genes in response to low-and moderate-volume resistance training, three normalization models were used and compared. The first model contained effective library sizes as a covariate, as (See figure on next page.) Fig. 3 Global mRNA expression and transcriptome profiles in response to low and moderate volume resistance training. The amounts of muscle tissue used during cDNA synthesis varied over the course of the study and between volume conditions (a low-volume, LOW; moderate-volume, MOD). Library sizes increased during the course of the intervention, with a tendency towards a greater increase in the low-volume condition (b). Difference in library sizes between volume conditions when expressed per-unit tissue weight were diminished, though increases from baseline were maintained (c). The tissue offset-normalized model identified 21 genes with higher expression in the moderate volume condition (d, e), ten of which was shared with the effective library-size normalized model at week 2 (e), and none of which was shared with the naïve model. No volume-dependent differences were found at Week 12 using the tissue-offset model. At this time point, library-size and naïve models both showed a marked skew towards augmented expression in the low-volume condition. At Week 2, functional annotation identified gene sets relating to extracellular matrix in response to higher training volume (tissue-offset model, orange and purple circles, f), all of which were more highly expressed in MOD, indicated by the positive enrichment score. Orange circles denote gene sets that were identified from rank-based enrichment tests based on the full data set. Purple circles denote gene categories that were also identified using over-representation analysis (ORA). Normalization strategies had global effects on enrichment analyses using rank tests, assessed using fold-changes and minimum significant differences scores (not shown), illustrated with the tissue-offset model leading to marked increases in genes associated with the "Collagen containing extracellular matrix" gene set (g) as well as a shift in the full distribution of Log2 fold-changes between volume conditions towards MOD (shown as density curves). Black bars represent genes that belong to the gene set identified as enriched (g). Genes symbols indicate genes identified as differentially expressed in each normalization scenario previously suggested [36], while also containing tissue weight as an offset to account for amounts of tissue used during RNA-seq library preparation (tissue-offset model).
The second model contained effective library sizes as a covariate only (library-size normalization), thus representing an effort to compare expression levels across training modalities while accounting for technical variation during library preparation [25,36].
The third model was a non-normalized model (naïve model, included for comparison) ( Table 1). At Week 2, 21 genes were identified as DE between low and moderate volume using the tissue-offset model, with all genes showing higher expression in the moderate volume condition ( Fig. 3d; rested-state biopsies sampled after four training sessions).   10 and seven genes were identified as DE between volume conditions using the library-size model and the naïve model, respectively ( Fig. 3d and e). For the library-size model, seven of these DE genes showed higher expression in the moderate-volume condition, overlapping completely with genes found using the tissue-offset model (Fig. 3e). For the naïve model, the seven DE genes all showed decreased expression in the moderate volume condition, with 3 genes still overlapping with the library-size model, thus resulting in a contra-intuitive decrease in expression. Using tissue-offset model-derived estimates for functional analyses (Rank-based enrichment tests of minimum significant differences, MSD), revealed enrichment of genes associated with extracellular matrix (ECM) gene ontology (GO) sets (Fig. 3f, Table 2). The top-ranked GO terms were also identified by over-representation tests (ORA) using DE genes (Fig. 3f, Table 2, detailed table in Additional files 1 and 2). Using library-size model estimates similar top-ranked GO sets were identified as with the tissue-offset model albeit with lower levels of significance and lower degrees of agreement between methods ( Table 2). The naïve model generally identified GO sets with negative enrichment scores indicating gene sets with lower expression in moderate volume compared to low volume, with a weak agreement between enrichment methods ( Table 2). The analytical consequences of using the different normalization strategies were particularly apparent in comparisons of rank metrics, such as fold-changes. Importantly, this analytical approach uses the entire gene set to identify enriched gene sets, rather than being confined to DE genes. After controlling for amounts of tissue used during preparation of RNA-seq libraries, the distribution of Log2 differences between volume conditions shifted markedly in favor of higher training volume (Fig. 3g), and robust gene sets appeared with higher expression in the moderate volume condition, such as genes belonging to the Collagen-containing ECM GO set (Fig. 3g). Accordingly, the number of DE genes identified to this GO set was highest using the tissue-offset model (n = 11), followed by the library-size model (n = 6), with no genes identified using the naïve model ( Fig. 3g). At Week 12, no genes were identified as differentially expressed between resistance training with low and moderate volume using the tissue-offset model ( Fig. 3h; restedstate biopsies sampled after finalization of the intervention). In contrast, a small number of genes were identified as DE between volume conditions using library-size and naïve models (n = 4 and n = 3, respectively; Fig. 3h and i), with all genes showing lower    Fig. 3h and i). Of these genes, two were shared between models (Fig. 3i). Using tissue-offset model-derived estimates for functional annotation analyses, revealed no consensus GO sets, using either of the two enrichment approaches (Table 2). In contrast, functional annotation based on estimates from the library-size and naïve model revealed GO terms related to cellular respiration and protein translation with enrichment scores (NES) indicating higher expression of genes in the low-volume condition ( Table 2). No DE genes contributed to ORA-identified GO terms among the top-ranked GO terms from the library-size model. In the naïve model, genes related to cellular respiration were identified as DE with subsequent contribution in ORA-identified GO terms ( Table 2).

Identification of genes with altered expression over time (0, 2 and 12 weeks)
At Week 2 and 12 (rested-state biopsies), we also investigated the overall effects of resistance training on transcriptome profiles: i.e. the time effect, assessed by combining data from the two training modalities. At Week 2, resistance training led to increased expression for 3923, 1609, and 3875 genes and decreased expression of 77, 289, and 100 genes using the tissue-offset, the library-size, and the naïve model, respectively (Fig. 4a). The majority of these DE genes were found in the intersection between all models (Fig. 4a lower panel). At Week 12, resistance training led to increased expression of 1733, 584, and 5108 genes and decreased expression of 2, 19, and 2 genes using the tissue-offset, the library-size, and the naïve model, respectively (Fig. 4b).
Here, the majority number of DE genes were found in the intersection between the tissue-offset model and the naïve model ( Fig. 4b lower panel). At both Week 2 and 12 (and using any normalization model), functional analyses of DE genes revealed enrichment of GO terms associated with ECM structure, organization, and synthesis, as well as stress responses ( Table 3).

Effects of acute exercise on transcriptome profiles (pre-to post-exercise in Week 2)
At Week 2, we also investigated the effects of acute bouts of resistance exercise with low and moderate volume on transcriptome profiles. As we did not expect changes in the total RNA-to-muscle mass ratio in this short time span [37] and rather fluid shifts [38] may have affected tissue weight and hence downstream analyses in an undesirable manner, transcriptome analyses were performed using the library-size model. First, we performed an analysis of the effects of resistance training per se on transcriptome profiles (combining data from the two training modalities). These analyses identified 1736 DE genes after acute resistance exercise, 707 of which showed increased expression and 1029 of which showed decreased expression (Fig. 5a). Genes that showed increased expression were generally associated with stress-related GO terms, including immune response ( Table 4). Genes that showed decreased expression were associated with ECMrelated GO terms (Table 4), contrasting observations made in rested-state biopsies at Week 2 and 12 (Table 3; detailed table in Additional files 1 and 2). We then compared the effects of low and moderate-volume conditions. These analyses identified one single DE gene (RFT1, Fig. 5b), which decreased to a greater extent in the moderate-volume condition. Despite this, rank-based enrichment tests with MSD identified five GO terms with significant enrichment. Among these five categories, three had genes with MSD > 0 (RNA splicing, RNA localization, and Covalent chromatin modification), indicating that the lower bound of 95% CI did not overlap no-change. However, as differences between volume conditions were both negative and positive, as indicated by the rug-plot in Fig. 5c, these categories were not identified in gene-set enrichment analysis based on fold-changes. Overall, these analyses do not provide evidence for pronounced volume-dependent regulation of mRNA expression in the acute recovery phase after resistance exercise (1-h).

Discussion
In the present study, we used a within-participant model to study the effects of low and moderate resistance training volumes on transcriptome responses. For these analyses, we used a subset of muscle biopsy material from a previously reported    study [3]. Training volume led to robust increases in muscle strength and limb lean mass, resembling observations made in the full study cohort [3], and previous studies [39]. Despite these benefits, few differences were detected in transcriptome profiles between volume conditions, with the most prominent exception being a selection of genes involved in extracellular matrix organization and biology in the early stages of resistance training. We have shown that disclosure of these differences was made possible by our systematic selection of bioinformatic tools.

Identification of a model-specific bioinformatics pipeline
The continued development of bioinformatic tools for RNA-seq analyses require continuous optimization of analytic pipelines to any specific study conditions [24,40].
To this end, we first sought to select a suitable read-trimming method. This is necessary to provide high-quality downstream alignment and k-mer search in reads [30]. Two commonly used algorithms were compared [30,31]. Trimmomatic provided data with higher quality than Trimgalore (Fig. 1G), which did not improve quality scores compared to non-filtered data. Second, we compared five mapping tools for performing transcript quantification of trimmed reads: two genome-based mapping tools (STAR [41] and HISAT2 [42]) and three transcript based mapping tools (RSEM [43], kallisto [44] and Salmon [45]). Transcript-based mapping tools resulted in stronger correlations between mRNA and protein profiles, measured as relationships between myosin heavy chain mRNA profiles and protein abundances in rested state biopsies [3] (Fig. 2a), which is known to correlate in resting human skeletal muscle [33][34][35]. This comparison was performed using gene-/protein-family normalization [34,46], allowing the deduction of mRNA-to-protein relationships without the need for other normalization assumptions. Notably, a marked skew was observed in the relationship between MYH1 proportions and its corresponding Type IIX fiber following the initial part of the training intervention. This coincided with robust changes in MYH1 expression, as is typically seen in response to mechanical loading in such short time frames [47]. Furthermore, the transcript-based mapping tools (RSEM, kallisto, and Salmon) resulted in transcriptome profiles with an expected bimodal distribution of counts and a larger subset of detected genes compared to genome-based tools (Fig. 1h). They were also associated with the less technical variation, evident as lower Log2-fold differences in expression for a selection of highly expressed genes between the two legs at baseline [23], assuming minimal biological variation between such paired samples. In these analyses, RSEM displayed slightly lower average variation between paired samples, thus outperforming kallisto and Salmon (Fig. 1j).  (RFT1, b). Three gene ontology categories were identified as significantly enriched using a rank-based test (minimum significant difference, MSD as the rank metric) as well as containing genes with unadjusted P-values < 0.05 (MSD > 0). Traces from the rank test are displayed in c

Comparison of normalization strategies
Transcriptome analyses often rely on the assumption that gene expression is counted and compared between conditions on a per-cell level [28]. This is implicitly assumed to be equivalent to measuring transcriptome data as ratios between mRNA and total RNA, as the input in sequencing or hybridization experiments usually is total RNA [13,15,48].
In the current study population, we previously showed that total RNA increases perunit-muscle tissue in a volume-dependent manner following initiation of resistance training [3]. Consequently, the preparation of cDNA libraries for RNA-seq experiments was unavoidably based on different amounts of tissue, as fixed amounts of total RNA (1000 ng) were used for this purpose. If unaccounted for, this would lead to a comparison of transcriptome data originating from different amounts of muscle tissue for the two volume conditions. We show that this leads to a contra-intuitive larger increase in global transcript counts in the low-volume condition compared to the moderate-volume condition (though without reaching statistical significance). In contrast, after adjusting for the difference in amounts of muscle tissue fed into RNA-seq experiments (i.e., tissueoffset normalization), this apparent difference in average library size disappeared. Overall, we thus observed an increase in global mRNA expression per-unit-muscle weight (43-53%) in response to resistance training that did not depend on training volume, contrasting the observed volume-dependency of total RNA expression [3]. This global change in mRNA expression was associated with substantial alterations in the expression of a multitude of genes, with as many as 26 and 12% of the total read-count pool showing increased expression using the tissue-offset model at Week 2 and Week 12, respectively. These genes were associated with biological processes such as ECM synthesis and organization corroborating with previous studies [13,18,22]. As the volume-dependent changes in muscle growth and total RNA levels in m. vastus lateralis arguably will affect downstream bioinformatics analyses and identification of DE genes, we aimed to compare three different normalization strategies during the subsequent analyses: tissue-offset, library-size, and naïve. At Weeks 2 and 12, the different normalization strategies resulted in marked shifts in global mRNA responses between the two volume conditions, having pronounced effects on identification of GO terms during enrichment analyses (Fig. 3g). In general, tissue-offset normalization (providing mRNA expression per-mg-muscle weight) resulted in a global shift in mRNA responses towards moderate volume. At Week 2, this was evident as more robust increases in mRNA expression in response to higher training volume for most genes, contrasting findings in library-size and naïve analyses. At Week 12, this was evident as a counterbalancing of mRNA expression profiles, with tissue-offset providing rather normally distributed responses to the two volume conditions, contrasting the skew towards larger mRNA expression in response to low training volume in library-size and naïve analyses. The utilization of generalized linear mixed models (GLMM) allowed convenient comparisons of normalization models, which could then be fitted using in the same statistical framework, as previously suggested [36]. GLMM also allowed the incorporation of random effects into the model to account for the repeated measures design. Although there are approaches to account for correlated observations in commonly used RNA-seq modeling frameworks [49], GLMMs provides a more robust and potentially more powerful framework for dealing with correlated data [36].

Training volume-dependent changes in transcriptome profiles.
Using the tissue-offset model, we were able to identify genes relating to ECM functions as volume-sensitive during the early stages of resistance training. This may indicate a role for ECM remodeling in the beneficial effects of higher training volumes on muscular adaptations and strength. As such, previous research has shown that ECM remodeling is induced by exercise training, both acutely [50][51][52][53], and after prolonged endurance and resistance training [13,18,54,55], at both mRNA and protein levels [13,18,54,55]. However, none of these previous studies have found ECM remodeling to be differentially affected by different exercise modalities. Rather, different resistance training modalities such as low-and high-load training have been associated with similar responses, measured as collagen synthesis [52]. Importantly, there seems to be a close association between training-induced changes in abundances of ECM-related mRNAs and their respective proteins, including collagen-organization proteins [54]. This suggests that ECM remodeling is primarily controlled at the transcriptional level, which arguably increases the biological relevance of the herein presented transcriptome analyses. However, this relationship seems to involve a complex time course dependency. For example, transcriptional regulation of COL1A2 shows a considerable lag from stimuli to transcription, as shown in fibroblasts [56]. In line with this, our data on the effects of acute resistance exercise on transcriptome profiles, suggests a counterintuitive reduction in the expression of e.g., collagen mRNA immediately after exercise, as has also been found by others [18]. Indeed, enrichment analyses confirmed these negative changes in ECM-related transcripts in response to acute exercise, contrasting the effects of chronic resistance training [13,18].
The physical properties of ECM are distorted by disuse and aging, resulting in increased stiffness and potentially decreased force transmission and muscle efficiency [57,58]. Training-induced ECM remodeling may thus constitute an effective measure to reverse these adversities [50][51][52][53], and has been suggested to exert a protective role against injury [51]. However, available studies are ambiguous in their conclusions [52,53], and the link between observed changes in ECM-related gene expression and muscle biology and functionality remain uncertain. Adding to this, ECM remodeling-responses to training seems to be age-dependent [50,59], and also shows a clear dependency of time. Indeed, in the present study, the volume-associated differences in ECM-related transcriptome profiles disappeared entirely after twelve weeks of resistance training, whereby no genes were identified as showing volume-dependent responses. From a general perspective, this indicates that after prolonged training, the biological state of the muscle may have reached a new equilibrium, with lowand moderate-volume training having led to similar muscle phenotypic traits. This would imply that the benefits of higher training volume are restricted to augmented increases in muscle mass, perhaps facilitated by increases in ribosomal biogenesis [3]. Notably, this is likely to be an oversimplification. Indeed, at the protein level, in the current study, higher training volume led to a more robust phenotypic switch from type IIX → type IIA fiber also after twelve weeks of training [3]. Taken together, our data provide valuable directions for future research. It suggests that ECM remodeling is volume-dependent, at least during the initial part of a training program. This needs to be confirmed by studies in other populations, and its biological and functional significance needs clarification. Such studies should take advantage of the increased biological resolution of within-subject contralateral models.
Analyses of transcriptome responses to acute bouts of low and moderate resistance training volume revealed one single gene with volume-dependent changes in expression (RFT1). RFT1 is associated with the GO terms lipid transport, carbohydrate transport, and endoplasmic reticulum membrane. RFT1 expression has previously been shown to decrease in muscle immediately after training [22]. Although this warrants more research, one single gene arguably provides limited information. As such, rank-based enrichment tests identified three GO terms with volume-dependent changes in gene expression in the acute data set (RNA splicing, RNA localization, and covalent chromatin modification). Upon closer examination, these gene sets consisted of genes with both increased and decreased expression in the moderate-compared to the low-volume condition. Only a small fraction of these transcripts showed actual positive MSDs, indicating changes with unadjusted P-values < 0.05. While this enrichment analysis supports acute volume-dependent regulation of gene sets after resistance training at the selected time point of biopsy sampling (1 h after sessions), it remains plausible that such regulation would have been more pronounced at later time points.

Conclusions
Transcriptomic analyses of skeletal muscle subjected to altered growth conditions should account for global changes in mRNA to total RNA and cell density to accurately reflect biologically meaningful events. When accounting for such aspects, ECM remodeling showed volume-dependent responses to resistance training. These recommendations could be applicable to studies of other cell types and model systems undergoing increased or arrested growth. Also, the optimized selection of bioinformatic tools increases the biological relevance of transcriptome analyses from resistance-trained skeletal muscle.

Participants and study overview
The full study design has been previously described in detail [3]. Thirty-four participants completed a 12-week training-intervention with legs allocated to either low-(one set per exercise) or moderate-volume (three sets per exercise) training (Fig. 1a). Muscle biopsies were obtained from each leg prior (Week 0) to and after the intervention (Week 12), as well as prior to (Week 2 Pre-ex) and 60-min after (Week 2 Post-ex) the last training session of week 2, as previously described [3]. Participants with a complete set of highquality RNA samples (RQI ≥ 7, n = 25) were selected for RNA-seq analyses (Fig. 1b). Training-induced changes in muscle size and strength were estimated for each leg using several methods (for a complete overview, see [3]). Herein, we present dual-energy X-ray absorptiometry (DXA) measurement of lean mass for the 25 participants eligible for RNA-seq, as well as a weighted combined measure of strength (combining data from different strength tests) ( Table 1).

Training protocol
The training protocol consisted of unilateral lower body exercises (leg-press, leg-curl, and knee-extension). Each participant's leg was randomly assigned to perform either one or three sets per exercise (low-vs. moderate-volume), ensuring within-participant comparisons. Rest periods between sets were 90-180 s. The single-set leg was always trained in the rest period between the second and third sets of the multiple-set protocol. Training protocols were performed in a progressive manner, whereby resistance was continuously adjusted to ensure that the targeted number of repetitions was reached at volatile fatigue. This was equivalent to 10 repetitions maximum (RM) in weeks one and two, followed by 8RM in weeks three to five and 6RM in weeks six to twelve. Each week consisted of either 2 or 3 training sessions. From week four, weeks with three sessions contained one session at a sub-maximal load (90% of previous session load). All sessions commenced with a standardized warm-up. After each session, participants were given a standardized milk-based drink [3].

Muscle strength and lean mass assessments
Muscle strength was assessed twice before and once after the intervention. A detailed description of strength outcomes resulting from the study has been previously reported [3]. For the purpose of the present analyses, we present a weighted average of strength gains for the 25 participants eligible for RNA-seq, based on data from unilateral isometric and isokinetic (60°, 120° and 240° × s −1 ) knee extension, and one-repetition maximum (1RM) in unilateral knee extension and leg press, as previously reported [3]. Isometric and isokinetic strength was assessed using an individually adjusted dynamometer (Cybex 6000, Cybex International, Medway USA). 1RM was defined as the maximum load lifted through the full range of motion. From pre-intervention tests, the highest values were used for change score calculations. Limb lean-mass was assessed from full-body dualenergy X-ray absorptiometry (DXA; Lunar prodigy, GE Healthcare, Oslo, Norway) scans performed prior to and after the intervention. Limb lean-mass was derived from a segment covering the full leg from collum femoris to the distal end of the foot defined in the analysis software (enCore, GE Healthcare, Oslo, Norway).

Muscle tissue sampling, immunohistochemistry and RNA extraction
Muscle tissue was obtained bilaterally from m. vastus lateralis using a 12-gauge needle (Universal-plus, Medax, San Possidonio, Italy) under local anesthesia (Xylocaine, 10 mg × ml −1 with adrenaline 5 μg × ml −1 , AstraZeneca AS, Oslo, Norge). Samples were obtained from the two legs within 10 min of each other at all time-points. All rested state samples were obtained in the morning after a standardized breakfast. Resting samples obtained in Week 2 were sampled approximately 48 h after the fourth session. After the training period (Week 12), samples were obtained six days after the last training session and three days after the last strength assessment. Participants were instructed to ingest standardized meals during the 24 h leading up to the sampling event, and to refrain from strenuous physical activity during the last 48 h. Samples for immunohistochemistry (~ 15 mg) were transferred to a 4% formalin solution for fixation 24 − 72 h, before further preparation. Samples for RNA analyses (~ 25 mg) were dissected in ice-cold sterile saline solution (0.9% NaCl), blotted dry and snap-frozen using −80 °C isopentane, before storage at −80 °C until further processing.
For RNA extraction, the frozen muscle was homogenized in 1 ml of TRIzol reagent (Invitrogen, Life technologies AS, Oslo, Norway) using a bead homogenizer (Bullet Blender, Next Advanced, Averill Park, NY, USA). After phase separation, 400 μl of the aqueous phase was used in isopropanol precipitation of RNA, and after three washing steps (70% ethanol) the pellet was eluted in TE buffer. All samples showed a 260/280 nm ratio > 1.95 assessed by a spectrophotometer (NanoDrop 2000, ThermoFisher Scientific, Oslo, Norway). RNA integrity scores (RQI) were determined by capillary electrophoresis (Experion Automated Electrophoresis Station using RNA StdSens Assay, Bio-Rad). For each participant, all samples were extracted in the same extraction session in a randomized order. Participants with complete sets of high-quality RNA samples had an average RQI score of 9.0 (0.4), [full data set, 8.1 (2.1), range 1-9.7] (Fig. 1c). Notably, to achieve accurate normalization of qPCR data (and potentially also RNA-seq data), a commercially available exogenous RNA control (λ polyA External Standard Kit, Takara Bio Inc., Shiga, Japan) was added at a fixed amount per extraction prior to homogenization (0.04 ng ml −1 of Trizol reagent), as previously described [3,60]. Unfortunately, at present, we do not have access to the sequence of this spike-in, prohibiting its identification in RNA-seq data and rendering its subsequent utilization for normalization purposes difficult.

Illumina library preparation and sequencing
For each muscle sample, mRNA sequencing libraries were prepared from the same amount of RNA (1000 ng, depending on the minimum amount available) using TruSeq Stranded Total RNA Library Prep (Illumina, San Diego, CA USA). Paired-end sequencing (150 bp) was performed using an Illumina HiSeq 3000 (Illumina) at the Norwegian Sequencing Centre.

Modeling of gene counts
Gene counts were modeled using negative binomial generalized linear mixed models (GLMM), as suggested in [36], with modifications. Three model formulations were used to represent three different normalization scenarios. First, to account for fluctuations in RNAto-tissue ratios, the amount of tissue used in cDNA synthesis was included as an offset term together with the effective library size and study conditions (time and volume condition), added as a fixed effects in the model (tissue offset model). A simplified model contained only the effective library size together with study conditions, included as fixed effects (Effective library-size model). And finally, a naïve model formulation, without any form of normalization term, was used for comparisons. For acute exercise effects (pre-to post-exercise in the last session of week 2), only the library size normalized model was used as we expected that fluid shifts [38] could influence the muscle weight measurement, and changes in total RNA were unlikely to occur in this short time span [37]. The effective library size was calculated by multiplying the total library size with the RNA composition normalization factor, calculated using the trimmed mean method [25], followed by dividing the value by the median effective library size, as suggested by Cui et al. [36]. The effect of resistance training on gene counts was assessed as (1) the effect of exercise volume and (2) the effect of time. For analyses of the effect of exercise volume, differential expression was evaluated using GLMMs containing the interaction between time and exercise volume. For analyses of the effect of time, differential expression was evaluated using GLMMs containing only the time factor, combining all data irrespective of volume condition. In all models, a single random effect was used, giving each participant an individual intercept. Models were iteratively fitted using glmmTMB [64]. Utilization of the negative binomial distribution was supported by comparing the full model with a Poisson model (not containing the dispersion term), providing likelihood-ratio tests P-values that were distributed primarily below 0.05 (0.37% of models showed P > 0.05). Heteroscedasticity was assessed using the uniformity test in the DHARMa package [65], which generally showed good agreement with model assumptions, providing P-values concentrated near 1 (98.51% of models showed P > 0.05). Genes were identified as differentially expressed when the absolute Log2 fold-change was greater than 0.5, and the adjusted P-value was below 5%. P-values were adjusted per-model coefficient to control the false discovery rate [66].

Functional annotation
Enrichment analyses of gene ontology (GO) gene sets were performed using three approaches. First, a non-parametric rank test (described in [67] and implemented in the tmod package [68], version 0.40) was performed based on gene-specific minimum significant differences (MSD). MSD was defined as the lower limit of the 95% confidence interval (CI, based on estimated standard errors) around the Log fold-change (FC) when Log(FC) > 0 and the negative inverse of the upper 95% CI when Log(FC) < 0. This metric has been shown to have lower false-positive rates compared to other metrics applied during enrichment analyses [69]. As the MSD metric is positive when the CI does not overlap 0 and negative when overlap occurs, the rank test does not discern between up and down-regulated gene sets. A second approach, gene set enrichment analysis (GSEA) [70], was used to quantify the directional regulation of the gene set. GSEA was performed using the fgsea package [71] with Log(FC) as the gene level metric. Thirdly, over-representation analysis (ORA) was performed to assess if genes identified as differentially expressed (|Log2 fold-change|> 0.5 and adjusted P-values < 0.05; DE-genes) belonged to specific gene sets. ORA was performed using the enrichGO function in the clusterProfiler package [72], (version 3.16.0). GO gene sets (biological process, cellular component, and molecular function) were retrieved from the molecular signature database (version 7.1) [73].

Statistical analysis
Descriptive data are presented as mean and standard deviation (SD). Changes in muscle strength and mass were estimated using linear mixed models on change scores with baseline values as covariates. The performance of alignment tools was assessed by comparing log-differences between biological replicates, as suggested by Teng et al. [74], with modifications. Briefly, a subset of genes previously shown to be stably expressed between tissues was selected [32], whereupon log-differences between paired biopsy samples were calculated (i.e., using biopsies collected from each of the two legs prior to the training intervention). In addition, alignment tools were assessed by comparing relationships (Pearson's correlation coefficient) between gene family profiling of myosin heavy chains (MYH1, MYH2, and MYH7; muscle-specific) and the corresponding myosin heavy chain protein expression (measured using immunohistochemistry as fiber types IIX, IIA, and I). These mRNA and protein profiles were expressed as a fraction of total counts, thus removing the need for normalization of the RNA-seq data, as previously described for qPCR data [34]. Notably, these data also provided insight into the overall biological validity of the RNA-seq data, linking gene counts to protein phenotypes.
Additional file 1: Gene count estimates. Model (negative binomial generalized linear mixed models) based estimates of gene counts from models using different normalization strategies. ensemblid: Ensembl gene identifiers; normalization_model: Indicates normalization model: lib_size_normalized, gene counts are modeled with study conditions as fixed effects and participant id as random effects and with the addition of library size as a fixed effect. tissue_offset_lib_size_normalized, same as lib_size_normalized but with the addition of tissue weight used in library prep as an offset. The non-normalized model is not included in the data set; interaction_model: Indicates if the coefficient is estimated in a model containing the interaction between volume condition and time. If FALSE, the model only contains time as a fixed effect (in addition to any normalization, see above), coefficients should in this case be interpreted as averages over volume conditions; coefficient: Names of coefficients (fixed effects, time and volume conditions). timew2pre and timew12 indicate differences from timew0 (intercept) in rested state models. timew-2post indicates differences from timew2pre (intercept) in acute exercise models. setsmultiple indicates interaction effects, the difference between setssingle (reference level) and setsmultiple; estimate: Estimates on the natural log scale; se: Standard errors (SE) on the natural log scale; zvalue: Z-values; pvalue: Un-adjusted P-values; pvalue_adjust: Adjusted P-values (FDR) per model/normalization method and coefficient.
Additional file 2: Functional annotation using gene ontology terms. Significance tests for functional annotation using rank-based and over-representation based analysis. ID: Gene ontology (GO) id; go_category: GO category, bp, Biological process; cc, cellular component; mf, molecular function; name: Descriptive GO name; normalization_ model: Indicate normalization model used to test enrichment: lib_size_normalized, gene counts are modeled with study conditions as fixed effects and participant id as random effects and with the addition of library size as a fixed effect, tissue_offset_lib_size_normalized, same as lib_size_normalized but with the addition of tissue weight used in library prep as an offset. The non-normalized model is not included; coefficient: Names of coefficients used to test enrichment (fixed effects, time and volume conditions). timew2pre and timew12 indicates differences from timew0 (intercept) in rested state models. timew2post indicates differences from timew2pre (intercept) in acute exercise models. setsmultiple indicates interaction effects, the difference between setssingle (reference level) and setsmultiple; cerno_statistic: Test statistic from the rank-based cerno-test; cerno_auc: Area under the curve from the Cerno test (see tmod documentation for details, https ://CRAN.R-proje ct.org/packa ge=tmod); cerno_pval: Un-adjusted P-values from cerno test, cerno_padj: Adjusted P-values (default settings in tmodCERNOtest, see tmod documentation); fgsea_pval: Un-adjusted P-values from gene set enrichment tests performed with the fgsea package; fgsea_ padj: Adjusted fgsea P-values; NES: Normalized enrichment scores from gene set enrichment tests; set_size: Size of gene sets using genes expressed in the present data set; ora_geneRatio: Gene ratio in over-representation analysis of genes identified as differentially expressed; ora_padj: Adjusted p-values from over-representation analysis.