MicroRNA-mRNA Co-sequencing Identifies Transcriptional and Post-transcriptional Regulatory Networks Underlying Muscle Wasting in Cancer Cachexia

Cancer cachexia is a metabolic syndrome with alterations in gene regulatory networks that consequently lead to skeletal muscle wasting. Integrating microRNAs-mRNAs omics profiles offers an opportunity to understand transcriptional and post-transcriptional regulatory networks underlying muscle wasting. Here, we used RNA sequencing to simultaneously integrate and explore microRNAs and mRNAs expression profiles in the tibialis anterior (TA) muscles of the Lewis Lung Carcinoma (LLC) model of cancer cachexia. We found 1,008 mRNAs and 18 microRNAs differentially expressed in cachectic mice compared with controls. Although our transcriptomic analysis demonstrated a high heterogeneity in mRNA profiles of cachectic mice, we identified a reduced number of differentially expressed genes that were uniformly regulated within cachectic muscles. This set of uniformly regulated genes is associated with the extracellular matrix (ECM), proteolysis, and inflammatory response. We also used transcriptomic data to perform enrichment analysis of transcriptional factor binding sites in promoter sequences, which revealed activation of the atrophy-related transcription factors NF-κB, Stat3, AP-1, and FoxO. Furthermore, the integration of mRNA and microRNA expression profiles identified post-transcriptional regulation by microRNAs of genes involved in ECM organization, cell migration, transcription factors binding, ion transport, and the FoxO signaling pathway. Our integrative analysis of microRNA-mRNA co-profiles comprehensively characterized regulatory relationships of molecular pathways and revealed microRNAs targeting ECM-associated genes in cancer cachexia.


INTRODUCTION
Cancer cachexia is a multifactorial syndrome characterized by an ongoing loss of skeletal muscle mass that affects up to 80% of cancer patients, depending on the cancer type (Dhanapal et al., 2011), and it represents the direct cause of at least 20% of cancer deaths (Loberg et al., 2007). Muscle wasting in cancer patients leads to a significant weight loss that affects the quality of life, response to radio-and chemotherapy, tolerance to treatment, and survival (Dewys et al., 1980;Fearon et al., 2013;Martin et al., 2013;Vaughan et al., 2013;von Haehling et al., 2016). Loss of skeletal muscle mass, the most apparent symptom of cachexia, is mainly related to metabolic dysregulation characterized by resistance to anabolic signals, increased energy expenditure, and catabolism (Porporato, 2016).
Several studies have linked cachexia to increased plasma levels of proinflammatory cytokines such as interleukin (IL) 1β, IL6, tumor necrosis factor-alpha (TNF-α), and interferongamma (IFN-γ) (Fearon et al., 2012). These molecules alter the ubiquitin-proteasome, insulin-like growth factor 1 (IGF1)-Akt, and autophagy-lysosome pathways, which ultimately lead to an imbalance between muscle protein synthesis and degradation (Fearon et al., 2012). However, the molecular mechanisms by which these cytokines mediate muscle wasting are not entirely understood.
Recent advances in mRNA transcriptome analysis have helped to unveil new molecular pathways, potential therapeutic targets, and biomarkers in cancer cachexia (Monitto et al., 2001;Stephens et al., 2010;Bonetto et al., 2011;Gallagher et al., 2012;Fontes-Oliveira et al., 2014;Judge et al., 2014;Shum et al., 2015;Fukawa et al., 2016). Genome-wide expression microarrays studies have also highlighted the importance of microRNAs (miRNAs) as an additional level of complex regulatory networks that posttranscriptionally control gene expression in muscle-wasting conditions (Eisenberg et al., 2007;Agarwal et al., 2013;Shen et al., 2013;Soares et al., 2014). While these studies have initially shown that global miRNA transcriptional profiles in wasting muscles seem to be peculiar to different catabolic conditions (Soares et al., 2014), it was subsequently demonstrated that the miRNA miR-29b contributes to multiple types of muscle atrophy, including those found in cardiac and cancer cachexia Moraes et al., 2017). Skeletal muscle miRNA transcriptome profiling using RNA sequencing (RNA-Seq) is also highly informative in cancer cachexia and has identified novel miRNAs associated with the syndrome in patients with pancreatic and colorectal cancer (Narasimhan et al., 2017). This same study also reanalyzed two independent mRNA microarrays datasets and identified that the miRNA-targeted transcripts, in a tissue-specific context, are involved in myogenesis and inflammation pathways. In a mouse model of cancer cachexia, miRNA sequencing combined with bioinformatics prediction analyses revealed that wasting muscles change the expression of miRNAs targeting transcripts that are essential for determining muscle size .
Although informative, these previous transcriptomic studies addressing muscle wasting in cancer cachexia did not explore different levels of gene expression regulation by integrating mRNAs and miRNAs RNA-Seq data from the same set of muscle samples. Nevertheless, miRNA-mediated regulation of gene expression in cellular networks involves complex interactions among various miRNA targets through several mechanisms (Flynt and Lai, 2008;Dragomir et al., 2018), which can be better addressed by simultaneous analysis of intrinsic interactions of biological entities across multiple omics layers. Thus, the identification of multi-omics features observed on the same set of samples provides a unique possibility for elucidating transcriptional and post-transcriptional regulatory networks during muscle wasting in cancer cachexia.
In the present study, we aimed to explore and integrate paired miRNA and mRNA co-profiles during skeletal muscle wasting in a mouse model of cancer-induced cachexia. This comprehensive analysis allowed the identification of new miRNA targets and regulatory strategies controlling gene expression in skeletal muscle atrophy. Although the atrophic muscles of cachectic mice showed high heterogeneity in the transcriptional profile, we successfully identified a set of differentially expressed genes uniformly regulated within these cachectic samples. We also characterized the regulatory relationships of molecular pathways, including miRNAs targeting extracellular matrix (ECM)-associated genes. Taken together, our results show that ECM-associated genes are dependent on inflammatory signaling and reveal miRNA-mRNA molecular networks that may play a role in the development of muscle wasting in cancer cachexia.

