Processes Underlying a Reproductive Barrier in indica-japonica Rice Hybrids Revealed by Transcriptome Analysis1[OPEN] Yanfen Zhu, Yiming Yu, Ke Cheng, Yidan Ouyang, Jia Wang, Liang Gong, Qinghua Zhang, Xianghua Li, Jinghua Xiao, and Qifa Zhang 2 National Key Laboratory of Crop Genetic Improvement and National Center of Plant Gene Research (Wuhan), Huazhong Agricultural University, Wuhan 430070, China ORCID IDs: 0000-0003-1436-810 (X.L.); 0000-0003-2693-2685 (Q.Z.).

In rice (Oryza sativa), hybrids between indica and japonica subspecies are usually highly sterile, which provides a model system for studying postzygotic reproductive isolation. A killer-protector system, S5, composed of three adjacent genes (ORF3, ORF4, and ORF5), regulates female gamete fertility of indica-japonica hybrids. To characterize the processes underlying this system, we performed transcriptomic analyses of pistils from rice variety Balilla (BL), Balilla with transformed ORF5+ (BL5+) producing sterile female gametes, and Balilla with transformed ORF3+ and ORF5+ (BL3+5+) producing fertile gametes. RNA sequencing of tissues collected before (MMC), during (MEI), and after (AME) meiosis of the megaspore mother cell detected 19,269 to 20,928 genes as expressed. Comparison between BL5+ and BL showed that ORF5+ induced differential expression of 8,339, 6,278, and 530 genes at MMC, MEI, and AME, respectively. At MMC, large-scale differential expression of cell wall-modifying genes and biotic and abiotic response genes indicated that cell wall integrity damage induced severe biotic and abiotic stresses. The processes continued to MEI and induced endoplasmic reticulum (ER) stress as indicated by differential expression of ER stress-responsive genes, leading to programmed cell death at MEI and AME, resulting in abortive female gametes. In the BL3 +5+/BL comparison, 3,986, 749, and 370 genes were differentially expressed at MMC, MEI, and AME, respectively. Large numbers of cell wall modification and biotic and abiotic response genes were also induced at MMC but largely suppressed at MEI without inducing ER stress and programed cell death , producing fertile gametes. These results have general implications for the understanding of biological processes underlying reproductive barriers.

Reproductive isolation, a phenomenon widely existing in natural organisms, promotes speciation and maintains the integrity of species over time (Coyne and Orr, 2004). It is divided into two general categories depending on whether it occurs before or after fertilization: prezygotic reproductive isolation and postzygotic reproductive isolation (Seehausen et al., 2014). The former prevents the formation of hybrid zygotes, while the latter results in hybrid incompatibility, including hybrid necrosis, weakness, sterility, and lethality in F1 or later generations (Ouyang and Zhang, 2013). The Dobzhansky-Muller model suggests that hybrid incompatibility is caused by the negative interactions between functionally divergent genes in the parental species (Dobzhansky, 1937; Muller, 1942). 1

This work was supported by grants from the National Key Research and Development Program (2016YFD0100903) and the National Natural Science Foundation (31130032) of China. 2 Address correspondence to [email protected]. Y.Z. and Qifa Zhang designed the research; Y.Z., Y.Y., K.C., Y.O., J.W., L.G., and Qinghua Zhang performed the research; X.L. and J.X. contributed reagents; Y.Z. and Qifa Zhang analyzed data and wrote the article. The author responsible for distribution of materials integral to the findings presented in this article in accordance with the policy described in the Instructions for Authors (www.plantphysiol.org) is: Qifa Zhang ([email protected]). [OPEN] Articles can be viewed without a subscription. www.plantphysiol.org/cgi/doi/10.1104/pp.17.00093

Rice (Oryza sativa), a major food crop and model plant for genomic study of monocotyledon species, also provides a model system for studying reproductive isolation. The Asian cultivated rice is composed of two subspecies, indica and japonica. Their hybrids are usually highly sterile, which is a barrier for utilization of the strong vigor in indica-japonica hybrids. A number of loci responsible for hybrid sterility have been identified, including ones that cause embryo-sac sterility, pollen sterility, and spikelet sterility (Ouyang et al., 2009). So far, six loci responsible for hybrid incompatibility have been cloned. Among them, DPL1/DPL2 (Mizuta et al., 2010), S27/S28 (Yamagata et al., 2010), and Sa (Long et al., 2008) are responsible for pollen fertility, while hsa1 (Kubo et al., 2016), S7 (Yu et al., 2016), and S5 (Chen et al., 2008; Yang et al., 2012) control female fertility; all of them cause spikelet sterility. Three evolutionary genetic models have been proposed to summarize the findings: parallel divergence, sequential divergence, and parallel-sequential divergence (Ouyang and Zhang, 2013). DPL1/DPL2 and S27/S28 in rice, HPA1/HPA2 in Arabidopsis (Arabidopsis thaliana), and JYAlpha in Drosophila melanogaster conform to the parallel divergence model. Sa, hsa1, and many loci in other species can be explained by the sequential divergence model. However, S5 in rice is the only case identified so far to follow the parallel-sequential divergence model. S5 is a major locus controlling indica-japonica hybrid sterility identified in many studies (Ikehashi and Araki,

Plant PhysiologyÒ, July 2017, Vol. 174, pp. 1683–1696, www.plantphysiol.org Ó 2017 American Society of Plant Biologists. All Rights Reserved.

1683

Zhu et al.

1986; Liu et al., 1992, 1997; Wang et al., 1998, 2005; Qiu et al., 2005; Song et al., 2005; Chen et al., 2008) and is also associated with the divergence and domestication of rice (Du et al., 2011; Ouyang et al., 2016). In a typical indicajaponica hybrid, the preferential abortion of female gametes with japonica allele is regulated cooperatively by three adjacent genes, ORF3, ORF4, and ORF5, at the S5 locus (Yang et al., 2012). Typical japonica and indica varieties carry ORF3-ORF4+ORF5- and ORF3+ORF42ORF5+, respectively, where the functional allele of each gene is indicated with “+” while the nonfunctional allele as “2.” ORF5 encodes an aspartic protease with cell wall localization. The indica allele ORF5+ differs from the japonica allele ORF52 by two nucleotides, both of which cause amino acid changes (Chen et al., 2008). The japonica allele ORF4+ encodes an unknown transmembrane protein, while the indica allele ORF42 has an 11-bp deletion resulting in an endoplasmic reticulum (ER) localization. ORF3 is predicted to encode a heat shock protein (HSP) resident in ER functioning to resolve ER stress. The japonica allele ORF32 has a frameshift in the C terminus relative to the indica ORF3+ due to a 13-bp deletion. In this system, extracellular ORF5+ presumably catalyzes its substrate, producing a substance that is sensed by the transmembrane protein ORF4+ that induces ER stress and subsequently leading to premature programmed cell death (PCD). It is speculated that ORF3+, but not ORF32, can prevent/resolve the ER stress caused by ORF5+/ORF4+. Based on the results of genetic and cytological analyses, Yang et al. (2012) proposed a killerprotector model for the roles of these three genes in this system: ORF5+, together with ORF4+, kills the gametes by inducing ER stress that causes premature PCD. ORF3+ protects the gametes from being killed and thus restores the hybrids fertility by preventing/resolving ER stress. In recent years, transcriptomic analysis based on RNAsequencing technology has become a powerful tool to characterize regulatory networks and pathways to enhance understanding of biological processes (Vogel et al., 2005; Digel et al., 2015; Groszmann et al., 2015; Coolen et al., 2016). In this study, we performed a comparative transcriptomic analysis of pistils from Balilla (a typical japonica variety), BalillaORF5+ (Balilla with a transformed ORF5+), and BalillaORF3+ORF5+ (Balilla with transformed ORF3+ and ORF5+) before, during, and after the megaspore meiosis. The comparisons detected thousands of genes that were differentially regulated in these three genotypes at various stages. Functional classification of the differentially regulated genes identified major patterns of differential expression, which suggested the biological processes underlying the killerprotector system of the S5 locus. These results provide insights into the functional roles of ORF3+, ORF4+, and ORF5+ genes as well as the cellular and biochemical nature of their interactions in regulating indica-japonica hybrid sterility. RESULTS Ovule Abortion Process of BL5+

Previous study (Yang et al., 2012) showed that the S5 reproductive barrier is composed of three tightly linked 1684

