To generate the LLC model of cancer cachexia, we used 8-week old, healthy C57BL/10 male mice obtained from a breeding colony, maintained by our institutional animal care facility at the São Paulo State University (UNESP, Botucatu, São Paulo, Brazil). The animals were housed in cages under a 12 h light/dark cycle with food and water ad libitum. Before inoculation into mice, LLC cells (ATCC^® CRL-1642^TM) were cultured in Dulbecco’s modified Eagle’s medium (DMEM, Thermo Fisher Scientific, United States) supplied with 10% of fetal bovine serum (FBS, Thermo Fisher Scientific, United States) and 1% penicillin/streptomycin (Thermo Fisher Scientific, United States) and maintained in a 5% CO2, 37°C humidified incubator. After 3 days of acclimation to their environment, mice were randomly assigned into two groups. Twenty mice (LLC group) were inoculated subcutaneously with a total of 1.5 × 10⁶ LLC cells (7.5 × 10⁵ cells in 0.1 mL PBS in each flank). Ten mice (control group) were injected with equal volumes of 1X PBS.

Lewis Lung Carcinoma and control mice were weighed daily and studied 22 days after LLC cells inoculation when the LLC group had developed overt cachexia. The animals were euthanized upon anesthesia with intraperitoneal ketamine and xylazine (100 and 10 mg/kg, respectively), and tumor diameter, tumor weight, and body weight were measured. The tibialis anterior (TA), soleus (SOL), and gastrocnemius (GAS) muscles were collected, weighed, and then snapped frozen in liquid nitrogen for further analyses. These muscles were chosen based on their fiber composition: TA muscle has a higher proportion of glycolytic fast-twitch fibers, SOL muscle has a higher proportion of oxidative slow-twitch fibers, and the GAS muscle, which comprises both fiber types.

The study was performed following the guidelines of the Control of Animal Experimentation and Ethical Principles in Animal Research (CONCEA - National Council for Control of Experimental Animals), under the approved protocol n° 702, emitted by the Institute of Bioscience of Botucatu Ethics Committee on Animal Use, from the São Paulo State University (UNESP, Brazil).

RNA extraction was performed using the TRIZOL reagent (Thermo Fisher Scientific, United States), according to the manufacturer’s instructions. RNA was quantitated by spectrophotometry (NanoVue; GE Healthcare Life Sciences, United States), and its integrity was ensured by obtaining RNA Integrity Number (RIN) > 8 (Agilent 2100 Bioanalyzer; Agilent Technologies, Germany). RNA samples were treated with DNA Free Kit (Thermo Fisher Scientific, United States) to remove genomic DNA contamination.

We randomly selected samples (four control and six LLC mice) to construct the RNA sequencing libraries with the TruSeq Stranded Total RNA Sample Prep Kits (Illumina, United States), using 1000 ng of total RNA. These six LLC mice are those that survived the entire period of the experiment (22 days). Samples were indexed with adaptors and submitted for paired-end 2 × 100-bp sequencing using a HiSeq 2000 instrument (Illumina, United States). Sequencing was performed on ten RNA samples in one lane of the flow cell, following the manufacturer’s instructions. The lane produced ∼ 600 million raw paired reads. The data output in the fastq file format contained sequence information, including the sequencing quality (Phred quality score). Average Phred scores ≥20 per position were used for alignment.

The TruSeq Small RNA Sample Preparation kit (Illumina, United States) was used to prepare the miRNA-Seq libraries for the same set of samples used for mRNA-Seq, following the manufacturer’s instructions. miRNA libraries were 50 bp single-end sequenced using a HiSeq2000 instrument (Illumina, United States). Sequencing was performed on ten RNA samples in one lane of the flow cell.

