The liverwort oil body is formed by redirection of the secretory pathway

Eukaryotic cells acquired novel organelles during evolution through mechanisms that remain largely obscure. The existence of the unique oil body compartment is a synapomorphy of liverworts that represents lineage-specific acquisition of this organelle during evolution, although its origin, biogenesis, and physiological function are yet unknown. We find that two paralogous syntaxin-1 homologs in the liverwort Marchantia polymorpha are distinctly targeted to forming cell plates and the oil body, suggesting that these structures share some developmental similarity. Oil body formation is regulated by an ERF/AP2-type transcription factor and loss of the oil body increases M. polymorpha herbivory. These findings highlight a common strategy for the acquisition of organelles with distinct functions in plants, via periodical redirection of the secretory pathway depending on cellular phase transition.

E ukaryotic cells originated from prokaryotes by expanding their endomembrane network during evolution, with the last eukaryotic common ancestor (LECA) likely possessing a complex set of organelles 1 . New membrane trafficking pathways were added to the LECA endomembrane network, some of which were secondarily lost in a lineage-specific manner, resulting in divergent and organism-specific membrane trafficking systems and organelle compositions of extant eukaryotes [2][3][4] . Comparative genomics has proposed that the emergence of a novel membranebounded organelle was accompanied by the development of a novel membrane trafficking pathway, which was accomplished by expansion and functional differentiation of machinery components through gene duplication followed by neo-and/or subfunctionalization of those machineries by accumulated mutations 5 . For example, coat-protein complexes, RAB GTPases, and SNARE proteins, which act in transport vesicle formation, tethering of the vesicles to target membranes, and membrane fusion between the two membranes, respectively, were shown to expand and functionally differentiate during evolution to develop the specialized membrane trafficking system in each organism 6,7 . However, empirical support for this hypothesis, i.e. the organelle paralogy hypothesis, still remains to be provided.
In the plant lineage, several organelles and organelle functions have been uniquely acquired during evolution. For example, the plant vacuole harbours a unique function that is not shared with the animal lysosome and the yeast vacuole: storage of proteins 8 . This vacuolar function is fulfilled through the plant-unique vacuolar trafficking system, which comprises multiple vacuolar transport pathways involving plant-unique machinery components acquired during plant evolution 9 . The cell plate, which is formed during the mitotic phase to accomplish cytokinesis in land plants, is also a prominent example of plant-specific cellular structures/organelles 10 .
One of the basal-most land plant lineages, liverworts, also possesses a unique organelle, the oil body, the existence of which is a synapomorphy of this lineage. The liverwort oil body contains bioactive compounds, such as sesquiterpenoids and cyclic bisbibenzyl compounds, and is not related to the oil body that stores neutral lipids in storage organs like seeds and fruits (i.e. lipid body or oleosome). About 90% of liverwort species have this organelle; however, its origin, biogenesis, and physiological function remain unclear with controversial origins proposed from microscopic observations [11][12][13][14][15] , although the first description of the liverwort oil body dates back to 1834 16 . Through the systematic analysis of SNARE proteins in the liverwort, Marchantia polymorpha (hereafter referred to as Marchantia), we identified an oil body-resident protein, MpSYP12B 17 , which is a homolog of animal syntaxin-1 that acts in the final step of exocytosis. Functions of SYP1 members to which MpSYP12B belongs dramatically diversified during plant evolution, suggesting expanded secretory trafficking systems in plants. Arabidopsis SYP1 consists of three subgroups, SYP11, SYP12, and SYP13. The SYP13 group mediates constitutive secretion, whereas the other groups are involved in plant-specific, higher-ordered functions; for example, KNOLLE (also known as SYP111) is specifically involved in membrane fusion at forming cell plates, and PEN1 (also known as SYP121) is required for intact penetration resistance against powdery mildew fungi and regulation of the ion channel activity [18][19][20][21][22][23] . Tip growth of root hairs and pollen tubes also involves SYP1 members; distinctive SYP1 members exhibit different expression and localization patterns in tip-growing cells, indicating diverged and specialized functions of SYP1 members in polar secretion in Arabidopsis 20,[24][25][26] .
In this study, for more insight into the functional diversification of SYP1 members in Marchantia and the mechanism of oil body biogenesis, we further characterized four SYP1 members in this organism 17,27 . MpSYP12A and MpSYP12B are distinctly targeted to the forming cell plate and oil body membrane, respectively, representing the cell-phase specific redirection of the secretory pathway. Furthermore, we identified a master transcription factor of oil body formation. Our findings indicate that the cell plate and the oil body are paralogous cellular structures/ organelles acquired during plant evolution, and also provide a strong support for the organelle paralogy hypothesis.

