Systematic integrated analyses of methylomic and transcriptomic impacts of early combined botanicals on estrogen receptor-negative mammary cancer

Dietary botanicals such as the cruciferous vegetable broccoli sprouts (BSp) as well as green tea polyphenols (GTPs) have shown exciting potential in preventing or delaying breast cancer (BC). However, little is known about their impact on epigenomic aberrations that are centrally involved in the initiation and progression of estrogen receptor-negative [ER(−)] BC. We have investigated the efficacy of combined BSp and GTPs diets on mammary tumor inhibition in transgenic Her2/neu mice that were administered the diets from prepubescence until adulthood. Herein, we present an integrated DNA methylome and transcriptome analyses for defining the early-life epigenetic impacts of combined BSp and GTPs on mammary tumors and our results indicate that a combinatorial administration of BSp and GTPs have a stronger impact at both transcriptome and methylome levels in comparison to BSp or GTPs administered alone. We also demonstrated a streamlined approach by performing an extensive preprocessing, quality assessment and downstream analyses on the genomic dataset. Our identification of differentially methylated regions in response to dietary botanicals administered during early-life will allow us to identify key genes and facilitate implementation of the subsequent downstream functional analyses on a genomic scale and various epigenetic modifications that are crucial in preventing ER(−) mammary cancer. Furthermore, our realtime PCR results were also found to be consistent with our genome-wide analysis results. These results could be exploited as a comprehensive resource for understanding understudied genes and their associated epigenetic modifications in response to these dietary botanicals.

www.nature.com/scientificreports/ epidermal growth factor 2 (HER2) gene. The treatment of TNBC is more challenging than ER(+) BC and the tumor progresses at a faster rate in comparison to ER(+) BC 4 . The absence of the ER gene renders it even more challenging to control, thereby confining the treatment options. Consumption of dietary botanicals such as green tea and cruciferous vegetables demonstrates promising results in cancer prevention. Epidemiological studies have shown that higher consumption of various fruits and vegetables in countries like Asia, leads to a decline in BC incidence 5 . Another study in Mediterranean and Greek populations has provided strong evidence supporting the hypothesis that a diet high in fiber may reduce the risk of ER(−) BC 6 . This dietary pattern contains multiple protective substances such as high amounts of fiber, essential fatty acids and vitamins E and C, which may be associated with lowering the risk of BC 7 . Mechanistically, dietary polyphenols exhibit an active involvement in various cancer pathways such as the RTK/RAS, PI3K, p53 and cell signaling pathways, which may contribute to their chemopreventive effects on cancers 8 . Amongst various dietary botanicals, isothiocyanates in cruciferous vegetables such as cabbage, kale, cauliflower and broccoli sprouts (BSp) as well as green tea polyphenols (GTPs) play an active role in preventing cancer 9 . Our previous studies have suggested that the combined administration of epigallocatechin-3-gallate in GTPs and the isothiocyanate, sulforaphane in BSp led to synergistic inhibition of cellular proliferation 10 . Furthermore, this combinatorial dietary treatment also resulted in inhibition of tumor development in a BC xenograft mice model 10 . However, even with the recent advancement towards analyzing various anti-cancerous properties of these dietary phytochemicals, little is known about their impact on the epigenomic machinery in ER(−) BC.
Herein, we present an integrated analyses of reduced-representation bisulfite sequencing (RRBS) methylome analyses coupled with RNA-seq transcriptome analyses to evaluate whether early life consumption of combined BSp and GTPs is highly effective in neutralizing epigenomic aberrations and epigenetic control of key genes which can lead to inhibition of ER(−) BC formation and progression. Using a Her2/neu transgenic mammary cancer ER(−) mouse model, our results not only exhibited the combinatorial impacts of BSp and GTPs-induced global DNA methylomic changes, but also illuminated the impacts of genomic regulatory RNA changes on gene expression at specific life stages. Our DNA methylome analyses based on p value ≤ 0.05, identified 2874 (2252 hypermethylated vs. 622 hypomethylated) differentially methylated regions (DMGs) in tumors from mice treated with BSp, 4074 (3538 hypermethylated vs. 536 hypomethylated) DMGs from GTPs treatment and 4181 (3639 hypermethylated vs. 542 hypomethylated) DMGs in mice that received the combination (BSp + GTPs) treatment. Similarly, transcriptomic analyses based on p value ≤ 0.05, identified a total of 146 differentially expressed genes (DEGs) in the tumors from mice receiving BSp, no DEGs in the tumors from mice receiving GTPs and 895 DEGs in the tumors from mice receiving the combination treatment group.
Subsequently, 8 DMGs-DEGs correlated pairs were identified in the mammary tumors of mice receiving BSp treatment and 39 DMGs-DEGs correlated pairs were identified in the tumors of mice receiving the combinatorial (BSp + GTPs) diet. Our pathway analyses further revealed association of DEGs and DMGs with various epigenetics modifications such as DNA methylation and histone modifications (such as Histone H3K4 methylation, Histone H3-K4 trimethylation and histone acetylation). These results indicate that a combinatorial administration of BSp and GTPs has a stronger impact at both the transcriptome and methylome levels in comparison to BSp or GTPs alone. Thus, our dataset including both DNA methylomic and RNA-seq transcriptome analyses is very exclusive and could be of importance to future investigation by revealing the impact of these dietary combinations as well as other novel dietary combinations.

