ANRIL, a long noncoding RNA (lncRNA), has recently been reported to have a direct role in recruiting polycomb repressive complexes PRC2 and PRC1 to regulate the expression of the p15/CDKN2B-p16/CDKN2A-p14/ARF gene cluster. Expression analysis of ANRIL, EZH2, SUZ12, EED, JARID2, CBX7, BMI1, p16, p15, and p14/ARF genes was evaluated in a large cohort of invasive breast carcinomas (IBC, n = 456) by qRT-PCR and immunohistochemistry (IHC) was performed on CBX7, EZH2, p14, p15, p16, H3K27me3, and H3K27ac. We observed significant overexpression in IBCs of ANRIL (19.7%) and EZH2 (77.0%) and an underexpression of CBX7 (39.7%). Correlations were identified between these genes, their expression patterns, and several classical clinical and pathologic parameters, molecular subtypes, and patient outcomes, as well as with proliferation, epithelial–mesenchymal transition, and breast cancer stem cell markers. Multivariate analysis revealed that combined EZH2/CBX7 status is an independent prognostic factor (P = 0.001). In addition, several miRNAs negatively associated with CBX7 underexpression and EZH2 overexpression. These data demonstrate a complex pattern of interactions between lncRNA ANRIL, several miRNAs, PRC2/PRC1 subunits, and p15/CDKN2B-p16/CDKN2A-p14/ARF locus and suggest that their expression should be considered together to evaluate antitumoral drugs, in particular the BET bromodomain inhibitors.

Implications: This study suggests that the global pattern of expression rather than expression of individual family members should be taken into account when defining functionality of repressive Polycomb complexes and therapeutic targeting potential. Mol Cancer Res; 14(7); 623–33. ©2016 AACR.

This article is featured in Highlights of This Issue, p. 587

Invasive breast carcinomas (IBC) are highly heterogeneous tumors, and triple-negative subtype (TNC) remains a major cause of death in women. Although more than 90% of breast tumors are clinically localized, about 50% of them relapse within 5 years. Thus, identification of new prognostic factors and potential targeted therapies is crucial to improve outcome of patients with IBC.

Changes in the epigenetic landscape, including histone modifications, DNA methylation, and noncoding RNAs (ncRNA,) altered expression are now considered as a new hallmark of cancer (1). Polycomb repressive complexes (PRC) play important roles in chromatin remodeling and inhibition of transcription, suppressor genes silencing, and interconnection with major signaling pathways (2). Polycomb group protein subunits are essential regulators of embryonic stem cell pluripotency and early developmental cell fate decisions that are often deregulated in cancer. Two multiprotein complexes, PRC2 and PRC1, act concertedly and sequentially in transcriptional repression via two distinct histone modifications, trimethylation of lysine 27 on histone H3 (H3K27me3) and monoubiquitination of histone H2A (H2AK119ub). PRC2 contains mainly four subunits (EZH2, SUZ12, EED, and JARID2) and is involved in initiation of silencing by catalyzing methylation of histone H3 (3). PRC1 also comprises mainly four subunits (BMI1, PHC, CBX7, and RING1) and is implicated in maintenance of silencing by catalyzing monoubiquitination of histone H2A (ref. 4; Supplementary Fig. S1).

PRC2 and PRC1 have been described to interact with long noncoding RNAs (lncRNA). LncRNAs are defined by a length ranged from 200 bp to 100 kbp and their number is in constant increase: 48,680 lncRNAs are at present identified (5–7). LncRNAs are important regulators coordinating expression of protein-coding genes nearby (cis-regulation) or at distance (trans-regulation). LncRNAs function as guide, tethers, scaffolds, and decoys and their mechanisms of action are located at epigenetics, translational, and posttranslational levels (8). Accumulating data suggest that deregulation of lncRNAs is pivotal in cancer initiation, progression, and metastatic spread. Among important lncRNAs deregulated in cancer, ANRIL (antisense noncoding RNA in the INK4 locus) is encoded in the chromosome 9p21 region (9) and has been reported to have a direct role in recruiting PRC2 and PRC1 complexes to specific loci and in repressing gene expression (10). ANRIL is transcribed as a 3.8-kb lncRNA in the opposite direction from the p15/CDKN2B-p16/CDKN2A-p14/ARF gene cluster. Common disease genome-wide association studies (GWAS) have identified the ANRIL gene as a shared genetic susceptibility locus to coronary disease, intracranial aneurysm, type 2 diabetes, and numerous cancers (11). Increased ANRIL expression levels were observed in prostate carcinomas and involved in repression of the CDKN2B/CDKN2A/ARF gene cluster in cis by directly binding to PRC complexes (10). EZH2 (PRC2) and CBX7 (PRC1) act as transcriptional repressors of many genes and are particularly implicated in silencing via ANRIL of the CDKN2B/CDKN2A/ARF locus (12). However, the functional role and underlying mechanism of action of ANRIL in breast carcinogenesis still remain unclear (13).

The aim of the present study was to analyze at RNA and protein levels, the ribonucleoproteic network composed of the lncRNA ANRIL, the two polycomb complexes PRC2/PRC1, and the downstream CDKN2B/CDKN2A/ARF gene cluster, as well as several miRNAs targeting PRC units, in IBCs and to evaluate their clinical significance.

Patients and samples

