Relevance of cancer cachexia models – muscle whole genome gene expression in human and animal cachexia

Background Cancer cachexia is a complex and multi-factorial syndrome. As currently available therapeutic options are limited, more in-depth knowledge on cachexia pathophysiology and the underlying molecular mechanisms remains warranted. Studies with animal models provide useful insights but they only mimic the human situation to a certain degree. Furthermore, there is heterogeneity in the design of published animal studies and outcomes. To further address this issue, we performed a comparative study analysing muscle whole genome gene expression of different cachexia studies in mice and human. Methods We selected data sets from the NCBI Gene Expression Omnibus database containing muscle gene expression data measured by micro-array or RNA-sequencing, at least comprising a cachectic/tumour bearing group (n>3) and a non-cachectic/control group (n>3). This provided 12 datasets; 9 from mouse models and 3 human datasets. All datasets were quality checked, normalised and annotated. Datasets were merged and compared at different levels. General similarity and differences in gene expression were determined using ordered list analysis and principal component analysis (PCA). Moreover, similarities and differences at pathway level were studied by applying gene set enrichment analysis (GSEA) of KEGG pathways. Results Animal models displayed similarities to each other and to human datasets at different levels and with different processes. At the gene level, a similarity analysis indicated little similarity between the animal models and the human datasets, while animal models showed high similarity. Only one of the C26 mice models (GSE121972) showed significant similarity to more than one human dataset. Moreover, one human dataset comparing cachectic and non-cachectic cancer patients showed no similarity to any of the other datasets. PCA results indicated that a xenograft model showed most different expression from the other datasets and the Lewis lung carcinoma model to be least different from the human datasets. GSEA results showed four pathways clearly standing out across experiments with downregulation of oxidative phosphorylation and thermogenesis pathway, and upregulation of the proteasome and RNA transport pathway. However, these pathways were not consistently changed in the human datasets. Conclusions Our comparative analysis showed that there is currently no basis to define a preferred animal model for human cachexia. More human datasets containing proper controls are needed. Repetition of the current analysis upon publication of additional human datasets is warranted.


Introduction
Cachexia is a multifactorial syndrome characterised by disease-induced progressive loss of muscle and/or fat mass. It is often present in patients with chronic kidney disease, heart failure, COPD or cancer. In many patients cachexia contributes to reduced treatment efficacy and quality of life. Experts generally agree that cachexia should be treated with a multimodal approach addressing dietary intake, systemic inflammation and physical activity [1,2]. These different treatment components are often studied using different mouse models. The current study aimed to evaluate the relevance of different animal models mimicking human cancer cachexia by comparing muscle gene expression in different mouse models with that in human cancer cachexia. To this end, we used a big data approach by analysing publicly available datasets present in the GEO database containing whole genome gene expression measured by either microarray or sequencing.
Different types of tumours can cause cachexia both directly, by secreting cachexia-inducing factors, and indirectly, by triggering a systemic inflammatory process. To simulate human cancer cachexia and investigate treatment options, different types of animal models exist [3][4][5][6][7][8][9][10]. The traditional and most common models, i.e. C26 or Lewis Lung Carcinoma (LLC), use direct injection of in vitro cultured tumour cells [3,4,[6][7][8]10]. These are usually injected in the flank or muscle of the mouse where they grow into a solid tumour. More recent are models with genetic modifications causing spontaneous tumours [5], models using adenoviruses to induce tumour development [9], and xenograft models where tumour cells are harvested from human tumours and implanted in mice [8]. The difference between the models not only lies in the type of tumour induction. Different models have different time-spans ranging from 14 days to 8.5 weeks after tumour induction. In addition to the different types of models, experimental conditions may also vary between experiments using the same model. Despite these differences, all these models aim to mimic human cancer cachexia.
To assess animal model relevance, meta-analyses and systematic reviews have been performed on differences and commonalities of molecular processes between mice and men [11,12]. However, no studies directly compared data of molecular processes driving the cachexia in animals or humans, which makes the relevance of the animal models used a subject for debate. Therefore, our aim was to elucidate possible differences and commonalities in molecular processes in muscle tissues from different cachexia samples.
Whole genome gene expression of muscle tissue has been measured in several cachexia studies using micro array analysis or RNA sequencing. Upon publication of these whole genome gene expression studies, raw data is often made available via the NCBI Gene Expression Omnibus database. By comparing the different datasets using a big data approach, differences and similarities between models and experimental conditions can be examined. Moreover, possible knowledge can be gathered on underlying processes causing muscle wasting and on possible differences between animal models for cancer cachexia and human cancer cachexia.

Tumour model Dataset inclusion
We selected data series present in the NCBI Gene Expression Omnibus database using the keywords cachexia or (cancer AND skeletal muscle). From these, we selected those series containing muscle gene expression data of a cachectic/tumour bearing group (n>3) and a non-cachectic/control group (n>3). This provided 12 data series; 9 from mouse models and 3 human datasets (figure 1 and table 1).