Paired-end reads for mRNA were mapped to the mm10 genome using TopHat2 (Kim et al., 2013) with the following options: –mate-inner-dist 200 –mate-std-dev 100 –no-novel-juncs –min-intron-length 40. Single-end reads for miRNA were mapped to the miRBase version 21 using Bowtie (Langmead et al., 2009) with the following options: -n 0 -l 8 -a –best –strata –phred33-quals. Counts for RefSeq genes were obtained using HTSeq (Anders et al., 2015) with the default settings, and DESeq2 (version 1.4) (Love et al., 2014) was used to normalize expression counts. Fold change was calculated as the ratio of normalized counts for each sample against the average of the reference group. Genes were considered differentially expressed if | fold change| (FC) ≥ 1.5, and p-values ≤ 0.05. For assess mRNA transcript abundance, reads were converted to reads per thousand base pairs peak per million mapped reads (RPKM). For assessing miRNA abundance, the reads were converted to counts per million mapped reads (CPM). Finally, we used dot plots to demonstrate the expression of selected genes either previously associated with muscle atrophy in cachexia or that we hypothesized may be associated with muscle atrophy, based on the literature.

We standardized the normalized counts of each gene and applied k-means clustering using Euclidian distance and with random initialization and 10000 executions, and finally selecting four clusters.

Pscan web interface (Zambelli et al., 2009)¹ was used to detect DNA motifs overrepresented in the promoter of the differentially expressed genes. Gene promoters were considering between nucleotides −300 and +50 relative to the Transcription Start Site (TSS). Significance was tested against CpG-content-matched promoters as background. Binding sites were considered significantly overrepresented when the p-value <0.01.

Gene Ontology enrichment was performed using the ClueGO Cytoscape plugin (Bindea et al., 2009), using a hypergeometric test with a Benjamini-Hochberg False Discovery Rate correction (Benjamini and Hochberg, 1995). A p-value cut-off of 0.05 was used to identify enriched terms.

Candidate miRNA-mRNA targets relationships were predicted by at least one or more of the following target prediction algorithms (union set) extracted from mirDB (Wang, 2008), TargetScan 5.1 (conservation and non-conservation sites) (Garcia et al., 2011), DIANA-microT (Maragkakis et al., 2009), and PicTar (4-way, and 5-way) (Krek et al., 2005). Additionally, we used validated targets deposited in miRTarBase (Hsu et al., 2011). We also filtered our data using differentially expressed genes (mRNA and miRNA) identified by RNA-Seq, considering that mRNA and miRNA expression levels should be inversely correlated.

Based on the differentially expressed genes, protein-protein interaction networks were generated using the STRING database (Snel et al., 2000; Szklarczyk et al., 2017)², which also detects functional interactions among the corresponding genes. Network visualization was performed using the open-source software platform Cytoscape (Shannon et al., 2003) v. 3.6.1³.

Cryostat transverse sections of the TA muscle (10 μm thick) were collected from control and LLC tumor-bearing mice. Collected samples were placed on the same slide to minimize staining differences; sections were incubated with a saturated picric acid solution followed by Picrosirius red (0.1% Sirius red in saturated picric acid) for 3 min, dehydrated, and mounted in Permount. Eight color pictures per sample were captured using a light microscope (Olympus, Japan). The light intensity parameters used were the same for all samples. Picrosirius staining areas were assessed using the Image J software. As previously described, picrosirius staining areas were normalized for cell density to account for individual fiber atrophy (Kirby et al., 2005). To quantify muscle tissue disorganization, we employed fractal dimension analysis by binarizing photographs using ImageJ software, as previously described (Pacagnelli et al., 2016). Briefly, the fractal dimension was estimated using the box-counting tool (ImageJ software), which quantifies pixel distribution in the space without considering image texture. The fractal dimension value is expressed from 0 to 2, where values close to 2 represent higher tissue disorganization.

