Next-Generation Sequencing Reveals the Role of Epigallocatechin-3-Gallate in Regulating Putative Novel and Known microRNAs Which Target the MAPK Pathway in Non-Small-Cell Lung Cancer A549 Cells

Lung cancer constitutes 85% of non-small cell lung cancer diagnosed cases. MicroRNAs are novel biomarkers that are capable of modulating multiple oncogenic pathways. Epigallocatechin-3-gallate (EGCG) is a potent chemopreventive and chemotherapeutic agent for cancer. We aimed to identify important known and putative novel microRNAs modulated by EGCG in A549 cells using next-generation sequencing and identify their gene targets. Preliminary analysis revealed an IC50 value of 309 μM with G0/G1 phase arrest at 40 μM EGCG treatment. MicroRNA profiling identified 115 known and 4 putative novel microRNAs in 40 μM and 134 known and 3 putative novel microRNAs in 100 μM EGCG-treated A549 cells. The top 10 up-expressed microRNAs were similar between the untreated control and EGCG-treated A549 cells. An up-expression in oncogenic microRNAs, which belong to broadly conserved seed families, were observed in untreated control and EGCG-treated A549 cells. Kyoto Encyclopedia of Genes and Genomes and Protein Analysis Through Evolutionary Relationships pathway analyses of the validated microRNA targeting genes strengthened the hypothesis that EGCG treatment can modulate microRNAs that play a significant role in the MAPK signaling pathway. Expression profile of microRNAs was validation by quantitative real time PCR of randomly selected microRNAs. This study identified signature microRNAs that can be used as novel biomarkers for lung cancer diagnosis.


Introduction
Green tea, brewed from the unfermented dry leaves of the plant Camellia sinensis, is the most consumed non-alcoholic beverage in Asian countries and is now gaining popularity in western countries as well. It contains a wide range of phytochemicals which exhibits anti-cancer, anti-oxidative, and anti-inflammatory properties [1]. Among these phytochemicals, (−)-epigallocatechin-3-gallate (EGCG) accounts for 18-36% of the total phenolic compounds and 70% of the catechins present in green tea. Animal and cell line studies have demonstrated an important role of EGCG in promoting apoptosis and reducing cancer growth. EGCG is a potential chemopreventive and chemotherapeutic compound against skin [2], lung [3], breast [4], colon [5], prostate [6], and other cancers [7,8].
Lung cancer dominates among all the cancers with the leading mortality rate worldwide, 85% of which are non-small cell lung cancer (NSCLCs) cases [9]. Over 50% of the lung cancer patients die within a year of diagnosis with an estimated 1.6 million deaths per year [10]. The prevalent

EGCG Induced G0/G1 Phase Arrest in A549 Cell Line
EGCG treatments significantly increased the percentage of cells at the G0/G1 phase of the cell cycle ( Figure 1). The percentage of cells at the G0/G1 phase of the cell cycle increased from untreated control (55.58%) to 40 μM EGCG treatment (79.15%). Furthermore, we observed that 79. 15, 70.76, 78.76, and 72.04% of cells were persistent at the G0/G1 phase of the cell cycle at 40, 60, 80, and 100 μM EGCG treatments, respectively. A rational cell percentage at the G0/G1 phase of the cell cycle beyond 40 μM EGCG treatment attests to G0/G1 phase arrest at 40 μM EGCG treatment ( Figure 1). In addition, no apoptosis was observed in EGCG-treated A549 cells. EGCG treatment (f) the cell cycle distribution analysis. The data are shown as mean ± SD. Significant (p < 0.05) difference between the groups are indicated by: " * " between untreated control and 40 μM EGCG treatment, "@" between untreated control and 80 μM EGCG treatment, and "$" between untreated control and 100 μM EGCG treatment.