Data processing
All datasets were quality checked, normalised and annotated with the latest annotations available using R statistical computing software (https://www.r-project.org/). The flowchart in figure 1 summarizes data selection and processing work flow. Three different measurement modalities were distinguished; (1) Affymetrix arrays, (2) Illumina bead-chip arrays and (3) Illumina HiSeq RNA-sequencing results. (1) Data obtained with Affymetrix arrays were quality checked and normalized using the robust multi-array analysis (RMA) algorithm [13] as implemented in the Bioconductor package AffyPLM. Probe sets were identified with genome information according to Dai et al. [14] based on annotations provided by the Entrez Gene database (custom CDF v23). (2) Data collected with Illumina beadchip arrays were background corrected and normalized using quantile normalization (neqc) [15]. Subsequently, unspecific probes were removed [16] (3) Illumina HiSeq RNA-sequencing results were filtered for genes having an expression level greater than 10 counts.
Next to that, the library size and the experimental design was taken into account and corrected for, using the Bioconductor package edgeR [17,18]. Subsequently, library size differences were adjusted using the trimmed mean of M-values normalization method [19], implemented in the Bioconductor package edgeR [18]. Counts were then logtransformed and the observed mean-variance trend was converted into precision weights by the voom function [20] in the Bioconductor package limma [21].
After normalization, all different data types were analysed similarly. Differential expression of probe sets (genes) was determined using linear models (package limma) and an intensity-based moderated t-statistic [21,22]. For each data set, samples from tumour bearing animals/patients were compared to a healthy control group, with one exception: in GSE85017 a comparison between cachectic and non-cachectic cancer patients was made, because this dataset did not contain a non-cancer control group. Genes with a p-value of p < 0.01 were considered to be significantly

Data analysis
All data analysis was performed in R. Data visualisation was done using ggplot2 and corrplot [23,24]. For similarity analysis we used the Bioconductor R library OrderedList [25]. Here, we used ranked listed based on the t-test statistic to compare all experiments in a pairwise fashion. Similarity scores with p<0.05 were normalized and visualized in an association matrix using corrplot and in a network plot using Cytoscape [26]. To assess differences between experiments, a sparse principal component analysis (sPCA) was performed available in the mixOmics package [27,28]. Using this sPCA, we reduced the dimensions and identified the top genes responsible for variation between the different datasets. Based on the number of different experimental types we selected the number of components to be (n-1) and we specified to keep the 100 genes most responsible for the variation within each component. Changes in individual genes were related to changes in pathways by gene set enrichment analysis (GSEA) [29] using the subset of metabolic and signalling pathways retrieved from the expert-curated Kyoto Encyclopedia of Genes and Genomes (KEGG) database [30]. For each comparison, genes were ranked on their t-value that was calculated by the empirical Bayes method. Statistical significance of GSEA results was determined using 10,000 permutations. GSEA and visualization was performed using the Bioconductor package clusterProfiler [31]. To assess variation between separate samples of all datasets, an sPCA with 6 components of each 25 genes was performed on gene expression, normalized per experiment, of individual samples of all experiments (n=168).

Included datasets
We identified data sets from five different animal models that met our criteria for inclusion: two inoculation models; Next to this, we included 3 human datasets, one of which containing whole genome gene expression data of m. rectus abdominis samples of cachectic or non-cachectic pancreatic and colorectal cancer patients, one containing quadriceps samples of upper gastrointestinal cancer (UGI) patients and healthy controls, and one containing m. rectus abdominis samples of UGI patients and weight stable (WS) patients undergoing surgery for benign, non-inflammatory conditions.

Number of significantly changed gene transcription activities
Analysis of the number of significantly changed gene transcription activities in each experiment showed that there were large differences in the magnitude of tumour effects on the muscle (figure 2). The C26 model was found to be overall the most invasive model affecting expression of up to 56% (GSE56555) of all genes. However, this model also showed considerable between-study variation, as in the apparently least invasive study (GSE63032) transcription of only 16% of all genes was changed. The other mouse models showed on average a lower number of changed gene transcription activities. The adenovirus-induced non-small cell lung cancer (NSCLC, GSE107470) affected 31%, the Lewis Lung Carcinoma (LLC, GSE114820) affected 29%, the human xenograft model (GSE80081) affected 16%, and the spontaneous pancreatic cancer (GSE81931) affected transcription of 12% of all genes. Compared to these mouse models, human datasets showed only very mild effects with 6%, 2% and 1% of all genes being affected in transcription (respectively GSE34111, GSE85017 and GSE18832).

Similarity analysis
Subsequently, we created a similarity network and matrix (figure 3 and 4) based on our pair-wise similarity analysis.
Results showed that the highest similarity is among the C26 models and that only very minor similarities are found between the mouse models and the human models. Surprisingly, the similarity of two skeletal muscle samples from the

Sparse Principal Component Analysis Expression Differences
To assess differences between datasets, we performed a sparse Principal Component Analysis with 6 components of

Gene Set Enrichment Analysis
When looking at the KEGG pathway enrichment, four pathways seemed to be important in all experiments; the proteasome and RNA transport pathway were both strongly upregulated in most experiments, while the oxidative phosphorylation and thermogenesis pathway were strongly downregulated in most experiments ( figure 6).
Interestingly, the effects on these pathways were quite uniform along all animal models, but differed to some extent between the human experiments with the Proteasome and RNA transport pathway being downregulated in the human m. quadriceps (GSE34111). When looking at the change in expression of the genes in these separate pathways ( figure   7), it became clear that different models had differential effects on these four important pathways. In the most

Discussion
Different models are in use to evaluate possible therapeutic strategies in cachexia, like those based on nutritional, exercise and pharmacological intervention. The present study aimed to determine the relevance of mouse models used to mimic human cachexia, by comparing their effects on muscle gene expression, both between models and with data from human cancer cachexia studies. Publicly available datasets like those present in the GEO database are increasingly offering opportunities for such meta-analyses without the need to perform additional (animal) studies.
Overall, our analysis revealed some similarities between animal models and, to lesser extent, to human datasets, in particular with respect to biological process level.
On gene level, analysis indicated only little similarity between the data from animal models and those from human datasets, while more similarities were found between the animal models. Only the m. extensor digitorum longus and m. soleus data from one of the C26 models (GSE121972) showed significant similarity to more than one human dataset. Moreover, the human dataset without healthy controls showed no similarity to any of the other datasets.
However, principal component analysis to examine the specific differences between datasets revealed that the xenograft model showed the most different expression pattern compared to those of the other datasets. Here, the LLC model seems the least different from the human datasets based on clustering in the sPCA heatmap ( figure 5.5 D). These results give a good impression of general similarities and differences between data obtained from animal models and those from human cachectic patients. However, for analysis of relevant molecular processes, GSEA results are preferred.
With GSEA, four pathways were found to be clearly standing out. A clear downregulation of oxidative phosphorylation was observed in all animal models, which is in line with literature data [32]. However, this downregulation was only seen in the human m. quadriceps, while in the other human datasets, oxidative phosphorylation was not significantly enriched. For the thermogenesis pathway, a clear downregulation was seen in all animal models, whereas also here only the human m. quadriceps data showed significant downregulation. The downregulation in muscle thermogenesis is surprising since literature suggests an important role for increased thermogenesis during cachexia, not only in brown adipose tissue [33] but also in muscle [34]. When considering upregulated pathways, the proteasome pathway appeared to be strongly upregulated which is in line with literature [35]. However, the term proteasome is most frequently used in combination with the ubiquitin system [12,36,37]. In our analysis, we do not see the ubiquitin mediate proteolysis to be significantly enriched in any of the datasets. This is in line with previous suggestions that the ubiquitin-proteasome system is not as important in cancer-induced cachexia as previously thought [38]. Another relevant pathway is the strongly upregulated RNA transport pathway which is the only pathway where human datasets cluster together with some animal models (NSCLC and the LLC model). In this pathway, Eukaryotic translation initiation factor 4E binding protein 1 [EIF4EBP1] seems to be an important factor being highly downregulated (figure 7 B).
Interestingly, this gene is mostly referred to in the context of its influence on the mTOR and AKT pathways in cachexia, whereas we obtained no indication of importance of either the mTOR nor the AKT pathway in our analysis (supplemental figure 1) [39]. Finally, some pathways came out of the human datasets without apparently being relevant in the animal models. For example, the RNA degradation was upregulated in both human m. rectus abdominis samples, while not significantly found to be enriched in any of the animal models (supplemental figure 1).
In line with literature, we see that the measurement platform markedly influenced clustering based on individual gene expression of separate samples [40]. This occurred despite the fact that we used the mixOmics analysis package correcting for differences in platform and experiment [27,28]. However, we do not see platforms clustering together in analyses on differential expression on gene or pathway level making these comparisons still highly relevant.
Another limitation of this analysis is that the merge of all data from different platforms causes some data loss. This is the first study on muscle gene expression in cachexia directly comparing whole genome gene expression data from different experiments. This type of analysis is getting more important due to the large increase in raw data made available and is already used for analysis of tumour gene expression in different types of cancer [41][42][43][44]. However, in case of muscle gene expression in cancer, the number of available data sets are still limited, making a repetition of this analysis upon publication of new data highly relevant. Unfortunately, phenotypic (meta-) data was often not present or insufficient with the datasets, making integration with other data like body weight loss and gene expression not possible for this analysis. Moreover, integration of different omics data like protein expression might be highly relevant to also gain insight into post-transcriptional processes important in cancer cachexia [45].
In conclusion, this study allows to do some suggestions, based on the currently available data, on the relevance of different models mimicking human cancer cachexia. Most importantly, we do not see one model out-performing other models. Each model shows its own differences and similarities to the publicly available human datasets. Thus, the

Declarations
Ethics approval and consent to participate: not applicable.

Consent for publication : not applicable
Availability of data and material: R scripts have been added to make the analysis used available Author contributions: RP, GH, RW and KvN did the concept and design of the study. Acquisition, analysis and interpretation of the data were carried out by RP, GH and KvN. Drafting of the manuscript was performed by RP, GH, RW and KvN. All authors revised the manuscript critically for important intellectual content and approved the final manuscript.

Conflict of interest:
The authors declare no competing interests.
Funding: research was funded by the Wageningen University