Induced changes to DNA methylation and associated perturbations to gene expression has been reported for genes associated with organogenesis (De-La-Peña et al., 2012) and other developmental processes in plants (Nic-Can et al., 2013). In comparison, *in vitro* regeneration protocols that preserve meristem anatomy and function, and which generate new plants from the activation of previously dormant meristems are known as micropropagation systems and are widely viewed as being genetically analogous to field cuttings. These systems are thought to generate daughter ramet plants that are morphologically and genetically faithful replicates of the original explant material (Kahn, 2012). To date, however, little is known about the degree to which the epigenetic profiles (and so associated cell regulatory processes) of regenerated plants from micropropagation represent faithful replicates of the original mother plant. In this study we combine methylation-sensitive amplified polymorphism (MSAP) and methylation-sensitive genotype by sequencing (msGBS; Xia et al., 2014) to assess the epigenetic fidelity of meristem micropropagation and to seek specific methylation signatures associated to *in vitro* propagation in cassava (*M. esculenta*).

## Materials and Methods

### Plant Material

Five varieties of cassava (*M. esculenta* Crantz) namely Kiroba, Kizimbani, Kibandameno, Mfaransa, and Mzungu were used in this study. *In vitro* micropropagated samples were obtained from the tissue culture facility of the Mikocheni Agriculture Research Institute (MARI) Dar es Salaam, Tanzania and were propagated as described by (Konan et al., 1997). Field samples of the same genotypes were grown at the Sugarcane Research Institute-Kibaha (SRI-KIBAHA), Tanzania. Samples were collected from young leaves (last leaf emerged from bud), newly mature leaves (first fully expanded leaf) and primary root tips from three individual plants representing each cassava variety both from field cutting and meristem culture samples. All samples were kept in dry ice in the field and stored at −80◦C until required for DNA extraction.

### DNA Isolation

DNA was extracted from all samples at MARI using a DNeasy plant mini kit (Qiagen) according to the manufacturers' instructions. DNA concentration and quality was estimated using a Nano-Drop 1000 Spectrophotometer (Thermo Scientific). DNA was lyophilised prior to transport to the Plant Research Centre in Adelaide, Australia for use in subsequent MSAP or msGBS analyses. Upon arrival, all DNA samples were re-suspended in nuclease free water (Sigma), and re-quantified using the Thermo Scientific NanoDropTM 1000 Spectrophotometer. DNA concentrations were standardized to produce working solutions of 10 or 20 ng/ul.

## Methylation-Sensitive Amplification Polymorphism Procedure

A modification of the MSAP technique (Reyna-Lopez et al., 1997) was used as described by Rodríguez López et al. (2012). In brief, genomic DNA was digested with a combination of the methylation insensitive restriction enzyme *EcoR*I and one of two isoschizomer enzymes that exhibit differential sensitivity to DNA methylation (*Hpa*II and *Msp*I; **Table 1**). Adapters were ligated to the digested gDNA and then used as template for the first of two consecutive selective PCR amplifications in which the primers were complementary to the adaptors but possessed unique 3 overhangs (**Table 1**). *Hpa*II/*Msp*I selective primers were end labeled using a 6-FAM reporter molecule for fragment detection using capillary electrophoresis. A total of six primer combinations (**Table 1**) were tested in a pilot study using eight randomly selected DNA samples.

## Sample Fractionation by Capillary Electrophoresis

Single base resolution separation of the MSAP products was achieved by capillary electrophoresis on an ABI PRISM 3130 (Applied Biosystems, Foster City, CA, USA) housed at the Australian Genome Research Facility Ltd, Adelaide South Australia. Sample fractionation was performed as follows: 2 μl of the labeled MSAP products were combined with 15 μl of HiDi formamide (Applied Biosystems, Foster City, CA, USA) and mixed with 0.5 μl of GeneScanTM 500 ROXTM Size Standard (Applied Biosystems, Foster City, CA). Samples were heatdenatured at 95◦C for 5 min and snap-cooled on ice for 5 min. Samples were fractionated at 15 kV for 6 s and at 15 kV for 33 min at 66◦C.

### Methylation Sensitive Genotyping by Sequencing

We performed the methylation-sensitive modification of the genotype by sequencing (GBS) technique (Poland et al., 2012) as described by Xia et al. (2014). In brief, a two-enzyme approach was used to generate restriction products. In this experiment, only one enzyme combination was used *(Msp*I with *EcoR*I). The

∗*Indicates the selective primers used to analyze all samples.*

selected enzyme combination was based on the results obtained using the MSAP approach. Two hundred nanogram of genomic DNA from each of the 86 selected samples [comprising three replicate per tissue/variety and growing condition (i.e., *in vitro* or field) see Supplementary Table S1] were used in a reaction volume of 20 μl containing 2 μl of NEB Smart cut buffer, 8 U of HF-*EcoR*I (High-Fidelity) and 8 U of *Msp*I (New England BioLabs Inc., Ipswich, MA, USA). Reactions were prepared in a 96 well plate containing 87 reactions (86 DNA samples plus one Negative control water sample) and conducted on a BioRad 100 thermocycler at 37◦C for 2 h and then 65◦C for 20 min for enzyme inactivation. A set of 96 barcoded adapters with an *Msp*I overhang and a common Y adapter with an *EcoR*I overhang were designed for the ligation reaction using barcode script made by Thomas P. van Gurp1 . Adapters were annealed prior to ligation as described by Poland et al. (2012). A full list of adapters for *Msp*I (with corresponding barcodes and cassava samples) and *EcoR*I is listed in Supplementary Table S1. The ligation reaction (40 μl) was carried out on the same PCR plate adding to the restriction products T4 Ligase (200 U) and T4 Ligase buffer (NEB T4 DNA Ligase #M0202), 0.1 pmol of the respective barcoded *Msp*I adapter and 15 pmol of the common Y-adapter. Ligation was completed at 22◦C for 2 h followed by an enzyme inactivation step of 20 min at 65◦C. Five micro liter from each ligation reaction were pooled into a single tube and then divided into two equal volumes for column clean-up using PureLinkR PCR Purification Kit (Life Technologies) according to manufacturer's instructions. Samples were re-suspended in 60 μl of nanopure water. Both clean-ups were then combined and divided again into eight samples for PCR amplification. Each 25 μl PCR consisted of 10 μl of DNA digested/ligated library), 5 μl of 5x NEB MasterMix, 2 μl of 10 uM Forward and Reverse primers at 10 uM (Supplementary Table S1). Reactions were performed in a BioRad T100 thermocycler for eight cycles consisting of 95◦C (30 s), 62◦C (30 s), 68 ◦C (30 s). All eight PCR products were pooled and then purified first using a PureLinkR PCR Purification Kit (Life Technologies) according to manufacturer's instructions (resuspended in 30 μl). Excess adaptor was finally removed using Ampure XP magnetic beads (Beckman) by mixing 30 μl of the pooled PCRs with 22.5 μl of beads. Captured fragments were eluted in 30 μl of water. Next, 125 bp paired-end sequencing was performed in one Illumina HiSeq 2000 v4lane (Illumina Inc., San Diego, CA, USA) by QBI Centre for Brain Genomics.

### Statistical Analysis

Analysis of Genetic/Epigenetic Variability using MSAP

The MSAP technique uses MspI or HpaII as isoschizomers; both can cleave the motif CCGG in the absence of methylation. MspI can also cleave hemi-methylated dsDNA (mC in one DNA strand only) or fully methylated DNA sequences where the internal cytosine is methylated CmCGG. However, it cannot digest hemi-methylated and fully methylated at the external cytosine site, viz: mCCGG and mCmCGG motifs (Walder

1www.deenabio.com/gbs-adapters

et al., 1983; Reyna-Lopez et al., 1997). In contrast, HpaII is more sensitive to methylation but can cleave hemimethylated DNA at the external cytosine position (mCCGG; Mann and Smith, 1977; Reyna-Lopez et al., 1997). Direct comparison of MspI profiles with those generated by the more methylationsensitive HpaII therefore does not provide a definitive contrast between genetic variation and that attributable to changes in methylation (Fulnecek and Kova ˇ ˇrík, 2014). For these reasons, simple comparisons were made between profiles generated from various tissues of plants grown in the two settings (micropropagation and field cuttings) under the reasonable assumption that consistent differences will arise from differential methylation (driven by RdDM) rather than by repeated chance mutations.

MSAP profiles were visualized using GeneMapper Software v4 (Applied Biosystems, Foster City, CA, USA). Two matrices containing allelic information were generated. First, a qualitative analysis was carried out in which epiloci were scored as "present" (1) or "absent" (0) to form a presence/absence binary matrix. In this case, the selection of MSAP fragments was limited to allelic sizes between 100 and 580 bp to reduce the potential impact of size homoplasy (Caballero et al., 2008). Profile polymorphisms between DNA samples from the same cassava variety but extracted from different tissues (young leaves, newly mature leaves, and primary root tips) were retained as inter-tissue methylation differences. Polymorphisms between DNA samples from *in vitro* culture plants and from field grown plants were considered as *in vitro* culture induced methylation differences. Second, a matrix containing peak heights of fragments with allelic size between 50 and 550 bp was created for quantitative analysis (Rodríguez López et al., 2012). In both cases, different levels of hierarchy were generated to group the samples. Samples were first grouped according to cassava variety. Then, samples were divided into field grown and *in vitro* grown. Finally, samples were separated into the three different tissues of origin (young leave, mature leave, and roots).

For the analysis of the MSAP qualitative data, GenAlex v6.4 software (Peakall and Smouse, 2006) was used to infer pairwise epigenetic PhiPT distances (estimation of genetic/epigenetic distances) between different cassava samples. Analysis of molecular variance (AMOVA) was then performed using the same software to test the significance of the estimated PhiPT between tissues (Michalakis and Excoffier, 1996). An allele frequency table was generated using GenAlex 6.4 to find *in vitro*/field specific qualitative markers for each cultivar and for all cultivars. Finally, the visualization of the patterns of tissue epigenetic variations in this study was done by constructing a Principal Coordinates Analysis (PCoA).

For each variety, the peak height intensities of the epiloci generated using MSAP were compared between field grown and *in vitro* tissue samples as well as compared between samples from different tissue origins within field and *in vitro* groups. Preceding comparative analysis the data was filtered by removing epicloci containing excessively low peak height intensities across the complete set of samples. From this reduced set of epiloci the peak height libraries were normalized using the model based weighted trimmed mean method derived in Robinson and Oshlack (2010). For each pair of tissue groups being investigated, the normalized peak heights were extracted and compared using the approach described in Robinson and Smyth (2007, 2008). This approach initially assumes the normalized peak heights are distributed as a negative binomial with a common dispersion calculated across the complete set of epiloci for the two groups. From this, individual epiloci dispersions were calculated using the empirical Bayes methods of Robinson and Smyth (2007). An exact statistical test was then conducted for each epiloci to determine differences in peak heights between the two groups (Robinson and Smyth, 2008). The *p*-values obtained from these tests were then appropriately adjusted for multiple comparisons using the false discovery rate (FDR) method of Benjamini and Hochberg (1995). All analyses were performed using the differential expression analysis R package edgeR (Robinson et al., 2010) available in the R statistical computing environment (R Development Core Team Foundation, 2015).

### Analysis of Genetic/Epigenetic Variability using GBS Data

For the processing of Illumina HiSeq 2000 v4 data, the sequences from the unfiltered fastq Illumina output were separated into samples using the barcode sequence and trimmed to 64 bp using the software TASSEL (sourceforge.net/projects/tassel/). Only sequences with one of the exact used barcodes followed by the expected sequence of three nucleotides remaining from an *Msp*I cut-site (5 -CGG-3 ) were retained for analysis. Sequences present in the negative water control were also removed from the analysis. Finally only sequences present in three or more different samples were kept for analysis. A matrix of sequence abundance of was then generated for further analysis.