Analysis of MicroRNAs
The miRBase-21 database was used for known miRNA detection using sequence similarity approach (ncbi-blast-2.2.30). Novel microRNA sequences were predicted using Mireap_0.22b [36]. In total, 958, 944 and 935 known microRNAs were detected in the untreated control, 40, and 100 μM EGCG treatments, respectively. MicroRNAs with ≥50 read count was counted as 206, 194, and 199 for untreated control, 40, and 100 μM EGCG treatments, respectively. In addition, on an average, maximum microRNAs were predicted from chromosome 1 followed by chromosome 17 and 14 ( Figure 2). About 105, 108, and 111 microRNAs were predicted from chromosome 1 followed by 82, 77, and 75 from chromosome 17 in the untreated control, 40, and 100 μM EGCG treatments, respectively. EGCG treatment (f) the cell cycle distribution analysis. The data are shown as mean ± SD. Significant (p < 0.05) difference between the groups are indicated by: "*" between untreated control and 40 µM EGCG treatment, "@" between untreated control and 80 µM EGCG treatment, and "$" between untreated control and 100 µM EGCG treatment.

Prediction of Putative Novel MicroRNAs
We identified top ten up-expressed predicted putative novel microRNA sequences in control, 40, and 100 μM EGCG treatments and named as TPP-A549-1, TPP-A549-2 and so on. The chromosomal location, precursor, and mature sequences of the predicted putative novel miRNAs are presented in Table 1. We observed six putative novel microRNA sequences which were expressed in more than one sample. The predicted secondary structure of these putative novel microRNAs is presented in Figure 4. The putative novel microRNA TPP-A549-7 was expressed in all the samples with the MFE (Minimum free energy) value of -29 Kcal/mol (Table 1). In the untreated control, the putative novel microRNA TPP-A549-1 was highly expressed with the read count of 3852. In addition, the putative novel microRNAs TPP-A549-11 and TPP-A549-17 were up-expressed in 40 and 100 μM EGCG treatments with the read count of 998 and 1334 respectively (Table 1).

Prediction of Putative Novel MicroRNAs
We identified top ten up-expressed predicted putative novel microRNA sequences in control, 40, and 100 µM EGCG treatments and named as TPP-A549-1, TPP-A549-2 and so on. The chromosomal location, precursor, and mature sequences of the predicted putative novel miRNAs are presented in Table 1. We observed six putative novel microRNA sequences which were expressed in more than one sample. The predicted secondary structure of these putative novel microRNAs is presented in Figure 4. The putative novel microRNA TPP-A549-7 was expressed in all the samples with the MFE (Minimum free energy) value of -29 Kcal/mol (Table 1). In the untreated control, the putative novel microRNA TPP-A549-1 was highly expressed with the read count of 3852. In addition, the putative novel microRNAs TPP-A549-11 and TPP-A549-17 were up-expressed in 40 and 100 µM EGCG treatments with the read count of 998 and 1334 respectively (Table 1).

Differential Expression Analysis of Known MicroRNAs
A complete microRNA profiling is depicted in Figure 5, indicating the effect of EGCG on A549 cells. MicroRNA expression with greater than 2 log2 fold change was determined in the untreated control vs. 40 µM and the untreated control vs. 100 µM EGCG treatments. The heat maps of the top 100 differentially expressed microRNAs between the samples are presented in Figure 6a,b.
A complete list of up-and down-regulated microRNAs is presented in Table 2. We found 115 microRNAs differentially expressed in the untreated control vs. 40 µM EGCG treatment, and 134 microRNAs differentially expressed in control vs. 100 µM EGCG treatment (Figure 6c). Out of the 115 differentially expressed microRNAs reported in the untreated control vs. 40 µM EGCG treatment, 53 were up-and 62 were down-regulated. Furthermore, in the untreated control vs. 100 µM EGCG treatment, we reported 69 up-and 65 down-regulated microRNAs (Figure 6d).
By comparing the data with all the reported up-regulated microRNAs, hsa-miR-125a-3p showed the highest change of log2 fold expression in the untreated control vs. 40 µM EGCG treatment (7.12 log2 fold change) and untreated control vs. 100 µM EGCG treatment (7.47 log2 fold change). Furthermore, hsa-miR-548o-3p was down-regulated by -9.12 and -8.12 log2 fold change in the untreated control vs. 40 µM EGCG and the untreated control vs. 100 µM EGCG treatments, respectively. We observed 21 up-and 24 down-regulated microRNAs in the untreated control vs. 40 and the untreated control vs 100 µM EGCG treatments ( Figure 7).