Samples of 456 primary unilateral invasive primary breast tumors excised from women managed at Institut Curie (Saint-Cloud, France) from 1978 to 2008 have been analyzed. All patients cared in our institution before 2007 were informed that their tumor samples might be used for scientific purposes and had the opportunity to decline. Since 2007, patients treated in our institution have given their approval by signed informed consent. This study was approved by the local ethics committee (Breast Group of Rene Huguenin Hospital). Samples were immediately stored in liquid nitrogen until RNA extraction. A tumor sample was considered suitable for this study if the proportion of tumor cells exceeded 70%. The patients (mean age, 61.7 years; range, 31–91 years) all met the following criteria: primary unilateral nonmetastatic breast carcinoma for which complete clinical, histologic, and biologic data were available; no radiotherapy or chemotherapy before surgery; and full follow-up at Institut Curie. Treatment consisted of modified radical mastectomy in 283 cases (63.9%) and breast-conserving surgery plus locoregional radiotherapy in 160 cases (36.1%). The patients had a physical examination and routine chest radiography every 3 months for 2 years and then annually. Mammograms were done annually. Adjuvant therapy was administered to 369 patients, consisting of chemotherapy alone in 91 cases, hormone therapy alone in 176 cases, and both treatments in 102 cases. The histologic type and the number of positive axillary nodes were established at the time of surgery. The malignancy of infiltrating carcinomas was scored according to Scarff Bloom Richardson (SBR) histoprognostic system. Estrogen receptor (ERα), progesterone receptor (PR), and human EGF receptor 2 (ERBB2) status was determined at the protein level by using biochemical methods [dextran-coated charcoal method, enzyme immunoassay, or immunohistochemistry (IHC)] and confirmed by real-time quantitative RT-PCR assays (14, 15). The population was divided into 4 groups according to HR (ERα and PR) and ERBB2 status, as follows: two luminal subtypes [HR+ (ERα+ or PR+)/ERBB2+ (n = 54)] and [HR+ (ERα+ or PR+)/ERBB2− (n = 290)]; an ERBB2+ subtype [HR− (ERα− and PR−)/ERBB2+ (n = 45)]; and a triple-negative subtype [HR− (ERα− and PR−)/ERBB2− (n = 69)]. Standard prognostic factors are shown in Supplementary Table S1. Patients (n = 169) developed metastasis during a median follow-up of 8.9 years (range, 6 months to 29 years). Ten specimens of adjacent normal breast tissue from patients with breast cancer (n = 7) and normal breast tissue from women undergoing cosmetic breast surgery (n = 3) were used as sources of normal RNA.

RNA extraction

Total RNA was extracted from breast tumor samples by using acid-phenol guanidium as previously described (16). RNA quality was determined by electrophoresis through agarose gels, staining with ethidium bromide, and visualization of the 18S and 28S RNA bands under ultraviolet light.

Real-time RT-PCR

LncRNA and protein-coding genes expression analysis.

ANRIL, EZH2, SUZ12, EED, JARID2, CBX7, BMI1, P16/CDKN2A, P15/CDKN2B, and P14/ARF gene mRNA expression levels were quantified by using real time RT-PCR in a retrospective series of 456 well-characterized tumors. Quantitative values were obtained from the cycle number (Ct value) at which the increase in the fluorescence signal associated with exponential growth of PCR products started to be detected by the laser detector of the ABI Prism 7900 Sequence Detection System (Perkin-Elmer Applied Biosystems), using PE Biosystems analysis software according to the manufacturer's manuals. The precise amount of total RNA added to each reaction mix (based on optical density) and its quality (i.e., lack of extensive degradation) are both difficult to assess. We therefore also quantified transcripts of the TBP gene (Genbank accession: NM_003194) encoding the TATA box–binding protein (a component of the DNA-binding protein complex TFIID) as an endogenous RNA control and normalized each sample on the basis of its TBP content. We selected TBP as an endogenous control because the prevalence of its transcripts is moderate and because there are no known TBP retropseudogenes (retropseudogenes lead to coamplification of contaminating genomic DNA and thus interfere with RT-PCR, despite the use of primers in separate exons; ref. 14). Results, expressed as N-fold differences in target gene expression relative to the TBP gene and termed "Ntarget," were determined as Ntarget = 2ΔCt_sample, where the ΔCt value of the sample was determined by subtracting the average Ct value of the target gene from the average Ct value of the TBP gene. The Ntarget values of the samples were subsequently normalized such that the median of the Ntarget values for the 10 normal breast tissues was 1. The primers for TBP, ANRIL, and the 9 protein-coding genes were chosen with the assistance of the Oligo 6.0 program (National Biosciences). We scanned the dbEST and nr databases to confirm the total gene specificity of the nucleotide sequences chosen for the primers and the absence of SNPs. To avoid amplification of contaminating gDNA, 1 of the 2 primers was placed at the junction between 2 exons or on 2 different exons. Agarose gel electrophoresis was used to verify the specificity of PCR amplicons. The conditions of cDNA synthesis and PCR were as described (17).

miRNAs expression analysis.

miRNAs were isolated with the RNA extraction procedure used for the protein-coding transcripts (total RNA extraction). Reverse transcription was performed with the QIAGEN miScript Reverse Transcription Kit, according to the manufacturer's protocol (QIAGEN, GmbH). Specific miRNAs were quantified by real-time PCR with the QIAGEN miScript SYBR Green PCR Kit (QIAGEN). The small nucleolar RNU44 was used as an internal control. The relative expression level of each miRNA, expressed as N-fold difference in target miRNA expression relative to RNU44, and termed “Ntarget”, was calculated as follows: Ntarget = 2ΔCt_sample. The ΔCt value of a given sample was determined by subtracting the Ct value of the target miRNA from the Ct value of RNU44. The Ntarget values of the samples were subsequently normalized so that the median Ntarget value of samples with low level of EZH2 and CBX7 transcripts was 1.

IHC

We performed IHC assay by using CBX7 (Novus, polyclonal, rabbit, NBP1-79042, 1/250, pH6), EZH2 (Novocastra, monoclonal mouse, 6A10, 1/100, pH6), p15 (Abcam, polyclonal, rabbit, ab53034, 1/50, pH6), p16 (CINtec p16, Ventana 705-4713, mouse, monoclonal, E6H4 prediluted, 1 μg/mL), H3K27me2/3 [Abcam, anti-histone H3 (di- + tri-methyl K27), ab6147, mouse, monoclonal, 1/200, pH6], and H3K27ac [Abcam, ab 4729, Rabbit polyclonal to Histone H3 (acetyl K27), 1/400, pH6] antibodies in a series of 80 IBCs. Paraffin-embedded tissue blocks, obtained at the time of the initial diagnosis, were retrieved from the archives of the Department of Biopathology, Institut Curie. Sections of 3 μm were cut with a microtome from the paraffin-embedded tissue blocks of normal breast tissue, preinvasive lesions, and IBCs. Tissue sections were deparaffinized and rehydrated through a series of xylene and ethanol washes. Briefly, key figures included: (i) antigen retrieval in 0.1 mol/L citrate buffer, pH 6 (Biocare) in a pressure cooker (4 minutes); (ii) blocking of endogenous peroxidase activity by immersing sections in 3% hydrogen peroxide in methanol for 15 minutes and subsequently rinsing them in water and PBS; (iii) incubation with primary antibodies against the targeted antigen; (iv) immunodetection with a biotin-conjugated secondary antibody formulation that recognizes rabbit and mouse immunoglobulins, followed by peroxidase-labeled streptavidin and linking with a rabbit biotinylated antibody against mouse immunoglobulin G (DAKO SA); and (v) chromogenic revelation with AEC and counterstaining with Mayer's hematoxylin. All immunostaining was processed by using a DAKO automated immunostaining device. The specificity of the antibodies was confirmed by doing IHC studies with the same protocol on paraffin-embedded human tissue sections containing lymphocytes. A semiquantitative histologic score (H-score = intensity × frequency) was performed (score 0 = negative staining, score 1 = weak staining, score 2 = moderate staining, score 3 = strong staining).

