Possible epigenetic regulatory effect of dysregulated circular RNAs in Alzheimer’s disease model

As circular RNAs (circRNAs) regulates the effect of micro RNAs (miRNAs), circRNA–miRNA-mRNA network might be implicated in various disease pathogenesis. Therefore, we evaluated the dysregulated circRNAs in the Tg2576 mouse Alzheimer’s disease (AD) model, their possible regulatory effects on downstream target mRNAs, and their pathomechanistic role during the disease progression. The microarray-based circRNA expression analysis at seven- and twelve-months of ages (7 M and 12 M) returned 101 dysregulated circRNAs at 7 M (55 up-regulated and 46 down-regulated) and twelve dysregulated circRNAs at 12 M (five up-regulated and seven down-regulated). For each dysregulated circRNA, potential target miRNAs and their downstream target mRNAs were searched. Dysregulation of circRNAs was associated with increased frequency of relevant dysregulation of their downstream target mRNAs. Those differentially expressed circRNA–miRNA-mRNA regulatory network included 2,275 networks (876 for up-regulated circRNAs and 1,399 for down-regulated circRNAs) at 7 M and 38 networks (25 for up-regulated circRNAs and 13 for down-regulated circRNAs) at 12 M. Gene ontology (GO) and pathway analyses demonstrated that the dysregulated mRNAs in those networks represent the AD pathomechanism at each disease stage. We concluded that the dysregulated circRNAs might involve in the AD pathogenesis by modulating disease relevant mRNAs via circRNA–miRNA-mRNA regulatory networks.

circRNA might complete with that of the corresponding linear mRNAs. Second, circRNAs regulates the expression of gene transcription machineries. Most importantly, circRNAs contains multiple binding sites for miRNAs that enable to sequestrate or buffer the effect of the target miRNA. In this regard, circRNA might modulate the expression of their target genes by participating in circRNA-miRNA-mRNA regulatory network [25][26][27] .
Some properties of circRNA indicate that circRNAs might have a particular role in the pathogenesis of AD. First, the regulation of circRNA expression shows a time-and region-specific pattern and is independent from that of the corresponding linear-form mRNAs 28,29 . Second, circRNAs are highly abundant and stable in the brain due to their closed loop structures making them resistant to RNA exonucleases or RNase R-mediated degradation 21,29 . Therefore, circRNAs can consistently buffer the highly fluctuating effect of miRNA and their altered expression might sufficiently direct the overall gene expression profile into a certain disease process via circRNA-miRNA-mRNA networks 22 . In this context, studying the altered profile of circRNA expression in different disease stages and its implication on the downstream target mRNAs might help elucidate a novel epigenetic pathomechanism of AD.
In this study, we hypothesized that changes in circRNA expression during the AD progression might be implicated in the pathogenesis of AD by epigenetic regulation via circRNA-miRNA-mRNA network. Therefore, we investigated the association among the dysregulated circRNAs and the altered expression profiles of their downstream target miRNAs and mRNAs in different disease stages of AD by constructing circRNA-miRNA-mRNA network based on the microarray database and evaluated their potential role in the pathogenesis of AD by bioinformatics analyses.

