**3. Discussion**

DNA methylation can provide additional heritable information beyond that of the DNA sequence in plant genomes [2]. Tissue culture is considered a stressful environment and thus trigger epigenetic changes in plants [34]. Culture-induced DNA methylation has been found in different species, including rice, maize, and barley [35–37]. Furthermore, many studies have shown that regeneration under various selection stresses or from various donor tissues induced changes in methylation patterns [38–40]. When uniform callus donor tissue was used in an *Agrobacterium*-mediated transformation procedure in rice [41,42], there was no difference in methylation among donor tissues between experiments. It has previously been reported that tissue culture reduces mC in rice, and this reduction in mC is stable from T2 to T6 generation of regenerated plants [18]. Our sequence data confirmed that the tissue culture process leads to a reduction in global DNA methylation in CK and that this reduction was maintained in the *pOsNAR2.1-OsNAR2.1* line at least until the T4 generation (Figure 1a). Moreover, we showed that the global methylation level returned to a level similar to that of WYJ7 after removal of the transgene by a segregation process in the SDWT line. Since Stroud et al. [18] confirmed that the loss of methylation in regenerated plants is stable across generations, our results suggest that the global methylation status in the SDWT is more unstable and is unable to maintain the mC decrease across the generations. However, the insertion of the transgene stabilized the massive loss of mC in *pOsNAR2.1-OsNAR2.1*. The tissue culture process can cause the loss of methylation in both WYJ7 and NP (Figure 1a) according to both our results and those of Stroud, et al. [20], but *pOsNAR2.1-OsNAR2.1* in WYJ7 showed stabilization of the loss rather than *pUbi-OsNAR2.1* in NP. Therefore, we consider that overexpression of *OsNAR2.1* is not the reason behind the stabilization. We prefer the explanation that the unstabilized methylation of SDWT was caused by the double genetic change, gaining and losing the transgene during generation. Furthermore, the way insertion of the transgene alters global methylation appears to depend on the rice genetic background, since global methylation reduction was observed in the *pOsNAR2.1-OsNAR2.1* WYJ7 line whereas it was increased in the *pUbi-NAR2.1* and RNAi NP lines.

We found global methylation decreased in the WYJ7 *OsNAR2.1-*overexpression line, while increased *OsNAR2.1*-overexpression, and RNAi lines in the NP. We suggest that rice varieties have different sensitivities to DNA methylation in the transgenic process. It has been reported that there are extensive variations in DNA methylation among plant inbred lines, and that DNA methylation can provide unique information in explaining variation of phenotype in maize [43]. Vilperte et al. [5] reported that the methylation status of genes showed significant differences in four different maize backgrounds with the same transgene. In our results, although in both the *pOsNAR2.1-OsNAR2.1* and *pUbi-OsNAR2.1* lines, plant height and yield of per plant were significantly higher than in the wild-type, but the plant height increased 4.69% and 17.36%, and yield of per plant increased 49.07% and 34.97%, respectively, in WYJ7 and NP (Table 1, Supplementary Tables S2 and S3). Even overexpression of the same gene was able to cause different phenotypes in different backgrounds. This result could be caused by several possible factors. The first of these is the original traits of the two wild types: WYJ7 is photo-sensitive late-maturing japonica rice and NP is photo-insensitive early-maturing japonica rice [44,45]. The two wild-types are, therefore, affected differently by the circadian clock. The circadian clock regulates NO3 - uptake and usage, and thus the expression of *OsNAR2.1* [30,46]. A second factor is the promoter of the transgene. It is well known that different promoters have different effects on transgenes [47]. In our research, the expression of *OsNAR2.1* increased around 2.5 times with the *OsNAR2.1* native promoter and increased four times with the *Ubi* promoters. Wang et al. [48] reported that overexpressed auxin-inducible gene (*ARGOS*) increased plant height in *Arabidopsis* under the 35*s* promoter but showed no phenotypic change under the *Ubi* promoter. The third factor is background. There have been many reports of the same transgene showing different phenotypes in plants from a variety of genetic backgrounds. It has been reported that, while the *DREB3* transgene was detected in wheats from four different genetic backgrounds, only three lines expressed the transgene, and only two showed phenotypes of higher yield [24]. Knockout of *OsNramp5*, a member of the natural

resistance-associated macrophage protein (NRAMP) family, decreased yield in Xidao 1, a japonica rice cultivar, but did not alter yield in indica hybrid rice [26,49]. Even knockout of the same gene in different backgrounds with same genome also resulted in phenotypic variation, as shown by the findings of Yang et al [50]: knockout of *OsNramp5* using CRISPR/Cas9 in two japonica varieties, Nanjing 46 (NJ46) and Huaidao 5 (HD5), resulted in similar plant height, grain number, and seed setting rate, but with increased panicle number in NJ46 but not in HD5. Researchers surmised that damage to the recipient genome caused by insertion fragments, efficiency of the transgene promoter, metabolic imbalances and environmental constraints could be the reasons behind undesirable phenotypes of transgenic plants [24,27,51], while we suspected that the insertion of the transgene causing different methylation change in different backgrounds could be another reason.

