Diet preparation and study design of the rat feeding trials were described in a previous publication (Zeljenkova et al., 2014) following OECD test guideline 408. Briefly, intestinal tissue was collected after two repeated-dose 90-day oral toxicity studies in rats (referred to as studies A and B). The animal trials were performed in compliance with GLP in the experimental animal house at the Department of Toxicology of the Slovak Medical University in Bratislava and approved by the ethics commitee of the Veterinary State Administration, Slovak Republic Ro-4372/12-221. The two studies used maize materials with the MON810 event but with different varieties having different genetic backgrounds (study A = MON810 from Monsanto; study B = MON810 from Pioneer). For RNA analysis, we used intestinal tissue (ileum, ascending colon) from high dose study groups, which were fed diets containing 33% GM or 33% near-isogenic control maize. For study A, male and female animals were sampled on four consecutive days, each day starting with animals from the control group and finishing with animals from the GM-maize fed group. Sampling was performed over a period of 10 h. A similar sampling scheme was applied in study B, but the interval between sacrifice of the first and last animal was 8 h. Due to the high number of animals per group a variation in time of sacrifice was given as a total of 160 rats were introduced in a randomized block design.

Tissue samples from 10 animals per group and sex were used (n = 10). In total, 20 male and 20 female rats were examined, one half being fed the 33% GM maize feed and one half being fed the 33% isogenic maize feed. The tissue samples (ileum, colon) were dissected after sacrifice, immediately frozen in liquid nitrogen and stored at -80°C. Total RNA from total ileum and colon macro-dissections (∼250 mg wet-weight) was isolated using the mirVana^TM miRNA Isolation Kit (Ambion), according to the manufacturer’s protocol. The RNA quantity was determined with the Nanodrop 1000 Spectrophotometer (Thermo) whereas the quality was controlled as previously described via the 2100 Bioanalyzer Instrument (Agilent Technologies) (Sharbati et al., 2010). Quantification of mRNA expression was performed by RT-qPCR as described earlier (Sharbati et al., 2011), with some modifications. One microgram of total RNA was reverse-transcribed using the RevertAid^TM M-MuLV Reverse Transcriptase (Fermentas GmbH) in 20 μl total volume using random hexamers. SYBR Green qPCR was performed using the SensiMix DNA Kit (Quantace Ltd.) and 0.2 μM of gene specific primers, synthesized by Sigm–Aldrich (electronic Supplementary Table S1). The first step of amplification was denaturation at 95°C for 2 min, followed by 40 cycles for 15 s at 95°C, for 10 s at 60°C and for 10 s at 72°C. Melting curve analysis allowed testing for specificity of each qRT-PCR. The PCR reactions were optimized for each primer pair by the use of serial diluted PCR products for standard curves with a dynamic range between 5 and 8 orders of magnitude and 80–110% efficiency at 60°C annealing temperature. All targeted genes were verified by sequencing of purified amplicons. All transcripts were quantified by triplicate measurements of 1 μl 1:5 diluted cDNA in 10 μl final reaction volume using either the StepOnePlus^TM Real-Time PCR System (Life Technologies) or the PikoReal^TM Real-Time PCR System (Thermo Scientific). Normalization of expression data was performed using three stably expressed reference genes. The RefFinder algorithm (Xie et al., 2012) was used to determine the most useful reference genes of the gene set: ACTB, B2M, HPRT1, LDHA1, and RPLP1. The relative gene expression was obtained by calculating the ΔΔCq values between the control and treatment animals.

For RNA sequencing eight animals per diet and sex were used (n = 8). In total 16 females and 16 males were included, one half being fed the 33% GM maize feed and one half being fed the 33% isogenic maize feed. TruSeq libraries were prepared using the Illumina “TruSeq Stranded Total RNA with Ribo-Zero Gold” library kit from 500 ng of intact whole RNA samples with RIN (RNA Integrity Number) values of at least 5.6 and OD260/280 of at least 1.8. All RNA samples from both male and female biological replicates were prepared and subsequently analyzed using the Agilent 2200 Tapestation. Ribosomal depletion was performed using library preparation kit. Sequencing was carried out on Illumina HiSeq 2500 instruments as a 2 × 125 bp run. RNAseq mass data have been submitted to the NIH resource via SRA submission¹ under the ID SUB2976951 Bioproject PRJNA400118 (accession numbers SAMN07556813 – SAMN07556876).

