EAT1 transcription factor, a non-cell-autonomous regulator of pollen production, activates meiotic small RNA biogenesis in rice anther tapetum

The 24-nucleotides (nt) phased secondary small interfering RNA (phasiRNA) is a unique class of plant small RNAs abundantly expressed in monocot anthers at early meiosis. Previously, 44 intergenic regions were identified as the loci for longer precursor RNAs of 24-nt phasiRNAs (24-PHASs) in the rice genome. However, the regulatory mechanism that determines spatiotemporal expression of these RNAs has remained elusive. ETERNAL TAPETUM1 (EAT1) is a basic-helix-loop-helix (bHLH) transcription factor indispensable for induction of programmed cell death (PCD) in postmeiotic anther tapetum, the somatic nursery for pollen production. In this study, EAT1-dependent non-cell-autonomous regulation of male meiosis was evidenced from microscopic observation of the eat1 mutant, in which meiosis with aberrantly decondensed chromosomes was retarded but accomplished somehow, eventually resulting in abortive microspores due to an aberrant tapetal PCD. EAT1 protein accumulated in tapetal-cell nuclei at early meiosis and postmeiotic microspore stages. Meiotic EAT1 promoted transcription of 24-PHAS RNAs at 101 loci, and importantly, also activated DICER-LIKE5 (DCL5, previous DCL3b in rice) mRNA transcription that is required for processing of double-stranded 24-PHASs into 24-nt lengths. From the results of the chromatin-immunoprecipitation and transient expression analyses, another tapetum-expressing bHLH protein, TDR INTERACTING PROTEIN2 (TIP2), was suggested to be involved in meiotic small-RNA biogenesis. The transient assay also demonstrated that UNDEVELOPED TAPETUM1 (UDT1)/bHLH164 is a potential interacting partner of both EAT1 and TIP2 during early meiosis. This study indicates that EAT1 is one of key regulators triggering meiotic phasiRNA biogenesis in anther tapetum, and that other bHLH proteins, TIP2 and UDT1, also play some important roles in this process. Spatiotemporal expression control of these bHLH proteins is a clue to orchestrate precise meiosis progression and subsequent pollen production non-cell-autonomously.


