Differentially expressed alternatively spliced genes in skeletal muscle from cancer patients with cachexia.

BACKGROUND: Alternative splicing (AS) is a post-transcriptional gene regulatory mechanism that contributes to proteome diversity. Aberrant splicing mechanisms contribute to various cancers and muscle-related conditions such as Duchenne muscular dystrophy. However, dysregulation of AS in cancer cachexia (CC) remains unexplored. Our objectives were (i) to profile alternatively spliced genes (ASGs) on a genome-wide scale and (ii) to identify differentially expressed alternatively spliced genes (DASGs) associated with CC.
METHODS: Rectus abdominis muscle biopsies obtained from cancer patients were stratified into cachectic cases (n = 21, classified based on International consensus diagnostic framework for CC) and non-cachectic controls (n = 19, weight stable cancer patients). Human transcriptome array 2.0 was used for profiling ASGs using the total RNA isolated from muscle biopsies. Representative DASG signatures were validated using semi-quantitative RT-PCR.
RESULTS: We identified 8960 ASGs, of which 922 DASGs (772 up-regulated and 150 down-regulated) were identified at ≥1.4 fold-change and P < 0.05. Representative DASGs validated by semi-quantitative RT-PCR confirmed the primary findings from the human transcriptome arrays. Identified DASGs were associated with myogenesis, adipogenesis, protein ubiquitination, and inflammation. Up to 10% of the DASGs exhibited cassette exon (exon included or skipped) as a predominant form of AS event. We also observed other forms of AS events such as intron retention, alternate promoters.
CONCLUSIONS: Overall, we have, for the first time, conducted global profiling of muscle tissue to identify DASGs associated with CC. The mechanistic roles of the identified DASGs in CC pathophysiology using model systems is warranted, as well as replication of findings in independent cohorts.

Introduction

Alternative splicing (AS, also used synonymously with the term ‘alternatively spliced’) is a crucial post‐transcriptional gene regulatory mechanism that is involved in generating gene isoforms from a single precursor mRNA, thereby increasing proteome diversity.1, 2 More than 90% of the human genes are alternatively spliced, but such complexity is not observed in lower organisms.3 AS is regulated in a tissue‐specific manner, with skeletal muscle exhibiting the highest number of alternatively expressed exons,4, 5 and is recognized to contribute to a wide range of physiological and cellular processes.6, 7 Aberrant splicing mechanisms due to splice site mutations also contribute to tumourigenesis.8, 9 Aberrant splicing mechanisms were shown to be associated with various muscle‐related pathologies such as Duchenne muscular dystrophy.10 However, the contribution of AS dysregulation in human cancer cachexia (CC) is unknown.

Cancer cachexia, a debilitating condition seen in advanced cancer patients is associated with involuntary weight loss and loss of lean body mass with or without loss of fat mass.11 CC can be seen as a consequence of complex host‐tumour interactions eventually leading to a state of energy imbalance.12, 13 Degeneration of skeletal muscle and impaired myogenesis are prominent features in CC.14, 15 In vitro models suggested that AS plays a major role in myogenic differentiation in a spatiotemporal manner.16 Dysregulation of splicing mechanisms is observed in polygenic muscular disease such as sporadic inclusion body myositis.17 It is therefore very likely that AS dysregulation contribute to CC pathophysiology.

Traditional gene expression arrays (at transcriptional level) have enabled researchers to identify many differentially expressed genes to understand the underlying biology in a disease‐specific context, especially in cancers.18 However, it is now possible to address gene regulatory mechanisms at a finer level (post‐transcriptional) with the availability of global microarrays and massive parallel sequencing technologies.19 In vitro experiments have shown that isoform‐specific expressions delineate tumour associated signatures from non‐tumour signatures more accurately.20 Given its importance in a disease context, it is imperative to understand the effect of AS in CC. Here, we chose to profile muscle tissue as skeletal muscle atrophy is a characteristic feature of CC. Recognizing the plasticity exhibited by the muscle tissue and the dynamic mechanism exhibited by AS, identifying tissue‐specific isoforms may shed more light on CC pathophysiology, which remains as a gap in the literature. The study design is cross‐sectional, and the aims of the current sub‐study are (i) to profile alternatively spliced genes (ASGs) in CC on a genome‐wide scale from human skeletal muscle biopsies and (ii) to identify differentially expressed alternatively spliced genes (DASGs) associated with CC. In this current study, we have shown that DASGs are associated with CC pathophysiology. The identified DASGs may potentially play a role in myogenesis, inflammation, and ubiquitination, which are classically associated with skeletal muscle wasting.

Materials and methods

Procurement of muscle biopsies