Muscle proteins were extracted using Tris-Triton buffer (10 mM Tris pH 7.4, 100 mM NaCl, 1 mM EDTA, 1 mM EGTA, 1% Triton X-100, 10% glycerol, 0.1% SDS, 0.5% deoxycholate) containing Protease Inhibitor Cocktail (Sigma-Aldrich, United States) and quantified by the Bradford method (Bradford, 1976). Subsequently, the Lammeli buffer (Sigma-Aldrich, United States) was added to each sample and boiled at 100°C for 10 min. Proteins were subjected to SDS-PAGE on 10% polyacrylamide gels. After electrophoresis, proteins were electrotransferred to nitrocellulose membranes (Bio-Rad, United States) for 2 h at 120V. The membranes were blocked with 5% non-fat dry milk diluted in TBS-Tween for 2 h, and then incubated overnight at 4°C with collagen I (1:100 dilution, sc-25974, Santa Cruz, CA, United States) or β-actin (1:1000 dilution, sc-81178, Santa Cruz, CA, United States) antibodies. Secondary antibodies conjugated with horseradish peroxidase (HRP) and ECL chemiluminescent detection (GE Healthcare, United States) system were used for visualization of the blots in ImageQuantTM LAS 400 (GE Healthcare, United States). We quantified the blots by densitometry using ImageJ software, and collagen I values were normalized to β-actin.

Data are expressed as mean ± standard deviation (SD). Statistical analysis was performed using the GraphPad Prisma software v 6.07 (GraphPad Software, Inc., United States). For all statistical analyses not described elsewhere, Student’s t-test was applied to compare the groups. Statistical significance was considered achieved when the p-value was <0.05.

As expected, all mice subcutaneously inoculated with Lewis Lung Carcinoma (LLC) cells developed cancer cachexia (LLC group) compared to mice injected with PBS (control group). The tumor was detected by palpation after 7 days of cell injection in LLC, after 15 days, the tumor site was visually identified as a skin projection; and after 22 days, when the animals were euthanized, tumor mass was observed under the skin. After euthanasia, the surgically exposed tumor was solid, vascularized, roughly spherical, measuring ∼2 cm in diameter, and weighing ∼4 g (Supplementary Figures S1A,B).

Although the LLC cell line is highly tumorigenic, metastases were not visually identified. Five out of twenty mice died during the experiment (25%) (Figure 1A). Consistent with cachexia syndrome, the LLC group exhibited 14.4% body weight (BW) loss after 22 days of LLC cell injection compared to the control group (Figure 1B and Supplementary Table S1). This BW reduction was associated with a loss of adipose tissue (Figure 1C and Supplementary Table S1) and skeletal muscle mass (Figure 1D and Supplementary Table S1). Cachexia was further confirmed by splenomegaly (Figure 1E, Supplementary Table S1, and Supplementary Figure S1C). Finally, among the studied muscles – tibialis anterior (TA), soleus (SOL), and gastrocnemius (GAS) – we selected TA for further analysis because its weight presented strong correlations between tumor weight (inverse correlation) and BW (positive correlation) (Figure 1F).

Transcriptome analysis of TA muscle revealed 11,436 genes out of the nearly 45,000 mouse RefSeq genes. Principal Component Analysis (PCA) was able to discriminate LLC and control samples (Figure 2A). High heterogeneity in mRNA profiles of the LLC group was evidenced by spatial dispersion. We found 1008 differentially expressed genes (DEGs), of which 487 and 521 were up- and down-regulated, respectively (Supplementary Table S2). We validated by RT-qPCR the expression of some genes involved in myogenesis, sarcomere, proteasome, and ECM (Supplementary Figure S2). Unsupervised hierarchical clustering analysis showed that DEGs were organized in four clusters (I–IV), according to their direction and variability of expression within LLC samples (Figure 2B): clusters I (n = 386) and II (n = 101) contain up-regulated genes, while clusters III (n = 157) and IV (n = 364) include down-regulated genes. Considering the gene expression variability within LLC samples, clusters II and IV contain genes that are uniformly regulated (fold change% CV: 32.62 ± 15.50 and 32.49 ± 15.90, respectively Figure 2C); whereas clusters I and III contain genes with high variability in expression levels among LLC samples (fold change% CV: 65.60 ± 38.82 and 54.97 ± 15.42, respectively, Figure 2C).

