Small RNAome profiling from human skeletal muscle: novel miRNAs and their targets associated with cancer cachexia.

BACKGROUND: MicroRNAs (miRs) are small non-coding RNAs that regulate gene (mRNA) expression. Although the pathological role of miRs have been studied in muscle wasting conditions such as myotonic and muscular dystrophy, their roles in cancer cachexia (CC) are still emerging.
OBJECTIVES: The objectives are (i) to profile human skeletal muscle expressed miRs; (ii) to identify differentially expressed (DE) miRs between cachectic and non-cachectic cancer patients; (iii) to identify mRNA targets for the DE miRs to gain mechanistic insights; and (iv) to investigate if miRs show potential prognostic and predictive value.
METHODS: Study subjects were classified based on the international consensus diagnostic criteria for CC. Forty-two cancer patients were included, of which 22 were cachectic cases and 20 were non-cachectic cancer controls. Total RNA isolated from muscle biopsies were subjected to next-generation sequencing.
RESULTS: A total of 777 miRs were profiled, and 82 miRs with read counts of ≥5 in 80% of samples were retained for analysis. We identified eight DE miRs (up-regulated, fold change of ≥1.4 at P < 0.05). A total of 191 potential mRNA targets were identified for the DE miRs using previously described human skeletal muscle mRNA expression data (n = 90), and a majority of them were also confirmed in an independent mRNA transcriptome dataset. Ingenuity pathway analysis identified pathways related to myogenesis and inflammation. qRT-PCR analysis of representative miRs showed similar direction of effect (P < 0.05), as observed in next-generation sequencing. The identified miRs also showed prognostic and predictive value.
CONCLUSIONS: In all, we identified eight novel miRs associated with CC.

Introduction

Cancer cachexia (CC) is a multifactorial syndrome characterized by severe depletion of skeletal muscle with or without fat loss. The pathophysiology includes systemic inflammation, reduced food intake and negative energy and protein balance.1 Patients with CC have poor treatment response, reduced survival and a severe compromise in their quality of life.2 CC is characterized by complex host–tumour interactions that remain to be fully elucidated. Tumour derived mediators lead to aberrant host tissue responses. Several proinflammatory and procachectic cytokines released from tumour cells contribute to systemic inflammation in the host and lead to metabolic alterations.3 Genes (mRNA) involved in the pathophysiology of CC have been fairly well studied.4, 5 However, the role of finer post‐transcriptional gene regulatory mechanisms and its implications on CC at the whole genome level has not been comprehensively explored.

MicroRNAs (miRs, 18‐25 nt) are a class of small non‐coding RNAs that are considered as global regulators of gene expression. They primarily bind to 3' untranslated region (3'UTR) of mRNA and cause either translational repression or mRNA degradation, depending on the degree of complementarity shared between the two molecules.6 Myo‐miRs, a suite of miRs highly enriched in skeletal muscle, are known to play a role in myogenesis.7 For example, miR‐1 and miR‐133 regulate skeletal muscle proliferation and differentiation.8 Similarly, diverse roles of miRs have been reported in physiological process such as adipogenesis, exercise and in general, health and disease states9, 10, 11 and have also been identified as promising biomarkers for many diseases.12, 13 Recent evidence has suggested that miR‐21 and miR‐378 promote muscle cell apoptosis and adipolysis, respectively, in CC.14, 15 While these studies have focused on animal models and human adipose tissues, it remains to be explored if miRs identified from human muscle associate with CC pathophysiology and can also act as potential biomarkers. Biomarkers for CC have previously been identified from preclinical models.16, 17 Nonetheless, many of these molecules are far from being universally accepted as these are yet to be validated in human subjects.

Whole genome miR profiling using next‐generation sequencing (NGS) are increasingly being utilized, as NGS offers better sensitivity and specificity compared with hybridization techniques.18 However, comprehensive profiling of miRs for CC using NGS platform has not yet been attempted. We hypothesized that deregulation of miRs contributes to the aetiology of CC. The study objectives were (i) to profile human skeletal muscle expressed miRs; (ii) to identify differentially expressed (DE) miRs between cachectic and non‐cachectic cancer patients; (iii) to identify mRNA targets for the DE miRs to gain mechanistic insights; and (iv) to explore the prognostic and predictive potential of miRs. We report eight novel miRs associated with CC pathophysiology. We have also identified gene targets (mRNA) for the miRs and potential regulation of canonical pathways by these miR‐mRNA pairs.

Methods

Recruitment of study subjects and acquisition of muscle samples

Skeletal muscle biopsies from patients who underwent laparotomy at the Foothills Medical Centre between 2006 and 2013 were obtained from the University of Calgary Hepatopancreaticobiliary/Gastrointestinal Tumour Bank. Patients were approached for consents in pre‐admission clinics (out‐patient) prior to surgery. All patients provided written informed consent for study participation, and the study was approved by the Conjoint Health Research Ethics Board at the University of Calgary (Ethics ID E‐17213). Biopsies of rectus abdominus were taken at the start of the surgery using sharp dissection, immediately frozen in liquid nitrogen to minimize ischemic shock post‐devitalization, and stored at −80°C until further use. Post‐biopsy molecular profiling of samples for NGS and collection of patient information was conducted under ethics protocol ETH‐21709, which was approved by Health Research Ethics Board of Alberta (HREBA)–Cancer Committee.