Differential Expression Analysis of Putative Novel MicroRNA Sequences
Heat maps were plotted to study the differential expression pattern of putative novel microRNAs. A complete putative microRNA profiling is shown in Figure 8. Significantly differentially expressed putative novel microRNAs are plotted as heat maps as shown in Figure 9a,b. Comparison of differential expression

Differential Expression Analysis of Putative Novel MicroRNA Sequences
Heat maps were plotted to study the differential expression pattern of putative novel microRNAs. A complete putative microRNA profiling is shown in Figure 8. Significantly differentially expressed putative novel microRNAs are plotted as heat maps as shown in Figure 9a,b. Comparison of differential expression of putative novel microRNAs in the untreated control vs. 40, the untreated control vs. 100, and 40 vs. 100 µM EGCG treatments showing greater than 2 log2 fold change revealed 4 putative novel microRNAs differentially expressed in control vs. 40, 3 in control vs. 100, and 4 in 40 vs. 100 µM EGCG treatment. In the untreated control vs. 40 µM EGCG treatment, the putative novel microRNA TPP-A549-26 was up-regulated and three others namely TPP-A549-24, TPP-A549-27, and TPP-A549-28 were down-regulated. In addition, in the untreated control vs. 100 µM EGCG treatment, putative novel microRNAs TPP-A549-32 and TPP-A549-29 were up-regulated and TPP-A549-30 and TPP-A549-31 were down-regulated ( Table 3). The chromosomal location, precursor, and mature sequence details of greater than 2 log2 fold change of putative novel microRNAs in the untreated control vs. 40 and the untreated control vs. 100 µM EGCG treatments are presented in Table 3. Three putative novel microRNAs namely TPP-A549-23, TPP-A549-24, and TPP-A549-25 were commonly differentially expressed in the untreated control vs. 40 and the untreated control vs. 100 µM EGCG treatments (Figure 9). It was interesting to observe that putative novel microRNA TPP-A549-23 was up-regulated in the untreated control vs. 40 µM EGCG treatment and was down-regulated in the untreated control vs. 100 µM EGCG treatment (Figure 9c).