Lewis Lung Carcinoma (LLC) Model of Cancer Cachexia
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 R 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 6 LLC cells (7.5 × 10 5 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).

Total RNA Isolation
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.

Preparation and Processing of mRNA-Seq Libraries
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 pairedend 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.

Preparation and Processing of miRNA-Seq Libraries
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.

Read Alignment and Differential Gene Expression Analysis
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 -minintron-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.

Clustering Analysis of the RNAseq Expression Data
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.

Motif Analysis
Pscan web interface (Zambelli et al., 2009) 1 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-contentmatched promoters as background. Binding sites were considered significantly overrepresented when the p-value <0.01.

Gene Ontology (GO) Enrichment Analysis
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.

miRNA Target Prediction
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.

Interaction Network
Based on the differentially expressed genes, protein-protein interaction networks were generated using the STRING database (Snel et al., 2000;Szklarczyk et al., 2017) 2 , 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 3 .

Picrosirius Red Staining
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.

Western Blotting Analysis
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.

Statistical Analysis
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 2 http://string-db.org/ 3 https://cytoscape.org/ 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.

Lewis Lung Carcinoma Cells Induced Cachexia in Mice
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).

Comprehensive Transcriptome Characterization of Muscle Wasting
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 downregulated 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.

Relevant Genes and Regulatory Pathways Associated With Muscle Wasting in Cancer Cachexia
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 downregulated ( 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).

Transcriptional Factors Motifs Enriched During Muscle Wasting in Cancer Cachexia
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).