Results
The combinatorial treatment demonstrated higher efficacy in preventing ER(−) mammary tumor development in comparison to BSp and GTPs administered singly. To observe tumor development in response to BSp, GTPs and combination (BSp + GTPs) treatments, we evaluated the effects of these treatment approaches in Her2/neu transgenic mouse model that develops spontaneous ER(−) mammary tumors driven by overexpression of an oncogene. We found that the combinatorial treatment group exhibited the strongest inhibitory effect on tumor growth (Fig. 1). The BSp and the combination treatment groups displayed a significant decline in tumor incidence from 21 weeks of age and thereafter (p value < 0.05) (Fig. 1a). Further, the weight of tumors was significantly decreased in the BSp, GTPs and combination (BSp + GTPs) treatment groups with p value < 0.001 (Fig. 1b). Early life BSp and/or GTPs administration had no effect on mouse body weight, food and water consumption of Her2/neu mice at the same point in the time ( Supplementary  Fig. S1). Subsequently, the BSp treatment group was more efficacious than the GTPs treatment group. Overall, the combination treatment group was most effective in suppressing the tumor development, as demonstrated by its lowest tumor incidence and smallest tumor weight in comparison to the BSp and GTPs treatment groups.

Differential gene expression profile analyses induced by combined administration of BSp and
GTPs. Our current studies and previous publications indicate that combinatorial treatment with BSp and GTPs exhibited the most significant preventive and inhibitory effects on BC in vitro compared to these two compounds that were administered singly. In order to better understand the global influence of combination dietary treatment compared to the individually administered BSp or GTPs, we evaluated the RNA expression using RNA sequencing (RNA-seq) analyses in mammary tumors of Her2/neu mice in the different treatment groups (Control-BSp, Control-GTPs, Control-Combination) as done previously 11 . The mammary tumor tissues were harvested at 31 weeks (endpoint) (N Control = 5, N BSp = 5, N GTPs = 5, N Combination = 5) (Fig. 2) and further sent for RNA-seq pair-end library preparation. Initially, the RNA-seq data were transformed for linear modeling and samples outliers were identified by generating a boxplot across all the samples in different treatment groups ( Supplementary Fig. S2a). As a result, GTP1 displayed abnormal distribution and was removed from further processing/analyses ( Supplementary Fig. S2b) Fig. S2c). Furthermore, unsupervised principal component analyses (PCoA) were conducted on gene expression profiles for individual samples among different treatment groups to reveal gene expression shift (Fig. 3a). Based on the spatial arrangements in a PCoA plot, GTPs had no significant overlap; however, we noticed some overlaps among BSp and the combination treatment groups (Fig. 3a,b). The latter result was verified by differential expression analyses using qRT-PCR. Based on the spatial arrangements of GTPs, as expected, there were no transcripts which were differentially expressed (DE). However, gene expression was slightly changed among the BSp treatment group compared to the control group. Among the total 14,157 transcripts (genes) detected, 146 (1.03%) genes were DE in BSp treatment with 90 (61.64%) genes up-regulated and 55 (38.35%) downregulated using a 5% false discovery rate (FDR) and fold-change (log 10 FC) cutoff. In comparison to BSp and GTPs treatment alone, gene expression profiling was drastically changed among the combination treatment group. Out of 14,157 genes in the combination treatment, 895 (6.32%) genes were DE amongst which 575 (64.25%) genes were up-regulated and 320 (35.75%) genes were down-regulated using a 5% FDR and fold-change cutoff (Fig. 3b). The list of all the transcripts that are DE across BSp, GTPs and combination (BSp + GTPs) treatment groups is displayed in Supplementary File 1: Table S1-S5. The top 30 down-regulated and up-regulated genes ranked by statistical significance with the combination treatment are shown in Tables 1 and 2, respectively. Additionally, out of 146 DEGs in the BSp treatment group and 895 genes in the combination (BSp + GTPs) diet group, we identified 125 overlapping transcripts that were DE. Out of 125 overlapping DEGs, 75 genes were up-regulated, and 50 genes were down-regulated.
To further elucidate the transcriptional profiles changes across different dietary treatments, we generated a heatmap with the top 30 up-regulated and down-regulated genes in the combination group between rows corresponding to DE genes and columns corresponding to biological replicates in the control and combination treatment groups (Fig. 3c). We also generated a heatmap of DE genes that were up-regulated and downregulated in the BSp treatment group (Supplementary Fig. S3). In the combination treatment group, Pkd1 and Bdp1 genes were the top two up-regulated genes. These genes are found to be down-regulated in many different types of cancers. Pkd1 is a tumor suppressor gene that plays a crucial role in many different types of cancers including BC. Bdp1 is a protein coding gene found in TFIIIB that can regulate JNK1 expression in c-jun N-terminal kinases (JNKs), which displays both oncogenic and tumor suppressor properties 12,13 . Furthermore, the most down-regulated www.nature.com/scientificreports/ genes were related to biological processes such as DNA-templated transcription positive regulation of "DNA binding" "protein phosphorylation", "intracellular signal transduction" and top up-regulated genes were related to "histone acetylation", "covalent chromatin modifications", "apoptosis" and "DNA-methylation".
To gain a holistic understanding of DEGs with combination and BSp treatments, we used WebGestalt software to perform GO SLIM subset analyses which can identify major clusters within various biological processes, cellular components and molecular functions. In the combination treatment group, among the biological processes category, the major subsets were "metabolic processes" and "biological response to stimulus" (Fig. 4a). Similarly, in the cellular component category, the most abundant terms were "membrane' , "nucleus" and "cytosol" (Fig. 4b) and in the molecular function category and the most frequent terms were related to "protein binding", "ion binding" and "nucleotide binding" (Fig. 4c). Due to lesser transcripts that were DE in the BSp treatment group, we noticed fewer abundance at these subcategories level ( Supplementary Fig. S4a-c).
Additionally, we used DAVID to perform a broader analysis of specific biological function related to DEGs in BSp and combination treatment groups. In the BSp treatment group, 33 GO terms (p value ≤ 0.05) mainly related to regulation of transcription" (GO:0006355 and GO:0006351), "oxidative-reduction process" (GO: 0055114) and "DNA methylation" (GO:0006306) ( Supplementary Fig. S4d). Unlike the BSp treatment group, combination treatment group exhibited higher gene ontology functions with 50 GO terms which were statistically significant www.nature.com/scientificreports/ DE transcripts. These transcripts were primarily related to various epigenetics mechanisms and other biological pathways such as "cellular response to DNA damage stimulus" (GO: 0006974), "covalent chromatin modifications" (GO:0016569), "DNA repair" (GO:0006281), "mitotic nuclear division" (GO: 0007067), "methylation" (GO:00032259), "histone H3-K4 methylation" (GO:0051568) and "histone H3-K4 trimethylation" (GO:0080182) (Fig. 4d). Although the combination treatment group exhibited more associations with gene ontology functions, there were many biological processes overlapping between the treatment groups. Based on these gene expression results, we concluded that the combination (BSp + GTPs) treatment can impose greater effects on expression level changes by identifying a greater number of DE genes (DEG Combination = 895) in comparison to BSp and GTPs administered alone with less DE genes (DEG BSp = 145 and DEG GTPs = 0). Further, these transcripts in the combination treatment group elicited greater changes at the functional level by identifying various epigenetics modifications which might eventually result in delaying or inhibiting ER(−) mammary cancer. Consistent with the phenotypic changes that the combination treatment can lead to the most promising inhibitory effects on mammary cancer development as compared to single compound treatment alone, these dramatic transcriptomic pattern changes may contribute to chemopreventive effects induced by combination treatment against mammary cancer.

