Genome-wide miRNA response to anacardic acid in breast cancer cells

MicroRNAs are biomarkers and potential therapeutic targets for breast cancer. Anacardic acid (AnAc) is a dietary phenolic lipid that inhibits both MCF-7 estrogen receptor α (ERα) positive and MDA-MB-231 triple negative breast cancer (TNBC) cell proliferation with IC50s of 13.5 and 35 μM, respectively. To identify potential mediators of AnAc action in breast cancer, we profiled the genome-wide microRNA transcriptome (microRNAome) in these two cell lines altered by the AnAc 24:1n5 congener. Whole genome expression profiling (RNA-seq) and subsequent network analysis in MetaCore Gene Ontology (GO) algorithm was used to characterize the biological pathways altered by AnAc. In MCF-7 cells, 69 AnAc-responsive miRNAs were identified, e.g., increased let-7a and reduced miR-584. Fewer, i.e., 37 AnAc-responsive miRNAs were identified in MDA-MB-231 cells, e.g., decreased miR-23b and increased miR-1257. Only two miRNAs were increased by AnAc in both cell lines: miR-612 and miR-20b; however, opposite miRNA arm preference was noted: miR-20b-3p and miR-20b-5p were upregulated in MCF-7 and MDA-MB-231, respectively. miR-20b-5p target EFNB2 transcript levels were reduced by AnAc in MDA-MB-231 cells. AnAc reduced miR-378g that targets VIM (vimentin) and VIM mRNA transcript expression was increased in AnAc-treated MCF-7 cells, suggesting a reciprocal relationship. The top three enriched GO terms for AnAc-treated MCF-7 cells were B cell receptor signaling pathway and ribosomal large subunit biogenesis and S-adenosylmethionine metabolic process for AnAc-treated MDA-MB-231 cells. The pathways modulated by these AnAc-regulated miRNAs suggest that key nodal molecules, e.g., Cyclin D1, MYC, c-FOS, PPARγ, and SIN3, are targets of AnAc activity.

For miRNA RNA-seq The Truseq Small RNA kit (Illumina) was used to prepare miRNA libraries from 1 μg total RNA. Each Library was individually gel purified on a Novex TBE 6% gel and resuspended in 10uL 10mM Tris-Cl, pH 8.5. Libraries were validated and quantitated by running 1μL on the Agilent Technologies 2100 Bioanalyzer DNA High Sensitivity Chip. 36-cycle single sequencing reads were generated on the Illumina NextSeq500 instrument utilizing the 500 Mid-output v2 (75 cycle) sequencing kit. The resulting samples were divided into 48 FASTQ [15] single-end raw sequencing files representing four conditions: MCF-7 control, MCF-7 treated with AnAc 24:1n5 (MCF-7 AnAc), MDA-MB-231 control, and MDA-MB-231 treated with AnAc 24:1n5 (MDA MB-231 AnAc). These raw data of our RNA-seq are available at Gene Expression Omnibus (GEO) database: accession number GSE78011.

Differential miRNA expression analysis
A total of three biological replicates for each treatment were analyzed, with four flow cell lanes per replicate. Raw sequence data files were downloaded from Illumina's BaseSpace (https:// basespace.illumina.com/) onto the KBRIN server for analysis the miRDeep2 [16] and edgeR [17]. Each of the four single-end raw. FASTQ files for each replicate (representing the four flow cells) was concatenated into one single-end. FASTQ file using the unix cat command. Quality control (QC) of the raw sequence data was performed using FastQC (version 0.10.1) [18]. The FastQC results indicated sequence trimming was not necessary since the minimum quality value for all samples was well above Q30 (1 in 1000 error rate) (data not shown).
Given that this is a miR sequencing project, preliminary adapter trimming was performed on each of the samples using a custom file adaptersToTrim.fa which contains a subset of the Illumina TruSeq Small RNA adapter and primer sequences taken from https://support. illumina.com/content/dam/illumina-support/documents/documentation/chemistry_ documentation/experiment-design/illumina-adapter-sequences_1000000002694-00.pdf Sequences were trimmed of the adapters with Trimmomatric v0.33 [19]. The trimmed sequences were directly aligned to the human hg19 reference genome assembly using the mapper.pl wrapper of the miRDeep2 package (v 0.0.7) [16]. This script used bowtie (version 1.1.1) [20], generating alignment files in arf format. The aligned sequences were then used as inputs into the miRDeep2 package and the script quantifier.pl. In addition, this script used the mirBase release 21 [21] mature miRNA and miRNA hairpin sequences downloaded from ftp://mirbase.org/pub/mirbase/CURRENT/. The result was a file containing the number of reads mapping to each of the 2,822 human (hsa) miRs for the specific sample. After quantification, the resulting counts for each miR in each sample were combined into a reads matrix. This was accomplished using a custom perl script, createReadMatrix.pl. Differentially expressed miRs were determined using edgeR [17] and a customized R script, Schultz-Klinge. miRNA.R. Using a p-value cutoff of 0.05, the number of differentially expressed miRs in each comparison is shown in Table 1.