Results
MpSYP12A is important for cell plate formation. Our phylogenetic analysis suggested that MpSYP1 members belong to two distinct subgroups: the SYP13 group (MpSYP13A and 13B) and the SYP11/12 group (MpSYP12A and 12B) that contains SYP11 and SYP12 groups in seed plants ( Supplementary Fig. 1a) 26 . MpSYP12A, 13A, and 13B were ubiquitously expressed and localized to the plasma membrane (PM) in thalli, whereas MpSYP12B exhibited specific expression in a subpopulation of thallus cells ( Supplementary Fig. 1b-e). Using ubiquitously expressed mCitrine-MpSYP1 protein fusion constructs driven by their own regulatory elements, we detected localization of MpSYP12A at forming cell plates associated with the phragmoplast that rapidly stained with FM4-64 in addition to the PM, whereas MpSYP13A, 13B, and 12B were not detected at forming cell plates ( Fig. 1a-f, Supplementary Fig. 1f, g). This localization suggested that MpSYP12A could be a functional counterpart of KNOLLE in Arabidopsis, which was further verified genetically.
We tried to generate a complete Mpsyp12a knockout mutant by genome editing (Supplementary Fig. 1h) but were not successful. Although we could not completely rule out the possibility of impaired functions of other genes due to off-target effects, it was likely because MpSYP12A is an essential gene. However, we succeeded in isolating chimeric plants comprising mutated and wild-type cells, in which we frequently observed enlarged cells with cell wall stabs, suggesting that MpSYP12A is required for cell plate formation during cytokinesis similar to the function of KNOLLE (Fig. 1g). The KNOLLE promoter is sufficient to target non-cytokinesis-specific SYP132 to forming cell plates 28 . Similarly, mCitrine-tagged MpSYP13A and 13B were localized to cell plates and the PM when expressed under the MpSYP12A promoter (Fig. 1h, i). Furthermore, a similar localization was also observed when mCitrine-MpSYP13 was expressed by the MpCYCB1 promoter ( Fig. 1j and Supplementary  Fig. 1i), indicating that the cell cycle-dependent transcriptional regulation is also critical for cell plate targeting of SYP1 proteins in Marchantia. We also found that clathrin light chain (CLC) tagged with Citrine localized to the PM and trans-Golgi network in non-dividing thallus cells, as well as on forming cell plates in Marchantia as reported in other plants ( Supplementary  Fig. 1j, k) [29][30][31] . These results strongly suggested that fundamental mechanisms of cell plate formation were conserved during land plant evolution.
The oil body and plasma membranes share common properties. Distinct from the other MpSYP1 members, MpSYP12B exhibits specific expression in oil body cells and localizes to the oil body membrane (Supplementary Fig. 2b) 17 . Oil body cells expressing a pro MpSYP12B:2×YFP construct indicated that MpSYP12B distributed around meristematic regions in young thalli, which was observed using light-sheet microscopy ( Fig. 2a; Supplementary Movie 1). These results suggest that oil body formation occurs around meristematic regions and into thalli during the growth. mCitrine-MpSYP12B expressed under its own regulatory elements localized to the oil body membrane with a faint signal on the PM (Fig. 2b), which was also confirmed by  immune-electron microscopy using an anti-GFP antibody (Fig. 2c, d). We then tested whether organelle markers for the endoplasmic reticulum (MpSEC20, MpUSE1A, MpSEC22, mRFP-HDEL, and GFP-HDEL), Golgi apparatus (MpGOS11 and MpSFT1), trans-Golgi network (MpSYP6A and MpSYP4), and tonoplast (MpSYP2 and MpVAMP71) were targeted to the oil body membrane, none of which were detected ( Fig. 2e and Supplementary Fig. 2c, d). Using transmission electron microscopy, we found that clathrin-coated vesicles formed from the oil body membrane, and emergence and disappearance of clathrinpositive foci at the oil body membrane was also observed in transgenic plants expressing Citrine-tagged MpCLC1 ( Fig. 2g; Supplementary Movie 2).
Given that clathrin-mediated endocytosis occurs at the PM, and SYP1 members generally function on the PM, the oil body membrane could share similar characteristics with the PM. This hypothesis was supported by dual localization of the PM proteins, PM-type aquaporin MpPIP2 (Fig. 2h) and PM-resident SNARE MpSYP13A (Fig. 3a) driven by their own regulatory elements. The plasma membrane-like nature of the oil body membrane and the existence of the SYP1 member that is homologous to KNOLLE suggested that the oil body could be formed by the fusion of secretory vesicles similar to the cell plate. Thus, the luminal space of the oil body should be topologically equivalent to the extracellular space, which was tested by expressing a general secretion marker sec-mRFP, which is composed of the signal peptide for ER translocation and mRFP under the constitutive MpEF1α promoter. In non-oil body cells, mRFP fluorescence was detected only in the extracellular space ( Supplementary Fig. 2e). However, mRFP accumulated in the oil body in addition to the extracellular space in oil body cells, demonstrating equivalent topology between the lumen of the oil body and extracellular space (Fig. 2i). Unlike MpSYP12A and KNOLLE, MpSYP12B loss of function did not result in a detectable abnormality in oil body formation and transport of MpSYP13A or sec-mRFP (Supplementary Fig. 3), probably reflecting a functional redundancy between MpSYP12B and 13A, similar to the partial functional redundancy of KNOLLE and SYP132 32 .
The oil body is formed by redirecting the secretory pathway. Distinct from MpSYP13A with the dual localization to the plasma and oil body membranes, the close homolog MpSYP13B was only localized to the PM in oil body cells (Fig. 3a, b). Surprisingly, both of these proteins were targeted predominantly to the oil body membrane when expressed under the regulation of the MpSYP12B promoter (Fig. 3c, d). This effect was not restricted to SYP1 homologs; the PM-resident protein MpPIP2 and a general secretion cargo sec-mRFP were also targeted almost exclusively to the oil body when driven by the MpSYP12B promoter (Fig. 3e, f). These results indicated that the secretory pathway is redirected inward and vesicles fuse with each other to form the oil body during the phase in which the MpSYP12B promoter is active. In contrast, when the MpSYP13B promoter is active, the secretory pathway is directed to the PM and extracellular space, which was further confirmed by the accumulation of sec-mRFP predominantly in the extracellular space when the construct was driven by the MpSYP13B promoter (Fig. 3g). Other organelle markers did not change their localization even when expressed under the MpSYP12B promoter ( Supplementary Fig. 2c, d). These results indicate that the redirections of the secretory pathway in oil body cells is regulated at the transcription level. The reorientation of secretory directions should take place periodically and repeatedly through oil body cell development, because we never observed oil body membrane localization for MpSYP13B at any developmental stages of oil body cells (Fig. 4a-c), and increase in the size of the oil body cell along with the increase in the oil body size was observed (Fig. 4d). Based on these findings, we propose the "oil body cycle hypothesis", which states that Marchantia oil body cells cycle between two distinct cellular phases: the "PM phase" in which the secretory pathway is oriented to the PM and extracellular space, and the "oil body phase" when the secretory pathway is oriented to the oil body, with the phase transition under the regulation of a transcriptional regulatory system (Fig. 4e, f).
MpERF13 transcription factor regulates oil body formation. To investigate the regulatory system of oil body formation and the oil body cycle, we screened mutants defective in oil body formation ( Supplementary Fig. 4a). From 48,825 T-DNA insertion lines, we identified a mutant (Mperf13 GOF ) with an increased number of oil bodies. The mutant gemmae contained 334.4 ± 85.6 oil bodies (mean ± s.d.), whereas wild-type gemmae possessed 51.9 ± 8.7 oil bodies under our experimental conditions ( Fig. 5a-e). The T-DNA was inserted 2674-bp upstream of the start codon of Mapoly0060s0052.1 (MpERF13), which encodes a putative ERF/ AP2 transcription factor in the subgroup containing DREB1A and TINY in Arabidopsis ( Supplementary Fig. 4b). Moreover, MpERF13 and MpSYP12B transcripts accumulated in this mutant, compared to the wild type, as detected by RNAsequencing (RNA-Seq) and quantitative RT-PCR analyses (Supplementary Fig. 5b). We then generated knockout mutants in which the MpERF13 gene was deleted by genome editing (Supplementary Fig. 4c). Two independent mutants (Mperf13-1 ge and Mperf13-2 ge ) exhibited no detectable abnormalities in thallus development and reproductive growth; however, these mutants completely lacked oil bodies in the gemma and thalli (Fig. 5c-e and Supplementary Fig. 4f, g). These phenotypes of gain-and loss-of-function mutations suggested that MpERF13 is a major transcription factor regulating oil body formation. Consistently, the MpERF13 promoter was highly active in oil body cells (Fig. 5f).
As MpERF13 is homologous to ERF/AP2 transcription factors, this protein was expected to be involved in the transcriptional regulation of oil body formation and/or the oil body cycle. To identify genes downstream of MpERF13, we analysed transcriptomes in Tak-1 (wild type) and Mperf13 GOF and Mperf13-1 ge mutant lines by RNA-Seq. Through the comparison between the three groups using an ANOVA-like test, we identified 136 differentially expressed genes (DEGs) other than MpERF13 whose expression was higher than two log 2 -fold change (FC) in Mperf13 GOF compared with Tak-1, and lower than -2 log 2 FC in Mperf13-1 ge than Tak-1 (FDR < 0.01) (Supplementary Data 1). To verify the RNA-Seq result, we also performed the quantitative RT-PCR analysis for selected 11 DEGs including MpERF13 and MpSYP12B, which confirmed that all of these DEGs were highly expressed in Mperf13 GOF , but their expression was significantly lower or not detectable in Mperf13-1 ge and Mperf13-2 ge compared with Tak-1 ( Supplementary Fig. 5). These data strongly suggested that these 136 genes function downstream of MpERF13 in oil body cells, which was further verified for the two genes Mapoly0083s0014.1 and MpMYB02 (Mapoly0006s0226). Mapoly0083s0014.1 was positively regulated by MpERF13 and encodes a member of the G subfamily of ABC transporters 27 , which we named MpABCG1. Although ABCG members are generally targeted to the PM and mediate efflux of secondary metabolites from intracellular to extracellular spaces 33,34 , mCitrine-MpABCG1 was targeted to the oil body membrane when expressed under its own regulatory elements (Fig. 5g); however, this protein was targeted to the PM in non-oil body cells when expressed under the MpSYP13B promoter (Fig. 5h). MpMYB02, whose expression was also positively regulated by MpERF13, is reported to be a transcription factor responsible for the production of marchantin A and its derivatives that accumulate in oil bodies 15,35 . MpABCG1 and MpMYB02 are specifically expressed in oil body cells and when sec-mRFP was driven by the promoters of these genes, it was specifically targeted to the lumen of oil bodies. These results indicate that MpABCG1 and MpMYB02, presumably together with other DEGs we identified, act during oil body formation downstream of MpERF13.
Lastly, we asked why liverworts possess oil bodies. Using fresh and alcohol-washed liverworts, the oil body was proposed to be an effective chemical protection from herbivores 36 , which was also recently supported by a feeding assay using Armadillidium vulgare (pill bug) 37 . We fed starved pill bugs with wild-type liverwort and genetically established mutants that contain extra or no oil bodies. Thalli of Tak-1 and Mperf13 GOF remained almost intact after 24 h of herbivory, while the areas of Mperf13-1 ge and Mperf13-2 ge thalli were significantly reduced after herbivory ( Fig. 6 and Supplementary Fig. 6). This result further supports that the oil body is effective in protecting liverworts from herbivores.