Statistical analysis

The distributions of target mRNA levels were characterized by their median values and ranges. Relationships between mRNA levels of the different target genes and between mRNA (and protein) levels and clinical parameters were identified by using nonparametric tests, namely, the χ2 test (relation between 2 qualitative parameters), the Mann–Whitney U test (relation between 1 qualitative and 1 quantitative parameter), and the Spearman rank correlation test (relation between 2 quantitative parameters). Differences were considered significant at confidence levels greater than 95% (P < 0.05). Metastasis-free survival (MFS) was determined as the interval between initial diagnosis and detection of the first metastasis. Survival distributions were estimated by the Kaplan–Meier method, and the significance of differences between survival rates was ascertained with the log-rank test. The Cox proportional hazards regression model was used to assess prognostic significance, and the results are presented as HRs and 95% confidence intervals (CI). Hierarchical clustering was performed with GenANOVA software (18).

ANRIL and EZH2 are upregulated, whereas CBX7 is downregulated in IBCs

We measured mRNA levels of the lncRNA ANRIL and 6 genes belonging to the two repressive polycomb PRC2 and PRC1 complexes (PRC2: EZH2, SUZ12, EED, and JARID2; PRC1: BMI1 and CBX7) in a series of 456 patients with unilateral invasive breast tumors and in a series of 10 normal breast tissues (including adjacent normal breast tissue from patients with breast cancer and normal breast tissue from women undergoing cosmetic breast surgery) by using quantitative RT-PCR assays (Table 1). We did not observe significant differences concerning the expression profiles of the 7 genes tested between the adjacent normal breast tissue from patients with breast cancer (7 samples) and normal breast tissue from women undergoing cosmetic breast surgery (3 patients; data not shown).The mRNA values of the breast cancer samples were normalized such that the median of the 10 normal breast tissue mRNA values was 1. To determine the cutoff point for altered target genes expression in breast cancer tissues, the Ntarget value, calculated as described in Materials and Methods, was determined for the 10 normal breast RNA samples. As these normal values were consistently between 0.39 (CBX7) and 2.45 (ANRIL), values of 0.33 or less were considered to represent underexpression, and values of 3 or more to represent overexpression of these genes in tumor samples. We have previously used the same approach to determine cutoff points for tumor gene altered expression (14, 19). In our series of 456 breast carcinomas, we observed overexpression of ANRIL in 19.7% of IBCs (median = 1.58; minimum = 0.05; maximum = 13.5). Concerning Polycomb subunits, we observed a marked overexpression of EZH2 in 77.0% of IBCs (median = 5.09; minimum = 0.76; maximum = 48.3) and to a lesser extent of BIM1 (median = 1.23; minimum = 0.00; maximum = 14.0) in 10.1% of IBCs and an unexpected underexpression of CBX7 in 39.7% of IBCs (median = 0.38; minimum = 0.00; maximum = 2.59). SUZ12, EED, and JARID2 were overexpressed only in 2.6%, 1.8%, and 5.1% of IBCs, respectively.

We performed an IHC analysis by using anti-EZH2 and anti-CBX7 antibodies in a series of 80 IBCs (standard prognostic factors are shown in Supplementary Table S2). We identified a moderate-to-intense nuclear staining of tumor cells (H-score = 2–3) with EZH2 antibodies in 75.9% of IBCs (n = 63 of 80) and an absence or a slight nuclear and cytoplasmic staining of tumor cells (H-score = 0–1) with CBX7 antibodies in 63.8% of IBCs (n = 51 of 80; Fig. 1).

The p15/CDKN2B-p16/CDKN2A-p14/ARF cluster is unexpectedly not silenced in IBCs

We measured mRNA levels of p15/CDKN2B, p16/CDKN2A, and p14/ARF in the series of 456 patients. We showed marked overexpression of p16, p14/ARF, and p15, respectively, in 15.8%, 39.9%, and 21.7% of IBCs (Table 1).

We studied protein expression of p14, p16, and p15 in IBCs by using IHC in our series of 80 breast tumors. We observed nuclear and cytoplasmic overexpression of cancer cells (H-score = 2–3) with anti-p14, p16, and p15 antibodies in, respectively, 41.2% (n = 33 of 80), 51.3% (n = 41 of 80), and 37.5% (n = 30 of 80) of IBCs (Supplementary Fig. S2).

PRC2/PRC1 activity status revealed decrease of H3K27me3 expression levels and increase of H3K27ac expression levels in IBCs

To evaluate PRC2/PRC1 repressive activity status in breast tumors, we performed an IHC analysis in our series of 80 IBCs by using anti-H3K27me3 and anti-H3K27ac antibodies. We could identify an absence or a slight nuclear staining of cancer cells (H-score = 0–1) with H3K27me3 antibody in 76.3% of IBCs (n = 61 of 80) and a moderate-to-intense nuclear staining of tumor cells (H-score = 2–3) with H3K27ac antibody in 73.8% of IBCs (n = 59 of 80; Supplementary Fig. S3). An expected negative correlation was observed between staining of H3K27me3 and H3K27ac (P = 0.0000073).

ANRIL, PRC2/PRC1 subunits and p15/CDKN2B-p16/CDKN2A-p14/ARF cluster expression levels are associated with clinicopathologic parameters, molecular subtypes, and patient outcome

Clinicopathologic parameters.

ANRIL mRNA expression level was exclusively and weakly correlated with ERα (P = 0.038) and PR (P = 0.011) status (Supplementary Table S3).