The Venn diagram of genes in the methylation region in WYJ7, SDWT, CK, and *pOsNAR2.1-OsNAR2.1* showed that all three had over 1000 unique methylation genes compared with WYJ7 (Figure 1f), and all had 200–400 unique methylation genes compared to each other. The results suggest that, compared with WYJ7, all three samples (CK, SDWT, and *pOsNAR2.1-OsNAR2.1*) have abundant genes with different degrees of methylation, but not all three lines showed phenotypic changes (Figure 4a). In spite of having these differentially methylated genes, SDWT still showed a similar phenotype to WYJ7. These results suggest that the global methylation status change had a stronger influence on plants than changes in the methylation of individual genes.

Our results showed that the regenerated line CK had more hypo-DMRs than hyper-DMRs and that all three transgenic lines had more hyper- than hypo-DMRs, compared to the wild-type. However, in SDWT, the percentage of hyper- and hypo-DMRs are similar, 51% and 49%, respectively. The total number of DMRs in SDWT (2732 DMRs) was much less than that in CK (5122 DMRs) and the three transgenic lines (3539 DMRs in *pOsNAR2.1-OsNAR2.1*, 5232 DMRs in *pUbi-OsNAR2.1* and 4239 DMRs in RNAi) (Figure 2a–c,f,g). Stroud et al. [18] reported that all of their 12 regeneration lines had much more hypo-DMRs than hyper-DMRs, confirming our results. Therefore, it is possible that the reprogramming process could cause an abundance of hypo-DMRs, while insertion of the gene could cause hyper-DMRs instead, regardless of the function of the inserted gene. Furthermore, for the SDWT lines, WT went through both reprogramming and transgene processes but after segregation, it was found that this double genetic change resulted in nearly 3000 DMRs with the numbers of hyper- and hypo-DMRs returning to the baseline level of WYJ7. Since the phenotype of SDWT did not show significant differences compared to WYJ7, this indicated that the balance of hyper- and hypo-DMRs could be an important factor for the stability of the epigenetic status.

We found that the methylation status in the CK line regenerated from the tissue culture process gave rise to a more variable DMR methylation level clustering pattern than in the CK replicates. This is similar to the observations of Hsu et al. [9]. These authors reported that methylomes in different stages of callus showed a high level of variability and the global methylome demonstrated a prominent genetic/cultivar-related impact. Kaeppler et al. [4] showed that epigenetic modifications of genomic DNA are less stable in culture. Interestingly, both the transgenic line (*pOsNAR2.1-OsNAR2.1*) and the wild-type lines derived from the segregation process of *pOsNAR2.1-OsNAR2.1* (SDWT) showed a more consistent DMR methylation level clustering pattern. This suggests that the transgenic process may play a role in stabilizing the global methylation status of tissue cultured plants.

However, it is possible that the CK line has a genetic change in the genome in addition to epigenetic changes. For example, there have been many reports of the occurrence of transposons, such as *Tos17,* in regenerants [34,52]. To clarify the effects of epigenetic changes on the CK phenotype, we have listed the methylation sequence data and agronomic traits of each CK line from three independent calli in Supplementary Table S4. We found that, despite the differences between individuals, there is an overall consistency in the data of CK replicates (Supplementary Table S4). The Table shows that normalized mC number, mCHG ratio of epigenetic and plant height, panicle length, grain number per panicle, and yield of phenotype are all significantly lower in CK. Even if we cannot exclude possible genetic change in CK, the changes, at least, did not influence the phenotypes of the individual CK replicates

in our study. It has also been reported that genetic variation during tissue culture is related to the CCGG target, which suggests that genetic and epigenetic changes in regenerants are relevant [34,53,54]. These changes could be the reason for the clustering of variable DMR methylation levels in CK but not causing phenotypic changes.

Previous research has reported that phenotypes are influenced by both genetic and epigenetic mechanisms [20,55,56]. Epigenetic modifications caused by the tissue culture process occur in an apparently random manner and it is thus difficult to predict phenotypic changes resulting from these modifications [4,5,7,57,58]. Here, we summarized global methylation, genetic and phenotype changes of five plant materials from two varieties, as shown in Table 2. The results indicate that, in the transgenic lines (*pOsNAR2.1-OsNAR2.1*, *pUbi-OsNAR2.1,* and RNAi), the phenotype is mainly determined by the nature and function of the transgene, but the global methylation change in transgenic lines is determined by the rice background rather than the function of the transgene, which may cause the difference between the phenotypic changes in the two overexpression lines. However, the CK lines showed fewer phenotypic changes, with the methylation status being significantly decrease from WYJ7 as result of the tissue culture process and showing DMR variation between replicates (Figure 3).