Discussion
Plant SYP1 family generally localizes to the PM and execute membrane fusion between the PM and secretory vesicles. The SYP1 family, which has expanded during land plant evolution, consists of two subgroups, the SYP13 and SYP11/12 groups. The SYP11/12 group in seed plants is further divided into two subclades with distinct functions, SYP11 and SYP12 groups, suggesting sub-and/or neofunctionalization in the SYP11/12 group in an ancestor of seed plants ( Supplementary Fig. 1a) 17,26 .
Another key machinery component of membrane trafficking, RAB GTPase, is also shown to have been expanded during land plant evolution, in which secretory RAB GTPases including RABA/RAB11 exhibit remarkable functional diversification 38,39 . Cell plate formation during cytokinesis involves distinctive subgroups of the SYP1 and RAB11 families 28,40,41 , probably reflecting that diversified SYP1 and RAB11 members were coopted to the trafficking pathway to the cell plate. In this study, we found that an organelle specific to liverworts, the oil body, should also be acquired through paralogous expansion followed by neofunctionalization of the SYP1 group. Although the functions of the cell plate and oil body are totally different, they share several important traits that include specific SYP11/12 members that reside at their membranes, clathrin-coated vesicles generated from the plasma and oil body membranes, and both the cell plate and oil body are formed by reorientation of secretory pathways. In both cases, the orientation of the secretory pathway is redirected in response to transitioning cellular states during the cell cycle for cell plate formation and the hypothetical oil body cycle for oil body formation. The orientation of the plant secretory pathway is also modulated during plant-microbe interactions; symbiotic arbuscular mycorrhizal fungi redirects the plant secretory pathway to form the arbuscule in host plant cells that acts as an interface for sugar and mineral exchange 42,43 . Thus, the redirection of the secretory pathway is a common strategy for plants to acquire novel organelles/cellular structures that have resulted in maximizing fitness for land plants during evolution. A lineage-specific acquisition of secretion-related organelles associated with duplication and neofunctionalization of trafficking machineries, such as RAB GTPase, was also reported in protozoa 6,44 . Thus, acquisition of new organelles by redirection of the secretory pathway through the paralogous expansion of trafficking machinery components followed by their neofunctionalization would be a commonly adopted evolutionary path, concordant with the organelle paralogy hypothesis.
The oil bodies in liverworts contain a variety of secondary metabolites such as sesquiterpenoids and bisbibenzyls 11,15 , whose biological functions remain to be determined. In this study, we demonstrated that the oil body reduces Marchantia herbivory using mutant plants with excess or without oil bodies. Consistently, a mutant of the MpC1HDZ transcription factor with a reduced number of the oil body was also recently shown to exhibit increased herbivory 37 . The secondary metabolites stored in the liverwort oil body were also reported to exhibit a diverse range of bioactivities, e.g., anti-cancer and anti-virus activities 11 . Thus, modulation of the oil body formation using MpERF13 would provide a promising methodology for the pharmaceutical application of these compounds. We also demonstrated in this study that a general secretion marker sec-mRFP was targeted to the luminal space of the oil body when expressed during the hypothetical oil body phase, suggesting that any soluble proteins could be targeted and stored in the oil body just by adding the signal peptide for translocation to the ER lumen. The knowledge obtained from the study of the Marchantia oil body, in combination with attempts to modulate organelle biogenesis using other tools 45 , would open a new road to engineering organelle functions to maximize and/or design plant cellular functions.  17 and Bowman et al. 27 , were used for phylogenetic analyses of SYP1 and ERF/AP2 members, respectively. Multiple sequence alignments were performed using the MUSCLE program version 3.8.31 46,47 with the default parameter. The alignment gaps were removed by Gblocks program version 0.91b 48,49 or manually. The maximum likelihood phylogenetic analyses were performed using PhyML 3.0 50 under the LG model. Bootstrap analyses were performed by resampling 1000 sets.
Vector construction. Open reading frames (ORFs) and genomic sequences of Marchantia genes were amplified by polymerase chain reaction (PCR) from cDNA and genomic DNA prepared from the accession Tak-1; the amplified products were subcloned into pENTR/D-TOPO (ThermoFisher, Waltham, USA) according to the manufacturer's instructions. The primer sequences and sizes of amplified products are listed in Supplementary Data 2. For the construction of pENTR pro MpERF13, the 5′ sequence (promoter + 5′ UTR) was amplified and subcloned into pENTR/D-TOPO. The resultant sequence was introduced into pMpGWB316 51 using the Gateway LR Clonase™ II Enzyme Mix (ThermoFisher) according to the manufacturer's instructions. For the construction of pENTR genomic mGFP-MpSYP12A and mRFP-MpSYP13B constructs, the cDNA for mGFP and mRFP were inserted into SmaI and BamHI sites in the pENTR genomic XFP(SmaI)-MpSYP12A and pENTR genomic XFP(BamHI)-MpSYP13B, which were previously prepared 17 using the In-Fusion HD Cloning System (Clontech, Shiga, Japan). For the construction of pENTR genomic mCitrine-MpPIP2, pENTR genomic mCitrine-MpSYP4, pENTR genomic mCitrine-MpSEC22, and pENTR genomic mCitrine-MpABCG1, genomic sequences of the protein-coding regions and 3′ flanking sequences were amplified with a SmaI site at the 5′ end and subcloned into the pENTR vector. The 5′ sequences (promoter + 5′ UTR) and cDNA of mCitrine were amplified and inserted into NotI and SmaI sites, respectively, of the resultant pENTR vectors using the In-Fusion HD Cloning System. The resultant chimeric genes were then introduced into pMpGWB301 51 using the Gateway LR Clonase™ II Enzyme Mix. For the construction of pENTR pro MpSYP12B:mCitrine, the pro Mp-SYP12B:mCitrine sequence was amplified from the pENTR genomic mCitrine-MpSYP12B vector, and subcloned into pENTR/D-TOPO. The resultant sequence was introduced into pMpGWB307 51 using the Gateway LR Clonase™ II Enzyme Mix to create the pMpGWB307 pro MpSYP12B:2×YFP construct. For the construction of pENTR genomic MpCLC1-mCitrine, the genomic sequence comprising the 5′ flanking (promoter + 5′ UTR) and protein-coding sequences were amplified with a SmaI site at the 3′ end and subcloned into the pENTR vector. The 3′ flanking sequence including 3′ UTR and cDNA for mCitrine were amplified and inserted into the AscI and SmaI sites, respectively, of the resultant pENTR vectors using the In-Fusion HD Cloning System. The resultant chimeric genes were then introduced into pMpGWB301 51 using the Gateway LR Clonase™ II Enzyme Mix. For the construction of pENTR sec-mRFP, sequences for the signal peptide and mRFP were independently amplified by PCR, and the sequence of sec-mRFP was then amplified by PCR using the mixture of the amplified products as templates.
The amplified products were subcloned into pENTR/D-TOPO. For the construction of pENTR mRFP-HDEL, the sequence of SP-mRFP-HDEL was amplified by PCR using pENTR sec-mRFP as template, and subcloned into pENTR/D-TOPO. For the construction of pENTR pro MpSYP12B:sec-mRFP, pENTR pro MpSYP13B:sec-mRFP, pENTR pro MpABCG1:sec-mRFP, and pENTR pro MpMYB02:sec-mRFP, the 5′ sequences (promoter + 5′ UTR) were amplified and were inserted into NotI sites of the pENTR sec-mRFP vector using the In-Fusion HD Cloning System. The resultant sequences were then introduced into pMpGWB301 51 using the Gateway LR Clonase™ II Enzyme Mix. To construct pMpGWB301-derived Gateway vectors ( Supplementary Fig. 2a), the promoter:mCitrine sequences were amplified from the genomic constructs described above, and inserted at a HindIII site of pMpGWB301 using the In-Fusion HD Cloning System. pENTR MpTUB2 52,53 was kindly provided from Dr. R. Nishihama and introduced into pMpGWB301 pro MpSYP2:mCitrine-GW ( Supplementary Fig. 2a) using the Gateway LR Clonase™ II Enzyme Mix. For construction of the genome-editing vectors, the target sequences were selected using CRISPR direct (https://crispr.dbcls.jp/) 54 ; double-stranded oligonucleotides of the target sequences were inserted into the pMpGE_En03 vector 55 . The resultant gRNA cassettes were introduced into the pMpGE010 or pMpGE011 vectors 55 using the Gateway LR Clonase II Enzyme Mix. For the construction of the homologous recombination-mediated gene targeting vector, the 3.5-kb homologous genomic sequences were amplified from the Tak-1 genome and inserted at PacI and AscI sites of the pJHY-TMp1 vector 56 using the In-Fusion HD Cloning System.
Plant materials and transformation. Marchantia accession Takaragaike-1 (Tak-1, male) and Takaragaike-2 (Tak-2, female) 57 were used throughout the study. The growth conditions and transformation methods were previously described 17,57,58 . Thalli were grown asexually and maintained on 1/2× Gamborg's B5 medium containing 1.0% (w/v) agar at 22°C under continuous white light. The transition from the vegetative to reproductive growth was induced by supplementing far-red light, and spores were generated by crossing male and female lines. For transformation using regenerating thalli, 11-day-old thalli, from which apical regions including meristems were removed, were cultured on 1/2× Gamborg's B5 medium containing 1.0% (w/v) sucrose and 1.0% (w/v) agar at 22°C for three days under continuous white light for regeneration. The regenerating plantlets were co-cultured with the agrobacterium strain GV2260 harbouring a binary vector in the liquid 0M51C medium containing 2% (w/v) sucrose and 100 μM acetosyringone at 22°C under continuous white light with agitation at 130 rpm for three days. Transformants were selected on 1/2× Gamborg's B5 medium plates containing 10 mg L −1 hygromycin B and 100 mg L −1 cefotaxime for pMpGWB101, pMpGWB103, pMpGE010, and pJHY-TMp1 vectors, and on the medium plates plus 0.5 μM chlorsulfuron and 100 mg L −1 cefotaxime for pMpGWB301, pMpGWB307, pMpGWB316, and pMpGE011, respectively. The method of transformation using sporelings is described in the forward genetic screening for the oil body formation mutants section. M. polymorpha expressing mCitrine-MpSYP12A, mCitrine-MpSYP12B, mCitrine-MpSYP13A, mCitrine-MpSYP13B, mCitrine-MpSYP2 or mCitrine-MpVAMP71 under the regulation of their own regulatory elements (e.g. promoter, UTRs, and introns), and M. polymorpha expressing mRFP-MpSYP6A or GFP-HDEL driven by the MpEF1α promoter were previously reported 17,51,59 .  Tak-1 After Before a Fig. 6 Loss of the oil body increases Marchantia herbivory. a Ten-day-old thalli in indicated genotypes before they were fed to starved pill bugs (Before) or 24 h after feeding (After). b Ratios of thallus areas of After: Before pill bug feeding were calculated to indicate thallus area change ratios (n = 36 thalli for each genotype). Bars indicate means ± s.d. Statistical analyses between Tak-1 and each genotype were conducted by two-tailed Welch's t-test. The pvalues were 1.07 × 10 −4 for Mperf13 GOF and 4.10 × 10 −17 for Mperf13-1 ge . Bars = 1 cm.
Genotyping. Genomic DNA was extracted from the thalli in an extraction buffer containing 1 M KCl, 100 mM Tris-HCl (pH 9.5), and 10 mM EDTA. This DNA was used as templates for PCR. Genome fragments with putative mutation sites were amplified by PCR using KOD FX Neo (TOYOBO, Osaka, Japan) according to the manufacturer's instructions. The primers used in PCR-based genotyping are listed in Supplementary Data 2.
Electron microscopy. For immunoelectron microscopy, five-day-old M. polymorpha thalli from Tak-1 and thalli expressing mCitrine-MpSYP12B were used. The sandwiched samples with copper disks were frozen in liquid propane at −175°C. After the samples were frozen, they were freeze-substituted with 0.5% tannic acid in acetone and 3% distilled water at −80°C for 48 h. The samples were then kept at −20°C for 3 h followed by incubation at 4°C for 1 h. Next, the samples were dehydrated in anhydrous ethanol 3 times for 1 h each, followed by infiltration with a 50:50 mixture of ethanol and LR white resin at 4°C for 3 h. The samples were transferred to new 100% LR white resin and incubated at 4°C for 1 h, and this process was repeated three times. Next, the samples were transferred to fresh 100% resin and polymerized at 50°C overnight. The blocks were sliced into ultra-thin (80-nm) sections using an ultramicrotome, and the sections were placed on nickel grids. The grids were then incubated with the primary antibody (Anti-GFP polyclonal antibody, ab6556, abcam) in PBS plus 1% BSA for 90 min at room temperature, followed by three rinses with PBS plus 1% BSA for 1 min each. The grids were then floated on drops of secondary antibody conjugated with gold particles (10 nm) (anti-rabbit IgG polyclonal antibody, EM.GAR10, BBI Solutions) for 1 h at room temperature. The primary and secondary antibodies were used at 1:50 and 1:20, respectively. After rinsing with PBS, the grids were placed in 2% glutaraldehyde in 0.1 M phosphate buffer, dried, stained with 2% uranyl acetate for 15 min, rinsed with distilled water, and then subjected to secondary staining with lead stain solution at room temperature for 3 min. The grids were observed using a transmission electron microscope (JEM-1200EX; JEOL Ltd.) at an acceleration voltage of 80 kV. Digital images were obtained with a CCD camera (VELETA, Olympus Soft Imaging Solutions GmbH).
For morphological observation of the oil body, five-day-old Tak-1 thalli were fixed with 2% paraformaldehyde and 2% glutaraldehyde in 0.05 M cacodylate buffer pH 7.4 at 4°C overnight. The fixed samples were washed three times with 0.05 M cacodylate buffer for 30 min each and were then post-fixed with 2% osmium tetroxide in 0.05 M cacodylate buffer at 4°C for 3 h. The samples were dehydrated in graded ethanol solutions (50% and 70% ethanol for 30 min each at 4°C , 90% for 30 min at room temperature, four times with 100% for 30 min each at room temperature, and 100% overnight at room temperature). The samples were infiltrated with propylene oxide (PO) two times for 30 min each, and then placed into a 70:30 mixture of PO and resin (Quetol-651; Nisshin EM Co.) for 1 h. The caps of tubes were opened overnight to volatilize PO. The samples were transferred to fresh 100% resin and polymerized at 60°C for 48 h. Ultra-thin sample sections were mounted on copper grids, stained with 2% uranyl acetate and lead stain solution (Sigma-Aldrich), and observed under a transmission electron microscope (JEM-1400Plus; JEOL Ltd) at an acceleration voltage of 80 kV. Digital images (2048 × 2048 pixels) were obtained using a CCD camera (VELETA; Olympus Soft Imaging Solutions GmbH).
Light-sheet microscopy. Five-day-old thalli expressing 2×Citrine driven by the MpSYP12B promoter were embedded in low-melt agarose gel and observed using a Light-sheet Z.1 microscope (Carl Zeiss) equipped with a water immersion lens (×5, numerical aperture = 0.16). The samples were excited at 488 nm (Argon 488 laser). Acquisition and construction of three-dimension images from multi-angle images were conducted using the ZEN2014 software (Carl Zeiss). The images were processed digitally with Imaris 8.2 (Bitplane) and Photoshop software.
Forward genetic screening for the oil body formation mutants. The workflow of the mutant screening is described in Supplementary Fig. 4a. Liverworts with transfer DNA (T-DNA) insertions were generated by co-culture of M. polymorpha sporelings with the agrobacteria strain GV2260 harbouring the binary vector pCAMBIA1300 (https://cambia.org/welcome-to-cambialabs/cambialabs-projects/ cambialabs-projects-legacy-pcambia-vectors-pcambia-legacy-vectors-1/ cambialabs-projects-legacy-pcambia-vectors-list-of-legacy-pcambia-vectors-3/). T 1 plants were first selected on 1/2× Gamborg's B5 plates containing 10 mg L −1 hygromycin and 100 mg L −1 cefotaxime. Hygromycin-resistant T 1 plants were then transferred and incubated on hygromycin-free medium for about one week, followed by BODIPY 493/503 staining to visualize oil bodies. Visible screening for mutants defective in oil body formation was performed using a fluorescent stereoscopic microscope (M165 FC, Leica). To identify flanking sequences of T-DNA insertions, thermal asymmetric interlaced-PCR (TAIL-PCR) was performed as described previously 57,60,61 with minor modifications using crude-extracted DNA as templates. Flanking sequences were amplified using KOD FX neo DNA polymerase and T-DNA-specific (TR1-3 and TL1-3 for right and left borders of T-DNA, respectively) and universal adaptor (AD1-6) primers. The reaction cycles are shown in Supplementary Data 3. After agarose gel electrophoresis of the final TAIL-PCR products, DNA bands were excised and purified using the Wizard SV Gel and PCR Clean-Up System (Promega). Purified products were directly sequenced using TR3 or TL3 primers. The T-DNA insertion sites were identified using the genome sequence registered in MarpolBase, genome version 3.1.
RNA extraction and RNA-sequencing analysis. Total RNA from Tak-1, Mperf13 GOF and Mperf13-1 ge thalli was extracted by RNeasy (QIAGEN) according to the manufacturer's instructions. The eluted total RNA samples were treated with DNase I (Takara) to remove DNA contamination. The quality and quantity of the total RNA were evaluated with a NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific), Qubit 2.0 Fluorometer (Thermo Fisher Scientific) and a Bioanalyzer RNA6000 Nano Kit (Agilent Technologies). The sequence libraries were prepared with the TruSeq RNA Sample Prep Kit v2 (Illumina) according to the manufacturer's low sample protocol. The quality and quantity of each library were determined using a Bioanalyzer with the High Sensitivity DNA kit (Agilent Technologies) and the KAPA Library Quantification Kit for Illumina (Illumina). Equal amounts of each library were mixed to make the 2 nM pooled library. Illumina sequencing was performed using a HiSeq 1500 platform (Illumina) to produce 126-bp single-end reads. Three biological replicates were prepared for the library construction and RNA-Seq analysis. All reads are available through the Sequence Read Archive (SRA) under the accession number DRA009193. The abundance of transcripts from Illumina RNA-Seq data were quantified with the kallisto program (version 0.43.1) 62 with the default parameters using the primary transcript dataset from the Marchantia genome version 3.1. To estimate DEGs, an ANOVA-like test was performed by edgeR (version 3.18.1) 63 . Genes with FDR values lower than 0.01 and an absolute log 2 FC > 2 were considered to be differentially expressed. Among all combinations of DEGs (log 2 FC(Mperf13 GOF /Tak-1) > 2 and log 2 FC(Mperf13-1 ge /Tak-1) > 2, log 2 FC(Mperf13 GOF /Tak-1) > 2 and log 2 FC(Mperf13-1 ge /Tak-1) < -2, log 2 FC(Mperf13 GOF /Tak-1) < −2 and log 2 FC (Mperf13-1 ge /Tak-1) > 2, log 2 FC(Mperf13 GOF /Tak-1) < −2 and log 2 FC(Mperf13-1 ge /Tak-1) < −2), we selected the set log 2 FC(Mperf13 GOF /Tak-1) > 2 and log 2 FC (Mperf13-1 ge /Tak-1) < −2 as candidate genes regulated downstream of MpERF13, which included both MpERF13 and MpSYP12B.
Quantitative reverse transcription (qRT)-PCR. For qRT-PCR, cDNA was synthesized from total RNA using SuperScript III Reverse Transcriptase (Invitrogen) and an oligo dT(18) primer according to the manufacturer's instructions. qRT-PCR was performed with LightCycler 480 (Roche) using FastStart SYBR Green Master (Roche) according to the manufacturer's protocol. Sequences of primers used are listed in Supplementary Data 2. MpAPT (Mapoly0100s0027.1/Mp3g25140.1) was used as a housekeeping reference for normalization 64 . Three biological replicates were prepared, and three technical replicates were performed for each reaction.
Pill bug feeding assay. The feeding assay was performed according to Nakazaki et al. 65 with modification. Gemmae were cultured on 1/2× Gamborg's B5 medium containing 1% (w/v) agar and 1% (w/v) sucrose for five days at 22°C under continuous white light (50 μmol m −2 s −1 ). Five-day-old thalli were transferred onto 1/2× Gamborg's B5 medium containing 1% (w/v) agar without sucrose and cultured for additional five days under the same conditions. Pill bugs (Armadillidium vulgare) were collected in the Myodaiji area of NIBB (Nishigonaka 38, Myodaiji, Okazaki 444-8585 Aichi, Japan). Pill bugs were maintained on Prowipe (Daio Paper, Tokyo, Japan) moistened with sterilized water for 48 h at 22°C in the dark without feeding before the assay. Six pill bugs were introduced into each medium plate containing 10-day-old thalli and kept for 24 h in the dark at 22°C. The thallus areas were calculated using ImageJ (National Institute of Health, https://imagej.nih.gov/ij/). Gene accession numbers. Accession numbers of genes used in this study are listed in Supplementary Data 4. We followed the nomenclature of genes, proteins, and mutants of Marchantia reported in Bowman et al. 66 . Gene IDs were taken from MarpolBase (http://marchantia.info/), genome version 3.1 and version 5.1 27,67 .
Statistics and reproducibility. For confocal and electron microscopic images, we presented representative images among biological and technical replicates. The numbers of generated transgenic lines and observed cells are listed in Supplementary Data 5. Statistics were calculated with Excel 2016, and graphs were drawn using Excel 2016 or GraphPad Prism 8.4.3.