Next, we explored the identity of the genes found within the clusters. Cluster I include up-regulated genes that have variable expression. The genes in cluster I are associated with the proteasome complex (e.g., Trim63, Fbxo32, Ubc, Ubb, Psmd4, Psma7, Psmc2, Fbxo31, and Ube4a), autophagosome (e.g., Ctsl and Retreg1), and the translation inhibitor Eif4ebp1 (Figure 2D and Supplementary Figure S3A). Remarkably, cluster I includes the interleukin-6 receptor (Il6ra), which presented a 5-log range in expression variability within LLC samples (Figure 2D). Cluster III contains down-regulated genes that have variable expression; this cluster includes genes associated with ECM (e.g., Has3, Col15a1, Col22a1, Col9a1, Cpq, and Mmp15) and muscle metabolism and contraction (e.g., Myom3, Myh2, Myl3, Myoz2, Myh7, Lmod3, Synpo2, and Myl1) (Figure 2D and Supplementary Figure S3A). Cluster II comprises up-regulated genes that are uniformly regulated within LLC samples, which are associated with the immune system (e.g., Selp, S100a9, Il1b, Cxcr2, and Csf3r), ECM organization (e.g., Mmp9, Mmp8, and TLL1), and apoptotic process (e.g., Cd300lf, Scarn1, Lcn2, Gadd45g, and Chac1) (Figure 2D and Supplementary Figure S3A). Cluster IV contains down-regulated genes that are uniformly regulated within LLC samples, which are associated with ECM and sarcomere (Figure 2D and Supplementary Figure S3A). Remarkably, both cluster III and IV contain genes related to the ECM (collagens) and sarcomere (myosins), but with differences in gene expression stability (low and high, respectively) within LLC samples. Considering the variability in gene expression profiles within the LLC samples, we determined a reduced set of DEGs able to differentiate cachectic and control groups. For this, we use Euclidian distances, and we found that samples L2, L3, and L4 are part of a cluster that we call cluster A which is the closest to the control group (Supplementary Figure S3B). Next, we determined the number of DEGs between the subgroup of samples. We found 443 DEGs that were sufficient to effectively segregate both LLC and control groups (Supplementary Figure S3D), presenting high intragroup stability (fold change% CV: 28.95 ± 10.02 and 27.30 ± 11.06, for up- and down-regulated genes, respectively) (Supplementary Figures S3E,F). This set of DEGs was determined using samples from cluster A and comprises genes associated with the ECM, proteolysis, and inflammatory response (Supplementary Figure S3G). This set of 443 genes represents a reduced number of deregulated genes in all cachectic muscle samples, regardless of the variability in the gene expression within LLC samples.

We considered as biologically relevant for the muscle wasting those most abundant transcripts with the highest degree of regulation. Initially, we used a scatter plot that integrated the degree of regulation (fold change; FC) and abundance (Reads Per Kilobase Million; RPKM) of all transcripts (Figure 3A). Abundant transcripts presented subtle changes in gene expression when compared to rare transcripts. Notably, the up-regulated genes Trim63, Fbxo32, and Ubc - associated with the proteasomal degradation pathway - presented the highest abundance and degree of change in expression (Figure 3A). Additionally, we identified the up-regulated genes Ddit4 and Eif4ebp1 (Figure 3A), which have been previously implicated in protein synthesis during skeletal muscle atrophy (Bodine et al., 2001). We also found up-regulation in the expression of the antioxidant genes Gpx3, Mt1, and Mt2 (Figure 3A), which have been described with a role in muscle repair in atrophic conditions (Nuoc et al., 2017; Summermatter et al., 2017). We identified the down-regulation of the sarcomere genes Myl1, Myh1, Fhl1, Gsn, Myl3, and Actc1 in cachectic muscle samples (Figure 3A). Interestingly, we found that the muscle-specific myoglobin (Mb) transcript was highly abundant and is down-regulated cachexia (Figure 3A). Also, a high degree of down-regulation of the ECM genes Col3a1, Thbs4, Col6a1, Col1a1, Col1a2, and Col6a2 was identified. Finally, we detected low abundance transcripts with a high degree of regulation, including Il6ra and the ECM remodeling genes Mmp9 and Mmp8 (Figure 3A).