Combination of dietary botanicals induces genome-wide differential DNA methylation patterns.
To further assess the changes in DNA methylation patterns on combinatorial dietary treatment in comparison to BSp and GTPs treatment alone, we applied RRBS analyses across harvested mammary tumor samples. A total of twenty libraries were constructed, and each of them produced a minimum of 5 Gb clean www.nature.com/scientificreports/ reads, which were sequenced and aligned to the reference genome of Mus musculus (mm10) using Bismark (https:// www. bioin forma tics. babra ham. ac. uk/ proje cts/ bisma rk/). The reads of individual samples mapped to the reference genome, thereby producing the relevant BAM files for different samples within each group. Furthermore, these files were used for final downstream analyses resulting in CpGs methylation levels for each treatment group (5 samples/treatment group). Additionally, we identified differentially methylated regions (DMRs) within a minimum range of 500 bp across different treatment groups. We performed empirical optimization of the methylation regions generated using Bismark across different treatment groups. Additionally, we performed dependency adjustment for DMRs within each treatment group and identified statistically significant differentially methylated genes (DMGs) in the BSp, GTPs and combination treatment group.
In the BSp and GTPs groups, a total of 2874 genes and 4074 genes associated with DMGs were identified (p value ≤ 0.05), respectively. Out of 2874 genes in the BSp group, 622 genes were hypomethylated and 2252 genes were hypermethylated, and out of 4074 genes in the GTPs group, 536 genes were hypomethylated and 3538 genes were hypermethylated (Fig. 5a). Comprehensive lists of all the transcripts that are DM in the different treatment groups are provided in Supplementary File 2: Tables S6-S12. Based on the annotation results, the DMRs in the BSp treatment group were mainly distributed in the intronic regions (35%), exonic regions (18%), intergenic regions (31%) and promoter regions (16%) (Supplementary Fig. S5a). The DMRs in the GTPs treatment group were found in the intronic regions (38%), exonic regions (18%), intergenic regions (32%) and promoter regions (17%) (Supplementary Fig. S5b). Unlike BSp and GTPs treatments administered individually, the combination treatment exhibited greater variation at the methylation level with a total of 4181 genes associated with DMGs. Out of 4181 genes in the combination treatment group, 542 genes were hypomethylated and 3639 genes were hypermethylated (Fig. 5a). In the combination treatment group, the DMRs were distributed in the intronic regions (35%), exonic regions (16%), intergenic regions (32%) and promoter regions (17%) (Fig. 5b).
This list served (Supplementary File 2: Table S10) as a reference list for identification of DMRs in the combination treatment group and was further used for correlation with transcripts that were DE in the combination treatment group.  Fig. S6a). The correlation of DMGs with DEGs was not performed in the GTPs treatment group as there we no DEGs identified in GTPs treatment group. Subsequently, a heatmap between DM and DE in BSp treatment group was generated in order to visualize the transcription and methylation level changes between these overlapping genes ( Supplementary Fig. S6b). Out of 8 genes which were DM and DE, 4 genes were up-regulated and hypomethylated and 4 genes were down-regulated and hypermethylated (Supplementary File 2: Table S8). Unlike the BSp treatment group, the combination treatment group exhibited a higher correlation among DEGs and DMGs. In the combination treatment group, out of 895 DEGs with 575 up-regulated genes and 320 downregulated genes, and 4181 DMRs with 542 hypomethylated and 3639 hypermethylated regions, we identified 39 overlapping genes (Fig. 6a). Table 3 provides a comprehensive depiction of overlapping genes in the combination treatment group which were DE and DM simultaneously. Subsequently, 2 genes (Pdx1 and Tmem132d) were identified with a positive association of up-regulated DEGs and hypomethylated genes, while 3 genes (Get4, Rpl13 and Ndufa1) were found to be down-regulated DEGs and hypermethylated genes. To better visualize the relationship between the 5 transcripts that were DM (meth.diff) and DE, we generated a scatter plot between methylation difference and FC (Fig. 6b). The heatmap in Fig. 6c shows DE and DM between rows and columns corresponding to biological replicates in the combination treatment group. Table 4 provides a reference list of unique transcripts that showed significantly differential changes (p value ≤ 0.05) with positive correlation of gene expression and DNA methylation patterns (up-regulated gene expression with hypomethylated DNA and down-regulated gene expression with DNA hypermethylated) in the combination treatment group.
A pathway enrichment analyses using ConsensusPathDB was performed to explore the pathways associated with the overlapping genes. The up-regulated and hypomethylated genes were enriched for "transcriptionrelated", "DNA binding", "regulation of gene expression" and "cell division" pathway. The down-regulated and www.nature.com/scientificreports/ hypermethylated genes were enriched with "apoptosis", "oxidative-phosphorylation", "translation" and "DNAmethylation" pathway (Table 5). These results suggest that the molecular mechanisms that reinforce the biologically beneficial effects of lifestyle modification by consumption of BSp + GTPs together drive transcriptome level changes significantly and have a substantial effect on methylation level changes.

