Primer Pairs for Rice (Oryza sativa L.) Bisulfi te Sequencing Studies

Many agriculturally important properties such as heterosis, inbreeding depression, phenotypic plasticity, and resistance for biotic and abiotic stresses are thought to be affected with epigenetic components. New discoveries related with epigenetics are likely to have a major impact on strategies for crop improvement in rice breeding. However, assessing the contribution of epigenetics to heritable variation in plant species still poses major challenges. Methylation of cytosine in DNA is one of the most important epigenetic mechanisms in plants. DNA methylation not only plays signifi cant roles in the regulation of gene activity, but also it is related with genomic integrity. Although most of next generation sequencing (NGS) technologies do not require the use of target specifi c primer pairs to identify and study DNA cytosine methylation, validation studies of NGS uses selective primer pairs. Bisulfi te sequencing technique is a gold method for DNA cytosine methylation studies. However, bisulfi te sequencing requires the development of bisulfi te primer pairs to selectively study DNA sequences of interest. In this study 9 bisulfi te specifi c primer pairs were identifi ed and validated. These primer pairs successfully amplifi ed bisulfi te converted and unconverted genomic DNA extracted from radicle and plumule of rice (Oryza sativa L.) seedlings. Results of the present study clearly revealed the occurrence of CG, CHG and CHH (H stands for C. T or A nucleotides) contents in studied DNA sequence targets were different indicating potential role of DNA cytosine methylation in these genes. Primer pairs reported in this study could be used to detect DNA methylation which is one of the most important epigenetic mechanisms affecting the development, differentiation or the response to biotic and abiotic stress in rice (Oryza sativa L.). Research Article Primer Pairs for Rice (Oryza sativa L.) Bisulfi te Sequencing Studies Mehmet Karaca1* and Ayse Gul Ince2 1Department of Field Crops, Akdeniz University, 07070 Antalya, Turkey 2Vocational School of Technical Sciences, Akdeniz University, 07070 Antalya, Turkey *Address for Correspondence: Mehmet Karaca, Department of Field Crops, Akdeniz University, 07070 Antalya, Turkey, Email: mkaraca@akdeniz.edu.tr Submitted: 22 October 2018 Approved: 05 November 2018 Published: 06 November 2018 Copyright: © 2018 Karaca M, et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited


Introduction
Paddy rice (Oryza sativa L.) is grown in more than one hundred countries in the world with a total harvested area of approximately 158 million hectares, producing more than 700 million tons, 470 million tons of which are milled rice. Nearly 640 million tons of paddy rice are grown in Asia, followed by Sub-Saharan Africa that produces about 19 million tons and Latin America that produces about 25 million tons. Currently China is the world's largest rice producer, accounting for 35% of total world's rice production with a planting area of 29 million hectares [1].
Turkey is not among the larger rice producing countries. The Marmara-Thrace region of Turkey is the largest rice-growing area and makes the greatest contribution to production, followed by the Black Sea region. Paddy rice (rice in the husk) production of Turkey is about 750000-850000 tons from approximately 95000-100000 hectares planting area [2]. In Turkey, average yield of paddy rice per hectare ranges from 6.5 tones (Baldo variety) to 10 tones (Luna variety) while average milling rate ranges from 55% (Baldo variety) to 68% (Luna variety). Edirne province in the Thrace region produces more than 50 percent of Turkey's paddy rice production while Central Anatolia ("Çankırı" and "Çorum" regions) produces nearly a quarter of paddy rice production. Rice varieties such as Osmancık 97, Baldo, Cammeo, Galileo, Ronaldo, and Luare have popularity in Turkey because of their high yields and high milling rates along with their better cooking and lavor characteristics. Rice in Turkey is irreplaceable for making pilaf, rolling stuffed, stuffed vegetables and several desserts such as rice pudding. However, rice production of the Turkey is not suf icient to counteract its consumption which is about 790000 tons annually [3].
In many parts of the world including Turkey low and high temperatures limit and damage rice growth during the different stages of plant development. Low temperatures during germination and at the seedling negatively affect panicle formation, lowering and pollination, resulting in yield losses. On the other hands, high temperatures in many parts of rice production areas in the world especially in the southeastern part of world cause spikelet sterility. Also drought or water shortage, salinity and alkalinity negatively affect the rice growth and development. In Turkey, blast disease (Pyricularia oryzae), brown leaf spot (Helminthosporium oryzae) and root rot (Fusarium moniliforme) are among the most important rice diseases [2,3].
Many agriculturally important properties such as heterosis, inbreeding depression, phenotypic plasticity, and biotic and abiotic stresses are thought to be affected and regulated with epigenetic components. New epigenetic discoveries are likely to have a major impact on strategies for crop improvement in plant breeding. Most abiotic and biotic growth limiting factors are known to be under epigenetic regulations [4]. Aberrant DNA methylation is known associated with loss of DNA homeostasis and genomic instability leading to development of human diseases such as cancer. Epigenetic variation due to the DNA cytosine methylation and histone levels in the genome results from enzymemediated chemical modi ications. Methylation of cytosine in DNA is one of the most studied and important epigenetic modi ications in plants. In plant genomes, DNA cytosine methylation occurs at the C5 position of cytosines in the three sequence contexts or islets, CG, CHG and CHH where H is A, T, or C. In plants, Domains Rearranged Methyltransferase 2 (DRM2) plays major roles in de novo cytosine methylation. There are speci ic enzymes for methylation and maintenance among the three sequence contexts CG, CHG and CHH. For instance, DNA methyltransferase 1 (MET1) and chromomethylase 3 (CMT3) catalyze methylation of CG and CHG while CHH methylation is maintained through de novo methylation by DRM2 [5,6]. It is well known that DNA methylation not only plays signi icant roles in the regulation of gene activity, but it is also related with maintenance of genomic integrity. Therefore, the studies of DNA cytosine methylation offer important inding for gene expression and regulation [7,8].
Currently there exist several techniques for assaying and detecting global and genespeci ic cytosine DNA methylation in plant species. Among these techniques bisul ite sequencing is the gold standard for assessing DNA methylation. Bisul ite sequencing uses bisul ite-mediated deamination before sequencing. It could reveal the methylation status of every cytosine nucleotide of a genome under study. It was developed by Frommer and colleagues and recognized as a revolution in DNA methylation analysis. It is based on conversion of genomic DNA by using sodium bisul ite application [8]. This method takes the advantage of the selective and complete conversion of unmethylated cytosine to uracil by sodium bisul ite application. In this method converted cytosine nucleotides are ampli ied using polymerase chain reactions (PCR) as thymine nucleotides [7]. PCR ampli ied products are cloned (or directly sequenced) and, the sequencing experiments reveal the initial methylation pro ile of the region of interest [8,9].
Reliability and success of bisul ite-mediated deamination and sequencing method is strongly dependent on the availability of primer pairs that could uniquely and successfully amplify bisul ite treated and untreated genomic DNA region. By selecting speci ic primer pair any gen structure and regulatory sequences could be studied. The aim of this study was to develop a set of bisul ite speci ic primer pairs using DNA sequences of rice deposited in NCBI GenBank database and validate these primers on rice cultivar Osmancik 97.