Gene ontology analysis revealed 33 terms that were clustered in eight main modules that are relevant to muscle wasting, which include genes associated with the scaffold, cellular metabolism, cellular signaling, cellular differentiation, cellular immune system process, cellular respiration, proteolysis, and ion regulation (Figure 3B and Supplementary Table S3). Notably, some novelties were found, such as the negative regulation of cell junctions (e.g., gap and tight junctions), carbohydrate metabolism (e.g., glycolytic process), cell differentiation (e.g., axonogenesis, angiogenesis, and PDGF signaling), and positive regulation of the immune cell system (e.g., neutrophil and leukocyte chemotaxis) (Figure 3B). Moreover, previously cachexia-associated terms such as negative regulation in the sarcomere, cell migration, and ECM genes, as well as positive regulation of genes involved in the proteasome complex, autophagy, IL-6 signaling, and cell differentiation were detected (Figure 3B). We also identified the percentage of up- and down-regulated genes in each ontology term (Figure 3B). Consistent with the atrophic phenotype, this analysis demonstrated that all deregulated genes related to proteasome complex were up-regulated, while most of the deregulated genes related to myofibril and ECM were down-regulated (Figure 3B). Considering the transcriptome variability across the LLC samples, we also asked which pathways were enriched in LLC transcriptome when compared with the reduced set of genes enriched explicitly in the LLC subgroup A. This analysis confirmed changes in the expression of genes associated with protein degradation such as proteasome complex, autophagosomes, and macroautophagy (Figure 3C).

The transcriptional profile can provide a step toward the identification of key transcription factors that regulate gene expression. We performed an enrichment analysis of transcriptional motifs in the promoter sequences of 1008 differentially expressed genes (Figures 3D,E). The promoters of the up-regulated genes revealed motif enrichment for the Forkhead transcription factor (FoxO) (Figure 3D). However, when we analyzed the changes in expression and abundance of FoxO family members genes, only FoxO1 and FoxO6 were up- and down-regulated, respectively (Figure 3F). The promoters of the up-regulated genes also revealed binding sites for transcriptional factors within the NF-κB and STAT families (Figure 3D). Additionally, components of the NF-κB and STAT families were up-regulated: Rela, RelB, Nfkb2, and Stat3 (Figure 3F). We found enrichment of the AP-1 transcription factor, as well as the up-regulation of Junb and Fosl2, which are translated into proteins that constitute the AP-1 heterodimer (Figure 3F). We also found enrichment of other transcriptional factors without a change in their expression (Figures 3D,F); among these are transcription factors related to cell cycle regulation (E2f3, Yy1, and Creb3), unfolding protein response (Xbp1), and the SMAD family (Smad2, Smad3, Smad4, and Smad5).

Interestingly, promoters of down-regulated genes also revealed motif enrichment of transcriptional factors related to myogenesis (Myf6, Myod1, Myog, Tcf12, Pbx1, lbx1, Nfix, and Nfic), lipid homeostasis (Rora and Rorc), energy metabolism (Med1), and muscle fiber-type specification (Six1, Six2, Tead1, Tead3, Tead4, Egr1, Klf3, Hsf1, Hsf2, and Hsf4) (Figure 3E). However, only genes coding for the transcription factors Myf6, Myog, Rorc, Rora, and Egr1 changed their expression (Figure 3F).