Validation of target genes with unique differentially expressed (DE) and differentially methylation (DM) patterns in combination treatment by qRT-PCR.
To further verify the identified target genes that showed significant changes in DE and DM patterns in response to the combination treatment, we evaluated the expression of Pdx1, Tmem132d, Get4, Rpl13 and Ndufa1 genes in mouse mammary tumors using qRT-PCR. We found that combinatorial treatment induced different mRNA expression levels in all tested genes as shown in Fig. 7a- (Fig. 7c-e). Consistent with our genome-wide analysis, these results indicate that combined BSp + GTPs treatment are highly effective in inhibiting early BC in comparison to BSp treatment group and GTPs treatment group alone, which may contribute to epigenetic regulation of these key genes.

Discussion
Novel therapeutic approaches and targets have been developed for different types of cancer 19 . For instance, a recent study on molecular signatures identified various mutant proteins related to BC such as BRCA1, BRCA2 and PTEN. These molecular signatures can potentially serve as a potential target for understanding various The height in the bar plot represents the total number of differentially expressed genes. (d) Gene ontology enrichment terms and REACTOME pathways analyses using DAVID web-based tool. The plot is sorted based on decreasing FC wherein Y-axis represents specific GO terms related to biological pathways and X-axis represents log 10 (FC) associated with each GO term. www.nature.com/scientificreports/ changes at the genomic level during breast tumor progression 20 . Significant evidence indicates the role of various dietary phytochemicals or compounds in cancer prevention by elucidating changes at the molecular and biological levels 21,22 . Previous studies have also suggested that several mechanisms may lead to various preventive and therapeutic effects of dietary botanicals on different types of cancers such as apoptosis 23 , cell cycle arrest 24 and regulation of various signaling pathways such as MYC-WWP1 inhibitory pathway targeting PTEN 25 .
In our previous studies, we found that BSp 26,27 and GTPs 28,29 can potentially act as therapeutic and preventive agents against BC. However, the molecular mechanisms through which these dietary botanicals (BSp and GTPs alone and in combination) affect different stages of life remain unknown.
In this study, we evaluated the effect of BSp and GTPs alone and in combination in prevention of ER(−) mammary cancer. Herein, we administered BSp, GTPs and a combination of BSp and GTPs to transgenic Her2/neu mice. As a result, the combination treatment group showed stronger efficacy towards inhibiting tumor growth. www.nature.com/scientificreports/ www.nature.com/scientificreports/ Furthermore, the tumor incidence rates demonstrated a significant decrease in the BSp and combination treatment groups along with a decline in tumor weight. To better understand the underlying mechanisms behind the impacts of these dietary botanicals on BC prevention or delay, we characterized the global DNA methylation and gene expression patterns in mammary tumors in response to BSp, GTPs or in their combination (BSp + GTPs). This is the first report describing the impacts of BSp, GTPs or combination (BSp + GTPs) dietary treatment at both transcriptomic and methylomic levels in mammary tumors. Although high throughput analyses revealed various changes across different treatment groups, the combinatorial administration has a greater affect in comparison to BSp and GTPs administered individually. Gene set enrichment analyses identified different epigenetics pathways that were modulated due to these dietary treatments.
Firstly, we identified 145 DEGs and 895 DEGs that were obtained in response to the BSp treatment group or the combination treatment, respectively and no DEGs in GTPs treatment group. We identified DEGs between the control and the GTPs group, as well as other treatment groups, with the lower FDR-adjusted p value of 0.1 and higher fold change (higher than 2, in either direction). Although no DEGs were screened out from the GTPs group, numerous genes were up-or down-regulated with a fold-change less than 2 (Supplementary File 1: Table S6). Thus, significantly lower tumor biomass in GTPs group should be primarily attribute to an accumulative effect of those changed expressed genes. Additionally, posttranscriptional and posttranslational regulation, like noncoding RNA and histone modification, are proven to remarkably affect protein expression, leading to efficient regulation of biological processes, without noticeable changes at the transcriptional level. Future studies are warranted to investigate underlying mechanisms of mammary tumor suppression due to our dietary treatments at different regulatory levels. Subsequently, we identified a group of target genes (such as Pdx1, Tmem132d, Get4, Ndufa1 and Rpl13) that could be regulated by epigenetic mechanisms such as DNA methylation in the combination treatment group. Based on our results from DNA methylation profiles, the breast tumors treated with the combination diet (BSp + GTPs) tend to be more hypermethylated compared to the control group, which provides strong evidence that most of the CpGs are unmethylated in the control group. These DNA modifications at the global level may lead to changes in gene expression that are related to upregulation and down-regulation. Based on our analyses, we found considerably less hypomethylated regions in comparison to hypermethylated regions across different dietary treatment groups (Combination hypo = 542 and Combination hyper = 3639). Subsequently, in our previous studies, we also found that dietary BSp treatment potentially increased DNA methylation levels globally 11 .
We identified several target genes that show positive correlation of gene expression and DNA methylation changes through integrated analyses (Table 4). For example, the combination treatment with BSp and GTPs up-regulated gene expression of Pdx1 and Tmem132d genes. Pancreatic and duodenal homeobox 1 (Pdx1) is a transcription factor expressed in the promoter regions of Mus musculus and human olfactory epithelium (MOE). Pdx1 primarily contributes to morphogenesis during the development of mouse embryonic pancreas 30 . A plethora of studies have suggested that Pdx1 is a tumor suppressor gene which is primarily involved in different types of cancers such as BC, pancreatic cancer, gastric cancer and pseudopapillary tumors [31][32][33] . Further, we found that the Pdx1 gene was up-regulated and hypomethylated after administration of the combination (BSp + GTPs) treatment. However, Pdx1 was neither DE or DM after the administration of BSp or GTPs treatment alone. This could potentially be due to a synergistic effect in combination treatment in comparison to BSp or GTPs alone. Similarly, transmembrane protein 132d (Tmem132d) is a tumor suppressor gene actively involved in different types of cancer such as lung cancer 15 , ovarian cancer 34 and BC 35 . Tmem132d was up-regulated and hypomethylated in the combination treatment with BSp and GTPs but not in BSp or GTPs treatment alone. These gene expressions were found to be positively correlated with their DNA methylation changes with a methylation difference of − 27.08 and − 31.00. Functional annotations of these genes identified key pathways that were Table 5. Pathway enrichment analyses of overlapped genes between RRBS and RNA-seq data in the combination (BSp + GTPs) treatment group.