Rectus abdominis muscle biopsies were obtained from University of Calgary Hepatopancreaticobiliary/Gastrointestinal Tumor Bank from pancreatic cancer and colorectal cancer patients with liver metastasis who underwent laparotomy at the Foothills Hospital from 2006 to 2013. Tumour stage is reported according to American Joint Committee on Cancer v7. Standard procedures were followed for tissue procurement and storage. Specifically, biopsies of rectus abdominis muscle were taken by the operating surgeon within 30 min of the start of the surgery using sharp dissection, immediately flash frozen in liquid nitrogen to minimize ischaemic shock post‐devitalization, and stored at −80°C until further use. Written informed consent was obtained from all study participants, which was approved by Conjoint Health Research Ethics Board at the University of Calgary (Ethics ID E‐17213). Health Research Ethics Board of Alberta (HREBA)–Cancer Committee approved the current study protocol for transcriptome profiling and access to the patient's clinical information (protocol number ETH‐21709).

The stratification of patients was based on the International consensus diagnostic framework for CC.11 Patients were classified into cachectic cases, defined as those with either (i) >5% pre‐illness weight loss (WL) over a period of 6 months, (ii) >2% WL with either body mass index (BMI) <20, (iii) or sarcopenia (defined by skeletal muscle index (SMI) cut‐points using computed tomography, CT) with >2% WL. Non‐cachectic controls are defined as those who were weight stable (WS) cancer patients over a period of 6 months, compared with their pre‐illness weight. Physician documented weight loss at first presentation of the patient in the clinic was used for the study. We excluded patients with no clinical chart, or no recorded WL information, who were below 18 years of age and were unable to provide written consent.

Body composition analysis using computed tomography

Of the 42 patients in the study, 34 patients had a CT prior to surgery (71.51 ± 45.9 days). CT‐based body composition analysis was carried out using lumbar vertebrae (L3) as a standard landmark, as described elsewhere.21 Cross‐sectional muscle area (cm2), SMI [cross‐sectional muscle area normalized to their stature (cm2/m2)], total adipose tissue, and muscle radiation attenuation (MA) were measured. Muscle radiation attenuation was measured in Hounsfield units, and the ranges for these measurements are described elsewhere.21 Patients were classified sarcopenic based on the previously described SMI values.21, 22

RNA isolation and integrity of RNA

Total RNA was isolated using TRIzol and QIAGEN RNeasy maxi kit (Mississauga, ON, Canada); 260/280 ratio was measured using Nanodrop and RNA integrity number was assessed using Agilent Bio‐analyzer 2100 for all the samples.

Human transcriptome array 2.0, hybridization methods, and data analysis

The human transcriptome (HTA) array 2.0 has been designed to capture all known and putative coding exons and coding gene exon‐intron boundaries in the human genome. Ten probes span each exon, and four probes span each exon‐intron boundaries (splice junctions). The array also contained non‐protein coding transcripts, but these were not mined in the current study. Instead, we focused primarily on the protein coding gene transcripts (isoforms and exon level). The entire protocol was carried out as per the manufacturer's instructions (http://www.affymetrix.com). Briefly, 100 ng of total RNA was used as a starting material for labelling and hybridization. Hybridization was performed for 16 h in Genechip hybridization oven 645 using the standard procedures (http://www.affymetrix.com). The washing and staining protocol was carried out using Genechip Fluidics Station 450 according to manufacturer's protocol. The HTA 2.0 arrays were scanned using Affymetrix GCS 3000 7G scanner to generate the raw intensity CEL files for downstream analyses. Identification of ASGs and Gene Set Enrichment Analysis (GSEA) were carried out using Partek Genomics Suite 6.6 (PGS 6.6).

Analysis for alternatively spliced genes

Exon level intensity estimates were generated using RMA method, which includes background correction, quantile normalization, and log2 transformation. Exons that were not expressed in all of the samples were excluded from further analysis. In the AS analysis, PGS generates two results—at the transcript level and at the exon level. Initially, at a transcript level, differentially expressed ASGs were identified between cachectic cases and WS cancer patients at P‐value of <0.05 (PGS defined this as ‘alt‐splice P‐value’). At the exon level, the expression between cases and WS cancer patients with P < 0.05 and FC ≥ 1.4 were identified using one‐way analysis of variance, and these exons were mapped to their corresponding transcripts. Final representation of data would therefore reflect a composite signature of the ASG transcripts overlapped with exon level results (P < 0.05, FC ≥ 1.4). These overlapped transcripts were called DASGs. The identified DASGs were used for subsequent GSEA and for ingenuity pathway analysis (IPA). Select DASGs were validated using semi‐quantitative RT–PCR. The raw files and normalized counts have been submitted to Gene Expression Omnibus database (GEO accession ID‐GSE85017).