qRT-PCR Analysis of MicroRNAs
The qRT-PCR analysis was performed to validate the NGS dataset. Hsa-miR-21-5p, hsa-miR-548o-5p, hsa-miR-181c, and hsa-miR-212-5p microRNAs were randomly selected for the study. qRT-PCR analysis showed 3.2, and 2.17 log2 fold change in expression in the untreated control vs. 40 μM EGCG treatment and 2.5 and 1.28 log2 fold change in expression in the untreated control vs. 100 μM EGCG treatment in hsa-miR-548o-5p and hsa-miR-181c, respectively. Furthermore, minimal differential log2 fold change expression of 0.31 and 0.03 in the untreated control vs. 40 μM EGCG treatment and 0.76 and 0.6 in the untreated control vs. 100 μM EGCG treatment was observed in hsa-miR-21-5p and hsa-miR-212-5p respectively ( Figure 10). In microRNA sequencing data, an up-regulation of hsa-miR-181c by 3.51 and 1.5 log2 fold was observed in the untreated control vs. 40 and the untreated control vs. 100 μM EGCG treatments, respectively. Furthermore, a significant down-regulation of hsa-miR-548o-5p by 9 and 8 log2 fold change was noted in the untreated control vs. 40 and the untreated control vs. 100 μM EGCG treatments respectively ( Figure 11). About 0.1 and 0.5 log2 fold change of hsa-miR-212-3p was observed in untreated control vs. 40 μM EGCG treatment. In hsa-miR-181c, the log2 fold change of 0.15 and 1.5 was noted in untreated control vs. 100 μM EGCG treatment. A comparative fold change between the qRT-PCR and sequencing dataset supports each other. Therefore, q-RT PCR results validate the present NGS dataset. In microRNA sequencing data, an up-regulation of hsa-miR-181c by 3.51 and 1.5 log2 fold was observed in the untreated control vs. 40 and the untreated control vs. 100 µM EGCG treatments, respectively. Furthermore, a significant down-regulation of hsa-miR-548o-5p by 9 and 8 log2 fold change was noted in the untreated control vs. 40 and the untreated control vs. 100 µM EGCG treatments respectively (Figure 11). About 0.1 and 0.5 log2 fold change of hsa-miR-212-3p was observed in untreated control vs. 40 µM EGCG treatment. In hsa-miR-181c, the log2 fold change of 0.15 and 1.5 was noted in untreated control vs. 100 µM EGCG treatment. A comparative fold change between the qRT-PCR and sequencing dataset supports each other. Therefore, q-RT PCR results validate the present NGS dataset.  Figure 12). This fold change was compared with the log2 fold change obtained from the NGS data analysis. This sequence was up-regulated by 1.45 log2 fold change in the untreated control vs. 40 μM EGCG treatment and was down-regulated by 4.21-fold change in the untreated control vs. 100 μM EGCG treatment ( Figure 12).
It was observed that log2 fold change obtained in sequencing dataset was similar to the log2 fold change obtained from the qRT-PCR analysis. Furthermore, the present study also validates a novel microRNA sequence obtained in the NGS dataset. This qRT-PCR analysis attests to our computational analysis of NGS data.  (Figure 12). This fold change was compared with the log2 fold change obtained from the NGS data analysis. This sequence was up-regulated by 1.45 log2 fold change in the untreated control vs. 40 µM EGCG treatment and was down-regulated by 4.21-fold change in the untreated control vs. 100 µM EGCG treatment ( Figure 12).
It was observed that log2 fold change obtained in sequencing dataset was similar to the log2 fold change obtained from the qRT-PCR analysis. Furthermore, the present study also validates a novel microRNA sequence obtained in the NGS dataset. This qRT-PCR analysis attests to our computational analysis of NGS data.

KEGG and PANTHER Pathway Enrichment of Targets of Validated microRNAs
The high precision target prediction for hsa-miR-548o-5p, hsa-miR-181c, hsa-miR-212-5p, and hsa-miR-21-5p was carried out using TargetScan and miRDB target computational prediction software. Default cut off values were used for gene target prediction. KEGG and PANTHER pathway analysis were carried out using the Database for Annotation, Visualization and Integrated Discovery (DAVID) and the pathways were shortlisted. Pathway analysis for microRNAs was evaluated and common pathways predicted between TargetScan and miRDB are presented in Figure 13. Pathways in cancer, MAPK, regulation of actin cytoskeleton, wnt signaling, ErbB signaling, B-cell and T-cell receptor signaling, and long-term potentiation were the most significant pathways obtained in KEGG pathway analysis (Figure 13 a,c,e, and g). In addition, Ras signaling, angiogenesis, FGF signaling, wnt signaling, FGF, and PDGF (Platelet-derived growth factor) signaling pathway genes were reported in PANTHER pathway analysis (Figure 13b,d,f,h).