Gene symbol
Pathways related to gene symbol p value for differential expression (combination vs control) www.nature.com/scientificreports/ primarily related to "DNAtemplated transcription", "positive regulation of DNA binding", "regulation of gene expression" and "cell division". In addition, our analyses also revealed that combinatorial administration of BSp and GTPs resulted in downregulated and hypermethylated changes in the Ndufa1, Rpl13 and Get4 genes. The Ndufa1 gene (MWFE protein) is the first component of the electron transport chain that accepts electrons from NADH oxidation 36 . It is primarily responsible for NADH-ubiquinone oxidoreductase (complex-I) [37][38][39] . This gene is also responsible for direct effects on reactive oxygen species (ROS) 39 , functional complex-I assembly 40 and possess a wide variety of roles in aerobic activity 41 . These functionalities associated with Ndufa1 eventually play a significant role in genetic instability and tumorigenesis. In our analyses, we identified Ndufa1 to be downregulated and hypermethylated with a methylation difference of 28.91. Another gene, Rpl13 was found to be down-regulated in our study. Studies have shown Rpl13 plays a crucial role in tumor development in gastrointestinal carcinoma 42 , colorectal cancer 43 , prostate cancer and BC 44 with greater proliferative capacity and attenuated chemoresistance 45 . Additionally, functional ontology of these down-regulated genes revealed pathways that were primarily associated with various biological pathways such as "apoptosis", "oxidative-phosphorylation process", "DNA methylation activity" and "translation related pathway". Furthermore, a study in colorectal cancer reported that Golgi to ER traffic www.nature.com/scientificreports/ protein 6 (Get4) has been considered to promote tumorigenesis or tumor progression by demonstrating clinical significance of Get4 expression in colorectal cancer 16 . Subsequently, we validated the gene expression of these uniquely identified target genes by real time RT-PCR assay. Our results indicated that the combinatorial treatment with BSp and GTPs can induce significant gene expression in most of the tested genes, which are also consistent with our RNA-seq and RRBS results. Thus, combinatorial treatment with BSp and GTPs may induce epigenetic regulations of these key tumor-related genes with important biological functions, which may contribute to its tumor inhibitory effects on mammary cancer development consistent with the cancer phenotypic results that we have observed in this study.