Concerning the PRC2 complex, EZH2 mRNA expression level was strongly associated with SBR histologic grade (P < 10−7), ERα status (P = 3 × 10−7), and PR status (P = 3 × 10−7), moderately associated with ERBB2 status (P = 0.00020) and tumor size (P = 0.0062; Supplementary Table S4). Expression of SUZ12 was exclusively associated with lymph node status (P < 10−7; Supplementary Table S5), EED with SBR histologic grade (P = 0.0011), and PR status (P = 0.0028; Supplementary Table S6) and JARID2 with ERα status (P = 0.0027) and ERBB2 status (P = 0.0035; Supplementary Table S7).

Concerning the PRC1 complex, CBX7 mRNA expression level was highly significantly correlated with SBR histologic grade (P = 0.000025), ER status (P = 0.000011), and PR status (P = 0.0000013; Supplementary Table S8). BMI1 was weakly associated with ERα status (P = 0.036) but not with any of other clinical predictive factors (Supplementary Table S9).

Finally, ARF/p14 and p16 mRNA expression levels were highly correlated with SBR histological grade (P < 10−7 and P = 0.000019, respectively), ERα status (P = 1.8 × 10−7 and P < 10−7, respectively), and PR status (P = 1.3 × 10−5 and P = 1.8 × 10−6, respectively; Supplementary Table S10), whereas p15 was moderately correlated with ERα status (P = 0.009) and PR status (P = 0.0008), but not with SBR histologic grade, like ARF/p14 and p16 (Supplementary Tables S11 and S12).

It is also noteworthy that none of these 10 gene expressions showed significant link with the adjuvant treatment status including hormone therapy alone (176 patients), chemotherapy alone (91 patients), both hormone therapy and chemotherapy (102 patients), and no adjuvant therapy (87 patients; Supplementary Tables S3–S12).

Breast cancer molecular subtypes.

By using HR (ERα and PR) and ERBB2 status, we subdivided our total population of IBCs (n = 456) into 4 subgroups: HR+/ERBB2+ (n = 54), HR+/ERBB2− (n = 289), HR−/ERBB2+ (n = 45), and HR−/ERBB2− (n = 68). ANRIL, PRC1/2 genes, and p14/p15/p16 mRNA levels of expression in the different breast cancer molecular subtypes are shown in Supplementary Table S13. A strong correlation was observed between molecular subgroups and ANRIL (P = 0.00013), EZH2 (P < 10−7), JARID2 (P = 0.00090), CBX7 (P = 0.000020), and p14/ARF and p16 (P < 10−7) mRNA levels. High frequencies of ANRIL, p14/ARF, and p16 deregulated expressions (under- or overexpressions) were mainly observed in triple-negative subtypes, EZH2 overexpression, and CBX7 underexpression in the two HR− subtypes, and JARID2 overexpression in HR−/ERBB2+ subtypes (Supplementary Table S13).

Patient outcome.

We assessed the impact of variations of ANRIL, EZH2, SUZ12, EED, JARID2, CBX7, BMI1, p14/ARF, p15, and p16 mRNA levels on patient outcome. We used a log-rank test to identify relations between MFS and mRNA levels of these 10 genes (Table 2). EZH2 overexpression and CBX7 underexpression were significantly associated with shorter MFS (P = 0.000045 and P = 0.0051, respectively). We observed also a light significant difference in MFS among patients with p16 overexpression (P = 0.025). MFS was not influenced by expression status of ANRIL and other genes analyzed (i.e. SUZ12, EED, JARID2,BMI1, ARF, and p15).

Two molecular markers may provide a more accurate prediction of patient outcome when combined than when considered in isolation. By combining EZH2 and CBX7 status, we identified four separate prognostic groups with significantly different MFS curves (P = 0.00011; Fig. 2). The patients with the poorest prognosis had EZH2 overexpression and CBX7 underexpression [5-year MFS: 61.7% (57.7%–65.7%); 10-year MFS: 51.4% (47.2%–55.6%)], whereas those with the best prognosis had EZH2 and CBX7 normal expression [5-year MFS: 86.4% (82.4%–90.4%); 10-year MFS: 81.2% (76.5%–85.9%)].

Multivariate analysis (Cox proportional hazards model) showed that lymph node status, tumor size, and combined EZH2/CBX7 status were independent variables predictive of MFS (P = 0.00021, P = 0.014, and P = 0.001, respectively; Supplementary Table S14).

There is a complex interplay between ANRIL, PRC1/PRC2 subunits and p15/CDKN2B-p16/CDKN2A-p14/ARF genes cluster

By using the Spearman rank correlation test for continuous variables, we could find a strong positive correlation between mRNA expression levels of the following genes belonging to the ANRIL pathway: (i) ANRIL and both ARF (r = +0.594, P < 10−7), p15 (r = +0.381, P < 10−7), and p16 (r = +0.323, P < 10−7); (ii) EZH2 and both EED (r = +0.295, P < 10−7), JARID2 (r = +0.278, P < 10−7), and BMI1 (r = +0.229, P = 0.0000031); and (iii) SUZ12 and CBX7 (r = +0.317, P < 10−7), reflecting very complex interplay across this gene set (Supplementary Table S15).

Hierarchical clustering of ANRIL pathway genes according to their expression patterns in the series of 456 IBCs represented in dendrogram format, subdivided genes into two subgroups (Supplementary Fig. S4). Group I included ANRIL, p14/ARF, p15, p16, EZH2, EED, JARID2, and BMI1 and contained 4 genes, ANRIL and the p15/CDKN2B-p16/CDKN2A-p14/ARF gene cluster, highly related. Group II included SUZ12 (from the PRC2 complex) and CBX7 (from the PRC1 complex).

Correlations of ANRIL and PRC2/PRC1 genes expression levels with epithelial–mesenchymal transition, proliferation, and breast cancer stem cell markers

Data from the literature have suggested that Polycomb complexes are implicated in proliferation, epithelial–mesenchymal transition (EMT) induction, and stem cell differentiation. In consequence, we tested the possible links between ANRIL pathway gene expression and mRNA levels of various markers of proliferation (i.e., MKI67), of EMT (i.e., VIM and TWIST1), and of breast cancer stem cell (i.e., ALDH1A1, ALDH1A3, and CD133) in our series of 456 IBCs (Table 3). Marked positive associations were observed between CBX7 and two EMT markers (VIM and TWIST1; P < 10−7 and P < 10−7, respectively). A significant positive correlation was also identified between CBX7 and one breast cancer stem cell marker (ALDH1A1; P < 10−7). EZH2 was highly linked to MKI67 (P < 10−7), whereas EED was linked to two breast cancer stem cell markers: ALDH1A3 and CD133 (P < 10−7 and P = 0.000082, respectively).