In the case of paraffin embedded tissues a standard rehydration protocol for 6 μm sections was applied as described earlier (Scholven et al., 2009). To reduce auto-fluorescence, sections were immersed in 0.1% Sudan Black diluted in 70% ethanol for 20 min at RT, followed by washing 5 min in 50% ethanol and rehydration in PBS (PBS including 1% BSA) 10 min at RT. Antigen retrieval was achieved by boiling in 10 mM citrate buffer, pH 6, for 40 min. Non-specific binding of the primary antibody was blocked by incubation in 10% goat serum in PBS (1% BSA in PBS) for 1 h at RT. Mucous and cell membranes were stained using an Alexa Fluor 594 conjugate of wheat germ agglutinin (Life Technologies) at a concentration of 5.0 μg/mL HBSS for 15 min at RT, followed by washing for 20 min in PBS. Rabbit monoclonal anti-NR1D1 antibody (Abcam, EPR10376, 0.086 mg/ml) was used at a 1:50 dilution (PBS with 1% BSA). A goat anti-rabbit IgG conjugated to Alexa Fluor 488 (pre-absorbed, Abcam) was used as the secondary antibody at a dilution of 1:400 for 1 h at RT. Incubations without the primary antibody, or with rabbit IgG1, were performed as negative controls. Nuclei were counterstained using 200 ng/ml DAPI (Roche) for 10 min, and sections were embedded in 50% glycerol in PBS.

Consistent with the RNA analysis, we analyzed protein samples of eight rats (same individuals) per study group and sex of study A (except for only seven males of the control group, due to limited tissue material). Intestinal tissues were homogenized in RIPA buffer (Cell Signaling Technology) supplemented with 1 mM PMSF and 1% (v/v) protease inhibitor cocktail (Calbiochem). Protein concentration was determined using a 2-D Quant Kit (GE Healthcare). Twenty micrograms of total protein were separated on a 12% (w/v) polyacrylamide gel under reducing conditions. Subsequently, proteins were transferred onto a nitrocellulose blotting membrane 0.45 μm (Sartorius) by semi-dry electroblotting for 70 min at 1.2 mA/cm². The membrane was blocked for 1h with 5% skim milk (Carl Roth) in TBST. Overnight incubation was performed at 4°C with a monoclonal anti-rabbit NR1D1 antibody (Abcam, EPR10376, 0.086 mg/ml) at a 1:1000 dilution or a HPRT monoclonal anti-rabbit antibody (Abcam) at a 1:7500 dilution in TBST 3% BSA. Donkey anti-rabbit IgG horseradish peroxidase-conjugated secondary antibody (Amersham Biosciences) was used at a dilution of 1:20000 for 2 h. Immunoreactive proteins were visualized using enhanced chemiluminescence (ECL Select, Amersham), and the obtained bands were visualized using the documentation system Fusion SL (Vilber Lourmat). Signals were quantified using the program Bio1D (Vilber Lourmat). For statistical analysis, the means (triplicate individual measurements) of treated samples and non-treated controls were calculated and were subjected to the Mann–Whitney-U-Test, a p < 0.05 being required for significance.

For RT-qPCR statistical analysis, the log two ratio of means (triplicate RT-qPCR measurements) of treated samples and non-treated controls were calculated and were subjected to the Kruskal–Wallis-test followed by the Mann–Whitney U-test as post hoc test, p < 0.05 being required for significance using IBM SPSS Statistics 20. Volcano plots as well as the comparison of NGS and RT-qPCR data (Spearman’s rank correlation) were performed by using Microsoft Excel 2010. To evaluate the RNAseq data, we applied Tophat (version 2.0) and Cufflinks (version 2.2) as described (Trapnell et al., 2012) with default parameters. Cufflinks also contains the programs Cuffdiff and CummeRbund to compute and visualize the differential expression between the groups. Cuffdiff computes variance estimates given multiple replicates and uses them to calculate the significance of the observed expression differences between the groups. Besides the raw computed p-value, Cuffdiff also outputs p-values corrected for multiple testing (FDR method (Benjamini and Hochberg, 1995)) and log2 fold changes, which we used to extract the differentially expressed genes for downstream analyses. Gene Ontology (GO) biological process analysis of genes with highest fold changes was performed using DAVID (Huang da et al., 2009). As the threshold for gene enrichment we applied an EASE Score (modified Fisher Exact P-Value) of 0.1. We ranked genes according to modified fisher exact p-values (the smaller the more enriched). In Figures 2–4 we showed 10 genes with lowest modified fisher exact p-value. To explore the biological context of differentially expressed genes, we performed the GO classification, making use of the category BP_Fat (biological process). Network analysis was applied to predict protein-protein interaction. We imported genes with highest fold changes to STRING database and visualized the networks by Cytoscape software (Shannon et al., 2003).