Results
Overall expression profile of mRNA, miRNA, and circRNA. Among a total of 39,429 mRNAs analyzed, the number of differentially expressed mRNAs in the brain of 7 M Tg2576 mouse was 1,108 (up-regulated, 417; down-regulated, 691) and 12 M mouse was 1,534 (up-regulated, 1,121; down-regulated, 413). Scatter plots for the mRNA expression pattern were demonstrated in Fig. 1A. In most cases, dysregulation of mRNAs was time specific, not persisting or being reversed between 7 M and 12 M (2186/2414, 90.6%, Fig. 2A). Top five most In the volcano plot, red squares represent the differentially expressed circRNAs (P-value < 0.05). Panel D shows hierarchical cluster analysis of differentially expressed circRNAs. The log2 signal intensity is reflected in the color scale, which runs from green (low intensity) to red (strong intensity). dysregulated mRNAs in 7 M and 12 M Tg2576 mice were listed in Table 1. The full data of mRNAs expression is available in Supplemental Data S1.
The total number of analyzed miRNAs was 326. Among them, the number of differentially expressed miRNAs in the brain of 7 M Tg2576 mouse was 163 (up-regulated, 135; down-regulated, 28) and 12 M mouse was 121 (up-regulated, 110; down-regulated, 11). The majority of the miRNA dysregulation was time specific (154/219, 70.3%, Fig. 2B). The full data of miRNA expression is demonstrated in Supplemental Data S2.
The total number of analyzed circRNAs was 14,119. Among them, the number of differentially expressed circRNAs in the brain of 7 M Tg2576 mouse was 101 (up-regulated, 55; down-regulated, 46) and 12 M mouse was 12 (up-regulated, 5; down-regulated, 7). Scatter plots, volcano plots, and hierarchical clustering for the circRNA expression pattern were demonstrated in Fig. 1B-D. Most of the dysregulation of circRNAs was also time specific (109/111, 98.2%, Fig. 2C). Top five most highly dysregulated circRNAs in 7 M and 12 M Tg2576 mice were listed in Table 2. The full data of circRNA expression is in the Supplemental Data S3.  PCR analysis of circRNA and mRNA expression. Three RNAs were randomly selected from each of the top ten (five up-regulated and five down-regulated) most highly dysregulated circRNAs and mRNAs in 7 M and 12 M Tg2576 mice. Quantitative polymerase chain reaction (PCR) analysis of those 12 RNAs (six circRNAs and six mRNAs) was performed using the primers described in the Supplemental Data S4. As the result, 5/6 (83.3%) circRNAs and 5/6 (83.3%) mRNAs were confirmed to be differentially expressed relevantly to the microarray data ( Fig. 3, see Supplemental Data S5 for the PCR raw data).
Analysis of circRNA-miRNA-mRNA regulatory interaction. First, we evaluated the regulatory effect of miRNAs on their target mRNAs, separately for the up-and down-regulated miRNAs in each time point (7 M and 12 M). 10,243 mRNAs were identified as the potential targets for the 135 up-regulated miRNAs in the 7 M brain, 3,100 mRNAs for the 28 down-regulated miRNAs in the 7 M brain, 9,475 mRNAs for the 110 up-regulated miRNAs in the 12 M brain, and 1,587 mRNAs for the 11 down-regulated miRNAs in 12 M the brain. However, no significantly higher frequency of dysregulation was observed in the downstream target mRNAs of the differentially expressed miRNAs (Table 3A). Rather, the frequencies of down-regulation in the target mRNAs of up-regulated miRNAs in the 12 M brain and up-regulation in the target mRNAs of down-regulated miRNAs in the 7 M brain were significantly low. Second, we evaluated the regulatory effect of differentially expressed circRNAs on the expression of their target miRNAs. From the MRE sequence analyses, 205 miRNAs were identified as the potential targets for the 55 up-regulated circRNAs in the 7 M brain, 184 miRNAs for the 46 down-regulated circRNAs in the 7 M brain, 27 miRNAs for the 5 up-regulated circRNAs in the 12 M brain, and 32 miRNAs for the 7 down-regulated circRNAs in the 12 M brain. From the microarray data of 326 miRNAs, the expression of 23/205, 21/184, 3/27, and 3/32 of those potential target miRNAs, respectively, were identified. However, no significantly altered frequency of dysregulation was observed in the downstream target miRNAs of the differentially expressed circRNAs, supporting that the regulatory effect of circRNAs on their target miRNAs does not result in the alteration of the target miRNA expression level (Table 3B).
Third, we evaluated the regulatory effect of differentially expressed circRNA-miRNA-mRNA regulatory networks on their target mRNAs. 5,406 mRNAs were identified as the potential targets for the 23 target miRNAs of the 55 up-regulated circRNAs in the 7 M brain, 5,266 mRNAs for the 21 target miRNAs of the 46 down-regulated circRNAs in the 7 M brain, 787 mRNAs for the 3 target miRNAs of the 5 up-regulated circRNAs in the 12 M brain, and 65 mRNAs for the 3 target miRNAs of the 7 down-regulated circRNAs in the 12 M brain. Notably, the frequency of dysregulation in the relevant direction of the downstream target mRNAs was significantly high in every group of different time and direction of upstream circRNA expression (Table 4A). After adjusting the effect of dysregulated upstream miRNAs, the altered expression of upstream circRNAs was significantly associated with the higher frequency of relevant dysregulation in downstream target mRNAs (P < 0.001 for up-regulated   Table 4B).
Differentially expressed circRNA-miRNA-mRNA regulatory network. Construction of the cir-cRNA-miRNA-mRNA regulatory networks was performed separately for the up-and down-regulated circRNAs in each time point (7 M and 12 M). For the 55 up-regulated circRNAs and 23 downstream target miRNAs in 7 M brain, 242 downstream mRNAs were relevantly up-regulated, consisting a total of 876 circRNA-miRNA-mRNA regulatory network (Fig. 4A). For the 46 down-regulated circRNAs and 21 downstream target miRNAs in 7 M brain, 388 downstream mRNAs were relevantly down-regulated, consisting a total of 1,399 circRNA-miRNA-mRNA regulatory network (Fig. 4B).
For the 5 up-regulated circRNAs and 3 downstream target miRNAs in 12 M brain, 20 downstream mRNAs were relevantly up-regulated, consisting a total of 25 circRNA-miRNA-mRNA regulatory network (Fig. 4C). For the 7 down-regulated circRNAs and 3 downstream target miRNAs in 12 M brain, 9 downstream mRNAs were relevantly down-regulated, consisting a total of 13 circRNA-miRNA-mRNA regulatory network (Fig. 4D). The full list of differentially expressed circRNA-miRNA-mRNA regulatory network is in the Supplemental Data S6).
Additionally, a substantial complexity was observed in the differentially expressed circRNA-miRNA-mRNA regulatory networks, in which each miRNA was targeted by up to three upstream circRNAs (median 1, interquartile range, IQR 1−1), targeted up to 123 downstream mRNAs (median 32, IQR 18−54.25), and each mRNA was targeted by up to 9 upstream miRNAs (median 2.5, IQR 1−7).   www.nature.com/scientificreports www.nature.com/scientificreports/ Gene ontology and pathway analysis. Gene ontology and pathway analyses were performed separately for the up-and down-regulated circRNA-miRNA-mRNA regulatory networks in each time-point. In each analyses, the top five enriched GO terms and KEGG pathways were summarized in Tables 5 (up-regulated mRNAs) and 6 (down-regulated mRNAs, full list of the enriched GO terms and KEGG pathways in Supplemental Data S7).
The most enriched terms of the dysregulated target mRNAs were indicative of that mRNA regulation has directed pattern which are related with the AD pathogenesis. In the 7 M brain, the most enriched GO and pathway terms for the up-regulated mRNAs are related with the immune activation ('Macrophage colony-stimulating factor signaling pathway' in BP, 'Fc-gamma receptor I complex binding' in MF, and 'CSF1-CSF1R complex' in CC), activation of inflammatory cascade ('Complement and coagulation cascades' , 'Staphylococcus aureus infection' , and 'Hematopoietic cell lineage' in Pathway analysis), cellular adhesion ('Integrin complex' in CC), and production of reactive oxygen species ('Superoxide-generating NADPH oxidase activity' in MF and 'NADPH oxidase comple'x in CC) which are known as the major responses to overproduced Aβ in brain 2,3,[9][10][11]30 . Also in the 12 M brain, most of the enriched terms for the up-regulated mRNAs are related with the immune acitivation ('T cell activation via TCR contact with antigen bound to MHC molecule on APC' and 'Negative regulation of natural killer cell activation' in BP, 'C-X3-C chemokine binding' in MF, and 'Iimmunological synapse' and 'Phagocytic vesicle membrane' in CC), inflammatory response ('Positive regulation of tumor necrosis factor secretion' and 'Positive regulation of interleukin-1 production' in BP and 'Graft-versus-host disease' and ' Allograft rejection' in pathway analysis), and cellular adhesion ('Lipoteichoic acid binding' , 'Fibronectin binding' , and 'Integrin binding' in MF, and 'Integrin complex' in CC).
In contrast, the most enriched GO and pathway terms for the down-regulated mRNAs in the 7 M brain are commonly related with progenitor self-renewal and neuronal differentiation ('Muscle cell fate commitment' , 'Negative regulation of muscle hyperplasia' , and 'Schwann cell proliferation' in BP) and maintenance ('Regulation www.nature.com/scientificreports www.nature.com/scientificreports/ of store-operated calcium entry' and 'Mitochondrial genome maintenance' in BP, 'Structural constituent of bone' and 'cAMP response element binding protein binding' in MF, 'Trans-Golgi network membrane' and 'Transcription elongation factor complex' in CC, and 'Pantothenate and CoA biosynthesis' , ' Arginine and proline metabolism' , and 'Glucagon signaling pathway' in pathway analysis). Additionally, the down-regulated mRNAs in 12 M brain were enriched in synapse function and preservation of neuronal networks ('Positive regulation of action potential' in BP, 'Inward rectifier potassium channel activity' and 'Voltage-gated ion channel activity' in MF, 'Voltage-gated potassium channel complex' and 'Receptor complex' in CC, and 'Retrograde endocannabinoid signaling' and 'Cholinergic synapse' in pathway analysis), which is a competent finding with the overt manifestation of AD phenotypes at twelve months 31 .