As previous studies suggested a role of ANRIL in EMT (13), we further investigated the correlations between the expression levels of ANRIL and 9 different biomarkers of EMT in a series of 60 breast carcinomas (30 tumors with high expression level for ANRIL and 30 tumors with low expression level for ANRIL). We confirmed the absence of significant associations between ANRIL and various EMT-associated markers including TWIST1, ZEB1, ZEB2, SNAI1/Snail, SNAI2/Slug, VIM/Vimentin, CDH2/N-cadherin, CDH1/E-cadherin, and ZO-1 (Supplementary Table S16).

miRNAs are implicated in posttranscriptional regulation of PRC1/PRC2

Previous studies suggested complex regulations between Polycomb subunit genes and miRNAs in carcinomas (20, 21). To investigate roles of miRNAs in posttranscriptional deregulation of the ANRIL/PRC pathway in IBCs, we selected 17 miRNAs from the literature. Among them, 4 miRNAs (miR-26A1, miR-101, miR-125B2, and miR-214) are known to target directly EZH2 transcript (22–24) and 2 miRNAs (miR-181B1, miR-181A2) to target directly CBX7 transcript (25). The 11 others miRNAs also target putatively EZH2 and CBX7 (26). We analyzed expression of these 17 miRNAs, using real-time quantitative RT-PCR, in IBCs with normal (n = 20) and marked low (n = 20) levels of CBX7 transcripts and with normal (n = 20) and marked high (n = 20) levels of EZH2 transcripts. We showed that (i) miR-181B1, miR-181A2, and miR-203 levels were marked significantly increased (P < 0.01) in CBX7-underexpressed tumors (Table 4) and that (ii) miR-26A1 and miR-125B1 levels were significantly decreased in EZH2-overexpressed tumors (Table 5).

These results of microRNA levels, shown in Tables 4 and 5, were standardized versus RNU44 as endogenous microRNA control. Similar results (significantly links between microRNA and CBX7 and EZH2 expressions) were obtained using a second well-known endogenous microRNA control (i.e., miR-191).

Among the major discovered lncRNAs deregulated in carcinogenesis, ANRIL has been shown to regulate its neighbor tumor suppressor genes CDKN2B, CDKN2A, and ARF by epigenetic mechanisms and thereby control cell proliferation and senescence (27). ANRIL, acting in cis, specifically binds various Polycomb proteins, in particular EZH2 (PRC2) and CBX7 (PRC1), resulting in histone modification in the CDKN2A/CDKN2B/ARF locus. EZH2 is a histone methyltransferase and a member of PRC2 that not only methylates histone H3 on Lysine 27 but also interacts with and recruits DNA methyltransferases to methylate CpG of certain EZH2 target genes to establish firm repressive chromatin structures, contributing to tumor progression. CBX7 is a Chromobox family protein and a member of the PRC1 and is directly involved in regulation of genes frequently silenced in cancer.

To assess expression of the lncRNA ANRIL and the Polycomb subunits of PRC2/PRC1 as well as CDKN2A/CDKN2B/ARF locus in IBCs, we used qRT-PCR in a large series of 456 IBCs from patients with known clinical and pathologic status and long-term outcome. Evaluation at mRNA level was combined with an IHC study in a series of 80 IBCs. Indeed, these two techniques of expression analysis (qRT-PCR and IHC) give complementary information, the technique of qRT-PCR with an extraction in homogeneous solution quantified transcripts levels both in tumor cells and in stroma cells, while IHC study allows only a semiquantitative analysis of protein levels, but on an individual cell basis (allowing distinction between epithelial tumor cells and stroma cells).

In the present study, we observed ANRIL overexpression in 19.7% of IBCs. ANRIL depletion has been associated with reduced proliferation (28), suggesting a procarcinogenic role supported by overexpression of this lncRNA in various others cancer types (29–31). In our IBC series, ANRIL overexpression was observed more particularly in the aggressive TNC subtype.

We also observed a marked overexpression of EZH2 and to a lesser extent of BMI1, respectively, in 77.0% and 10.1% of IBCs and CBX7 underexpression in 39.7% of the tumors. We confirmed these data at the protein level in the series of 80 IBCs by identifying a moderate and intense nuclear staining with an anti-EZH2 antibody in 82% of IBCs and an absence or a slight nuclear and cytoplasmic staining of cancer cells with a CBX7 antibody in 43% of IBCs. In malignant tumors, ANRIL is a natural antisense transcript which usually recruits the two Polycomb repressor complexes PRC2 and PRC1, resulting in chromatin reorganization with silencing of the INK4A/ARF/INK4B locus. For example, positive links were observed between CBX7 overexpression, ANRIL overexpression, and CDKN2A/CDKNB/ARF underexpression in prostate cancer (10). Conversely, in our series of breast tumors, we showed marked CBX7 underexpression and a nonexpected high positive link between ANRIL overexpression and CDKN2A/CDKNB/ARF overexpression, and more particularly with p14/ARF. Although these results seem to be contradictory, it can be hypothesized that CBX7 can exhibit both oncosuppressive and oncogenic functions, depending on type of cancer, microenvironment, and presence of interacting proteins (32–35).

CBX7 underexpression was strongly correlated with high SBR histologic grade, negative ERα and PR status, and EMT and breast cancer stem cell markers. Apart from few exceptions (as prostate cancer), CBX7 expression is lost in human malignant neoplasias, and a clear correlation between its downregulated expression and a cancer aggressiveness and poor prognosis has been observed (36).

EZH2 overexpression was also strongly associated with high SBR histologic grade and negative ERα and PR status. But, unlike CBX7, EZH2 overexpression was highly linked to the proliferation marker MKI67, suggesting a different role of these two proteins in breast tumorigenesis (37). Numerous studies point on the widespread roles of PRC2 (as well as H3K27 methylation) in developmental and differentiation processes of multicellular organisms, and on their implication in fundamental chromatin mechanisms that underlie stem cell regulatory circuits and cancer progression (38). Thus, it is not surprising that increasing evidences indicate that EZH2 deregulation is frequently observed in a variety of cancers (39–41). EZH2 overexpression is mainly found in solid tumors (as in our present study in breast cancer), whereas activating or inactivating mutations are identified in hematologic malignancies (42).