Validation of representative differentially expressed alternatively spliced genes using semi‐quantitative RT–PCR

A total of 1 μg RNA was converted into cDNA using high‐capacity cDNA Reverse Transcription Kit (Applied Biosystems, ON, Canada) using the manufacturer's protocol. Reverse transcription was performed using the following thermal cycler conditions: 37°C for 60 min and 95°C for 5 min.

A total of 100 ng cDNA was used for semi‐quantitative RT–PCR. Go Taq G2 Hotstart green Mastermix (Promega, Madison, Wisconsin, USA) was used. Primers were designed using primer3 software (v 0.4.0) (http://bioinfo.ut.ee/primer3‐0.4.0/) and oligocalc (http://biotools.nubic.northwestern.edu/OligoCalc.html). Six DASGs (IFRD1, KCNQ5, DEPDC1, UBA3, FNDC1, and CNNM3) were validated in representative cases and WS cancer patients. The amplified products were separated using 2.5–3% agarose gel and stained using Redsafe stain solution (Intron Biotechnology) for visualization. Densitometric scans were used to quantify gel bands using image j software and to calculate the ratios between the splice variants.23

Validation of AS events can also be carried out using other methods. Use of boundary spanning primer captures the relative abundance of transcript with high specificity.24 An alternative method has also been established using RT‐qPCR, which can differentiate smaller expression differences between two transcripts.25 These aforementioned methods capture the exponential amplification phase of the detected signals. We used semi‐quantitative PCR methods to explain the relative abundance of transcript and quantified the relative abundance of the expressed transcript by utilizing gel electrophoresis with imaging software, as described by others.26

Functional annotation of the identified differentially expressed alternatively spliced genes and identification of canonical pathways

Gene Set Enrichment Analysis was carried out using PGS 6.6 to understand the potential functions of the identified DASGs in CC pathophysiology. IPA was used for identifying the canonical pathways and upstream regulators for the DASGs.

Statistical analyses

Patient demographics data were represented as mean ± standard deviation. Independent t‐test and chi‐squared test were used for continuous and categorical variables, as appropriate. For the AS analysis, one‐way analysis of variance test was used to identify DASGs. To understand the association between DASGs and clinical factors, Pearson correlation test was carried out. For all the analyses, P < 0.05 was considered to be statistically significant.

Results

Patient demographics

In all, 42 age‐matched patients with pancreatic cancer and colorectal cancer with liver metastasis were selected for the study. Of these 42 patients, 22 were cachectic cases (hereafter referred to as cases), and 20 were non‐cachectic controls (hereafter referred to as WS cancer patients). Eight patients (four cachectic cases and four WS cancer patients) had completed a course of neo‐adjuvant chemotherapy but had not received any chemotherapy 4 weeks prior to surgery. The remaining study participants did not receive any chemotherapy before surgery.

The 260/280 ratio for all the samples ranged from 1.6 to 1.8. The RNA integrity number values for samples were between 5.9 and 8.9. Two samples had poor single cRNA yield (one of the intermediate steps leading to hybridization) and was not processed further in study, leaving 40 samples for further analysis; 21 patients belonged to cases and 19 patients belonged to WS cancer patients.

No significant difference was observed in age, gender, tumour stage, and tumour type, whereas BMI was significantly different between cases and WS cancer patients (n = 40, Table 1A). Results from the body composition analysis are represented in Table 1B. SMI, z‐score, and MA were found to be significant between cases and WS cancer patients. Abnormally low MA has been observed in pathological conditions such as cancer, where there is an excess infiltration of fat in muscle against the normally observed levels.27

Table 1A

Patient demographics

Characteristic	Cachectic cases (n = 21)	Non‐cachectic controls (n = 19)	P‐value
Weight loss (% mean)	11.4 ± 6.5	—	—
Age (mean, in years) a [range]	65.7 ± 10.5 [39–84]	64.2 ± 8.1 [46–77]	0.67
Tumour typeb
Pancreatic	12	12	0.69
Colorectal	9	7
Gender b			0.55
Male	8	9
Female	13	10
Body mass indexa (mean, in kg/m2) [Range]	24.2 ± 3.6 [19–29]	26.9 ± 3.9 [21–40]	0.02
Tumour Stage c			0.64
I	2	1
II	3	3
III	2	0
IV	14	15

Unpaired t‐test.

Chi‐square test.

Fisher's exact test.

Statistical analyses were carried out using SPSS v16, and the data were represented as mean ± standard deviation. P < 0.05 was considered statistically significant.

Table 1B

CT derived body composition analysis

Characteristics	Cachectic cases (n = 19)	Non‐cachectic controls (n = 15)	P‐value
Cross sectional skeletal muscle area (cm2)a
Male	137.8 ± 15.7	158 ± 12.9	0.17
Female	96.2 ± 14.6	103.5 ± 14.1
Skeletal muscle index (cm2/m2)a
Male	44.6 ± 5.9	49.1 ± 3.1	0.05
Female	36.3 ± 5.6	41.52 ± 7.43
Sarcopeniab
Yes	16	8	0.06
No	3	7
Z‐scorec	−0.75 ± 0.7	−0.09 ± 0.9	0.04
Total adipose tissue a
Male	215.6 ± 84.5	266.29 ± 77.2	0.9
Female	328.8 ± 126.4	302.2 ± 122.7
Muscle attenuation (HU)a
Male	32.9 ± 8.6	39.8 ± 6.9	0.02
Female	29.5 ± 7.3	36.4 ± 8.6

Unpaired t‐test.

Fisher's exact test.

Z‐score indicates how many standard deviations an element is from the mean value for cancer patients of the same age and sex.

CT quantification was carried out for patients who had CT prior to surgery. Skeletal muscle index, z‐score, and muscle attenuation were found to be significant between cachectic cases and non‐cachectic controls. Sarcopenia status based on Martin et al.21 classification was trending towards significance (P = 0.06).

Identification of differentially expressed alternatively spliced genes associated with cancer cachexia

The HTA array data identified a total of 8960 ASGs. At ≥1.4‐fold change and P < 0.05, 922 DASGs were identified, of which 772 DASGs were up‐regulated (Table ), and 150 were down‐regulated (Table ).

Representative DASGs from HTA array were selected (Table 2A) for semi‐quantitative RT–PCR. The thermal cycle profiles and the primers sequences used are given in Table . The primer sequence for β‐actin (used as an internal control) was available in literature.28 Exons of four up‐regulated DASGs (IFRD1, KCNQ5, DEPDC1, and UBA3) and two down‐regulated DASGs (FNDC1 and CNNM3) were validated (Figure 1 and Tables 2A and 2B). IFRD1 and KCNQ5 were found to be associated with skeletal muscle differentiation (see GSEA results later). However, the roles of DEPDC1 and FNDC1 in skeletal muscle or CC are unknown. CNNM3 is vital for magnesium transport, and magnesium is known to regulate muscle contraction.29 From the HTA 2.0 results, exon 3 of UBA3 (ENST00000415609) was found to be up‐regulated in cases. However, designing the primer for this exon (ENST00000415609) was not feasible because of the presence of an upstream cassette exon event in another UBA3 transcript (ENST00000361055). Because both these transcripts were identified in muscle tissues, we chose to validate a cassette exon event (reported in UCSC genome browser), present in exon 4 (ENST00000361055).

Table 2A

DASGs validated using semi‐quantitative PCR and detected from HTA 2.0 array

Gene	Gene ID	Probe set ID	AS event	Fold change	P‐value
IFRD1	NM_001197079	PSR07011131	Exon 5 is a CE (up)	1.76	0.007
KCNQ5	NM_001160132	PSR06008716	Exon 9 is a CE (up)	1.51	0.005
DEPDC1	NM_001114120	PSR01043217	Exon 8 is a CE (up)	1.55	0.0002
FNDC1	ENST00000297267	PSR06013583	Exon 10 is a CE (down)	−1.32	0.03
CNNM3	ENST00000305510	PSR02009044	Exon 2 is a CE (down)	−1.58	0.01

AS event for all the validated DASGs were identified using UCSC genome browser (hg19). The direction of effect of the identified exons for DASGs are represented, with respect to the cachectic cases relative to non‐cachectic controls, as indicated by: ‘up’ means up‐regulated and ‘down’ means down‐regulated.

Figure 1

Semi‐quantitative RT–PCR validation of AS event. The forward and reverse primers are represented by arrows, which are designed from flanking constitutive exons. The five exons identified in human transcriptome array (indicated in orange) exhibited similar direction of effect when validated in representative cachectic cases and non‐cachectic controls. UBA3 (exon 4, ENST00000361055) also showed up‐regulation in semi‐quantitative RT–PCR. All six validated DASGs exhibited cassette exon property. NTC is a non‐template control and ‘M’ is a DNA ladder marker. β‐actin was used as an internal control.

Table 2B

Densitometry analysis for DASGs

Gene	Isoform inclusion ratio	Cachectic cases_mean	Non‐cachectic controls_mean	Fold change (densitometric analysis)
IFRD1	Incl Ex 5/ Skip Ex 5	6.41	5.5	1.16
KCNQ5	Incl ex 9/ Skip Ex 9	1.66	1.39	1.2
DEPDC1	Incl ex 8/ Skip ex 8	3.82	2.44	1.57
FNDC1	Incl Ex 10/Excl Ex 10	0.53	0.59	0.89
CNNM3	Incl Ex 2/Skip Ex 2	4.34	5.40	0.8
UBA3	Incl Ex 4/ Skip Ex 4	3.74	3.03	1.23

Representative DASGs with cassette exon event were validated using semi‐quantitative RT–PCR. Similar direction of effect was observed in both microarray and PCR experiments.

If an exon of a particular DASG is up‐regulated in cases, it means increased inclusion of an exon is observed in cases, relative to WS cancer patients. For example, exon 5 of IFRD1 was found to be up‐regulated in cases and was identified as a cassette exon (exon included or skipped). Based on semi‐quantitative RT–PCR, an increased inclusion of exon 5 was observed in cases, relative to WS cancer patients. Similarly, if an exon was down‐regulated in CC cases, an increased inclusion of an exon was observed in WS cancer patients, relative to CC cases. Exon inclusion rate was calculated based on the ratio of long isoform/short isoform.30 Exon inclusion rate (densitometric analysis of RT–PCR amplicons) was computed from the semi‐quantitative RT–PCR data, and these independent results were consistent with the findings from transcriptome array (Table 2B). All six validated DASGs exhibited cassette exon property, one of the most common AS events.

Up to 10% of the DASGs in the current study were identified to have cassette exons (exon present in some transcripts and absent in other transcripts). The next common splicing event observed in our study was alternate promoter usage (4.5% of identified DASGs had alternative transcription initiation sites, leading to isoform diversity). Other events identified were (i) variations in the sequence at the 3′ and 5′ positions in the splice donor‐acceptor sites termed alternate 3′ and 5′ transcripts, (ii) alternate finish (transcription ending at multiple sites and hence multiple isoforms), and (iii) intron retention. The definitions for each of the above AS events are adopted from the UCSC genome browser (https://genome.ucsc.edu/). Validation of AS events other than cassette exons are quite challenging, and an optimal way to validate diverse AS events is still emerging.

Differentially expressed alternatively spliced genes and its biological relevance to cancer cachexia using gene set enrichment analysis

The identified DASGs were associated with muscle structure and function (ENAH and ANKRD1), skeletal muscle differentiation (MEF2C, IFRD1, and KCNQ5), circadian rhythm (MKL1, QKI, and ARNTL), inflammation (ADIPOR1 and IL18), magnesium transport (CNNM3), protein ubiquitination (ANAPC1, UBB, and UBC), and signalling pathways (DEPDC1 and MAP2K5). A list of DASGs with their respective functions is summarized in Table 3.

Table 3

Functional annotation of DASGs

Function	Associated genes
Muscle
Muscle structure and function	DAAM1, PFN2, TRPM7, TPM3, MYL6B, PDE4D, SLMAP, ANKRD1, ENAH
Skeletal muscle cell differentiation	IGF1, LEMD3, MEF2C, PAXBP1,CYR61, IFRD1, FLRT3, KCNQ5, ACVR2A, ROCK2
Extracellular matrix protein	ADAM10, ADAM32, COL18A1, COL20A1
Lipid Biosynthesis	INPP4B , PIGX, PLA2G2A, PLA2G15, SYNJ1, FABP4, FAR1, GK, B4GALT4, ST8SIA5
Inflammation
Cytokine signalling and B cell activation	ADIPOR1, PTPN2, IL18, CLCF1, PPM1B, BCL6, PTPRC, TGFBR2, CD47
Protein ubiquitination and proteolysis	UBE4B, UBE2B, USP45, PSMA5, HERC4, HUWE1, UBR5, UBA3, UBE2G2, FBXO11, USPL1, UBA6, UBLCP1
Signalling pathways	DNM1L, MAP2K5, MKLN1, NCF2, DEPDC1, ENAH, FYB, PIK3R1, JAK2, KIDINS220, NEDD9
Transcription factors	INPP5F, CDCA7
Circadian rhythm	QKI, SIK1, NCOR1
mRNA‐splicing	CSTF3, SRSF10, SRSF3, SRSF4, SRSF5, SRSF6, TRA2A, CSTF1

Gene set enrichment analysis was carried out using PGS 6.6. Clusters with P < 0.05 and enrichment score of more than 1.2 were considered significant. Majority of the DASGs were associated with muscle structure and function, which might be directly associated with skeletal muscle dysregulation. Other functions include inflammation, energy homeostasis, and protein ubiquitination—dysregulation of all of these processes are well known to be associated with CC pathophysiology. The DASGs in bold are up‐regulated and the italicized genes are down‐regulated in the study.

Canonical pathways identified from ingenuity pathway analysis

A total of 83 canonical pathways were identified at P < 0.05 (Table ). Increased and decreased activity of a particular pathway is inferred based on the z‐score. If the z‐score is positive, then a pathway has an increased activity and vice versa for the decreased activity. The highly activated pathways include FLT3 signalling, IGF‐1 signalling, IL‐8 signalling, CNTF signalling, and CXCR4 signalling (z‐score range: 1.7–3.5). PTEN signalling was found to be decreased (z‐score −1.732). While IGF‐1 and IL‐8 have been studied in the context of CC, FLT3 signalling is an emerging new pathway with role in myogenic differentiation.31 Other significant pathways (P < 0.05) identified were protein ubiquitination pathway, glucocorticoid signalling, and IL‐4 signalling, findings that are consistent with CC literature on pathway‐based gene expression analysis.

Many of the identified DASGs have not previously been associated with CC. However, the pathways that the DASGs belonged to are known to be associated with CC.32, 33 The complete lists of canonical pathways are given in Table . While canonical pathways are informative, our interest was also to search for upstream modulators, which affect several of the identified downstream effectors. In this search, we identified TGFB1 as one of the upstream regulators with activation z‐scores of 3.1 with overlap P‐value of 0.004 (Figure 2). The overlap P‐value is an estimate of overlap between the DASGs identified in this study and the upstream regulator. Many of the up‐regulated DASGs were predicted to be activated than inhibited. A recent study had shown the effect of muscle dysfunction in CC through TGFB1 signalling. Activation of TGFB1 signalling in model systems led to an up‐regulation of NOX‐4 among other molecules, which eventually leads to defective muscle contraction.34 We also observed NOX‐4 to be up‐regulated in the human skeletal muscle tissues, in support of previous observations.

Figure 2

TGFB1 as an upstream regulator identified by IPA along with its downstream targets. Many of the up‐regulated DASGs were predicted to be (i) activated by TGFB1 (orange lines); (ii) inhibited by TGFB1 (blue lines); IPA identified molecules with no set predictions (grey lines) or those which could not be fit to a pattern of downstream molecules (yellow lines) are also illustrated. The DASGs highlighted in green are the down‐regulated, and the remaining DASGs are up‐regulated.

Pearson correlation analysis for DASGs with BMI and body composition measurements

The DASGs identified in GSEA (Table 3) were subjected to Pearson correlation test. DASGs identified with muscle structure and function, ubiquitination and inflammation were correlated to SMI; DASGs associated with lipid biosynthesis were correlated to muscle radiation attenuation. As CT was available for 34 patients, expression values from those samples were considered for correlations with body composition measurements. ENAH (r = −0.38, P = 0.03), KCNQ5 (r = −0.36, P = 0.04), and ROCK2 (r = −0.46, P = 0.006) were negatively correlated with SMI and MYL6B (r = 0.44, P = 0.009) was positively correlated to SMI. DASGs associated with inflammation such as BCL6 (r = −0.47, P = 0.004) and TGFBR2 (r = −0.35, P = 0.04) were negatively correlated to SMI, and ADIPOR1 (r = 0.34, P = 0.05) was positively correlated to SMI. INPP4B (r = 0.37, P = 0.04) was positively correlated with muscle radiation attenuation, and PLA2G2A (r = −0.47, P = 0.005) and B4GALT4 (r = −0.37, P = 0.03) were negatively correlated to muscle radiation attenuation. As BMI was significant between cachectic cases and non‐cachectic controls, we also correlated BMI with DASGs identified in GSEA (Table 3). MYOZ2 (r = −0.41, P = 0.01), TPM3 (r = −0.34, P = 0.04), and UBE2G2 (r = −0.34, P = 0.05) were negatively correlated to BMI.

Discussion

This is the first study to identify DASGs associated with CC. It is recognized that different isoforms generated from same pre‐mRNA are expressed in a tissue‐specific manner,35 and these isoforms are known to perform diverse functions.36 Therefore, an understanding of the isoform‐specific expression in muscle tissue may explain the hidden complexity hitherto not revealed at the conventional gene level studies. We identified several DASGs, which are involved in protein ubiquitination, skeletal muscle differentiation, and inflammation. These aforementioned mechanisms have been well documented for their role in CC pathophysiology. Circadian rhythm is another pathway that is slowly gaining prominence in CC pathophysiology,37 an indication of the pleotropic nature of the genes in diverse phenotypes.38 Many of the DASGs identified in this study either have not been previously reported in CC literature or were shown to be associated with other tissue types. For example, FNDC1 is associated with apoptosis of cardiomyocytes under hypoxic conditions.39 DEPDC1 gene with role in bladder cancer40 has not been reported to be associated with skeletal muscle functions. Independent confirmation of the described DASGs by semi‐quantitative RT–PCR gives us the confidence in the study findings that these are indeed expressed in skeletal muscle and may have a role in CC. It is currently known that >90% of all genes express multiple isoforms.3 With muscle having highest number of alternatively expressed exons4, 5 and skeletal muscle atrophy being a hall mark of CC, our study addresses a critical gap in literature by identifying the DASGs in CC pathophysiology using rectus abdominis muscle.

Many of the identified DASGs were associated with muscle structure and development, lipid biosynthesis, extracellular matrix, inflammation, and protein ubiquitination (Table 3). Some of the identified DASGs have been reported to be associated with CC pathophysiology at the gene expression levels41 but have not been explored at isoform levels. The ensuing discussion explains the potential role of representative DASGs identified from correlation analysis as well as from other DASGs identified in the study.

Extracellular matrix protein

Collagen and its family members are one of the most abundant ECM proteins and are associated with a range of muscle diseases.42 Collagen gene expression levels have been shown to be down‐regulated in muscle tissue in various catabolic states, including CC.41 DASGs of collagen such as COL18A1 and COL20A1 were down‐regulated in our current study. Matrix metalloproteinases, an ECM remodelling enzymes, play an important role in the breakdown of ECM components in normal physiological process and also in tissue turnover.43 ADAM 10, a cell surface protein identified in this study, has been shown to be a critical player in maintenance of satellite cells.44 Up‐regulation of this isoform in CC context needs to be elucidated in future experiments using model systems.

Inflammation

Systemic inflammation is a hallmark of CC in which there is an imbalance between the levels of pro‐inflammatory and anti‐inflammatory cytokines.45 While interleukin‐1 (IL‐1) and IL‐6 have been implicated in pathogenesis of CC,46 enhanced IL‐18 levels have been associated with fat loss and cachexia.47 Evidence suggests that mRNA levels of IL‐18 are increased in the skeletal muscle of COPD patients when compared with controls and may potentially play a role in muscle wasting.48 Based on this premise, it could be conjectured that IL‐18 isoforms may also play a role in cancer‐induced muscle wasting but this needs to be validated in future studies. ADIPOR1 is an adipokine molecule that is abundantly expressed in skeletal muscle. Its spliced isoforms have been shown to play a role during myogenesis and also in insulin sensitivity.49 Up‐regulation of ADIPOR1 in the skeletal muscle of CC may potentially lead to impaired myogenesis. Both ADIPOR1 and IL‐18 are involved in cytokine signalling, deregulation of which may contribute to CC pathophysiology. Other DASGs that are associated with cytokine signalling, B cell activation are listed in Table 3.

Lipid biosynthesis and fatty infiltration in skeletal muscle

Abnormal deposition of fat in organs such as liver, muscle, and bone are being recognized as pathogenic. While abnormal accumulation of fat in liver (hepatosteatosis) is well documented,50 molecular mechanisms involved in fatty infiltration of skeletal muscle (myosteatosis) have not been addressed. Myosteatosis has been found to be associated with insulin resistance and is also known to reduce the survival period in cancer patients.27 In the present study, cachectic cases have lower muscle radiation attenuation relative to WS cancer patients, which likely indicates the presence of fat accumulation in muscle. FABP4 has been reported to be expressed in skeletal muscle and play a role in fatty acid transport.51 FAR1 is studied for its role in β oxidation of fatty acids and acetyl‐CoA translocation.52 Up‐regulation of both these DASGs may potentially lead to accumulation of fat in skeletal muscle. Evidence also suggests that sphingolipid accumulation contributes to increased fat accumulation in skeletal muscle.53 The current study has identified DASGs such as B4GALT4 and ST8SIA5, which are associated with sphingolipid biosynthesis. With the potential role of the above‐mentioned DASGs leading towards lipogenesis and adipogenesis, it may be inferred that they may play a role in conferring the fatty muscle characteristics.

Protein ubiquitination

In CC, increased protein degradation or decreased protein synthesis or sometimes, both are observed. One of the well‐studied pathways for muscle atrophy is the ubiquitin proteasome pathway (UPP).13 Some of the identified DASGs such as UBA3, UBE2B, UBE2G2, PSMA5, and UBR5 have been implicated to play a role in protein ubiquitination. PSMA5 was shown to be up‐regulated in various catabolic states in animal models,41 and our study also found its isoform to be up‐regulated. UBE2B was also shown to be up‐regulated in vitro leading to myofibrillar protein loss.54 The current study has also identified UBE2B to be up‐regulated. Results from in vivo and in vitro models suggest that UBR5 acts as an activator of smooth muscle differentiation by stabilizing myocardin protein.55 While not much has been reported on spliced isoforms of UBR5, sequencing studies in mantle cell lymphoma have identified a high number of lethal mutations in mantle cell lymphoma patients, which include splice site mutations.56 However, defective splicing of UBR5 on muscle and muscle‐related conditions has not been investigated yet. Other up‐regulated DASGs associated with ubiquitination may also play a role in protein degradation pathways but have not been reported in the CC literature.

Skeletal muscle function and differentiation

Approximately, 5% of the DASGs identified were associated with skeletal muscle function in some capacity, ranging from muscle structure to extracellular matrix protein (Table 3). ROCK2 identified in this study plays a role in myoblast fusion during myogenesis by interacting with an RNA‐binding protein.57 ACVR2A, up‐regulated in the current study, is involved in BMP signalling and plays a role in maintaining muscle mass. Their functional role in CC pathophysiology needs to be elucidated in future studies.58 FLRT3 is a cell surface protein that is expressed during the somite development and plays a role in cell adhesion and in FGF signalling.59 However, their role in adult skeletal muscle and in disease states remains to be established. Our IPA analysis also predicted CXCR4 signalling pathways to be associated with CC. Martinelli et al.,60 recently showed that down‐regulation of genes involved in this pathway led to muscle wasting in Yoshida hepatoma‐bearing rodents and subsequent activation‐reduced muscle wasting. The study also showed that the levels of SDF1 isoforms were reduced in muscle of cachectic mice.60 Although SDF1 did not meet the predefined cut‐off for DASG, a similar trend was also observed in our study (down‐regulated in cachectic cases when compared with WS cancer patients, data not shown). It would be interesting to further investigate if this pathway can be targeted for potential interventions for human CC.

This study has a few limitations. Firstly, replication of findings using independent samples is required, and secondly, the identified signatures need to be interrogated for their putative biological roles in the context of CC in appropriate model systems. Our association study premise was based on cachectic cases and non‐cachectic controls. Weight stability for 6 months, high BMI, and the absence of sarcopenia in the non‐cachectic cases would appear to be a reasonably robust set of criteria for establishing the absence of cachexia. However, we acknowledge the limitation that we cannot prove that these individuals did not have precachexia. Indeed, some of the non‐cachectic controls may eventually go on to gain or lose skeletal muscle or adipose tissue in their cancer trajectory. A comparison of non‐weight losing cancer patients with matched healthy control subjects would be valuable. Patient consent and ethical approval for muscle biopsies from healthy controls are an issue; healthy controls have been infrequently included in prior studies and where present very small numbers (i.e. n = 6)38 have been obtained.

Potential influence of polymorphisms (single nucleotide polymorphisms or SNPs) on promoter, enhancer, exonic regions, or splice acceptor/donor sites could affect expression of transcript levels or encoded protein functions.61 Focused studies are needed to characterize the influence of SNPs on splicing mechanisms and their contribution to CC pathophysiology.62 Therapeutic potential of exon skipping mechanism has been demonstrated for Duchenne muscular dystrophy.63 Such studies are encouraging as some of the DASGs identified in this study exhibiting exon skipping mechanism may also be explored for potential interventions for CC in future.

To put the current study into perspective, this is one of the largest sample sizes used to date for profiling AS in human skeletal muscle and CC literature. We recognize that independent validation of findings is needed and international collaborations are sought to gain access to the rare and precious source of skeletal muscle biopsies from cancer‐affected patients. Being the first study in CC and AS, we chose to profile all of the samples available (n = 40) and not split the samples into subsets for use in discovery and validation stages, as such a stratification attempt may weaken the statistical power. Although we recognize that gender and age may have an impact on expression studies,64 we have refrained from such stratified analyses due to sample size limitations. It would be interesting to conduct these studies in future using well‐powered cohorts and validate in independent cohorts. If muscle biopsies from patients affected by different tumour types are accrued, this can also help us identify tumour‐specific cachexia signatures, which may aid in therapeutic interventions in future, to support the premise of personalized medicine.

Conflict of interest

The authors declare that they have no conflict of interest.

Supplementary information Table S1: List of up‐regulated DASGs identified at 1.4 FC, p < 0.05

Click here for additional data file.

Supplementary information: Table S2: List of down‐regulated DASGs identified at 1.4 FC, p < 0.05

Click here for additional data file.

Supplementary information: Table S3: The thermal cycle profiles and the primers sequences used for validating representative DASGs using semi‐quantitative RT‐PCR

Click here for additional data file.

Supplementary information: Table S4: List of significant (p < 0.05) canonical pathways identified using IPA for DASGs.

Click here for additional data file.