Introduction
Small noncoding RNAs are 20-30 nucleotides (nt) long and associate with Argonaute family proteins to serve as guide molecules for RNA silencing in various biological processes, such as cell type specification, cell proliferation, cell death, metabolic control, transposon silencing and antiviral defense [1]. Plant genomes encode precursors of microRNA (miRNA) and small interfering RNA (siRNA), as do animal genomes [2]. miRNA is produced from a hairpin structure of a single precursor RNA molecule, and siRNA is derived from a precursor RNA that is either naturally double-stranded or is formed by RNA-dependent RNA polymerases.
The third class of animal small RNAs is Piwi-interacting RNA (piRNA). The piRNA is abundantly expressed in the germline and acts in silencing of transposable elements (TEs) [3], massive elimination of paternally derived mRNAs [4], systemic recognition of self and nonself mRNAs [5,6], and so on. piRNA associates with Piwi family proteins, a distinct subgroup of Argonaute proteins. In contrast, plants have no Piwi family Argonautes [7,8], and consequently lack piRNA species. In place of piRNA, trans-acting siRNA (tasiRNA) and phased secondary siRNA (phasiRNA) are identified as plant-specific small RNA subgroups. In monocot model plants, rice and maize, phasiRNAs are abundantly expressed in the male reproductive organs, and in this study, the term "phasiRNA" will be used for monocot reproductive phasiR-NAs derived from protein-noncoding regions. Both tasiRNA and phasiRNA are produced via miRNA-dependent primary processing, and characterized by phased alignment on both sense and antisense strands in genomic regions. However, they are distinct in several points. First, phasiRNAs are abundantly expressed in developing reproductive organs [9][10][11][12][13], while 21-nt tasiRNAs are expressed in both vegetative and reproductive phases [14]. Second, phasiRNAs are transcribed from hundreds or thousands of unique, namely nonrepetitive, intergenic regions [9,[11][12][13], while a few tasiRNA-producing (TAS) loci are conserved in the plant genome [15][16][17]. Finally, no phasiRNA targetting a protein-coding gene has been identified, whereas tasiRNAs are complementary to particular genes important for defense and developmental events [14]. In plant reproduction, 24-nt unphased siRNAs or 21-nt epigenetically activated siRNAs (easiRNAs) are thought to maintain genome integrity by programmed DNA methylation of TEs [18][19][20]. The roles of phasiRNAs during plant reproduction largely remain elusive.
In rice, a single-stranded PHAS precursor RNA is primarily processed with 22-nt miRNA triggers; miR2118 for 21-PHASs and miR2275 for 24-PHASs [10]. PHAS and TAS RNA members each have one or two conserved complementary sequences to miRNAs, and are cleaved via the one-hit or two-hit processing pathway; the one-hit mode is mediated by the AGO1-miRNA complex for 5'-end cleavage of precursor RNAs [15] to generate the 3' fragment that becomes double-stranded, and the two-hit mode depends on AGO1-or AGO7-miRNA, which potentially associates with both ends and cleaves either end or both [14]. The processed RNA is made double-stranded by RNA DEPENDENT RNA POLYMERASE6 (RDR6) [21], and chopped into 21-and 24-nt lengths by DICER-LIKE4 (DCL4) and DCL5 (previous DCL3b in rice), respectively [10].
The anther is a four-lobed male reproductive organ in angiosperms. Each anther lobe is composed of central sporogenous cells and four concentric somatic layers; the epidermis, endothecium, middle layer and tapetum, from outward to inward ( Fig 1A) [22][23][24]. Sporogenous cells undergo several rounds of mitosis and mature into pollen mother cells (PMCs) to prepare for meiosis [22][23][24]. Maize OUTER CELL LAYER4 (OCL4), an HD-ZIP IV transcription factor (TF), expressed in the anther epidermis and MALE STERILE23 (MS23), a basic helix-loop-helix (bHLH) TF expressed in the tapetum are required for 21 and 24-nt phasiRNA  Table 1  biogenesis, respectively [12,25]. Small RNA-mediated intercellular signaling is proposed in various steps of plant reproduction, for example, between sperm and vegetative cells in the pollen [19,26] and between megaspore mother cells and somatic nucellar cells in the ovule [18]. The intercellular movement of reproductive phasiRNAs has been proposed in maize [12,13], while there is yet no decisive evidence. The underlying mechanism to determine the spatiotemporal expression of reproductive phasiRNAs in anthers has largely remained elusive.
In this study, we focused on the rice bHLH TFs, because they are key transcriptional regulators for differentiation and development of anther somatic layers. TDR INTERACTING PRO-TEIN2 (TIP2)/bHLH142 is expressed in several undifferentiated cell layers to form the middle layer and tapetum [27,28]. TAPETUM DEGENERATION RETARDATION (TDR)/bHLH5 makes a heterodimer with TIP2 to promote tapetal differentiation [29]. ETERNAL TAPE-TUM1 (EAT1)/bHLH141, 41% similar to TIP2, also dimerizes with TDR, and activates transcription of aspartic protease-encoding genes to promote programmed cell death (PCD) of postmeiotic tapetal cells [30,31]. UNDEVELOPED TAPETUM1 (UDT1)/bHLH164 [32] is expected to function upstream of the regulatory cascade for anther wall development. However, downstream targets of these bHLH TFs are largely unknown.
In addition to its role in tapetal PCD, we found that EAT1 is required earlier in tapetal development to support meiosis, while the loss of EAT1 function has little impact on the tapetum morphology. EAT1 shows a bimodal expression at both early meiosis and postmeiosis. Interestingly, EAT1 expressed during early meiosis promoted both transcription and processing of 24-PHAS precursor RNAs to produce 24-nt phasiRNAs in tapetum. This study demonstrates that EAT1 is one of key regulators triggering meiotic phasiRNA biogenesis in anther tapetum, and that other bHLH proteins, TIP2 and UDT1, also play important roles in this process.

EAT1 is expressed in anther tapetum during early meiosis
To determine the impact of bHLH proteins in communication between somatic tapetal cells and PMCs in rice anthers, we first performed quantitative reverse-transcription PCR (qRT-PCR) of four bHLH genes: UDT1, TDR, TIP2 and EAT1, all of which are involved in tapetal cell-fate decision [27][28][29][30][31][32]. In this study, we separated anther developmental processes into six stages to characterize spatiotemporal expression of these genes (ST.1 to ST.6; Table 1). qRT-PCR of meiotic anthers demonstrated that UDT1, TDR and TIP2 were expressed as expected from previous reports (S1 Fig, S1 Data). However, EAT1 expression was bimodal, both at early meiosis (ST.2) and postmeiosis (ST.5), whereas it was previously thought to function only in postmeiotic tapetal PCD [30,31]. To investigate EAT1 expression during early meiosis, an EAT1pro-EAT1-GFP transcriptional fusion construct (Fig 1B) was introduced into male-sterile eat1-4 plants homozygous for a putative null allele with a Tos17-retrotransposon insertion (S2A- S2F Fig). The transgenic plants recovered male fertility (Fig 1C), indicating that the EAT1-GFP protein is functional in planta. EAT1-GFP expression was bimodal at ST.2 and ST.5, as was mRNA expression, and the two expression peaks were clearly separated by the silent ST.4 ( Fig 1D). Transcription of AP25, an aspartic protease gene required for tapetal PCD initiation [30], was fully dependent on EAT1 at ST.5 (S3 Fig, S1 Data), while no AP25 transcript was detected at ST.2 or ST.3. These results confirm that the role of meiotic EAT1 is distinct from its postmeiotic role in tapetal PCD and further suggest that the EAT1 bHLH TF has distinct bHLH partners at these two developmental stages. The eat1-4 mutant phenotype was remarkable in postmeiotic ST.5 and ST.6 anthers, in which tapetal cells were unusually degenerated at ST.5, concurrent with abortive microspores On the other hand, no morphological phenotype was found in earlier stages, ST.1 to ST.4 by light microscopy (Fig 2A and S4F-S4H Fig). Degradation of beta-1,4-glucan on cell walls of tapetal cells and PMCs occurred normally in eat1-4 anthers at ST.2-ST.3 stages (Fig 2B). These observations were largely consistent with previous results [30].

Delayed and asynchronous male meiosis in the eat1-4 mutant
We detected an unreported defect in male meiosis of eat1-4 mutants: PMCs harbor aberrantly decondensed bivalent chromosomes frequently, 74.2% at diakinesis (n = 70) and 69.4% at metaphase I (n = 36) (Fig 2C). In addition, two out of 27 eat1-4 PMCs at anaphase I harbored lagging chromosomes or chromosomal bridges, which were not found in the wild-type (n = 42) (Fig 2C). Another 5.5% eat1-4 PMCs exhibited interphase-like nuclei with fully decondensed chromosomes (n = 163), in contrast to wild-type PMCs (n = 192, Fig 2C). In addition, meiotic division timing was retarded in mutant anthers, with asynchronous progression within an anther lobe (S5 Fig, S1 Data). Despite these meiotic defects, male meiosis could complete, but resulting microspores were aborted most likely by the aberrant tapetum, which normally secretes nutrients and exine components required during post-meiotic pollen development (S4J Fig). These results suggest that non-cell-autonomous signaling or some nutrient delivery between somatic tapetal cells and PMCs is mediated by EAT1 during meiosis, in addition to post-meiosis.

EAT1 activates transcription of 101 loci encoding 24-PHAS RNAs
To identify genes under the control of meiotically expressed EAT1, we conducted mRNA-seq experiments using whole anther samples and compared the data between wild-type and eat1-4 plants. The data were obtained from three different meiotic stages: premeiosis (ST.1), early meiosis (ST.2) and late meiosis (ST.4), each with three biological replicates. 142,048,793 reads from wild-type and 146,928,874 reads from eat1-4 anthers (S1 Table) Table). The ontology terms for 7 of 115 genes were enriched in lipid metabolism based on the agriGO algorithm [33] (S3 Table), implying that they function in pollen coat formation [34].
mRNA-seq also identified 6,097 regions generating long intergenic noncoding RNAs (lincRNAs), and 248 showed ST.2-enriched and EAT1-dependent expression (Fig 3A, S4  Table). Next, we conducted small RNA-seq (sRNA-seq) to ask whether these lincRNAs are small RNA precursors or not. 52,726,712 reads of total small RNAs extracted from wild-type and 62,364,061 from eat1-4 anthers were mapped onto the rice genome (S1 Table). As a result, the 93 lincRNAs were defined as 24-PHAS RNAs, because a large number of 24-nt small RNAs were mapped in a 24-nt phasing manner on the lincRNA loci (see below for details). Of 44 24-PHAS loci previously reported [9,10], 24 were included in the loci identified in this study. Another 8 loci, which were left out of our first selection by their length or overlapping coding genes, generated EAT1-dependent and ST.2-enriched 24-nt phasiRNAs (S4 Table), while the remaining 12 loci did not. Thus, adding the 8 loci, a total of 101 loci were specified as ST.2-enriched and EAT1-dependent 24-PHAS loci and analyzed hereafter.
Median FPKM values of 24-PHAS transcripts detected at the 101 loci in wild-type ST.2 anthers were 688-fold and 24-fold higher than those in ST.1 and ST.4 anthers, respectively. In addition, the values were 55-fold higher than in eat1-4 anthers at ST.2 (Fig 3B and 3C, S4 Table). This result reconfirmed the EAT1-dependent and early meiosis-enriched nature of 24-PHAS transcripts. This trend was reproducible in qRT-PCR of five 24-PHASs (Fig 3D, S1  Data). In contrast, most 24-nt RNAs from the corresponding PHAS loci were abundant not only in ST.2, but also in ST.4 anthers (Fig 3B and 3C, S4 Table), implying slower turnover of small RNAs than precursor transcripts.
The 101 PHAS loci were unevenly distributed in the genome as reported previously [9], except for chromosomes 1 and 9, and many loci formed several clusters on each chromosome ( Fig 4A, S4 Table). Sequence comparison by the MEME program [35] demonstrated that 93  Methods), and frequencies of repetitive sequences including TEs (gray charts). The horizontal length of each box corresponds to the physical distance of respective rice chromosomes. (B) A conserved sequence logo found in upstream of ninetythree 24-PHAS loci detected by MEME program [35], which are potentially targeted by miR2275. The arrow indicates the predicted cleaved position by DCL1 and miR2275 complex [10]. (C) Frequency of repetitive sequence (grey), gene coding region (blue) and miR2275 targeted site (red) around 24-PHAS loci. The data was examined in 93 24-PHAS transcripts with conserved miR2275 targeted sites. The reason why a small peak of miR2275 target site appeared at the 3' end of 24-PHAS is that some 24-PHAS loci were relatively small in length (~500 bp). (D) Characterization of three 24-PHAS loci. From the top to the bottom, the graphs indicate the mapping results of mRNA-seq and 24-nt sRNA-seq reads (gray histograms), the 24-nt phasing pattern (green and orange charts), and the plot of read counts from the degradome-seq using young panicles of indica variety, 93-11 [38]. The degradome analysis revealed that the cleavage of three 24-PHAS transcripts frequently occurs at the position shown in (B), within the predicted miR2275 sites (red dots), while few degradome-seq reads were mapped onto both sense and antisense strands of other regions (gray dots). Reads were depicted by IGV [78]. (E) An example of distribution of EAT1-dependent 24-PHAS-loci cluster (green boxes) on the long arm of chromosome 12, with the context of surrounding genes (blue) and repetitive sequences (black).  Table). The miR2275 sites were conserved at the 5'-region in 92 loci (Fig 4C, S4  Table), consistent with previous results that 22-mer miRNA triggers one-hit processing [36,37]. The phased pattern tended to start at the 13th position in the 22-mer miR2275 site in most of 24-PHAS loci ( Fig 4D). This position corresponded to the cleavage site of the AGO1/ miR2275 complex reported previously [10]. Consistent with this, the degradome data from the indica rice variety [38] demonstrated that the cleavage actually occurred at the same position relative to the miR2275 complementarity in 62 of 93 24-PHAS loci (Fig 4B and 4D, S4 Table), and that almost of lincRNAs detected here were the unprocessed, primary 24-PHAS RNAs.
Of 24-nt small RNAs mapped on 93 24-PHAS loci, the 77.1% reads from wild-type ST.2 and ST.4 anthers showed a 24-nt phased pattern which starts from putative AGO1/miR2275 cleavage site (S6 Fig, S1 Data), indicating that 24-nt small RNAs produced from these loci were processed by DCL5.
Most 24-PHAS loci were mapped to unique or low copy regions (Fig 4C and 4E, S1 Data).  Table). We concluded that meiotic 24-nt phasiRNAs originate from 101 intergenic 24-PHAS loci and that they have a role distinct from TE silencing.
The expression of 24-PBR genes other than DCL5 was examined. DCL1 and RDR6 are respectively required for processing of miR2275 precursors and RNA double-strand formation [10,21]. DCL1 and RDR6 transcripts were abundant in ST.2 anthers; however, both were also abundant in ST.1 and ST.4 anthers and were unaffected by the eat1-4 mutation (S8A Fig, S1 Data), indicating that expression of DCL1 and RDR6 is EAT1 independent and not restricted to meiotic stages. Transcripts of pri-miR2275a/b, the precursors of mature miR2275, were enriched in ST.2 anthers. In contrast to 24-PHASs and DCL5, the amount of pri-miR2275 transcripts was elevated in the eat1-4 mutant (S8A Fig, S1 Data). pri-miR2275b promoter sequences were not enriched in ChIP of EAT1-GFP-expressing anthers, despite containing E-box motifs (S8C and S8D Fig, S1 Data).
To investigate the EAT1 ability to promote the transcription of 24-PHAS and DCL5 loci, we performed the transient expression assay. The bHLH proteins have homo-and heterodimerization ability [42]. Thus, the effector construct encoding any two of EAT1, TIP2, UDT1 and TDR was cotransfected with the 24-PHAS or DCL5 promoter (pPHAS, pDCL5)-Luciferase fusion reporter into rice protoplasts (S10A Fig), and the promoter activity was measured. The activity of two pPHASs was significantly 4.46 (chr5-20) and 3.99-fold (chr6-97) elevated in EAT1-UDT1 cotransfection, compared to the no effector control (Fig 6A). However, contrary to expectations, the same combination displayed insignificant effects on the pDCL5 (Fig 6A). Little effect on pPHASs nor pDCL5 was observed in the transfection of EAT1 alone and EAT1-TIP2, while the EAT1-TDR cotransfection slightly affected the activity of pPHASs (1.85 and 2.17 fold) and pDCL5 (1.95 fold) ( Fig 6A). Interestingly, EAT1-UDT1 cotransfection induced the pEAT1 activity by greater 7.61 fold (S10B Fig), while it was slightly upregulated by the EAT1-TDR cotransfection (1.58 fold). Cotransfection of EAT1 with TIP2, TDR or UDT1 displayed no significant effect on the pDCL3a (S10B Fig). To examine the protein-protein interaction between EAT1 and UDT1, we performed the bimolecular fluorescence complementation analysis (BiFC) in rice protoplasts. EAT1 fused with the C-terminal split of YFP (EAT1-cYFP and cYFP-EAT1) gave positive BiFC signals when coexpressed with UDT1-nYFP (Fig 6B, S11A and S11B Fig), while they tended to be detectable faintly in the nucleus (Fig 6B, S11A Fig arrows) or intensely in the cytoplasm ( Fig  6B, S11A Fig arrowhead). In both cases, the positive signals were always more intense compared to negative controls (Fig 6B, S11A-S11C Fig).
The above results demonstrate that the meiotic EAT1 TF promotes the transcription of 24-PHAS precursors and the EAT1 gene itself by interacting with UDT1 at the molecular level. EAT1 also promotes the DCL5 transcription, but likely with an unknown bHLH partner.
Collectively, these results suggest that TIP2 has the potential to activate transcription of both 24-PHASs and DCL5 by interacting with UDT1 at the molecular level in early meiosis.
The mapping mode of 24-nt masiRNAs was shown in two 24-PHAS loci for example (chr12-82 and chr12-85, Fig 7C). On the chr12-82 locus, 165 and 207 reads of only a 24-nt masiRNA species (masiRNA_u_0815) were mapped at the third phase of the sense strand in ST.2 and ST.4 anthers, respectively (Fig 7C left). A significant reduction of the masiR-NA_u_0815 in male-sterile eat1-4 plants (Fig 7C, S5 Table) confirmed their origin in anthers, not in pistils, although MEL1 is expressed in both male and female cells [7]. A similar tendency was found in the chr12-85 and masiRNA_u_1708 (Fig 7C right).
Collectively, above results indicate that a subset of EAT1-dependent 24-nt phasiRNAs, at least the versions retaining 5'-terminal cytosine, was bound by MEL1.
This study gave new insights in the relationship of tapetal bHLH proteins during early meiosis. First, the EAT1 expression is bimodal, not only in post-meiosis, but also in early meiosis https://doi.org/10.1371/journal.pgen.1007238.g007 (Fig 1). Second, the transient expression assay suggests a possibility that the transcription of EAT1 gene during early meiosis is activated by the TIP2/UDT1 heterodimer, and reinforced by the EAT1/UDT1 (Fig 6, S10 Fig). Third, both EAT1 and TIP2 can activate transcription of 24-PHAS lincRNAs and the DCL5 gene in tapetum during early meiosis (Figs 3, 5 and 6A, S10B Fig). The activation by EAT1 is thought to be independent of that by TIP2, because of no interaction between two proteins as previously reported [27,28] and shown in this study ( Fig  6A, S10 Fig). In these two pathways, UDT1 is a strong candidate for the dimerization partner of EAT1 and TIP2 (Fig 6B, S11 Fig), while dimerization of unknown bHLH proteins with EAT1 is supposed in the DCL5 transcription (Fig 6A). In the udt1 mutant, the tapetum is aberrantly vacuolated and the tetrads are degenerated during meiosis [32]. This observation is consistent to the idea that UDT1 acts with TIP2 and EAT1 in 24-nt phasiRNA biogenesis in rice anther tapetum during meiosis. The temporal replacement of binding partners from UDT1 to TDR may enable EAT1 and TIP2 to switch downstream targets from meiotic phasiRNA production to postmeiotic tapetal PCD induction.
In this study, we performed mRNA-seq and sRNA-seq to estimate 24-nt phasiRNA production only in the eat1-4 (Fig 3), but not in the tip2-2. This is because in the tip2 mutant, tapetum is replaced by undifferentiated cell layers [27,28] (S4U-S4X Fig), and the absence of 24-PHAS and DCL5 transcripts is possibly a by-product of the missing tapetum. However, the results that at least two 24-PHAS transcripts enriched at ST.2 were transcribed EAT1-independently (green spots in Figs 3A and 5C), and that non-negligible amounts of 24-PHAS and DCL5 transcripts are expressed still in eat1-4 anthers at ST.2 (Fig 3B left, Fig 5C). Taken together with the results of ChIP-qPCR and transient expression assay, it is obvious that TIP2 has an indispensable role in 24-nt phasiRNA production.
The maize (Zm) bHLH122, the EAT1 ortholog, also shows bimodal expression [25], and MALE STERILE23 (MS23), the TIP2 ortholog, promotes the expression of bHLH122/ ZmEAT1, DCL5, 24-PHAS transcripts and meiotic 24-nt phasiRNAs [12,25]. A positive interaction in the yeast two hybrid analysis (Y2H) is reported between MS32/ZmUDT1 and bHLH122/ZmEAT1, consistent to the results of this study (Fig 6, S10 and S11 Figs). Thus, the bHLH TF-mediated mechanism underlying specification and development of tapetum is well conserved in rice and maize, and commonly coupled with meiotic small RNA production. A contradiction between maize and rice is in the relationship of TIP2 and UDT1. In maize, a negative Y2H interaction of MS23/ZmTIP2 and MS32/ZmUDT1 is reported [25], whereas rice TIP2 and UDT1 interact with each other at the molecular level ( Fig 6B) and promote the activity of pPHASs, pDCL5 and pEAT1 (Fig 6A, S10B Fig). Further analyses will be necessary for conservation and differentiation of tapetal bHLH protein functions in these monocot model plants.