miRNAs Associated With Muscle Wasting in Cancer Cachexia
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 downregulated, 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 FIGURE 3 | Relevant genes and regulatory pathways associated with muscle wasting in cancer cachexia. (A) Scatterplot comparing abundance (RPKMs, x-axis) and their degree of expression (log2-fold change, y-axis). Red and blue dots represent up-and down-regulated genes (fold-change > 1.5 and FDR < 0.05; Wald test), respectively. These genes were either previously associated with muscle atrophy in cachexia or that we hypothesized may be associated with muscle atrophy, based on the literature. Vertical lines represent thresholds for low, medium, and high abundance genes as defined by quartiles (Q1 and Q3) (B) Gene-term enrichment analysis of differentially expressed genes (DEGs) in the tibialis anterior of Lewis Lung Cancer (LLC) tumor-bearing mice showing the top canonical pathways. The colored horizontal bars represent the percentage of genes presented in the dataset compared to the total number of genes in each term (% Genes/Term). The fraction of up-and down-regulated genes (horizontal bars) in each term are shown in red and blue, respectively. The vertical colored bars (y-axis) represent major gene terms modules. (C) Gene ontology analysis of DEGs from LLC vs. LLC subgroup A samples. Each horizontal black bar represents the ontology term fold enrichment compared to the total number of genes in each term. De novo motif analysis was performed on promoters (-300 and +50 relative to Transcription Start Site, TSS) of up-(D) and down-regulated (E) genes. Motifs were compared using the transcription factor JASPAR database to determine the closest annotated matches. Percentage (%) represents a fraction of foreground (Fg) and background (Bg) sequences that contain at least the occurrence of one motif. (F) Scatterplot comparing abundance (Basal RPKMs, x-axis) and their degree of expression (log2-fold change, y-axis). Each dot represents a differentially expressed genes (DEGs; fold-change > 1.5 and FDR < 0.05; Wald test) encoding for transcription factors (red and blue dots, up-and down-regulated genes, respectively). Vertical lines represent thresholds for low, medium, and high abundance genes, as defined by quartiles (Q1 and Q3). 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.

Integrative Analysis Revealed a Set of ECM mRNAs Regulated by miRNAs in Cancer Cachexia
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.

ECM Remodeling 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).