genes, ORF3, ORF4, and ORF5, in which ORF5+ and ORF4+ together function as the killer of the female gametes and ORF3+ protects the gametes from killing. A typical indica variety has the haplotype of ORF3+, ORF42, and ORF5+, and japonica variety is of ORF32, ORF4+, and ORF52. Balilla (abbreviated BL) is a typical japonica variety. Transgenic plants BalillaORF5+ (abbreviated BL5+) carrying a transformed ORF5+ under its native promoter showed extremely low spikelet fertility (6.50%) compared with wild-type BL (75.52%). In contrast, BalillaORF3+ORF5+ (abbreviated BL3+5+), carrying both transformed ORF5+ and transformed ORF3+ (also under its native promoter), had normal spikelet fertility (77.35%) due to the protective role of ORF3+ (Supplemental Fig. S1, A–C). Paraffin section results showed that the ovule development process in both BL and BL3+5+ was normal. In brief, the archesporial cell enlarged initially and developed into a megaspore mother cell (Fig. 1, A and K). The megaspore mother cell then underwent meiotic division, forming four megaspores in a line (Fig. 1, B, C, L, and M). The three megaspores near the micropyle degenerated, while the megaspore nearest the chalaza enlarged and developed into a functional megaspore (FM; Fig. 1, D and N). The FM underwent three rounds of mitotic division to produce an eight-nucleate embryo-sac (Fig. 1, E and O). In the aborted ovule of BL5+, the megaspore mother cell, dyad, and tetrad could be produced as in BL (Fig. 1, F–H). However, the nucellus cells showed vacuolization and even degeneration at the meiosis stage (Fig. 1, G and H). After meiosis, morphological abnormalities of both nucellus cells and FM were observed (Fig. 1I). The abnormal FM did not enlarge and failed to develop into a mature embryo-sac. The aborted embryo-sac accumulated unknown substances that were deeply stained (Fig. 1J). Thus, the megaspore meiosis process in BL5+ was normal, but the nucellus cells around the megaspores showed serious defects at the megaspore meiosis stage. The nucellus cell degradation and FM abortion caused the sterility of BL5+. ORF5 and ORF4 were neighboring genes overlapping in their promoter region (Yang et al., 2012). These kinds of gene pairs are frequently coexpressed (Wang et al., 2009; Kourmpetli et al., 2013). ORF5 and ORF4 had a similar expression profile and relatively high expression in reproductive organs (Yang et al., 2012). The in situ hybridization showed that ORF5+ had localized expression in developing ovules (Chen et al., 2008). Similarly, ORF4+ also displayed high and specific expression in ovule tissues (Supplemental Fig. S2, G–K). The specific expression of ORF5+ and ORF4+ in ovules was in accordance with the specific abortion of ovules in BL5+. However, the protector ORF3+ showed a wide range of expression in many tissues including ovules, indicating its possible function in many processes (Supplemental Fig. S2, A–F). Based on the ovule development process, we divided the ovule abortion process into three stages (Fig. 1). The first was the megaspore mother cell (MMC) stage when the ovule had megasporocyte and no obvious abnormalities were observed in BL5+. The second was the Plant Physiol. Vol. 174, 2017

Transcriptomic Basis of a Rice Reproductive Barrier

Figure 1. The ovule development process of BL, BL5+, and BL3+5+. A to E, Ovule development of BL (gametes fertile). F to J, Ovule development of BL5+ (gametes sterile). K to O, Ovule development of BL3+5+ (gametes fertile). A, F, K, The MMC stage. B, G, L, The dyad stage. C, H, M, The tetrad stage. D, I, N, The FM stage. E, J, O, The mature embryo sac. Dy, dyad; Te, tetrad; DM, degenerated megaspore; AFM, aborted functional megaspore; ANC, abnormal nucellus cells; ES, embryo sac; AES, aborted embryo sac. Bars = 30 mm. The gray bars at the top indicate the corresponding stage in transcriptome data.

meiosis (MEI) stage when meiosis was in progress and the nucellus cells of BL5+ were vacuolated. The third was after meiosis (AME) when some unknown substances accumulated in embryo-sac (Fig. 1). Processes Underlying Gamete Abortion in BL5+

We collected the pistils of wild-type BL, sterile BL5+, and fertility-restored BL3+5+ at MMC, MEI, and AME, and performed transcriptomic analysis. A data matrix of three stages from three genotypes was generated. A total of 19,269 to 20,928 genes were detected as expressed in these tissues at different stages representing 34.53 to 37.50% of total annotated genes of the Nipponbare reference genome (Supplemental Table S1; Kawahara et al., 2013). A total of 8,339, 6,278, and 530 genes were differentially expressed in BL5+ compared with BL at MMC, MEI, and AME, respectively, including 4,234, 2,430, and 410 genes that were up-regulated at the three stages and 4,105, 3,848, and 120 genes that were down-regulated (Supplemental Table S2; Supplemental Data Set 1). Among the differentially expressed genes (DEGs), 2,067 were up-regulated and 3,325 were down-regulated at both MMC and MEI (Fig. 2A). These 5,392 common DEGs accounted for 64.66 and 85.89% of total DEGs at MMC and MEI, respectively, indicating that MMC and MEI involved substantially overlapping cellular and molecular responses. Plant Physiol. Vol. 174, 2017

An Overview of the Differentially Induced Processes in BL5+ Revealed by the Transcriptomes

We performed MapMan and Gene Ontology (GO) enrichment analyses to classify the DEGs. MapMan is an effective tool to analyze metabolic pathways and other biological processes of the transcriptome systematically (Thimm et al., 2004; Ramsak et al., 2014). To get a global view of the 8,339 DEGs in BL5+ at MMC, we used the “Metabolism overview” and “Cellular response” of MapMan to outline the transcriptional changes of genes with putative functions in metabolism and cellular response (Fig. 2, B and C). In “Metabolism overview,” primary metabolic processes, including cell wall (354 genes), amino acid (226 genes), lipids (316 genes), nucleotide (81 genes), and major carbohydrate (128 genes) metabolism, accounted for 13.25% of total DEGs. The cell wall localization of ORF5+, which was previously hypothesized to initiate the premature PCD (Yang et al., 2012), led us to focus our attention on the cell wall metabolic genes. We found that 187 (128 up-regulated and 59 down-regulated) of the 354 genes were involved in cell wall degradation and modification, which may affect the cell wall structures in one way or another. The secondary metabolism processes, including terpenes (64 genes), flavonoids (105 genes), and phenylpropanoids (120 genes) metabolism, which are important in plant stress response (Dixon and Paiva, 1995; Winkel-Shirley, 2002; Tholl, 1685

Zhu et al.

Figure 2. Analysis of genes differentially expressed in BL5+ compared to BL. A, Numbers of the up- and down-regulated DEGs at MMC and MEI. B, Numbers of DEGs in different MapMan functional categories at MMC, MEI, and AME, respectively. C, The heat map of “Metabolism overview” and “Cellular response” in BL5+ at MMC based on the MapMan results. Each small square represents a gene differentially expressed in BL5+ versus BL. The color of the square represents the level of up-regulation (red) or down-regulation (green) based on the log2 fold change of the DEGs.

2006), constituted 3.47% of all DEGs. In the term of “Cellular response,” 8.65% (721 genes) of DEGs were responsive to biotic and abiotic stress; 1.94% (162 genes) of DEGs were related to redox, a phenomenon usually triggered by stress; and 6.67% (667 genes) of DEGs functioned in developmental processes. GO analysis for up-regulated genes revealed that three GO terms related to cell wall metabolism, seven GO terms involved in stress or stimulus response, and eleven GO terms associated with gene expression, protein synthesis, and translation were significantly enriched (false discovery rate [FDR] , 0.005; Fig. 3). The change in cell wall structure by cell wall modification genes may cause cell wall integrity damage, which 1686

was similar to the situation of pathogen or pest attack. We thus explored the “Biotic stress” in BL5+ at MMC and found that DEGs in various functional categories of biotic stress were enriched, including R genes, pathogen-related proteins, signaling, hormone signaling, proteolysis, ROS regulation, and biotic stress-responsive transcription factors (Supplemental Fig. S3A). Thus, both the MapMan and GO results indicated that a severe stress response including both biotic and abiotic stresses was induced in BL5+ at MMC, which may be caused by the damage in cell wall integrity, even though no abnormalities were observed at this stage (Fig. 1F). At MEI, 163 of the 300 DEGs in the “cell wall” were related to the cell wall degradation and modification Plant Physiol. Vol. 174, 2017

Transcriptomic Basis of a Rice Reproductive Barrier