A possibility of intercellular mobilization of EAT1-dependent phasiRNAs in anthers
The observation that a subset of tapetum-originating phasiRNAs was sorted to MEL1 Argonaute, which is abundantly expressed in PMCs but not in tapetal cells (Fig 7). Though the possibility that 24-nt phasiRNA functions mainly in tapetum cannot be excluded, the result of this study suggests another possibility that the 24-nt phasiRNA is mobile between somatic and reproductive cells in rice anthers. This idea is attractive and proposed previously [12,13], but should be considered carefully. It is difficult to exclude the possibility that 24-nt phasiRNAs are produced cell-autonomously in PMCs by EAT1 and/or TIP2-independent pathways, for example, DNA double-strand break (DSB)-induced small RNAs [44,45]. However, we think this unlikely, because mel1 mutant anthers with few meiotic DSBs in male meiocytes [46] produce a robust level of 24-nt phasiRNAs (S1 and S4 Tables). In addition, few amounts of 24-nt phasiRNAs are detectable in eat1-4 anthers (Fig 3C, S3 Table). A recent study unveiled that 24-nt phasiRNA and miR2275 expression is depleted in two rice mutants, multiple sporo-cytes1 (msp1) and tpd1-like gene in rice1a (tdl1a), in which a subset of inner anther-wall cells turn into PMCs [24,47]. In maize, the ms23 anther lacking the tapetum fails to produce 24-nt phasiRNAs, but the ocl4 anther developing the tapetum succeeds [12]. These results suggest that 24-nt phasiRNA production occurs exclusively in tapetum, consistent to the conclusion of this study.
An alternative possibility is that precursor PHAS transcripts or their processed intermediates are transferred from tapetum, and processed into mature 24-nt phasiRNAs by 24-PBR components in PMCs. TIP2 and EAT1 are detectable in somatic companions, but hardly in PMCs (Fig 2D, S9C Fig), implicating that most of DCL5-mediated 24-PHAS processing is completed in anther tapetum. However, to answer the above question, further analyses for tissue-specific expression of precursor transcripts and 24-PBR components are required.
Another question for the intercellular small-RNA movement is whether the undetectable level of MEL1 proteins accumulates in tapetal cells during meiosis and associates with tapetum-expressing 24-nt phasiRNAs. However, MEL1 mRNA expression is ranked at the top 1.7 percentile (the 629th highest) of all protein-coding transcripts expressed in ST.2 anthers (S13 Fig), and as reflecting the higher mRNA level, the MEL1-GFP signal in male meiocytes made a striking contrast to undetectable signals in somatic anther cells in transgenic plants ( Fig 7A). Thus, small RNAs immunoprecipitated with somatic MEL1 are, if any, hard to be detected in the RIPseq analysis of anther samples, that is, the MEL1 RIPseq data of this study largely comes from the masiRNA population derived from male meiocytes. In any case, rigorous verification requires some breakthrough technologies for live-imaging of small RNAs or sequestering them into the particular cell type, such as tapetal cells.
Molecular transport in plants occur either symplastically through plasmodesmata, or apoplastically across the cell membrane, cell walls and intercellular space [48]. Tapetal cells and PMCs are connected with plasmodesmata and form symplastic continuity by the onset of meiotic leptotene (ST.2 in this study) [49,50], when EAT1-dependent meiotic 24-nt phasiRNAs are produced in tapetal cells (Fig 3). This interconnection is broken by callose accumulation [49]. Callose is the highly impermeable polysaccharide distinct from cellulose [50], and can be a barrier for apoplastic molecular movement. However, in ST.2 anthers, cellulosic components still remain between tapetum and PMCs (Fig 2B), in turn, callose accumulation is absent or less. Thus, both symplastic and apoplastic movements are currently possible mechanisms underlying meiotic phasiRNA movement in anthers.
Taking previous findings into consideration, we propose that considerable amounts of 24-nt meiotic phasiRNAs are imported from tapetum to PMCs during early meiosis in rice. If it is true, not only the phasiRNAs with the 5'-teminal cytosine (C-terminal phasiRNAs), but also non-C-terminal ones are supposed to move together in the intercellular movement, because the enrichment of C-terminal phasiRNAs in MEL1-RIPseq in this study is simply due to the selectivity of MEL1 [11]. The analysis of other Argonautes expressed in PMCs will be beneficial to trace tapetum-originating non-C-terminal phasiRNAs.
Arabidopsis AGO4 plays important roles in chromosome condensation and segregation during the first meiotic division [52], comparable to rice EAT1 function in male meiosis (Fig 2C). ZmAGO104, orthologous to Arabidopsis AGO9, is also required for meiotic chromosome condensation [53]. In either case, the relationship of Argonaute/small RNA complexes to the nuclear RdDM and histone modification will be one of the most important questions regarding epigenetic regulation of plant meiosis.
Meiosis is a special type of cell division to transmit new haplotypes to the next generation, and additionally, to survey incompatibilities in ploidy levels and chromosomal structures between both parents. This process must be strictly regulated by complicated mechanisms genetically and epigenetically. Recent genome-wide studies have revealed that small RNAmediated and non-cell-autonomous regulation is likely general in reproduction of eukaryotic species. Further analyses of tapetum-expressing bHLH TFs and meiotic phasiRNAs in anthers will bring new epigenetic insights into plant reproduction systems.