KEGG and PANTHER Pathway Enrichment of Targets of Validated microRNAs
The high precision target prediction for hsa-miR-548o-5p, hsa-miR-181c, hsa-miR-212-5p, and hsa-miR-21-5p was carried out using TargetScan and miRDB target computational prediction software. Default cut off values were used for gene target prediction. KEGG and PANTHER pathway analysis were carried out using the Database for Annotation, Visualization and Integrated Discovery (DAVID) and the pathways were shortlisted. Pathway analysis for microRNAs was evaluated and common pathways predicted between TargetScan and miRDB are presented in Figure 13. Pathways in cancer, MAPK, regulation of actin cytoskeleton, wnt signaling, ErbB signaling, B-cell and T-cell receptor signaling, and long-term potentiation were the most significant pathways obtained in KEGG pathway analysis (Figure 13a,c,e,g). In addition, Ras signaling, angiogenesis, FGF signaling, wnt signaling, FGF, and PDGF (Platelet-derived growth factor) signaling pathway genes were reported in PANTHER pathway analysis (Figure 13b,d,f,h). Furthermore, KEGG pathway analysis showed the MAPK signaling pathway as the target for hsa-miR-21-5p, hsa-miR-548o-5p, hsa-miR-181c, and hsa-miR-212-5p microRNAs (Table 4). Our analysis with PANTHER pathway did not show any common pathway among the microRNAs. Ras signaling pathway, wnt pathway, angiogenesis, p53, and EGF (Epidermal growth factor) receptor signaling pathway were the most significant pathways predicted by TargetScan and miRDB target list.
MiRanda software was used for target prediction and pathway analysis of putative novel microRNA sequences. The common pathways found in the target prediction were wnt, angiogenesis, p-53, PI3K, and MAPK signaling pathways. The putative novel microRNA TPP-A549-23 targeted  19) were noted to be potential target genes for the putative novel microRNA TPP-A549-25. The pathway analysis showed that the putative microRNAs could play an important role in cell cycle proliferation, MAPK, Hedgehog, FOXO, and TGF-beta signaling cascade. We believe that these predicted putative novel microRNA sequences play a major role in cancer proliferation and metastasis.

