Identification of Differentially Expressed miRNA by Next Generation Sequencing in Locally Advanced Breast Cancer Patients of South Indian Origin

Purpose: miRNAs are known to be aberrantly expressed in the serum, tissue, and Peripheral Blood Mononuclear Cells (PBMC) of cancer patients and could serve as potential noninvasive diagnostic markers for breast cancer. The aim of this study was to identify the differentially expressed miRNA using next-generation sequencing (NGS) from the paired PBMC samples from breast cancer patients and age-matched healthy individuals and explore their functional significance. Methods: In this study, PBMCs were employed for the detection of miRNAs by NGS in locally advanced breast cancer (LABC) women of South Indian origin who were divided into three age groups, (a) 40yrs-50yrs (b) 50yrs-60yrs and (c) 60yrs-70yrs, compared with age-matched control groups. Results: Four miRNAs (hsa-miR-192-5p, hsa-miR-24-2-2p, hsa-miR-3609, and hsa-miR-664b-3p) were found to be differentially expressed among LABC patients compared with age matched healthy women of the South Indian population. While miR-24-2-5p, miR3609, and miR-664b-3p were down-regulated, miR-192-5p was up-regulated. Gene Ontology (GO) annotations implicated miRNA with signaling pathways in peripheral nerve synapses, glutamatergic synapse, and cell morphogenesis, all of which play a pivotal role in the manifestation of cancer. Conclusion: Four miRNAs- 3 (While miR-24-2-5p, miR3609, and miR-664b-3p) downregulated and one upregulated (miR-192-5p) were identified as potential biomarkers for patients with locally advanced breast cancer. These markers could be validated in studies with a larger sample size.


Introduction
Breast cancer (BC) is the most common cancer in women worldwide, affecting 12% of all women and leading to 450,000 deaths each year (Torre et al., 2015) which accounts for 15% of cancer related deaths among women.
Genomic profiling of tumors to predict response to therapies has become an approach of precision medicine (Stover et al., 2015). Identification of biomarkers, recognition of potential targets for therapy and plausible predictions using bioinformatics platforms have accelerated the research. Such multiple approach leads the discovery of microRNAs, a large family of non-coding RNAs that play crucial roles in cellular processes and disease manifestations.
Micro RNAs or miRNAs are non-coding single stranded RNAs of 20-22 nucleotides in length. They exert regulatory functions on the expression of multiple genes by initiating translational inhibition or degradation of their complementary mRNA targets. The varied

Identification of Differentially Expressed miRNA by Next Generation Sequencing in Locally Advanced Breast Cancer
Patients of South Indian Origin characteristics of miRNAs allows for regulation of various cellular processes at different developmental stages (Loh et al., 2019). Various studies indicate that expression of miRNAs are tissue-specific and that are reflected in pathophysiological processes like tumorigenesis and metastasis. miRNAs can be detected not only in tissue but also in blood, serum, urine and other sources that are accessible with minimal invasiveness. Hence identification of miRNAs and the related target mRNAs has become the paramount focus of various clinical research studies, particularly cancer studies (Mallanna et al., 2010).
Expression levels of a few miRNAs were found to be associated with the initiation, development and manifestation of breast cancer, its prognosis and treatment at cellular and subcellular levels (Ohzawa et al., 2017). Dysregulation of miRNA expression and its role in the progression of mammary cancer was first described in 2005 by Iorio et al., (2005) and since then many miRNAs with aberrant expression in conjunction with breast cancer progression have been reported. An et al., (2017) have reported the plausible role of miRNA dysregulation in Having perceived the role of miRNAs and their role in breast cancer development, the race is on for the discovery of novel miRNAs and their target mRNAs so as to identify a novel target or therapeutic means to fight cancer. However, the discovery and profiling of miRNAs using microarrays and RT-qPCR has limitations including low sensitivity to miRNA expression, difficulty deducing structural changes and the inability to identify novel miRNAs. To overcome this limitation massive parallel sequencing technology like NGS would be a better alternative to discover novel mi-RNA which occurs at low frequency. The objective of this pilot study was to identify differentially expressed mi-RNA in the PBMC of patients with locally advanced breast cancer compared with age and gender-matched healthy control by using NGS. In particular, we have focused on the mechanism of profiling miRNA expression associated with breast cancer through examining the expression of their targets, followed by pathway analyses. Hence aim of the present study was to identify the differentially expressed miRNA using next-generation sequencing (NGS) from the paired PBMC samples of locally advanced breast cancer patients and age-matched healthy individuals that have the potential to affect breast cancer tumorigenesis and progression.