Plant materials
The eat1-4 mutant is a Tos17 insertion line produced from the rice variety, cv. Nipponbare [60], NF9876, kindly provided by the Rice Genome Resource Center, Japan. The mel1-1 mutant [7], another Tos17 line with the Nipponbare background, was kindly provided by the National Bioresource Project (NBRP) Rice, conducted by the Japan Agency for Medical Research and Development (AMED). The tip2-2 mutant is a T-DNA tag line with the genetic background of cv. Dongjin [61,62], 1B-24309, kindly provided by Dr. G. An (POSTECH, Korea). All plants were grown in moist chambers, greenhouses, and/or open paddy fields at the National Institute of Genetics (NIG), Mishima, Japan. Plant genotypes were determined by PCR using GoTaq Green Master Mix (Promega) and gene-specific and T-DNA/Tos17-internal primers (S6 Table).

Histology
Rice spikelets were fixed in PMEG buffer (50 mM PIPES, 10 mM EGTA, 5 mM MgSO 4 , and 4% glycerol, pH 6.8) containing 4% paraformaldehyde (PFA) for 3 h and washed six times in PMEG buffer for 2 hours. After dehydration using ethanol series, they were embedded in Technovit7100 resin (Heraeus Kulzer), sectioned in 2 μm thick slices using a LM2255 microtome (Leica Microsystems), stained with 0.1% toluidine blue O (Wako Pure Chemicals) and photographed using a BX50 light microscope (Olympus) and a DP50 camera system (Olympus). Cellulosic cell wall staining was conducted according to the method described previously [63]. Fluorescent images were captured using a Fluoview FV300 CLSM system (Olympus), and pseudo-colored and merged using Photoshop CS4 (Adobe Systems Inc.).
In MEL1pro-GFP-MEL1 construction, the GFP sequence was inserted just in front of MEL1 TIS in pKN16, a binary vector containing the 18 kbp MEL1 genomic fragment [7]. Two DNA fragments, corresponding to 5 0 upstream and 3 0 downstream regions of MEL1 TIS, were amplified from pKN16 with primer pairs up_nf/up_nr and up_atgf/up_r, respectively. Linkerattached sGFP coding sequence was amplified from CaMV35S-sGFP(S65T)-nos3 0 with ngfp_f/ ngfp_r primers. The PCRs were conducted using a PrimeSTAR Max DNA polymerase (TaKaRa). The three amplified DNA fragments were mixed with the NruI-AscI-digested pKN16 and incubated with an In-Fusion HD enzyme premix (TaKaRa) to assemble MEL1pro-GFP-MEL1, following manufacturer's instructions. All the primer sequences for the construction were listed in S6 Table. The constructs were transformed into rice calli using agrobacterium-mediated transformation [67], in which Hygromycin B (50 mg/L in media; Wako Pure Chemicals) or glufosinateammonium PESTANAL (5 mg/L in media; Sigma-Aldrich) was used for a positive selection.