Gene expression studies addressing the physiologically important pathways of apoptosis, DNA damage and repair and unfolded protein response (UPR) have been introduced in order to evaluate if GM feed has an impact on fundamental signaling pathways associated with cellular stress response and intestinal dysfunction. In the present study, samples were taken from two 90-day studies (the so-called trials A and B; (Zeljenkova et al., 2014)) and 79 marker genes were quantified using qPCR. These marker genes split into 26 genes whose products are involved in apoptosis and UPR pathways and 27 genes whose products are involved in DNA damage and repair (DDR) pathways. Intestinal tissue samples from ileum (Figure 1A) as well as colon (Figure 1B) of the 33% GMO maize and 33% near-isogenic maize control group of study A were analyzed. We applied the generally accepted threshold of twofold difference for considering qPCR results for differential expression. In both intestinal regions, we observed statistically significant fold changes between the groups. However, those marker genes did not exceed a twofold up- or down-regulation cut-off threshold (Figures 1A,B). This applied to 6 genes in ileum of males, 26 genes in female ileum, 11 genes in male colon, and 4 genes in female colon. Individual profiles of marker genes and statistical significances are shown in Supplementary Figures S1, S2 (electronic Supplementary Material). In the case of the apoptosis pathway, the profiles demonstrate that essentially no distinct regulatory pattern was observed. One pro-apoptotic gene (Fas) exceeded the cut-off point by a fold change of 2.17 in the female ileal group and a fold change of 2.25 in the male colon group without a statistical significance in both cases. For some marker genes, we observed significant tendencies (less than twofold difference) toward up- or downregulation. However, these genes were well below the twofold cut-off and did not point to co-ordinately differential expression pattern of pro- and anti-apoptotic genes. The analysis of DDR marker genes did not reveal any regulation with respect to the fold change cut-off, while 13 genes showed statistically significant tendencies, which again did not reveal informative patterns of differential expression (Supplementary Figures S1, S2). The UPR pathway analysis revealed no differential expression with respect to the applied fold change cut-off. In female ileum 54% of the studied UPR marker genes showed statistically significant 1.2 – 1.82-fold differences between the study groups, whereby seven of these genes encode heat shock proteins (Figure 1A and electronic Supplementary Material 1, Supplementary Figure S1F). However, in male animals as well as in colon tissue no such tendencies were observed (Figure 1B and electronic Supplementary Material 1, Supplementary Figures S1E, S2E,F). Taken together, we conclude that only minor tendencies, which do not reflect distinct perturbations of the studied pathways, were observed.

To assess if feeding GM-maize affects global gene expression in the rat small intestine, we compared RNAseq based expression profiles of rat ileum from both study groups receiving 33% GM or near-isogenic control maize. In total, we detected 34.276 (study A) and 25.570 (study B) transcripts by the RNAseq profiling. In study A, we identified 102 transcripts to be at least twofold up-or down regulated in GM maize-fed male animals, while in female intestinal tissues, 63 transcripts showed altered concentrations comparing both study groups (Log2 fold-change > 1 or < -1, q-value < 0.05, false discovery rate [FDR] < 5%; see electronic Supplementary Material 1, Supplementary Figures S3A,C). In study B, we observed a differential expression of 30 and 18 transcripts in males and females, respectively (Log2 fold-change > 1 or < -1, q-value < 0.05, false discovery rate [FDR] < 5%; see electronic Supplementary Material 1, Supplementary Figures S3B,C). In a comparative analysis, we identified 18 commonly differentially expressed transcripts in males after matching both studies A and B, while 12 transcripts overlapped in datasets of differentially expressed genes comparing females of both studies A and B (electronic Supplementary Material 1, Supplementary Figure S3C).