Study Population and Characteristics
Fifteen patients who were diagnosed with LABC during the year 2018 were selected from a regional cancer center in South India for the current study. Patients who had not undergone chemotherapy or radiotherapy but had histological evidence of disease with South Indian descent were selected and grouped into three age groups of five patients each-(a) 40yrs-50yrs (b) 50yrs-60yrs and (c) 60yrs-70yrs. Patients diagnosed with triple-negative breast cancer were excluded from the study. Simultaneously, fifteen age and gender matched healthy individuals with no history of chronic disease or acute illnesses were identified as controls for the study. The institutional ethics committee approved the study (JIP/IEC/2017/0267). All the details were explained and written informed consent was taken before including the subject into the study. Blood samples were collected from the study and control groups and stored for further analysis.

RNA isolation and sequencing
For the identification of differentially expressed miRNAs, 5 mL of blood was collected in EDTA tubes from breast cancer patients and their corresponding controls. The samples were then subjected to density gradient centrifugation and peripheral blood mononuclear cells (PBMC) were isolated. RNA was isolated from the PBMCs using miRNA easy kit (Qiagen, Gmbtt Germany) in accordance with the manufacturer's instructions. The yield of RNA was checked by Nanodrop (Bio photometer plus) and Fluorimeter (Qubit4) and the samples were stored at -80°C.The samples were divided into three groups of five each based on the age of patients and controls for expression studies as shown in Table 1.
RNA samples were pooled for each age group and a single sample per group was obtained. This was done as an attempt to eliminate biological variations as a result of differential gene expressions between samples from the same age group.
The six pooled RNA samples were quantified by using the Qubit BR Assay Kit (Invitrogen cat# Q10211). RNA purity check was done by using QiAXPERT and RIN (RNA integrity number) was assessed on Tape station using RNA screen tape (Agilent, Cat# 5067-5576) to detect purity, concentration, and integrity of RNA sample. The RIN number is based on the numbering system from 1 to 10 with 1 being most degraded and 10 the most intact. Next-generation sequencing (NGS) with Illumina (HiSeq2500) technology was used for the miRNA study as explained below.

Statistical Analysis
Differential expression analysis was performed with read count of the known miRNA using DESeq2 to compare miRNA expression between LABC patient and aged matched healthy individual. Log2FC≥ (1) for upregulated and log2FC≤ (-1) for downregulated miRNAs with p-adj ≤0.05 was considered.

Library preparation and Sequencing
The extracted RNA was used to prepare small RNA libraries using the NEB Next Multiplex Small RNA Library Prep kit protocol (NEB, #E7560S) and in accordance with the manufacturer's protocol. Qubit High Sensitivity Assay (Invitrogen, Cat# Q32852) was employed for the quantification of the same. The libraries were then pooled and sequencing was carried out on the Illumina HiSeq 2500.
Bioinformatics analysis: Figure 1 The sequencing of samples resulted in generation of 0.9 to 1.7 GB raw data from each sample. The raw data was subjected to quality check and further analysis.

Pre-processing of data
Raw fastq sequences (PRJNA615345) thus generated were de-multiplexed with CASAVA v1.8pipeline. The raw reads were quality checked with FastQC. FastQC tool analyzes the base quality score distribution, sequence quality score distribution, average base content per read, GC distribution in the reads and read-length distribution. Unwanted sequences such as adapter sequences, primers, poly-A tails and other similar sequences are trimmed using Cutadapt tool (http://www.bioinformatics.babraham. ac.uk/projects/fastqc). Structural RNA contamination including ribosomal RNAs and transfer RNA sequences were removed using Bowtie2 (version2.2.4) with default parameters ( Table 1 in Supplementary data).

Prediction of miRNAs
To identify both known and novel miRNAs, and to determine differential gene expression profiles, the quality controlled reads were aligned to the reference human genome (Human genome -hg19) and to a miRNA database, such as miRBase (http://www.mirbase.org/). The Identification of Deferentially Expressed miRNA in Breast Cancer by Next Generation Sequencing HiSeq2500. More than 15 million reads were obtained for each sample. The reads were subjected to quality control checks using FastQC and trimmed where necessary to retain high quality sequences for further analysis using Cutadapt. In addition, using Bowtie, low-quality sequence reads were discarded (Q < 20). Also, other RNA contamination and sequences <17 bases were quality filtered using the same,1-2% data loss was observed at every step of pre-processing.

miRNA prediction, identification and quantification
Using miRDeep2, the reads were aligned to the reference human genome [Human genome -hg19 from Illumina iGenomes (http://support.illumina.com/ sequencing/sequencing_software/igenome.html)] and to a miRNA database, miRBase (http://www.mirbase.org/). This was done to identify both known and novel miRNAs. The sample wise distribution of known as well as novel miRNAs is summarized in ( Table 2 in supplementary data) while the structures for the same with the highest read count in each of the samples are shown in (Figure  1a and 1b in supplementary data). In the known miRNA top 10 predominant miRNA were found across the 6 sample are shown in (Figure 2 in supplementary data), of which miR-21 identified as the key regulator of oncogenic processes.
DESeq2 was used to compute the differential expression of miRNAs for the control and cancer groups. A total of 640 miRNAs were found to be differentially expressed (Table 3 in supplementary data) in PBMCs of breast cancer patients as compared to that in gender and age-matched healthy women. Out of which only four reads which map to multiple positions within a genome and/or map to known small RNA coordinates (e.g., snRNA, rRNA, tRNA), along with any reads that do not map to the reference genome, were discarded. miRDeep2 is a software package which allows both expression profiling and identification of novel miRNAs. The quality checked reads were mapped to the reference human genome (hg19) from Illumina iGenomes (http:// support.illumina.com/sequencing/sequencing_software/ igenome.html). miRDeep2 is used to identify known as well as novel micro-RNA.

Differential expression analysis
The sRNA sequences were aligned against mature miRNA sequences from miRbase22. Reads were aggregated into counts and the expression values were quantified. Differential expression analysis was done based on normalized deep sequencing counts in RPM (Reads Per million) against known miRNAs using DESeq2. Log2 fold changes with positive (Log2fc≥1) values were taken as up-regulated genes and the same with negative values (Lo2fc≤-1) as down-regulated genes.

Functional analysis of target genes
The potential functions of the differentially expressed miRNAs were identified using the Amigo tools in which the targets of the up and down-regulated miRNAs were used as the input for the gene set enrichment analysis. This include GO (Biological process, cellular component, and molecular function) and Reactome Pathway. In each category, 10 prime values were considered for further analysis. A p-value of less than 0.05 was considered statistically significant.

Results
The present study was aimed to identify potential miRNA signatures related with breast carcinoma in PBMCs on an NGS platform. Expression profiles for the paired samples were generated on the Illumina  Table 2). The numbers of targets for the four differentially expressed miRNAs were predicted using two online web servers-miRDB (https://mirdb.org) and miRWalk.

Functional Annotations of mRNAs (Targets of miRNA)
Functional annotation is a method of recognizing functional elements along the sequence of a given genome. Functional insights into microRNAs are gained from the pathways or biological processes that are overrepresented in their target genes. Moreover, numerous genes could be regulated by a single miRNA and the vice versa is also possible. Functional annotation was performed using the miRNet web server for selected genes that were predicted as miRNA targets. GO annotations in terms of molecular functioning and Reactome pathway were not observed significantly for miRNA-192. Significant up-regulation was observed with biological processes and specifically in the regulation of the localized translations. A similar pattern was observed with cellular components wherein the transferase complexes were highly up-regulated. However, the functional analysis for the down-regulated miRNAs-miR-24-2-5p, miR3609 and miR-664b-3p, demonstrated that the target genes were involved in a variety of positive regulation processes like cell morphogenesis, nervous system development, neuron differentiation. Cellular component analysis identified a significant down regulation of the cell leading edge and cytoplasmic vesicles. Down regulation of glutamatergic synapse was also observed. GO molecular functioning demonstrated association of the down-regulated miRNA with purine ribonucleoside binding, cytoskeletal protein binding, and protein kinase activity while the Reactome pathway analysis demonstrated association with signaling pathways involving Fc-gamma Receptor activation, antimicrobial peptides and gene expression. Further analysis found that their target genes were enriched in pathways involved in many cancers including prostate, colorectal cancer, and breast cancer (Figure 3a, 3b and 3c in supplementary data).
Of the four differentially expressed miRNAs, miR-192 was downregulated in breast cancer patients when compared to healthy individuals. This tumor suppressor property is supported by other studies in different types of cancer where miR192 was found to be downregulated. The study done by Zhang et al reported mir-192 acts as a tumor suppressor by inhibiting cell proliferation and inducing apoptosis in breast cancer (Zhang et al., 2019). In contrast, there is some evidence to suggest that mir-192 act as oncomir by promoting tumor growth in breast cancer patients. (Bhardwaj et al., 2019;Zhou et al., 2018).
The present study documented the downregulation of miR-24-2-5p in healthy individuals compare to LABC patients which indicate its oncomir nature. According to  enhances invasion and progression of tumour by targeting PTPN9 (phosphatases tyrosineprotein phosphatase non-receptor type 9) and PTPRF (receptor-type tyrosine-protein phosphatase F) thereby promoting EGF (Epidermal growth factor) (Du et al., 2013). Its expression, however, was down regulated after chemotherapy (Camps et al,. 2014;Sochor et al., 2014). miR-3609 was found to be downregulated in healthy individual as compared with LABC patient. But in contrast to our study, the other studies reported low miR-3609 expression was involved in the malignancy of breast and other type of cancer Slattery et al., 2014). Report from Fitzpatrick C et al. confirmed, the upregulation of miR3609 sensitized breast cancer to Adriamycin (Fitzpatrick et al., 2019).
According to our finding, miR-664b-3p was downregulated in healthy individual as compared with LABC patients which confirm its oncogenic nature. miR-664b-3p was downregulated after the treatment  (Khraman et al., 2018). In contrast, study by Wu et al confirmed tumour suppressor nature of miR-664 by suppressing the invasion and proliferation in breast cancer by targeting Insulin Receptor Substrate 1 (Wu et al., 2019). Hence further studies are needed to validate the role of miRNA-664b-3p as oncomir or tumour suppressor in breast cancer.
These conflicting evidences about the roles of miRNAs (miR-664b-3p, miR-3609, miR-192) in breast cancer pathogenesis illustrates these miRNAs exert both oncogenic and tumour suppressive effects. (Khraman et al., 2018;Wu et al., 2019) .The overall net oncogenic or tumour suppressive effect of a miRNA depends partly upon the balance between miRNA-mediated upregulation or downregulation of oncogenic and tumour suppressive pathways, as well as tumour immune system interactions and tumour-modifying extrinsic factors. There is a need for further in -vitro or in-vivo experiments, to confirm the relevance of these miRNAs.
Functional annotation from this study implicate miRNA with the manifestation of breast cancer through signaling pathways associated with glutamatergic synapse, cell morphogenesis and neurological synapses. Glutamate regulation has been previously reported in association with cancer progression and its potential as a biomarker for therapeutic purposes has been identified (Stepulak et al., 2014;Prickett and Samules, 2012). Through the peripheral sympathetic nervous system a physiological stress can be exerted which in turn enhances the metastasis of breast cancer (Pimentel et al., 2013). Also, dysregulation of transcription and translation have been time and again reported in association with the metastasis of various cancers (Imai et al., 2005). Based on the functional annotations and Reactome pathway analysis it can be assumed that the differentially expressed miRNA identified herein manifest tumor progression by acting on signaling pathways associated with peripheral synapses and cell morphogenesis.
In conclusion, four differentially expressed mi-RNAs: hsa-miR-192-5p, hsa-miR-24-2-2p, hsa-miR-3609, and hsa-miR-664b-3p were identified in the blood sample of patients with locally advanced breast cancer compared with that of age matched healthy controls. These mi-RNAs were elucidated through massive parallel sequencing and their functions have been mapped. The bioinformatics analysis showed their role in breast cancer development and may act as prognosis or diagnostic biomarker. Further analysis is required to validate these findings in more number of subjects for their plausible role in breast cancer both as diagnostic and prognostic biomarkers.

Author Contribution Statement
Jayanthi Mathaiyan: Conceptualization, Methodology, manuscript reviewing and editing. Shyam Kumar Tripathi: Data curation, manuscript writing and editing. Rajesh NG, Smita Kayal: Supervision, guidance in doing and writing manuscript. Availability of data and materials: Authors consent to the availability of data and materials.