In silico network analysis
We performed pathway and network analysis of differentially expressed genes in MetaCore™ version 6.27 (GeneGO, Thomson Reuters, New York, N.Y.). MetaCore™ is a web-based software suite for multiple applications in systems biology including RNA-seq analysis as used here. MetaCore™ analyses are based on MetaBase (http://metadatabase.org/), a 100% manually-curated integrated database of mammalian biology that contains over 6 million experimental findings on protein-protein, protein-DNA, protein-RNA, and protein-compound interactions; metabolic and signaling pathways; and other information [22]. Generation of heatmaps: Files of miRNAs significantly altered by AnAc treatment in each cell line were imported into Partek software Version 6.6 (Partek Inc., St Louis, MO.) and Partek Genomic Suite™ was used to generate heatmaps (Fig 1, S1 and S2 Figs). Each hierarchical clustering was created using Euclidean distance as similarity measure for genes and samples. We noted that one of the three MCF-7 AnAc samples appeared to behave as a hybrid between the other two AnAc treated and three control (EtOH)-treated samples (S2 Fig). RNA isolation, RT-PCR and quantitative real-time PCR (qPCR) of miRNAs and mRNAs Cell growth, treatment and RNA isolation and quantification/quality assessment were performed as described above. For miRNA, RNA was converted to cDNA using the Taqman 1 miRNA Reverse Transcription kit (PE Applied Biosystems). For mRNA, RNA was converted to cDNA using the High Capacity cDNA Reverse Transcription kit (PE Applied Biosystems). Primers for hsa-miR-268g, hsa-miR-612, hsa-miR-20b-5p, and hsa-miR-20b-3p were purchased from TaqMan (Advanced miRNA assays) and RNU48 (TaqMan) was used as the reference for normalization [23]. Primers for VIM (Vimentin) [24]: Forward 5'-GACAATGCGTC TCTGGCACGTCTT-3'; Reverse 5'-TCCTCCGCCTCCTGCAGGTTCTT-3'; for ZFP36L1 (ZFP36 Ring Finger Protein Like 1, aka ERF1 and BRF1) [25]: Forward, 5 0 -AGGATGACCAC CACCCTCGTGTCT-3 0 , Reverse, 5 0 -CCC CCTGCACTGGGAGCACTA-3 0 , and for GAPDH [26] were purchased from IDT. qPCR was performed using ABI Viia 7 (Life Technologies) with each reaction run in triplicate. The comparative threshold cycle (Ct) method (2 -ΔΔCT ) was used to determine fold change relative to vehicle treated or control transfected cells [27].

Transient transfection
MCF-7 and MDA-MB-231 cells were transiently transfected for 24 h with miR-612 mimic, miR-612 inhibitor, Anti-miR ™ negative control #1, or mirVANA™ miRNA mimic negative control #1 (all from Ambion, Life Technologies, Thermo Fisher Scientific, Carlsbad, CA, USA), using Lipofectamine RNAiMAX transfection reagent (Invitrogen, Thermo Fisher Scientific) and Opti-MEM 1 Reduced Serum Medium (Invitrogen, Thermo Fisher Scientific). After 24 h of transfection, cells were treated with ethanol (EtOH, vehicle control) or 13.5 or 35 μM AnAc, for MCF-7 and MDA-MB-231 respectively, in phenol red-free IMEM medium containing 5% DCC-stripped FBS for 48 h prior to MTT assay (CellTiter 96, Promega, Madison, WI, USA). Two separate experiments were performed with quadruplicate wells within each experiment. For analysis of miR-612 expression in transfected cells, the medium was changed 24 h Table 1. Differentially expressed miRNAs (DEmiRs). The log2-fold change with zero value in the control conditions was arbitrarily set to one and the maximum log2-fold change value and those with zero value in the treatment conditions were arbitrarily set to the minimum log2-fold change value of minus one. The number of differentially expressed genes in each comparison is shown and the number of upregulated genes indicated with the upward arrow and downregulated genes indicated by downward arrow. after transfection as above, without any treatment and RNA was harvested (see above) a total of 72 h post transfection, i.e., at the same time the MTT assay was performed for qPCR of miR-612 using RNU48 as a control (see above).

Results and discussion
RNA-seq analysis of AnAc-regulated miRNAs  Table 1. Tables 2-5 list the AnAc-regulated miRNAs in MCF-7 and MDA-MB-231 cells, their genomic location and host gene (if applicable), information about their relevance in breast or other cancers and their experimentally verified, i.e., bona fide, targets. The expression of more miRNAs was significantly changed in response to AnAc in MCF-7 cells vs MDA-MB-231 cells (Figs 1 and 2). The heatmap shows that MCF-7 and MDA-MB-231 cells have different responses to AnAc with MDA-MB-231 cells showing less change in response to AnAc compared with MCF-7 cells (Fig 1). These data suggest that AnAc selectivity alters miRNA transcript expression in these two cell lines through mostly non-overlapping mechanisms.
As shown in the Venn Diagrams of  Table 2). The common GO Processes for upregulated miR-20b and miR-612 were identified by Meta-Core™ analysis and listed in Fig 2; however no matches between genes/proteins for miR-20b and miR-612 were identified in Pathway Maps by MetaCore analysis. Interestingly, AnAc increased miR-20-3p in MCF-7 and miR-20-5p in MDA-MB-231 cells. This suggests that distinct miR-20b targets would be expected to be regulated in response to AnAc upregulation of miR-20b-3p versus miR-20-5p in the two cell lines. The selection of which mature miRNA 5p or 3p arm is dominant is determined by thermodynamic and structural properties of the processed pre-miR-duplex AGO protein (reviewed in [32]). The functional consequences of arm selection are therefore distinct. The exact mechanism of miRNA Induced Silencing Complex (miRISC) assembly remains elusive and includes a human miRNA loading complex containing the ds-pre-miRNA, DICER1, TRBP2 and miRNA-free AGO protein as its components. [33]. Recent studies in Huh7 human hepatoma cells showed that an increase in target genes, i. e., SLC7A1 (CAT-1), increased the processing of pre-miR-122 to miR-122, implying that increases in target mRNA levels can promote miRNA biogenesis [33]. Whether this is true for other cells and miRNAs remains to be examined. The MetaCore network enrichment analysis of the miRNAs upregulated in AnAc-treated MCF-7 vs. MDA-MB-231 cells identified "Cellular response to inorganic substance" as the top GO process (S3A Fig). The network analyses for miR-20b and miR-612 are shown in S3A and S3B Fig. There is only one previous examination of miRNAs, mRNAs, and lncRNAs in MCF-7 and MDA-MB-231 cells, but that study used a microarray expression profiling [167] rather than an unbiased RNA-sequencing approach. None of the AnAc-regulated miRNAs was among the miRNAs more highly expressed in MCF-7 compared with MDA-MB-231 cells [167]. In contrast, miR-4284 was more highly expressed in MDA-MB-231 cells [167] and we observed that AnAc decreased miR-4284 in MDA-MB-231 cells (

miRNAs downregulated by AnAc in MCF-7 cells
Twenty-one miRNAs were downregulated by AnAc in MCF-7 cells (Table 3). miRNAs are encoded within a gene (intronic or exonic) or are intergenic (reviewed in [168]). miRNAs can be regulated independently or are cotranscribed with their host gene (reviewed in [8]). To examine if the miRNA host gene was downregulated by AnAc in MCF-7 cells we searched miR-520a-5p Chr19, intergenic. miR-520a-3p inhibits proliferation by targeting HOXD8 in non-small cell lung cancer None experimentally validated for 5p. For 3p: CCND1 and CD44 [74] miR-520d-5p Chr19, intergenic. involved in HER2-receptor-related differentiation through undefined mechanisms [75]. Overexpression by lentiviral-miR-520d infection of human HLF and Huh7 hepatoma cells converted the cells to non-tumorigenic and less differentiated normal stem cells, but no miRNA target genes were validated [76]. Acts as a tumor suppressor in colorectal cancer [77].
(Continued) Chr11, host gene MTRNR2L8. Is transported into mitochondria and inhibits 16S rRNA processing and mitochondrial protein synthesis [112]. Acts as a tumor suppressor in MCF-7 cells in vitro and in MDA-MB-231 cells in xenograft studies in mice [112].

miR-1247-5p
Chr14, in the DLK1-DIO3 genomic imprinted microRNA cluster [149]. Downregulated in aromatase-resistant MCF-7 breast cancer cells [103] and lung adenocarcinomas [150]. Acts as a tumor suppressor in pancreatic cancer [151]. Silenced by DNA methylation in lung adenocarcinomas and cell lines and overexpression promotes apoptosis and inhibits cell invasion and migration [152]. Overexpressed in castrationresistant prostate cancer [153]. GSE78011. In AnAc-treated MCF-7 cells, six downregulated host genes for downregulated miRNAs were identified: MiR-548j host gene HMGB1P10; miR-597 host gene TNKS; miR-1915 host gene CASC10; miR-3146 host gene TWISTNB; miR-5187 host gene TOMM40L; and miR-6814 host gene RIPK4. Whether AnAc selectively inhibits the transcription of these genes via its p300/PCAF histone acetyltransferase (HAT) inhibitory function [169] remains to be examined. Inhibition of HAT activity would be expected to increase gene expression. Interestingly, AnAc inhibits p300/PCAF histone acetyltransferase (HAT) activity [169] and thus could coordinately downregulate this set of miRNAs and host genes by promoting a more condensed genomic state, but experimentally examining the veracity of the supposition is outside this current study and remains to be examined fully. MetaCore transcription factor (TF) network analysis identified CREB1, FosB, SOX4, TCF7L2 (TCF4), PRDM14, JunD, GATA-3, FRA-1, cFos, JunB, FOXp3, and YY1 as significantly associated with these genes. The ability of AnAc to inhibit the activity of these TFs will also need to be experimentally verified.
A decrease in a miRNA would be expected to result in an increase its target transcript expression. Validated targets of each miRNA were identified in the literature. An important note in searching the literature for miRNA targets is that often, whether the miRNA# is the 3p or 5p arm is not stated. However, if the miRNA sequence is provided in a diagram along with the seed match site in a target mRNA's 3'-UTR, the miRNA sequence can be identified as either 3p or 5p by entering the miRNA sequence in miRBase.org. Clearly, a miRNA-3p and miRNA-5p will have different targets, and thus potentially different cellular effects. When identified in our RNA seq study, the 3p or 5p arm is indicated. AnAc reduced miR-378g that targets VIM (vimentin) [51] and VIM mRNA transcript expression was increased in AnAc-treated MCF-7 cells (GSE78011), suggesting a reciprocal regulation. None of the other validated targets of decreased miRNAs (Table 3) were found among the upregulated mRNA transcripts identified in GSE78011. MetaCore network enrichment analysis did not match any of the downregulated miRNAs and Pathway Maps, GO processes, or Process Networks. Networks identified were 1) miR-509: positive regulation of macromolecule metabolic process; 2) miR-584: regulation of gene expression; 3) miR-509, miR584, MDM2, ERK1/2: positive regulation of gene expression (S4 Fig). Based on their CSC and tumor-promoting activities the AnAc downregulation of miR-378g, miR-548, miR-548j, miR-548l (Table 3) would be expected to contribute to the anti-proliferative activity of AnAc.

miRNAs oppositely regulated by AnAc in MCF-7 and MDA-MB-231 cells
In contrast, miR-6873 showed opposite AnAc regulation in the two cell lines: it was downregulated in MCF-7 and upregulated in MDA-MB-231 cells (Tables 2 and 5). There are no publications in PubMed on miR-6873 and miR-6873 was not listed in microRNA.org or miRTarBase. Thus, its relevance to AnAc responses in these two cell lines is unknown.

miRNAs downregulated by AnAc in MDA-MB-231 cells
Twenty-two miRNAs were downregulated by AnAc in MDA-MB-231 cells and none of these overlapped with miRNAs downregulated by AnAc in MCF-7 cells ( Table 5). The chromosome location and host gene, if warranted, of each of the AnAc-downregulated miRNAs are identified in Table 5. To examine if the miRNA host gene was downregulated by AnAc in MDA-MB-231 cells, we searched GSE78011. miR-1277 host gene WDR44 was downregulated by AnAc in MDA-MB-231 cells. WDR44 encodes a protein that interacts with the small GTPase rab11 and is involved in endosome recycling [171]. There are no validated targets for miR-1277 in miRTarBase.
Downregulation of a miRNA would be expected to increase the expression of its targets; hence, we searched our data of mRNAs upregulated by AnAc in MDA-MB-231 cells (550 genes, GSE78011) for the validated targets in Table 5, but none were reciprocally upregulated. This may be because the miRNA and mRNA for RNA seq were extracted at the same time, i.e., after 6 h of AnAc treatment, or that these mRNAs are not expressed or targeted in MDA-MB-231 cells. Given their roles as putative oncomiRs the downregulation miR-23b and miR-1247 may play a role in the anti-proliferative activity of AnAc in in MDA-MB-231 cells.
Analysis of the data identified ZFP36L1 as a putative target of miR-3614 in MDA-MB-231 cells. Interestingly, AnAc downregulated miR-3614 and upregulated ZFP36L1 transcript expression in MDA-MB-231 cells, suggesting an inverse correlation. ZFP36L1 has been identified as a cancer gene due to mutations in breast cancer and acts in a recessive manner [172]. ZFP36L1 is a member of the TTP family of tandem zinc finger proteins that bind AU-rich elements (AURE) in the 3 0 -end of target gene transcripts and promote target degradation, e.g. STARD1 [173], VEGFA [174], NR4A2 [175], BCL2 [176], LDLR [177], STAT5B [178], and CDK6 [179]. Of these genes, only VEGFA and LDLR were identified as differentially expressed genes in AnAc-treated cells.    (Table 6). We have described miR-20b-5p and miR-612 upregulation in the context of similar results in AnAc-treated MCF-7 cells (Table 2, Fig 2, S2 Fig). The chromosome location and host gene, if warranted, of each of the AnAc-upregulated miRNAs are identified in Table 6. Interestingly, most of the downregulated miRNAs were intergenic. miR-1298 is in encoded in HTR2C, but HTR2C was not among the AnAc-regulated genes in MDA-MB-231 cells in GSE78011. An increase in a miRNA would be expected to result in a decrease of its target transcript. miR-20b-5p target EFNB2 (ephrin B2) expression was downregulated in AnAc-treated MDA-MB-231 cells, but none of the validated targets of the upregulated miRNAs (Table 6) were found among the AnAc-downregulated mRNA transcripts identified in RNA seq (GSE78011). Given their roles as 'tumor suppressor' miRNAs (see Table 6), the increases in miR-29b, miR-612, and miR-1298 may contribute to the antiproliferative activity of AnAc in MDA-MB-231 cells.
MetaCore analysis of these upregulated miRNAs identified "cellular response to inorganic substance" as the top GO process (S7A Fig). MetaCore analysis identified two networks: 1) Table 6. miRNAs upregulated by AnAc in MDA-MB-231 cells. The genomic location of each miRNA was identified in miRAD http://bmi.ana.med.uni-muenchen.de/miriad/ [34]. Verified targets are those experimentally validated targets of the indicated miRNA as demonstrated by 3'-UTR luciferase reporter assay. Since many publications do not include whether the 5p or 3p arm of the miRNA was studied, if the sequence of the miRNA was provided, it was searched in miRBase.org to identify which arm was used in the target gene 3'-UTR luciferase reporter assay.