Discussion
To evaluate the changes of EGCG induced microRNA expression in A549 cells, an established in vitro model of human lung adenocarcinoma, and next-generation sequencing analysis were employed. The dose of EGCG used in the present study was decided based on cell cycle analysis and previous literature. EGCG is capable of causing G0/G1 phase arrest in many cancer cell lines including A549 [37][38][39][40]. The present cell cycle data analysis showed G0/G1 phase arrest at 40 µM EGCG treatment. Kweon et al. [41] reported that out of eight cell lines they studied, A549 showed no sign of apoptosis and was highly resistant even at 100 µM EGCG treatment. As much as 85% of the cell viability was sustained for 72 h at 40 µM EGCG treatment. We observed browning of the medium above 100 µM EGCG treatment. Previous data suggested that O 2 − and quinones generation in the cell culture medium occurs due to the auto-oxidation properties of EGCG [42,43]. Hence, in the present study, we have chosen 40 µM and 100 µM EGCG concentrations. Nine hundred and fifty-nine microRNAs out of 1881 microRNAs (50.9%) reported in miRbase (as per Feb 2018) were expressed in A549 cells. Some microRNAs did not show differential expression while certain other microRNAs were significantly influenced by EGCG treatment. Therefore, EGCG dependent microRNA profiling was studied according to the three independent criteria namely log2 expression analysis, differentially up-and down-regulated microRNAs and putative gene targets of the microRNAs.
Interestingly, the top 10 up-expressed microRNAs namely, hsa-miR-21-5p, hsa-let-7i-5p, hsa-miR-100-5p, hsa-miR-27b-3p, hsa-miR-151a-3p, hsa-miR-148a-3p, hsa-miR-30a-5p, hsa-miR-192-5p, hsa-miR-3529-3p, and hsa-miR-30d-5p were similar in untreated control and EGCG treatments. The hsa-miR-21 is the most commonly up-expressed microRNA in human cancers [44][45][46]. This microRNA is the first one to be named as "oncomir" [47] and was found to highly up-expressed in 540 clinical samples from cancer patients [48]. Yang et al. [49] reported significant up-regulation of hsa-miR-21 in NSCLC patients. In addition, they also reported that inhibition of hsa-miR-21 expression reduced cell proliferation, migration, and invasion in A549 cell. Up-expression of hsa-miR-21 in untreated control A549 cells further attests to our NGS data. Furthermore, due to treatment with EGCG, no significant down-regulation of hsa-miR-21 was observed in differential expression analysis of our NGS data (Figure 11a) which is validated by qRT-PCR analysis (Figure 10a). It is well known that hsa-miR-21 is an oncogene which plays an important role in programmed cell death and targeting apoptosis [50] but, in the present study no such apoptosis was observed on EGCG treatments. Therefore, we hypothesize that EGCG treatment has no modulatory effect on hsa-miR-21. Therefore, the present data support the non-apoptotic effect of EGCG in A549 cells. Another clinical study validated the chemoresistance role of hsa-mir-21 on platinum-based chemotherapy [51]. The persistent expression of hsa-miR-21 in untreated control and EGCG-treated A549 cells in the present study is supported by the clinical microarray data validated by Fujita et al. [51].
Let-7i, a member of the let-7 family, is an oncogenic driver of NSCLC. Increased expression of let-7i in the present study supports the statistically significant clinical study conducted on NSCLC [52]. Furthermore, a microarray analysis of primary lung cancer tumors and non-cancerous lung tissues revealed down-regulation of hsa-miR-181c [53]. In addition, Fujita et al. [51] reported the significance of hsa-miR-181c. It was observed that patients which did not respond to platinum-based chemotherapy had elevated hsa-miR-181c expression as compared to the well responsive patients. This study thus validated the chemoresistant role of hsa-miR-181c [51]. In the present study, EGCG treatment elevated the expression of hsa-miR-181c (Figure 9c). Moreover, this data was validated by qRT-PCR analysis (Figure 8c), where we found a significant up-regulation of hsa-miR-181c by 2.17 and 1.28 log2 fold change by 40 and 100 µM EGCG treatments, respectively.
The direction of log2 fold change (up-or down-regulation) was analyzed in microRNA expression. About 53 and 69 microRNAs were up-regulated after 40 and 100 µM EGCG treatments, respectively. Hsa-miR-125a-3p was highly up-regulated with a log2 fold change of 7.12 and 7.47 in 40 and 100 µM EGCG treatments, respectively. Jiang et al. [54] showed an inverse relationship between hsa-miR-125a-3p with invasion and metastasis. Furthermore, it was marked that hsa-miR-125a-3p was down-regulated in NSCLC [54]. Significant up-regulation of hsa-miR-125a-3p was observed with EGCG treatment in A549 cells. As it was established that EGCG inhibits cell proliferation in A549 cells by causing G0/G1 phase arrest (Figure 1), up-regulation of signature microRNA hsa-miR-125a-3p attests to the role of EGCG in inhibiting cell proliferation and invasion.
Hsa-miR-548o-3p was significantly down-regulated by 9.12 and 8.12 log2 fold change in 40 and 100 µM EGCG treatments, respectively. Mir-548 family consist of 68 microRNAs and is the largest, and poorly conserved primate-specific gene family [55]. The reports on one of the members of mir-548 family, hsa-miR-548c-3p, revealed it as a functional biomarker in prostate and gastric cancer progression [54,56]. However, the exact function of hsa-miR-548o-3p in A549 cells is still unknown, but evidence supports the up-regulation of the mir-548 family in cancer progression. We noted down-regulation of hsa-miR-548o-3p by EGCG treatments (Figures 8b and 9b) which indicated it as a potential biomarker for cancer progression.
Hsa-miR-212 is a tumor suppressor microRNA which negatively regulates anti-apoptotic protein PED/PEA-15 [57]. A study by Jiang et al. [58] validated the biological role of hsa-miR-212/132 in A549 cells. Up-regulation of hsa-mir-212 blocked proliferation and migration and was observed to cause cell cycle arrest by modulating p21 and cyclin D1 expression [58]. An analysis of our NGS data supported the biological role of hsa-mir-212. We observed no significant expression of hsa-mir-212 in untreated control A549 cells. The lower expression of hsa-miR-212 signifies the cancerous property of A549 cells. In addition, negligible differential expression was observed in the EGCG-treated A549 cells. This indicated the fact that A549 cells are resistant to EGCG treatment.
Notably, a significant modulation in hsa-miR-146b-3p was observed with EGCG treatment. Down-regulation of hsa-miR-146b-3p was observed with 40 µM EGCG treatment and a notable up-regulation was noticed in 100 µM EGCG treatment ( Figure 6). A clinical study elucidated a down-regulation of hsa-miR-146b-3p in breast cancer tissues. Further findings revealed that hsa-miR-146b-3p's over-expression suppressed migration, invasion, metastasis, and growth in breast cancer cell lines [59]. The principal findings on the role of hsa-miR-146-3p in untreated control A549 cells revealed the association of its lower expression with early-stage NSCLC recurrence [60]. Functional role of hsa-miR-146b-3p is still unexplored, but our data indicated its important role in tumor suppression.
The goal of this study was to identify putative novel microRNA sequences in A549 cells and their potential differences in expression between untreated control and EGCG treatments. We reported, 4, 3, and 4 putative novel microRNA sequences in untreated control vs. 40, untreated control vs. 100, and 40 vs. 100 µM EGCG treatments respectively (Table 3). Only three putative microRNA sequences were persistently present between the untreated control vs. 40 and the untreated control vs. 100 µM EGCG treatments (Figure 9c). To validate the log2 fold change of the putative microRNA sequences, qRT-PCR analysis was done. The expression analysis of the putative novel microRNA TPP-A549-23 was done using qRT-PCR analysis. We observed similar log2 fold change between the untreated control vs. 40 and the untreated control vs. 100 µM EGCG treatments (Figure 8). The NGS dataset and validation confirmed the significant expression of this putative novel microRNA. Further validation of the sequence is required.
TargetScan and miRDB databases identified approximately 1200 potential target genes of microRNAs (data not shown). Next, we have used DAVID web-based gene ontology and pathway prediction software for further analysis. KEGG pathway and PANTHER pathway analysis were carried out with both the gene targets obtained by TargetScan and miRDB.
In summary, we demonstrated the important role of EGCG in modulating microRNA expression profiling. Substantial evidence from previous studies have provided an important role of EGCG in inhibiting cell proliferation in lung cancer, but its mechanistic insights on microRNA profiling remain incompletely understood. However, our results provided a useful approach for better understanding of EGCG induced microRNA modulation. Furthermore, exploration with transfection and human lung cancer tissue should be performed to validate the microRNA profiling and their predicted targets.