Discussion
This study demonstrated the altered expression of circRNAs and their possible epigenetic regulatory effect at different time points in the brain of an AD model. Although a very small number of circRNAs had a significantly altered expression (0.7% in 7 M and 0.1% in 12 M), they targeted large number of mRNAs by circRNA-miRNA-mRNA interaction. Dysregulated circRNAs were associated with a significantly higher frequency of the relevant dysregulation of the downstream target mRNAs, after adjusting the effect of dysregulated upstream miRNAs. Additionally, gene ontology and pathway analyses demonstrated that the mRNAs of the differently expressed circRNA-miRNA-mRNA interactions might be closely involved in the pathomechanism of AD. Numerous studies have evaluated the miRNAs' epigenetic regulatory effect on disease relevant genes 1,2,[8][9][10][11][12][13][14][15] . However, this study demonstrated that circRNAs might have more significant regulatory effect on the disease relevant genes than miRNAs.
In this study, we found no significant change in the frequency of dysregulation in the downstream target mRNAs of the differentially expressed miRNAs. This negative association might be due to the high complexity of miRNA-mRNA interactions and the high temporal variability of miRNA expression levels in the brain [16][17][18][19] . Additionally, no significant change was observed in the frequency of dysregulation in the downstream target  www.nature.com/scientificreports www.nature.com/scientificreports/ miRNAs of the differentially expressed circRNAs. This might be attributable to that the regulatory effect of circR-NAs on miRNAs might not necessarily result in an alteration of target miRNA expression level 25 . Some properties of circRNAs might explain the regulatory effect of circRNA on target mRNAs that exceeds the effect of miRNA. First, due to substantial complexity in the circRNA−miRNA−mRNA interactions, a single dysregulation of a miRNA can hardly induce a significant dysregulation in downstream target genes, as they are simultaneously regulated by hundreds of other upstream RNAs 18,19 . Second, as circRNAs are more stable than miRNAs in CNS 21,29 , their altered expression can sufficiently buffer the fickle changes of miRNAs' effect and modulate the overall gene expression profile in response to a certain disease process [20][21][22] .
No circRNAs except for the mmu_circRNA_29980 was consistently dysregulated in 7 M and in 12 M brains, suggesting a high time specificity in the regulation of the circRNA expression, which is consistent with the findings of the previous studies. This might indicate that circRNA expression is dynamically regulated according to each stage of the disease progression to exert an epigenetic regulatory effect on their downstream targets relevantly to the disease pathomechanism 28,29 .
Gene ontology and pathway analyses demonstrated that the mRNAs of the differently expressed circRNA-miRNA-mRNA regulatory networks might be closely involved in the pathomechanism of AD. The enriched terms of up-regulated target mRNAs in both 7 M and 12 M brains were commonly associated with the immune activation, activation of inflammatory cascade, and cellular adhesion. Considering that Tg2576 mice show progressive accumulation of amyloid-b42 (Aβ42) in brain which becomes evident at seven months, these up-regulated mRNAs might be involved in the major responses to overproduced Aβ in brain 2,3,[9][10][11]30 . In contrast, the enriched terms of down-regulated target mRNAs in both time points were associated with the progenitor self-renewal, neuronal differentiation and maintenance. These findings represent the ongoing deterioration of metabolic homeostasis and subsequent neuronal degeneration during the disease progression 2,3,8,10,32 . Additionally, the enriched terms of down-regulated target mRNAs 12 M were associated with synapse function and preservation of neuronal  www.nature.com/scientificreports www.nature.com/scientificreports/ networks, which is competent with the overt manifestation of AD phenotypes at twelve months 3,4,8,31 . Therefore, the results of the gene ontology and pathway analyses suggests that the altered expression of circRNAs might direct the overall gene expression profile relevantly to a certain disease process in cellular, network, and phenotypical (cognitive) levels via circRNA-miRNA-mRNA networks. In this regard, the dysregulated circRNAs might involve in the AD pathogenesis in disease stage specific manners.
There are several limitations in this study that should be acknowledged. First, because of the large number of genes analyzed, validation of their expression ratios with quantitative PCR analyses was performed for the limited number of the differentially expressed circRNAs and mRNAs. Second, we only evaluated the association between the dysregulated circRNA, miRNA, and mRNAs and did not directly validate the regulatory interactions among them. Especially, quantitative analyses of the association between the changes of the MRE expression in dysregulated circRNAs and the expression levels of their target mRNAs were not performed, due to the high complexity of circRNA-miRNA-mRNA interactions. Third, the role of the circRNA−miRNA−mRNA regulatory network in the AD pathomechanism was not validated but only predicted by bioinformatics methods. Future studies should endeavor to confirm the regulatory interaction among the circRNA, miRNA, and mRNA and network verify a specific circRNA−miRNA−mRNA regulatory network with a sufficient pathomechanistic role in AD, which might serve as a disease marker and potential therapeutic target.