In both studies A and B 0.3–0.8% of all transcripts detected showed significant differences between the GMO and control groups (Supplementary Figures S3A,B). To determine its biological meaning, we performed an enrichment analysis of GO annotated biological processes using DAVID (Huang da et al., 2009) and applied a network analysis using STRING (Franceschini et al., 2013) and Cytoscape (Shannon et al., 2003). In study A, it was apparent, that the list of genes with the highest fold changes included various members of the circadian clock pathway (Figures 2A,B). For example, in males of study A, all transcripts being more than sixfold down regulated in the GMO group (Per3, Nr1d1, Tef, and Dbp) are connected with the pathway of circadian rhythms, while genes showing the highest up-regulation point to the relevance of metabolic processes as revealed by GO-analysis. For both male and female animals of study A, the top 50 transcripts with a more than twofold difference between the study groups are shown in Figures 2A,B. Throughout both animal trials and both sexes, the top ten biological processes associated with differentially expressed genes (Log2 fold-change > 1 or < -1, q-value ≤ 0.05) pointed to circadian rhythms, and in all cases, the highest significance in the enrichment of genes was observed for “rhythmic process” (Figures 3, 4C,D). Consequently, further biological processes included metabolic processes and nutrient transport, which are well known to be strongly affected by circadian rhythms (Green et al., 2008; Hussain and Pan, 2009; Hussain, 2014). In line with the potential functional interactivity between transcripts that showed differential expression between study groups, we performed a network analysis, which revealed connected elements that are known to interact (based on Text mining, Experiments and Databases). The interaction network nodes with most direct interactions always pointed to the circadian clock pathway, as shown in Figure 2 for study A males (e) and females (f), as well as study B males and females (Figures 3E,F). In Figure 4, the commonly differentially expressed transcripts comparing studies A and B are shown, thereby underlining the consistency of observed perturbations in circadian pathways in both animal trials.

In order to confirm our hypothesis that differences in gene expression between study groups in both studies A and B are related to the time of sampling, we performed RT-qPCR analysis of seven genes centrally involved in circadian rhythms. First, we technically compared fold differences observed using RNAseq and RT-qPCR analyses and observed a good correlation of differential expression values (Spearman’s rank correlation coefficient: 0.94), while the RT-qPCR approach tended to generate higher total fold differences (electronic Supplementary Material 1, Supplementary Figure S4). Using the RT-qPCR approach, we observed strong up regulation of the circadian clock gene Bmal1 (Arntl) as well as a strong down regulation of Dbp, Nr1d1 (rev-Erbα), and Per3 in male and female ileal tissues and in both studies (Figures 5A,B). Likewise, we detected similar effects in the male and female colonic tissues (electronic Supplementary Material 1, Supplementary Figure S5). Here, the mRNAs of bmal, dbp, nr1d1, and per3 showed the same regulation as in the ileum, and cry1 was in addition significantly upregulated. Nr1d1 exists in two slightly different isoforms: rev-Erbα1 and rev-Erbα2 (Rambaud et al., 2009). A detailed analysis of these two isoforms in the female ileum via RT-qPCR revealed a comparable circadian regulation. However, we observed a notable difference in the expression of the isoforms. Rev-Erbα1 showed a 40–60-fold higher expression compared to rev-Erbα2 (data not shown). Quantification of Nr1d1 protein in female ileal samples resulted in a statistically non-significant 21% down-regulation between the GMO and the isogenic control group (p = 0.117) (Figure 5C). Furthermore, Nr1d1 expression was localized in the nuclei of small intestinal epithelial cells as assessed by immunostaining of paraffin sections of small intestine tissue samples (Figure 5D). Accordingly, a tendency toward a lower protein expression was observed in the GMO groups.

We compared our data sets with a gene set from TissGDB, which is tissue-specific gene annotation database in cancer² (Kim et al., 2017). The RNAseq based expression profiles from both studies A and B included 65 of the 80 reference genes for colorectal cancer. In both male as well as in female animals, there was a near-complete absence of significantly altered gene expression levels in the GMO-group compared to controls (Log2 fold-change > 1 or < -1, q-value < 0.05, false discovery rate [FDR] < 5%; see electronic Supplementary Material 1, Supplementary Figure S6). However, there was an exception in male rats in study A. Here, we observed a significantly 2.8- and 4.3-fold higher expression of Cyp2s1 and Agp8 in the GMO-group compared to controls, respectively.