Using the differential expression analysis procedure outlined in Section "Analysis of Genetic/Epigenetic Variability Using MSAP," the sequence abundances were compared between field grown and *in vitro* tissue samples and also compared between samples from different tissue origins within field and *in vitro* groups. Sequences presenting significantly different number of reads between all *in vitro* and all field grown samples for each variety were isolated. Finally, only those sequences that presented the same variation (increase or reduction of number of reads in all varieties when comparing *in vitro* against field grown samples) in at least four of the five studied varieties were considered micropropagation induced markers. Due to the extremely low probability of a mutational event leading to the generation of these markers happening in all plants from all varieties during culture we can consider that such markers are differentially methylated regions (DMRs) induced by the micropropagation procedure. Detected DMRs were then selected for blast analysis against the cassava (BLASTN, nucleotide query to cassava nucleotide database genome blast tool in Phytozome. Top hits indicating differential methylation of a genic region were sought by comparing exons, introns, and flanking sequences (5 kb upstream of the Transcription Start Site and 5 kb downstream of the Transcription Termination Site).

## Results

### Analysis of Genetic/Epigenetic Variability using MSAP

### Estimation of Genetic/Epigenetic Differences Based on Qualitative Analysis

Methylation-sensitive amplified polymorphism profiles generated a total of 164 loci (13 unique to *Hpa*II, 22 unique to *Msp*I, and 129 common to both enzymes) for the 86 samples of five cassava cultivars used in this study. PCoA analysis created from a simple presence/absence similarity matrix of the combined MSAP profiles revealed clear separation between *in vitro* propagated cultivars and their field counterparts for all cultivars (**Figure 1**). Calculated genetic/epigenetic distances between field and *in vitro* samples were significant for all cultivars (**Table 2**). In general, calculated distances between *in vitro* and field samples where higher for samples restricted using *Msp*I (**Table 2**). All pairwise PhiPT between *in vitro*-grown samples and those grown in the field tissues were significant for all varieties using both *Hpa*II and *Msp*I (**Table 2**). In general, genetic/epigenetic distances were higher between tissues recovered from plants grown in the field than between those taken from *in vitro-*grown plants (**Table 3**). The reduced divergence between tissues taken from *in vitro*grown material was also evident from the PCoA analysis, with samples from different tissues of *in vitro* propagated plants occupying less eigen space than those of the same tissues obtained from field grown plants (**Figure 1**). However, the level of variability observed within tissue types did not differ significantly between *in vitro-*grown and field-grown material.

We further analyzed the differences existing between each tissue derived from the field grown plants and all samples from plants grown *in vitro.* The aim here was to investigate which of the *in vitro*-grown tissues generated MSAP profiles were most similar to the field samples. Distance estimates were significant for all pairs, but it was consistently smaller between young leaves from field grown plants and bulked *in vitro* tissues (**Table 4**).

### Estimation of Genetic/Epigenetic Differences Based on Quantitative Analysis

We selected 106 markers for quantitative analysis of MSAP profiles based on peak height data. In general, both enzymes, yielded more markers separating between tissues from fieldgrown plants than those taken from *in vitro* material (62 vs. 44 for *Msp*I and 44 vs. 42 for *Hpa*II; **Table 5**; For a list of all fragments and their levels of significance see Supplementary Table S2 for *Msp*I and Supplementary Table S3 for *Hpa*II). However the number and scale of these differences varied between cultivars.

A total of 14 and 15 markers for *Msp*I and *Hpa*I respectively were found to be significantly different (*p* < 0.005) between all *in vitro*- and field-derived material. Two of these markers generated using *Hpa*II, epiloci 55 and 101 bp, were able to diagnose *in vitro* from field samples of three varieties, Kiroba, Kibandameno, and Kizimbani (**Figure 2**). The epilocus 55 bp generated using *Msp*I

was also significantly different between propagation systems for the same varieties but not epilocus 101 bp (**Figure 2**).

### Analysis of Epigenetic Variability using GBS Data

In total, we generated 236,624,193 raw reads from the HiSeq 2000 v4 lane, of which 71,723,843 (32.3%) passed quality filter and contained the expected exact matches to sequences of the barcode adapter, *Msp*I restriction product site and the *EcoR*I adapter, and which appeared in at least three biological samples but were absent from the negative (water) control. On average, 754,980 high quality reads were produced per DNA sample. Collectively, this included 357,271 unique sequence tags across all samples. The number of these markers that differed significantly in abundance between *in vitro*-grown and field-grown samples varied considerably between varieties: 3,298 for Kiroba; 25,683 for Kizimbani; 2,029 for Mfaransa; 34,098 for Mzungu and 17,702 for Kibandameno. Most of these sequences were more abundant in the field-grown samples (**Table 6**) and the overwhelming majority was variety-specific responses. We next sought to identify candidate generic epimarks that differentiate between propagation systems across all varieties. When the most stringent filter for differential abundance was applied (i.e., reads with an FDR lower than 0.05; the phase of differential abundance being conserved across all genotypes and tissues, and the absence of variety-specific SNPs) just 105 (0.03%) of unique differential sequences featured in the profiles of all varieties and showed a common pattern of phasing (Supplementary Table S4). There was also a marked difference in the phase of these marks, with just four tags being more abundant among *in vitro-*grown samples compared with 101 that were more numerous in the field-grown plants.

When compared against the cassava nucleotide database genome using the blast tool in Phytozome2 , 102 differentially methylated sequences generated one or multiple hits against the cassava genome. Eighty nine top hits were associated to a gene (i.e., mapped within a window of 5 kb from the gene). BLAST results indicate that the homologs of these genes are involved in many processes, including cell differentiation, plant development, sugar metabolism, nucleic acid methylation, stress response, photosynthesis, signaling and transposon activation (Supplementary Table S5). Of the 89 differentially methylated genes 45 have been previously mentioned in the literature as having homologous genes that are: regulated by DNA methylation (14) or other epigenetic mechanisms (10), methylated or differentially methylated under different growing conditions (4), implicated in the regulation of DNA methylation (9) or other epigenetic mechanisms (1) and DNA binding proteins affected by methylation of their target sequences (3).

## Discussion

Micropropagation via nodal cuttings relies on the regeneration of pre-existing meristems and so is widely considered to be

2http://www.phytozome.net

### Kitimu et al. Micropropagation-induced epigenetic changes in cassava

### TABLE 2 | Summary of calculated genetic/epigenetic distance (PhiPT) between *in vitro-* and field-grown samples.

*Pairwise genetic/epigenetic distance between samples grown in vitro and those grown in the field from each variety were calculated using MSAP profiles obtained by restricting gDNA from five cultivars of cassava with EcoRI and HpaII or MspI using GenAlex 6.5.1 AMOVA. p-values were calculated based on 1000 permutations.*

TABLE 3 | Effect of *in vitro* culture on epigenetic differentiation between tissues.

*Calculated pairwise PhiPT values (epigenetic distances) between tissues from field cuttings and in vitro propagated lines. Distances were calculated using GenAlex 6.5.1 AMOVA on MSAP profiles obtained by restricting genomic DNA with MspI and HpaII. Ylv, young leaves, Mlv, mature leaves, Rt, roots.*

analogous to field cuttings because they do not pass through a state of disorganized (dedifferentiated) tissue state (Dale and McPartlan, 1992). Nevertheless, ramets recovered from nodal micropropagation can still exhibit signs of increased morphological variability (somaclonal variation) when compared to those recovered from field cuttings (Debnath, 2005). Us-Camas et al. (2014) suggested that such observations might be explained by the stressful environment experienced by *in vitro*grown plants (i.e., high relative humidity, low ventilation rate, high concentrations of sugars and plant growth regulators, and low light availability). Under these conditions, cultured plants cells are also forced to change their molecular make ups in order to generate different cell types. Cell division to generate tissues and organs require a precise coordination of genetic and epigenetic processes (Miguel and Marum, 2011; Smulders and de Klerk, 2011). For micropropagation systems that rely on dedifferentiation and *de novo* organogenesis, *in vitro* culture can often yield occasional regenerants that are phenotypically off-type (e.g., Lakshmanan and Taji, 2000; Rout et al., 2000; Da Silva et al., 2015). In contrast, those recovered from meristem micropropagation are widely reported to remain more faithful to the phenotype of the parental plant in range of species (e.g., Villordon and LaBonte, 1996) including cassava (e.g., Santana et al., 2009). There is nevertheless a large body of evidence indicating that changed growing conditions often induces moderations in global methylation patterns in culture (for review see Pastor et al., 2013) and this leads to the plausible expectation of epigenetic divergence between plants cloned by meristem micropropagation and field cuttings. Evidence supporting this assertion came from a study by Baranek et al. (2010), who used MSAP profiles to compare daughter plants recovered from field cuttings and micropropagated nodal segments of two grape vine varieties. The authors found consistent differences between the two systems in their clustering on resultant dendrogams. However, the work failed to further characterize the variation in terms of tissue type or to provide sequence identity for the differential epimarks. Characterizing such epigenetic differences may prove useful not only for the mere detection of putative somaclonal variants (Causevic et al., 2006) but for use in epiallele discovery, and as a tool for directed crop epigenetic improvement. In this study we combine MSAP and msGBS (Xia et al., 2014) to survey for *C*-methylation perturbations associated with the micropropation of elite clones of cassava (*M. esculenta*).

### Analysis of Genetic/Epigenetic Variability using MSAP

Both quantitative (**Figure 1**) and qualitative analysis (**Table 5**) of MSAP generated profiles showed clear separation of all tissues in all five varieties studied. Higher levels of diversity and divergence were observed when using *Msp*I than *Hpa*II (**Tables 2** and **4**). Care should be exercised before tentatively assigning this variability as likely to have arisen through genetic or epigenetic causes. Moreover, polymorphic markers between propagation systems that appear in the profiles generated of both isoschizomer restriction enzymes (*Hpa*II and *Msp*I; **Figure 2**) could be caused by either a genetic or an epigenetic change. Conversely, a polymorphic marker detected by only one of the enzymes can only be epigenetic in nature (Pérez-Figueroa, 2013). Application of this reasoning implies that variation at the propagation system diagnostic 55 epilocus could be explained by both genetic and epigenetic changes whereas that of epilocus 101 was due to differential methylation arising from the tissue culture conditions (**Figure 2**). However, since the chance of a genetic mutation occurring at exactly the same location on more than one occasion is extremely low (Rodríguez López et al., 2010b) combined with the fact that these two markers polymorphic were found in three different cultivars implies that both markers probably have an epigenetic origin rather than one caused by genetic mutation.

The variability in MSAP profiles seen between DNA extracted from the same tissue type was both modest and consistent,

*Calculated pairwise tissue PhiPT values (epigenetic distances) obtained from MSAP profiles obtained by restricting genomic DNA from ten cultivars of cassava with MspI using GenAlex 6.5.1 AMOVA. The values show the distance between individual tissues of field cutting lines from bulks of all tissues of in vitro propagated lines. Young lv: distance between field derived young leaves and bulked in vitro tissues, Mature lv: distance between field derived mature leaves and bulked in vitro tissues, and Root: distance between field derived root and bulked in vitro tissues. Probability values based on 9999 permutations are shown in parenthesis. Asterisks indicate lower PhiPT values.*

*Quantitative markers were generated from MSAP profiles peak heights obtained by restricting genomic DNA from five cultivars of cassava (in vitro and in the field) with MspI and HpaII. Column Field vs. in vitro shows the number of significantly different (p* > *0.005) epimarkers between samples of the same cultivar grown in vitro or in the field. Columns Tissues (in vitro) and Tissues (Field) show the number of significantly different (p* > *0.005) epimarkers between tissues grown either in the field or in vitro.*

regardless of the propagation system used to produce the plants (**Figure 1**). This finding suggests that these DNA methylation changes induced by micropropagation are not random, as would be expected for genetic somaclonal variation (Bairu et al., 2011) and so more likely to be associated with methylation events associated with cell and tissue differentiation. Circumstantial evidence in support of this inference can be taken from the PhiPT distance estimate, which showed that samples from *in vitro* nodal micropropagation ramets were always (epigenetically) closer to young leaves of their field counterparts (**Table 4**).

Genetic variation induced during *in vitro* nodal micropropagation cannot be ruled out in this study. However, our results suggest that the majority, if not all, the variability detected using MSAPs is epigenetic in nature. This is supported by the lack of higher levels of variation between micropropagated samples than in field grown samples and the fact that the observed variability seems to be conserved between different plants and between different varieties (**Figure 1**) and also by previous studies that show that micropropagated plants using this approach present high levels of genetic stability (Lata et al., 2010) but measurable levels of epigenetic variability (Baranek et al., 2010).

Our MSAP results suggest that (1) *in vitro* nodal micropropagation introduces *de novo* variability in the global methylation patterns; (2) micropropagation induced epigenetic variability does not seem to be random.

### Analysis of Epigenetic Variability using GBS Data

Most studies of the epigenetic basis of somaclonal variation have used MSAPs to characterize culture-induced epigenetic variation. This technique is reliable and does not require previous knowledge of the studied organism. Conversely, it presents the disadvantage that the generated markers are anonymous. It is possible to isolate and sequence the differential markers (Massicotte et al., 2011), although the process can be cumbersome, expensive and time-consuming (Schrey et al., 2013), especially when many markers and samples are involved. The use of Next-generation sequencing can significantly reduce the cost of epiallele sequence characterization. The recent development of GBS (Elshire et al., 2011; Poland et al., 2012) and its methylation-sensitive version (ms-GBS; Xia et al., 2014) has allowed for a simple, time and cost effective system for the sequencing of multiple DMRs in non-model organisms.

Our study uncovered 105 unique sequences (Supplementary Table S4; 0.03% of those generated) with different levels of methylation between propagation systems. Although total sequence reads were similar between systems, the vast majority of differential tags (101/105) were more abundant among ramets recovered from field cuttings, suggesting again lower global levels of methylation in field grown plants. This is in contradiction with previous evidence suggesting that in vitro culture is related to low DNA methylation (Valledor et al., 2007). However, deciphering global hyper/hypomethylation from restriction products is not a reliable approach (Fulnecek and Kova ˇ ˇrík, 2014). What is more, other studies have shown that methylation levels during in vitro propagation are related to the donor tissue (Fang et al., 2009; Wang et al., 2012), to the length of the culture (Diaz-Sala et al., 1995; Rodríguez López et al., 2010a,b), and the media components (LoSchiavo et al., 1989; Arnholdt-Schmitt, 1993).

FIGURE 2 | Differential peak intensity of informative epiloci (55 and 101 bp) between *in vitro* tissues and field tissues for three cassava cultivars (Kibandameno, Kiroba, Kizimbani). Peak intensities were obtained from MSAP profiles generated restricting genomic DNA from three different tissues (i.e., young leaf, mature

leaf, and roots) in three plants grown either *in vitro* or in the field with *Msp*I (A) and *Hpa*II (B) and amplifying using primer combination *Hpa*II2.2/*EcoR*I5. Box plots show the average normalized intensity scores for a 55 and a 101 bp MSAP fragment selected using. Asterisk indicates *p* < 0.005.

TABLE 6 | Number of significantly different msGBS sequences between *in vitro* and field grown samples.

*"In vitro" and "Field" columns indicate the number of sequences with a higher number of reads from samples grown under each condition.*

BLAST analysis against the cassava genome of the micropropagation induced DMRs generated in this study yielded significant hits for 102 sequences of which the 89 top hits were each associated to a gene (i.e., mapped within a window of 5 kb from the gene; Supplementary Table S5). BLAST results indicate that the homologs of these fragments are involved in many processes, including cell differentiation, plant development, sugar metabolism, nucleic acid methylation, stress response, photosynthesis, cell wall modifications, signaling and transposon activation (Supplementary Table S5). However, it is important to remember that the mere presence of differential methylation in or around a gene is not sufficient evidence to infer that expression of the gene is actually regulated by methylation. There are nevertheless enticing hints to suggest that these candidate methylation markers for propagation system may indeed also play a role in metabolic divergence between field cuttings and meristem micropropagated plants.

A series of studies have implicated DNA methylation in the regulation of genes controlling pathways in plant developmental progression or tissue differentiation (Messeguer et al., 1991), during embryogenesis, seed formation (Xiao et al., 2006), apical dominance regulation, flowering, and floral and leaf formation (Finnegan et al., 1996). Several differentially abundant loci identified in the present study showed high sequence homology to loci in cassava that have been previously implicated in cell differentiation and development: CWF19, XPMC2, EXO70, TAP42-like, Sterile alpha motif (SAM) domain-containing protein, AP3M, Enhancer of polycomb-like transcription factor protein, cassava protein containing a transcription factor UCC1, GFS9, GT-2, EMB71, and ARF2 (Supplementary Table S5). Even if it were shown that the changes in methylation among loci identified are causally linked to changes in gene expression, further work would still be required to establish whether such changes are sufficient to cause a biological meaningful change in cell metabolism and phenotype. Once again, however, there are some grounds to reason that at least some loci might.

A number of studies have shown that DNA methylation plays a central role in gene expression and plant development under stress (for extensive reviews see Chinnusamy and Zhu, 2009; Kinoshita and Seki, 2014). Not surprisingly perhaps, abiotic stresses like those encountered under *in vitro* culture conditions have been found to impose an effect on DNA methylation and have been correlated with subsequent organogenesis (Us-Camas et al., 2014). However, the comparative paucity of marks that appear at higher abundance in meristem culture (just 4 of the 105 generic marks; Supplementary Table S4) suggests that changes of this type lay in the minority. Explanation is therefore required for the far more commonly encountered appearance of marks among the field cutting samples only. Perhaps the most plausible hypothesis for this divergence lay in the more variable living environment experienced by field cutting plants when compared with the more homogeneous environment in culture.  It is therefore, not difficult to presume that *in vitro* grown plants will not experience the same pathogen or abiotic stress-induced DNA methylation marks as will those grown exposed to pathogens in the field. In our study, 13 of the 105 differentially methylated loci mapped to locations of the cassava genome associated to known stress response genes [i.e., Calcium-dependent lipid-binding (CaLB domain) family protein, Plastocyanin-like domain, Leucine-Rich Repeat Receptor-Like Protein Kinase (LRR-RK), Disease resistance protein (TIR-NBS-LRR class) family, ATBCB, ATHCHIB, ATTTM2, MLP-Like Protein 28, ATMGL, ATBZIP1, XTH1, Peroxidase superfamily protein, ATATG18F and HT1; Supplementary Table S5]. Plants have evolved two different strategies involving LRR proteins to perceive microbial pathogens. LRR-RKs are transmembrane host-encoded patternrecognition receptors that directly recognize pathogens while NBS-LRR indirectly recognize pathogen effectors by sensing their effects on plant target proteins (Yu et al., 2013). The latest has been shown to be regulated by DNA methylation (Yu et al., 2013). Previous studies have shown that overexpression of NBS-LRRs induces a severe drop in fitness (Tao et al., 2000). Our results show several of the 105 candidate loci associated to both types of LRRs, so it would be tempting to speculate that the agronomic performance of micropropagated plants could be potentially affected if the observed changes on DNA methylation led to the overexpression of such genes. Calcium lipid-binding domain (CaLB domain) proteins are repressors of abiotic stress response in plants (de Silva et al., 2011) and have been shown to be regulated by environmental conditions through DNA methylation (Dubin et al., 2015). Curiously, our results showed one differentially abundant sequence matched to the xyloglucan endotransglucosylase/hydrolase1 (XTH1) gene. Previous studies have shown that XTHs have a function on cell wall modifications and that changes on their DNA methylation levels are associated to colonization of potato plants by beneficial bacterial endophytes (Da et al., 2012).

In all, 45 of the 89 differentially abundant sequences matched to genes that have been previously reported to be: regulated by DNA methylation (14) or other epigenetic mechanisms (10), methylated or differentially methylated under different growing conditions (4), associated to the regulation of DNA methylation (9) or other epigenetic mechanisms (1) and DNA binding proteins affected by methylation of their target sequences (3) (For references see Supplementary Table S5). Which (if any) of these genes is playing a role in a possible divergence in cell metabolism and phenotype between plants replicated by meristem-propagation and field cutting warrants further attention. Looking further ahead, identifying the developmentally important genes whose expression is sensitive to culture growth conditions may ultimately allow for the development of new culture regimes that yield regenerants with the lowest possible incidence of off-types. In the shorter term, however, the provision of methylation marks that consistently diverge in abundance between plants propagated by meristem culture and those recovered by field cuttings could have utility in the optimization of *in vitro* meristem propagation protocols and also in the diagnosis of the origin of clonal stocks.

## Acknowledgments

This study was carried out as part of a biosafety capacitybuilding project for sub-Saharan Africa, being implemented by the ICGEB. This report is based on research funded by the Bill & Melinda Gates Foundation. The findings and conclusions contained within are those of the authors and do not necessarily reflect positions or policies of the Bill & Melinda Gates Foundation nor the ICGEB.

## Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2015.00590

## References

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2015 Kitimu, Taylor, March, Tairo, Wilkinson and Rodríguez López. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Maize miRNA and target regulation in response to hormone depletion and light exposure during somatic embryogenesis

Elva C. Chávez-Hernández, Naholi D. Alejandri-Ramírez, Vasti T. Juárez-González and Tzvetanka D. Dinkova\*

Departamento de Bioquímica, Facultad de Química, Universidad Nacional Autónoma de México, Mexico City, Mexico

### Edited by:

Clelia De-la-Peña, Centro de Investigación Científica de Yucatán, Mexico

### Reviewed by:

Abu Hena Mostafa Kamal, National Agriculture and Food Research Organization, Japan Soumitra Paul, Krishnagar Government College, India

### \*Correspondence:

Tzvetanka D. Dinkova, Departamento de Bioquímica, Facultad de Química, Conjunto E, Universidad Nacional Autónoma de México, 04510 Mexico City, Mexico<EMAIL_ADDRESS>
### Specialty section:

This article was submitted to Plant Biotechnology, a section of the journal Frontiers in Plant Science

Received: 27 April 2015 Accepted: 07 July 2015 Published: 22 July 2015

### Citation:

Chávez-Hernández EC, Alejandri-Ramírez ND, Juárez-González VT and Dinkova TD (2015) Maize miRNA and target regulation in response to hormone depletion and light exposure during somatic embryogenesis. Front. Plant Sci. 6:555. doi: 10.3389/fpls.2015.00555 Maize somatic embryogenesis (SE) is induced from the immature zygotic embryo in darkness and under the appropriate hormones' levels. Small RNA expression is reprogrammed and certain miRNAs become particularly enriched during induction while others, characteristic to the zygotic embryo, decrease. To explore the impact of different environmental cues on miRNA regulation in maize SE, we tested specific miRNA abundance and their target gene expression in response to photoperiod and hormone depletion for two different maize cultivars (VS-535 and H-565). The expression levels of miR156, miR159, miR164, miR168, miR397, miR398, miR408, miR528, and some predicted targets (SBP23, GA-MYB, CUC2, AGO1c, LAC2, SOD9, GR1, SOD1A, PLC) were examined upon staged hormone depletion in the presence of light photoperiod or darkness. Almost all examined miRNA, except miR159, increased upon hormone depletion, regardless photoperiod absence/presence. miR528, miR408, and miR398 changed the most. On the other hand, expression of miRNA target genes was strongly regulated by the photoperiod exposure. Stress-related miRNA targets showed greater differences between cultivars than development-related targets. miRNA/target inverse relationship was more frequently observed in darkness than light. Interestingly, miR528, but not miR159, miR168 or miR398, was located on polyribosome fractions suggesting a role for this miRNA at the level of translation. Overall our results demonstrate that hormone depletion exerts a great influence on specific miRNA expression during plant regeneration independently of light. However, their targets are additionally influenced by the presence of photoperiod. The reproducibility or differences observed for particular miRNA-target regulation between two different highly embryogenic genotypes provide clues for conserved miRNA roles within the SE process.

Keywords: hormone depletion, maize, miRNA, photoperiod response, polyribosomes, somatic embryogenesis

## Introduction

Maize (Zea mays L.) is one of the most important widely cultivated cereal crops, used as valuable source for human food, livestock feed, and raw material for the industry (Huang et al., 2002). Besides its agricultural and economic relevance, it has been a major model system in plant genetics and improvement. Common methods for maize transformation involve a process known as somatic embryogenesis (SE). SE begins with callus induction, a process characterized by re-organization and re-structuring of somatic cells in the presence of the synthetic auxin 2,4-D (Zimmerman, 1993). The second step of SE is plant regeneration. Plants can be regenerated when somatic embryos are depleted of external hormones and are exposed to photoperiod (Garrocho-Villegas et al., 2012). SE has been advantageously used in many plant species for clonal propagation, plant transformation and genetic improvement (Stasolla and Yeung, 2003).

Plant regeneration through maize SE was first reported by Green and Phillips (1975) using immature embryos as the initial explant. Various conditions for callus induction and plant regeneration have been tested since then, encountering the embryogenic potential highly dependent on explant and maize genotype used (Obert et al., 2009; Shen et al., 2012). The immature embryo's developmental stage and size, usually from 12 to 18 days upon pollination, appear as particularly relevant to generate embryogenic callus type II associated with high plant regeneration frequency over long subculture periods (Armstrong and Green, 1985). A few inbred lines have been reported with good embryogenic potential (Armstrong and Green, 1991; Jakubeková et al., 2012; Shen et al., 2012). Creole varieties VS-535 and H-565 (Costeño mejorado), derived from the germplasm of Mexican landrace Tuxpeño, are highly embryogenic and their plant regeneration frequency could be maintained for over 2 years of subculture (Garrocho-Villegas et al., 2012). However, while VS-535 has been successfully used for more than 15 years in tissue culture, H-565, released in 2008, has shown variable behavior in callus proliferation and plant regeneration responses through SE.

The mechanisms underlying gene activation/repression in the SE process have been poorly characterized (Salvo et al., 2014). Recently, a growing number of reports evidenced essential regulatory roles for microRNAs (miRNAs) in plant developmental and differentiation processes, including zygotic embryogenesis (Nodine and Bartel, 2010; Willmann et al., 2011), hormone signaling (Guo et al., 2005; Reyes and Chua, 2007) and stress response (Sunkar and Zhu, 2004). MicroRNAs are small 21–22 nt RNAs derived from longer precursors by Dicer-like (DCL) endonuclease activity and recruited to protein complexes by Argonaute (AGO) to target specific mRNA repression, either through degradation or translation inhibition (Llave et al., 2002; Li et al., 2013b).

 Most of these reports used a high-throughput sequencing technology to compare the presence of conserved and species-specific miRNAs in the embryogenic callus (EC) before and during different stages of SE. Overall, an up-regulation of specific miRNAs was observed at particular stages of the somatic embryo differentiation. For most of the analyzed plant species, the proliferative embryogenic callus is characterized by low expression of many miRNAs related to flowering and leaf development, while stress-related miRNAs are increased in comparison to the original explant. Although miRNA expression patterns are altered upon callus induction, between embryogenic and non-embryogenic callus, as well as during SE and differentiation, the enriched or decreased miRNA species are particular for each of the plant species analyzed. Therefore, it is important to approach the miRNA regulation in the context of each plant SE system, taking into account the specific conditions and genotypes used in dedifferentiation as well as in plant regeneration induction.

Shen et al. (2013) analyzed by deep sequencing miRNA expression patterns upon callus induction from the maize inbred line 18-599R. They identified miR528, miR156, miR166, miR168, miR390, miR164, miR167, miR398, miR397, miR408, and miR319 as the most abundant during dedifferentiation. These miRNAs increased as the embryos were dedifferentiated into calli, with the exception of miR166 and miR167 that decreased. A degradome analysis indicated that most of their targets are involved in hormone signaling transduction pathways. In a previous study performed on VS-535-derived EC, we found that development-related miRNAs such as miR156, miR159, miR164 and miR168 decreased as the length of subculture increased, while stress-related miRNAs such as miR397, miR398, miR408, and miR528 remained highly expressed (Dinkova and Alejandri-Ramirez, 2014). However, the regulation on miRNA expression and their target mRNAs has not been explored in maize plant regeneration through SE. In this study, we addressed the impact of light exposure and hormone depletion as plant regeneration cues on maize miRNA levels and target gene expression. Taking into account the relevance of maize genotype on plant regeneration success through SE, we investigated how general the observed changes were by comparing two closely related cultivars, VS-535, and H-565.

## Materials and Methods

### Callus Induction and Subculture

To induce embryogenic callus type II, immature embryos were collected at 15–18 days after pollination from two closely related maize cultivars, VS-535 and H-565. These cultivars have been derived from the Tuxpeño landrace germplasm (Márquez-Sánchez, 2008). Tuxpeño (VS-535) was previously shown to display high embryogenic potential during in vitro culture (Garrocho-Villegas et al., 2012), while Costeño mejorado (H-565) was introduced more recently (Márquez-Sánchez, 2008) and its behavior in tissue culture has been variable (unpublished data). The ears from the middle part of the husk (similar developmental conditions) were gently washed, first with 70% ethanol (1 min); then with 50% bleach solution (15 min); and three times with sterile deionized water. Next, the immature embryos were dissected and placed on a Petri dish with sterile deionized water and 1 g L−1 Cefotaxime. Thirty embryos (embryo axis side down) were placed per Petri dish on N6I medium (Supplementary Material, Data sheet 1) and maintained for 3 weeks under darkness at 25 ± 2 ◦C. Upon this time, proembryogenic masses were selected for subculture on proliferation medium N6P (Supplementary Material, Data sheet 1). Every 3 weeks the embryogenic callus was subcultured on fresh N6P medium.

### Plant Regeneration

To test the effects of hormone depletion and light on miRNAmediated regulation during plant regeneration, the embryogenic callus was subjected to stage hormone depletion under darkness or light photoperiod (16 h light/8 h dark). During the first stage, the 2,4-D and kinetin concentrations were half-reduced and during the second stage (2 weeks after the first regeneration subculture) hormones were omitted from N6P. Samples were collected 1 week upon each subculture and stored at −70◦C until used. Three biological samples were collected from each stage for RNA extraction. Regenerating callus was maintained on N6P devoid of hormones with every 2 weeks-subcultures until plantlets appeared under photoperiod. Plantlets were subcultured on MS medium (Murashige and Skoog, 1962; Supplementary Material, Data sheet 1).

### Total and Polysomal RNA Isolation

Total RNA was isolated from triplicates with Trizol reagent (Invitrogen, USA) according to the manufacturer's instructions. For polyribosome fractionation 5 g of embryogenic callus were pulverized in liquid nitrogen with mortar and pestle. The powder was suspended in 25 mL of lysis buffer (200 mM Tris-HCl pH 8.5, 50 mM KCl, 25 mM MgCl2, 2 mM ethylene glycol tetraacetic acid and 0.05 mg mL−1 cycloheximide) and clarified by centrifugation at 15,000 rpm for 15 min. The supernatant was layered onto 4 mL sucrose cushion buffer (50 mM Tris-HCl pH 8.5, 25 mM KCl, 10 mM MgCl2, 60% sucrose and 0.05 mg mL−1 cycloheximide) and centrifuged at 45,000 rpm in a 75Ti rotor (Beckman Coulter, Mexico City, Mexico) for 3 h to concentrate ribosomes. The ribosomal pellet was suspended in 0.5 mL of DEPC water, layered onto 15–60% continuous sucrose gradient and centrifuged in SW-40 rotor (Beckman) at 36,000 rpm for 2.0 h. Fractionation and absorptivity at 260 nm of the gradient was performed in an Auto Densi-flow system (Labconco, Kansas City, MO, USA) connected to Econo UV Monitor EM-1 (BioRad). RNA was isolated from each fraction as described previously (Martinez-Silva et al., 2012). The RNA quality was tested by agarose gel electrophoresis and the concentration was measured with Nanodrop.

### Northern Blot

Ten micrograms of total RNA were separated by electrophoresis on 15% polyacrylamide gels including 7 M urea, transferred to Hybond-N+ membranes (GE Healthcare Life Sciences, USA) and hybridized in ULTRAhyb R -Oligo hybridization buffer (Ambion, USA) with oligonucleotide probes for each miRNA (Supplementary Material, Table S1) end-labeled with [γ-32P]ATP (PerkinElmer, USA) by T4 polynucleotide kinase (NEB, USA).

### miRNA Target Prediction

Targets of zma-miRNAs were predicted by psRNATarget program (http://plantgrn.noble.org/psRNATarget/) with default parameters (Dai and Zhao, 2011), except for the Maximum expectation parameter which was set at five to get a higher prediction coverage. The selected genomic library for the target search was: "Zea mays (maize), transcript, NSF-funded Maize Genome Sequencing Project, Release 5a, filtered set." Target gene descriptions were retrieved from biomart (Kasprzyk, 2011). Some of the predicted miRNA targets (Supplementary Material, Table S2) were chosen to evaluate their expression levels during plant regeneration (**Table 1**) according to the following criteria: (1) appropriate miRNA:target pairing in the seed region; (2) annotation of the target as protein coding transcript; and (3) experimental evidence as miRNA target in maize and/or in other plant species.

## qRT-PCR

Total RNA was extracted from two independent biological samples, treated with RQ1 RNase-Free DNase (Promega, USA) and reverse-transcribed using the ImProm-II™ reverse transcription system (Promega, USA). Each RNA sample was reverse-transcribed in two replicate reactions. Quantitative PCR (qPCR) was performed on the two biological samples for each genotype with three technical replicates per cDNA, using the Express GreenER qPCR reagents (GE Healthcare Life Sciences, USA) in a 7500 Real-time PCR System (Applied Biosystems, USA). Specific primers for qPCR were designed using Primer3Plus (Untergasser et al., 2007) and are available in Supplementary Material, Table S3. Relative abundance was calculated using the 2−11Ct method (Livak and Schmittgen, 2001). Target levels were normalized by rRNA 18S as internal housekeeping control and then compared to the levels found for the initial tissue (100% hormones, darkness). qRT-PCR foldchange data were summarized as Mean + Standard Error. The results obtained for each condition were compared performing a two-way (light × hormone)-ANalysis Of VAriance (ANOVA). The significance of mean difference within and between the groups was retrieved using Tukey Honestly Significance Difference (HSD) at P < 0.05.

## Results

### Plant Regeneration through Somatic Embryogenesis

Staged hormones depletion combined with light photoperiod promotes plant regeneration through SE in maize (Jakubeková et al., 2012; Garrocho-Villegas et al., 2012). Both signals exert important effects at molecular, biochemical, and physiological levels leading to a developmental switch from highly proliferating dedifferentiated tissues to fully differentiated plantlets. In agreement with this, clear morphological changes were observed upon hormone depletion in both, light photoperiod and darkness (**Figure 1**). Particularly, globular compact structures appeared at 50% hormone (2,4-D and kinetin) reduction in light and 0% hormones in darkness (**Figure 1E** and **Figure 1C**, respectively). The light presence rapidly (< 24 h) induced pigment deposition in the callus. For both conditions, the globular structures became elongated, but while under light they eventually resulted in plantlets, in darkness only organogenesis was appreciated (**Figure 1G** vs. **Figure 1D**). Fully developed plantlets were derived from 2 years-subcultured EC at 6–7 weeks upon hormones depletion for both, VS-535 and H-565. However, qualitative differences during plant regeneration were evident

### TABLE 1 | Analysis of miRNA targets.

miRNA targets were predicted with psRNAtarget (http://plantgrn.noble.org/psRNATarget). The Expectation (E) value refers to the score of complementarity between each miRNA and target transcript. The UnPair Energy (UPE) value refers to the energy required for secondary structure melting around the target site. A lower UPE implies higher possibility for establishing a contact between miRNA and target mRNA, as well as for AGO-mediated cleavage.

between these genotypes. For instance, the greenish color was more intense for VS-535 callus and a greater proportion of plantlets were obtained per gram of tissue for this variety (data not shown).

### Specific miRNA Expression Is Mostly Affected by Hormones Depletion

Previous studies in maize long-term subcultured EC indicated that miR156, miR159, miR164, miR168, and miR319 importantly reduce their levels in subcultures maintained for more than 18 months (Dinkova and Alejandri-Ramirez, 2014). On the other hand, the stress-related miR397, miR398, miR408 and miR528, become enriched upon callus induction and remain at high levels once the proliferation is established. This is in accordance with the proposal that specific miRNA expression in the undifferentiated EC associates with their proliferation maintenance, suggesting a switch from this expression pattern during plant differentiation (Luo et al., 2006). To evaluate which stimulus promotes miRNA expression changes, we tested the effect of hormone depletion in darkness or light for two independent maize 2 year-subcultured EC, owing that only in the presence of light plant regeneration could be achieved.

In VS-535-derived EC, most of the analyzed miRNAs increased in response to hormone half reduction (50%), regardless the presence of light (**Figure 2A**; lanes b and d). Although similar behavior was observed in H-565-derived callus, fold changes with respect to 100% hormones were distinct if compared to VS-535 (**Figure 2B**). For instance, at least two-fold increase in miR156, miR164, miR168 and miR408, was observed in VS-535, 50% hormones, while a modest 1.2-1.4-fold increase under the same conditions was evident for H-565 (**Figure 2** and Supplementary Material, Figure S1). The next stage, 0% hormones, implied a sharp reduction for most miRNAs in VS-535, while in H-565 the reduction was lower or not observed depending on the miRNA. The levels of miR156 and miR164 decreased by three-fold (with respect to 50% hormones) for VS-535, but showed no change for H-565 (**Figure 2**, lanes c vs. b and e vs. d). Curiously, miRNA levels in the initial tissue (100% hormones) at 2 year-subcultures were also different between maize genotypes (data not shown). The effect of hormone depletion on miRNA expression was greater in the presence of light, at least for VS-535. miR398 decreased almost by 40 fold between 50 and 0% hormones in VS-535 against a modest two-fold decrease in H-565.

An exemption to most miRNA expression patterns was miR159, showing little or no change in response to hormones depletion in both cultures. In addition, while miR168 importantly increased (around two-fold) in 50% hormones for VS-535, it remained at similar levels for H-565. Regarding the differences observed between genotypes, it is important to stress out that

subculture.

tissues cultured in vitro are highly heterogeneous. Although for RNA analysis, EC was always sampled a week upon subculture, and visibly embryogenic differentiating tissues were selected, the heterogeneity inherent to each cultivar could not be avoided. The miRNA expression levels achieved during hormone depletion were maintained or further decreased in plantlets regenerated from VS-535 (**Figure 2**, lane f). Particularly, the stress-related miR397, miR398, miR408, and miR528 showed about two-fold reduction in VS-535 fully developed plantlets with respect to dedifferentiated tissues (100% hormones). However, in H-565 derived plantlets miR397 and miR398 remained expressed at higher levels.

(C,F). Further subcultures were performed every 2 weeks until plant

### Light Is a Major Stimulus Affecting Development-related miRNA Target Levels during Maize SE

Several of the miRNA targets analyzed in this study (SBP23, GAMYB, CUC2) encode for transcription factors known to participate in plant developmental switches including zygotic embryogenesis (**Table 1**). Others represent enzymes involved in plant stress response (targets of miR397, miR398, miR408, miR528) or the miRNA biogenesis pathway itself (miR168). According to this, we separated the results from qRT-PCR in development-related and stress-related miRNA targets (**Figures 3**, **4**). The levels of GRMZM2G039455\_T01, an AGOlike transcript also termed AGO117 or AGO1c (Qian et al., 2011) were considered within the development-related targets, but according to its proposed function in miRNA biogenesis it actually corresponds to any of the subdivisions.

Development-related target mRNA levels displayed mostly contrasting behavior between darkness and photoperiod. The Squamosa Promoter Binding protein (SBP)-like transcript GRMZM2G126018\_T01 (SBP23) targeted by miR156 showed significantly higher levels upon hormone half reduction under darkness than in the presence of light (**Figure 3**, b–c vs. d–e). Under darkness, SBP23 inversely mirrored miR156 changes due to hormone depletion (**Figure 3A**, lanes b–c; and **Figure 3B**, lanes a–b). However, its levels were higher in this condition

FIGURE 2 | Changes in miRNA levels occur upon hormone depletion under both, darkness and photoperiod conditions. Northern Blot assays were performed to evaluate miRNA abundance during plant regeneration from somatic embryos in two maize cultivars: VS-535 (A) and H-565 (B). Northern blot signals were acquired by densitometry and normalized according to 5S rRNA (Supplementary Material, Figure S1). The value corresponding to EC in 100% hormones (2,4-D and kinetin) was set as 1. Fold changes in 50% hormones (b, d); 0% hormones (c, e) and plantlet (f) with respect to 100% hormones were plotted on charts.

than in the presence of light. Such behavior correlates with the differentiation status of the tissue and with major influence of light on SBP23 expression. Since miRNA changes were essentially similar between darkness and light, but levels of the target displayed significant differences between these conditions, additional levels of regulation are likely operating to downregulate SBP23 during differentiation. Although we only tested one SBP-like target of miR156, 11 out of 26 SBP-like maize genes were predicted as targets for this miRNA, all of them exhibiting the same target sequence as SBP23 (Supplementary Material, Table S2). Targeting of SBP-like transcripts by miR156 is highly conserved in plants (Rhoades et al., 2002; Xie et al., 2006; Gandikota et al., 2007) and has been reported as relevant for vegetative to reproductive phase changes and plastochron length (Wu and Poethig, 2006; Wang et al., 2008). For maize, it has been shown that miR156 gene family overexpression renders plants with increased number of leaves and delayed flowering (Chuck et al., 2007).

Similar to SBP23, a miR164 previously validated target, GRMZM2G393433\_T01 (NAC-domain transcription factor CUC2 or NAC107, Zhai et al., 2013; Liu et al., 2014) was more abundant in darkness than light for either genotype (**Figure 3**, lanes a–c). However, while for H-565 the miRNA/target levels followed an inverse tendency with a consistent decrease of CUC2 upon hormone removal and light exposure, for VS-535 a correlation was observed only under light (**Figure 3A**, lanes d–f). Global transcriptome analyses indicated that CUC2 expression was restricted to immature inflorescence and ear in maize (Sekhon et al., 2011; Fan et al., 2014). Therefore, it is interesting that during SE this transcript was particularly abundant at 100% (H-565) or 50% (VS-535) hormones in darkness, both stages characterized by undifferentiated tissues (**Figure 1**).

Although miR159 and miR168 displayed slight changes in response to plant regeneration stimuli, their target levels (GRMZM2G139688\_T01, Giberellic Acid-responsive MYB, GA-MYB for miR159 and GRMZM2G039455\_T01, AGO1c for miR168) showed significant fluctuations at least for 0% hormones in darkness and photoperiod (**Figure 3**, lanes c and e). For both genotypes GA-MYB levels increased at 0% hormones in darkness, while strongly decreased for the same condition in light photoperiod. This correlates with the observation of drastic morphological changes in somatic embryos (**Figure 1**). miR159 is one of the most abundant miRNAs in maize EC (Shen et al., 2013) and its levels change the least (**Figure 2**). Nevertheless, an inverse correlation between GA-MYB and slight miR159 fluctuations were usually observed for H-565 and to a lesser extent for VS-535, suggesting that miR159 expression during SE might be important to keep GA-MYB levels low when not needed. Regarding normal maize plant development, GA-MYB expression has been detected in 24 h germinating seed and during seed development in the endosperm and pericarp, in meiotic tassel and in anthers (Sekhon et al., 2011). In addition, miR159-mediated regulation on MYB transcripts might depend on the tissue and developmental stage, at least in Arabidopsis (Woodger et al., 2003; Alonso-Peral et al., 2012).

miR168 is known as master regulator of the general miRNA pathway since it is required for fine-tuning the AGO1 levels in Arabidopsis (Vaucheret et al., 2006). In maize, at least five AGO1 putative genes were previously identified (Qian et al., 2011) in contrast to only one AGO1 encoding gene in Arabidopsis. AGO1c shows the greatest identity with the Arabidopsis AGO1 homolog and has conserved the miR168 target sequence (**Table 1**). Similar to miR156 and SBP23, there was poor correlation between miR168 and AGO1c transcript levels. While the greatest increase in miR168 was observed for VS-535 at 50% hormones under photoperiod, AGO1c levels reached a two-fold increase in the same sample (**Figure 3A**, lane d). It was then decreased by four-fold at 0% hormones under the same condition when also miR168 decreased (**Figure 3A**, lane e). A greater reduction in this AGO1c transcript was observed for H-565, 0% hormones and photoperiod (25-fold) without accompanying changes in miR168 (**Figure 3B**, lane e). This suggests a complex regulation on miRNA biogenesis in maize, probably involving additional AGO1 isoforms during plant regeneration.

### Targets of Stress-related miRNAs Display Essentially Contrasting Patterns between Maize Genotypes during SE

Stress-related miRNAs are highly expressed in the process of EC induction for one Chinese maize genotype (Shen

et al., 2013). Stress has been also regarded as one of the main stimuli promoting dedifferentiation as well as SE. Major targets for stress-related miRNAs are copper proteins, multi-copper oxidases, superoxide dismutases and laccases (Jones-Rhoades and Bartel, 2004; Yamasaki et al., 2007; Abdel-Ghany and Pilon, 2008). Here we tested the mRNA levels for maize GRMZM2G146152\_T01 (laccase-like, LAC2) as target for miR397, GRMZM2G058522\_T01 (SOD9) as validated target for

through SE. Transcript levels of miRNA stress-related targets were analyzed by qRT-PCR in the same samples as miRNAs for VS-535 (A) and H-565 (B). The18S rRNA was used as reference control and plotted values represent the expression of each mRNA relative to 100% hormones (a, set to 1). Fold changes are shown accordingly for 50% hormones (b, d); 0% hormones (c, e) and plantlet (f). Filled bars

target: GRMZM2G146152\_T01 (LAC2); miR398 target: GRMZM2G058522\_T01 (SOD9); miR408 target: GRMZM2G384327\_T03 (GR1); miR528 targets: GRMZM2G106928\_T01 (SOD1A) and GRMZM2G107562\_T01 (PLC). Error bars represent + Standard Error; n = 6. Statistically significant differences were identified using Two-Way ANOVA and Tukey as described in Methods at \*P < 0.05.

miR398 (Shen et al., 2013), GRMZM2G106928\_T01 (SOD-1A) and GRMZM2G107562\_T01 (plastocyanin-like protein, PLC) as targets for miR528, and GRMZM2G384327\_T03 (Gamma Response 1 protein, GR1), as previously predicted miR408 target (Li et al., 2013a). It is worth mentioning that although stressrelated miRNAs have been identified as very highly expressed during germination, immature embryo dedifferentiation and stress-response, their target identification has remained elusive in maize. Our selection included both, validated and predicted targets according to the properties indicated in **Table 1**.

The level of miR397 target, LAC2, showed contrasting changes between VS-535 and H-565 genotypes throughout hormone depletion. For VS-535, it increased upon 50% hormone reduction under darkness or light (to a lesser extent), and remained unaltered after complete hormone removal (**Figure 4A**). On the contrary, for H-565 it was strongly reduced in 0% hormones, darkness or light (**Figure 4B**). For other stress-related miRNA targets, such as SOD9 (miR398), SOD1A (miR528) and GR1 (miR408), lower levels were also detected in H-565 than in VS-535, particularly in differentiating tissues under light (**Figure 4A** vs. **Figure 4B**, lanes e and f). There was poor correlation between miR397 and LAC2 levels for either of the embryogenic cultivars, suggesting this transcript is regulated by additional mechanisms during plant regeneration. We observed an inverse correlation between miR398 (increase) and SOD9 (decrease) upon 50% reduction in hormones for both cultivars in darkness, but only for H-565 in light (**Figure 4A** vs. **Figure 4B**). This indicates even for a validated target it is difficult to find strong miRNA/target correlation under light, particularly upon full hormone depletion. Similarly, SOD1A showed a clear inverse correlation with miR528 levels for both, VS-535 and H-565, in darkness but not in light (**Figure 4**, upper panels corresponding to miR528). It is interesting to notice that although SOD9 and SOD-1A catalyze the same reaction, their high and differential transcript expression supports a non-overlapping function in maize plant regeneration through SE. Similar to other stress-related miRNA targets, GR1 levels were different between the tested maize genotypes (**Figure 4**). It increased upon hormone depletion in darkness or light for VS-535, while strongly decreased in the presence of light for H-565 (**Figure 4A** vs. **Figure 4B**). A correlation between GR1 and miR408 was observed for H-565, but not for VS-535, suggesting this target regulation during SE could be genotype-specific.

miR528 is one of the most abundant, stress-related miRNAs in SE for maize (Shen et al., 2013), rice (Luo et al., 2006; Chen et al., 2011), and citrus (Wu et al., 2015). In our miR528 target prediction analysis, several copper-binding protein encoding transcripts, many uncharacterized genes and even transcription factors were found (Supplementary Material, Table S2). Taking into account the expected relevance for miR528-mediated regulation in SE, we tested, in addition to SOD1A, the expression levels of GRMZM2G107562\_T01 that codes for a plastocyaninlike protein (PLC). Unlike SOD1A, PLC showed a very weak correspondence with miR528 changes. However, consistent with other stress-related miRNA targets, its levels differed between VS-535 and H-565 throughout the hormone depletion stages (**Figure 4A** vs. **Figure 4B**, lower panels corresponding to miR528). Curiously, PLC levels were high in darkness, in spite of its predicted photosynthesis-related function (Abdel-Ghany and Pilon, 2008). In this regard, a light-independent induction of photosynthesis genes has been documented throughout embryogenesis stages (from globular to torpedo) in Arabidopsis. Genes involved in energy production comprised the largest upregulated functional group, strongly biased toward components of the photosynthetic apparatus (Spencer et al., 2007). Taking this into account, high expression of PLC is probably required at early stages of SE.

### Polyribosomal Distribution of miR528 in Maize EC

According to the lack of consistent correlation between levels of conserved development- or stress-related miRNAs and their targets (predicted or confirmed) during maize SE and plant regeneration, the regulation by additional mechanisms was evident. Primary regulation at transcription levels was already described for several of the analyzed targets (SBP23, GA-MYB, CUC2), particularly in response to developmental cues and hormones' presence. However, another level of regulation might be exerted by miRNAs at translation without affecting the transcript levels (Brodersen et al., 2008; Lanet et al., 2009; Iwakawa and Tomari, 2013). To assess the possible role of some conserved miRNAs at translation level, polyribosomal fractions were obtained from 2 years-subcultured VS-535 EC and the

RNA abundance and integrity across fractions (F) was measured by absorbance at 254 nm and agarose gel electrophoresis. Northern blotting of selected miRNAs was performed upon RNA isolation from every two-fraction pools resulting in nine fractions. The 5S ribosomal RNA was used as control for the Northern blotting.

presence of miR159, miR168, miR398, and miR528 was analyzed (**Figure 5**). Only miR528 was detected in polyribosomal fractions (F 10–18) suggesting this miRNA could regulate some of its targets at translation level. In Arabidopsis, the presence of miR168 and miR398 in polyribosomes was previously demonstrated in 10 days-old seedlings (Lanet et al., 2009). However, we could not detect these miRNAs in maize polyribosomes from EC, although they had detectable levels in total RNA obtained from this tissue (**Figure 2**, Supplementary Material, Figure S1).

 These studies found that miRNA patterns change upon callus induction, between embryogenic and non-embryogenic callus, as well as during SE and differentiation in a plant species-dependent fashion. For example, while miR171, miR390, and miR398 are preferentially expressed in EC before induction of plant differentiation in rice (Luo et al., 2006); these miRNAs are increased during the differentiation process in citrus (Wu et al., 2011). Even miRNAs from the same family, i.e., miR156a and miR156b, may display differential expression patterns during the SE differentiation process, as demonstrated in a recent study on Larix leptolepsis (Zhang et al., 2012). Hence, the current knowledge based on global miRNA approaches highlights the relevance of exploring particular miRNA landscapes and their target regulation in the context of species-specific SE conditions.

### miRNA Expression in Maize SE

Separated analysis on hormone depletion and light effects during maize SE indicated that miRNA expression patterns are affected mostly by hormones, rather than the light presence. By analyzing the process in two independent embryogenic genotypes we expected to confirm whether the observed miRNA changes are inherent to the process rather than the cultivar. However, we noticed a great influence of the genotype, first on the initial miRNA level appreciated in the 2 years-subcultured EC, and second on the degree of changes registered during hormone depletion. In this sense, miRNA expression regulation resulted more dramatic in VS-535 than in H-565. Strikingly lower levels were detected for all tested miRNAs in VS-535 at 0% hormones and regenerated plantlets under light, whereas H-565 maintained the presence of some of the stress-related miRNAs (miR397 and miR398) higher.  The authors suggested miR164, miR169, miR528, and miR529 might be primarily participating in the process of EC induction through the regulation of targets involved in auxin and gibberellin signaling. However, other miRNAs, significantly up-regulated in the dedifferentiation process, were miR156k, miR168, miR397, miR398, and miR408. We found the same miRNAs transiently increased by the reduction of hormones concentration in half during plant regeneration. However, in the absence of hormones their levels were reduced. These observations support the occurrence of specific miRNA expression readjustments in embryogenic tissues in response to hormone changes in the environment.

One important question is whether the high concentration of specific miRNAs associates with the embryogenic potential of the callus. Shen et al. (2013) found that the initial increase in miRNAs was either enhanced or maintained during the dedifferentiation process. However, in maize EC subcultured for long periods (up-to 2 years) we have found a gradual reduction in miR156, miR164 and miR168 levels without impairment on the callus embryogenic potential (Dinkova and Alejandri-Ramirez, 2014). In agreement, the data presented here support an initial burst on certain miRNA levels preceding their further decrease during maize SE. Hence, miRNA expression response to hormone changes could be a major factor impacting on the embryogenic potential of maize cultivars during both, dedifferentiation and plant regeneration.

### Correlation between miRNA and Target Levels in Maize SE

miRNAs are known to regulate their target mRNAs by degradation, translation inhibition or both. In plants, an inverse correlation between miRNA and target levels is commonly observed, suggesting that mRNA degradation is the preferred regulatory mechanism. However, global miRNA and degradome sequencing data have shown that not always the degradation products of predicted or even validated targets could be detected in the libraries (Shen et al., 2013; Wu et al., 2015). Conversely, these studies have identified novel targets for conserved miRNAs such as miR156 and miR164, or known targets with novel miRNA sites. An increasing number of reports have also revealed the lack of inverse relationship between miRNAs and their targets, depending on the tissue or the process analyzed (Brodersen et al., 2008; Wu et al., 2011; Alonso-Peral et al., 2012).

Here we found the expected inverse correlation between a miRNA and its predicted and/or validated target is highly dependent on environmental (light presence) and internal (genotype) signals during maize plant regeneration through SE (**Figure 6**). Major differences were observed for developmentrelated miRNA target levels between genotypes under darkness, but not under light where plant regeneration took place. On the other hand, inverse relationship between miRNA/target levels was more easily found in darkness than in light. Therefore, for most of the analyzed targets we propose there is a major influence of additional regulation under a photoperiod, making it difficult to appreciate the effect of miRNAs. On the other hand, increments of both, miRNA and target, upon hormone 50% reduction (e.g., miR156, miR164, miR168) suggest miRNA upregulation might be required to control the levels of transcripts induced during SE.

Accessing miRNA-mediated regulation through quantitative transcript evaluation has been widely used in plants, particularly during SE (Li et al., 2012; Yang et al., 2013; Wu et al., 2011, 2015). However, an inverse miRNA/target relationship has not always been observed, even in the presence of miRNA-mediated

degradation evidence (degradome or 5′RACE fragments). Our finding on miR528 distribution in maize EC polysomal fractions supports the notion that target regulation is probably exerted at multiple levels depending on the developmental process, miRNA, and analyzed target. In agreement, we observed such expected inverse correlation between miR528 and SOD1A target, but not for PLC target. Therefore, for future studies, it would be relevant to include miRNA-target evaluations at protein in addition to transcript levels.

### Physiological Relevance of Development-related miRNA Regulation in Maize SE

miR156 and miR164 have been found as SE-abundant miRNAs in several species, including maize (Li et al., 2012; Shen et al., 2013; Dinkova and Alejandri-Ramirez, 2014; Wu et al., 2015). The miR156-mediated suppression of SBP transcripts is probably required for early SE, as demonstrated for Arabidopsis zygotic embryogenesis (Nodine and Bartel, 2010). Similarly, miR164 initial increase during SE is consistent with maintaining low levels of its CUC2 target during plant regeneration under light. In plants, conserved miR164 targets are NAC transcription factors (Mallory et al., 2004). The NAC family of proteins includes NAM, ATAF1-2, and CUC2. Proteins belonging to the NAC family are involved in many plant developmental processes, such as flowering, embryogenesis, senescence, auxin signaling, secondary wall thickening and others (Fan et al., 2014). NAM was reported as related to the shoot apical meristem and primordium formation in Petunia hybia (Souer et al., 1996) and CUC2 has been involved in Arabidopsis shoot apical meristem development (Aida et al., 1997). This context is consistent with the observed CUC2 reduction in the absence of hormones and in fully differentiated tissues. Therefore, although miR164 and miR156 targets might display contrasting behavior in undifferentiated tissues (darkness, hormones' presence) between maize genotypes, their expression regulation is apparently required for plant regeneration through SE, regardless the genotype (**Figure 6**).

One possibility underlying the contrasting behavior of miRNA targets between darkness and light could be their transcription responsiveness to photoperiod. Another is a differential function of miRNA-mediated silencing pathways depending on light presence/absence during SE. It was previously reported that AGO1 expression is highly induced early in carrot SE and further decreased after the globular-staged embryo (Takahata, 2008). In maize, the presence of several putative AGO1 genes, not all of them exhibiting the miR168 target sequence, further adds a complexity level to the interpretation of miRNAmediated regulation. According to a global transcriptome analysis for the maize B73 line, AGO1c transcript (target of miR168) is highly expressed in differentiating tissues, such as the shoot apical meristem, immature cob, and tassel (Sekhon et al., 2011). It remains high for few days upon pollination, while in the mature and germinating maize seed AGO1c and miR168 are both decreased. Another AGO1 isoform also exhibiting the miR168 target sequence (GRMZM2G441583\_T01, AGO113 or AGO1a) remains highly expressed during seed maturation and germination (Sekhon et al., 2011). The behavior of AGO1c in maize SE is consistent with the transient increase observed in carrot and in early zygotic embryogenesis. We did not analyze other AGO1 transcripts during maize plant regeneration, but it remains possible that they were differentially expressed depending on the genotype and regulated during SE.

### Physiological Relevance of Stress-related miRNA Regulation in Maize SE

Stress-related miRNAs and their targets have been associated with sweet orange callus embryogenic potential (Wu et al., 2011, 2015) and SE in other species (Li et al., 2012). Highly proliferating tissues are thought to produce an excess of reactive oxygen species (ROS) making transcript accumulation of stress-related genes such as superoxide dismutase, cupredoxin, and multicopper oxidases necessary to minimize cell damage and promote SE. Members of the copper superoxide dismutase (CSD, SOD) family are miRNA targets in several plant species (Jones-Rhoades and Bartel, 2004; Jovanovic et al., 2014; Naya et al., 2014 ´ ). These enzymes are in charge of destroying oxygen reactive species accumulating during fast plant growth (Ravet and Pilon, 2013). On the other hand, an Arabidopsis ortholog of the predicted maize miR408 target, GR1, is expressed in mitotically active tissues, such as the shoot apical meristem and floral primordium from unstressed plants in a similar to cell cycle-related gene expression profiles (Deveaux et al., 2000). AtGR1 is responsive to genotoxic stress, such as gamma radiation, and has been proposed to block mitotic cell divisions in irradiated cells to prevent premature entry into mitosis before completion of DNA repair. Therefore, the control of this transcript by miRNA makes sense in highly proliferating undifferentiated callus (Luo et al., 2006; Shen et al., 2013). Here we found that during maize hormone depletion, targets of stress-related miRNAs

## References

Abdel-Ghany, S. E., and Pilon, M. (2008). MicroRNA-mediated systemic downregulation of copper protein expression in response to low copper availability in Arabidopsis. J. Biol. Chem. 283, 15932–15945. doi: 10.1074/jbc.M8014 06200

often displayed higher expression under darkness. In addition, their expression patterns under light showed striking differences between genotypes contrary to what observed for developmentrelated targets. Accordingly, we propose that while stress-related miRNAs are required to control their targets in undifferentiated tissues, their participation during the plant regeneration process could be dependent on the genotype and the physiology of the EC used to initiate SE (**Figure 6**).

## Conclusion

The present study provides important novel information about the separate effects of hormone depletion and light presence on miRNA patterns and their target regulation during plant regeneration from maize embryogenic callus. While development or stress-related miRNAs are responsive to hormone concentration, their targets are additionally influenced by the presence of photoperiod. The reproducibility or differences observed for particular miRNA-target regulation between two different highly embryogenic genotypes provide clues for conserved miRNA roles within the SE process. Future work should aim to approach the mechanism underlying regulation of SE-related miRNAs on their targets and its relevance for the plant regeneration success.

## Author Contributions

Conceived and designed the experiments: EC, NA, TD. Performed the experiments: EC, NA, VJ. Analyzed the data: EC, NA, VJ. Contributed reagents/ materials/ analysis tools: EC, NA, VJ. Wrote the paper: EC, TD. Revised the manuscript: EC, TD. Steered the whole study: TD.

## Acknowledgments

The authors acknowledge the financial support from Programa de Apoyo a Proyectos de Investigación e Innovación Tecnológica, UNAM, IN210912 and IN211215. We wish to thank M.Sc. Ma. Teresa de J. Olivera-Flores for the support and advice in maize in vitro tissue culture and Dr. Carolina Peña-Montes for the help in qRT-PCR experiments.

## Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2015. 00555

Aida, M., Ishida, T., Fukaki, H., Fujisawa, H., and Tasaka, M. (1997). Genes involved in organ separation in Arabidopsis: an analysis od the cup-shaped cotyledon mutant. Plant Cell 9, 841–857. doi: 10.1105/tpc.9.6.841

Alonso-Peral, M. M., Sun, C., and Millar, A. A. (2012). MicroRNA159 can act as a switch or tuning microRNA independently of its abundance in Arabidopsis. PLoS ONE 7:e34751. doi: 10.1371/journal.pone.0034751

embryogenesis reveals altered expression of stress factors and embryogenesisrelated genes. PLoS ONE 9:e111407. doi: 10.1371/journal.pone.0111407

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2015 Chávez-Hernández, Alejandri-Ramírez, Juárez-González and Dinkova. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# The rice enhancer of zeste [E(z)] genes SDG711 and SDG718 are respectively involved in long day and short day signaling to mediate the accurate photoperiod control of flowering time

### *Xiaoyun Liu1, Chao Zhou1, Yu Zhao1, Shaoli Zhou1, Wentao Wang1 and Dao-Xiu Zhou1,2 \**

1 National Key Laboratory of Crop Genetic Improvement, Huazhong Agricultural University, Wuhan, China

2 Institut de Biologie des Plantes, UMR8618, Université Paris-Sud 11, Orsay, France

### *Edited by:*

Raúl Alvarez-Venegas, Centro de Investigación y de Estudios Avanzados del Instituto Politécnico Nacional, Mexico

### *Reviewed by:*

Hairong Wei, Michigan Technological University, USA Donna Glassop, Commonwealth Scientific and Industrial Research Organisation, Australia

### *\*Correspondence:*

Dao-Xiu Zhou, Institut de Biologie des Plantes, UMR8618, Université Paris-Sud 11, 90405 Orsay, France e-mail<EMAIL_ADDRESS>
Recent advances in rice flowering studies have shown that the accurate control of flowering by photoperiod is regulated by key mechanisms that involve the regulation of flowering genes including Heading date1 (Hd1), Early hd1 (Ehd1), Hd3a, and RFT1. The chromatin mechanism involved in the regulation of rice flowering genes is presently not well known. Here we show that the rice enhancer of zeste [E(z)] genes SDG711 and SDG718, which encode the polycomb repressive complex2 (PRC2) key subunit that is required for trimethylation of histone H3 lysine 27 (H3K27me3), are respectively, involved in long day (LD) and short day (SD) regulation of key flowering genes. The expression of SDG711 and SDG718 is induced by LD and SD, respectively. Over-expression and downregulation of SDG711 respectively, repressed and promoted flowering in LD, but had no effect in SD. By contrast, down-regulation of SDG718 had no effect in LD but delayed flowering in SD. SDG711 and SDG718 repressed OsLF (a repressor of Hd1) respectively in LD and SD, leading to a higher expression of Hd1 thus late flowering in LD and early flowering in SD. SDG711 was also found to be involved in the repression of Ehd1 in LD. SDG711 was shown to directly target to OsLF and Ehd1 loci to mediate H3K27me3 and gene repression. The function of the rice E(z) genes in LD repression and SD promotion of flowering suggests that PRC2-mediated epigenetic repression of gene expression is involved in the accurate photoperiod control of rice flowering.

**Keywords:** *Oryza sativa***, heading date, epigenetics, epigenomics, histone methylation, PRC2,** *OsCLF***,** *OsiEZ1*

### **INTRODUCTION**

The control of flowering time is a critical step for successful grain production in rice. Day length is a key factor controlling rice flowering. Most rice cultivars recognize 13.5 h light/10.5 h dark as a critical photoperiod to separate long day (LD) from short day (SD) periods, day length shorter than 13.5 h will greatly induce rice flowering (Itoh et al., 2010; Tsuji et al., 2013). Despite some genes are shared between *Arabidopsis* and rice flowering regulatory pathways, there are considerable differences between the regulation of flowering of both species such as absence of the vernalization pathway in rice (Shrestha et al., 2014)*. Heading date1 (Hd1)*, the rice ortholog of *Arabidopsis CONSTANS* (*CO*), has a dual role in flowering time control. It promotes flowering at SD and represses flowering in LD (Yano et al., 2000). (*Hd3a* is a rice ortholog of the *Arabidopsis* florigen gene *FLOWERING LOCUS T* (*FT*; Yano et al., 2000; Kojima et al., 2002; Tamaki et al., 2007). Hd1 activates *Hd3a* in SD, but in LD Hd1 is converted by phytochrome B to a repressor of *Hd3a* (Yano et al., 2000; Hayama and Coupland, 2004). *Hd1* expression is controlled by *Oryza sativa* GIGANTEA (OsGI), a key factor of circadian rhythms in rice flowering control (Yano et al., 2000; Hayama et al., 2003). There is a second, *Hd1*-independent, photoperiod inductive pathway in rice. *Early*

*hd1* (*Ehd1*), a B-type response regulator that is activated by a SD flowering promoter (i.e., *OsMADS51*), activates the expression of *Hd3a*, *RFT1* (another Rice *FT* gene), and *OsMADS14,* and mainly confers SD-dependent flowering promotion (Doi et al., 2004; Kim et al., 2007). In LD, *Ehd1* is activated by *Ehd2/OsId1/RID1* (*Rice Indeterminate1*) and *OsMADS50*, but repressed by *Ghd7* (*Grain number, plant height and heading date7*; Lee et al., 2004; Matsubara et al., 2008; Park et al., 2008; Wu et al., 2008; Xue et al., 2008; Itoh et al., 2010; Tsuji et al., 2013). *Ehd1* activates the expression of *Hd3a* and *RFT1* which activates *MADS14* and *MADS15* (Komiya et al., 2009). The tight control of expression of the flowering promoter *Ehd1* and the flowering repressor *Ghd7* allows to measure the slight differences in day lengths to control*Hd3a* and*RFT1* transcription with a critical day length threshold (Itoh et al., 2010; Tsuji et al., 2013). Unlike *Hd1* that is conserved with *Arabidopsis CO*, *Ehd1,* and *Ghd7* are evolutionarily acquired rice-specific genes. Therefore, the control of florigen expression is regulated by key mechanisms that involve the regulation of *Hd1* expression, the conversion of Hd1 function on *Hd3a* expression, and modulation of *Ehd1/Ghd7* expression.

Recent studies have established a close relationship of epigenetic regulation with flowering. For instance, during the process of vernalization in *Arabidopsis*, polycomb repressive complex 2 (PRC2)-mediated trimethylation of histone H3 lysine 27 (H3K27me3) represses the expression of *Flowering Locus C* (*FLC*) to stimulate flowering (Sung and Amasino, 2004; De Lucia et al., 2008). In rice, SDG724 mediates H3K36me2/3 deposition at *OsMADS50* and *RFT1*, promoting flowering and establishing a difference of functionality between paralogs *RFT1* and *Hd3a* under LD or SD conditions (Sun et al., 2012). SDG725 also mediates H3K36me2/3 deposition at *Ehd2*, *Ehd3*, *OsMADS50*, *Hd3a,* and *RFT1*, promoting flowering under LD or SD conditions (Sui et al., 2013). The PRC2 complex was first discovered in *Drosophila*, which has four core proteins: ENHANCER OF ZESTE [E(z)], SUPPRESSOR OF ZESTE 12 [Su(z)12], EXTRA SEX COMBS (ESC) and P55 (Schuettengruber and Cavalli, 2009). The E(z) protein has the H3K27 methyltransferase activity (Cao et al., 2002; Czermin et al., 2002). *Arabidopsis* contains three E(z) genes [*CURLY LEAF* (*CLF*), *SWINGER* (*SWN*), and *MEDEA* (*MEA*)], three Su(z)12 genes [*FERTILIZATION INDEPENDENT SEED 2* (*FIS2*), *VERNALIZATION 2* (*VRN2*) and *EMBRYONIC FLOWER 2* (*EMF2*)], and one ESC gene [*FERTILIZATION INDE-PENDENT ENDOSPERM* (*FIE*)] (Luo et al., 1999; Ohad et al., 1999; Gendall et al., 2001; Yoshida et al., 2001; Hennig et al., 2003), which form three PRC2-like complexes. The FIS complex, which contains MEA/SWN, FIS2, FIE, and MSI1 (Multicopy Suppressor of IRA1), functions during gametogenesis and seed development. The others are involved in flowering control: the EMF complex that is comprised of CLF/SWN, EMF2, FIE, and MSI1 and is involved in the suppression of early flowering and the VRN complex that plays critical roles in the vernalization pathway by maintaining the high level of H3K27me3 on the *FLC* locus after vernalization (Hennig and Derkacheva, 2009). The VRN complex is associated with VERNALIZATION INSENSITIVE 3 (VIN3, a PHD-domain containing protein) and VIN3-like proteins to form PHD-VRN PRC2 complexes (Wood et al., 2006; De Lucia et al., 2008). The VIN3 protein enhances H3K27me3 throughout the target loci to a level sufficient for stable silencing.

The rice genome contains two genes for E(z) (*OsiEZ1* and *OsCLF*), Su(Z)12 (*OsEMF2a* and *OsEMF2b*) and ESC (*OsFIE1* and *OsFIE2;* Luo et al., 2009). No morphological changes are observed in *osclf* and *osfie1* mutants, while *osfie2* and *osemf2b* mutants display earlier flowering at LD and abnormal flower organs (Luo et al., 2009). Recent results have shown that rice VIN3-like proteins OsVIL2, OsVIL3, or RICE LEAF INCLINATION 2 (LC2, hereafter referred to as LC2) promote rice flowering through the photoperiod pathway (Wang et al., 2013; Yang et al., 2013). These results suggest that PRC2 and PRC2-associated genes are involved in photoperiod regulation of flowering in rice. However, how PRC2-mediated gene repression is involved in accurate photoperiod control of rice flowering is not clear.

In this work we show that the two rice E(z) genes, *OsCLF* (or *SDG711,* Os06g16390, here after referred to as *SDG711*) and *OsiEZ1* (or *SDG718,* Os03g19480, here after referred to as *SDG718*), displayed distinct function in photoperiod regulation of flowering in rice. *SDG718* is induced in SD and represses *OsLF*, a repressor of Hd1 (Zhao et al., 2011), leading to a higher expression of *Hd1* (that activates *Hd3a* in SD) and thus early flowering.

*SDG711* is induced in LD and represses *OsLF*, *Ehd1*, and other flowering-promoting genes leading to late flowering in LD. The data suggested that the two E(z) genes are involved respectively, in LD and SD signals to differentially control key flowering genes expression and contribute to the accurate photoperiod control of flowering time in rice.

### **MATERIALS AND METHODS**

### **PLANT MATERIALS AND GROWING CONDITIONS**

Rice cultivar (*Oryza sativa* spp. *japonica*) "Zhonghua 11" (ZH11) and "DongJin" (DJ) were used for genetic transformation in this study. T-DNA insertion line of *SDG711* (3A-60654.R) was obtained from the Postech rice mutant database (http://signal.salk.edu/cgi-bin/RiceGE). Insertion was confirmed by PCR using the specific primers F and R and a T-DNA left side primer RB2. The primers used for genotyping and real-time PCR analysis are listed in Table S1. The rice plants were grown either in a paddy field in summer in Wuhan (day length >13.5 h) or in controlled growth chambers for 6 week-old under SD (10 h light at 30◦C/14 h dark at 25◦C) or 8 week-old LDs (14 h light at 30◦C/10 h dark at 25◦C) conditions as described previously (Yang et al., 2013).

### **EXPRESSION ANALYSIS BY NORTHERN BLOTS AND RT-PCR**

Total RNA was isolated from rice callus, stems, leaves, flag leaves, shoots, panicles, endosperm, and roots using TRIzol reagent (Invitrogen). Three μg of total RNA were reverse-transcribed in a reaction of 20 μl by using DNase I and SuperScript III (Invitrogen) according to the manufacturer's instruction to obtain cDNA. For northern blotting analysis, fifteen micrograms of total RNA extracted from field-grown rice leaves was separated in 1.2% (w/v) formamide-denaturing agarose gels, then transferred to nylon membranes. Gene-specific probes were labeled with 32P-dCTP using the Random Primer kit (Invitrogen) and hybridized to the RNA blots. The probe of *SDG711* was amplified from *SDG711* cDNA using primers Insitu-SDG711-F and Insitu-SDG711-R (Table S1), resulting in a fragment of 505 bp of the cDNA.

Real-time PCR was performed in an optical 96-well plate that included SYBR Premix EX Taq and 0.5 μl of Rox Reference Dye II (Takara), 1 μl of the reverse transcription reaction, and 0.25 μM of each gene-specific primer in a final volume of 25 μl on a PRISM 7500 PCR instrument (Applied Biosystems). The reactions were performed at 95◦C for 10 s, 45 cycles of 95◦C for 5 s, and 60◦C for 40 s. Disassociation curve analysis was performed as follows: 95◦C for 15 s, 60◦C for 20 s, and 95◦C for 15 s. Data were collected using the ABI PRISM 7500 sequence detection system following the instruction manual. The relative expression levels were analyzed using the 2–--CT method (Livak and Schmittgen, 2001). The rice *ACTIN1* gene was used as the internal control. Accession numbers of genes analyzed in this study: *SDG711*: LOC\_Os06g16390; *SDG718*: LOC\_Os03g19480; *OsLF*: LOC\_Os05g46370; *Ehd1*: LOC\_Os10g32600; *RFT1*: LOC\_Os06g06300; *Hd1*: LOC\_Os06g16370; *Hd3a*: LOC\_ Os06g06320; *RID1*: LOC\_Os10g28330; *OsGI*: LOC\_Os01g08700; *Ghd7*: LOC\_Os07g15770; *OsMADS14*: LOC\_Os03g54160; *OsMADS15*: LOC\_Os07g01820; *OsMADS50*: LOC\_Os03g03070; *OsMADS51*: LOC\_Os01g69850. The primers for real-time PCR are listed in Table S1.

### **VECTOR CONSTRUCTION AND PLANT TRANSFORMATION**

For over-expression (OX) vector, the *SDG711* full-length cDNA was amplified from DJ leaf mRNA using primer set OXSDG711- F and OXSDG711-R, then inserted into the OX vector pU1301 under the control of the maize ubiquitin gene promoter within *Kpn*I sites (Sun and Zhou, 2008). For RNAi vectors, gene-specific sequences of *SDG711* and *SDG718*, spanning from 2352 bp to 2860 bp relative to the translation start site and from 918 bp to 1599 bp relative to the translation start site, respectively, were amplified from cDNA using primer sets RiSDG711-F/RiSDG711- R and RiSDG718-F/RiSDG718-R, then inserted into the RNAi vector pDS1301 (Huang et al., 2007). Sequence amplified using the primers sets are listed in Table S1. The constructs were transformed into DJ (SDG711 OX and RNAi) and ZH11 (*SDG718* RNAi) plants by *Agrobacterium tumefaciens* (strain *EHA105) –* mediated transformation as previously described (Huang et al., 2007).

### **WESTERN BLOTTING ANALYSIS**

For Western blot analysis, histone-enriched fractions were extracted from wild type (WT), mutant, and transgenic leaves as described previously (Huang et al.,2007). Antibodies used inWestern blot are: anti-H3K27me3 (07-449,Millipore), anti-H3K27me2 (ab24684, Abcam), anti-H3K27me1 (ab113671, Abcam), anti-H3 (ab1791, Abcam), anti-H3K4me3 (07-473, Millipore), anti-H3K4me2 (07-430, Millipore), anti-H3K4me1 (07-436, Millipore), anti-H3K9me3 (ab8898, Abcam), anti-H3K9me2 (07-441, Millipore), anti-H3K9me1 (ab9045, Abcam) H3K36me1 (ab9048, Abcam), anti-H3K36me2 (ab9049, Abcam), and anti-H3K36me3 (ab9050, Abcam). Anti-SDG711 was prepared by immunizing rabbits with SDG711 protein produced in *Escherichia coli* (in pET-28a vector) and purified with His-tag protein purification beads (V8550, GE Healthcare). The anti-serum was affinity-purified with protein-A agarose beads purchased from Millipore (16-157).

### **CHROMATIN IMMUNOPRECIPITATION (ChIP)**

The chromatin immunoprecipitation (ChIP) experiment was performed as described (Huang et al., 2007). Rice leaves (0.8–1.0 g fresh weight) were harvested at the end of the dark period of 8 week-old LD grown and 6 week-old SD grown plants and crosslinked in 1% formaldehyde under vacuum. Chromatin was extracted and fragmented to 200–750 bp by sonication, and ChIP was performed using the following antibodies: H3K27me3 (07- 449, Millipore), H3K4me3 (07-473, Millipore), and Anti-SDG711. The precipitated and input DNA samples were then analyzed by real-time PCR with gene-specific primers listed in Table S1. All assays were performed at least three times from two biological replicates.

### **RESULTS**

### **THE RICE E(z) HOMOLOGOUS GENES DISPLAY DISTINCT PHOTOPERIODIC EXPRESSION PATTERNS**

Phylogenetic analysis has previously shown that rice E(z) homologous genes, *SDG711* and *SDG718*, are closely related to *Arabidopsis* *CLF* and *SWN*, respectively (Luo et al., 2009; Figure S1). The two rice genes share 51% amino acid sequence identity, with much higher conservations in the catalytic SET domain and the protein interaction domains (SAND, Cys-rich). Analysis of expression of the genes by real-time PCR indicated that *SDG711* was widely expressed in different tissues/organs, whereas *SDG718* was more expressed in leaves than other tested tissues/organs (**Figure 1A**). Because *OsEMF2b* and PRC2-interacting PHD domain protein genes are involved in photoperiod regulation of flowering time in rice (Luo et al., 2009; Wang et al., 2013; Yang et al., 2013), we studied whether the rice E(z) genes were photoperiod-responsive. We analyzed mRNA isolated from leaves of 6 week-old SD (10 h light/14 h dark)- and 8 week-old LD (14 h light/10 h dark)- grown plant leaves harvested at intervals of 4 h during a 36 h period. The analysis revealed that *SDG711* mRNA levels were higher in LD than in SD, whereas that of *SDG718* were higher in SD than LD (**Figure 1B**).

### *SDG711* **IS INVOLVED IN LD REPRESSION OF RICE FLOWERING**

To study the function of *SDG711*, we produced transgenic rice plants in DJ which contains functional *Hd1* (Naranjo et al., 2014) to knockdown the gene by RNAi and to over-express the gene by using the maize ubiquitin gene promoter (Sun and Zhou, 2008). Analysis of the transgenic plants revealed several lines with reduced expression and six lines with increased expression of the gene (**Figure 2A**). To check the protein level of SDG711, we performed Western blot analysis of protein extracts from the leaves of the WT and transgenic plants grown in LD by using antibodies generated against *E. coli*-produced SDG711 protein. The analysis confirmed the OX and RNAi of the gene in the transgenic lines (**Figure 2B**). Three single copy T3 homozygous lines (offspring of single insertion T2 lines that did not segregate transgene-negative

### **FIGURE 1 | Expression profiles of** *SDG711* **and** *SDG718***. (A)** Real-time PCR detection of transcript levels of the genes in indicated organs or developmental stages. **(B)** transcript levels of the two genes in young leaves harvested at intervals of 4 h during a 36 h period from plants grown under long day (LD) or short day (SD) conditions. Dark periods are indicated by black bars. Bars = mean ± SD from two biological repeats. Values are shown as relative to ACTIN transcript levels.

individuals) per transgene were selected for phenotypic analysis. During vegetative growth the transgenic plants did not display any visible morphological defects but exhibited altered HD (or flowering time) compared to the WT plants. In LD (14 h light/10 h dark), the HD of the OX plants was largely delayed (20–25 days, *p* < 0.001, student's *t*-test), while that of the RNAi plants was significantly earlier (10–16 days, *p* < 0.001, student's *t*-test) than WT (**Figure 3**). However, in SD (10 h light/14 h dark) the HD of *SDG711* OX and RNAi plants was not significantly different from WT (**Figure 3**). These data suggested that *SDG711* may have a function to repress flowering in LD. In addition to the effect on flowering time, the OX plants produced a higher number of stamens (Figure S2). The pollen viability was reduced in both the OX and RNAi plants of *SDG711* (Figure S2). These observations suggested that *SDG711* might play a role in fertility. To confirm the transgenic results, we characterized a T-DNA mutant that had the insertion located in the 5 -UTR of the gene (Figures S3A,B). One single insertion was identified by Southern blotting (Figure S3C). Real-time PCR revealed a higher level of transcripts of the gene (Figure S3D). Northern blotting experiments revealed that the increased transcripts were about the same size as in WT (Figure S3E). The mutant displayed the same phenotype as that of the OX plants (**Figure 3**). After three backcrosses the phenotype co-segregated with the presence of the insertion, indicating that the mutation was a gain-of-function mutation.

To study the effect of *SDG711* OX or RNAi on histone methylation, we performed Western blotting analysis of histones isolated

from leaves of the gain-of-function mutant, an RNAi and an OX line. The analysis revealed that compared to WT the levels of H3K27me3 were lower in the RNAi line, but higher in the gainof-function mutant and OX plants (Figure S4), suggesting that *SDG711* was required for the overall H3K27me2/3 in rice.

were surveyed. Significance of differences between WT and transgenic

plants are indicated by (∗∗p < 0.01, Student's t-test).

The flowering time phenotype suggested that *SDG711* had a function to suppress flowering in LD while without a clear effect in SD. We therefore analyzed the expression of flowering regulatory genes in both pathways. mRNA were isolated from 8 week-old LD or 6 week-old SD plant leaves at intervals of 4 h during a 24 h period. In LD, the mRNA levels of LD flowering activators including *Ehd1*, *Hd3a*, *MADS14*, and *MADS15* were clearly decreased in the OX, but increased in the RNAi plants compared to WT (**Figure 4A**). By contrast, the transcript levels of the LD flowering repressor *Ghd7* were not clearly altered in the transgenic plants (Figure S5). The expression of *RID1* and *OsMADS50* also appeared unchanged in the transgenic plants (**Figure 4A**; Figure S5). Interestingly, *Hd1* that acts as a repressor of *Hd3a* in LD was highly induced by *SDG711* OX and repressed by *SDG711* RNAi (**Figure 4A**), while the expression level of *OsGI* that activates *Hd1* was not changed (Figure S5). Because E(z) proteins are supposed to be transcriptional repressors, the activation of *Hd1* in *SDG711* OX plants might be due to an indirect effect. Recent data have shown that OsLF, a bHLH protein that directly represses *Hd1* (Wang et al., 2013). Examination of *OsLF* transcripts in the transgenic plants indicated that *OsLF* was repressed by *SDG711* OX, but activated by *SDG711* RNAi in LD (**Figure 4A**). The expression of flowering regulatory genes supported the flowering time phenotype observed in the transgenic plants (**Figure 3**). Analysis of the gain-of-function mutant confirmed the above observed effects of *SDG711* OX on flowering time gene expression (Figure S6). In SD, the expression levels of the flowering genes were not clearly affected by *SDG711* OX or RNAi, except the peak expression of *OsLF* in the RNAi plants was about 4 h earlier than the WT (**Figure 4B**). These results suggested that *SDG711* is mostly involved in LD repression of rice flowering time.

**(A)** Key flowering genes expression in SDG711 transgenic plants compared to WT in LD. **(B)** Key flowering genes expression in SDG711 transgenic plants

from plants of 3 T3 lines grown for 8 week-old under LDs and 6 week-old under SDs. Bars = mean ± SD from three biological repeats.

To study whether SDG711 protein levels were differentially regulated in LD and SD, Western blotting analysis of protein extracts from WT, RNAi and OX plants grown in LD and SD harvested at 12:00 (mid day, after 7 h in light for LD condition and 5 h in light for SD condition) and 24:00 am (mid night, after 5 h in dark for LD condition and 7 h in dark for SD condition) was performed by using anti-SDG711. The results shown in **Figure 5** revealed that SDG711 protein levels were lower in SD compared in LD and that the high SDG711 level of the OX line was reduced in SD and in the dark, suggesting that the stability of SDG711 protein may be regulated by day length and light/dark conditions.

### *SDG718* **PROMOTES FLOWERING IN SD**

To study whether *SDG718* had a function in flowering time control, we obtained several RNAi lines in ZH11 variety (**Figure 2C**), which also has a functional Hd1 (Luan et al., 2009; Naranjo et al., 2014). These plants showed a similar HD as the WT in LD, but flowered later (>20 days) than WT in SD (**Figure 3B**). Analysis

of flowering time gene expression revealed that *SDG718* RNAi clearly induced *OsLF*, but repressed the SD activators including *Hd1*, *Hd3a*, *OsMADS14,* and *RFT1* (**Figure 4C**), suggesting that *SDG718* played a role in promoting flowering in SD. However, the *SDG718* RNAi plants did not show any clear defect in flower organ and pollen viability (Figure S2).

### *SDG711* **AND** *SDG718* **REGULATE H3K27me3 ON FLOWERING GENES**

To study whether knockdown and OX of *SDG711* altered chromatin modification of flowering regulatory genes, we analyzed H3K27me3 on *RID1*, *Ehd1*, *RFT1*, *Hd3a*, *MADS14*, *MADS15*, *Hd1*, and *OsLF* in WT, OX and RNAi plants grown under LD conditions by ChIP assays. The 8 week-old plant leaves were harvested at the end of the darkness period under LD conditions. Because H3K27me3 is located in gene body with enrichment on the 5 end of the gene in rice (Hu et al., 2012; Li et al., 2013), we analyzed the ChIP by real-time PCR using two primer sets, one corresponding to the 5 region, the other to the gene body (**Figure 6**). Among these genes, *Hd1* appeared to be unmodified by H3K27me3 in either the 5 region or the gene body, further supporting an indirect effect of *SDG711* on its expression. *RID* displayed a very low level of H3K27me3 (**Figure 6A**), consistent with the observation that the expression of the gene was not affected in the *SDG711* transgenic and the gain-of-function mutant plants (**Figure 4**; Figures S5 and S6). The remaining genes displayed H3K27me3 in the 5 region and/or gene body, suggesting that regulation of these genes might involve PRC2 function. H3K27me3 on these genes (except the gene body of *RFT1* which displayed a relatively weak level of H3K27me3) was clearly reduced in the RNAi but increased in the OX plants (**Figure 6A**), which reversely correlated with their expression changes in the transgenic plants. These data suggested that SDG711-mediated H3K27me3 was involved in the regulation of the genes.

Because H3K27me3 is antagonistic to H3K4me3 on gene activity, we therefore analyzed whether alteration of H3K27me3 affected H3K4me3 on the flowering genes in the transgenic plants.

performed, one repetition is shown. One other biological replicate is shown in Figure S7. Bars = mean ± SD from three technical repeats. **(C)** ChIP analysis of H3K27me3 on the indicated flowering genes in WT and SDG718 RNAi plants pooled from 3 T3 lines. Leaves of SDG711 transgenic plants were harvested from 8 week-old LD plants grown, and leaves of SDG718 transgenic plants were harvested from 6 week-old SD plants grown.

The analysis revealed relatively higher levels of H3K4me3 in *Hd1*, *Hd3a*, *MADS15*, *Ehd1,* and *OsLF* than in *MADS14*, *RFT1,* and *RID* in WT plants (**Figure 6B**). H3K4me3 levels in *Hd1*, *MADS15,* and *Hd3a* were not affected by the *SDG711* transgenes except some increases on *Hd1* in the OX plants (**Figure 6B**). The increased H3K4me3 on *Hd1* may be due to increased expression of the gene in the OX plants, as H3K4me3 is thought to be associated with active genes (Hu et al., 2011). However, H3K4me3 on *Ehd1* and *OsLF* was increased in *SDG711* RNAi but decreased in *SDG711* OX plants (**Figure 6B**), which conversely correlated to that of H3K27me3 and suggested that *SDG711*-mediated H3K27me3 might affect H3K4me3 in the two loci, which may be linked to the repression of the genes. Analysis of *SDG718* RNAi plants also revealed a clear decrease of H3K27me3 in *OsLF* in SD (**Figure 6C**).

### **SDG711 BINDS TO Ehd1 AND OsLF LOCI**

To further assess the function of *SDG711* on the flowering gene regulation, we performed anti-SDG711 ChIP assays and analyzed by real-time PCR using the same primer sets as for the histone methylation ChIP. Non-immunized serum was used as control. The analysis revealed that in LD-grown WT plants SDG711 binding was clearly enriched in the 5 end and the gene body of *Ehd1* and *OsLF* compared to the other tested genes (**Figure 7A**). However, some enrichment was also observed in the gene body of *RID1*. The SDG711 enrichment on the three genes was sensibly reduced in RNAi and clearly enhanced in the OX plants (**Figure 7A**), suggesting that SDG711 may directly target to the genes in LD. In SD-grown plants, the enrichment of SDG711 on the three genes was much weaker than that in SD (**Figure 7B**), which was consistent with lower accumulation of SDG711 in SD (**Figure 5**) and suggested that SDG711 enrichment in the flowering genes was regulated by day length.

### **DISCUSSION**

Our data suggest that *SDG711* and *SDG718* repress *OsLF* in LD and SD, respectively, leading to the activation of *Hd1* that inhibits *Hd3a* and flowering in LD but activates *Hd3a* and *RFT1* and flowering in SD (**Figure 8**). Therefore, the two E(z) genes likely contribute to the accurate photoperiod control of flowering in rice. The increased level of *SDG711* expression in LD and that of *SDG718* in SD may be critical for *OsLF* repression. The SD-induced *SDG718*

expression is consistent with the results showing that *OsVIL* genes that promote flowering in SD are also induced in SD (Wang et al., 2013). The *SDG711*-mediated LD repression of flowering is supported by previous results showing that mutants of other PRC2

**genes. (A)** SDG711 protein enrichment on the indicated flowering genes in LD-grown WT, SDG711 RNAi and OX seedlings tested by ChIP with anti-SDG711. Non-immunized rabbit serum (IgG) was used as control. Three biological repeats were performed, one repetition is

shown. The two other repeats are shown in Figure S8. Bars = means ± SD from three technical repeats. **(B)** SDG711 protein enrichment on the indicated flowering genes in SD-grown plants. Leaves were harvested from plants grown for 8 week-old LD and 6 week-old SD plants grown.

genes such as *osfie2* and *osemf2b* also display an early flowering phenotype in LD (Luo et al., 2009). The observations that the *SDG711* OX only enhanced LD repression of flowering but did not affect flowering in SD (**Figure 3**), suggest that *SDG711* may be also regulated at posttranscriptional levels by day length. This hypothesis is supported by the observation that accumulation of SDG711 protein in the OX lines in LD was reduced in SD and that SDG711 binding to key flowering genes was largely reduced in SD (**Figure 7B**). This day length-dependent stability of the protein may allow SDG711 to regulate flowering only in LD. The changes of H3K27me3 levels on the *OsLF* locus in *SDG711* transgenic plants imply that PRC2-mediated repression may involve the deposition of H3K27me3 on the gene. As SDG711 was shown to be directly associated with the *OsLF* locus in LD (**Figure 7**), *OsLF* may be a primary target for rice PRC2-mediated LD-repression of flowering. The results are in agreement with recent data showing directly repression of *OsLF* by LC2 that interacts with the OsEMF2b protein (Wang et al., 2013; Yang et al., 2013). However, it is not excluded that *SDG711* or *SDG718* represses other genes that mediate *Hd1* repression.

Previous studies have shown that the flowering promoter *Ehd1* and the flowering repressor *Ghd7* could enable manipulation of slight differences in day lengths to control *Hd3a* transcription with a critical day length threshold (Itoh et al., 2010). *Ehd1* and *Ghd7* provide a gating mechanism to set critical day length for *Hd3a* expression in SD, in which *Ehd1* is repressed by the morning activation of *Ghd7* in LD (Itoh et al., 2010; Tsuji et al., 2013). The observations that *SDG711* transgenic and gain-of-function mutant plants affected the expression of *Ehd1* but not that of its upstream regulators (i.e., *Ghd7*, *MADS50* or *RID1*) and that SDG711 was associated with the *Ehd1* locus (**Figure 6**; Figures S5 and S6), indicate that SDG711 is a direct repressor of *Ehd1* in LD. In addition, H3K27me3 on the *Ehd1* locus was modulated by the expression level of *SDG711* (**Figure 6A**). Therefore, besides *Ghd7* and *OsMADS50*, *SDG711* represents an additional LD repressor of *Ehd1* (**Figure 8**). The observations that *Ehd1* expression was not clearly affected by down-regulation of *SDG718* and *SDG711* or up-regulation of *SDG711* in SD (**Figures 4B,C**), suggest that PRC2 may mainly target to the OsLF-Hd1 pathway instead of that of *Ehd1* in SD. Alternatively, the SD-induced expression of *Ehd1* (mainly due to the repression of *Ghd7* in SD) may be overwhelming, which may mask the effect of *SDG718* RNAi. In addition, the chromatin analysis data indicated that key flowering genes displayed different histone modification patterns (**Figure 6**). The changes of H3K27me3 levels on the marked genes caused by *SDG711* RNAi and OX and *SDG718* RNAi, which are correlated with their expression change in the transgenic plants in LD or SD, suggest that *SDG711*/*SDG718* may be also involved in the deposition of the mark on these loci. Collectively, the data demonstrating that the *SDG711* and *SDG718* are involved respectively, in the LD and SD signaling to promote LD repression and SD activation of flowering, suggest the involvement of PRC2 in the accurate photoperiod control of flowering in rice.

*SDG711* and *SDG718* are closely related to *Arabidopsis CLF* and *SWN*, respectively. It is suggested that *CLF* and *SWN* act redundantly to regulate vegetative growth. No vegetative phenotype observed in *SDG711* and *SDG718* transgenic plants support the idea that the two rice genes may also have a redundant function in the vegetative growth. *Arabidopsis* genome contains a third E(z) gene, *MEA*, which is mainly involved in the regulation of gene imprinting and embryo and endosperm development. The counterpart of *MEA* is not found in rice. Therefore the question arises as that whether either *SDG711* or *SDG718* plays a role in reproductive development. Our data showing that besides flowering time, *SDG711* RNAi and OX affects stamen number, pollen viability, and fertility (Figure S2) suggests that this E(z) gene may be involved in reproductive development in addition to flowering time control in rice.

### **AUTHOR CONTRIBUTIONS**

Xiaoyun Liu and Dao-Xiu Zhou designed the experiments and wrote the manuscript. Xiaoyun Liu, Chao Zhou, and Yu Zhao performed the experiments. Shaoli Zhou performed some images processing. Wentao Wang performed some materials cultivation. All authors have read and approved the manuscript.

### **ACKNOWLEDGMENTS**

We thank Prof. Caiguo Xu and Prof. Xianghua Li for help in field experiments and management. This work was supported by grants from National Science Foundation of China (31371241 and 31371468), the 863 Key Project "Rice Functional Genomics" of the Chinese Ministry of Science grants (2012AA10A303), the special transgenic program of the Chinese Ministry of Agriculture (2014ZX0800902B), and the Melinda and Bill Gates Foundation.

### **SUPPLEMENTARY MATERIAL**

The Supplementary Material for this article can be found online at: http://www.frontiersin.org/journal/10.3389/fpls.2014.00591/ abstract

### **REFERENCES**

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

### *Received: 02 September 2014; accepted: 10 October 2014; published online: 31 October 2014.*

*Citation: Liu X, Zhou C, Zhao Y, Zhou S, Wang W and Zhou D-X (2014) The rice enhancer of zeste [E(z)] genes SDG711 and SDG718 are respectively involved in long day and short day signaling to mediate the accurate photoperiod control of flowering time. Front. Plant Sci. 5:591. doi: 10.3389/fpls.2014.00591*

*This article was submitted to Plant Biotechnology, a section of the journal Frontiers in Plant Science.*

*Copyright © 2014 Liu, Zhou, Zhao, Zhou, Wang and Zhou. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited andthatthe original publication inthis journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# An efficient method for miRNA detection and localization in crop plants

### *Flor de Fátima Rosas-Cárdenas1,2 , Rocío Escobar-Guzmán1, Andrés Cruz-Hernández 3 , Nayelli Marsch-Martínez 4 and Stefan de Folter 1\**

1 Laboratorio Nacional de Genómica para la Biodiversidad, Centro de Investigación y de Estudios Avanzados del Instituto Politécnico Nacional, Guanajuato, México

2 Centro de Investigación en Biotecnología Aplicada del Instituto Politécnico Nacional, Tlaxcala, México

3 Facultad de Ingeniería, Universidad Autónoma de Querétaro, Querétaro, México

4 Departamento de Biotecnología y Bioquímica, Centro de Investigación y de Estudios Avanzados del Instituto Politécnico Nacional, Guanajuato, México

### *Edited by:*

Clelia De-la-Peña, Centro de Investigación Científica de Yucatán, Mexico

### *Reviewed by:*

Serena Varotto, University of Padova, Italy Norberto Daniel Iusem, Universidad de Buenos Aires, Argentina

### *\*Correspondence:*

Stefan de Folter, Laboratorio Nacional de Genómica para la Biodiversidad, Centro de Investigación y de Estudios Avanzados del Instituto Politécnico Nacional, Km. 9.6 Libramiento Norte, Carretera Irapuato-León, CP 36821 Irapuato, Guanajuato, México e-mail: sdfolter@langebio. cinvestav.mx

microRNAs are a class of non-coding small RNAs (sRNAs) that are important regulators of gene expression at the post-transcriptional level by mRNA cleavage or translation inhibition. Another class of sRNAs are siRNAs, which also regulate gene expression but by causing DNA methylation. This epigenetic regulatory role has been observed for some miRNAs as well.The use of sRNAs allows the development of biotechnological applications in plants.To develop these types of applications, and to better understand the natural roles they play, it is important to be able to detect and to localize these sRNAs at the plant tissue level. Sometimes, in crop plants this can be challenging.Therefore, we developed a tissue printing hybridization protocol for easy and efficient detection of sRNAs and demonstrate this by the analysis of the spatio-temporal expression patterns of the miRNAs miR159 and miR164 in fruits of various crop plants. Moreover, we show the possibility to also detect the expression of miRNAs in fruit juice using a dot blot hybridization approach.

**Keywords: tissue-printing, miRNA detection, hybridization, fruits, miR164**

### **INTRODUCTION**

microRNAs (miRNAs) are a class of small RNAs (sRNAs) that are fundamental regulatory elements of eukaryotic genomes (Voinnet, 2009) and their widespread conservation and divergence in the plant kingdom has been demonstrated (Chavez Montes et al., 2014). Detection of the spatio-temporal expression of miRNAs is critical to understand their function (van Rooij, 2011). Several techniques such as Northern-blot hybridization (Pall and Hamilton, 2008), microarray analysis (Yin et al., 2008), stem-loop RT-PCR analysis (Chen et al., 2005), or sRNA-seq analysis (Chavez Montes et al., 2014) permit the identification of miRNA expression patterns, which may suggest their involvement in certain biological processes. Besides these advantages, a drawback of these techniques is the lack of information about the actual localization of the miRNA in the tissue or organ itself, which is important for understanding the biological function of the miRNA.

A technique overcoming this drawback is *in situ* hybridization, which also allows the detection of miRNAs in thin tissue sections using a labeled complementary probe against the miRNA of interest (Kidner and Timmermans, 2006;Várallyay and Havelda, 2011). Nevertheless, the generation of thin tissue sections from large tissues using paraffin-embedding or cryosections may be challenging. Tissue printing in nitrocellulose or nylon membranes is a technique employed to study the localization of proteins, nucleic acids, and soluble metabolites from freshly cut tissue slices. Tissue printing has been defined as 'the art and science of visualizing cellular material and information that are transferred to