functional categories of biotic stress were identified (Fig. 2B; Supplemental Fig. S3B), suggesting that the stress response also continued at MEI. By the AME stage, DEG numbers in all the metabolic processes and cellular responses were largely reduced compared to those at MMC and MEI (Fig. 2B). The ovule had aborted with the accumulation of unknown substances. In contrast, however, a GO term “cellulose biosynthetic process” was significantly enriched (Fig. 3), indicating misregulation of cellulose biosynthesis in BL5+. Large-Scale Transcriptional Changes of Cell Wall-Modifying and -Degrading Genes in BL5+ at MMC and MEI

Figure 3. Partial significantly (FDR , 0.005) enriched GO terms for up-regulated genes in BL5+ compared to BL. The color in each cell represents the significance of enrichment based on the FDR value. The cells without color indicate not significantly enriched. The full list of GO terms is given in Supplemental Table S3.

including 102 up-regulated and 61 down-regulated DEGs (Fig. 2, B and C). Cell wall-related GO terms such as “glucan metabolic process,” “polysaccharide metabolic process,” and “polysaccharide biosynthetic process” were enriched (Fig. 3). The continuous up-regulation of cell wall-degrading and -modifying genes at MMC and MEI might be related to the observed cell wall abnormalities such as vacuolization and degradation at MEI (Fig. 1, G and H). In addition, 599 DEGs (9.54% of total DEGs) in “biotic stress” and “abiotic stress” and genes in various Plant Physiol. Vol. 174, 2017

The MapMan and GO results showed that cell wall metabolism, especially cell wall modification and degradation, were active in BL5+ (Fig. 2). Cell wall-modifying and -degrading genes, also called cell wall-remodeling genes, affect cell growth by changing the extensibility and physiochemical properties of the cell wall (Cosgrove, 1999; ZhiMing et al., 2011; Gou et al., 2012; Hepler et al., 2013; Xiao et al., 2014). We analyzed the expression of genes from seven main enriched families: expansins (EXPs), xyloglucan endotransglucosylase/hydrolases (XTHs), pectin methyl esterases (PMEs), pectin acetylesterases (PAEs), polygalacturonases (PGases), glycoside hydrolase family 9 (GH9), and arabinogalactan proteins (AGPs), which affect cell wall structure and extension by enzymatic or nonenzymatic ways (Bosch et al., 2005; Coimbra et al., 2008; Maris et al., 2009; Gou et al., 2012; Che et al., 2016). We found that 95, 85, and 23 genes from these families were differentially expressed at MMC, MEI, and AME, respectively, of which 73 genes were differentially expressed at both MMC and MEI, including 59 up-regulated and 14 down-regulated (Fig. 4, A and B). The 59 up-regulated genes consisted of 13 EXPs, 7 XTHs, 5 PMEs, 4 PAEs, 8 PGases, 7 GH9s, and 15 AGPs. OsEXP4 was elevated in BL5+ at MMC and MEI. In rice, overexpression of OsEXP4 results in increased cell length, while repression of OsEXP4 reduces cell length (Choi et al., 2003). Three up-regulated XTHs (Os06g48160, Os03g01800, and Os10g39840) have been characterized to utilize xyloglucan, one of the cell wall components, as substrate (Hara et al., 2014). OsGLU1, a GH9 member that regulates cell elongation in rice (Zhou et al., 2006; Zhang et al., 2012), was found up-regulated as well. The up-regulation of genes involved in cell wall degradation and cell size regulation was in accordance with cell wall abnormalities in BL5+ (Fig. 1, G and H). Thus, we inferred that the up-regulation of the cell wall-remodeling genes induced by ORF5+ in BL5+ led to cell structure change, damaging the cell wall integrity. Cellulose Deposition in BL5+

In BL5+, unknown substances were accumulated in the aborted ovule at AME (Fig. 1, I and J). The enrichment of the GO term “cellulose biosynthetic process” suggested a misregulation of cellulose biosynthesis and that the deposited substances in the aborted ovule 1687

Zhu et al.

Figure 4. Differential regulation of cell wall genes in BL5+ and BL3+5+ compared with BL. A, The number of differentially expressed genes involved in cell wall modification and degradation in BL5+ and BL3+5+. The x axis indicates the gene numbers. The y axis indicates the gene families at different stages. B, Heat map of expression change of genes involved in cell wall modification and degradation. The color in each cell represents the level of expression change based on the log2 fold. The cell without color indicates the gene was not differentially expressed. Full list of genes and their expressions is given in Supplemental Table S5. C, Expression change of cellulose biosynthetic genes.

might overaccumulate cellulose. To validate this, we used Pontamine fast scarlet 4B and calcofluor white to stain cellulose (Anderson et al., 2010; Thomas et al., 2012), and observed strong fluorescence signals in BL5+ at the site where the substances accumulated, indicating that cellulose was indeed in large quantity (Fig. 5, B and E), which was also supported by the much higher cellulose content in panicles of BL5+ than in BL (Fig. 5G). In contrast, normal mature embryo-sac structures were easily identified in the ovule of BL and BL3+5+, where low fluorescence intensity was detected (Fig. 5A, C, D, and F). Cellulose is a major cell wall polysaccharide composed of unbranched b-1, 4-linked glucan chains. In rice, OsMYB46 and OsSWNs are secondary cell wall biosynthetic regulators. OsSWN1 acts on the upstream of OsMYB46 by binding to the SNBE element in its promoter region. Overexpression of OsSWN1 or OsMYB46 in Arabidopsis results in the ectopic deposition of cellulose, xylan, and lignin (Zhong et al., 2011). In BL5+, OsSWN1, OsSWN4, OsSWN6, and OsMYB46 were up-regulated at MMC and MEI. Three 1688

cellulose synthases, OsCESA4, OsCESA7, and OsCESA9, which are regarded as subunits of cellulose synthase complex (Tanaka et al., 2003), were continuously up-regulated from MMC to AME. In addition, expression of five cellulose synthase-like genes (OsCslA1, OsCslA3, OsCslA5, OsCslA9, and OsCslF6) and two brittle culm genes (BC1 and BC3) were also elevated at MMC or MEI (Fig. 4C). Upstream transcription factors such as OsSWN1, OsSWN6, and OsMYB46 had peak expression at MEI, while their downstream genes such as OsCESAs and two brittle culm genes had the highest expression at AME (Supplemental Fig. S4), in agreement with the overaccumulation of cellulose at AME (Fig. 5, B and E). These results showed that the cellulose biosynthetic pathway was highly activated in BL5+, resulting in the accumulation of cellulose in the aborted embryo-sac. Callose Deposition in BL5+

Callose is a dynamic component of the plant cell wall that plays important roles in megasporogenesis. The Plant Physiol. Vol. 174, 2017

Transcriptomic Basis of a Rice Reproductive Barrier

Figure 5. Cellulose deposition in the mature ovules of BL, BL5+, and BL3+5+. A to F, Cellulose deposition in BL (A and D), BL5+ (B and E), and BL3+5+ (C and F). A to C, Pontamine fast scarlet 4B staining. D to F, Calcofluor white staining. The strong red or blue fluorescence indicates where the cellulose accumulated. Bars = 50 mm. G, Cellulose content in the panicles of BL and BL5+. The data are means 6 SE (n = 3). ***Significant difference at P , 0.001.

uneven distribution of the callose wall around the megaspores decides the polarization of megaspores, the formation of a functional megaspore, and the abortion of the other three nonfunctional megaspores (Ünal et al., 2013). Since many genes encoding b-1,3-glucanases, which have putative functions in callose metabolism, were differentially expressed in BL5+ at MMC and MEI (Supplemental Fig. S5), we investigated the dynamic callose distributions in ovule of BL, BL5+, and BL3+5+. During meiosis, callose preferentially accumulated at the micropylar end of megaspores. Two callose bands were observed in all three genotypes at the dyad stage (Fig. 6, A–F), and at the tetrad stage at least three callose bands could be recognized also in all three genotypes although slightly different in distribution patterns (Fig. 6, G–L). This apparent normal callose distribution in megaspores of BL5+ indicated likely normal meiosis. Unlike BL and BL3+5+, however, weak callose deposition was detected in nucellus cells at the tetrad stage in BL5+ (Fig. 6I). At the functional megaspore (FM) stage, callose accumulated in the three degenerated

gametes but not in the FM (Fig. 6, M–R). The survived FM subsequently developed to embryo-sac (Fig. 6, S, T, W, and X). However, the ovule of BL5+ displayed irregular callose deposition around the FM and in nucellus cells (Fig. 6, O and P), and finally, the aborted embryo-sac was stuffed with a large amount of callose (Fig. 6, U and V). Thus, ORF5+ resulted in a misregulation of callose deposition during megasporogenesis such that callose began to accumulate abnormally from MEI and was deposited in a large quantity at AME in the ovule of BL5+. ER Stress and Premature PCD Induced in BL5+