Cell Culture and Treatment
Human non-small cell lung cancer A549 cells (purchased from National Centre for Cell Science, Pune, India) was cultured in F12 Ham-K medium containing 1X antibiotic-antimycotic solution and 10% FBS (Himedia, India) at 37 • C with 95% humidified atmosphere and 5% CO 2 . For cell cycle analysis, A549 cells were seeded in 60 mm cell culture dishes for 24 h in the F12-Ham-K medium. The cells were further treated with 20, 40, 80, 100, and 150 µM EGCG (Sigma Chemical Co. St. Louis, Missouri, United States, USA). After 24 h, cells were harvested and re-suspended in 100 µL PBS. The cells were fixed by drop-wise addition of 8-10 mL 75% ice-cold ethanol, continuous stirring, and incubating on ice for 15 min. After centrifugation, the cells were re-suspended in 200 µL PBS containing 0.04 mg/mL propidium iodide and 0.1 mg/mL RNase and was incubated for 30 min at 37 • C in dark. Cell cycle analysis was carried out on BD FACSVerse™ flow cytometer (ThermoFisher Scientific, Massachusetts, United States.

Next-Generation Sequencing
A549 cells were treated with EGCG at the concentration of 40 and 100 µM. Untreated cells were used as a control. Cells were harvested for RNA extraction after 24 h of culture. Total RNA was extracted from the untreated and treated A549 cells using TRizol reagent (Takara Bio Inc., Kusatsu, Japan) following the manufacturer's protocol. The small RNA library construction and deep sequencing were carried out at AgriGenome Labs Pvt. Ltd., Kochi, Kerala, India. NEBNext®Multiplex Small RNA Library Prep Set for Illumina was used for library construction. The raw counts and the normalized files were submitted in NCBI (National Center for Biotechnology Information) by accession number GSE110514.

Classification and Differential Expression Analysis of MicroRNAs
The human genome (GRCh38) was used for the mapping of microRNA reads using Bowtie. Known miRNAs were identified using miRBase-21 database based on miRNA sequence similarity approach. The sequences were checked for other ncRNA contamination (rRNA, tRNA, snRNA, snoRNA, and piRNA). The novel microRNA prediction was evaluated by Mireap_0.22b [36]. Furthermore, secondary hairpin structures were predicted using mfold [61]. Expressed reads for each microRNAs were calculated and DESseq R software package was used for differential expression analysis. Differentially expressed microRNAs in control vs. 40, control vs. 100 and 40 vs. 100 µM EGCG treatments were determined by their expression in each sample [62]. The expressed reads in untreated control, 40, and 100 µM EGCG treatments were used to calculate the log2 fold change of expression between untreated control vs. 40, untreated control vs. 100, and 40 vs. 100 µM EGCG treatments.

Validation of MicroRNAs
To validate the expression of some of the significant known and putative novel microRNAs, total microRNA was isolated using miRNeasy Mini Kit (Qiagen Sciences, USA) from untreated control and EGCG-treated (40 and 100 µM) A549 cells. First-strand cDNA was synthesized using miScript PCR starter kit (Qiagen Sciences, Germantown, MD, USA) and SYBR green (nucleic acid strain) was used for qRT-PCR (ABI Prism 7000, ThermoFisher Scientific, Massachusetts, United States). The log2 fold change was calculated using the ∆∆CT method [63].
The comparative analysis of qRT-PCR and NGS was done to validate the microRNA expression profile obtained by NGS. The log2 fold change obtained by qRT-PCR of the known microRNAs namely miR-548o-3p, miR-212, miR-125a, and miR-181c was calculated and compared with the log2 fold change obtained in the NGS sequencing data. The log2 fold change of the putative novel microRNA TPP-A549-23 obtained by qRT-PCR and NGS was also compared for the sequence validation.

KEGG and PANTHER Pathway Enrichment of Targets of Validated MicroRNAs
The target prediction for known microRNAs was performed using TargetScan and miRDB [56] and the pathway analysis was done by The Database for Annotation, Visualization and Integrated Discovery (DAVID V 6.7). In addition, target prediction of the putative novel microRNA sequences was carried out using miRanda software [64].

Statistical Analysis
Statistical analysis was performed with one-way analysis of variance (ANOVA). The experimental data was represented as mean ± SD. The results were considered significant when p < 0.001, p < 0.01 or p < 0.05.

Conclusions
We used NGS to study a complete microRNA profiling of A549 cells and identified 958, 944 and 935 known microRNAs in the untreated control, 40, and 100 µM EGCG treatments, respectively. The oncogenic microRNAs were highly up-expressed in the untreated control and EGCG-treated A549 cells which are the part of broadly conserved let-7, hsa-miR-21, and hsa-miR-30 seed families. The up-expressions of these oncogene microRNAs with EGCG treatment indicated resistance character of A549 cells for EGCG. The differential expression analysis of microRNAs identified highly up-regulated hsa-miR-125a-3p and highly down-regulated hsa-miR-548o-3p in the untreated control vs. 40 and the untreated control vs. 100 µM EGCG treatments, respectively. A similar log2 fold change in the comparative analysis between NGS and qRT-PCR of randomly selected known and putative novel microRNAs validated the NGS data. Our data indicated EGCG as an effective natural compound which regulates microRNA profile in A549 cells. This study also attested the modulation of microRNAs by EGCG which regulates cell cycle and inhibits cell proliferation and metastasis. KEGG and PANTHER pathway analysis revealed the MAPK pathway as the most potent targeted pathway by EGCG modulated microRNAs. The findings explored an important role of EGCG in microRNA regulation, which targets MAPK cascade. Furthermore, the putative novel microRNA sequences reported in this study can be a novel approach towards the microRNA targeting gene therapies.