Abstract
Background
Epigenetic modifications play important roles in the regulation of
plant development. DNA methylation is an important epigenetic
modification that dynamically regulates gene expression during
developmental processes. However, little studies have been reported
about the methylation profiles of photoperiod- and thermo-sensitive
genic male sterile (PTGMS) rice during the fertility transition.
Results
In this study, using methylated DNA immunoprecipitation sequencing
(MeDIP-seq), the global DNA methylation patterns were compared in the
rice PTGMS line PA64S under two different environments (different
temperatures and day lengths). The profiling of the DNA methylation
under two different phenotypes (sterility and fertility) revealed that
hypermethylation was observed in PA64S (sterility), and 1258
differentially methylated regions (DMRs) were found between PA64S
(sterility) and PA64S (fertility). Twenty differentially methylated
genes of them were further validated through bisulfite sequencing, and
four of these genes were analyzed by qRT-PCR. Especially, a
differentially methylated gene (LOC_Os08g38210), which encoded
transcription factor BIM2, is a component of brassinosteroid signaling
in rice. The hypermethylated BIM2 gene may suppress some downstream
genes in brassinosteroid signaling pathway, and thus affect the male
fertility in PA64S.
Conclusions
The results presented here indicated that hypermethylation was observed
in PA64S (sterility). Gene Ontology (GO) analysis and KEGG analysis
revealed that flavone and flavonol biosynthrsis, circadian rhythm,
photosynthesis and oxidative phosphorylation pathways were involved in
sterility-fertility transition of PA64S.
Electronic supplementary material
The online version of this article (doi:10.1186/s12864-015-1317-7)
contains supplementary material, which is available to authorized
users.
Keywords: Oryza sativa, DNA methylation, PTGMS, Fertility transition,
Gene expression
Background
DNA methylation is an indispensable epigenetic modification of plant
genomes. It occurs predominantly in the CG context in both plants and
animals [[31]1-[32]3]. The state of CpG methylation regulates and
stabilizes chromatin structure, and possibly regulates the
accessibility of these DNA regions to the transcriptional machinery
[[33]4]. Studies have shown that DNA methylation plays an important
role in many plant processes, including transposon silencing, virus
defence, and gene imprinting [[34]5-[35]7].
In plant, DNA methylation is mainly distributed in heterochromatic
regions that are composed of tandem or inverted repeats, i.e.,
transposons. Transposable elements (TEs) are usually highly methylated
along their entire length, but DNA methylation within genes is
distributed away from the 5′ and 3′ ends. It has been proposed that DNA
methylation is primarily used to inhibit the transcription and movement
of TEs to protect genomes against uncontrolled insertions. Some studies
have reported that transposons or retrotransposons are activated in
response to environmental stress and may generate true genetic or
epigenetic changes for adapting to stresses [[36]8]. In rice,
epigenetic regulation has been reported to affect the transposition
activity of Tos17 and modulates the activity of the neighboring genes
[[37]9]. Highly methylated transposons may affect transcription of
neighboring genes. For example, the patterns of DNA methylation in long
terminal repeat (LTR) transposable elements differed between rice
leaves and roots and affected the transcription of the neighboring
genes [[38]10].
Recently, three DNA treatment techniques have been typically employed
for DNA methylation detection: (i) methylation sensitive amplified
polymorphism (MSAP), (ii) methylated DNA enrichment (MeDIP), (iii)
bisulfite-sequencing. With the advent of new technologies, the
elucidation of the genome-wide methylation profiles is possible. In
this regard, the immunoprecipitation of methylated DNA by monoclonal
antibodies specific to 5-methylcytidine (5mC) combined with
high-throughput sequencing (MeDIP-seq) has been used as a valuable tool
to map methylated DNA at the genomic scale [[39]11,[40]12]. The
MeDIP-seq approach employs an antibody against 5-methylcytosine or
methyl-binding domain proteins to capture methylated DNA, which are
subsequently subjected to next generation sequencing [[41]11]. After
methylated DNA enrichment, the unmethylated DNA fragments are removed
from the genome samples for reducing the redundancy. Thus, a relatively
lower sequencing throughput is required, compared with
bisulfite-sequencing. In addition, MeDIP is more unbiased and less
limited than MSAP. Due to these advantages, MeDIP-seq is considered to
have great potential for the development of cost-effective and unbiased
strategies for whole genome DNA methylation profiling.
Rice (Oryza sativa L.) is one of the most important grain crops
worldwide and provides a staple food for almost half of the world’s
population [[42]13]. Recently, hybrid rice has made a tremendous
contribution to food security both in China and many other countries
which are losing arable land. Hybrid rice technologies are mainly based
on two well-known male sterility systems, namely CMS (Cytoplasmic Male
Sterile) and EGMS (Environmentally sensitive Genic Male Sterile).
Compared with the three-line system (CMS system), the two-line system,
is based on the discovery and application of EGMS lines, which serve as
both the male sterile lines and maintainer lines under different
environmental conditions [[43]14]. PTGMS (Photoperiod- Thermo-Sensitive
Genic Male Sterile) is the major type of EGMS germplasm resources and
has been widely used for breeding two-line hybrid rice. The PTGMS
system is advantageous for broad restoration ability, easy maintenance,
and multiplication [[44]15-[45]17]. Thus, the utilization of PTGMS
lines is an important approach for the better exploitation of heterosis
in rice.
Peiai 64S (PA64S) is one of the most important indica rice genic male
sterile lines in two-line system. Its paternal line is Peiai64, and its
maternal line is Nongken58, which is the first photoperiod- sensitive
genic male sterile line discovered by Mingsong Shi in Hubei, China in
1973 [[46]18]. Using PA64S as the maternal line, more than 10 two-line
hybrid rice varieties have been developed since 1996 in China. Among
these new varieties, LiangYouPei9 (LYP9) is an elite super hybrid rice
and increases grain yield by approximately 20% per hectare compared
with other hybrid rice lines in the past years [[47]7,[48]19]. The
fertility transition of PA64S was controlled by the photoperiod and
temperature: the sterility trait exhibited at temperatures higher than
23.5°C during the anther development and long-day (LD) (≥14 h)
conditions can suppress the fertility transition. Therefore, under
lower temperatures (~21-23°C) and short-day (SD) (<14 h) conditions
during anther development, PA64S can transform from male sterility to
fertility [[49]20,[50]21].
To date, some studies have reported that a substitution of C-to-G in
the male-fertility-associated noncoding RNA p/tms12-1 in PA64S produced
a mutant small RNA, namely osa-smR5864m. This mutated noncoding small
RNA gene may lead to PTGMS in PA64S [[51]21]. Furthermore, the
transcriptome of the PGMS rice Nongken 58S was significantly suppressed
under LD condition and the circadian rhythm pathway had been shown to
be involved in the male sterility transition [[52]22]. Another
noncoding RNA called LDMAR was also found in Nongken 58S and was
required for the fertility transition. A siRNA (Psi–LDMAR) was detected
to associate with the DNA methylation pattern of the LDMAR promoter in
Nongken 58S. And increasing methylation in the promoter of LDMAR of
Nongken 58S also reduced the expression of LDMAR, leading to male
sterility [[53]23]. These results suggested that DNA methylation or
RNA-dependent DNA methylation (RdDM) might be involved in the
regulation of PTGMS.
Brassinosteroids (BRs), a class of steroid hormones, have been reported
to control male fertility by regulating the expression of several key
genes involved in anther and pollen development, such as SPL/NZZ, TDF1,
AMS, AtMYB103 and MS1 [[54]24,[55]25]. The basic helix-loop-helix
(bHLH) protein BIM1 (BES1-interacting Myc-like1) is a brassinosteroid
signaling component involved in regulating BR-induced genes in
Arabidopsis [[56]26]. BRI1- EMSSUPRESSOR1 (BES1) is shown to interact
with BIM1 and together bind to the E-box of a BR-induced gene [[57]27].
In this study, we investigated the global DNA methylation alterations
in young panicles of the PTGMS line PA64S under two different
environmental conditions using the MedIP-seq method. The aim of this
work was to explore the DNA methylation patterns and their influences
on gene expression during the transition from sterile to fertile in
PA64S. The results will provide some information for a better
understanding of the role of DNA methylation in PTGMS rice.
Results
Cytological observation of pollen morphology
Under higher temperature and LD, anthers from PA64S were light yellow,
and its pollen wizened as well as unstained by I[2]-KI (Figure [58]1a,
b), displaying sterility (S). However, under lower temperatures
(~21-23°C) and SD condition, PA64S showed male fertility. The anthers
of PA64S (fertility, F) were dark yellow, and its spherical pollen
could be darkly stained by I[2]-KI, indicating much starch accumulation
(Figure [59]1c, d).
Figure 1.
Figure 1
[60]Open in a new tab
Cytological observation of pollen morphology. (a) and (c), the anther
phenotypes of PA64S under two environmental conditions. (b) and (d) 1%
I[2]-KI staining of the pollen grains of PA64S under two environmental
conditions.
Genome-wide DNA methylation profiles of PA64S (S) and PA64S (F)
To decipher the genome-wide DNA methylome of PA64S under two different
environmental conditions, we immunoprecipitated sheared genomic DNA
with an antibody which specifically recognizes 5-methylcytosine and
sequenced the enriched methylated DNA. Only the uniquely mapped reads
were used in the scanning of the methylation peak. A total of
24,489,796 raw reads were generated for the two samples, and more than
92% of the reads were mapped, and about 58% of the reads in each sample
were uniquely mapped to the rice genome in each sample (Table [61]1 and
Additional file [62]1).
Table 1.
Data generated by MeDIP-seq in PA64S (S) and PA64S (F)
Total MeDIP-Seq data Percentage of mapped reads in total reads Unique
mapped reads Unique mapped bases Percentage of unique mapped reads
PA64S(S) 24,489,796 92.05% 14,676,561 719,151,489 59.93%
PA64S(F) 24,489,796 92.20% 14,410,405 706,109,845 58.84%
[63]Open in a new tab
According the principle of MedIP method, analysis of the methylation
enrichment in different components of the genome showed that 2 kb
region upstream and 2 kb region downstream of the CpG Island and
intragenic had enriched more methylation reads, especially, when
departing from the CpG Island or intragenic regions (Figure [64]2).
Furthermore, the upstream 2k and 5’UTR or downstream 2k and 3’UTR were
enriched more reads (Figure [65]3a).
Figure 2.
Figure 2
[66]Open in a new tab
Distribution of reads around CpG Islands and gene body. (a)
Distribution of reads around CpG Islands; (b) Distribution of reads
around gene body. The upstream and downstream 2 kb regions were split
into 20 equal regions. In the gene body, each gene was split into 40
equal regions. For each region, the normalized number of reads was
calculated. The “Y” axis is the average of the normalized depth for
each region.
Figure 3.
Figure 3
[67]Open in a new tab
Different distributions of PA64S(S) and PA64S (F). (a) Peak
distribution in different components of the PA64S(S) and PA64S (F)
genome. (b) Distribution of methylated peaks of different repeat types
in PA64S (S) and PA64S (F). (c) Distribution of CpG (O/E) in PA64S(S)
and PA64S (F).
The CpG observed/expected (O/E) ratio of PA64S (S) and PA64S (F) was
mainly in a range of 0.4 to 1.2. There was no significant difference in
the CpG (O/E) ratio between PA64S(S) and PA64S (F). However, PA64S(S)
had more peaks distributed at CpG (O/E) ratios of 0 ~ 0.8, while the
peaks of PA64S (F) were inclined to CpG (O/E) ratios of 1.0 ~ 1.8
(Figure [68]3b). The CpG O/E ratio is a widely used parameter to
predict the DNA methylation level based on C → T transition mechanisms
resulting from deamination of mCs over the course of evolution
[[69]2,[70]28]. In our study, more peaks of PA64S (S) were distributed
in 0 ~ 0.8 (Figure [71]3c), and the methylation level of PA64S (S) was
significant higher than that of PA64S (F) (Figure [72]4). These results
were consistent with previous studies, showing that genes with higher
methylation levels tend to have lower CpG O/E ratios [[73]28,[74]29].
Figure 4.
Figure 4
[75]Open in a new tab
Boxplots showing the 5mC content (read count, y axis) in 1 kb tiled
windows for PA64S(S) and PA64S (F). (a) Total reads of PA64S(S) and
PA64S (F); (b) Reads for each elements of PA64S(S) and PA64S (F). The
asterisks indicate significant differences between PA64S(S) and PA64S
(F), as determined by Student’s t test (**P < 0.01, ***P < 0.001).
To get insight into the DNA methylation status of PA64S (S) and PA64S
(F), RPKMs (Reads per 1 kb window per million mapped reads) were
calculated for each sample and different elements (Figure [76]4). The
results showed that PA64S(S) had higher methylation level than that of
PA64S (F). And the DNA methylation levels between PA64S(S) and PA64S
(F) were observed to be significantly different on the upstream 2k,
5’UTR, intron, 3’UTR, and downstream 2k regions of the gene body
(p < 0.01, p < 0.001). However, there were not significant differences
on CDS (coding DNA sequence) of the differential methylated regions
(DMRs).
DNA methylation patterns in genes and TE regions
Comparison of the gene methylation status showed that there were 725
genes hypermethylated and 533 genes hypomethylated in PA64S (S). Most
of these DMRs were observed on the upstream 2k and downstream 2k of
genes and many of them were hypermethylated in PA64S (S) (Table [77]2
and Figure [78]4). The number of hypermethylated genes showed
significant difference between PA64S (S) and PA64S (F) in the
upstream2k, 5’UTR, 3’UTR, and downstream 2k regions comparing with the
hypomethylated genes (Table [79]2 and Figure [80]4). Interestingly, the
differentially methylated genes of the CDS regions between PA64S (S)
and PA64S (F) exhibited almost no differences (Table [81]2 and
Figure [82]4).
Table 2.
Numbers of differentially methylated genes between PA64S(S) and PA64S
(F) in different gene regions
Contrast PA64S(S) VS. PA64S(F)- Hyper PA64S(S) VS. PA64S(F)- Hypo
Upstream 2k 245 122
5’UTR 25 10
CDS 86 84
Intron 93 69
3’UTR 23 10
Downstream 2k 253 106
Total 725 533
[83]Open in a new tab
The significant enrichment of down-regulated genes in PA64S (S) was
annotated (Table [84]3 and Additional file [85]2). These differential
genes included: transposon, retrotransposon, phosphate carrier protein,
NB-ARC domain containing protein, ent-kaurene synthase, PsbP and so on.
They were involved in disease resistance, signal transduction,
transport, chloroplast development and photosynthesis, energy
metabolism, and pollen development. Further research on these genes may
elucidate the functions of hypermethylation in the fertility transition
of PA64S.
Table 3.
Significant GO enrichment of hypermethylated genes in PA64S (S)
Locus ID DMR location Annotation p-value
LOC_Os01g70080 Upstream 2k/ Intron NB-ARC domain containing protein
1.43E-42
LOC_Os07g04950 Upstream 2k Early nodulin 75-like protein 2.33E-29
LOC_Os01g70100 Upstream 2k/5'-UTR zinc finger DHHC domain-containing
protein 1.21E-22
LOC_Os09g38100 Upstream 2k phosphate carrier protein, mitochondrial
precursor 1.84E-15
LOC_Os09g39570 Upstream 2k beta-amylase 2.29E-15
LOC_Os03g25840 Upstream 2k Amino acid permease family protein 2.8E-14
LOC_Os12g30520 Upstream 2k pumilio-family RNA binding repeat containing
protein 1.06E-13
LOC_Os12g35730 Upstream 2k transposon protein, putative, CACTA, En/Spm
sub-class 1.36E-12
LOC_Os10g31780 Upstream 2k oxidoreductase 2.25E-12
LOC_Os01g04330 Upstream 2k Probable calcium-binding protein CML16
2.77E-11
LOC_Os06g42420 Upstream 2k transposon protein, putative, unclassified
1.08E-10
LOC_Os01g70360 Upstream 2k Long chain base biosynthesis protein 2b
1.21E-10
LOC_Os12g40130 Upstream 2k phosphoribosyl transferase 6.86E-10
LOC_Os01g18620 Upstream 2k transferase family protein 7.06E-10
LOC_Os12g30520 5'-UTR pumilio-family RNA binding repeat containing
protein 1.06E-13
LOC_Os04g30490 5'-UTR/ CDS/ Intron MATE efflux family protein 3.05E-11
LOC_Os01g62890 5'-UTR prenylated rab acceptor 3.75E-11
LOC_Os01g70540 CDS/ Intron retrotransposon protein 3.94E-31
LOC_Os01g20880 CDS OsWAK3 1.06E-15
LOC_Os01g70040 CDS/Intron transposon protein, CACTA, En/Spm sub-class
1.08E-14
LOC_Os01g70030 CDS/Intron transposon protein, CACTA, En/Spm sub-class
1.5E-14
LOC_Os04g02920 CDS/3'-UTR leucine-rich repeat family protein 2.2E-14
LOC_Os04g35210 CDS Leucine Rich Repeat family protein 9.11E-12
LOC_Os06g03850 CDS impaired sucrose induction 1 8.63E-11
LOC_Os07g41610 CDS/Intron retrotransposon protein, Ty1-copia subclass
1.86E-10
LOC_Os06g36540 CDS retrotransposon protein, Ty3-gypsy subclass 9.81E-10
LOC_Os07g04940 Intron PE-PGRS family protein precursor 2.33E-29
LOC_Os08g21920 Intron retrotransposon protein 5.37E-23
LOC_Os06g21950 Intron inorganic phosphate transporter 7.93E-15
LOC_Os12g19549 Intron telomerase reverse transcriptase 1.36E-12
LOC_Os09g23650 Intron FAM10 family protein 2.88E-11
LOC_Os06g03850 Intron impaired sucrose induction 1 8.63E-11
LOC_Os10g30090 Intron amino acid permease 4.79E-10
LOC_Os09g21770 3'-UTR ES43 protein 1.03E-12
LOC_Os04g47410 3'-UTR DHHC zinc finger domain containing protein
5.04E-12
LOC_Os01g70190 Downstream 2k exostosin family domain containing protein
9.18E-65
LOC_Os01g70550 Downstream 2k heparan-alpha-glucosaminide
N-acetyltransferase 3.94E-31
LOC_Os12g36410 Downstream 2k transposon protein, Pong sub-class
9.30E-24
LOC_Os01g70810 Downstream 2k homeobox domain containing protein
2.97E-23
LOC_Os01g70270 Downstream 2k auxin response factor 2.93E-21
LOC_Os04g12060 Downstream 2k retrotransposon protein 2.75E-19
LOC_Os01g67980 Downstream 2k cysteine proteinase EP-B 1 precursor
1.81E-16
LOC_Os01g71070 Downstream 2k xylanase inhibitor 1.73E-15
LOC_Os09g39560 Downstream 2k genetic modifier 2.29E-15
LOC_Os01g70020 Downstream 2k DEK C terminal domain containing protein
3.62E-15
LOC_Os05g23610 Downstream 2k protein phosphatase inhibitor 2 containing
protein 2.20E-14
LOC_Os01g70690 Downstream 2k Rapid ALkalinization Factor RALF family
protein 7.42E-14
LOC_Os02g17780 Downstream 2k ent-kaurene synthase, chloroplast
precursor 8.65E-14
LOC_Os05g12210 Downstream 2k chalcone synthase 1.15E-12
LOC_Os04g47420 Downstream 2k transmembrane amino acid transporter
protein 5.04E-12
LOC_Os07g17390 Downstream 2k PsbP 2.81E-11
LOC_Os10g34910 Downstream 2k secretory protein 3.37E-11
LOC_Os02g15690 Downstream 2k polygalacturonase 5.54E-11
LOC_Os01g01840 Downstream 2k helix-loop-helix DNA-binding domain
containing protein 6.27E-11
LOC_Os06g40200 Downstream 2k calcium-binding mitochondrial carrier
1.04E-10
LOC_Os01g70370 Downstream 2k serine palmitoyltransferase 2 1.21E-10
LOC_Os01g45750 Downstream 2k bile acid sodium symporter family protein
7.05E-10
LOC_Os01g18630 Downstream 2k aspartic proteinase oryzasin-1 precursor
7.06E-10
LOC_Os06g36550 Downstream 2k retrotransposon protein, Ty3-gypsy
subclass 9.81E-10
LOC_Os01g70900 Downstream 2k retrotransposon protein 9.84E-10
[86]Open in a new tab
Transposable elements (TEs) are frequently methylated in plant genomes.
The analysis of 26 repeat types revealed that LTR /Gypsy had the
highest 5mC content in PA64S and was increased in PA64S (F)
(Figure [87]3c). Further analysis of four types of LTR and LINE, SINE
showed that all of the LTR retroelements were highly methylated in
PA64S (F) (Table [88]4). In contrast, the two types of retroelements
(LINE, SINE) had higher methylation in PA64S (S) than in PA64S (F)
(Table [89]4).
Table 4.
Transposable element methylation of PA64S(S) and PA64S (F) based on
10 kb tiled genome windows
Type PA64S(S) PA64S(F)
reads % reads %
LTR/Copia 521,175 5.62 533,088 5.79
DNA/En-Spm 658,162 7.1 676,617 7.35
LTR/Gypsy 4,045,323 43.63 4,293,192 46.61
LTR 55,692 0.6 55,798 0.61
LINE/L1 174,518 1.88 171,516 1.86
SINE 167,785 1.81 154,945 1.68
[90]Open in a new tab
Validation of MedIP data by bisulfite sequencing
In this study, the DNA methylation patterns of 20 DMR-associated genes
were selected to carry out bisulfite sequencing for the validation of
MedIP data, including LOC_Os10g10560 (invertase/pectin methylesterase
inhibitor, IPMI), LOC_Os12g08770 (photosystem I reaction center subunit
N), LOC_Os03g51200 (H2A), LOC_Os08g38210 (transcription factor BIM2)
and so on (Figures [91]5 and [92]6). We found that many of the
DMR-associated genes showed hyper-methylation in PA64S (S), and the
results were almost in accordance with the MeDIP-seq data
(Figures [93]5 and [94]6, and Additional file [95]3). For example,
LOC_Os10g10560 (IPMI), LOC_Os03g01820 (vacuolar ATP synthase subunit
E), LOC_Os06g23980 (MADS-AGL16) and LOC_Os04g37619 exhibited a higher
methylation level in PA64S(S) compared with PA64S (F) (Additional file
[96]3). However, the DMR-associated genes LOC_Os03g51200 (H2A) and
LOC_Os03g27770 (heme oxygenase 2) were observed to be hypomethylation
in PA64S(S) (Figure [97]6a, b).
Figure 5.
Figure 5
[98]Open in a new tab
DNA methylation patterns of four selected DMR-associated genes
validated by bisulfite sequencing. (a) LOC_Os10g10560 (invertase/pectin
methylesterase inhibitor, IPMI); (b) LOC_Os05g12210 (chalcone
synthase); (c) LOC_Os02g17780(ent-kaurene synthase, chloroplast
precursor); and (d) LOC_Os12g08770 (photosystem I reaction center
subunit N). The red, green and blue columns in the histograms refer to
the collective methylation levels (in percentage) of CG, CHG, and CHH,
respectively.
Figure 6.
Figure 6
[99]Open in a new tab
DNA methylation patterns and gene expression of the four selected
DMR-associated genes. (a) Bisulfite sequencing and qRT-PCR analysis of
LOC_Os03g27770 (heme oxygenase 2); (b) Bisulfite sequencing and qRT-PCR
analysis of LOC_Os03g51200 (H2A); (c) Bisulfite sequencing and qRT-PCR
analysis of LOC_Os06g16270 (heat shock factor binding protein 2); and
(d) Bisulfite sequencing and qRT-PCR analysis of LOC_Os08g38210
(transcription factor BIM2). The red, green and blue columns in the
histograms refer to the collective methylation levels (in percentage)
of CG, CHG, and CHH, respectively.
Association of DNA methylation and gene expression
Analysis of the relationship between DNA methylation and gene
expression revealed that hypermethylation in some regions were related
to the expression of neighboring genes. DNA methylation level of
PA64S(S) on the upstream 2k, 5’UTR, intron, 3’UTR, and downstream 2k
regions were hypermethylated compared with PA64S (F) (p < 0.001)
(Figure [100]4b and Table [101]2). Correspondingly, many DMR-associated
genes may be down-regulated in PA64S (S).
The DNA methylation patterns and gene expression of four DMR-associated
genes, namely LOC_Os03g51200 (H2A), LOC_Os08g38210 (transcription
factor BIM2), LOC_Os06g16270 (heat shock factor binding protein 2) and
LOC_Os03g27770 (heme oxygenase 2), were further examined in this study
(Figure [102]6). Three of these had lower DNA methylation levels, but
higher expression in PA64S (S). Furthermore, The DMR-associated gene
LOC_Os08g38210 showed a higher DNA methylation level and a lower
expression in PA64S (S) (Figure [103]6d).
Expression of the genes involved in brassinosteroid signaling pathway
One of the hypermethylated gene OsBIM2 (LOC_Os08g38210), which is
homolog with BIM gene family in Arabidopsis, was thought to be involved
in BR signaling pathway. Comparison with the structure of BIM1
suggested that OsBIM2 had higher similarity to BIM1 of Arabidopsis,
with an extended N-terminal domain that was lacking in BIM2 and BIM3
[[104]30]. Comparing the respective genomic loci also supported the
idea that OsBIM2 and Arabidopsis BIM1 probably share a more ancestral
intron-exon structure (Figure [105]7). And BIM1 was reported to be
cooperated with BES1 to regulate BR-induced genes in Arabidopsis
[[106]26]. Therefore, we postulated that BR may participate in
regulating male fertility in PA64S. To test this hypothesis, we
examined the expression of some genes related to BR signaling pathway
and the genes that regulated by BR. Expression analysis of these genes
showed that some of them had lower expression in PA64S (S) (P < 0.01),
including OsBAK1, OsBZR1, OsSPL8 (Figure [107]8). The expression of
OsMS2, one of the regulating genes for BR signaling in the tapetum and
microspore development, was also dramatically reduced in PA64S (S) (P <
0.01) (Figure [108]8).
Figure 7.
Figure 7
[109]Open in a new tab
The genomic organization of the rice OsBIM2 (LOC_Os08g38210) loci is
compared to those of the BIM1 , BIM2 and BIM3 genes in Arabidopsis.
Note: This genomic organization is modified on the basis of Xing et
al., 2013 [[110]30].
Figure 8.
Figure 8
[111]Open in a new tab
Expression analysis of genes related to BR signaling pathway.
Gene ontology (GO) and KEGG pathway enrichment analysis
In this study, genes that overlapped with the methylation peaks in the
upstream 2k, downstream 2k or gene body regions were termed as
methylated genes. A total of 36,321 methylated genes were found in the
two samples: 32,678 in PA64S (S) and 32,990 in PA64S (F). Of them, 1258
methylated genes were identified to be differentially expressed genes
between PA64S (S) and PA64S (F), and these included 725 genes that were
hypermethylated and 533 genes that were hypomethylated in PA64S (S)
(Table [112]2).
We further used WEGO ([113]http://wego.genomics.org.cn) to functionally
categorize the methylated genes and observed significant differences
(Figure [114]9 and Additional file [115]4). Many genes enriched
hypermethylation in PA64S(S), such as translation regulator,
development process, multi-organism process and so on (Figure [116]9).
The results showed that some genes were significantly down-regulated
genes in PA64S (S), including 12 genes for “cellular component” (9), 16
genes for “molecular function” (21) and 64 genes for “biological
process” (50) (Table [117]5 and Additional file [118]5). Methylated
genes tend to be enriched in transferase activity, including
transcription regulator as well as electron transporter. With respect
to biological processes, the genes are enriched in functions associated
with photorespiration, cellular metabolism, oxidative phosphorylation
and mitochondrial ATP synthesis. Interestingly, the methylated genes of
“response to stimulus” were increased in PA64S (S) (Figure [119]9). In
contrast, for cellular component, the significantly differentially
down-regulated genes in PA64S (S) were enriched in mitochondrial
respiratory chain and mitochondrial inner membrane (Table [120]5).
Figure 9.
Figure 9
[121]Open in a new tab
Gene categories and distribution of differential methylated genes
between PA64S(S) and PA64S (F).
Table 5.
Overrepresented Gene Ontology (GO) categories in the hypermethylated
genes in PA64S (S)
GO term Ontology ^a Element Ontology description Number in input list
Number in BG/Ref p-value
GO:0005750 C Down2k_ Down mitochondrial respiratory chain complex III 4
14 2.74e-05
GO:0005746 C Down2k_ Down mitochondrial respiratory chain 3 22 0.00104
GO:0005743 C Down2k_ Down mitochondrial inner membrane 5 78 0.00223
GO:0016742 F 3’-UTR_ Down hydroxymethyl-, formyl- and related
transferase activity 4 4 1.61e-14
GO:0016741 F 3’-UTR_ Down transferase activity, transferring one-carbon
groups 4 84 3.12e-08
GO:0030528 F Intron_ Down transcription regulator activity 4 178
1.09e-05
GO:0045153 F Down2k_ Down electron transporter, transferring electrons
within CoQH2-cytochrome c reductase complex activity 4 13 4.58e-05
GO:0009853 P 3’-UTR_ Down photorespiration 4 4 3.65e-13
GO:0046653 P 3’-UTR_ Down tetrahydrofolate metabolic process 4 4
3.65e-13
GO:0006760 P 3’-UTR_ Down folic acid and derivative metabolic process 4
8 2.56e-11
GO:0042558 P 3’-UTR_ Down pteridine and derivative metabolic process 4
11 1.20e-10
GO:0043094 P 3’-UTR_ Down cellular metabolic compound salvage 4 19
1.41e-09
GO:0006732 P 3’-UTR_ Down coenzyme metabolic process 4 63 2.17e-07
GO:0009259 P 3’-UTR_ Down ribonucleotide metabolic process 4 81
6.08e-07
GO:0051186 P 3’-UTR_ Down cofactor metabolic process 4 102 1.55e-06
GO:0009117 P 3’-UTR_ Down nucleotide metabolic process 4 111 2.19e-06
GO:0006725 P 3’-UTR_ Down cellular aromatic compound metabolic process
4 129 4.02e-06
GO:0032787 P 3’-UTR_ Down monocarboxylic acid metabolic process 4 142
5.94e-06
GO:0046483 P 3’-UTR_ Down heterocycle metabolic process 4 189 1.88e-05
GO:0006122 P Down2k_ Down mitochondrial electron transport, ubiquinol
to cytochrome c 4 6 2.35e-06
GO:0042775 P Down2k_ Down mitochondrial ATP synthesis coupled electron
transport 4 18 0.00045
GO:0006119 P Down2k_ Down oxidative phosphorylation 4 28 0.00289
GO:0042773 P Down2k_ Down ATP synthesis coupled electron transport 4 28
0.00289
[122]Open in a new tab
^a, C: cellular component; F, molecular function; P, biological
process.
To generate further insight view of pathway perturbation in PA64S (S)
under long day length and high temperature conditions, we performed
KEGG pathway analysis of the differentially methylated genes. As a
result, 58 predicted pathways were enriched, including metabolism,
circadian rhythm, photosynthesis, porphyrin and chlorophyll metabolism,
peroxisome and oxidative phosphorylation (Additional files [123]6).
Through the GO term enrichment and KEGG pathway analysis of the
differentially regulated genes, we observed that certain genes were
hypermethylated in PA64S (S) and that they were involved in
photorespiration and mitochondrial ATP synthesis pathways. These
pathways were closely connected to photosynthesis and energy synthesis.
For instance, cytochrome c1-1 (LOC_Os05g23620) was involved in
mitochondrial electron transport.
Discussion
Although global DNA methylation surveys have been performed in some
plants [[124]3,[125]31,[126]32], little studies have been reported
about the methylation profiles of PTGMS rice (for example PA64S) during
the fertility transition. In this study, the objective was to perform a
genome-wide identification of the methylated genes that affect the
sterility-fertility transition in PTGMS rice PA64S. To confirm the
results from MeDIP-seq, the methylation patterns of twenty regions in
each sample were analyzed by bisulfite sequencing. And the methylation
levels obtained by the two methods were generally in accord with each
other. The read distribution analysis performed in this study found
that uniquely mapped reads were enriched in the upstream 2k, 5’UTR,
3’UTR and downstream 2k regions. In addition, the rice sterile line
PA64S(S) had higher methylation level than that of PA64S (F), and many
DMR-associated genes were down-regulated in PA64S (S).
Transposable elements (TEs) are abundant features in plant genomes that
are frequently marked by cytosine methylation. In this study, we
examined many categories of transposable elements from the RepBase
transposable element annotation database. LTR retroelements are the
predominant type of transposable element sequences in plant genomes
[[127]33]. Our results showed that the LTR categories had higher 5mC
content in PA64S (F) than in PA64S (S). And the LINE and SINE
retroelement categories had an increased in 5mC content in PA64S (S)
compared with PAS64 (F), although this difference was not significant.
Previous studies using methylation-deficient mutants have demonstrated
that hypomethylation can result in transposable element mobilization
[[128]34]. The increased 5mC content of TEs found in the present study
suggested that the operation of adaptive genomic mechanisms might help
protect genomes against transposable element spread. It would be of
interest to measure the relationship of DNA methylation and TE elements
during the fertility transition of PA64S.
In this study, the DNA methylation patterns of twenty DMR-associated
genes were validated through bisulfite sequencing and we found that
hypermethylation in some regions affected the expression of neighboring
genes (Figures [129]5 and [130]6, Additional file [131]3). Although
hypermethylated DMR-associated genes did not show significant
enrichment of down-regulated genes, there were still some
down-regulated genes with higher methylation levels (Figures [132]5 and
[133]6). Some flower/anther development related genes were
hypermethylated in PA64S (S), for example, LOC_Os06g23980 (putative
MADS-box protein AGL16) and LOC_Os10g10560 (IPMI) (Additional file
[134]3). Previous studies of a male sterile mutant (Ms-cd1) in Brassica
oleracea have demonstrated that IPMI was specifically expressed in
anther and filament and was not expressed in CMS NILs [[135]35]. The
increased 5mC content in IPMI of PA64S(S) found in the present study
might lead to a lower expression of IPMT and most likely disturbed the
anther cell wall metabolism. Furthermore, photosystem and mitochondrial
ATP synthesis associated genes, including LOC_Os12g08770 (photosystem I
reaction center subunit N), LOC_Os07g17390 (PsbP) and LOC_Os03g01820
(vacuolar ATP synthase subunit E) were also hypermethyalted in PA64S
(S) (Additional file [136]2). High temperature (>23.5°C) and LD (14 h)
conditions increased the methylation levels of these genes and
down-regulated their expression during anther development in PA64S,
which might in turn lead to PA64S sterility.
In contrast, some heat shock proteins, such as LOC_Os08g43490 (heat
shock protein DnaJ), LOC_Os06g16270 (heat shock factor binding protein
2) were found to be hypomethylated in PA64S (S) in our study, and the
expression level of LOC_Os06g16270 was also higher in PA64S (S)
(Figure [137]6c). Because heat shock proteins are always expressed in
response to stress, their higher expression in PA64S (S) suggested that
the adverse environment might affect anther development. These findings
suggested that different environments might regulate the gene
expression through DNA methylation in the gene body or flanking
sequence of fertility associated genes.
In Arabidopsis and rice, the BR signaling pathway is similar.
Generally, BRs bind the extracellular domain of BRI1, and then BRI1
heteromerizes with the BRI1-associated kinase 1 (BAK1). The receptor
complex then activates the BR-signaling kinases (BSKs) and inactivates
the glycogen synthase kinase-3-like kinase BIN2, resulting in nuclear
accumulation of unphosphorylated forms of two transcription factors
BZR1 and BES1, which activate the expression of BR-responsive genes
[[138]36-[139]38]. In this study, some of the components in BR
signaling pathway, such as OsBRI, OsBAK1, OsBZR and OsBIM2 were
dramatically reduced in male sterile line PA64S (S) (Figure [140]8).
Our results indicated that BR signaling involved in regulating the
sterility-fertility transition of PA64S under different environments.
On the other hand, OsBIM2 (LOC_Os08g38210) showed a higher DNA
methylation level in PA64S(S) and its expression dramatically reduced
(Figure [141]6d). In Arabidopsis, BIM1 and SPL8 act cooperatively in a
common complex or via synergistic pathways to promote male fertility
[[142]30]. In our study, both OsBIM2 and OsSPL8 were downregulated in
PA64S(S) (Figure [143]8). Since OsBIM2 is homolog with BIM1, the lower
expression level will affect BR signaling pathway and then regulated
the downstream genes which involved in anther and pollen development.
Therefore, DNA methylation of BR signaling related genes may affect the
expression of key genes (for example, OsMS2) for anther development,
leading to male sterility under long-day and high temperature
conditions. However, whether the hypermethylated OsBIM2 directly
suppress key genes need to be further investigated.
The GO term significant enrichment analysis showed that many genes,
including the transposons and retrotransposons, were hypermethylated in
PA64S (S) (Tables [144]3 and [145]5 and Additional files [146]3 and
[147]4). For example, LOC_Os09g38100 (phosphate carrier protein
mitochondrial precursor) and LOC_Os06g40200 (calcium-binding
mitochondrial carrier) were observed higher methylation levels in PA64S
(S) than PA64S (F) (Table [148]3 and Additional file [149]3). These two
genes participate in mitochondrial respiratory chain and facilitate
mitochondrial ATP synthesis [[150]39]. Thus, their hypermethylated
would affect the energy synthesis during pollen development. Another
hypermethylated gene LOC_Os02g17780 (ent-kaurene synthase, chloroplast
precursor) was an early intermediate in GA biosynthesis [[151]40].
Other significantly differentially expressed DMR-associated genes
included LOC_Os12g40130 (phosphoribosyl transferase), LOC_Os04g30490
(MATE efflux family protein) and LOC_Os05g12210 (chalcone synthase).
Some studies have reported that the adenine phosphoribosyl transferase
activity deficient mutants of Arabidopsis thaliana were male sterile
due to abortion of pollen development after the meiotic divisions of
pollen mother cells [[152]41]. The MATE (Multidrug And Toxic Compound
Extrusion) efflux protein family has also been reported to be
differentially expressed in sterile and fertile NILs of Brassica
oleracea [[153]35]. Chalcone synthase catalyzes the initial step of the
branch of phenylpropanoid pathway which leads to flavonoids. The
absence of chalcone synthase activity was found to have a pleiotropic
effect in maize and petunia mutants, and pollen fertility/pollen
germination and flavonoid synthesis was disrupted [[154]42]. Therefore,
these differentially expressed genes might be involved in fertility
transition.
To uncover the regulatory mechanisms of the fertility transition, the
regulatory network of differentially methylated genes was examined in
this study. Several important pathways were found, including flavone
and flavonol biosynthesis, circadian rhythm, photosynthesis, porphyrin
and chlorophyll metabolism, and oxidative phosphorylation (Additional
files [155]6). Previous studies have reported that flavone and flavonol
biosyntheses are associated with pollen development. Flavonoids are
required for pollen viability in many species [[156]42-[157]44]. In
addition, flavonols can promote pollen germination frequency and pollen
tube growth in tobacco [[158]45]. Many higher plants use the day length
as an environmental cue to switch from vegetative to reproductive
growth and plants with defects in their circadian regulation cannot
properly regulate the timing of the floral transition [[159]46]. In
PGMS rice, the circadian rhythm genes OsPRR37 in inflorescences have
been reported to directly affect the male sterility transition in
Nongken 58S [[160]22]. In this study, some DMR-associated genes, such
as LOC_Os05g12210 (chalcone synthase) and LOC_Os08g38210 (transcription
factor BIM2), were found to be involved in the circadian rhythm. The
other three pathways, namely photosynthesis, porphyrin and chlorophyll
metabolism, and oxidative phosphorylation were related to carbon or
energy metabolism [[161]47,[162]48]. The down-regulated genes of these
pathways in PA64S (S) may lead to insufficient energy and materials for
pollen development. Therefore, these five pathways were regarded as
potentially related to the fertility transition of PA64S in the present
study.
Conclusions
We systematically investigated the complete DNA methylome of
photoperiod- and thermo-sensitive male sterile rice Peiai 64S under two
different environments (different temperatures and day lengths). From
the whole-genome DNA metylation map, hypermethylation was observed in
PA64S (S) and many pathways such as circadian rhythm, photosynthesis
and oxidative phosphorylation were participated in the fertility
transition of PA64S. Furthermore, the hypermethylated OsBIM2 may be
involved in brassinosteroid signaling pathway and affected the male
fertility in PA64S.
Methods
Plant materials
The rice photo-thermo-sensitive genic male sterile line PA64S, which
was maintained at Wuhan University, was used in this study. The seeds
of PA64S were sown under two different natural ecological conditions:
in Hainan (18°48′ N, 110°02′ E) from December, 2010 to March, 2011
(fertility, F), and in Wuhan (30°30′ N, 114°18′ E) from May 2011 to
August 2011 (sterility, S). Young panicles of PA64S of the two
phenotypes at the meiosis stage were collected, and the pollen
fertility was assayed by staining with potassium iodide (1% I[2]-KI).
DNA preparation and methylated DNA immunoprecipitation
Total genomic DNA was extracted using a Genomic DNA Miniprep Kit
(Axygen) following the manufacturer’s instructions. The DNA quality was
then evaluated by agarose gel electrophoresis and a BioPhotometer Plus
spectrophotometer (Eppendorf, Germany).
Subsequently, the two samples were sonicated to produce DNA fragments
ranging from 100 to 500 bp. After end repairing, phosphorylating and
A-tailing with Paired-End DNA Sample Prep kit (Illumina, USA), the DNA
was ligated to an Illumina sequencing primer adaptor. Then the
fragments were used for methylated DNA immunoprecipitation (MeDIP)
enrichment using a Magnetic Methylated DNA Immunoprecipitation kit
(Diagenod, Belgium) following the manufacturer’s recommendations and
the qualifying DNA was used for PCR amplification. Then bands between
220 and 320 bp were excised from the gel and purified with a QIAquick
Gel Extraction Kit (Qiagen, Germany). The products were quantified with
a Quant-iT™ dsDNA HS Assay Kit (Invitrogen, USA) using an Agilent 2100
Analyzer (Agilent Technologies, USA). Following qPCR qualification, DNA
libraries were sequenced on the Illumina HiSeq 2000 (Illumina, CA, USA)
to generate paired-end 50-bp reads by the Beijing Genomics Institute
(BGI, China).
Bioinformatic analysis
The raw data obtained from Illumina sequencing were first processed to
filter out reads containing adapters, unknown or low quality bases and
were then were mapped to the rice reference genome
([163]ftp://ftp.plantbiology.msu.edu/pub/data/Eukaryotic_Projects/o_sat
iva) by SOAPaligner v 2.21 ([164]http://soap.genomics.org.cn/) with no
more than 2 bp mismatches [[165]49]. The uniquely mapped data were
retained for read distribution analysis, including the distribution in
rice chromosomes and the distribution in different components of the
genome.
Subsequently, the regions of differential methylation (DMRs) with DNA
methylation peaks were employed for differential DNA methylation
analysis. Methylated regions were deemed significant differentially
methylated between PA64S (S) and PA64S (F) with a p-value < 0.05, a
false discovery rate (FDR) < 0.05 and at least a 2.0-fold change in
sequence counts. The base data manipulations and statistical analysis
were performed using the R package ([166]http://www.r-project.org).
Significant differences of DNA methylation level between PA64S (S) and
PA64S (F) were evaluated using Student’s t test.
Gene Ontology and KEGG pathway analysis
Significantly DMRs in the whole gene body (DMGs) common to comparisons
between PA64S (S) and PA64S (F) with opposing directions of methylation
difference (for example, down-methylated or up-methylated in PA64S (S)
relative to PA64S (F)) were screened and annotated with Gene Ontology
(GO), and functionally annotated through KEGG Orthology (KO) enrichment
analysis. The GO enrichment analysis of the functional significance
applies a hypergeometric test to map all of the differentially
expressed genes (DEGs) to terms in the GO database, looking for
significantly enriched GO terms in the differentially expressed genes
compared with the complete genome. The formula used in this analysis is
the following:
[MATH: P=1‐∑i=0
m−1iM
n−iN−MnN
:MATH]
where N is the number of all genes with a GO annotation in the rice
genome, n is the number of those N genes that are differentially
expressed, M is the number of all of the genes that are annotated to
certain GO terms, and m is the number of differentially expressed genes
in M. The hypergeometric test was also performed using a user written
program in R (R Core Team, Vienna, Austria). The KEGG pathway
enrichment analysis identifies significantly enriched metabolic
pathways or signal transduction pathways in differentially expressed
genes comparing with the whole genome background. The calculating
formula and the applied program were the same as that used in the GO
analysis. Then, FDR (Q-value) ≤0.05 was used as the threshold to
determine the most significantly enriched pathways in DEGs.
Bisulfite sequencing
Genomic DNA was modified using an EpiTect Bisulfite kit (Qiagen,
Germany) according to the manufacturer’s instructions. Briefly, 1 μg of
DNA was first treated with bisulfite containing a C-T conversion
reagent and then incubated at 95°C for 5 min, and 64°C for 5 h. The
modified DNA was purified using a MinElute DNA spin columns (Qiagen,
Germany) and stored at -20°C until use. Primers for selected genic
sequences for bisulfite sequencing were designed by the Meth-Primer
Program ([167]http://www.urogene.org/methprimer/) and were given in
Additional file [168]7. For each PCR reaction, 1.0 μl of bisulfite
treated DNA was used in a 50 μl reaction system. These PCR products
were gel-purified using a Gel Extraction Kit (Omega, USA), and then
cloned into the pMD18-T vector (Takara, Dalian, China) and sequenced.
At least fifteen clones were sequenced for each sample. The methylation
levels which are expressed as percentage (%) per site for each of the
three types of cytosines, CG, CHG and CHH, were calculated by dividing
the number of non-converted (methylated) cytosines by the total number
of cytosines within the assay. The analyses of the bisulfite sequencing
results were conducted using the Kismeth website
([169]http://katahdin.mssm.edu/kismeth).
RNA isolation and real time RT-PCR
Total RNA was extracted from young panicles using an RNA isolation kit
(BioDev, Beijing, China) and treated with DNase I (Fermentas, Canada)
to remove the genomic DNA contamination. Approximately 2 μg of RNA from
each sample was used for reverse transcription using the RevertAid
First Strand cDNA Synthesis Kit (Fermentas, Canada). The reaction was
carried out in a total volume of 20 μl, including 2 μg of RNA, 1 μl of
oligo (dN)18 primers (100 μmol/L), 4 μl 5× reaction buffer, 2 μl 10
mmol/L dNTP mix, 1 μl of RiboLock RNase Inhibitor (20 U/μl), and 1 μl
of M-MuLV Reverse Transcriptase (Fermentas, Canada). The procedure was
as follows: the mixture was first denatured at 70°C for 10 min and then
incubated at 37°C for 5 min, followed by 42 °C for 60 min and at last
incubated at 70°C for 15 min.
Quantitative real-time PCR (qPCR) was performed using a StepOne Plus
real-time PCR system (Applied Biosystems) in a total volume of 10 μl
reaction mixture, which included 1 μl cDNA, 0.4 μl each primer
(5 μmol/L each), 0.2 μL ROX Reference Dye (50×) and SYBR premix Ex Taq
(Tli RNaseH plus) (2×) (Takara, Dalian, China). The primers used for
the qPCR experiments are given in Additional file [170]7. The qPCR
reactions were repeated three times for each sample from three
independent experiments. After the amplification steps, the melting
curve was determined for each primer pair to verify that only one
specific product had been amplified. Relative quantification was
performed using the 2^△△Ct method [[171]50]. Rice actin1 was used as a
reference for the mRNA expression analyses.
Availability of supporting data
All MeDIP-seq raw data supporting the results of this article were
deposited in the sequence read archive repository
([172]http://www.ncbi.nlm.nih.gov/sra/) under Sequence Read Archive
(SRA) database (accession # SRP050401).
Acknowledgements