Effect of altered miR-612 on cell viability
Since AnAc increased miR-612 in both MCF-7 and MDA-MB-231 cells (  (Table 4)  qPCR validation of AnAc-mediated changes in mRNAs targeted by miR-378g We selected VIM, a target of miR-378g downregulated by AnAc in MCF-7 cells, and ZFP36L, a target of miR-3614 downregulated by AnAc in MDA-MB-231 cells for validation by qPCR.
As anticipated from the decrease in miR-378g in RNA seq data (Table 3), we detected a slight increase in VIM transcript expression in MCF-7 as well as an increase in VIM in MDA-MB-231 cells (Fig 5). However, because qPCR indicated an increase in miR-378g levels in AnActreated MCF-7 cells (Fig 5), it is possible that VIM is upregulated by AnAc by mechanisms unrelated to miR-378g. In addition, miRNA and mRNA were extracted at the same time, i.e., after 6 h of AnAc treatment, and it may be that changes in VIM mRNA levels require a longer time to be degraded after miR-378g targeting. Transcript levels of ZFP36L were increased in AnAc-treated MDA-MB-231 cells (Fig 5), corresponding with the observed downregulation of miR-3614 (Table 5). These data confirm the reciprocal expression of these mRNA transcripts detected in RNA seq and their target miRNAs in the respective AnAc-treated cell line.

Pathways affected by DEGs and DEmiRs in AnAc-treated MCF-7 cells
MetaCore analysis of DEGs from both mRNA and miRNA data sets of AnAc-treated MCF-7 cells identified NETosis in SLE as the top pathway. The release of neutrophil extracellular traps (NETs) by dying cells (NETosis) was first described as the release of nuclear chromatin, nuclear histones and many granular antimicrobial proteins from neutrophils as one of the first lines of defense against pathogens (reviewed in [180]). The top GO processes were chromatin silencing, negative regulation of gene expression (epigenetic, nucleosome assembly, chromatin assembly, and nucleosome organization. The three gene networks identified were 1): PDEGF PDE6G, APOBEC3H, GGTF II beta, CDIP, p53; 2) miR-499, BMCC1, Histone H1, miR-20b, miR-23b; 3) UCHL1, Protein C, PDK4, EGR1, miR-1298 5p. Network #2 processes include anoikis, negative regulation of fat cell proliferation, regulation of DNA metabolic processes, which reflect the antiproliferative, pro-apoptotic, and NRAM activity of AnAc detected previously in MCF-7 cells [13].