Materials and Methods
Study design. Tg2576 AD transgenic mice (Taconic, Hudson, NY) were used in this study, which highly express the human 695-amino acid isoform of amyloid precursor protein containing the Swedish double mutation causing early-onset AD. Tg2576 mice show increased amyloid-b42 (Aβ42) expression in brain at seven months and manifests AD phenotypes at twelve months 1,31 . Hence, four male Tg2576 mice with seven months of age (7 M) and four male with twelve months of age (12 M) were chosen and compared with the normal male C57BL/6 mice. All animals were managed with standardized procedures approved by the Institutional Animal Care and Use Committee of Seoul National University Hospital.
Tissue preparation and microarray. Animals were sacrificed and their brain samples were taken. Control mice of the same number were sacrificed at the same age as the Tg2576 mice. The brains were obtained from each mouse and immediately stored at −80 °C. For microarray analysis, total RNAs were extracted using TRIZOL reagent (Invitrogen, NY, USA) and purified by an RNeasy Mini Kit (Qiagen, Hilden, Germany). RNA quantity was measured with a Nanodrop ND-1000 (Thermo Fisher Scientific, MA, USA) and quality checked by an Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA) 1 .
For circRNA microarray, the RNAs were treated with RNase R (Epicenter, WI, USA) to remove linear RNAs and enrich circRNA. The enriched circRNA samples were amplified and transcribed into fluorescent cRNA using a random priming method using Arraystar Super RNA Labeling Kit (Arraystar, Rockville, MD, USA). The labeled cRNAs were hybridized onto the Arraystar Mouse circRNA Array V2 (8 × 15 K, Arraystar). After the slides were washed, the arrays were scanned by an Agilent Scanner G2505C (Agilent Technologies). The miRNA and mRNA microarray data were also obtained using the Agilent Mouse miRNA Microarray 8 × 15 K kit and Agilent Mouse Gene Expression Microarray 4 × 44 K kit respectively, according to the manufacturer's protocol (Agilent Technologies) 1 .
Agilent Feature Extraction software (Agilent Technologies) was used to analyze the acquired array images. Quantile normalization and subsequent data processing were performed with the Agilent's GeneSpring Software (Agilent Technologies). Mann-Whitney test was used to detect differentially expressed circRNAs and mRNAs between the Tg2576 mice and the controls by fold-changes of ≥1.5 and P-values of < 0.05 with relatively small false discovery rates (FDRs). Scatter plots (for circRNAs and mRNAs), volcano plot, and hierarchical clustering (for circRNAs) were used to demonstrate the expression pattern of circRNAs and mRNAs.
Analysis of circRNA-miRNA-mRNA regulatory interaction. The regulatory effect of miRNAs on their target mRNAs was separately evaluated for the up-and down-regulated miRNAs in each time point (7 M and 12 M). For each dysregulated miRNA, potential target mRNAs were predicted by the combination of TargetScan and miRanda 14,22 . Expression profiles of those potential target mRNAs were obtained from the microarray data and the change in the frequency of relevant dysregulation (down-regulation of target mRNAs for up-regulated miRNAs or up-regulation of target mRNAs for down-regulated miRNAs) in those target mRNAs were evaluated.
To evaluate whether the dysregulate circRNAs alter the expression level of their target miRNAs, up to five target miRNAs for each differentially expressed circRNA were identified by microRNA response element (MRE) analysis for the circRNA sequence. MRE sequence analysis was performed using miRNA target prediction software (Arraystar) based on TargetScan (http://www.targetscan.org) and miRanda (http://www.microrna.org) algorithms 18,26 . Using the microarray data, the change in the frequency of relevant dysregulation in those target miRNAs were evaluated.
To demonstrate whether the dysregulated circRNAs alter the expression level of their downstream target mRNAs by circRNA-miRNA-mRNA interactions. We found that the dysregulated circRNAs does not alter the frequency of relevant dysregulation in their target miRNAs [25][26][27] . However, as circRNAs inhibit or buffer the miRNA's effect of suppressing its target mRNAs, we speculated that an up-regulated circRNA might exert a disinhibitory effect on its downstream targets mRNAs and a down-regulation of a circRNA might induce a hyper-inhibition on its downstream target mRNAs 17,18 . Therefore, we evaluated whether dysregulated circRNAs can induce a significant change in the frequency of relevant dysregulation in the downstream target mRNAs 25,26 . www.nature.com/scientificreports www.nature.com/scientificreports/ PCR analysis. Based on the microarray data, three circRNAs and three mRNAs were randomly selected from each of the top ten (five up-regulated and five down-regulated) most highly dysregulated circRNAs and mRNAs in 7 M and 12 M. For PCR analysis, two brain samples were pooled into one RNA sample as a unit. cDNAs were synthesized from 0.5 μg of total RNA of brain tissues by reverse transcription. Standard curves were prepared using 2 × SuperArray PCR master mix (Arraystar) according to the manufacturer's protocol. The relative expression ratio of each circRNA and mRNA was calculated with the Rotor-Gene Real-Time Analysis Software 6.0 (Qiagen), using the housekeeping gene, Gapdh, expression for normalization. All real-time reactions were performed in triplicate 27 .
Specific circRNA-miRNA-mRNA regulatory network. We found that the regulatory effect of circR-NAs on their target miRNAs does not changes the frequency of the target miRNA dysregulation but is associated with a higher frequency of relevant dysregulation in the downstream target mRNAs 25 . Therefore, we constructed differentially expressed circRNA-miRNA-mRNA regulatory networks comprised of a dysregulated circRNA, its target miRNA, and miRNA's downstream target mRNA which is also dysregulated in the same direction with the upstream circRNA dysregulation, without considering whether the miRNA expression level was altered [25][26][27] . Additionally, visualizations of the differentially expressed circRNA-miRNA-mRNA regulatory networks were performed using Cytoscape 3.4.0 33 .

Gene ontology and pathway analysis.
To demonstrate the pathophysiologic role of the differentially expressed circRNA-miRNA-mRNA regulatory network in AD pathophysiology, gene ontology and pathway analyses were performed for the target mRNAs. The gene ontology domains included Biological process (BP), Molecular function (MF), and Cellular component (CC) and were obtained using Database for Annotation, Visualization and Integrated Discovery (DAVID; http://www.david.abcc.ncifcrf.gov/) 34,35 . Pathway analysis was performed using the Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg) database 36 . In both analyses, Fisher's exact or chi-squared test with FDR were used, where a GO term or KEGG pathway with P-value < 0.05 and FDR < 0.05 was considered statistically significant. The top five enriched GO terms and pathways of the differentially expressed target mRNAs were ranked by fold enrichment score.
Statistical analysis. Data were reported as number (percentage) or mean ± standard deviation. Mann-Whitney U test was used to detect differentially expressed circRNAs, miRNAs, and mRNAs between the two groups by fold-changes of ≥1.5 and P-values ≤ 0.05. Changes in the frequency of relevant dysregulation in the downstream target RNAs were evaluated using Pearson's chi-square test. To evaluate factors associated with significant dysregulation of mRNA, logistic regression analyses including parameters with P-values < 1.20 in univariate analyses were performed. SPSS (version 23.0; SPSS Inc., Chicago, IL, USA) was used to all statistical analyses. P-values < 0.05 were considered statistically significant.

Data Availability
Full microarray, gene ontology, and pathway analysis analysis data of this study is available in the Supplemental Datasets.