Meiotic chromosome observation
Spikelet (lemma) and anther lengths were measured under SMZ645 stereo microscopy (Nikon). 0.8-1.2 mm anthers were fixed with 4% PFA/PMEG and provided for chromosome observations as previously described [7]. Fluorescent images of DAPI were taken as described above.

RNA extraction and quantitative RT-PCR (qRT-PCR)
Anther or spikelet samples were separated by their lengths as corresponding to ST.1-ST.6 stages (S7 Table), immediately frozen with liquid nitrogen in microtubes, and stored at -80˚C until use. Total RNAs were extracted from the samples using TRIzol reagent as manufacturer's recommendation (Life Technologies), and treated with DNase I (TaKaRa). In qRT-PCR, 1 μg of total RNA was reverse-transcribed by oligo(dT) [12][13][14][15][16][17][18] primer (Life Technologies) and Super-scriptIII reverse transcriptase (Life Technologies). The products were 20-fold diluted and supplied for real-time qPCR using gene-specific primers (S6 Table), KAPA SYBR FAST universal qPCR Kit (KAPA Biosystems) and Thermal Cycler Dice Real Time System (TaKaRa). Rice Ubiquitine gene was used as an internal standard.

mRNA-seq, sRNA-seq and data analyses
Total RNAs were extracted from ST.1, ST.2 and ST.4 anthers of wild-type and eat1-4 plants, three biological replicates each. For mRNA-seq, 1 μg of total RNA was subjected to library construction using KAPA stranded mRNA-seq Kit Illumina Platforms (KAPA biosystems). Eighteen libraries differentially indexed by FastGene Adapter kit (Nippon Genetics) were multiplexed (9 per lane) and sequenced by HiSeq2500 (Illumina) with SR50 (single ended). Adapter sequences were removed in silico using R package QuasR [68].
mRNA-seq reads were mapped on the rice genome IRGSP1.0 using Tophat (v2.0.14) [69]. Differential expression analysis of annotated genes were conducted using Cuffdiff2 program [70]. The genes fulfilling all of the following conditions were regarded as EAT1-dependent and ST.  (3) genes with each standard deviation less than a half of the FPKM mean value of three replicates in wild-type ST.2 anthers. The lincRNAs were determined by Cuffdiff2 (merged.gtf), in which protein-coding genes were removed as referring to MSU7.0 annotation, and unannotated but transcribed genomic regions larger than 200 bp were extracted. FPKM values of lincRNAs were calculated by BEDtools [71]. Furthermore, EAT1-dependent and ST.2-enriched lincRNAs were extracted according to the same conditions described above for coding genes.
For sRNA-seq, 1 μg of total RNA was provided for library construction by NEBNext Multiplex Small RNA Library Prep Set for Illumina (New England BioLabs). The libraries were 9-plexed per lane and sequenced by HiSeq2500 (illumina) with SR52, a 2-bp extended version of SR50, for higher-quality sequencing. After trimming by QuasR, 24-nt long sRNA-seq reads were extracted by ShortRead [72], and mapped to the rice IRGSP1.0 genome using Tophat, in which reads having >50 multi-hits on rice genome or any mismatches were cut off (-N 0 -g 50). If 24-nt RNAs with >10 FPKM values were mapped on each of EAT1-dependent and ST.2-enriched lincRNA loci identified above, the loci were defined as 24-PHAS loci. Regional abundance of mRNA-seq and 24-nt sRNA-seq reads mapped on the rice genome (shown in Fig 4A) was calculated in a sliding window (window; 50 kbp, step; 25 kbp) by BEDtools. Conserved motifs were searched in each 24-PHAS locus, in addition to 200 bp regions both upstream and downstream sequences, by MEME SUITE program [35]. Phased scores were calculated as described by Howell et al. [73].