Using log-rank analysis to identify relations between MFS and mRNA levels of the 10 ANRIL pathway genes, only EZH2 overexpression and CBX7 underexpression were significantly associated with shorter MFS. Multivariate analysis showed that combined EZH2/CBX7 status was an independent prognostic factor.

Results from previous studies suggest complex interactions between Polycomb complexes and miRNAs in cancer progression. As no mutation or allelic loss of CBX7 located in 22q13.1 were identified in large genomic studies (43) and focal amplification of EZH2 located in 7q36 was only observed in 1% of IBCs (43), we suspected epigenetic mechanisms implicating miRNAs could account for CBX7 underexpression and EZH2 overexpression in IBCs. Among the 17 miRNAs selected, (i) miR-181B1, miR-181A2, and miR-203 expressions were significantly increased in IBCs with CBX7 underexpression and (ii) miR-26A1 and miR-125B2 expressions were significantly decreased in IBCs with EZH2 overexpression, in agreement with the literature (44, 45). Our results thus showed that EZH2 is targeted by 3 oncosuppressive miRNAs (miR-26A1, miR-125B, and miR-214) that are decreased in IBCs, suggesting an epigenetic mechanism of EZH2 overexpression in IBCs. Concerning CBX7, our results were also consistent with recent data in the literature, suggesting that the 2 oncomiRs miR-181B1 and miR-181A2 could be induced by the oncogene HMGA1, resulting in CBX7 underexpression (46). Indeed, we observed a positive link between HMGA1 mRNA levels and expression of both miR-181B1 (r = +0.340, P = 0.0038) and miR-181A2 (r = +0.264, P = 0.025) in 71 tumors of our IBC series (data not shown).

Recent data revealed that histone modifier complexes PRC2 and PRC1 have a bivalent role in cancer, bearing both oncogenic and tumor-suppressive properties depending on cell tumor type (47). Complex deficiency in PRC2 and PRC1 components could decrease expression levels of epigenetic repressive mark H3K27 trimethylation, thereby contributing to H3K27 acetylation and resulting in increased transcription and activation of oncogenic pathways. Thus, defining a PRC2/PRC1 activity status in IBCs may be pivotal in combining or reconsidering therapeutic options. We identified a moderate-to-intense nuclear staining of tumor cells with H3K27ac antibodies in 73.8% of IBCs associated to an absence or a slight nuclear staining of cancer cells with H3K27me3 antibody in 76.3% of IBCs. Collectively, these data suggest that loss of Polycomb repressive function in IBCs might be due to epigenetic alteration of PRC1 and PRC2 by several microRNA targeting CBX7 and EZH2, resulting in altered H3K27trimethylation and H3K27 acetylation. Development of therapies targeting epigenetic processes in cancer has gained increasing focus with study of HDAC and DNA methyltransferase inhibitors and more recently of BET bromodomain inhibitors (48). Indeed, recent studies of malignant peripheral nerve sheath tumors (MPNST) have showed that alterations of PRC2 subunits lead to the amplification of Ras-driven transcription and confer sensitivity to BET bromodomain inhibitors (49, 50).

In conclusion, our results revealed a complex pattern of interactions between ANRIL, PRC2/PRC1 and suppressive and oncogenic miRNAs in IBCs. Our findings also point to EZH2 and CBX7 as the most promising of the 10 genes investigated for clinical applications and therapeutic approaches to breast cancer. This study suggest that it is the global pattern of expression, rather than expression of individual family members, that should be taken into account when defining functionality of repressive Polycomb complexes and evaluating antitumor drugs targeting these genes. Thus, establishment of a PRC2/PRC1 activity status in IBCs could be contributive in targeting BET proteins by using bromodomain inhibitors when levels of H3K27me3 are decreased and H3K27ac increased. Moreover, inhibition of oncomiRNAs miR-181B1 and miR-181A2 (as well as restoration of oncosuppressive miR-26A1, miR-125B, and miR-214) may represent another attractive combined strategy for breast cancer–targeted therapy.

No potential conflicts of interest were disclosed.

Conception and design: D. Meseure, E. Pasmant, I. Bieche

Development of methodology: D. Meseure, S. Vacher, I. Bieche

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): D. Meseure, S. Vacher, K. Drak Alsibai, A. Nicolas, I. Bieche

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): D. Meseure, S. Vacher, K. Drak Alsibai, W. Chemlali, R. Lidereau, I. Bieche

Writing, review, and/or revision of the manuscript: D. Meseure, E. Pasmant, C. Callens, I. Bieche

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): K. Drak Alsibai, M. Caly, C. Callens, I. Bieche

Study supervision: I. Bieche

This work was supported by grant INCa-DGOS-4654.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