DNA Sequences utilized
Genomic DNA sequences of 8 different rice genes (Oryza sativa L.) deposited in NCBI GenBank databases (https://www.ncbi.nih.gov/) were downloaded. Downloaded DNA sequences were used as input iles for EpiOne software [8]. The occurrence of methylation context, CG, CHG and CHH were calculated and occurrences and densities were identi ied for each gene elements.
Sequences obtained from EpiOne were used to design primer pairs speci ic to bisul ite sequencing studies using Primer 3 software [10] based on the following main parameters: GC content value of a primer was set between 40% and 80%, annealing temperature (Tm) was set between 58°C and 62°C, and expected ampli ied product size was de ined as 300-800 bp. After designing the sequences, cytosine bases (C) in the forward primers were replaced with Y (C/T) and the guanine bases (G) in the reverse primers were replaced with R (A/G). In the primer sequences R stands for purine (A or G) and Y stands for pyrimidine (C or T).

Genomic dna isolation and bisulfi te treatment experiments
Genomic DNA samples of rice cultivar Osmancik 97 were extracted from plumule and radicle sections of seedlings using a protocol previously described in Karaca et al. [11]. Genomic DNA samples obtained from plumule and radicle sections were bisul ite treated using a bisul ite conversion kit (Invitrogen Corp., Carlsbad, CA, USA) as previously described in Karaca et al. [12].
The Td-PCR ampli ication pro ile used in the present study was as follows: initial denaturation at 94°C for 3 minutes, ten cycles with denaturation at 94°C for 30 seconds, annealing at 60°C (or 66°C) for 30 seconds in the irst cycle, diminishing by 0.5°C each cycle, and extension reactions at 72°C for 1 minute using a Veriti 96-well thermal cycler (Applied Biosystems, Foster City, CA, USA). Additional 40 more PCR cycles were run using the same cycling parameters mentioned above with constant annealing at 55°C (or 60°C). Denaturation and extension conditions were the same as indicated above. The ampli ication reactions ended with inal extension reactions at 72°C for 10 minutes PCR products in the presence of 1x DNA loading buffer were loaded in 3% (w/v) high-resolution agarose gels (SERVA Electrophoresis GmbH, Heidelberg, Germany) and electrophoresed at 5 V per cm at constant voltage for 4-6 hours [8].

Electrophoretic purifi cation, cloning, purifi cation and sequencing experiments
PCR products of bisul ite treated and untreated targets (control or reference DNA) ampli ied with four different primer pairs were added with 5 μL 6x DNA loading buffer, mixed well and loaded in 3% (w/v) high-resolution agarose gels and electrophoresed at 5 V per cm of gel at constant voltage for 4-6 hours. PCR products in the gels were transferred into 2 mL tubes and puri ied using a PureLink Quick Gel Extraction & PCR Puri ication Combo Kit (Invitrogen Corp. Carlsbad, CA, USA). Puri ied amplicons were ligated into pTZ57R/T, vectors and were transformed to E. coli bacteria strain JM107 using a Transform Aid Bacterial Transformation Kit (Thermo Scienti ic, Waltham, MA USA). Targets carrying colonies were selected and sub-cultured and plasmid DNA samples were extracted using a GeneJET Plasmid Miniprep Kit (Thermo Scienti ic). A total of 12 plasmids with PCR fragments from bisul ite treated genomic DNA and 8 plasmids with PCR fragments from untreated DNA were sequenced using M13R sequencing primers (Macrogen Inc., Amsterdam, The Netherland) [12].
Sequences of plasmids were downloaded from Macrogen Inc. and primer sequences of forward and reverse were identi ied from each clone sequence and trimmed off along with the vector sequences. Target sequences were formatted into fasta ile using Sequencher program and analyzed using the KisMeth software program [13]. The methylation status of each target was calculated for each cytosine sequence context (CG, CHG and CHH) using the percentage methylation (%), calculated as 100×C/(C+T). The statistical signi icance between and within the samples of the entire set of either CG, CHG or CHH sites was separately evaluated with the nonparametric Mann-Whitney U test. In this study, two-tailed p values, 0.05 were considered statistically signi icant. Methylation context (CG, CHG or CHH) was used as factor and the methylation percentage was use as response.

Results and Discussion
In this study, a total of seven genes and regulatory sequences with different functions were studied. These genes were further mined into gen bodies and regulatory sequences using EpiOne software program [8]. Based on the EpiOne programs, sequences were mined into gene structures or elements such as exons, introns, promoters and untranslated regions (UTRs). Also three Coding DNA Sequences (CDS) that were portions of mRNA transcripts translated by ribosomes were also studied ( Table 1). Introns and exons are usually called as gene bodies that do not contain regulatory regions such as promoters and UTRs. Studies of DNA methylation in gene bodies provide further advantages [14] since type, amount, and location of DNA cytosine methylation within a gene provides signi icant knowledge on the enzymatic pathway.
Occurrence of GC, GHG and GHH islets among sucrose phosphate synthase, 3-hydroxy-3-methylglutaryl-CoA reductase, BP-73 gene for SAP-like protein, tub A2 gene for tubulin alpha, promoter region and 5'-UTR of CDPK2, tub6 gene for beta-tubulin 6 exon 1-2 and granule-bound starch synthase (Wx) with 5'-UTR were different as calculated with EpiOne program. This indicated that occurrence and densities of these methylation speci ic islets within gene bodies and regulatory regions could play signi icant roles in regulation of these genes studied in rice gene expression [12][13][14].
We could design more primer for some genes while some had fewer number of primers (Table 1). We strictly avoid CG dinucleotides, potential hairpin structures and possible primer dimers. Ampli ication of targets used touch-down PCR conditions. This help us not to rely on a gradient PCR thermocycler to determine the appropriate annealing temperature. Primer pairs from intron, exons and regulatory sites were designed so that DNA methylation could be targeted at speci ic regions on a gene element (either exon, intron or UTRs). In the present study, we also used a strategy similar to that of primer walking to cover whole intronic or exonic regions for a gene under study. Gene targets used in this study are among genes that are important for rice yield and productivity. For instance, sucrose synthase gene is very important for sucrose which is the major photosynthetic product in many plant species. Also sucrose is main sugar for translocation of carbohydrate form leaves to tissues such as grain in rice. Furthermore, sucrose is also known a signal molecule that regulates gene expression in plants [15]. Therefore, information regarding the epigenetic and genetic regulation of sucrose synthase gene could be used in the studies of rice improvement plant protection and productivity.
Target genes also included 3-hydroxy-3-methylglutaryl-CoA reductase (rice HMGR like) and granule-bound starch synthase (Wx) gene. The gene HMGR like coded in rice has been found to increase expression not by wounding but pathogen, Magnaporthe grisea, a well-known causal agent of rice blast disease [16]. This suggested that this enzyme in rice might play a key role as a component of the inducible defense mechanism, therefore, epigenetic studies of this gene may provide information in understanding of epigenetic regulation. In rice as well as in other cereals, granulebound starch synthase gene is responsible for amylose synthesis. This enzyme is regulated at the transcriptional and post-transcriptional levels but its epigenetic regulation is not known. Since this enzyme affects the ratio of amylose to amylopectin in grain. Understanding the regulation of granule-bound starch synthase activity may lead to the improvement of rice eating quality [17].
A total of 18 bisul ite speci ic primer pairs were initially designed and they were irst tested using rice genomic DNA sample that was not bisul ite treated (called control or reference DNA). Genomic DNA samples were ampli ied using two different touch-down PCR pro iles. Primer pairs; OS-FR-11 and OS-FR-13 both targeting intronic and exonic region in BP-73 gene for SAP-like protein were failed to amplify targets in both touch-down ampli ication pro iles (Ann 60 and Ann 66). On the other hand, OS-FR-14 which was also failed to amplify target, could amplify target when annealing temperature reduced to 50 o C. Primer pairs OS-FR-15, OS-FR-19 produced an extra amplicon at Ann 60, but extra amplicons were vanished away when these primer pairs were used at Ann 66 ampli ication pro ile.
Among 16 primer pairs, which were found to successfully ampli ied genomic DNA samples of variety Osmancik 97, those OS-FR-01, OS-FR-02, OS-FR-07, OS-FR-09, OS-   Table 2. On the other hands, 7 primer pairs failed ampli ication of bisul ite treated genomic DNA probably due to the original C richness creates long stretches of uracil, which might result in DNA fragmentation during the bisul ite-treatment. We assumed that fragmentation of targets were encountered in those targets greater lengths. This is why bisul ite sequencing studies use short segments of genomic region in order to avoid failed ampli ication due to the bisul ite treatments. During the bisul ite treatment DNA fragmentation occurs and this may reduce or remove the number of available molecules for PCR ampli ication [12].
In order to validate bisul ite speci ic primer pairs developed in this study, primer pairs OS-FR-02 and OS-FR-10 were used. OS-FR-02 ampli ied exonic region of sucrose phosphate synthase gene while OS-FR-10 ampli ied intronic region of 3-hydroxy-3methylglutaryl-CoA reductase gene ( Figure 1). As shown in Figure 1, DNA cytosine methylation in the intronic region of 3-hydroxy-3-methylglutaryl-CoA reductase gene was not statistically different between radicle and plumule of rice variety Osmancik 97. On the other hands, methylation of CHH were different between radicle and plumule of rice variety Osmancik 97. Methylation of CG and CHG were not statistically differed between radicle and plumule of rice variety Osmancik 97. Based on the high level of occurrence and densities of methylation islets in 7 gene regions that we did not validate in the present study we expected that their use in DNA methylation studies in rice would reveal valuable information. Also primer pairs reported in this study could be used as molecular markers in genetic and breeding studies [18].
In the present study selected genes are important genes for rice yield and yield related parameters. Previous research identi ied a number of epigenetic regulators in rice [19]. For instance, it was found that CG and CHG hypomethylation associated with activation of genes preferentially expressed in rice seed development. Other epigenetic regulators in rice include involve in a wide range of developmental and stress responsive pathways. Epigenetic variations detected with the use of primer pair reported in the present study could be used to identify epialleles related to important agronomic traits and may provide new strategies to enhance agronomically favorable traits and grain productivity in rice [15,19].

Conclusion
Increased yield and biotic and abiotic resistance in rice is one of the most desired research goals. Yield of rice can be determined using data of panicle number per plant, the illed grain number per panicle, and the weight per grain. Genetic and epigenetic regulation in these grain related characteristics are important for increased rice yield. Results revealed that methylation islets in the contexts of CG, CHG and CHH among 7 genes or regulatory regions consisting of sucrose phosphate synthase, 3-hydroxy-3methylglutaryl-CoA reductase, BP-73 gene for SAP-like protein, tub A2 gene for tubulin alpha, promoter region and 5'-UTR of CDPK2, tub6 gene for beta-tubulin 6 exon 1-2 and granule-bound starch synthase (Wx) with 5'-UTR were different. We reported 9 gene speci ic primer pairs suitable in epigenetic studies in Osmancik 97 which is one of the most sown rice cultivars in Turkey. The use of these primer pairs could provide understanding of epigenetic regulation in rice, with an emphasis on the roles of key epigenetic regulators, the epigenomic landscape, epigenetic variation, transposon repression, and plant development. These primer pairs reported in the present study could be also used in the validation studies in high-throughput sequencing of rice genomes. Figure 1: Methylation status of 3-hydroxy-3-methylglutaryl-CoA reductase gene between radicle and plumule tissues of rice cultivar Osmancik 97. Each circle in the panels represents a cytosine residue, methylated cytosine residues are shown with full circles, and non-methylated cytosine residues are shown with empty circles. Cytosine residues are color-coded by their sequence context: red for CG, blue for CHG and green for CHH where H is C, T or A. Each row of panels represents a sequence of clone.