Degradome-seq data analysis
A degradome-seq dataset from young panicles of indica rice variety, cv. 93-11, was obtained from Sequence Read Archive of DNA Data Bank of Japan (DDBJ-SRA) under the accession code SRR034102 [38]. Adaptor sequence and low-quality reads were removed using FASTX-toolkit (http://hannonlab.cshl.edu/fastx_toolkit/) and the reads retaining 20-or 21-nt length were mapped onto rice IRGSP1.0 genome using Bowtie 2 [74]. The frequency of 5'-end of mapped reads were manually examined within and around 24-PHAS loci identified in this study (S4 Table).

0 RACE
To determine the TSS of two 24-PHASs and a pri-miR2275 (chr5-20, chr6-97 and pri-miR2275b), the standard 5 0 rapid amplification cDNA end (5 0 RACE) method was applied using a GeneRacer kit (Thermo Fisher Scientific), total RNA from ST.2 wild-type anther, and gene specific primers (S6 Table). Eight clones from each locus were sequenced using a BigDye Terminator v3.1 cycle sequencing kit (Applied Biosystems) and a PRISM 3130xl sequencer (Applied Biosystems) and the end of the longest read(s) was marked as TSS.

Chromatin immunoprecipitation (ChIP)-qPCR
Rice young panicles from transgenic derivatives were fixed, and the anthers at early meiosis (around 0.5 mm) were supplied for ChIP as described previously [75]. The anti-GFP antibody No.598 and the normal rabbit IgG (both from MBL International) were used for positive and negative ChIP experiments, respectively. The extracted DNAs were analyzed by real-time qPCR using region-specific primers (S6 Table). The 1/10 volume of chromatin-containing samples without IP treatment was prepared for the input samples.

Transient expression assay in rice protoplast
The 2-kbp upstream sequences from the translational start site of 24-PHAS (chr5-20, chr6-97), DCL5, EAT1 and DCL3a genes, all originated from the japonica rice cv. Nipponbare, were inserted in the upstream of the firefly Luciferase CDS and the nopaline synthase (nos) terminator. This reporter construct was cloned into pBSII-SK(-) plasmid (S10A Fig). For the effector construct, the cauliflower mosaic virus 35S (CaMV35S) promoter was fused with the cDNAs of EAT1, TIP2, TDR and UDT1 genes, originated from Nipponbare ST.2 anthers, with the nos terminator. and cloned into pBSII-SK(-) (S10A Fig). For normalization of the firefly Luciferase activity, the Luciferase cDNA of Renilla reniformis were fused with the CaMV35S promoter and the nos terminator, inserted into pBSII-SK(-) (S10A Fig), and cotransfected with the effector and reporter constructs as an internal control in all experiments. All PCR primers for the above constructions were listed in S6 Table. PrimeSTAR Max DNA polymerase (TaKaRa) was used for PCR amplification according to the manufacturer's instruction. Protoplast preparation from rice seedlings, transfection of plasmids, and protein extraction from protoplasts were according to the method previously described [76]. The Luciferase activity was detected using Dual-Luciferase Reporter Assay System (Promega) and Filter MAX F5 multi-mode microplate reader (Molecular Devices).

RNA immunoprecipitation (RIP)-seq and analysis of masiRNAs
RIP fractions from wild-type, mel1-1 and eat1-4 flowers at ST.1, ST.2 and ST.4, each of which included two biological replicates, were obtained using anti-MEL1 antibody as described previously [11]. Library construction, sequencing, adapter trimming, size filtration and mapping to rice genome were done as well as sRNA-seq methods described above. Reads per million (RPM) values were calculated in the respective 24-nt RNA sequences and compared among wild-type, mel1-1 and eat1-4 fractions. In this process, 24-nt masiRNAs were defined in 24-nt RNA sequences as having !15 RPM detected in wild-type ST.1, ST.2 or ST.4 flowers, and !RPM 4-fold enriched in wild-type compared to mel1-1.