1.
Piao
HL
,
Ma
LJ.
Non-coding RNAs as regulators of mammary development and breast cancer
.
J Mammary Gland Biol Neoplasia
2012
;
17
:
33
42
.
2.
Kerppola
TK
. 
Polycomb group complexes: many combinations, many functions
.
Trends Cell Biol
2009
;
19
:
692
704
.
3.
Margueron
R
,
Reinberg
D
. 
The Polycomb complex PRC2 and its mark in life
.
Nature
2011
;
469
:
343
9
.
4.
Luis
NM
,
Morey
L
,
Di Croce
L
,
Benitah
SA
. 
Polycomb in stem cells: PRC1 branches out
.
Cell Stem Cell
2012
;
11
:
16
21
.
5.
Iyer
MK
,
Niknafs
YS
,
Malik
R
,
Singhal
U
,
Sahu
A
,
Hosono
Y
et al 
The landscape of long noncoding RNAs in the human transcriptome
.
Nat Genet
2015
;
47
:
199
208
.
6.
Derrien
T
,
Johnson
R
,
Bussotti
G
,
Tanzer
A
,
Djebali
S
,
Tilgner
H
, et al
The GENCODE v7 catalog of human long noncoding RNAs: analysis of their gene structure, evolution, and expression
.
Genome Res
2012
;
22
:
1775
89
.
7.
Harrow
J
,
Frankish
A
,
Gonzalez
JM
,
Tapanari
E
,
Diekhans
M
,
Kokocinski
F
, et al
GENCODE: the reference human genome annotation for The ENCODE Project
.
Genome Res
2012
;
22
:
1760
1774
.
8.
Wang
KC
,
Chang
HY
. 
Molecular mechanisms of long noncoding RNAs
.
Mol Cell
2011
;
43
:
904
14
.
9.
Pasmant
E
,
Laurendeau
I
,
Héron
D
,
Vidaud
M
,
Vidaud
D
,
Bièche
I
. 
Characterization of a germ-line deletion, including the entire INK4/ARF locus, in a melanoma-neural system tumor family: identification of ANRIL, an antisense noncoding RNA whose expression coclusters with ARF
.
Cancer Res
2007
;
67
:
3963
9
.
10.
Yap
KL
,
Li
S
,
Munoz-Cabello
AM
,
Raguz
S
,
Zeng
L
,
Mujtaba
S
, et al
Molecular interplay of the noncoding RNA ANRIL and methylated histone H3 lysine 27 by polycomb CBX7 in transcriptional silencing of INK4a
.
Mol Cell
2010
;
38
:
662
74
.
11.
Pasmant
E
,
Sabbagh
A
,
Masliah-Planchon
J
,
Ortonne
N
,
Laurendeau
I
,
Melin
L
, et al
Role of noncoding RNA ANRIL in genesis of plexiform neurofibromas in neurofibromatosis type 1
.
J Natl Cancer Inst
2011
;
103
:
1713
22
.
12.
Popov
N
,
Gil
J
. 
Epigenetic regulation of the INK4b-ARF-INK4a locus: in sickness and in health
.
Epigenetics
2010
;
5
:
685
90
.
13.
Li
CH
,
Chen
Y
. 
Targeting long non-coding RNAs in cancers: progress and prospects
.
Int J Biochem Cell Biol
2013
;
45
:
1895
910
.
14.
Bieche
I
,
Onody
P
,
Laurendeau
I
,
Olivi
M
,
Vidaud
D
,
Lidereau
R
, et al
Real-time reverse transcription-PCR assay for future management of ERBB2-based clinical applications
.
ClinChem
1999
;
45
:
1148
56
.
15.
Bieche
I
,
Parfait
B
,
Laurendeau
I
,
Girault
I
,
Vidaud
M
,
Lidereau
R
. 
Quantification of estrogen receptor alpha and beta expression in sporadic breast cancer
.
Oncogene
2001
;
20
:
8109
15
.
16.
Chomczynski
P
,
Sacchi
N
. 
Single-step method of RNA isolation by acid guanidiniumthiocyanate-phenol-chloroform extraction
.
Anal Biochem
1987
;
162
:
156
9
.
17.
Bieche
I
,
Parfait
B
,
Le Doussal
V
,
Olivi
M
,
Rio
MC
,
Lidereau
R
, et al
Identification of CGA as a novel estrogen receptor-responsive gene in breast cancer: an outstanding candidate marker to predict the response to endocrine therapy
.
Cancer Res
2001
;
61
:
1652
8
18.
Didier
G
,
Brezellec
P
,
Remy
E
,
Henaut
A
. 
GeneANOVA-gene expression analysis of variance
.
Bioinformatics
2002
;
18
:
490
1
.
19.
Bieche
I
,
Noguès
C
,
Lidereau
R
. 
Overexpression of BRCA2 gene in sporadic breast tumours
.
Oncogene
1999
;
18
:
5232
8
.
20.
Cao
Q
,
Mani
RS
,
Ateeq
B
,
Dhanasekaran
SM
,
Asangani
IA
,
Prensner
JR
, et al
Coordinated regulation of polycomb group complexes through microRNAs in Cancer
.
Cancer Cell
2011
;
20
:
187
99
.
21.
Ye
S
,
Yang
L
,
Zhao
X
,
Song
W
,
Wang
W
,
Zheng
S
. 
Bioinformatics method to predict two regulation mechanism: TF-miRNA-mRNA and lncRNA-miRNA-mRNA in pancreatic cancer
.
Cell Biochem Biophys
2014
;
5
:
41
9
.
22.
Penna
E
,
Orso
F
,
Taverna
D
. 
miR-214 as a key hub that controls cancer networks: small player, multiple functions
.
J Invest Dermatol
2015
;
135
:
960
9
.
23.
Zhang
B
,
Liu
XX
,
He
JR
,
Zhou
CX
,
Guo
M
,
He
M
, et al
Pathologically decreased miR-26a antagonizes apoptosis and facilitates carcinogenesis by targeting MTDH and EZH2 in breast cancer
.
Carcinogenesis
2011
;
32
:
2
9
.
24.
Sun
YM
,
Lin
KY
,
Chen
YQ
. 
Diverse functions of miR-125 family in different cell contexts
.
J Hematol Oncol
2013
;
6
:
6
14
.
25.
Sepe
R
,
Formisano
U
,
Federico
A
,
Forzati
F
,
Bastos
AU
,
D'Angelo
D
, et al
CBX7 and HMGA1b proteins act in opposite way on the regulation of the SPP1 gene expression
.
Oncotarget
2015
;
6
:
2680
92
.
26.
O'Loghlen
A
,
Muñoz-Cabello
AM
,
Gaspar-Maia
A
,
Wu
HA
,
Banito
A
,
Kunowska
N
, et al
MicroRNA regulation of Cbx7 mediates a switch of Polycomb orthologs during ESC differentiation
.
Cell Stem Cell
2012
;
10
:
33
46
.
27.
Congrains
A
,
Kamide
K
,
Ohishi
M
,
Rakugi
H
. 
ANRIL: molecular mechanisms and implications in human health
.
Int J Mol Sci
2013
;
14
:
1278
92
.
28.
Chen
D
,
Zhang
Z
,
Mao
C
,
Zhou
Y
,
Yu
L
,
Yin
Y
, et al
ANRIL inhibits p15 (INK4b) through the TGFβ1 signaling pathway in human esophageal squamous cell carcinoma
.
Cell Immunol
2014
;
289
:
91
6
.
29.
Hua
L
,
Wang
CY
,
Yao
KH
,
Chen
JT
,
Zhang
JJ
,
Ma
W
. 
High expression of long non-coding RNA ANRIL is associated with poor prognosis in hepatocellular carcinoma
.
Int J Clin Exp Pathol
2015
;
8
:
3076
82
.
30.
Lin
L
,
Gu
ZT
,
Chen
WH
,
Cao
KJ
. 
Increased expression of the long non-coding RNA ANRIL promotes lung cancer cell metastasis and correlates with poor prognosis
.
Diagn Pathol
2015
;
10
:
14
21
.
31.
Qiu
JJ
,
Lin
YY
,
Ding
JX
,
Feng
WW
,
Jin
HY
,
Hua
KQ
. 
Long non-coding RNA ANRIL predicts poor prognosis and promotes invasion/metastasis in serous ovarian cancer
.
Int J Oncol
2015
;
46
:
2497
505
.
32.
Shinjo
K
,
Yamashita
Y
,
Yamamoto
E
,
Akatsuka
S
,
Uno
N
,
Kamiya
A
, et al
Expression of chromobox homolog 7 (CBX7) is associated with poor prognosis in ovarian clear cell adenocarcinoma via TRAIL-induced apoptotic pathway regulation
.
Int J Cancer
2014
;
135
:
308
18
.
33.
Zheng
X
,
Zhou
J
,
Zhang
B
,
Zhang
J
,
Wilson
J
,
Gu
L
, et al
Critical evaluation of Cbx7 downregulation in primary colon carcinomas and its clinical significance in Chinese patients
.
BMC Cancer
2015
;
15
:
145
55
.
34.
Karamitopoulou
E
,
Pallante
P
,
Zlobec
I
,
Tornillo
L
,
Carafa
V
,
Schaffner
T
, et al
Loss of the CBX7 protein expression correlates with a more aggressive phenotype in pancreatic cancer
.
Eur J Cancer
2010
;
46
:
1438
44
.
35.
Pallante
P
,
Terracciano
L
,
Carafa
V
,
Schneider
S
,
Zlobec
I
,
Lugli
A
, et al
The loss of the CBX7 gene expression represents an adverse prognostic marker for survival of colon carcinoma patients
.
Eur J Cancer
2010
;
46
:
2304
13
.
36.
Pallante
P
,
Forzati
F
,
Federico
A
,
Arra
C
,
Fusco
A
. 
Polycombprotein family member CBX7 plays a critical role in cancer progression
.
Am J Cancer Res
2015
;
5
:
1594
601
.
37.
Reijm
EA
,
Timmermans
AM
,
Look
MP
,
Meijer-van Gelder
ME
,
Stobbe
CK
, et al
High proteinexpression of EZH2 is related to unfavorable outcome to tamoxifen in metastatic breast cancer
.
Ann Oncol
2014
;
25
:
2185
90
.
38.
Wu
J
,
Crowe
DL
. 
The histone methyltransferase EZH2 promotes mammary stem and luminal progenitor cell expansion, metastasis and inhibits estrogen receptor-positive cellular differentiation in a model of basal breast cancer
.
Oncol Rep
2015
;
34
:
455
60
.
39.
Chen
S
,
Huang
L
,
Sun
K
,
Wu
D
,
Li
M
,
Li
M
, et al
Enhancer of zeste homolog 2 as an independent prognostic marker for cancer: a meta-analysis
.
PLoS One
2015
;
10
:
e0125480.
40.
Morin
RD
,
Johnson
NA
,
Severson
TM
,
Mungall
AJ
,
An
J
,
Goya
R
, et al
Somatic mutations altering EZH2 (Tyr641) in follicular and diffuse large B-cell lymphomas of germinal-center origin
.
Nat Genet
2010
;
42
:
181
5
.
41.
Sudo
T
,
Utsunomiya
T
,
Mimori
K
,
Nagahara
H
,
Ogawa
K
,
Inoue
H
, et al
Clinicopathological significance of EZH2 mRNA expression in patients with hepatocellular carcinoma
.
Br J Cancer
2005
;
92
:
1754
8
.
42.
Völkel
P
,
Dupret
B
,
Le Bourhis
X
,
Angrand
PO
. 
Diverse involvement of EZH2 in cancer epigenetics
.
Am J Transl Res
2015
;
7
:
175
93
.
43.
Koboldt
DC
,
Fulton
RS
,
McLellan
MD
,
Schmidt
H
,
Kalicki-Veizer
J
,
McMichael
JF
, et al
Comprehensive molecular portraits of human breast tumours
.
Nature
2012
;
490
:
61
70
.
44.
Derfoul
A
,
Juan
AH
,
Difilippantonio
MJ
,
Palanisamy
N
,
Ried
T
,
Sartorelli
V
. 
Decreased microRNA-214 levels in breast cancer cells coincides with increased cell proliferation, invasion and accumulation of the Polycomb Ezh2 methyltransferase
.
Carcinogenesis
2011
;
32
:
1607
14
.
45.
Xiangxiang
W
,
Xiaoyun
D
,
Shengcan
C
,
Haojun
S
,
Haizhong
J
,
Ying
F
, et al
The functional sites of miRNAs and lncRNAs in gastric carcinogenesis
.
Tumor Biol
2015
;
36
:
521
32
.
46.
Mansueto
G
,
Forzati
F
,
Ferraro
A
,
Pallante
P
,
Bianco
M
,
Esposito
F
, et al
Identification of a new pathway for tumor progression: microRNA-181b up-regulation and CBX7 down-regulation by HMGA1 protein
.
Gene Cancer
2010
;
1
:
210
24
.
47.
Lee
W
,
Teckie
S
,
Wiesner
T
,
Ran
L
,
Prieto Granada
CN
,
Lin
M
, et al
PRC2 is recurrently inactivated through EED or SUZ12 loss in malignant peripheral nerve sheath tumors
.
Nat Genet
2014
;
46
:
1227
32
.
48.
Fu
LL
,
Tian
M
,
Li
X
,
Li
JJ
,
Huang
J
,
Ouyang
L
, et al
Inhibition of BET bromodomains as a therapeutic strategy for cancer drug discovery
.
Oncotarget
2015
;
6
:
5501
16
.
49.
Baude
A
,
Lindroth
AM
,
Plass
C
. 
PRC2 loss amplifies Ras signalingin cancer
.
Nat Genet
2014
;
46
:
1154
5
.
50.
De Raedt
T
,
Beert
E
,
Pasmant
E
,
Luscan
A
,
Brems
H
,
Ortonne
N
, et al
PRC2 loss amplifies Ras-driven transcription and confers sensitivity to BRD4-based therapies
.
Nature
2014
;
514
:
247
51
.

Supplementary data