Under stress conditions, the number of misfolded proteins will increase and may exceed the capacity of the ER quality control system, thus inducing ER stress (Howell, 2013; Wan and Jiang, 2016). In BL5+, 78.95% of genes in “RNA transcription,” 66.23% of genes in “RNA processing,” 57.69% of genes in “protein synthesis initiation,” 58.33% of genes in “protein synthesis Figure 6. Callose deposition (green fluorescence) in the developing ovule of BL, BL5+, and BL3+5+ by aniline blue staining. A to F, The dyad stage. The red arrowheads indicate the callose bands around the dyad. G to L, The tetrad stage. The red arrowheads indicate the callose bands around the tetrad. M to R, The functional megaspore stage. The red arrowheads indicate the degenerated megaspores. S to X, The mature embryo-sac. I, O, and U, The gray arrowheads indicate abnormally deposited callose. Bars = 25 mm.

Plant Physiol. Vol. 174, 2017

1689

Zhu et al.

elongation,” 82.46% of genes in “ribosomal protein synthesis,” and 60.91% of genes in “ribosome biogenesis” were up-regulated at MMC (Supplemental Fig. S6), suggesting a greatly increased level of protein synthesis, which may increase the protein folding burden on ER. To detect whether ER stress occurred in BL5+ at MMC, we examined the expression of ER stress-responsive genes in BL5+ and found that 10 such genes were significantly up-regulated, including OsbZIP50, which is a key ER stress regulator in rice (Fig. 7A). The mRNA of OsbZIP50 was spliced when ER stress occurred (Hayashi et al., 2012). We then detected the unconventional splicing of OsbZIP50 mRNA in different materials and found that OsbZIP50 mRNA was spliced in BL5+ at MMC, but not in BL or BL3+5+ (Fig. 7B). The spliced OsbZIP50 mRNA encodes a transcription factor OsbZIP50-S, which binds to the pUPRE-II cis-element in the promoter regions of OsBIP4, ORF3, OsEBS, and OsPDIL2-3 (Hayashi et al., 2013; Takahashi et al., 2014; Wakasa et al., 2014). Accordingly, expression of these four genes was elevated (Fig. 7A). Some other genes including Fes1-like, CRT2-2, CRT1, and SAR1B-like, which were induced in ER stress, were also

up-regulated. These results showed that ER stress occurred in BL5+ at MMC relative to BL and BL3+5+. In response to ER stress, genes involved in protein modification (303 up-regulated and 358 down-regulated) and targeting (70 up-regulated and 53 down-regulated) were differentially expressed, presumably to help protein processing and transportation. In addition, 793 DEGs (339 upregulated and 454 down-regulated) involved in “protein degradation” were found, likely as a remedy to degrade the misfolded or unfolded proteins (Supplemental Fig. S6). The splicing of OsbZIP50 mRNA and up-regulation of ER stress-responsive genes were also detected at MEI and AME (Fig. 7, A and B), indicating that ER stress proceeded continuously in BL5+. The unresolved ER stress could induce premature programmed cell death (PCD) (Cai et al., 2014). In BL5+, PCD signals were detected from the meiosis stage (Yang et al., 2012), corresponding to the MEI stage in this study. At MEI, we detected up-regulation of two aspartic proteaseencoding genes, OsAP25 and OsAP37 (Fig. 7, A and C), both of which participate in the tapetal PCD process and induce PCD in yeast and plants (Niu et al., 2013).

Figure 7. Differential regulation of ER stress and PCD genes. A, Log2 fold change of the expression level of ER stress and PCD genes in BL5+ against BL. The color in each cell represents the levels of expression change based on the log2 fold. The cell without color indicates the gene was not differentially expressed. B, Splicing of OsbZIP50 mRNA in BL, BL5+, and BL3+5+ at each stage. C, qPCR verification for expression of tapetal PCD-related genes in BL and BL5+. The x axis represents different stages. The y axis represents the expression level relative to profilin (Os06g05880). The data are means 6 SE (n = 3). 1690

Plant Physiol. Vol. 174, 2017

Transcriptomic Basis of a Rice Reproductive Barrier

Figure 8. Analysis of genes differentially expressed in BL3+5+ compared to BL. A, Number of DEGs detected in the BL3+5+/BL and BL5+/BL comparisons at different stages. The overlapping areas indicate the genes differentially expressed in both BL3+5+ and BL5+. B, The number of DEGs in different functional categories of “Metabolism overview” and “Cellular response” by MapMan at MMC, MEI, and AME, respectively. C, Partial significantly (FDR , 0.005) enriched GO terms for up-regulated genes in BL3+5+ compared with BL. The color in each cell represents the significance of enrichment based on the FDR value. The cells without color indicate not significantly enriched. The full list of GO terms is given in Supplemental Table S4.

Overexpression of OsAP25 or OsAP37 in plants results in cell death, which is suppressed by aspartic proteasespecific inhibitor pepstatin A. Transformation of OsAP25 or OsAP37 into a PCD-deficient yeast strain restores the PCD ability of the strain. OsAP25 showed the highest expression at MEI (Fig. 7C), consistent with the observation of PCD signals at this stage. In the tapetal PCD process, transcription factor EAT1 directly regulates the expression of OsAP25 and OsAP37 (Niu et al., 2013). However, EAT1 was down-regulated at MEI, indicating that some other regulators might have replaced EAT1 to activate OsAP25 and OsAP37. TDR and TIP2, which work upstream of OsAP25 and OsAP37 in pollen development (Niu et al., 2013), were up-regulated in BL5+ at MMC and MEI (Fig. 7, A and D). The coordinated up-regulation of TDR, TIP2, OsAP25, and OsAP37 suggested that the ovule and tapetum may share some common genes/pathways in PCD. Taken together, an Plant Physiol. Vol. 174, 2017

OsbZIP50-dependent ER stress pathway and a premature PCD pathway involving the tapetal PCD genes were activated in BL5+. These results suggested that the ORF5+-induced up-regulation of cell wall degradation and modification genes in BL5+ caused cell wall damage, which induced both biotic and abiotic stresses. Lasting stress responses induced ER stress and unresolved ER stress led to premature PCD, eventually resulting in the ovule abortion.

Processes Detected in the BL3+5+/BL Comparison

BL3+5+ carrying both the killer ORF5+ and the protector ORF3+ showed no phenotypic difference from the wild-type BL with respect to fertility (Supplemental Fig. S1). However, 3,986 DEGs were still identified between 1691

Zhu et al.

integrity damage was suppressed. The cell wall compensation pathways such as cellulose biosynthesis and callose deposition were not activated (Figs. 4C and 6, K and Q), and ER stress and PCD did not occur (Fig. 7). As a result, the cell wall maintained integrity and the ovule meiosis was undisturbed (Fig. 1, L and M). At AME, the number of DEGs was further reduced to 370 in BL3+5+ compared with BL (Supplemental Table S2; Supplemental Data Set 2). These results indicated that ORF5+ induced similar but weaker damage processes in BL3+5+ as in BL5+ at MMC, but these processes were soon suppressed by ORF3+. Consequently, the embryo-sac developed normally, producing fertile female gametes (Fig. 1O).

DISCUSSION

Figure 9. A model for S5 locus-mediated ovule abortion.

BL3+5+ and BL at MMC (Supplemental Table S2; Supplemental Data Set 2), of which 3,843 DEGs were shared with the BL5+/BL comparison at this stage, indicating that ORF5+ induced similar processes in BL3+5+ (Fig. 8A). As expected, the process of “cell wall” was enriched with 217 DEGs (Fig. 8B), of which 64 (55 up-regulated and 9 down-regulated) were cell wall-modifying and -degrading genes from the seven main families (Fig. 4, A and B). GO terms “glucan metabolic process” and “polysaccharide metabolic process” were significantly enriched in up-regulated genes (Fig. 8C). Terms “biotic stress” and “abiotic stress” were also enriched involving 339 DEGs. However, the number of DEGs in each of the enriched functional categories was smaller than that in the BL5+/BL comparison (Fig. 4, A and B). Moreover, the splicing of OsbZIP50 mRNA and the up-regulation of its downstream genes were not detected (Fig. 7, A and B). These results suggested that ORF5+ induced similar processes in BL3+5+ at MMC, especially with respect to cell wall integrity damage, but to a much lesser extent, which did not trigger ER stress. At MEI, the difference between BL3+5+ and BL was greatly narrowed such that only 749 DEGs were identified, compared with 6,278 genes in the BL5+/BL comparison (Supplemental Table S2; Supplemental Data Set 2). Thus, numbers of DEGs, especially those in cell wall, biotic stress, and abiotic stress, were greatly reduced (Fig. 8B). Only seven DEGs from the seven cell wall families were found (Fig. 4, A and B), suggesting that cell wall 1692