Pathways affected by DEGs and DEmiRs in AnAc-treated MDA-MB-231 cells
MetaCore analysis of DEGs from both mRNA and miRNA data sets of AnAc-treated MDA-MB-231 cells identified "Immune response, IL-3 signaling via JAK/STAT, p38, JNK, and NFkB" as the top pathway. The top GO processes were "Positive regulation of biological process; cellular response to oxygen-containing compound, positive regulation of cellular process, response to oxygen-containing compound, regulation of developmental process, and response to lipid". The three gene networks identified were Network #1: Axin, Frizzled, cMyc, WNT, PI3K reg classIA: canonical Wnt signaling pathway, beta-catenin destruction complex disassembly, regulation of cell proliferation, cell surface receptor signaling pathway involved in cell-cell signaling, cell-cell signaling by wnt. Network #2: C/ EBPbeta, SOS, NGFR, H-Ras, NGF: positive regulation of cellular metabolic process, positive regulation of MAPK cascade, positive regulation of metabolic process, positive regulation of macromolecule metabolic process, and positive regulation of intracellular signal transduction. Network #3: GALNT4, Keratin80, BCMP101, HEXIM1, PNRC1: translational elongation, translation, amide biosynthetic process, peptide biosynthetic process, peptide metabolic process.

Conclusions
In summary, we describe the first comprehensive assessment of miRNA expression in response to anacardic acid in ERα+, luminal A MCF-7 and MDA-MB-231 TNBC breast cancer cells. The pathways modulated by these miRNAs suggest that key nodal molecules, e.g., Cyclin D1, SMAD, SP1, MYC, c-FOS, PPARγ, BCL2, FOXO3A, MDA2, and SIN3, are targets of AnAc activity. In agreement with the pathway analysis, we previously reported that AnAc reduced CCND1 transcript expression in MCF-7 and MDA-MB-231 cells [13]. The roles of the other proteins and pathways in AnAc responses remains to be investigated.