DISCUSSION
Despite advances in the study of cancer cachexia, its pathogenesis is complex and remains incompletely understood (Baracos et al., 2018). Thus, it is necessary to identify signaling pathways as well as transcriptional and post-transcriptional events underlying skeletal muscle atrophy in this condition. We used paired microRNA-mRNA co-profiles in wasting muscles of cachectic mice, which unveiled ECM remodeling events potentially regulated by miRNAs. Although our transcriptomic analysis demonstrated a high heterogeneity in mRNA profiles of cachectic mice, we successfully identified a reduced number of differentially expressed genes that were uniformly regulated with low variability within cachectic samples. Thus, in addition to the well-known down-regulation of sarcomere proteins genes, other genes encoding ECM structural proteins that are potentially regulated by miRNAs may contribute to the development of skeletal muscle wasting in cancer cachexia.
Our RNAseq data significantly expands previous genomewide studies in cancer cachexia by integrating mRNA and microRNA transcriptome profiling from the same set of muscle samples. This strategy allowed us to identify miRNA targets with higher accuracy. Moreover, we used a high number of biological replicates (six cachectic and four controls), which added higher precision and sensitivity for the identification of transcriptional and post-transcriptional events. Our transcriptomic data reliably differentiated muscle samples from cachectic and control mice. Notably, inter-individual variations in mRNA expression were evidenced in cachectic mice, which may be linked to individual genetic factors and stochastic tumor growth events that may determine the development and progression of muscle wasting. Studies in rodent models show heterogeneity in the occurrence of cachexia (Matsuyama et al., 2015;Norden et al., 2015). This characteristic is also found in human neoplasms, which shows variability in the prevalence and severity of cachexia among patients with the same diagnosis and cancer stage (Prado et al., 2013;Baracos et al., 2018). Importantly, we found some genes commonly associated with cachexia, which presented high variability in expression levels among cachectic mice. These data may help to explain the variability in the prevalence and severity of the syndrome. For example, the highly variable upregulated genes are related to protein catabolism (e.g., Trim63 and Fbxo32, Figure 2), while the highly variable down-regulated genes are mainly associated with sarcomere and ECM. This variability in the expression profile of these sets of genes within LLC samples may help to explain why protein catabolism genes in human studies have failed to recapitulate the findings from murine models (Gallagher et al., 2012;Johns et al., 2017). This variability in gene expression, also known as noise, has been The percentage of the variance of each principal component (PC1 and PC2) for the control (C1-C4) and LLC (L1-L3, L5-L6) samples. Sample L4 did not pass the quality filters, and it was removed from the analysis (B) Heatmap of 18 Z-score normalized differentially expressed miRNAs of control (C1-C4) and LLC (L1-L3, L5-L6) samples by unsupervised hierarchical clustering analysis. Down-regulated and up-regulated miRNAs with absolute values of fold-change > 1.5 and FDR < 0.05 (Wald statistics) are shown in red and blue, respectively. (C) Scatterplot comparing abundance (Counts per Million, CPM; x-axis) and their degree of expression (log2-fold change, y-axis). Each dot represents differentially expressed miRNAs (fold-change > 1.5 and FDR < 0.05; Wald test), and red and blue dots (up-and down-regulated miRNAs, respectively) highlight potentially relevant miRNAs associated with muscle wasting. Gray diamonds represent the muscle-specific miRNAs (MyomiRs). Vertical lines represent thresholds for genes with low, medium, and high abundance, as defined by quartiles (Q1 and Q3).
described as an essential feature of any biological system (Eldar and Elowitz, 2010), and previous gene expression studies have been able to classify different states of diseases based on this variability (Mar et al., 2011;Ecker et al., 2015;Zhang et al., 2015;Guan et al., 2016).
To generate a complete picture of transcriptome content and dynamics, we explored differential expression data to identify transcriptional factor-binding motifs. Our motif analysis for up-regulated genes identified enrichment for several transcription factors, including FoxO, NF-κB, AP-1, Stat3, and Smad. These factors have already been described in the regulation of muscle atrophy by activating genes related to the proteolytic ubiquitin-proteasome, autophagylysosomal, metabolic adaptation, myogenesis, differentiation, and immune-modulation (Cai et al., 2004;Costelli et al., 2005;Choi et al., 2012;Gerstein et al., 2012;Zhang et al., 2013;Chen et al., 2017). Similarly, motif analysis also showed that downregulated genes were associated with myogenic transcription factors (Berkes et al., 2004;Watanabe et al., 2007;Rossi et al., 2016), metabolism (Chen et al., 2010), fiber transition (Grifone et al., 2004;Tsika et al., 2008), and muscle contraction (Pacini et al., 2013). Most importantly, our motif analysis revealed transcription factors that have not yet been reported or identified previously in cancer cachexia. These factors are related to cell cycle and myogenesis (E2f3, Yy1, and Creb3) (Asp et al., 2009;An et al., 2014;Zhou et al., 2015), unfolding protein response (Xbp1) (Jheng et al., 2018), and muscle fiber metabolism (ESRRA) (LaBarge et al., 2014). Also, we identified increased expression of the myogenic regulatory factors Myf6 and Myog. Interestingly, it has been demonstrated that the FIGURE 5 | Integrative analysis revealed a set of extracellular matrix mRNAs regulated by miRNAs. (A) The network generated is consisted of 171 interactions between 18 miRNAs and 131 target genes transcripts. Up-and down-regulated miRNAs (fold-change > 1.5 and FDR < 0.05; Wald test) are represented by blue and red nodes, respectively. Node size indicates the number of miRNA-target gene transcripts, and the gray edge width denotes overlapping miRNA-target gene transcripts measured by the Jaccard coefficient (JC). (B) Gene ontology analysis of the mRNAs predicted and validated as regulated by miRNAs. Each horizontal black bar represents the ontology term fold enrichment compared to the total number of genes in each term. Regulatory network displaying predicted (dashed lines) and validated (solid lines) interactions between the miRNAs (rectangle) and target mRNAs (circles), and physical and pathway protein-protein interactions (solid lines) for extracellular matrix (C), cell migration (D), and transcription factors (E). Up-and down-regulated miRNAs (fold-change > 1.5 and FDR < 0.05; Wald test) are represented by blue and red nodes, respectively. Gray nodes represent non-regulated genes.
over-expression of Myf6 inhibits the transcription of sarcomeric proteins by inhibiting Mef2 (Moretti et al., 2016). On the other hand, Myog induces the expression of the atrogenes Trim63 and Fbxo32 (Moresi et al., 2010). Together, the differential expression profile of transcriptional factors suggests that muscle atrophy in cancer cachexia can only be understood in the context of simultaneous signaling pathways activated by different transcription factors. This result is in contrast with a previous report indicating that Foxo is a single master regulator controlling muscle wasting during cancer cachexia (Judge et al., 2014). The combined action of transcriptional factors is supported by the long-standing view that specific combinations of transcriptional factors act cooperatively or in sequential steps in the regulation of gene expression (Reiter et al., 2017). Also, it has been recently demonstrated that the combination of the transcriptional factors NF-κB, SRF, and IRF controls gene expression through logical OR gates in macrophages (Cheng et al., 2017). Moreover, the cooperative interaction of the transcriptional factors STAT3 and NF-κB promotes muscle atrophy (Ma et al., 2017), further demonstrating cooperative actions of transcription factors controlling gene expression in muscle-wasting conditions. Even though transcriptional factors are essential to understanding gene regulation, it is also crucial to identify post-transcriptional regulation mediated by miRNAs, which induce mRNA decay and inhibit translation (Bartel, 2018). To better understand the miRNA-mRNA interactome in cancer cachexia, we combined the miRNA-mRNAs expression coprofiles in the same set of muscle samples. This analysis allowed us to identify biological processes such as ECM, cell migration, and transcription as regulated by miRNAs during muscle wasting in cancer cachexia. Within these interactions, we highlight the up-regulation of miR-29b-3p, which regulates a plethora of processes such as migration, ECM, and myogenesis. The miRNA miR-29b-3p has been described as a critical regulator in a variety of processes, such as myogenesis  and muscle atrophy Moraes et al., 2017), and it has a multiplicity of targets in skeletal muscle cells, including collagen transcripts (Zhou et al., 2012b), YY1 (Zhou et al., 2012a), p85 (Park et al., 2009), and IGF-1 (Gao et al., 2016).
It is worth noting that the same cancer cachexia mouse model as ours has been previously used to identify the global miRNA expression profile in wasting muscles , but only miR-223-3p was shared as deregulated with our data. Additionally, our set of deregulated miRNAs does not overlap with the muscle miRNA expression profile from other mice models of cancer cachexia or cachectic patients (Soares et al., 2014;Narasimhan et al., 2017). These discrepancies may be a consequence of the highly dynamic post-transcriptional regulation mechanisms exerted by miRNAs, which may act to buffer fluctuations in gene expression and confer robustness in signaling outcomes for specific regulatory networks (Ebert and Sharp, 2012). Another plausible explanation is that different miRNAs sharing the same seed regions target the same gene transcripts or pathways during the development of muscle wasting (Ebert and Sharp, 2012). Consequently, target contextual features determine miRNA target recognition and regulatory outcome as well as RNA interaction networks in specific biological contexts (Carroll et al., 2014), reinforcing the importance of the use of parallel expression profiles of mRNA and miRNA, as used herein, to gain further insights into post-transcriptional controls underlying muscle wasting in cancer cachexia.
Collectively, our integrated mRNA and miRNA data pointed out the remodeling of the ECM components in the wasting muscles. These data are in line with previous studies showing down-regulation of ECM genes in cachectic muscles of tumorbearing mice and cancer patients (Gallagher et al., 2012;Judge et al., 2014;Moraes et al., 2017;Talbert et al., 2019). Muscle fibrosis in cachectic patients based on histological analysis has also been described in cancer cachexia (Judge et al., 2018), but these data remain to be validated at the RNA and protein levels. Therefore, determining the molecular mechanisms that lead to remodeling of the ECM and its relation to atrophy in cachexia justify further investigations.
While our study identified transcriptional and posttranscriptional regulatory networks underlying muscle wasting in cancer cachexia, our in silico approach is limited in several aspects. First, most studies, if not all, focusing on rodent models, suffer from limitations. Although the LLC mouse model has been extensively useful in elucidating several mechanisms of tissue wasting, it does not fully recapitulate the phenotype of human cancers, either by not forming spontaneous tumors or because of their inability to reconstitute a tumor microenvironment (Talbert et al., 2019). Second, although our work describes the variability in the expression of specific sets of genes within LLC samples, further studies are needed to establish how interindividual variation in gene expression affects the prevalence and severity of the syndrome. Finally, our study reveals several transcription factors potentially involved in the regulation of cachexia genetic program, and experiments are still necessary to quantify their sequence-specific DNA-binding activity. In the same way, we predicted miRNA-target interactions that deserve experimental verification.