Out of 1915 mature miRNAs, 302 were expressed in skeletal muscle (mapped reads >32 in at least one of the sequenced samples). Eighteen miRNAs were differentially expressed (FDR ≤ 0.05 and |fold change| ≥ 1.5) in muscle wasting during cancer cachexia in comparison to controls (13 up and five down-regulated, Figure 4A). PCA and clustering analysis showed that these 18 miRNAs were poorly clustered samples according to their experimental groups (Figure 4B) when compared to the clear clustering found in PCA of the mRNA expression data (Figures 2A,B). Furthermore, 44% of differentially expressed miRNAs in LLC muscle samples were expressed at low levels and the degree of regulation (Figure 4C). Notably, miR-10b-5p was regulated at high levels in the atrophying muscles (Figure 4C). The differentially expressed miRNAs also included miR-29b-3p, miR-146a-5p, miR-146b-5p, and miR-181c-3p, which have been previously studied in the skeletal muscle context (Naguibneva et al., 2006; Li et al., 2017; Sun et al., 2017). Interestingly, the MyomiRs mir-208a, mir-208b, mir-499, miR-133a, miR-133b, and miR-1 were not differentially expressed in the LLC group.

To improve the accuracy of our in silico mRNA target prediction used to identify potential mRNA targets of the differentially expressed miRNAs, we considered opposite directions of deregulated expression between miRNA and target mRNAs in the same set of samples. We found a network with 171 interactions between 18 miRNAs and 131 target genes (Figure 5A). This analysis revealed that the upregulated miRNA miR-350-3p has a higher number of target genes (n = 47). Interestingly, miR-29b-3p presented 22 potential targets, including many genes that encode proteins related to the ECM. Additionally, we found that repressed miRNAs do not share target genes. Furthermore, some genes such as Map2k6, Ptpn3, Mettl21c, Plxdc2, Ppargc1b, Rgs5, and Vegfa were found to be co-regulated by up to three upregulated miRNAs (Supplementary Table S4).

Based on the integrative miRNA-mRNA analysis, we identified enriched pathways for deregulated genes targeted by differentially expressed miRNAs (Figure 5B). Gene ontology analysis revealed miRNA interactions affecting genes regulating mainly the ECM, but also cell migration, transcription factor binding, ion transport, and FoxO signaling. To elucidate the functions of these complex interactions between mRNAs and miRNAs in cancer cachexia, we constructed a regulatory network displaying predicted and validated interactions between the miRNAs and target mRNAs, considering physical and pathway protein-protein interactions. We found sub-networks such as those related to ECM organization (Figure 5C), cell migration (Figure 5D), and transcription factors (Figure 5E). The ECM organization network (Figure 5C) contains a set of nine collagen genes, including validated targets of the miR-29b-3p. Furthermore, we found predicted interactions for miR-1843a-3p, miR-350-3p, miR-223-p, and miR-3535 with ECM components such as Col6a1, Timp2, Mmp15, Dcn, and Actn2. These identified networks share the miRNAs mir-29b-3p, mir350-3p, and miR-3535, suggesting a pleiotropic effect of these miRNAs on the ECM of atrophying muscles in cancer cachexia.

Together, our results point out to a crucial role of ECM remodeling in skeletal muscle atrophy in cancer cachexia. TA muscle cross-sectional area stained with Picrosirius red in LLC tumor-bearing mice presented reduced total collagen deposition (Figures 6A,B) and ECM disarrangement (Figure 6C), as revealed by picrosirius red staining area and fractal dimension analysis, respectively. The reduction of ECM was further confirmed by the reduced protein levels of collagen alpha-1 type I collagen (COL1A1), one of the main structural components of the ECM in skeletal muscle tissues (Figure 6D). We also found ECM genes deregulated in the gastrocnemius and soleus muscles (Supplementary Figure S4).