|Home | About | Journals | Submit | Contact Us | Français|
Conceived and designed the experiments: MD CR DC. Performed the experiments: SU SM IFN GC PO DD. Analyzed the data: LB EC. Contributed reagents/materials/analysis tools: EE SM. Wrote the paper: CR MD SM SU.
Trabectedin, a new antitumor compound originally derived from a marine tunicate, is clinically effective in soft tissue sarcoma. The drug has shown a high selectivity for myxoid liposarcoma, characterized by the translocation t(12;16)(q13; p11) leading to the expression of FUS-CHOP fusion gene. Trabectedin appears to act interfering with mechanisms of transcription regulation. In particular, the transactivating activity of FUS-CHOP was found to be impaired by trabectedin treatment. Even after prolonged response resistance occurs and thus it is important to elucidate the mechanisms of resistance to trabectedin. To this end we developed and characterized a myxoid liposarcoma cell line resistant to trabectedin (402-91/ET), obtained by exposing the parental 402-91 cell line to stepwise increases in drug concentration. The aim of this study was to compare mRNAs, miRNAs and proteins profiles of 402-91 and 402-91/ET cells through a systems biology approach. We identified 3,083 genes, 47 miRNAs and 336 proteins differentially expressed between 402-91 and 402-91/ET cell lines. Interestingly three miRNAs among those differentially expressed, miR-130a, miR-21 and miR-7, harbored CHOP binding sites in their promoter region. We used computational approaches to integrate the three regulatory layers and to generate a molecular map describing the altered circuits in sensitive and resistant cell lines. By combining transcriptomic and proteomic data, we reconstructed two different networks, i.e. apoptosis and cell cycle regulation, that could play a key role in modulating trabectedin resistance. This approach highlights the central role of genes such as CCDN1, RB1, E2F4, TNF, CDKN1C and ABL1 in both pre- and post-transcriptional regulatory network. The validation of these results in in vivo models might be clinically relevant to stratify myxoid liposarcoma patients with different sensitivity to trabectedin treatment.
Trabectedin (Yondelis®) is a tetrahydroisoquinoline molecule originally derived from the Caribbean marine tunicate, Ecteinascidia turbinata, with peculiar in vitro and in vivo cytotoxic activity in a wide range of tumor type: it has been approved by EMEA for the second line therapy of soft tissue sarcomas and ovarian cancer –.
From a mechanistic point of view, trabectedin binds to the minor groove of DNA, differently from the conventional alkylating agents, forming DNA adducts on N2 position of guanine, and bends DNA towards the major groove; two of the three fused rings of this molecule are involved in the minor groove binding while the third (ring C) plays an important role in the process of interaction with different DNA-binding proteins , . It is now widely accepted that trabectedin-DNA interaction affects transcription regulation although the level of selectivity is not completely elucidated , , .
Although trabectedin has shown activity against several types of sarcomas, its higher activity has been observed in patients with myxoid liposarcoma (MLS).
The high sensitivity displayed by MLS tumors is probably correlated to the ability of trabectedin to block the aberrant transcription activity of FUS-CHOP, the chimeric protein hallmark of MLS, generated by the chromosomal translocation t(12;16)(q13;p11) , .
Despite its activity, it is becoming a clinical experience that patients with long lasting trabectedin treatment develop resistance over time. To our knowledge there are only few studies focused on the mechanism of resistance and even fewer are based on a systems biology approach.
To investigate the molecular pathways involved in trabectedin resistance, the 402-91/ET cell line, resistant to trabectedin, was generated from the parental type I myxoid liposarcoma cell line 402-91 by stepwise increasing drug concentrations . Based on the assumption that trabectedin impairs mechanisms of transcription regulation, we exploited several high-throughput technologies to combine different layers of regulation to elucidate the molecular pathways involved in the mechanism of trabectedin resistance in a well defined in vitro model of myxoid liposarcoma, previously described .
The MLS type I 402-91 cell line was kindly supplied by Dr. P. Aman , while the trabectedin resistant 402-91/ET cell line was obtained in our laboratory as previously described . For the array experiments, different batches of exponentially growing 402-91 and 402-91/ET cells were seeded and then total RNA enriched in microRNA (miRNA) fraction and protein lysates were harvested according to standard protocols. Gene, miRNA and protein expression profiles were generated from three independent biological replicates. Each experiment was performed using four technical replicates.
Total RNA enriched in miRNA fraction was purified using the miRNeasy Mini Kit according to the manufacturer's instructions (Qiagen, Milan, Italy). RNA was quantified by spectrophotometry (Nanodrop; Thermo Scientific, Wilmington, USA) and quality assessed by on-chip electrophoresis (Bioanalyser RNA 6000 Nano, Agilent Technologies, Palo Alto, USA). Mature miRNA and mRNA expression levels were examined by qRT-PCR using Applied Biosystems 7900HT. Briefly, cDNA was generated from 1 µg of purified total RNA using the miScript® Reverse Transcription Kit, following the manufacturer's instructions (Qiagen). miRNA expression was quantified using QuantiTect SYBR Green PCR Master Mix and the specific miScript® Primer Assay (Qiagen); RU-6b was used as reference gene. mRNA expression was analyzed using QuantiFast SYBR Green PCR Master Mix (Qiagen) and dedicated primers that were designed as described ; cyclophilin A was used as reference gene. Experiments were run in triplicate for each case to assess technical variability. Analysis was performed by using the 2−ΔΔCt protocol and the median values were compared using the non parametric Mann-Whitney t-test. Differences were considered statistically significant with a two sided p-value<0.05. All tests and data plots were done using GraphPad Prism Version 5.01 (GraphPad Software, La Jolla, Ca, USA).
Total protein extracts were prepared by lysing cell extracts in lysis buffer (50 mM Tris-HCl, pH 8.0, 1% Triton X-100, 150 mM NaCl, 1 mM EGTA, 100 mM NaF, 10% glycerol (100%) and 1 mM MgCl2) with 2 mM phenylmethanesulfonyl fluoride (Sigma Aldrich, St. Louis, MO) as protease inhibitor, in agitation for 30 minutes at 4°C. Insoluble material was pelleted at 13 000 rpm at 4°C and the total protein concentration was determined by Bradford assay. The cellular proteins (40 µg) were separated using 8% and 15% (for detection of HMGA2 protein) SDS-PAGE and transferred electrophoretically to polyvinylidenedifluoride membranes (GE Healthcare, CA, USA). Blots were probed at 4°C overnight with primary antibodies in 5% milk/Tris Buffer Saline Tween20 (TBST). Immunoblotting was carried out with Cyclin D1 (Pharmingen, BD, La Jolla, CA, USA- 554181), E2F5 (abcam, Cambridge, UK - ab59769), SEMA4C (R&D Systems, Minneapolis MN, USA – AF6125), HMGA2 (kindly provided by Prof. Alfredo Fusco, Institute of Experimental Endocrinology and Oncology, CNR Naples, Italy) and PDCD4 (Sigma-Aldrich – SAB1407349) antibodies. Lamin B (Santa Cruz Biotechnology – sc6217) and Actin (Santa Cruz Biotechnology, CA, USA– sc130301) were used for loading control. Antibody binding was revealed by peroxidase labelled secondary antibodies and visualized using enhanced chemiluminescent (ECL) substrate (SuperSignal* West Pico, Thermo Scientific, Rockford IL USA). Experiments were repeated twice.
Gene expression (GE) was analyzed using Agilent slides on which 44,000 oligonucleotides were spotted (Agilent Technologies, Santa Clara, USA). Experiments were performed in four replicates, with dye-swapping duplication-scheme. The arrays were washed and scanned with a laser confocal scanner (G2565BA, Agilent Technologies) according to manufacturer's instructions. Oligo microarrays underwent standard post hybridization processing and the intensities of fluorescence were calculated by Feature Extraction software v.11 (Agilent Technologies). The microarray raw data are deposited in GEO (experiment number GSE20350) according to MIAME (Minimum Information About a Microarray Experiment) guidelines. GE signals were filtered according to the signal-to-background flag; features with more than 60% of signal below the local background were removed from the analysis. Then, lowess normalization was applied . The identification of differentially expressed genes was performed using Empirical Bayes test  implemented in Bioconductor package limma in R software (http://www.r-project.org) with Benjamini and Hochberg False Discovery Rate (FDR)  to control for multiple testing problem. Given the large number of differentially expressed genes, we decided to adopt a more conservative strategy with respect to that applied for differentially expressed miRNAs and proteins: FDR≤0.01 and log fold change greater/lower than 1/-1 were considered to select differentially expressed features. Enrichment analysis of Gene Ontology (GO) categories and KEGG pathways was performed with DAVID web tool . Cytoscape 2.7  was used to reconstruct miRNA-target network.
Transcription factor binding site enrichment analysis was performed using oPOSSUM web tool .
100 ng of total RNA enriched in miRNA fraction was Cy5-labelled and hybridized with a miRNA labelling and hybridization kit according to the manufacturer's instructions (Agilent Technologies). We used the commercially available G4470B human miRNA microarray kit (Agilent Technologies), which consists of 15 K features printed in an 8-plex format (8 of 15 K features printed in known human miRNAs :723 human and 76 human-viral miRNAs) sourced from the Sanger miRBASE public database, release 10.1. The arrays were washed and scanned (with a laser confocal scanner) according to the manufacturer's instructions. miRNA microarrays underwent standard post hybridization processing and the intensities of fluorescence were calculated by Feature Extraction software version 11 (Agilent Technologies). miRNA expression values were filtered as described in the previous section, and vsn algorithm has been applied to normalize the data . miRNA data have been submitted to Array Express (accession number E-MTAB-965).
The identification of differentially expressed miRNAs has been performed using Empirical Bayes test  implemented in Bioconductor package limma of R software with Benjamini and Hochberg FDR to control for multiple testing problem. FDR≤0.1 were considered to select significant miRNAs.
CHOP transcription factor binding site identification on putative miRNA promoter regions was performed using miRGen2.0  and CircuiDB database . Cytoscape 2.7  was used to reconstruct miRNA-target regulatory networks. KEGG database has been used to retrieve gene-gene interactions.
The Panorama Ab Microarray–XPRESS kit (Sigma-Aldrich) used in this study contains 725 different antibodies, spotted in duplicate on nitrocellulose-coated glass slides. These antibodies represent families of proteins known to be involved in a variety of different biological pathways (www.sigma.com/arrays). All reagents and solutions for protein purification were supplied within the kit and experiment was performed on four independent biological replicates, according to manufacturer's instructions. Briefly, total cell lysates were obtained from almost 10 million of exponentially growing 402-91 and 402-91/ET cells. 10 mg of total protein lysate, measured using Bradford protocol, were labeled with Cy3 and Cy5 fluorescent dyes (Amersham Biosciences, Buckinghamshire, UK) and their specific activity and concentration measured at the Nanodrop (Nanodrop, Thermo Scientific) after affinity column purification, provided within the kit. 150 µg of Cy3 and Cy5 labeled samples were co-hybridized on the same slides for 45′ at room temperature, accordingly to manufacturer's instructions. Washing and scanning (GENEPIX, version B) was performed according to the manufacturer's instructions. To increase the robustness of data analysis we performed technical replicates using dye-swap.
Scanned data files were first background subtracted using the normexp method . Positive and negative control spots were removed, and the data were normalized using a modified rank-invariant selection algorithm (in-vMod) .
Determination of differential proteins between the two cell lines was carried out with linear models for microarray analysis , correcting the resulting p-value for multiple testing with FDR. To balance statistical significance and biological relevance, given the relatively low number of spots on the arrays, proteins were selected as differential expressed when their adjusted p-value was less than or equal to 0.1, and they were subsequently annotated with the EntrezGene and RefSeq identifiers of the corresponding genes to facilitate data integration.
Network analysis andGO enrichment of differential proteins was performed using the MetaCore software , using both public data and internal MetaCore curated maps, and with DAVID web tool . Networks and GO terms were deemed significant if their p-value was less than or equal to 0.1.
In order to deeply understand the molecular differences between 402-91 and 402-91/ET cells, we performed multistep omics comparative analyses evaluating gene, miRNA and protein expressions.
A list of 2676 significant genes (corresponding to 3,083 non-unique probes) has been identified as differentially expressed between 402-91/ET and 402-91 cell lines (Figure 1A): 1473 (55%) underexpressed and 1203 (45%) overexpressed. Multiple probes in the microarray platform can be representative of a single gene; in our analysis more than 80% of probes belonging to the same gene behave concordantly. Only 505 genes represented by more than one probe in the platform have been identified as significant by a single probe. The complete list of differentially expressed probes is reported in Table S1. Remarkable is that a significant portion of the underexpressed (6%) and of the overexpressed (2%) genes are zinc finger proteins, i.e. proteins able to bind DNA (p=0.0003 with hypergeometric distribution).
Most of the genes overexpressed in 402-91/ET cell line are implicated in apoptosis, regulation of cell proliferation, angiogenesis, vascular development and negative regulation of developmental processes (enclosed the adipocyte differentiation category). On the other hand, pathways related to metabolism, gene expression and cell cycle progression are mainly populated by downregulated genes. Pathways associated to cancer (p53 signaling pathway and small cell lung cancer) are involved in trabectedin resistance as well as immunoregulatory mechanisms like PPARs and NOD-like receptors signaling. Some key genes involved in cell cycle and regulation of apoptosis have been validated with qRT-PCR (see Table S2 and section below).
Inferential analysis identified a list of 47 miRNAs that showed an altered regulation between 402-91/ET and 402-91 cell lines (Figure 1B): 23 overexpressed and 24 underexpressed (see Table S3). Among these, we found two miRNA clusters (defined as miRNAs whose chromosomal position is less than 1000 kb) on chromosome 13 (miR-17*, miR-19a, miR-19b, miR-20a, miR-92a) and 19 (miR-99b, let-7e, miR-125a-5p) and the miRNA families miR-320, miR-19 and miR-29.
Moreover, we performed a functional enrichment analysis on putative targets of significantly altered miRNAs in 402-91/ET cell line (see Table 3). The in silico predicted targets were cleaned out of those genes without differential expression. Gene targets of miRNAs regulation belong to categories involved in immune response, regulation of cell death, apoptosis and anti-apoptosis and G1/S transition of mitotic cell cycle. Among these genes, there were several known oncogenes playing a key role in tumor progression in melanoma, bladder and lung cancer. Some key miRNAs involved in cell cycle and regulation of apoptosis have been validated with qRT-PCR (see Table S2 and section below).
Apoptosis and cell cycle seem to be the most involved pathways in trabectedin resistance in 402-91/ET cell line either on transcriptional or on translational levels. Thus, to support the robustness of our analysis, differences in expression of genes and miRNA, belonging to the apoptosis and cell cycle pathways were validate in an independent batch of 402-91 and 402-91/ET cells. Representative results for let-7e and miR-21 are reported in Figure 2. As expected, let-7e resulted three folds downregulated (p<0.001) and miR-21 two folds upregulated (p<0.0001) in 402-91/ET compared to 402-91cells. qRT-PCR analysis of let-7e downstream targets confirmed four folds upregulation of CCDN1 (p<0.01), three folds upregulation for SEMA4C (p<0.01) and nine folds upregulation for E2F5 (p<0.001) in the resistant compared to the sensitive cell line (Figure 2A). Western blot analysis confirmed at protein level such differences. Analysis of PDCD4, a downstream target of miR-21, was confirmed to be downregulated in the resistant cells (1.5 times, p<0.0001) at both mRNA and protein levels (Figure 2B). These results confirmed our omics analysis even if the effect of let-7e on HMGA2 needs to be discussed separately. The three folds upregulation of HMGA2 at the mRNA level in 402-91/ET cells (p<0.01) was not further confirmed by western blot analysis (Figure 2A), suggesting that other mechanisms, not directly related to control of mRNA stability by miR-21, are involved in the HMGA2 protein stability.
Protein profiling identified 336 proteins that were significantly (FDR<0.1) different between sensitive and resistant cell lines (Table S4): 148 upregulated and 188 downregulated (see Figure 3). Cluster dendrograms of Figure 1A, 1B and and33 were performed on different data sources (genes, miRNAs and proteins expression levels) characterised by different sizes, means and range of values. Since hierarchical clusteringis highly sensitive to data structure, some differences in dendrograms are expected, given the intrinsic characteristics of the data. According to the experimental design, gene, miRNA and protein expression profiles were quantified using three different biological replicates, while within each batch technical replicates were performed. Thus, the difference in the dendrogram structures is due to the combination of these effects: biological and intrinsic data differences.
In order to assess the biological meaning of the differences in the protein expression, we performed enrichment analysis on GO terms, cellular networks and pathway maps using DAVID web tool and MetaCore software.
MetaCore analysis showed a prevalence of apoptosis-survival related functions for the upregulated proteins: the top three categories from enrichment analysis on GO terms showed prevalence of terms associated to cell death (p<0.0001). Of the top three pathways, two were associated with apoptosis (FAS signaling cascades and TNF family pathways), and “death domain receptors and caspases in apoptosis” was identified as significantly enriched in process networks. On the other hand, GO terms over-represented in downregulated proteins mainly involved cellular component reorganization and affected pathways and networks including AKT signaling, cytoskeleton remodeling and G1/S cell cycle regulation (Table 4). Other involved functions were related to cell cycle, cytoskeletal remodeling and growth factor regulation.
Likewise, the analysis performed with DAVID identified among the highest significant results pathways related to apoptosis, stress induction of HSP and the p53 pathway for upregulated proteins, and cell cycle, AKT signaling, and FAS signaling for downregulated proteins (Table 4).
It has been reported that FUS-CHOP is the main molecular target of trabectedin in myxoid liposarcoma cell line , . Under the assumption that the FUS-CHOP binding site is the same of CHOP (then of the heterodimer CHOP-C/EBPá), we tried to identify putative miRNAs transcriptionally regulated by the chimera and therefore having an enriched CHOP motif in their promoter regions.
Transcription factor binding sites enrichment analysis performed on separately up or downregulated genes (oPOSSUM web tool with Z-score>5 and p-value<0.5) identified respectively 10 and 11 binding sites enriched (see Table S5). CHOP binding site was not identified as enriched.
The identification of miRNA promoter regions is still a challenging task, since putative miRNA regulatory regions could be several thousands of bp upstream of the pri-miR. There are only a limited number of experimentally validated miRNA promoters. In particular only few of the differentially expressed miRNAs have an experimentally validated promoter region. In this context, TFBS enrichment analysis on the promoter regions of the differentially expressed miRNAs would be misleading. However, we could identify miRNAs among those differentially expressed with an experimentally validated promoter with a putative CHOP binding site. Using mirGen2.0 web tool, we found that miR-130a (downregulated in 402-91/ET cell line) and miR-21 (upregulated in 402-91/ET cell line) have a CHOP-C/EBPá binding site in their experimentally validated promoter regions, while using CircuiDB database  we found that miR-7 (upregulated in 402-91/ET cell line) has a CHOP-C/EBPá binding site in its computationally identified promoter region.
miRNA effect on mRNA degradation can be detected by the anti-correlation relationship between miRNA and gene expression. To reconstruct miRNA-mRNA network in this section we focused only on miRNAs having an inverse relationship with their respective mRNA targets.
Using only the list of differentially expressed miRNAs and genes we reconstructed the network in which nodes represent both genes and miRNAs (colored according to their expression levels) while edges represent computationally identified miRNA-target interactions. The entire resulting network is available in Figure S1 and Table S6.
We then focused on specific network modules (sub-networks), corresponding to pathways or specific miRNAs found to be interesting according to the biological problem under study (Figure 4). In particular we selected genes involved in negative regulation of apoptosis (Figure 4A) and cell cycle (Figure 4B). Then we selected the putative targets of those miRNAs (miR-7, miR-21 and miR-130a) with experimentally validated or computationally identified CHOP- C/EBPá binding sites in their promoter regions (Figure 5). Several miRNAs and genes identified as differentially expressed have been validated in qRT-PCR as previously reported (Figure 2; Table S2).
As the effect of miRNAs does not limit to mRNA degradation but also translational repression, we integrated miRNA and protein expression using the approach described previously (see above). As apoptosis and cell cycle were mainly involved in trabectedin resistance as reported in the previous enrichment results (see Table 4), we generated a miRNA-protein network for these two pathways (Figure 6, panels A and B respectively).
The simultaneous analysis of differentially expressed genes and proteins showed a core of 22 genes for which not only the mRNA expression but also the protein levels were significantly altered between 402-91/ET and 402-91 cell lines. Within this core it is possible to clearly identify two sets of genes: 54% of genes with mRNA positively correlated with protein levels (r=0.96) and 45% of genes with mRNA negatively correlated with protein levels (r=−0.6).
This result can be easily explained considering two points: i. mRNA and protein levels are not quantified from exactly the same biological sample but rather from two different biological replicates ii. mRNAs and proteins have different turnover.
miRNA, mRNA and protein can be combined in modules in which a miRNA acts on messenger targets and the resultant effect is measured on the protein.
In our study the combination of post-transcriptional and pre-translational regulatory networks allows the identification of triplets (miRNA-mRNA-protein) representing different miRNA action scenarios. Mapping the levels of expression on these modules allows us to identify putative different types of regulatory loops . The possible loops to be considered regard i) an up/down miRNA regulating the target mRNA through degradation (type A, Figure 7A), ii) an up/down miRNA regulating the protein level through translational repression (type B, Figure 7B), iii) an up/down miRNA whose expression does not support an action on the predicted targets or an up/down miRNA whose effect is not sufficient to overcome the effect of another external signal (and thus is impossible to hypothesize if miRNA acts as mRNA degradation or translational repressor) (type C, Figure 7C).
Depending on the signal, loops can be divided into two classes: coherent (loop of type A; characterized by mutually exclusive expression between miRNA and its targets) and incoherent (loops of type B and C; characterized by coexpression of miRNA and its target) . In the coherent loops both pathways from the signal to the target have the same effect (both repressing or activating target expression), while in the incoherent ones the two pathways have opposite effects . From a biological point of view, incoherent loops of type B have an easy interpretation while incoherent loops of type C have a less clear explanation and then usually require further investigations. Figure 7 shows in details these possibilities jointly with the loops found in our data. In particular we found 38 putative loops of type A, 34 of type B and 28 of type C. Looking at the coherent loops, we identified 3 loops involving differentially expressed miRNAs with roles in drug resistance and/or potentially regulated by FUS-CHOP chimera (let-7e, miR-130a, miR-98 and miR-7). These loops have as potential targets MDC1, TNFRSF21 (overexpressed in 402-91/ET), BACE1 and BAP1 (underexpressed in 402-91/ET). MDC1 activates the intra-S phase and G2/M phase cell cycle checkpoints, in response to DNA damage. It has been demonstrated that MDC1 is implicated in cellular apoptotic response . TNFRSF21 encodes a protein that is a member of the TNF-receptor superfamily. This receptor has been shown to activate NF-kappaB and MAPK8/JNK and induce cell apoptosis. The protein encoded by BAP1 localizes to the nucleus and it interacts with the RING finger domain of the breast cancer 1, early onset protein (BRCA1). This gene is thought to be a tumor suppressor gene that functions in the DNA repair. As far as the incoherent circuits of type B concerns we identified two loops with miR-130a that have PRKD1 and TNF as potential targets and one loop regarding a cluster of miRNAs (miR-320) with TNFRSF21. On the other hand, as far as circuits of type C concerns we identified 3 loops with let-7e and miR-98, FXR2, TUBB4 genes and maintenance of genomic stability pathway.
Despite its complex molecular heterogeneity, cancer results from the integration of at least ten common pathways that drive tumor progression and therapy response . It is now becoming clear that pathway integration rather than single gene manipulation could represent the basis to understand malignant transformation and therapy response.
Based on this hypothesis, in the current study we exploited different transcriptomic (miRNA and GE array) and proteomic (protein array) platforms to generate, in a well defined in vitro model of tumor cells, an integration map of the key pathways possibly involved in the mechanism of resistance to trabectedin.
In this perspective, our study is one of the first attempts towards the combination of the three layers of gene regulation (transcriptional, post-transcriptional and translational) for the identification of those pathways and networks playing a critical role in trabectedin action.
In particular, in this work, we looked for changes in molecular behavior due to an intensive and prolonged treatment with trabectedin that was performed to induce drug resistance. Then, using a systems biology approach we systematically screened the possible sources of variability in mRNAs, miRNAs and protein expressions. Our integration strategy started by extracting individual results from the three different layers and then combining them to increase their power.
Regarding gene expression data, a significant portion of the underexpressed (6%) and overexpressed (2%) genes are zinc finger (ZF) proteins. ZF proteins are multiple and play important role in the metabolism of normal cells. There are many reports on the inhibition of ZF transcription factors in early as well as advanced stages of oncogenesis, including the impairment of signal transduction . These results are also in agreement with the recent findings of Duan et al.  who performed gene expression analysis on chondrosarcoma cell lines resistant to trabectedin. Remarkably, the authors found some ZF proteins to be differentially expressed and suggest that ZF proteins are involved in the mechanism of trabectedin resistance. However, differently from Duan et al.  we did not find overexpression of ZNF93 and ZNF43 in the resistant 402-91/ET cell line.
Interestingly, we found BIRC2 (which encodes for c-IAP1 protein) overexpressed in 402-91/ET cells. c-IAP1 inhibits apoptosis by binding to TNF receptor-associated factors TRAF1 (that we found overexpressed in resistant cells) and TRAF2. It has recently been reported that BIRC2 is associated with resistance of esophageal squamous cell carcinomas to drug-induced apoptosis  and that c-IAP1 could be a novel predictive marker for resistance to radiotherapy in cervical squamous cell carcinomas .
The 402-91/ET cell line seems to be characterized by a reduced apoptosis and an increased cell proliferation. It is noteworthy that pathways involved in cell adhesion are enriched in downregulated genes, suggesting a loss of 402-91/ET intracellular communications and altered cellular processes such as proliferation, migration and differentiation.
miRNAs are involved in processes such as development, carcinogenesis, cell survival, and apoptosis, as well as cellular sensitivity to anticancer drugs . However, opposite effects can be seen towards the same compound in different tumor types, suggesting a complex relationship between miRNAs and drug resistance .
Among differentially expressed miRNAs in 402-91/ET we found miR-21, let-7e, miR-192, miR-130a and miR-98, whose ability to affect the potencies of a number of anticancer agents have been recently reported. Several studies recently reported a relation between the overexpression of miR-21 and the resistance to a variety of anticancer drugs. In particular, Shi and colleagues  demonstrated a role of miR-21 in temozolomide resistance in glioblastoma cells. Moreover, the overexpression of miR-21 has also been found to induce resistance against cytosine arabinoside (Ara-C) and arsenic trioxide. Targeting of PDCD4 (underexpressed in 402-91/ET cells) by miR-21 prevented Ara-C-induced apoptosis in the leukemic HL60 cell line , whereas in several myelogenous cell lines reduced expression of PDCD4 by miR-21 targeting prevented arsenic trioxide-induced apoptosis .
Let-7e, and miR-130a were found underexpressed in a panel of paclitaxel- and cisplatin-resistant cells lines . Furthermore, Herbert et al.  predicted HGMA2 (overexpressed in microarray and validated with qRT-PCR) as a target for miR-98, and showed their involvement in promotion of resistance to doxorubicin and cisplatin, while Boni and colleagues suggested a role of miR-192 in 5-fluorouracil resistance in colorectal cancer .
In general, pathway analysis on microRNA target genes highlights a subset of the pathways previously found in gene expression analysis. This means that in some specific pathways both levels of transcription are co-involved in the development of resistance to trabectedin. Cell cycle, angiogenesis, regulation of transcription and apoptosis-survival are the mostly over-represented pathways.
We investigated the presence of differentially expressed miRNAs, potentially transcriptionally regulated by FUS-CHOP chimera with an enriched CHOP motif in their promoter regions. We found three differentially expressed miRNAs (miR-21, miR-130a, miR-7) having either computationally or experimentally validated CHOP binding sites. Remarkably two of them, miR-21 and miR-130a, were previously discussed because of their documented involvement in drug resistance.
Among our differentially expressed genes we found 149 putative targets of miR-7 (with an anti-correlated expression profile), which are mainly involved in cell cycle, chromatin and cytoskeleton organization. Among them we found PIK3CD, RB1, CAMK2D, CUL5. There are emerging evidences about the involvement of PI3K signaling cascade in myxoid liposarcoma. The overexpression of receptor tyrosine kinases MET, IGF1R and IGF2 promote cell survival through both the PI3K/Akt and the Ras-Raf-ERK/MAPK pathways .
Mir-21 had 99 putative anti-correlated targets mostly involved in the regulation of transcription (CBX4, KLF13, PDCD4 and several ZF proteins).
Among the 118 anti-correlated putative targets of miR-130a, we identified genes involved in vasculature development (ROBO4, SIPR1, NRP1, ELK3, EDN1, THBS1) and cell motion, in particular FOSL1, PLCL2 and LYN. FOSL1 dimerizes with proteins of the JUN family (JUNB is overexpressed in resistant cell line) forming the transcription factor complex AP-1; FOS proteins have been implicated as regulators of cell proliferation, differentiation, and transformation. Furthermore, LYN has been recently shown to act as a mediator of tumor invasion.
Gene and miRNA expression profiles are able to separate 402-91/ET from 402-91 cell lines (Figure 1A and 1B), while protein expression dendrogram (Figure 3) has a more heterogeneous behaviour. 402-91 cell line is divided into two groups; one 402-91/ET sample clusters with 402-91 ones. It is much more difficult to quantify protein expression than gene expression in a multiplex manner than for gene expression, due to the larger variability in the physico-chemical properties of proteins . Therefore, this discrepancy could be due to the higher variability of protein arrays than that of mRNA and miRNA assessment.
GO enrichment analysis on differentially expressed proteins seems to be much more focused than the same analysis on genes. In particular, apoptosis-survival and cell cycle are categories capable of collecting the greatest number of differences between sensitive and resistant cells. Among the over-represented pathways we found p53, apoptosis and MAPK.
In particular we found CCND1, CCND2, BBC3, CDKN2A, PERP involved in p53 signaling, RB1, E2F2, E2F3, E2F4, SMAD4, CDC14, CDC25A, CDC6, CDC7 involved in cell cycle. Some of these proteins are involved in the regulation of transcription in G1/S phase of the mitotic cell cycle. G1/S transition is a rate-limiting step in cell cycle progression .
The cyclin D1 proto-oncogene is a key regulator of G1 to S phase progression in several cell types. Together with CDK4 and CDK6, cyclin D1 forms active complexes that promote cell cycle progression by phosphorylating and inactivating the retinoblastoma protein (RB1). RB1 in turn leads to the release of the E2F family of transcription factors driving the expression of several genes associated with S phase progression . A number of therapeutic agents have been observed to induce cyclin D1 degradation in vitro, indicating that such an induction may offer a useful avenue for therapeutic intervention.
After the analysis and the discussion of each regulatory aspect (mRNAs, miRNAs and proteins) separately, we combined all the biological evidences within a systems biology approach in order to have a wider perspective of the regulatory phenomena acting during trabectedin resistance.
miRNAs act as post-transcriptional regulators of gene expression. The reduction of protein levels is the final result of the regulation that can be achieved through two strategies: (i) mRNA degradation or (ii) translational repression. The effects of the two modes of action can be evaluated at different points of the transcription-translation pathway. Focusing on the two main highly correlated pathways that seem to be involved in trabectedin resistance (negative regulation of apoptosis and cell cycle), we reconstructed two different networks: i) the network composed of differentially expressed mRNAs targeted by differentially expressed miRNAs that allow us to study the regulation at the transcriptional level, hereafter the “post-transcriptional regulatory network” and ii) the network composed by differentially expressed miRNAs leading to differentially expressed proteins, hereafter “pre-translational regulatory network”.
It is clear the central role that CCND1, RB1, E2F4, TNF, CDKN1C, ABL1, TNFRSF21 play in both post-transcriptional and pre-translational regulatory networks. In particular, they act as a bridge between the two parts of the network with opposite regulations: upregulated genes/protein and downregulated miRNA and downregulated genes/protein and upregulated miRNA.
In some cases the involvement of miRNAs can explain the discordant relation between the transcriptome and the proteome. In a seminal review Inui et al.  describe the role of miRNAs in pathways as decision maker elements to discriminate real versus too weak or too transient signals. In this study we identified 3 loops involving differentially expressed miRNAs with roles in drug resistance; these loops have as potential targets MDC1 (overexpressed in 402-91/ET), BACE1 and BAP1 (underexpressed in 402-91/ET).
The present study is one of the first works that try to screen the behavior of genes, microRNAs (miRNAs) and proteins in order to reconstruct in silico the regulatory networks leading to resistance to anticancer agents such as trabectedin. Diverse sources of omic data have been used to describe altered regulatory mechanisms within a transcriptional, post-transcriptional and translational perspective. The combination of the three layers of regulation gave us the possibility to reconstruct putative gene-miRNA-protein circuits altered in trabectedin resistance. We found that transcriptome and proteome data agreed in recognizing anti-apoptosis and cell cycle proliferation as the two main biological processes characterizing the difference between 402-91 and 402-91/ET cell lines. All the hallmarks of cancer seem to be altered in the resistant cell line, raising the hypothesis that the mechanism of resistance is acting with more strength on the same molecular markers triggered by cancer development.
On the other hand, post-transcriptional analysis revealed two miRNAs (miR-21 and miR-130a) already identified as playing a role in drug resistance and putatively regulated by the FUS-CHOP chimera, which could be good markers for future functional studies. Most of miR-21 and miR-130a targets are known oncogenes or oncosuppressors highlighting their involvement in malignancies.
Once validated on in vivo models this approach might allow the identification of druggable targets for trabectedin-resistant malignancies.
miRNA and mRNA network reconstruction using mirSVR as target prediction tool. Colors represent expression (green and red for under and overexpressed genes, light green and light red for under and overexpressed miRNA) and edges between miRNA and gene represent mirSVR predictions, while gene-gene edges represent validated gene-gene interactions from KEGG database.
Complete list of probes identifying differentially expressed genes.
Complete list of validated mRNA and miRNA expression levels using qRT-PCR.
Complete list of differentially expressed miRNAs.
Complete list of differentially expressed proteins.
Complete list of predicted transcription factor binding sites enriched in the promoter region of up and downregulated gene in 402-91/ET cell lines.
The authors acknowledge Prof. Luca Tamagnone (Associate Professor, University of Torino, IRCC) for generously supporting set up SEMA4C western blot analysis.
Competing Interests: The authors have declared that no competing interests exist.
Funding: This work was supported by University of Padova (CPDR075919 to CR). Furthermore, the authors want to acknowledge the CINECA Award N. HP10BDJ9X8, 2010 for the availability of high performance computing resources and support. The authors acknowledge the generous contribution of the Italian Association for Cancer Research (AIRC). IFN is recipient of a fellowship from ACTO foundation. PO and GC are supported by Compagnia di San Paolo, Torino, Italy. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.