Study design and clinical annotation of cases

Weight loss (WL) information of patients (n = 42) was retrieved from medical charts. Computed tomography (CT) scans were available for majority of the samples (n = 35), and these were accessed for quantification of skeletal muscle and fat components. Patients were classified as cachectic based on the international consensus framework for CC; patients belonging to any one of the three diagnostic criteria were considered cachectic (more than 5% pre‐illness WL for a period of six months, body mass index (BMI) < 20 with any degree of WL >2% and lumbar skeletal sarcopenia, as defined by muscle index cut‐points using computed tomography).1 Non‐cachectic cancer patients had no WL over a period of 6 months compared with their pre‐illness weight. The exclusion criteria include as follows: (i) patients with no clinical chart or recorded WL information; (ii) below 18 years; and (iii) inability to give written informed consent. A subset of age matched patients in the bank meeting the aforementioned inclusion criteria and with available biopsies, and WL histories were accessed; of these 22 were cachectic cases (with mean WL of 11.8 ± 6.6%; hereafter referred to as cases) and 20 were non‐cachectic controls (hereafter referred to as controls). One patient had completed a course of chemotherapy (neo‐adjuvant) before surgery but had not received any chemotherapy 4 weeks prior to surgery. Remaining patients did not receive any chemotherapy before surgery.

Pancreatic and colorectal cancer (with liver metastasis) patients were considered for the present study from the collection of biopsies from the Hepatopancreaticobiliary/Gastrointestinal Tumour Bank to minimize the number of tumour types. Tumour stage is reported according to American Joint Committee on Cancer v7.

Computed tomography image quantification

Twenty cases and 15 controls had a single CT prior to surgery (70.75 ± 45.24 in days). The quantification protocol has been explained elsewhere in detail.19, 20 Briefly, the third lumbar vertebrae were used as a standard landmark to measure the cross‐sectional muscle area (cm2) and normalized to their stature (m2) to calculate the skeletal muscle index (SMI) (cm2/m2). Muscle attenuation (MA) (Hounsfield Units) was also captured for these patients. Sarcopenia status was assigned based on age and sex adjusted SMI values, as described earlier.21

RNA extraction and next‐generation sequencing profiling

Total RNA from muscle biopsies were isolated using Trizol (Sigma‐Aldrich, Oakville, ON, Canada) and purified using Qiagen's RNeasy Maxi kit (Mississauga, ON, Canada). Quantification of RNA was carried out using Nanodrop 1000 Spectrophotometer. All of the aforementioned protocols were carried out according to the manufacturer's instructions.