Overlapping up-regulated DEGs and hypomethylated DMGs
Altogether, these findings, based on an unbiased analysis of DE and DM genes along with functional characterization, identified that combined administration of BSp and GTPs dietary botanicals have shown greater impact on tumor suppression, gene expression and DNA methylation levels in comparison to BSp or GTPs treatment alone. We posit that the alteration in gene expression coupled with methylation changes modulates epigenetics pathways to a greater extent with the combinatorial treatment. Overall, our results indicate that the combination treatment consisting of BSp and GTPs may reverse epigenetic aberrations in BC leading to beneficial outcomes such as slow progression or delay of breast tumorigenesis. Therefore, these latter results could facilitate further understanding of the underlying etiology behind these epigenetic alterations.

Conclusion
In summary, we found DE and DM genes across different treatment groups and detected a subset of correlated genes that were DE and DM by using an unbiased approach. Additionally, we show that in comparison to BSp or GTPs administered alone, combinatorial treatment consisting of BSp and GTPs has a greater impact on transcriptomic and methylomic changes that may contribute towards reducing the tumor incidence rate of ER(−) BC consistent with the cancer phenotypic observations in this study. Identifying these key genes that potentially contribute towards BC-associated risk factors could serve as a key avenue in the area of translational precision medicine and eventually lead to improving risk assessment of BC. Further, the functional characterization of these transcripts may lead to the identification of novel therapeutic strategies against ER(−) BC.