CONCLUSION
In conclusion, our integrative analysis of miRNA-mRNA coprofiles comprehensively characterized regulatory relationships of molecular pathways, including miRNAs targeting ECMassociated genes, which may play a role in the development of muscle atrophy in cancer cachexia. The in silico analyses of the transcriptome data also revealed that these ECM genes are potentially regulated post-transcriptionally by miRNAs, such as miR-29a-3p, and transcriptionally by the NF-κB and Ciita transcription factors.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the Gene Expression Omnibus (GEO) DataSets (https://www.ncbi.nlm.nih.gov/gds) under the accession numbers GSE144567 (mRNAs) and GSE145393 (miRNAs).

ETHICS STATEMENT
The animal study was reviewed and approved by Institute of Bioscience of Botucatu Ethics Committee on Animal Use, from the São Paulo State University (UNESP, Brazil).

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene. 2020.00541/full#supplementary-material FIGURE S1 | (a) The exposed carcass of control and Lewis Lung Cancer (LLC) tumor-bearing mice, twenty-two days after subcutaneous injection of PBS or 1.5 × 10 6 LLC cells, respectively. (b) Tumor mass of LLC tumor-bearing mice. (c) Splenomegaly in LCC tumor-bearing mouse compared to a control.
FIGURE S2 | mRNA levels of the genes involved in myogenesis, sarcomere, proteasome, and ECM in LLC and control groups. RT-qPCR data are presented as Log2 fold change (2 − Ct ) relative to Rpl13a. Statistical difference was analyzed by Student's t-test. * P < 0.05.

FIGURE S3 | (a)
Gene-Ontology analysis of differentially expressed genes (DEG) from control vs. tumor-bearing mice (LLC) for individual clusters (I to IV) identified by the unsupervised hierarchical clustering analysis (shown in Figure 2B). Each horizontal gray bar represents the gene fold enrichment compared to the total number of genes in each ontology term. control (C1-C4) analyzed by unsupervised hierarchical clustering. Down-and up-regulated genes with absolute values of fold-change > 1.5 and FDR < 0.05 (Wald test) are shown in red and blue, respectively. Subgroup A-cumulative frequency distribution of the DEG (log2-fold change, x-axis) from LLC (L1-L6) vs. control samples, indicated as a percentage (%, y-axis) for up-and down-regulated genes in (e) and (f), respectively. (g) Gene-Ontology analysis of DEG from control vs. LLC samples from subgroup A. Each horizontal black bar represents the gene fold enrichment compared to the total number of genes in each ontology term.
FIGURE S4 | mRNA levels of the genes involved with ECM in different muscle types from LLC and control groups. RT-qPCR data are presented as Log2 fold change (2 − Ct ) relative to Rpl13a. Statistical difference was analyzed by Student's t-test. * P < 0.05.