Services from PlantBiosis Ltd (Lethbridge, Alberta, Canada; http://www.plantbiosis.com/) were used from small RNA library preparation up to the generation of alignment files (.bam files), as described earlier.22 All the samples were sequenced in a single batch in a single lane to avoid batch effects. Sequencing of small RNA was carried out using Truseq Small RNA sequencing kit (Illumina), TruSeq SR Cluster Kit v5‐CS‐GA (Illumina) and TruSeq SBS Kit v5‐GA (Illumina), according to manufacturer's instructions. The samples were sequenced using MiSeq platform with 36 bases single‐end protocol. Base calling and demultiplexing was performed using MiSeq Reporter FastQ workflow. Cutadapt 1.4.1 (https://pypi.python.org/pypi/cutadapt/1.4.1) was used to trim the adapters. Only reads with a quality score of more than 30 on the Sanger scale were considered. Quality control for sequenced reads before and after adapter trimming was performed using FastQC v0.11.3 software. Burrows–Wheeler aligner version 0.6.1 was used to align trimmed reads to the reference genome23 (human genome 19 build), downloaded from Illumina iGenome website. Samtools version 0.1.18 was used to convert sequence aligned data file (.sam files) to its binary format (.bam format), which was used for downstream analyses. The raw files and normalized counts of the NGS data have been submitted to Gene Expression Omnibus (accession ID GSE75473).

Identification of differentially expressed microRNAs

Partek Genomics Suite 6.6 (PGS 6.6) was used to analyze the .bam files generated from NGS experiment. Raw data were filtered to include miRs with more than or equal to five read counts in at least 80% of the samples and were normalized using reads per kilobase per million (RPKM) method.24 DE miRs were identified with a fold change (FC) of >1.4 and P < 0.05 using one‐way analysis of variance. The DE miRs were also subjected to permutation test (n = 10 000), and miRs with permuted P‐values < 0.05 were considered for subsequent analysis.

qRT‐PCR validation of the identified microRNAs

From the eight DE miRs identified from NGS, representative miRs: hsa‐miR‐3184‐3p, hsa‐let‐7d‐3p and hsa‐miR‐1296‐5p were selected and validated using qRT‐PCR. Total RNA was reverse transcribed using TaqMan MicroRNA reverse transcription kit. RNU6B was used as the internal control. qPCR was performed using 7900HT fast real‐time PCR system (Applied Biosystems). The 20 μL reactions for each sample were run in triplicates according to manufacturer's protocol. The following thermal cycler conditions were used for reverse transcription: 30 min. at 16°C, 30 min. at 42°C and 5 min. at 95°C; and the following conditions were used for qPCR: 50°C for 2 min, 95°C for 20 s and 40 cycles of 95°C and 60°C for 1 s and 20 s, respectively. Data analysis was carried out using 2‐∆∆Ct method.25

Target predictions and putative functional annotation for mRNA targets of differentially expressed microRNAs

Detailed description of the in‐house muscle transcriptome gene (mRNA) expression data generated using Agilent platform is elaborated elsewhere26 and has been deposited in Gene Expression Omnibus repository (GSE41726). Initially, independent muscle biopsies (non‐matched data sets) were used for miR profiling (this study, n = 42) and mRNA profiling (n = 90, a previous study GSE41726 from the same lab). Briefly, using the same definitions as in the miR study, 29 cachectic cases (WL > 5%) and 61 weight stable (WS) controls were used to identify DE mRNAs. PGS 6.6 was used for differential expression analysis. Raw intensity files were quantile normalized, log 2 transformed and DE mRNAs were identified at ≥ 1.4 FC and P < 0.05 using one‐way analysis of variance. The direct binding of miR to 3'UTR of the mRNA is one of the most commonly recognized mechanisms in post‐transcriptional silencing of genes.27, 28 However, miRs may exert influence on target mRNAs in an indirect manner29, and hence, we considered all miR‐mRNA correlations. We interrogated the mRNA data showing inverse and similar directional effect of expression between miR and mRNA (i.e. miR up‐regulated and its corresponding mRNA down‐regulated and both miR and mRNA being up‐regulated/down‐regulated) to understand the biological role of identified miRs in CC.22 Putative gene targets for DE miRs were predicted in silico using TargetScan 7.0,30 based on the complementary binding of the seed region of miRs with 3'UTR of mRNAs. These predicted targets were then overlapped with the DE genes identified from in‐house muscle transcriptome dataset. Ingenuity pathway analysis (IPA) was used to identify pathways for the overlapped gene targets of eight DE miRs. We have further interrogated the mRNA targets (GSE85017) for the miRs in matched data sets, albeit of lower sample size (n = 42), compared with the unmatched mRNA data set (n = 90). Independent confirmation of findings from two sources, matched and unmatched data sets as an approach for functional validation of targets in a tissue‐specific context would strengthen the study findings.

Statistical analysis

Data are presented as mean ± standard deviation. Independent t‐test and χ² test were used for continuous and categorical variables, as appropriate. Sample size estimation was carried out using the following parameters: α = 0.05, β = 80% and a fold difference of 1.4 or more in miR expression. Nineteen samples per group were required to conduct the study (http://bioinformatics.mdanderson.org/MicroarraySampleSize/) to meet the statistical power and identify DE miRs with confidence.

Identification of microRNAs as prognostic and predictive factors

One of the objectives of this study was to identify miRs associated with overall survival (OS). OS was defined as the time period between the date of surgery and the date of death. The median follow‐up period was 1060 days (Range: 6 to 2041 days), from the date of muscle biopsy accrual till the last follow‐up date (March 2014). Of the 42 patients, 19 died, and 23 were alive at the time of completion of this study.

Prognostic value: we considered DE miRs as continuous variables and subjected the miRs to univariate cox proportional hazards regression model. For OS, composite risk score was constructed for each sample using the parameter estimates (co‐efficient) and the normalized counts obtained for the eight DE miRs. Receiver operating characteristics curve was used to determine the optimal cut‐point for the composite risk score. This was carried out to dichotomize the 42 samples into high‐risk and low‐risk groups. Hazard ratios (HR) along with 95% confidence intervals were reported for univariate and multivariate analyses. For multivariate analysis, the risk score was adjusted for age, BMI and tumour type, where appropriate. Log‐rank test was carried out to identify if any significant difference existed between the survival curves of two‐risk groups.

Predictive value: the association between the weight change (loss vs. stable) group and the miRNA score (composite risk score) was tested using binary logistic regression. Univariate and multivariate logistic regression models are reported. For the multivariate model, we adjusted for age and BMI as continuous variables and tumour type (colorectal vs. pancreas) as dichotomous variable. Odds ratio (OR) and the corresponding 95% confidence interval are reported. P‐value < 0.05 was used to define statistical significance. Two‐sided tests were used for the comparisons. All statistical analyses were conducted using spss v16.0 and sas (SAS Institute Inc., Cary, NC) v 9.3.

Results

Patient demographics and body composition analysis

Forty‐two patients who met the study criteria with information on age, gender, BMI, tumour type and physician‐documented WL information were selected. There was no significant difference between cases and controls in age, gender and tumour type. BMI was found to be significant between the groups (P = 0.01) (Table 1a). We investigated if the CT‐derived body composition measurements also are in concordance with the phenotypes of cachexia in the study cohort (Table 1b). We observed expected trends in SMI values between cases and controls, (P = 0.07 for SMI). MA was found to be reduced in cases compared with controls and were significant at P = 0.04.

Table 1a

Patient demographics

Characteristics	Cachectic cases (n = 22)	Non‐cachectic controls (n = 20)	P‐values
Weight loss, % mean	11.8 ± 6.6	—
Age (mean, in years) [Range]	64.9 ± 10.1	63.6 ± 7.9	0.6a
	[40–83]	[45–76]
Tumour type (n)			0.2b
Pancreatic	12	7
Colorectal	10	13
Tumour stage (n)			0.6c
I	2	1
II	3	3
III	2	0
IV	15	16
Gender (n)			0.7b
Male	9	9
Female	13	11
Body mass index (mean, in kg/m2)	24.35 ± 2.5	27.02 ± 3.7	0.01a
[Range]	[19–29]	[20–40]

Unpaired t‐test.

χ² test.

Fisher's exact test.

Data represented as mean ± standard deviation. Statistical analyses were carried out using spss v16. Independent t‐test was conducted for age and body mass index. χ² test was carried out for tumour type and gender. P < 0.05 was considered statistically significant.

Table 1b

Body composition analysis for the study subjects

Characteristics	Cachectic cases (n = 20)	Non‐cachectic controls (n = 15)	P‐values
Cross‐sectional skeletal muscle area (cm2)a
Male	139.5 ± 15.4	158 ± 12.9	0.2
Female	96.2 ± 14.6	103.5 ± 14.1
Skeletal muscle index (cm2/m2)a			0.07
Male	45.3 ± 5.9	49.1 ± 3.1
Female	36.3 ± 5.6	41.52 ± 7.43a
z‐score	−0.7 ± 0.7	−0.08 ± 0.86	0.03
Total adipose tissuea			0.9
Male	226.8 ± 84.5	266.29 ± 77.2
Female	328.8 ± 126.4	302.2 ± 122.7
Muscle attenuation (HU)a			0.04
Male	34.14 ± 8.7	39.8 ± 6.9
Female	29.46 ± 7.3	36.42 ± 8.4

Unpaired t‐test.

Body composition was calculated for subset of patients (n = 35) who had computed tomography prior to surgery. Statistical analysis was conducted between cases and controls. Muscle attenuation values were significant in the overall comparison between cachectic cases and non‐cachectic controls. z‐score is the difference expressed in standard deviation of patients' values from age and gender‐specific mean values.21

Technical variation, data processing and identification of differentially expressed microRNAs using next‐generation sequencing

Earlier reports suggested that the technical variance arising from the NGS platform is minimal.31 This has been validated in the current investigation as well. Technical replicates of four samples (two cases and two controls) were sequenced. Excellent concordance was observed (r > 0.99) between the replicates, suggesting that technical variance in the data is absent or minimal because of experimental procedures adopted ( ).

A total of 7 926 299 and 9 155 808 reads were obtained from cases and controls (n = 42), respectively. A total of 97.12% (7 698 807) and 97.35% (8 913 914) of reads were retained after adapter trimming in cases and controls, respectively. A total of 87.78% (6 758 570) of reads in cases, and 87.64% (7 812 866) of reads in controls were aligned to the reference genome, of which 9.55% (646 069) of reads in cases and 8.14% (636 748) of reads in controls mapped to miRs (Table 2) and had read length distribution from 18 to 25 nucleotides (Figure 1).

Table 2

Descriptive statistics for the data obtained from next‐generation sequencing

Samples	Total reads	Reads retained after adapter trimming	Aligned reads	Reads mapped to miRNAs
Cachectic cases (n = 22)	7 926 299	7 698 807 (97.12%)	6 758 570 (87.78%)	646 069 (9.55%)
Non‐cachectic controls (n = 20)	9 155 808	8 913 914 (97.32%)	7 812 866 (87.64%)	636 748 (8.14%)

Step‐wise filtering of next‐generation sequencing data is shown starting from total reads obtained from next‐generation sequencing for both cases and controls to reads mapped to miRNA. Furthermore, the reads mapped to miRNA were subjected to differential expression analysis.

Figure 1

Length distribution of reads aligned to miRs: the aligned read length ranges from 17 to 27 nucleotides with the maximum distribution of reads captured between 18 and 25 nucleotides (reflecting miRNA read length).

A total of 777 miRs were profiled from the skeletal muscle tissues (out of a total of 2588 miRs annotated in miRbase v20, of which only a subset can be tissue specific). Eighty‐two miRs were retained after filtering and were subjected to DE analysis. A total of eight miRs were found to be up‐regulated with a FC of >1.4, P < 0.05 and with false discovery rate (FDR) ranging from 0.21 to 0.22. Permutation tests (n = 10 000) for the eight DE miRs were carried out and were also found to be significant with permuted P < 0.05. The eight up‐regulated DE miRs were hsa‐miR‐3184‐3p, hsa‐miR‐423‐5p, hsa‐let‐7d‐3p, hsa‐miR‐1296‐5p, hsa‐miR‐345‐5p, hsa‐miR‐532‐5p, hsa‐miR‐423‐3p and hsa‐miR‐199a‐3p (Table 3a). No down‐regulated miRs were identified in the current study under the stringent data filters applied and with the current sample size.

Table 3a

Differentially expressed miRNAs

miRNA	P‐value	FDR	Permutation P‐value	Fold change (Cachectic cases vs. non‐cachectic controls)	Direction of fold change
hsa‐let‐7d‐3p	0.01	0.21	0.01	1.48	Up
hsa‐miR‐345‐5p	0.02	0.21	0.02	1.47	Up
hsa‐miR‐423‐5p	0.009	0.21	0.009	1.42	Up
hsa‐miR‐532‐5p	0.02	0.21	0.02	1.48	Up
hsa‐miR‐1296‐5p	0.03	0.22	0.03	1.44	Up
hsa‐miR‐3184‐3p	0.009	0.21	0.008	1.42	Up
hsa‐miR‐423‐3p	0.01	0.21	0.01	1.43	Up
hsa‐miR‐199a‐3p	0.03	0.22	0.01	2.01	Up

FDR, false discovery rate.

All eight differentially expressed miRNAs were up‐regulated in cachectic cases with a fold change of ≥1.4 at P < 0.05. Permutation test was carried out (n = 10 000) for these seven miRNAs to rule out observation by chance. The permuted P‐value was significant for all the eight differentially expressed miRNAs. FDR for all eight differentially expressed miRNAs were also represented.

Validation of select microRNAs using qRT‐PCR

Representative miRs such as hsa‐miR‐3184‐3p, hsa‐let‐7d‐3p and hsa‐miR‐1296‐5p were validated in an independent platform – qRT‐PCR. All the three miRs showed similar direction of expression as observed in NGS with significant P‐values (P < 0.05) (Figure 2).

Figure 2

qRT‐PCR validation of miRs: representative miRNAs were validated using qRT‐PCR. Similar direction of effect was observed as seen in the next‐generation sequencing with statistical significance of P < 0.05.

Gene expression studies for identification of putative targets for differentially expressed microRNAs

TargetScan identified 20 988 mRNAs as potential targets for eight DE miRs. TargetScan in silico predicted targets included mRNAs from muscle and non‐muscle cell and tissue types. Therefore, to identify tissue‐specific gene (mRNA) targets for the DE miRs, an in‐house muscle transcriptome dataset was used.26 A total of 353 mRNAs with a FC of >1.4 and P < 0.05 were identified between cachectic cases (>5% WL) and non‐cachectic controls (WS), of which 127 were up‐regulated and 226 were down‐regulated (data not shown). Because we do not know a priori which among the 353 mRNAs serve as potential targets for the DE miRs in our study, we mapped the 353 mRNAs to the 25 348 mRNA targets predicted by TargetScan. In this analysis, we identified 191 mRNA targets as potentially regulated by the eight DE miRs (Table 3b). Up to 70% of the targets identified were down‐regulated, and is an expected consequence of the binding of miRs to 3'UTRs of mRNAs leading to gene silencing. To further confirm the direction of expression of these targets, an independent gene expression study for matched samples was conducted using Affymetrix Human Transcriptome Array 2.0 (GSE85017). We compared the direction of expression of the identified 191 mRNA targets in this independent dataset. We found 77% (147/191) of the mRNA targets expressed in the same direction, in these comparisons regardless of the P‐value and FC (data not shown). Further, the direction of expression of 10 representative mRNA targets chosen for discussion (see below) is indicated in . Although the P‐value and FC of these 10 mRNA targets were small in the independent dataset, perhaps due to the small sample size, we have confirmed the direction of expression in matching samples.

Table 3b

Summary of target identification

miRNA	Number of targets identified by target scan	Number of DE gene targets identified by in‐house datasets
hsa‐let‐7d‐3p	479	4
hsa‐miR‐345‐5p	3515	34
hsa‐miR‐423‐5p	4437	41
hsa‐miR‐532‐5p	3227	32
hsa‐miR‐1296‐5p	2007	13
hsa‐miR‐3184‐3p	3531	28
hsa‐miR‐423‐3p	779	14
hsa‐miR‐199a‐3p	3013	25

Differentially expressed (DE) miRs were mapped to TargetScan 7.0. The identified gene targets from TargetScan 7.0 were overlapped to in‐house muscle transcriptome dataset to identify tissue‐specific gene targets for DE miRNA. The identified targets were interrogated for pathway analysis using ingenuity pathway analysis to understand the potential biological roles of the miRNAs in cancer cachexia.

Pathway analysis using ingenuity pathway analysis

The 191 mRNA targets identified for eight DE miRs were subjected to IPA to identify the corresponding pathways and understand their biological relevance (Table 4). The identified targets are known to play a role in adipogenesis, myogenesis (SULF1, BMPR1B and DLK1), inflammation and innate immune response (RPS6KA6) and also in signal transduction pathways (SFRP4 and DKK2) that may directly or indirectly contribute to the phenotype of CC (See Discussion for details).

Table 4

List of significant pathways identified from ingenuity pathway analysis

Ingenuity canonical pathways	Molecules	miRNA‐ID
Actin cytoskeleton signalling	FN1	hsa‐miR‐199a‐3p
Adipogenesis pathway	DLK1	hsa‐miR‐345‐5p, hsa‐miR‐423‐5p, hsa‐miR‐3184‐3p
BMP signalling pathway	GREM1	hsa‐miR‐345‐5p, hsa‐miR‐3184‐3p, hsa‐miR‐199a‐3p
	BMPR1B	hsa‐miR‐3184‐3p
	SULF 1	hsa‐miR‐ 532‐5p
Calcium signalling	CAMK2A	hsa‐hsa‐miR‐423‐3p
Cholesterol biosynthesis I	SQLE	hsa‐miR‐3184‐3p
CNTF signalling	RPS6KA6	let‐7d‐3p, hsa‐miR‐1296, hsa‐miR‐199a‐3p, hsa‐miR‐532‐5p
Energy metabolism	NYP1R	hsa‐miR‐532‐5p
GDNF family ligand–receptor interactions	RET	hsa‐miR‐423‐5p, hsa‐miR‐3184‐3p
Glucocorticoid receptor signalling	PGR	let‐7d‐3p, hsa‐miR‐1296, hsa‐miR‐423‐5p
Glutamate receptor signalling	SLC1A7	hsa‐miR‐345‐5p
IGF‐1 signalling	CYR61, NOV	hsa‐miR‐345‐5p
IL‐6 signalling	COL1A1	hsa‐miR‐345‐5p, hsa‐miR‐423‐5p
IL‐8 signalling	EIF4EBP1	hsa‐miR‐423‐5p, hsa‐miR‐199a‐3p
Insulin receptor signalling	EIF4EBP1	hsa‐miR‐423‐5p
Integrin signalling	CAPN6	let‐7d‐3p, hsa‐miR‐1296, hsa‐miR‐345‐5p, hsa‐miR‐423‐5p, hsa‐miR‐3184‐3p
Mitochondrial dysfunction	SOD2	hsa‐miR‐345‐5p, hsa‐miR‐423‐5p, hsa‐miR‐3184‐3p, hsa‐miR‐199a‐3p
mTOR signalling	RPS6KA6	let‐7d‐3p, hsa‐miR‐1296, hsa‐miR‐199a‐3p
	EIF4EBP1	hsa‐miR‐423‐5p, hsa‐miR‐199a‐3p
NF‐κB signalling	BMPR1B	hsa‐miR‐3184‐3p
Oleate biosynthesis II	FADS2	hsa‐miR‐423‐5p
Phospholipase C signalling	BLNK	hsa‐miR‐345‐5p
Regulation of cellular mechanics by calpain protease	CAPN6	let‐7d‐3p, hsa‐miR‐1296, hsa‐miR‐345‐5p, hsa‐miR‐423‐5p, hsa‐miR‐3184‐3p
Serotonin receptor signalling	HTR2A	let‐7d‐3p, hsa‐miR‐1296, hsa‐miR‐3184‐3p, hsa‐miR‐199a‐3p
TGF‐β signalling	BMPR1B	hsa‐miR‐3184‐3p
Wnt/β‐catenin signalling	SFRP4	let‐7d‐3p, hsa‐miR‐1296, hsa‐miR‐3184‐3p

CNTF, ciliary neurotrophic factor; IGF‐1, insulin‐like growth factor gene 1; mTOR, mammalian target of rapamycin; TGF‐β, transforming growth factor beta.

All miRs selected for this analysis showed up‐regulation. Majority of the genes indicated in the table showed down‐regulation, as expected in the mRNA–miRNA correlations. Note that common molecules between pathways reflect in the redundancy of the genes due to multiple miRs regulating the same pathway.

microRNAs as potential independent prognostic factors

The eight DE miRs were considered as continuous variables and were subjected to univariate cox proportional hazards regression model. The parameter estimates thus obtained were used for constructing a risk score. A risk score cut‐point was estimated to dichotomize the patients into high‐risk and low‐risk groups. Cut‐point above 7.6 was considered as high risk and less than equal to 7.6 as low risk, following which the risk score was treated as categorical variable for univariate cox model. The risk score showed trend towards significance after adjusting for potential confounding factors (age, BMI and tumour type) in the multivariate analysis. The high‐risk group had a shorter OS when compared with the low‐risk group [HR: 2.32 (0.88–6.06), P = 0.08, Table 5a]. The log‐rank P‐value was significant in survival analysis between two‐risk groups (P = 0.001, Figure 3).

Table 5a

Univariate and multivariate results for overall survival

Parameter	Univariate analysis		Multivariate analysis
	HR (95% CI)	P‐value	HR (95% CI)	P‐value
Risk score	3.49 (1.51–8.04)	0.003	2.32 (0.88–6.06)	0.08
Age	0.99 (0.94–1.04)	0.81	0.95 (0.89–1.00)	0.07
BMI	0.97 (0.86–1.10)	0.67
Tumour type	0.21 (0.09–0.50)	0.0005	0.22 (0.07–0.63)	0.005

BMI, body mass index; CI, confidence interval; HR, hazards ratio.

Risk score and other clinical parameters were subjected to univariate cox proportional hazards model. In the multivariate analysis, risk score was marginally significant after adjusting for all potential confounders (age, BMI and tumour type).

Figure 3

Kaplan–Meier plot for risk score: Kaplan–Meier plot was constructed to assess the survival function of high‐risk group vs. the low‐risk group, based on the risk score. The high‐risk group had a shorter OS when compared to the low‐risk group. The log rank P‐value was significant in survival analysis between the two‐risk groups.

microRNAs as potential independent predictive factors

The association between the dichotomous risk score obtained from OS analysis for eight miR signature and weight change (cases vs. controls) was analyzed using binary logistic regression. Results from this analysis indicated that the composite risk score is an independent predictor of weight change. The odds of WL were seven times higher in patients with a high‐risk score compared with the low‐risk score patients. When adjusted for age, BMI and tumour type, the odds were still higher, about 7.95 times higher (1.44–43.97 at 95% confidence interval, Table 5b). Area under the curve was found to be 0.84 after adjusting for potential confounding factors (age, BMI and tumour type). The wider confidence interval observed may be attributed to a relatively small sample size used in the study.

Table 5b

Univariate and multivariate analysis for logistic regression

Parameter	Univariate analysis		Multivariate analysis
	OR (95% CI)	P‐value	OR (95% CI)	P‐value
Risk score (high score vs. low score)	7.00 (1.73–28.3)	0.006	7.95 (1.44–44.0)	0.01
Age	0.99 (0.92–1.06)	0.8
BMI	0.72 (0.53–0.96)	0.03	0.66 (0.45–0.96)	0.03
Tumour type (colorectal vs. pancreas)	0.45 (0.13–1.56)	0.21

BMI, body mass index; CI, confidence interval; OR, odds ratio.

The odds of weight loss was 7.95 times higher in patients with high‐risk score compared to low‐risk score after adjusting for potential confounders.

Discussion

We identified eight novel DE miRs from human skeletal muscle and also identified the potential contributions of miR‐mRNA correlations to the pathophysiology of CC. All eight DE miRs were up‐regulated, and none of the miRs reported here have been directly studied in the context of CC. This is the first study where NGS technology was used to comprehensively profile miRs on a genome‐wide scale for CC.

Forty‐two cancer patients were chosen and stratified as cases and controls for the study. BMI was less in cachectic cases when compared to controls (P = 0.01, Table 1a). In the current study, although cachectic cases have BMI in the overweight range, they have lost more than 10% of their weight in previous 6 months. In this population (Canada) considered for the study, WL regardless of BMI,32 as well as the combination of sarcopenia and WL, was shown to be associated with poor prognosis. Body composition analysis indicated that MA was reduced in cases compared with controls (P = 0.04). Previous studies suggest that cancer patients with reduced attenuation are associated with poor survival. In terms of attenuation values, we have observed a similar trend as reported in literature with large sample sizes.33 The reduced attenuation values in cases indicate that these muscles could have been infiltrated by fat.

Using in‐house muscle transcriptome dataset, we interrogated tissue (muscle) specific mRNA signatures to identify targets. The transcriptome dataset serves as a proxy for functional validation of eight DE miRs to gain mechanistic insights. A list of mRNA targets and inferred pathways for DE miRs is given in Table 4. However, further cell‐based assays would help confirm the specific targets and their regulation by their corresponding miRs. A brief discussion on eight DE miRs, their interacting partners (mRNAs) and their biological relevance to CC is described below. There is paucity of literature for the identified miRs associated with human CC, and hence, these comparisons are of high relevance.

We have identified Transferrin receptor (TFRC) as a putative target for let‐7d‐3p, in addition to three other targets. While let‐7d‐3p was up‐regulated in CC, its corresponding target, TFRC, was observed to be relatively down‐regulated in weight losing patients. Let‐7 members play a role in ageing; higher levels of let‐7 were observed in elderly subjects. They affect muscle cell proliferation by down‐regulating the target genes, eventually leading to reduced regenerative capacity of muscle.34 Transferrin, the ligand for TFRC, has been associated with myogenic differentiation, and down‐regulation of TFRC may impair myogenic differentiation. Furthermore, the expression of TFRC was observed to be down‐regulated in skeletal muscle and heart muscle in hereditary hemochromatosis and during muscle atrophy in conditions such as diabetes and uremia.35, 36 However, the role of TFRC in the context of CC remains to be elucidated.

Gene targets identified for miR‐345‐5p include NOV, COL1A1 and CYR61. While miR‐345‐5p was up‐regulated in CC, its targets NOV and COL1A1 were down‐regulated in the muscle transcriptome dataset, whereas CYR61 was found to be up‐regulated. NOV and CYR61 are known to be involved in insulin‐like growth factor 1 gene signalling via the AKT and mammalian target of rapamycin pathway,37 impairment of which is related to suppression of protein synthesis leading to muscle atrophy. COL1A1 is an extracellular matrix protein, and its down‐regulation in animal models is reported to be associated with skeletal muscle atrophy. let‐7d‐3p and miR‐345‐5p regulate genes such as TFRC and COL1A1 that are known to be involved in skeletal muscle atrophy in diverse disease states including CC.36

SQLE and FADS2 were the down‐regulated gene targets in muscle transcriptome dataset for miR‐423‐5p and miR‐3184‐3p; these miRs were up‐regulated and occur as clusters that are located within 10 kb distance.38 As they arise from the same clusters, their P‐value, FCs and direction of expression in the current study were found to be similar. Both these genes were found to be associated with lipid biosynthesis (Table 4) that is severely compromised in cachectic states. These genes contribute to sterol biosynthesis and fatty acid dehydrogenation. Additionally, miR‐423‐5p also regulates genes associated with energy homeostasis such as leptin. Perturbed levels of leptin in orexigenic and anorexigenic signals circuit and its impact on body weight regulation are fairly well understood. While increased leptin levels are observed in obesity, opposite effect is seen in CC.39 Leptin expression is low (but not lost) in mature skeletal muscle, but its expression is high in myocytes.40 Indeed, a similar trend was observed in our muscle transcriptome dataset, wherein leptin was down‐regulated in cases, relative to controls (data not shown). A similar argument may also be extended to miR‐423‐5p that is up‐regulated in CC (found in this study) and down‐regulated in morbid obese patients.41 miR‐423‐5p regulates DLK1 (down‐regulated in our study), which is known to be involved in skeletal muscle hypertrophy.42 However, the effect of this down‐regulation in CC needs to be further elucidated in future using model systems. miR‐423‐3p, which is derived from the same precursor miR as that of miR‐423‐5p, regulates CAMK2A (down‐regulated in this study), which is involved in calcium signalling. While CAMK2B is associated with inducing WL in CC,4 the role of CAMK2A remains to be elucidated.

miR‐3184‐3p, up‐regulated in the current study, is involved in Wnt/β‐catenin signalling, which plays a role in myogenic differentiation,43 and a defective signalling has been found to have an impact in muscle developmental defects.44 The targets of miR‐3184‐3p include BMPR1B and GREM1 that are up‐regulated and down‐regulated, respectively, in the muscle transcriptome dataset. They are involved in BMP signalling and transforming growth factor beta signalling, deregulation of which may contribute to CC.45

Targets identified for miR‐532‐5p (up‐regulated) in this study include SULF1, RPS6KA6 and NPY1R. While SULF1 plays a role in regulating BMP signalling and is also required for somite development,46 NPY1R is involved in appetite regulation. The action of neuropeptideY (NPY) in energy metabolism is carried out through NPY1R and other NPY receptors,47 and their roles have been widely implicated in the pathogenesis of cachexia.48 RPSKA6 is involved in ciliary neurotrophic factor (CNTF) signalling (see below).

The targets identified for miR‐1296‐5p are HTR2A and RPS6KA6; both these genes were down‐regulated in muscle transcriptome dataset. HTR2A plays a role in serotonin signalling. Altered serotonin signalling has been demonstrated to have an impact on food intake in cachectic mouse models.49 Serotonin was also shown to be involved in myogenesis by promoting longitudinal growth of muscle fibres.50 RPS6KA6 (also regulated by miR‐532‐5p) is involved in CNTF signalling. CNTF is expressed in both central and peripheral nervous systems and has been shown to induce cachexia in model systems, when injected systemically.51

EIF4EBP1 (up‐regulated) was one of the targets identified for miR‐199a‐3p. EIF4EBP1 has been demonstrated to be involved in mammalian target of rapamycin signalling, which regulates muscle protein synthesis. This signalling pathway is shown to be impaired in cachectic states by the action of IL‐6.52 The other targets of miR‐199a‐3p are HTR2A and RPS6KA6.

Overall, based on these insights and inferences, the pleiotropic and redundant roles of miRs are observed against their binding partners regulating different pathways.

Biomarkers identified to‐date in CC, although promising, are not ready for translation to clinic. One of the foremost concerns is the availability of well‐annotated specimens for the discovery stages of the studies in sufficient numbers. Sample sizes used for gene expression studies are critical for confidence in the study findings. While it is important to have large number of samples, it is equally imperative to recognize the difficulties in obtaining biopsy material with detailed clinical annotations for CC phenotype. Despite these challenges, the current study had sample size of 42 for miR profiling and a sample size of 90 for mRNA profiling, well above the sample sizes used in earlier gene expression profiling experiments.4 Nevertheless, validating the identified biomarkers in independent datasets with larger sample size is warranted to confirm the study findings. Although the study has identified several targets and pathways that have potential implications in CC, as judged from the miR‐mRNA correlations, further characterization of the identified targets using model systems is needed to validate the overall biological relevance of the pathways. This would enable us to delineate the role of these identified miRs in CC and also gain functional insights.

Recent studies have highlighted the role of miRs as biomarkers from serum in myotonic dystrophy,53 and it remains to be seen if this could be extended to CC in future. In this proof of principle study, we have demonstrated potential utility of miRs as prognostic and predictive factors from muscle biopsies. However, muscle biopsy is invasive for routine analysis. Therefore, clinical application of the current findings awaits development of alternative methodologies and study designs to make it feasible in clinical setting. To understand the disease trajectory of CC, longitudinal studies can be conducted in plasma and serum as the same may not be possible with muscle biopsies. However, till such time, the challenges and logistics involved to carry out such investigations are addressed, and muscle biopsy remains the best option for CC biomarker and discovery studies. Identification and characterization of miRs that drive CC followed by functional testing in preclinical models may eventually lead to clinical trials to determine the efficacy of certain miRs for CC treatment.

Conclusions

We identified eight novel miRs that could potentially contribute to the aetiology of CC and as promising biomarkers. Continuing efforts in validating the miR signatures in independent cohorts and characterization of the identified miRs to understand its biological relevance to CC may eventually lead to the use of miRs as therapeutics.

Conflict of interest

The author(s) declare that they have no conflict of interest.

Table S1: Technical replicates

Table S2: Comparison of mRNA targets regulated by the differentially expressed miRs from two independent transcriptome datasets

Supporting info item

Click here for additional data file.

Supporting info item

Click here for additional data file.