The killer-protector system of S5 is composed of three proteins: cell wall protein ORF5+, transmembrane protein ORF4+, and ER resident protein ORF3+. Based on the comparative transcriptome results of BL5+ and BL3+5+, we proposed a model to summarize the processes underlying the S5 system (Fig. 9). In BL5+, many cell wall-modifying and -degrading genes in the families EXPs, XTHs, PMEs, PAEs, PGases, GH9, and AGPs were up-regulated at MMC and MEI, and many genes of these families play important roles in cell elongation or cell expansion (Zhou et al., 2006; Chen and Ye, 2007; Yang et al., 2007; Coimbra et al., 2008; Maris et al., 2009; Zhang et al., 2010, 2012; ZhiMing et al., 2011; Gou et al., 2012; Ma et al., 2013; Miedes et al., 2013; Hara et al., 2014; Xiao et al., 2014; Che et al., 2016). It has been suggested that postsynthetic modifications of the cell wall by cell wall-modifying genes may affect cell wall integrity (CWI; Pogorelko et al., 2013). Thus, large-scale up-regulation of cell wall-remodeling genes induced by ORF5+ may lead to changes of cell wall structures, which affect CWI, thus inducing stress response. ORF5+ encodes an extracellular aspartic protease. In Pichia pastoris and Saccharomyces cerevisiae, extracellular aspartic protease yapsins are required for the CWI maintenance (Krysan et al., 2005; Guan et al., 2012). Two plant aspartic proteases, A36 and A39, participate in cell wall remodeling in pollen tube (Gao et al., 2017). The study in Candida albicans shows that secreted aspartic proteinases use a set of cell wall proteins as substrates (Schild et al., 2011). We speculated that one or more substrates of ORF5+ existed in the rice cell wall, which induced cell wall damage after activation by ORF5+. In Saccharomyces cerevisiae, CWI damage triggers unfolded protein response, a homeostatic response to alleviate ER stress, which in turn affects the CWI (Krysan, 2009; Scrimale et al., 2009). In plants, pathogen invasion or pest attack can induce CWI damage and subsequent stress response (Bellincampi et al., 2014; Pogorelko et al., 2013). CWI damage eventually leads to PCD (Paparella et al., 2015). In BL5+, the functional categories responsive to pathogen invasion or pest attack were enriched, and ER stress was accordingly Plant Physiol. Vol. 174, 2017

Transcriptomic Basis of a Rice Reproductive Barrier

induced, presumably by CWI impairment, which in turn exacerbated the cell wall damage and its subsequent stress responses, thus finally leading to PCD. To alleviate the adverse effects caused by CWI damage, cells usually trigger a cell wall compensatory pathway to strengthen the cell wall (Hamann, 2012). For example, disruption by aspartic protease Yps7p in Pichia pastoris reduces chitin content in the lateral cell wall, but the cell thickens the inner cell wall for compensation (Guan et al., 2012). Similarly, in BL5+, key genes in secondary cell wall biosynthesis were significantly up-regulated, resulting in the cellulose accumulation to repair cell wall. Callose, a polysaccharide usually accumulated under stress conditions to provide physical barrier to the adverse environments (Pirselova and Matusikova, 2012), was deposited in large quantity. Despite other roles of callose in ovule development, the large amount deposition of callose in BL5+ was most likely to strengthen cell wall. ORF4+ acted as a partner of ORF5+ with a transmembrane localization. We speculated that ORF4+ functioned as a sensor of CWI damage. In plants, a well-studied cell wall sensor is AtWAK1, a wall-associated kinase, which resides in the plasma membrane and responds to oligogalacturonides derived from the hydrolysis of cell wall component (Brutus et al., 2010). Xa4, a wall-associated kinase in rice, contributes to disease resistance by strengthening the cell wall via promoting cellulose synthesis and suppressing cell wall loosening (Hu et al., 2017). Unlike wall-associated kinase, no extracellular binding domain or inner kinase domain was found in ORF4+; thus there may exist a coreceptor that interacts with ORF4+ in the plasma membrane and/or a downstream signal protein that interacts with plasma membrane receptor to transmit signals. Based on these analyses, it was deduced that ORF5+ damaged CWI by cleaving its substrates, leading to up-regulation of the cell wall-modifying and -degrading genes. ORF4+ sensed the CWI impairment and induced biotic and abiotic stress responses, and consequently induced ER stress due to the accumulation of drastically increased misfolded proteins in ER. The occurrence of ER stress in turn deteriorated the CWI damage and led to a more serious stress response. Finally, the unresolved cell wall damage, stress response, and ER stress triggered PCD. Meanwhile, the cellulose and callose accumulated in the ovule as the cells’ unsuccessful remedy to repair the damaged cell wall. BL3+5+ showed no difference with BL with respect to ovule development and fertility, but a number of cell wall-modifying and -degrading genes were still up-regulated in BL3+5+ at MMC, which may induce CWI damage. However, transcriptional change of these genes was suppressed at MEI and AME. Thus, the subsequent negative effects caused by CWI damage, including ER stress, PCD, cellulose, and callose deposition, were prevented in BL3+5+. ORF3+ encodes a HSP 70 chaperone, which is a downstream target of OsbZIP50s and induced by ER stress. ORF3+ is believed to help protein folding in ER to resolve ER stress like OsBIP1-OsBIP5 in the HSP70 gene family Plant Physiol. Vol. 174, 2017

(Yang et al., 2012; Hayashi et al., 2013; Wakasa et al., 2014). BIP protein alleviates Cd2+ stress-induced PCD (Xu et al., 2013). The tobacco (Nicotiana tabacum) plants with overexpression of BIP show delayed leaf senescence (Carvalho et al., 2014). In BL3+5+, ORF3+ suppressed the ORF5+-induced ER stress and PCD. It was speculated that one role of ORF3+ was alleviating ER stress thus not to induce premature PCD. The chaperone activity of ORF3+ may also help the folding of the signal proteins downstream of ORF5+ and ORF4+ so that the stress response may be suppressed at MEI. In conclusion, we revealed the processes underlying a reproductive barrier composed of ORF3, ORF4, and ORF5 at the transcriptomic level and proposed a mechanistic model based on the results. Nevertheless, further studies are needed to support and enrich the model.

MATERIALS AND METHODS Sample Preparation and Library Construction Balilla was a typical japonica rice (Oryza sativa) variety. The transgenic plant BalillaORF5+ was generated by introducing ORF5+ from an indica variety Nanjing 11 into Balilla under its native promoter (Chen et al., 2008). The plant BalillaORF3 +ORF5+ was isolated from the hybrid offspring between BalillaORF5+ and BalillaORF3+, which contained a transformed ORF3+ from Nanjing 11 under its native promoter (Yang et al., 2012). All the rice plants were grown in the fields of the experimental station of Huazhong Agricultural University, Wuhan, China. We determined the relationship between the floret length and stage of ovule development by the paraffin section observation of florets from the variety Balilla. The pistils from florets of Balilla, BalillaORF5+, and BalillaORF3+ORF5+ ranging in length from 3.8 to 4.2 mm, 5.5 to 6.0 mm, and 7.0 to 7.5 mm, respectively, were collected as MMC, MEI, and AME stage samples with two biological replicates. RNA was extracted using TRIzol reagent (Invitrogen) according to the manufacturer’s manual. RNA yields and quality were measured using Nanodrop 2000 (Thermo Scientific) and further quantified using Qubit2.0 Fluorometer (Invitrogen). RNA integrity was evaluated using Experion RNA Analysis Kits (Bio-Rad) and Experion Automated Electrophoresis Station (Bio-Rad). Each library was constructed using 3 mg total RNA using the TruSeq RNA Library Preparation Kit v2 (Illumina). The libraries were sequenced using Illumina HiSEquation 2000 sequencing platform at National Key Laboratory for Crop Genetic Improvement in Huazhong Agricultural University, Wuhan, China.