Materials and methods
Animals. The 46 . In order to generate experimental colonies, the breeder mice at 4 weeks of age were obtained from the Jackson Laboratory (Bar Harbor, ME) and were bred from 10 weeks of age. The pups were weaned at 21 days after birth and were further tagged and genotyped, and a standard PCR analyses was performed on the tail of mice 47 48 . Out of 20 mice in each treatment group, 5 mice from each treatment group were randomly chosen for transcriptomic and methylomic analyses. The first group was comprised of control mice (N Control = 5) which were administered control AIN-93G diet. In group 2 (N BSp = 5), mice were fed with SFNrich broccoli sprouts powder (Natural Sprout Co.) which was added at 26% (w/w) into a modified AIN93G diet pellet from TestDiet, St. Louis, MO. 5.13-6.60 μM SFN per gram was evaluated for each batch of broccoli sprouts powder. In group 3 (N GTPs = 5), mice were orally fed with GTPs Sunphenon 90D (Sunphenon 90D, Taiyo International, Inc., Minneapolis, Minnesota) which was added at 0.5% in drinking water. Sunphenon 90D (Taiyo Inc., Minneapolis, MN, USA) is comprised of > 90% polyphenols, > 80% catechins, > 45% EGCG and < 1% caffeine, and is a decaffeinated extract of green tea encompassing purified polyphenols rich in green tea catechins. In group 4 (N Combination = 5), mice were fed with a combination diet of both BSp diet (Group 2) in pellets and GTPs (Group 3). Treatments were continued throughout the study until termination. Subsequently, the mice mammary tumor tissues were collected at 31 weeks during termination.
Tumor observation and sampling. The female Her2/neu WT mice model developed mammary tumors at an early age at around 20 weeks. Tumor incidence was measured and recorded weekly for tumor development throughout the experiment and the experiment was terminated when the average tumor diameter of mice in the control group exceeded 1.0 cm. Mouse body weight was recorded biweekly from 4 to 28 wk of age. Mouse food and water intakes were measured at 4,12, 20 and 28 wk of age respectively. At the end of the experiment, mice were euthanized by CO 2 and the mammary tumors were excised, weighed and further stored in liquid nitrogen for subsequent analysis. www.nature.com/scientificreports/ of RNA-seq and RRBS data to generate transcript level abundance and genome-wide DNA methylation profile for different phytochemical treatments is graphically demonstrated in Fig. 2. The details of the analytical steps related to analyses are further described in the sections with various subparts.
DNA and RNA extraction. Genomic DNA (gDNA) was extracted from frozen mammary tumor tissues of mice using the DNAeasy kit (Qiagen, Valencia, CA) based on manufacturer's guidelines. DNA concentration was determined using a NanoDrop spectrophotometer and the DNA quantity was measured using Qubit. Subsequently, total RNA was extracted from frozen mammary tumor tissues of mice using TRIzol reagent (Sigma-Aldrich, St. Louis, MO) based on the manufacturer's protocol. The RNA concentrations were determined with NanoDrop spectrophotometer and RNA Nano bioanalyzer chip was used to evaluate the integrity of RNA. Only the RNA samples with an integrity number greater than 7 were further used for sequencing.
Statistical analyses on animal experiments. Power calculation and sample size was measured using one-side 2-proportion comparison using an online power calculator (http:// power andsa mples ize. com/). In order to avoid any chances of single false positive during multiple comparison, Bonferroni adjustment with Power = 80%, Significance level = 0.01 and α = 0.05 was performed amongst four different treatment groups and statistical analyses was performed by SPSS software (v24.0). Furthermore, tumor incidence and significance were determined using Chi-square test. Two-tailed student's t test was performed to compare two treatment groups, and one-way independent ANOVA was performed to compare three or more groups. Tukey's post-hoc test was performed to determine significance between the groups. Error bars were standard error of the mean obtained from experiments. Statistically significant outcomes were represented as ** for p value < 0.01 and * for p value < 0.05.

Sequencing, processing and alignment of RNA-seq and RRBS libraries. RRBS pair-end libraries
were sequenced on Illumina NextSeq 500 platform. The quality of the raw Fastq were assessed using FastQC (v0.11.4). The adapter sequences were trimmed using trim_galore based on NuGEN Ovation RRBS system. The trimmed reads were aligned to the reference genome for mouse GRCm38/mm10 using Bismark alignment with default parameter settings. As a result, context-dependent methylation (CpG sites) call files were generated using bismark_methylation extractor. RNA-seq pair-end libraries were sequenced on Illumina NextSeq500. We performed fastQC to check the quality of the raw fastq data per samples. After FastQC, the RNA-Seq fastq reads were aligned to the mouse reference genome GRCm38/mm10 using Kallisto with their default parameter settings. The aligned BAM files were further processed using Kallisto (v 0.43.1-intel-2017a) 49 . A transcription level estimates (or abundance estimate) file was generated for each sample/treatment group. Further, a transcript-level estimates file was used as an input in tximport package 50 in R (v3.6.1) and the transcript level information was summarized to gene-level with their respective expression values.
Identification of differentially methylated (DM) CpG sites and differentially methylated regions (DMRs). Differential methylation analyses was performed using methylKit 51 package in R (version 3.6.1). Firstly, the methylation call files generated using Bismark 52 for treatment groups (BSp, GTPs, BSp + GTPs) were used as input files to generate CpG regions profiles. Subsequently, methylKit was used to identify DMG's and DMR's based on the false discovery rate (FDR) ≤ 0.05. In order to further explore the methylation profiles between the control and the treatment groups (BSp, GTPs and BSp + GTPs), hierarchical clustering was performed in R using hclust package.

Identification of differentially expressed genes (DEGs).
To characterize the transcriptional level changes in malignant tumors of Her2/neu mice (5 mice per treatment group) after BSp treatment, GTPs treatment, combined GTPs + BSp treatment, and control group (untreated group), we utilized R/Bioconductor package Limma (version 3.6.1). Limma package was used to evaluate differential gene expression by performing quantile normalization wherein an attempt is made to match gene count distributions across the samples in a dataset 53 . Subsequently, the significant threshold of DEGs was set as|log 2 (fold-change)|> 2 and false discovery rate (FDR) ≤ 0.01.

Integrated analyses of DNA methylation and gene expression.
To further elucidate the potential relationship between DNA methylation and mRNA expression, we correlated CGI-DMRs with their respective DEGs by selecting the significant DMRs-DEGs pair with p value < 0.05. In order to visualize the relationship between DEGs and DMRs, scatterplot with FC was generated using ggplot package in R (version 3.6.1) 54 .
Gene set function enrichment. Significant transcripts at the transcriptome level and methylome levels (p < 0.05) were considered to perform Gene Ontology (GO) enrichment. To identify associations of genes with GO terms, functional enrichment was conducted using web-based tool DAVID 55 and gsva package comprised of Reactome and KEGG database in R (version 3.6.1). Furthermore, we used a web-based tool WebGetSalt to perform an over-representation analyses with a significance level of 5% FDR 56 .
Quantitative real-time PCR (qRT-PCR). Total RNA was harvested from mouse mammary tissues by using QIAGEN RNeasy Plus Kit. Following extraction following the manufacturer's protocol. Following the extraction, RNA was reverse transcribed to cDNA by using iScript cDNA synthesis kit (BioRad) according to the www.nature.com/scientificreports/ manufacturer's instructions 57,58 . Gene expression for specific genes of interest was performed in triplicate and further analyzed using real-time PCR by SYBR GreenER qPCR Supermix (Invitrogen) in a Roche LC480 thermocycler. For the PCR setup, we used 1 µL of cDNA, 5 µL of iTaq SYBR green from Bio-Rad, 2 µL of nucleasefree water, 1 µL of forward and reverse primers for specific genes of interest with total volume of 10 µL. Once the www.nature.com/scientificreports/