Data Processing and DEG Identification The data of the two replicates of all samples were collected and analyzed following the workflow described previously (Trapnell et al., 2012). In brief, the raw sequencing reads were mapped to the rice reference genome (MSU7.0, http://rice.plantbiology.msu.edu/index.shtml) using TopHat, then the transcripts were assembled and quantified using Cufflinks, and finally, differentially expressed genes were analyzed using Cuffdiff. The raw data files can be found in the NCBI database under the accession number GSE95200. Genes with $2-fold expression up-regulation or down-regulation and adjusted P value # 0.05 (counted using Cuffdiff) were regarded as up-regulated or down-regulated genes.

GO and MapMan Analysis GO enrichment was performed in AgriGO, a web-based tool for GO analysis (Du et al., 2010; http://bioinfo.cau.edu.cn/agriGO/index.php). The species “Oryza sativa” was chosen. GO terms in biological process groups with FDR , 0.005 were identified and analyzed further. The MapMan software version 3.5.1.R2 and pathways were downloaded from the MapMen Site of analysis (http://mapman. gabipd.org/web/guest/mapmanstore). The mapping information of rice genes (osa_MSU_v7_2015-09-08_mapping.txt.tar.gz) that can be imported to the MapMan software was downloaded from GoMapMan (http://www.gomapman. org/export/current/mapman; Ramsak et al., 2014). 1693

Zhu et al.

In Situ Hybridization The methods for preparation of paraffin embedded sections followed the previously described protocol (Weng et al., 2014). Primer pairs (listed in Supplemental Table S4) ORF3insituF and ORF3insituR, ORF4insituF and ORF4insituR were used to prepare probes for ORF3 and ORF4, respectively. Products amplified from ORF3+ or ORF4+ cDNA were ligated to pGEM-T vector (Promega). The sense probe was transcribed with T7 RNA polymerase, and the antisense probe was transcribed with SP6 RNA polymerase. The probes were labeled with digoxigenin (Roche). Hybridization and detection were performed as described previously (De Block and Debrouwer, 1993).

PCR and RT-qPCR Total RNA of pistils was isolated with TRIzol reagent (Invitrogen) and treated with RNase-free DNaseI (Invitrogen) for 20 min at room temperature to digest the contaminating DNA. The first-strand cDNA was synthesized using SSIII reverse transcriptase (Invitrogen) with Oligo(dT)15 Primer (Promega) following the manufacturers’ instructions. For amplification of OsbZIP50s, specific primers around the splice site were used. PCR conditions and PCR product detection were performed as described before (Yang et al., 2012). For quantitative real-time PCR, we carried out the reaction in a total volume of 10 mL containing 5 mL FastStart Universal SYBR Green Master (Rox; Roche), 1 mL cDNA sample, and 3 mM each primer on ABI ViiA7 system. The primers for PCR and RT-qPCR are listed in Supplemental Table S4. The gene profilin (Os06g05880) was used as the internal control in RT-qPCR (Yang et al., 2012).

Supplemental Figure S5. Expression change of genes involved in callose metabolism. Supplemental Figure S6. The number of DEGs involved in RNA transcription, protein synthesis, and processing. Supplemental Table S1. The number of expressed genes in BL, BL5+, and BL3+5+ at different stages. Supplemental Table S2. The number of differentially expressed genes identified in different comparisons. Supplemental Table S3. Full list of significantly enriched GO terms (FDR , 0.005) for DEGs in BL5+ versus BL. Supplemental Table S4. Full list of significantly enriched GO terms (FDR , 0.005) for DEGs in BL3+5+ versus BL. Supplemental Table S5. The list of cell wall-modifying or -degrading genes and their expression changes based on log2 fold. Supplemental Table S6. The primers used in this study. Supplemental Data Set 1. The differentially expressed genes in BL5+ compared with BL. Supplemental Data Set 2. The differentially expressed genes in BL3+5+ compared with BL.

ACKNOWLEDGMENTS

Histological Staining The fresh florets were fixed with FAA solution (50% ethanol, 5% acetic acid, and 3.7% formaldehyde) at 4°C overnight, dehydrated with ethanol in a series of concentrations (50, 70, 85, 95, and 100%), infiltrated with xylene, and finally embedded in paraffin. The materials were cut into 10-mm-thick sections and then washed in a xylene/ethanol series (0–50% ethanol). For callose staining, the rehydrated sections were stained with 0.1% aniline blue solution for 60 min. For cellulose staining, the rehydrated sections were stained with 0.01% Pontamine Fast Scarlet 4B solution for 20 min or 0.1% calcofluor white (Sigma-Aldrich) for 5 min. The stained sections were scanned under fluorescence microscopy. The florets for observation of ovule development process were stained with ehrlich hematoxylin. After fixation, the materials were rehydrated, infiltrated with ehrlich hematoxylin for 4 d, then dehydrated with ethanol, infiltrated with xylene, and embedded in paraffin. The 10-mm-thick sections were dewaxed with xylene and sealed with Canada balsam (Sigma-Aldrich). The slides were scanned under light microscopy.

Cellulose Content Determination The cellulose content was determined as described previously (Guo et al., 2014). In brief, the materials were suspended with potassium phosphate buffer (pH 7.0), chloroform-methanol (1:1, v/v), dimethylsulphoxide-water (9:1, v/v), 0.5% (w/v) ammonium oxalate, and 4 M KOH containing 1.0 mg/mL sodium borohydride in turn. The remaining pellet was collected and extracted further with H2SO4 (67%, v/v), and the supernatants were collected for the determination of cellulose.

Accession Numbers Sequence data from this article can be found in the GenBank libraries under accession numbers JX138499.1 (ORF31), JX138498.1 (ORF32), JX138502.1 (ORF41), JX138503.1 (ORF42), EU889295 (ORF51), and EU889294.1 (ORF52).

Supplemental Data The following supplemental materials are available. Supplemental Figure S1. The fertility of BL, BL5+, and BL3+5+. Supplemental Figure S2. RNA in situ hybridization of ORF3+ and ORF4+. Supplemental Figure S3. DEGs in “Biotic stress” in BL5+ compared with BL at MMC and MEI. Supplemental Figure S4. qPCR verification for expression of the cellulose biosynthetic genes in pistils of BL, BL5+, and BL3+5+ at MMC, MEI, and AME, respectively. 1694

We thank Yihua Zhou (Institute of Genetics and Developmental Biology) for kindly providing the Pontamine Fast Scarlet 4B solution. Received January 25, 2017; accepted May 4, 2017; published May 8, 2017.

LITERATURE CITED Anderson CT, Carroll A, Akhmetova L, Somerville C (2010) Real-time imaging of cellulose reorientation during cell wall expansion in Arabidopsis roots. Plant Physiol 152: 787–796 Bellincampi D, Cervone F, Lionetti V (2014) Plant cell wall dynamics and wall-related susceptibility in plant-pathogen interactions. Front Plant Sci 5: 228 Bosch M, Cheung AY, Hepler PK (2005) Pectin methylesterase, a regulator of pollen tube growth. Plant Physiol 138: 1334–1346 Brutus A, Sicilia F, Macone A, Cervone F, De Lorenzo G (2010) A domain swap approach reveals a role of the plant wall-associated kinase 1 (WAK1) as a receptor of oligogalacturonides. Proc Natl Acad Sci USA 107: 9452–9457 Cai YM, Yu J, Gallois P (2014) Endoplasmic reticulum stress-induced PCD and caspase-like activities involved. Front Plant Sci 5: 41 Carvalho HH, Silva PA, Mendes GC, Brustolini OJ, Pimenta MR, Gouveia BC, Valente MA, Ramos HJ, Soares-Ramos JR, Fontes EP (2014) The endoplasmic reticulum binding protein BiP displays dual function in modulating cell death events. Plant Physiol 164: 654–670 Che J, Yamaji N, Shen RF, Ma JF (2016) An Al-inducible expansin gene, OsEXPA10 is involved in root cell elongation of rice. Plant J 88: 132–142 Chen J, Ding J, Ouyang Y, Du H, Yang J, Cheng K, Zhao J, Qiu S, Zhang X, Yao J, et al (2008) A triallelic system of S5 is a major regulator of the reproductive barrier and compatibility of indica-japonica hybrids in rice. Proc Natl Acad Sci USA 105: 11436–11441 Chen L, Ye D (2007) Roles of pectin methylesterases in pollen-tube growth. J Integr Plant Biol 49: 94–98 Choi D, Lee Y, Cho HT, Kende H (2003) Regulation of expansin gene expression affects growth and development in transgenic rice plants. Plant Cell 15: 1386–1398 Coimbra S, Jones B, Pereira LG (2008) Arabinogalactan proteins (AGPs) related to pollen tube guidance into the embryo sac in Arabidopsis. Plant Signal Behav 3: 455–456 Coolen S, Proietti S, Hickman R, Davila Olivas NH, Huang P-P, Van Verk MC, Van Pelt JA, Wittenberg AHJ, De Vos M, Prins M, et al (2016) Transcriptome dynamics of Arabidopsis during sequential biotic and abiotic stresses. Plant J 86: 249–267 Cosgrove DJ (1999) Enzymes and other agents that enhance cell wall extensibility. Annu Rev Plant Physiol Plant Mol Biol 50: 391–417 Plant Physiol. Vol. 174, 2017

Transcriptomic Basis of a Rice Reproductive Barrier

Coyne J, Orr H (2004) Speciation. Sinauer Associates, Sunderland, MA De Block M, Debrouwer D (1993) RNA-RNA in situ hybridization using digoxigenin-labeled probes: the use of high-molecular-weight polyvinyl alcohol in the alkaline phosphatase indoxyl-nitroblue tetrazolium reaction. Anal Biochem 215: 86–89 Digel B, Pankin A, von Korff M (2015) Global transcriptome profiling of developing leaf and shoot apices reveals distinct genetic and environmental control of floral transition and inflorescence development in barley. Plant Cell 27: 2318–2334 Dixon RA, Paiva NL (1995) Stress-induced phenylpropanoid metabolism. Plant Cell 7: 1085–1097 Dobzhansky T (1937) Genetics and the origin of species. Columbia University Press, New York Du H, Ouyang Y, Zhang C, Zhang Q (2011) Complex evolution of S5, a major reproductive barrier regulator, in the cultivated rice Oryza sativa and its wild relatives. New Phytol 191: 275–287 Du Z, Zhou X, Ling Y, Zhang Z, Su Z (2010) agriGO: a GO analysis toolkit for the agricultural community. Nucleic Acids Res 38: W64-70 Gao H, Zhang Y, Wang W, Zhao K, Liu C, Bai L, Li R, Guo Y (2017) Two membrane-anchored aspartic proteases contribute to pollen and ovule development. Plant Physiol 173: 219–239 Gou JY, Miller LM, Hou G, Yu XH, Chen XY, Liu CJ (2012) Acetylesterasemediated deacetylation of pectin impairs cell elongation, pollen germination, and plant reproduction. Plant Cell 24: 50–65 Groszmann M, Gonzalez-Bayon R, Lyons RL, Greaves IK, Kazan K, Peacock WJ, Dennis ES (2015) Hormone-regulated defense and stress response networks contribute to heterosis in Arabidopsis F1 hybrids. Proc Natl Acad Sci USA 112: E6397–E6406 Guan B, Lei J, Su S, Chen F, Duan Z, Chen Y, Gong X, Li H, Jin J (2012) Absence of Yps7p, a putative glycosylphosphatidylinositol-linked aspartyl protease in Pichia pastoris, results in aberrant cell wall composition and increased osmotic stress resistance. FEMS Yeast Res 12: 969–979 Guo K, Zou W, Feng Y, Zhang M, Zhang J, Tu F, Xie G, Wang L, Wang Y, Klie S, et al (2014) An integrated genomic and metabolomic framework for cell wall biology in rice. BMC Genomics 15: 596 Hamann T (2012) Plant cell wall integrity maintenance as an essential component of biotic stress response mechanisms. Front Plant Sci 3: 77 Hara Y, Yokoyama R, Osakabe K, Toki S, Nishitani K (2014) Function of xyloglucan endotransglucosylase/hydrolases in rice. Ann Bot (Lond) 114: 1309–1318 Hayashi S, Takahashi H, Wakasa Y, Kawakatsu T, Takaiwa F (2013) Identification of a cis-element that mediates multiple pathways of the endoplasmic reticulum stress response in rice. Plant J 74: 248–257 Hayashi S, Wakasa Y, Takahashi H, Kawakatsu T, Takaiwa F (2012) Signal transduction by IRE1-mediated splicing of bZIP50 and other stress sensors in the endoplasmic reticulum stress response of rice. Plant J 69: 946–956 Hepler PK, Rounds CM, Winship LJ (2013) Control of cell wall extensibility during pollen tube growth. Mol Plant 6: 998–1017 Howell SH (2013) Endoplasmic reticulum stress responses in plants. Annu Rev Plant Biol 64: 477–499 Hu K, Cao J, Zhang J, Xia F, Ke Y, Zhang H, Xie W, Liu H, Cui Y, Cao Y, et al (2017) Improvement of multiple agronomic traits by a disease resistance gene via cell wall reinforcement. Nat Plants 3: 17009 Ikehashi H, Araki H (1986) Genetics of F1 sterility in remote crosses of rice. In IRRI, SJ Banta, GS Argosino, eds. Rice Genetics. IRRI, Manila, pp 119– 130 Kawahara Y, de la Bastide M, Hamilton JP, Kanamori H, McCombie WR, Ouyang S, Schwartz DC, Tanaka T, et al (2013) Improvement of the Oryza sativa Nipponbare reference genome using next generation sequence and optical map data. Rice (NY) 6: 4 Kourmpetli S, Lee K, Hemsley R, Rossignol P, Papageorgiou T, Drea S (2013) Bidirectional promoters in seed development and related hormone/stress responses. BMC Plant Biol 13: 187 Krysan DJ (2009) The cell wall and endoplasmic reticulum stress responses are coordinately regulated in Saccharomyces cerevisiae. Commun Integr Biol 2: 233–235 Krysan DJ, Ting EL, Abeijon C, Kroos L, Fuller RS (2005) Yapsins are a family of aspartyl proteases required for cell wall integrity in Saccharomyces cerevisiae. Eukaryot Cell 4: 1364–1374 Kubo T, Takashi T, Ashikari M, Yoshimura A, Kurata N (2016) Two tightly linked genes at the hsa1 locus cause both F1 and F2 hybrid sterility in rice. Mol Plant 9: 221–232 Plant Physiol. Vol. 174, 2017

Liu A, Zhang Q, Li H (1992) Location of a gene for wide-compatibility in the RFLP linkage map. Rice Genet Newsl 9: 134–136 Liu K, Wang J, Li H, Xu C, Liu A, Li X, Zhang Q (1997) A genome-wide analysis of wide compatibility in rice and the precise location of the S5 locus in the molecular map. Theor Appl Genet 95: 809–814 Long Y, Zhao L, Niu B, Su J, Wu H, Chen Y, Zhang Q, Guo J, Zhuang C, Mei M, et al (2008) Hybrid male sterility in rice controlled by interaction between divergent alleles of two adjacent genes. Proc Natl Acad Sci USA 105: 18871–18876 Ma N, Wang Y, Qiu S, Kang Z, Che S, Wang G, Huang J (2013) Overexpression of OsEXPA8, a root-specific gene, improves rice growth and root system architecture by facilitating cell extension. PLoS One 8: e75997 Maris A, Suslov D, Fry SC, Verbelen JP, Vissenberg K (2009) Enzymic characterization of two recombinant xyloglucan endotransglucosylase/ hydrolase (XTH) proteins of Arabidopsis and their effect on root growth and cell wall extension. J Exp Bot 60: 3959–3972 Miedes E, Suslov D, Vandenbussche F, Kenobi K, Ivakov A, Van Der Straeten D, Lorences EP, Mellerowicz EJ, Verbelen JP, Vissenberg K (2013) Xyloglucan endotransglucosylase/hydrolase (XTH) overexpression affects growth and cell wall mechanics in etiolated Arabidopsis hypocotyls. J Exp Bot 64: 2481–2497 Mizuta Y, Harushima Y, Kurata N (2010) Rice pollen hybrid incompatibility caused by reciprocal gene loss of duplicated genes. Proc Natl Acad Sci USA 107: 20417–20422 Muller HJ (1942) Isolating mechanisms, evolution, and temperature. Biol Symp 6: 71–125 Niu N, Liang W, Yang X, Jin W, Wilson ZA, Hu J, Zhang D (2013) EAT1 promotes tapetal cell death by regulating aspartic proteases during male reproductive development in rice. Nat Commun 4: 1445 Ouyang Y, Chen J, Ding J, Zhang Q (2009) Advances in the understanding of inter-subspecific hybrid sterility and wide-compatibility in rice. Chin Sci Bull 54: 2332–2341 Ouyang Y, Li G, Mi J, Xu C, Du H, Zhang C, Xie W, Li X, Xiao J, Song H, et al (2016) Origination and establishment of a trigenic reproductive isolation system in rice. Mol Plant 9: 1542–1545 Ouyang Y, Zhang Q (2013) Understanding reproductive isolation based on the rice model. Annu Rev Plant Biol 64: 111–135 Paparella S, Tava A, Avato P, Biazzi E, Macovei A, Biggiogera M, Carbonera D, Balestrazzi A (2015) Cell wall integrity, genotoxic injury and PCD dynamics in alfalfa saponin-treated white poplar cells highlight a complex link between molecule structure and activity. Phytochemistry 111: 114–123 Pirselova B, Matusikova I (2012) Callose: the plant cell wall polysaccharide with multiple biological functions. Acta Physiol Plant 35: 635–644 Pogorelko G, Lionetti V, Bellincampi D, Zabotina O (2013) Cell wall integrity: targeted post-synthetic modifications to reveal its role in plant growth and defense against pathogens. Plant Signal Behav 8: e25435 Qiu SQ, Liu K, Jiang JX, Song X, Xu CG, Li XH, Zhang Q (2005) Delimitation of the rice wide compatibility gene S5 ( n ) to a 40-kb DNA fragment. Theor Appl Genet 111: 1080–1086  Baebler Š, Rotter A, Korbar M, Mozetic I, Usadel B, Gruden K Ramsak Z, (2014) GoMapMan: integration, consolidation and visualization of plant gene annotations within the MapMan ontology. Nucleic Acids Res 42: D1167–D1175 Schild L, Heyken A, de Groot PW, Hiller E, Mock M, de Koster C, Horn U, Rupp S, Hube B (2011) Proteolytic cleavage of covalently linked cell wall proteins by Candida albicans Sap9 and Sap10. Eukaryot Cell 10: 98–109 Scrimale T, Didone L, de Mesy Bentley KL, Krysan DJ (2009) The unfolded protein response is induced by the cell wall integrity mitogenactivated protein kinase signaling cascade and is required for cell wall integrity in Saccharomyces cerevisiae. Mol Biol Cell 20: 164–175 Seehausen O, Butlin RK, Keller I, Wagner CE, Boughman JW, Hohenlohe PA, Peichel CL, Saetre G-P, Bank C, Brännström A, et al (2014) Genomics and the origin of species. Nat Rev Genet 15: 176–192 Song X, Qiu SQ, Xu CG, Li XH, Zhang Q (2005) Genetic dissection of embryo sac fertility, pollen fertility, and their contributions to spikelet fertility of intersubspecific hybrids in rice. Theor Appl Genet 110: 205–211 Takahashi H, Wang S, Hayashi S, Wakasa Y, Takaiwa F (2014) Cis-element of the rice PDIL2-3 promoter is responsible for inducing the endoplasmic reticulum stress response. J Biosci Bioeng 117: 620–623 Tanaka K, Murata K, Yamazaki M, Onosato K, Miyao A, Hirochika H (2003) Three distinct rice cellulose synthase catalytic subunit genes required for cellulose synthesis in the secondary wall. Plant Physiol 133: 73–83 1695

Zhu et al.

Thimm O, Bläsing O, Gibon Y, Nagel A, Meyer S, Krüger P, Selbig J, Müller LA, Rhee SY, Stitt M (2004) MAPMAN: a user-driven tool to display genomics data sets onto diagrams of metabolic pathways and other biological processes. Plant J 37: 914–939 Tholl D (2006) Terpene synthases and the regulation, diversity and biological roles of terpene metabolism. Curr Opin Plant Biol 9: 297–304 Thomas J, Ingerfeld M, Nair H, Chauhan SS, Collings DA (2012) Pontamine fast scarlet 4B: a new fluorescent dye for visualising cell wall organisation in radiata pine tracheids. Wood Sci Technol 47: 59–75 Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, Rinn JL, Pachter L (2012) Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc 7: 562–578 Ünal M, Vardar F, Aytürk Ö (2013) Callose in plant sexual reproduction. In M Silva-Opps, ed, Current Progress in Biological Research. InTech, Croatia, pp 319–343 Vogel JT, Zarka DG, Van Buskirk HA, Fowler SG, Thomashow MF (2005) Roles of the CBF2 and ZAT12 transcription factors in configuring the low temperature transcriptome of Arabidopsis. Plant J 41: 195–211 Wakasa Y, Oono Y, Yazawa T, Hayashi S, Ozawa K, Handa H, Matsumoto T, Takaiwa F (2014) RNA sequencing-mediated transcriptome analysis of rice plants in endoplasmic reticulum stress conditions. BMC Plant Biol 14: 101 Wan S, Jiang L (2016) Erratum to: Endoplasmic reticulum (ER) stress and the unfolded protein response (UPR) in plants. Protoplasma 253: 765 Wang GW, He YQ, Xu CG, Zhang Q (2005) Identification and confirmation of three neutral alleles conferring wide compatibility in inter-subspecific hybrids of rice (Oryza sativa L.) using near-isogenic lines. Theor Appl Genet 111: 702–710 Wang J, Liu KD, Xu CG, Li XH, Zhang Q (1998) The high level of widecompatibility of variety ‘Dular’ has a complex genetic basis. Theor Appl Genet 97: 407–412 Wang Q, Wan L, Li D, Zhu L, Qian M, Deng M (2009) Searching for bidirectional promoters in Arabidopsis thaliana. BMC Bioinformatics (Suppl) 10: S29 Weng X, Wang L, Wang J, Hu Y, Du H, Xu C, Xing Y, Li X, Xiao J, Zhang Q (2014) Grain number, plant height, and heading date7 is a central regulator of growth, development, and stress response. Plant Physiol 164: 735–747 Winkel-Shirley B (2002) Biosynthesis of flavonoids and effects of stress. Curr Opin Plant Biol 5: 218–223

1696

Xiao C, Somerville C, Anderson CT (2014) POLYGALACTURONASE INVOLVED IN EXPANSION1 functions in cell elongation and flower development in Arabidopsis. Plant Cell 26: 1018–1035 Xu H, Xu W, Xi H, Ma W, He Z, Ma M (2013) The ER luminal binding protein (BiP) alleviates Cd(2+)-induced programmed cell death through endoplasmic reticulum stress-cell death signaling pathway in tobacco cells. J Plant Physiol 170: 1434–1441 Yamagata Y, Yamamoto E, Aya K, Win KT, Doi K, Sobrizal, Ito T, Kanamori H, Wu J, Matsumoto T, et al (2010) Mitochondrial gene in the nuclear genome induces reproductive barrier in rice. Proc Natl Acad Sci USA 107: 1494–1499 Yang J, Sardar HS, McGovern KR, Zhang Y, Showalter AM (2007) A lysine-rich arabinogalactan protein in Arabidopsis is essential for plant growth and development, including cell division and expansion. Plant J 49: 629–640 Yang J, Zhao X, Cheng K, Du H, Ouyang Y, Chen J, Qiu S, Huang J, Jiang Y, Jiang L, et al (2012) A killer-protector system regulates both hybrid sterility and segregation distortion in rice. Science 337: 1336–1340 Yu Y, Zhao Z, Shi Y, Tian H, Liu L, Bian X, Xu Y, Zheng X, Gan L, Shen Y, et al (2016) Hybrid sterility in rice (Oryza sativa L.) involves the tetratricopeptide repeat domain containing protein. Genetics 203: 1439–1451 Zhang GY, Feng J, Wu J, Wang XW (2010) BoPMEI1, a pollen-specific pectin methylesterase inhibitor, has an essential role in pollen tube growth. Planta 231: 1323–1334 Zhang JW, Xu L, Wu YR, Chen XA, Liu Y, Zhu SH, Ding WN, Wu P, Yi KK (2012) OsGLU3, a putative membrane-bound endo-1,4-betaglucanase, is required for root cell elongation and division in rice (Oryza sativa L.). Mol Plant 5: 176–186 ZhiMing Y, Bo K, XiaoWei H, ShaoLei L, YouHuang B, WoNa D, Ming C, Hyung-Taeg C, Ping W (2011) Root hair-specific expansins modulate root hair elongation in rice. Plant J 66: 725–734 Zhong R, Lee C, McCarthy RL, Reeves CK, Jones EG, Ye ZH (2011) Transcriptional activation of secondary wall biosynthesis by rice and maize NAC and MYB transcription factors. Plant Cell Physiol 52: 1856– 1871 Zhou HL, He SJ, Cao YR, Chen T, Du BX, Chu CC, Zhang JS, Chen SY (2006) OsGLU1, a putative membrane-bound endo-1,4-beta-D-glucanase from rice, affects plant internode elongation. Plant Mol Biol 60: 137–151

Plant Physiol. Vol. 174, 2017