CER16 Inhibits Post-Transcriptional Gene Silencing of CER3 to Regulate Alkane Biosynthesis1[OPEN]

Alkane biosynthesis in cuticular wax of Arabidopsis associated with stress responses is regulated through an RNA quality control mechanism. The aerial surfaces of land plants have a protective layer of cuticular wax. Alkanes are common components of these waxes, and their abundance is affected by a range of stresses. The CER16 protein has been implicated in alkane biosynthesis in the cuticular wax of Arabidopsis (Arabidopsis thaliana). Here, we identified two new mutant alleles of CER16 in Arabidopsis resulting in production of less wax with dramatically fewer alkanes than the wild type. Map-based cloning with genetic analysis revealed that the cer16 phenotype was caused by complete loss of AT5G44150, encoding a protein with no known domains or motifs. Comparative transcriptomic analysis revealed that transcripts of CER3, previously shown to play a principal role in alkane production, were markedly reduced in the cer16 mutants. To define the relationship between CER3 and CER16, we transformed the full CER3 gene into a cer16 mutant. Transgenic CER3 expression was silenced, and levels of small interfering RNAs targeting CER3 were significantly increased. Mutating two major components of the RNA-silencing machinery in a cer16 genetic background restored CER3 transcript levels to wild-type levels, with the stems restored to wild-type glaucousness. We suggest that CER16 deficiency induces post-transcriptional gene silencing of both endogenous and exogenous expression of CER3.

The surfaces of aerial organs of vascular plants are covered with a layer of cuticular wax that protects against environmental stresses (Bernard and Joubès, 2013;Lee and Suh, 2015). In Arabidopsis (Arabidopsis thaliana), cuticular wax is a mixture of very-long-chain (VLC) fatty acids and derived compounds synthesized by epidermal cells (Jetter and Kunst, 2008). These wax components include alkanes, secondary alcohols, and ketones, which are produced via the alkane-forming pathway, and primary alcohols and esters, which are generated from the so-called alcohol-forming pathway (Bernard and Joubès, 2013). Alkanes represent more than 50% of total wax on Arabidopsis stems and affect plant responses to biotic and abiotic stresses (Bourdenx et al., 2011).
The formation of alkanes is catalyzed by a complex of interacting proteins. This complex contains the VLC aldehyde decarbonylases CER1, CER1-LIKE1 and CER3, and cytochrome b5 (CYTB5; Chen et al., 2003;Bernard et al., 2012;Pascal et al., 2019). The CER1 and CER3 proteins are the principal contributors of VLC alkaneforming enzymatic activity. CER1 was previously reported as being regulated at transcriptional levels. For example, CER1 expression is negatively regulated by the AP2/ERF transcription factor DEWAX, which negatively regulates wax biosynthesis during daily dark and light cycles (Go et al., 2014). Additionally, a recent study in our laboratory revealed that CER1 is also regulated at the posttranscriptional level by the miR156-SPL9 (SQUAMOSA PROMOTER BINDING PROTEIN-LIKE 9) module in Arabidopsis (Li et al., 2019). Other transcription factors regulating CER3 have been identified. For example, MYB96, an abscisic acid-responsive R2R3-type MYB transcription factor, regulates CER3 expression during drought (Seo et al., 2011). CER3 is also controlled by several epigenetic modulating factors; for example, it was one target of the Arabidopsis histone methyl transferases SET DOMAIN GROUP8 (SDG8) and SDG25, which regulate wax biosynthesis through histone Lys methylation or H2B ubiquitination (Lee et al., 2016). Also, GENERAL CONTROL NON-REPRESSED PROTEIN5 (GCN5) encodes a histone acetyltransferase that modulates CER3 expression via H3K9/14 acetylation (Wang et al., 2018). Meanwhile, CER3 was also regulated at posttranslational level by SMALL AND GLOSSY LEAVES1 (SAGL1) mediated proteasome-dependent degradation .
In addition to these regulators, other studies also reported that the expression of CER3 is regulated posttranscriptionally, by RNA quality control (Hooker et al., 2007;Lam et al., 2012;Lam et al., 2015;Zhao and Kunst, 2016). RNA quality control is a surveillance mechanism in eukaryotes that eliminates selected endogenous aberrant RNAs, thereby preventing activation of the posttranscriptional gene silencing (PTGS) pathway (Liu and Chen, 2016). In plant, this surveillance mechanism involves two major exoribonucleolytic pathways that degrade RNAs in both directions (from 59 to 39 or from 39 to 59; Liu and Chen, 2016). There are several reports that expression of CER3 is regulated by the 39 to 59 degradation system (Lam et al., 2015;Zhao and Kunst, 2016). In addition, the core subunit of the RNA exosome, CER7/RRP45B, acts in the biosynthesis of wax on stems and siliques by regulating CER3 expression (Hooker et al., 2007). The silencing of CER3 in the cer7 mutant was attributed to the activation of PTGS pathway, as suppression of components involved in the PTGS pathway rescued CER3 expression and restored the wild-type phenotype in the cer7 background (Lam et al., 2012;Lam et al., 2015). The SUPERKILLER (SKI) family proteins, which are known to assist the exosome in the cytosol, are also found to participate in the CER7/RRP45B-mediated control of wax deposition (Zhao and Kunst, 2016). Although much progress was made in elucidating the regulation of wax biosynthesis mediated by RNA quality control, it is still unknown what other factors participate this process.
The Arabidopsis mutant cer16-1 in the Landsberg erecta (Ler-0) background is a wax-deficient mutant that was generated by fast neutron treatment (Koornneef et al., 1989). The cer16-1 mutant showed significant reduction on wax load in both stems and leaves, especially alkane production (Jenks et al., 1995). However, the gene is responsible for this phenotype has not been established. Using positional cloning, we located the mutation to chromosome 5. We also noted that six genes have been deleted at that location. To clarify which gene mutation caused the phenotype, we conducted specific genetic studies to identify the gene responsible for the cer16 phenotype as AT5G44150. Transcriptomic analysis reveals that CER3 transcripts were drastically reduced in cer16, and we also concluded that the CER3 transgene was silenced in cer16. Since the phenotypes associated with cer16 and cer7 are similar, we hypothesized that CER16 deficiency also triggered the PTGS pathway as it does in cer7. This led us to assess the levels of small interfering RNAs (siRNAs) targeting CER3. We then mutated components of the PTGS pathway in cer16, and found that CER16 deficiency induces the PTGS pathway. According to a recent bioRxiv preprint (Lange et al., 2019), CER16 is an exosome cofactor. We suggest here that CER16 regulates alkane biosynthesis through an RNA quality control mechanism.

Molecular Identification of CER16
The CER16 locus was previously rough-mapped to a 38.2 cM region on chromosome 2 (Rashotte et al., 2004). However, we found that this locus was linked not with the markers on chromosome 2, but with the markers on chromosome 5. We determined that CER16 is located in a region between two simple sequence length polymorphism (SSLP) markers Fo16.96M and Fo17.87M by analyzing 96 F2 plants with bright green glossy stems in a mapping population (Fig. 1A). To narrow down the region, we further analyzed 2,000 F2 individuals to locate the mutation to a 51-kb region between Fo17.731M and Fo17.7826M (Fig. 1A). DNA sequencing of this region revealed that one ;21 kb fragment between 17750256 and 17771918 bp was deleted in the genomic DNA of cer16-1. This deletion caused the disruption of AT5G44100 and complete loss of AT5G44110, AT5G44120, AT5G44130, AT5G44140, and AT5G44150 (Fig. 1B). The results were further confirmed by PCR analysis using the genomic DNA of the cer16-1 mutant as the template (Supplemental Fig. S1A).
To ascertain which gene mutation caused the cer16 phenotype, we transformed the intact genomic DNA containing the 2-kb upstream promoter, the full open reading frame, and the 1-kb downstream terminator of the six genes from AT5G44100 to AT5G44150 into cer16-1, respectively. Multiple cer16-1 transgenic complementation lines were generated for each gene. Only introduction of the genomic DNA of AT5G44150 completely rescued the cer16-1 glossy stem phenotype ( Fig. 1C; Supplemental Fig. S1, B and C). Two of the cer16-1 transgenic lines, whose phenotype was complemented by genomic DNA of AT5G44150, were selected for further analysis and designated as cer16-1R1 and cer16-1R2. Reverse transcription-PCR (RT-PCR) results showed that the transcripts of AT5G44150 were absent in cer16-1, but expressed at similar levels as in wild-type with cer16-1R1 and cer16-1R2 (Fig. 1D). We also generated two independent knockout lines of AT5G44150 in the Arabidopsis ecotype Columbia (Col-0) background with the clustered regularly interspaced short palindromic repeats (CRISPR)/Cas9 strategy, designated as cer16-2 and cer16-3. Sequence analysis revealed that the mutations were in the third exon of AT5G44150: a single nucleotide insertion in cer16-2 and a four nucleotides deletion in cer16-3, respectively (Supplemental Fig. S2, A and B). These caused frame shift of the AT5G44150 coding sequence. As was true for cer16-1, both cer16-2 and cer16-3 had a phenotype similar to wild type except for the visible glossy stem phenotype (Supplemental Fig. S2C).
The wax crystals on the surface of the inflorescence stem of plants carrying different alleles of cer16 and their corresponding wild type were analyzed by scanning electron microscopy (SEM). The wax crystal patterns on the stem surface of cer16-2 and cer16-3 were similar to that of cer16-1, including that the total number of visible crystals per area were similarly reduced, and that the crystal structures were more flattened than that of corresponding wild type ( Fig. 1E; Supplemental  Fig. S2D). The stem surfaces of complementation lines cer16-1R1 and cer16-1R2 were similar to wild-type Ler-0 (Fig. 1E). The wax chemical composition on inflorescence stems of the cer16 mutant lines, complementation lines, as well as corresponding wild-types were determined with gas chromatography with flame ionization detection. The total wax loads on inflorescence stems of cer16-1 were dramatically decreased to 53.4% of wildtype levels (Supplemental Table S1). This was mainly due to a general decrease of almost all major wax components, including the most abundant C29 alkane, C29 secondary alcohol, and C29 ketone, which were decreased to 40.8%, 53.7%, and 56.4% of wild-type levels (Supplemental Fig. S1D). However, the amount of C30 alcohol and C31 alkane in stem of cer16-1 greatly increased to 403.5% and 311.6% of wild-type levels (Supplemental Fig. S1D), which might be the result of shunting away from the impaired alkane forming pathway in cer16-1. Partial blockage of this pathway in cer16-1 likely causes the heavy accumulation of the C30-CoA wax precursor, which shunted to the alcohol forming pathway, or C32-CoA elongation used for Figure 1. Positional cloning of CER16 and complementation of cer16-1 mutant. A, CER16 fine mapping. The CER16 mutation was fine mapped to a 51-kb region between simple sequence length polymorphism markers Fo17.731MB and Fo17.7826MB. B, Genes in the deletion region of cer16-1 mutant. All genes in the deletion region are listed. C, Stems from 6-week-old Ler-0, cer16-1, cer16-1R1, and cer16-1R2. D, RT-PCR analysis of CER16 expression in stems of Ler-0, cer16-1, cer16-1R1, and cer16-1R2. E, Scanning electron micrographs of wax crystals on stems of Ler-0, cer16-1, cer16-1R1, and cer16-1R2. Scale bars 5 3 mm.
C31 alkane synthesis. The stem wax defects of cer16-1 were fully rescued in both the cer16-1R1 and cer16-1R2 complementation lines (Supplemental Fig. S1D; Supplemental Table S1). Consistent with the wax chemistry data for cer16-1, the total stem wax loads and compositions of wax constituents of cer16-2 and cer16-3 were similar to that of cer16-1 (Supplemental Fig. S2E; Supplemental Table S1). These results demonstrated that the mutation of AT5G44150 caused the wax-deficient phenotype in cer16-1 stem. Sequence analysis showed that CER16/AT5G44150 encodes a 355-amino acid peptide and does not contain any annotated domains or motifs.
The branched chain wax constituents, including branched chain alcohols and alkanes, are deficient in the total wax of leaves and flowers of cer16-1, respectively (Reisberg et al., 2013;Busta and Jetter, 2017). However, we found that the deficiency of these wax compounds was not rescued in two complementation lines cer16-1R1 and cer16-1R2 (Supplemental Fig. S3), suggesting that the AT5G44150 is not responsible for the branched chain wax-deficient phenotype in cer16-1. To further confirm this, we analyzed the wax chemical composition in leaves of cer16-2 and cer16-3 as well as a mixed sample from 20 F2 plants with the glossy stem phenotype in the cer16-1 mapping population. Levels of C30 and C32 branched chain alcohols were similar to wild-type levels (Supplemental Fig. S3; Supplemental Table S1), which further confirmed that the branched chain wax constituent-deficient phenotype of cer16-1 leaves was not caused by mutation of AT5G44150, but likely by the inactivity of another gene.

CER16 Is Present in Land Plants But Not in Alga
A homologs search for CER16 in Arabidopsis databases indicates that CER16 exists as a single copy gene in Arabidopsis. To further explore the evolutionary history of CER16 in green plants, BLAST searching was performed to recover the CER16 homologs from all major lineages of green plants (Supplemental Table S2). We found no positive blast hits from any of the algae genomes, even when we used a high E-value cutoff (E-value 5 10). Homologs were recovered from angiosperms, gymnosperms, lycophytes, ferns, and mosses. Phylogenetic analysis grouped all the homologs into several clades, which corresponded to the main lineages of land plants (Fig. 2). In most plant lineages investigated here, CER16 existed as a single copy, but multiple copies were present in several polyploid species, for example, flax (Linum usitatissimum), tomato (Solanum lycopersicum), and soybean (Glycine max; Supplemental Table S2). BLAST searching and phylogenetic analysis revealed that CER16 is present in all land plants lineages but not in algae (Fig. 2), and we propose that CER16 is unique to land plants.

CER16 Localizes to Endoplasmic Reticulum and Is Ubiquitously Expressed
Wax components are produced in the endoplasmic reticulum (ER), and many wax biosynthesis enzymes are localized in this organelle. Therefore, we suspected that CER16 might localize in the same compartment. We constructed Pro-35S::CER16-YFP and transiently coexpressed it with the ER membrane marker CD3-959 in Nicotiana benthamiana leaves. In these leaves, the signals of CER16-YFP colocalized with CD3-959, indicating that CER16 localizes to ER (Fig. 3A). To further confirm this result, we transformed the Pro-35S::CER16-YFP construct into cer16-2. This construct complemented the waxdeficient phenotype of cer16-2 (Supplemental Fig. S4), indicating that CER16-YFP protein was functional. Individual roots of transgenic line were photographed with confocal laser-scanning microscopy. Stable CER16-YFP signal was found in a reticulate network, the typical signal of the ER protein. When the roots were stained with ER dye (hexyl rhodamine B), the CER16-YFP showed a clear colocalization with the stain (Fig. 3B), suggesting that CER16 localizes to ER.
Many wax related genes are highly expressed in specific regions, particularly in epidermal cells. Thus, we wished to investigate if expression of CER16 is temporally and spatially regulated. We constructed ProCER16::GUS transgenic lines and assessed GUS activity at different development stages and different plant parts. In 3-, 5-, and 14-d-old seedlings, GUS gene was mainly expressed in taproots and lateral roots, with the strongest expression detected in root tips (Supplemental Fig. S5, A-D). In 4-week-old plants, GUS expression was mainly detected in leaf veins and in regions undergoing active cell division and elongation (Supplemental Fig. S5E). In mature plants, strong GUS activity was detected in all tissues including cauline leaves, inflorescence, flowers, and siliques (Supplemental Fig. S5, F-I). When the inflorescence stem was cross-sectioned, GUS expression was found in internal cell layers, and we conclude that it was not restricted to the epidermal cells (Supplemental Fig. S5J).

CER16 Deficiency Silenced CER3 Expression
Until now, the function of CER16 has not been investigated in detail, and its exact role in wax biosynthesis is unknown. To examine its function in cuticular wax biosynthesis, we performed comparative transcriptomic analysis. The profiles of gene expression were compared in 6-week-old inflorescence stem of wild-type Col-0 and cer16-2. We identified 149 differently expressing genes (fold change $2 and false discovery rate ,0.01). There were 105 upregulated and 44 downregulated genes in cer16-2 compared with wild type (Supplemental Table S3). A Gene Ontology term enrichment analysis revealed that the differentially expressed genes were involved in diverse biological processes (Supplemental Fig. S6). We screened for genes associated with wax biosynthesis pathways, and observed that, of these genes, only CER3 transcripts were significantly reduced in the cer16-2 mutant (Supplemental Table S3), which was further confirmed by quantitative RT-PCR analysis (RT-qPCR; Supplemental Fig. S7). To confirm that CER3 expression was associated with CER16, we also checked CER3 expression in cer16-1 and cer16-3 and found that CER3 expression was also reduced in these mutants (Fig. 4A). CER3 expression was also increased in cer16 complementation lines, and exceeded the wild-type expression (Fig. 4A). These results confirmed that inactivation of CER16 does indeed suppress expression of CER3.
To confirm that wax deficiency in cer16-2 was due to reduced CER3 expression, we overexpressed CER3 in cer16-2 mutant under the control of the epidermalspecific CER6 promoter, which we also introduced into Col-0. The levels of CER3 transcripts in Col-0 containing ProCER6::CER3 construct were significantly increased compared with wild-type (Supplemental Fig. S8), indicating that the ProCER6:CER3 is functional. However, CER3 transcripts were not up-regulated in cer16-2 mutants carrying the CER3 construct (Supplemental Fig.  S8), which produced similar amounts of stem wax as that of cer16-2 lacking the CER3 construct (Supplemental Fig. S9), suggesting that the exogenous CER3 transgene transcripts were silenced in cer16-2.

CER16 Deficiency Triggers the PTGS Pathway
We noticed that cer16 phenocopied cer7, and CER3 expression was reported to be regulated by CER7 via an siRNA-mediated PTGS pathway (Hooker et al., 2007;Lam et al., 2015). Based on this, we speculated that down-regulation of CER3 in cer16 mutants might also be caused by accumulation of siRNAs that targeted CER3. We examined the levels of siRNA-1 and siRNA-2 (Lam et al., 2015) that targeted CER3 in inflorescence stems by RT-qPCR, and found that the levels of these two siRNAs in cer16-1, cer16-2, and cer16-3 were much higher than that in Col-0, Ler-0, and the cer16-1 complementation lines ( Fig. 4B; Supplemental Fig. S10).
that targeted CER3, we mutated RDR1 or RDR6 in a cer16-2 background by CRISPR/Cas9 to generate cer16-2 rdr1 and cer16-2 rdr6 double mutants. RDR1 and RDR6 are two major components of PTGS pathway reported to suppress CER3 silencing in cer7 (Lam et al., 2012;Lam et al., 2015). Sequence analysis of cer16-2 rdr1 revealed that 521 nucleotides were deleted in RDR1 genomic region (Supplemental Fig. S11, A and B). The cer16-2 rdr1 exhibited a waxy stem, and the total wax load on cer16-2 rdr1 stems was increased to 113.2% of wild-type levels (Fig. 5A). This was mainly due to increased C30 aldehyde and C29 alkane (up to 167.3% and 111.4% of wild type, respectively; Supplemental Fig. S12). Sequence analysis of RDR6 genomic fragment in cer16-2 rdr6 revealed that 29 nucleotides were replaced by ACCCC near the Cas9 target (Supplemental Fig. S11, C and D). The cer16-2 rdr6 showed obvious downward curly leaf margins, a similar phenotype described in rdr6 mutants (Peragine et al., 2004). The wax load of stem in cer16-2 rdr6 was restored to approximate wild-type levels, with C30 aldehydes in cer16-2 rdr6 at 157.8% of wild-type levels ( Fig. 5A; Supplemental Fig. S12). The levels of CER3 transcripts, siRNA-1, and siRNA-2 were also examined in cer16-2 rdr1 and cer16-2 rdr6. As shown in Figure 5, B and C, the two siRNAs were accumulated to high levels in cer16-2, whereas they were decreased to wild-type levels in cer16-2 rdr1 and cer16-2 rdr6. The expression levels of CER3 in cer16-2 rdr1 and cer16-2 rdr6 were increased 2.7-and 2.2-fold compared with the wild type, respectively. This was consistent with the significant increase of C30 aldehyde observed in cer16-2 rdr1 and cer16-2 rdr6. These results indicated that the wax-deficient phenotype and CER3 silencing in cer16-2 was suppressed by the mutation of RDR1 or RDR6, implicating that the CER16 mutation serves as a trigger that activates the PTGS pathway in which RDR1 and RDR6 participate.

DISCUSSION
The cer16-1 wax deficient mutant was described long ago, but the gene that confers the phenotype was unknown. The cer16-1 mutation was induced by fast neutron treatments, and we observed that several genes are deleted in cer16-1. We found that alkanes are reduced in the cer16-1 mutant's stem and that the branched chain components (only present in leaf and flower) were also reduced in cer16-1 leaves. However, genetic analysis involving additional cer16 alleles revealed that the observed reduction in the branched chain components of the cer16-1 mutant is due to mutation of a different gene from the CER16. The reduced alkanes resulted from a mutation in AT5G44150, whereas the reduced branched component phenotype was caused by gene mutation roughly mapped to chromosome 1 (Xianpeng Yang, Shipeng Li, and Shiyou Lü, unpublished data). In this study, we sought to elucidate the function of the alkane associated AT5G44150 and characterized the phenotypes of all cer16 mutants representing multiple mutant Figure 5. RDR1 and RDR6 are required for CER16-mediated CER3 silencing. A, Total wax load of Arabidopsis inflorescence stems of Col-0, cer16-2, cer16-2 rdr1, and cer16-2 rdr6. Wax coverage is expressed as microgram per square decimeter stem surface area. Values shown are means 6 SD (n 5 4). **P , 0.01, ***P , 0.001, Student's t test. B, RT-qPCR analysis of CER3 expression levels in stems of Col-0, cer16-2, cer16-2 rdr1, and cer16-2 rdr6. ACTIN2 was used as the internal control, and wild-type Col-0 was normalized to 1. Data are mean 6 SD (n 5 4). ***P , 0.001, Student's t test. C, RT-qPCR analysis of levels of siRNA-1 and siRNA-2 in stems of Col-0, cer16-2, cer16-2 rdr1, and cer16-2 rdr6. U6 was used as the internal control, and wild-type Col-0 was normalized to 1. Data are mean 6 SD (n 5 4). *P , 0.05, **P , 0.01, ***P , 0.001, Student's t test.
alleles of AT5G44150. All cer16 mutants exhibited a major reduction in alkane levels compared with wildtype, results that clearly reveal this gene's function in alkane biosynthesis.
We identified the organellar and cellular expression patterns of CER16. It was ubiquitously expressed in all tissues of Arabidopsis, and was most highly expressed in regions undergoing active cell division and differentiation, such as the rapidly expanding tips of inflorescence stems, root tips, and young leaves (Supplemental Fig. S5). These locations typically have high expression of cuticle biosynthesis genes (Bernard and Joubès, 2013). CER16 was localized to the ER (Fig. 3), where previous studies demonstrate that wax biosynthesis originate in this compartment (Bernard and Joubès, 2013). It seems that the expression patterns of CER16 are closely aligned with its assigned function in cuticular wax synthesis. In addition, our homology search reveals that CER16 is present only in land plants, perhaps because its evolution was essential for the formation of the protective hydrophobic cuticle that would have been required for plants first colonizing dry land.
Two findings here provide strong evidence that the PTGS pathway is induced in cer16 mutants. First, CER16 deficiency silenced both endogenous and transgene CER3 expression. Second, a mutation of two of the main components of the PTGS pathways restored CER3 expression in a cer16 background. It is unknown, however, how CER16 deficiency activates the PTGS pathway. Also unknown is the basis for regulation of CER3 expression by both CER16 and CER7. Recent work by Lange et al. (2019; https://www.biorxiv.org/content/ 10.1101/617894v1) may provide an important clue to CER16 function. In a biochemical study, these workers identified that CER16, designated as RIPR in their report, is a RNA exosome cofactor that can physically and functionally link with the exosome, and another RNA exosome cofactor designated RESURRECTION 1 (RST1) whose mutation causes wax-deficiency (Chen et al., 2005). Impairment of the RNA quality control systems usually trigger RNA silencing pathways (reviewed by Liu and Chen, 2016). RNA quality control and RNA silencing together coordinate control of the correct partitioning of RNA substrates to regulate plant growth. It is not surprising, then, that mutations in the exosome subunits and cofactors cause the similar phenotypes exhibited by cer7, rst1, and cer16. Although involvement of the 39 to 59 RNA decay system in the regulation of wax biosynthesis has been clearly demonstrated, it is still unknown if XRN exonucleases that mediate 59 to 39 RNA degradation also participate in this process. Zhang et al. (2015) reported that suppression of both systems usually causes very severe phenotypes, suggesting cross talk between them. Future work on exploring the roles of XRN proteins in wax biosynthesis may resolve these questions, and lead to a more comprehensive understanding of the regulatory roles of RNA quality control in wax biosynthesis.
The appearance of cuticle is presumed critical for land plant adaptation to the highly desiccating, and highly solar radiated, terrestrial environment thought to have occurred some 450 million years ago (Edwards, 1993;Wood, 2005;Rensing et al., 2008). Characterization of PpABCG7 in moss revealed that genes related to cuticular wax transport were functionally conserved in land plant evolution (Buda et al., 2013). CER3 was the core component of alkane biosynthesis, and its homologs were found in early land plants, including liverworts, and mosses (spreading earthmoss [Physcomitrella patens], flat-topped bogmoss [Sphagnum fallax], and spikemoss [Selaginella rupestris]; Wang et al., 2019). Even in Klebsormidium flaccidum, CER3 was the only gene identified in alkane biosynthesis (Kondo et al., 2016), suggesting that CER3 might be an ancestral gene in the alkane production for cuticular wax. Likewise, RNA silencing is evolutionally conserved in eukaryotic organisms (Eamens et al., 2008;Martínez de Alba et al., 2013); therefore, it may be that the PTGS regulation mechanism of CER3 is also evolutionally conserved. Apart from CER3, CER1, another core member of alkaneforming complex, is also regulated at the posttranscriptional level by small RNA, miR156 (Li et al., 2019). How the two classes of small RNAs act on balancing the expression of CER1 and CER3 for alkane production is still unknown.
To our knowledge, no member of the CER16 family has been functionally characterized previously; therefore, CER16 represents a novel gene family involved in biosynthesis of cuticular wax. It is present in all land plants lineages, but not in algae (Fig. 2). This is unlike the CER7 and SKI proteins whose homologs were also previously found in yeast (Zinder and Lima, 2017). These findings suggest that the emergence of CER16, together with other cuticle wax related synthetic mechanisms, may have facilitated the territorial colonization of land plants. Therefore, future work to characterize the related CER16 homologs in other plant lineages, especially the early land plants, will be of great interest.

Map-based Cloning of CER16 and Complementation of cer16-1 Mutant
The positional cloning of CER16 was performed as previously described by Lukowitz et al. (2000). To complement the cer16-1 mutant, the genomic fragments of AT5G44100 to AT5G44150 were amplified from genomic DNA of Arabidopsis Col-0 and induced into pMDC83 vector by Hieff Clone Plus One Step Cloning Kit (Yeasen), respectively. Primers sequenced here are listed in Supplemental Table S4. All constructs were introduced into cer16-1 by Agrobacterium tumefaciens-mediated transformation (Clough and Bent, 1998). The transformed seeds were collected and screened for resistance to Hygromycin B, and plants of the T2 generation were used for phenotype analysis.

Cuticular Wax Analysis
The cuticular wax crystallization patterns on Arabidopsis inflorescence stems were imaged by scanning electron microscopy. The dried samples were coated with gold for 2 min in a sputter coater (MC1000, Hitachi) and viewed by a scanning electron microscope (TM3030, Hitachi). The inflorescence stems and rosette leaves of 6-week-old Arabidopsis were collected for wax analysis as previously described by Yang et al. (2017).

Phylogenetic Analysis
To explore the origin and diversification of CER16 in green plants, we performed blast search against the published genomes including algae, bryophytes, ferns, gymnosperms, and angiosperms (Supplemental Table S2). A high E-value cutoff (10) was used to maximize the sensitivity of searches in order to identify short and incomplete sequences. To infer the phylogeny of CER16, we conducted an iterative set of alignment and phylogenetic estimation steps. Initial alignments were performed with MAFFT v7.407 (Katoh and Standley, 2013) using default settings and trimmed with Trimal v1.2 (-gt 0.05 -st 0.0001; Capella-Gutiérrez et al., 2009). The trimmed sequences were realigned using MAFFT by using a variable scoring matrix (-ginsi,-allowshift-unalignlevel a max 5 0.8) to minimize the over aligned issues observed in previous MAFFT estimation (Katoh and Standley, 2016). The best-fit model was chosen according to the BIC (Bayesian Information Criterion) score (Schwarz, 1978) computed by modelFinder (Kalyaanamoorthy et al., 2017). Phylogenetic analyses were carried out by RAxML-NG v0.7.0b (Kozlov et al., 2019) with 10 parsimony starting trees, and 200 bootstrap replicates (-bs-metric tbe) were used to evaluate variation in the dataset.

Subcellular Localization of CER16 Protein
The Pro-35S::CER16-YFP construct was generated for subcellular localization of CER16 protein by expressing CER16-YFP in Nicotiana benthamiana leaves and stable Arabidopsis transgenic plants. The coding region of CER16 was first subcloned into pDONR207 by the BP (recombination of an attB substrate with an attP substrate) reaction and then transferred into pEarleyGate101 by LR (recombination of an attL substrate with an attR substrate) reaction. The construct was cotransformed with ER membrane Marker CD3-959 (Nelson et al., 2007) into N. benthamiana leaves (Sparkes et al., 2006). Photographs were obtained by confocal laser scanning microscopy (Leica, TCS SP8 MP) equipped with hybrid (HyD) detectors. This construct was then introduced into cer16-2 by A. tumefaciens-mediated transformation (Clough and Bent, 1998). The transformed seeds were collected and screening for resistance to herbicide BASTA, and the roots of T2 generation lines were dyed with hexyl rhodamine B for confocal microscope imaging as previously described by Haslam et al. (2015).

Promoter-GUS Reporter Gene Fusion and GUS Histochemical Staining
To generate ProCER16::GUS construct, 2-kb upstream fragments of CER16 was cloned from Arabidopsis genomic DNA and induced into pMDC162 by KpnI and AscI. The ProCER16::GUS was then introduced into Arabidopsis Col-0 by A. tumefaciens-mediated transformation (Clough and Bent, 1998), and T3 generation lines were then used for GUS histochemical analysis as previously described by Yang et al. (2017).
Ectopic Expression of CER3 in Col-0 and cer16-2 Driving by CER6 Promoter The 35S promoter in pMDC83 was replaced by CER6 promoter (Haslam et al., 2015), and then the coding region of CER3 was inserted in the downstream of CER6 promoter to generate ProCER6::CER3 construct, which was then introduced into Col-0 and cer16-2. The transformed seeds were collected and screening for resistance to Hygromycin B, and T3 generation lines were used for phenotype analysis.

RT-PCR and RT-qPCR
RT-PCR and RT-qPCR were performed to analyze gene expression levels in Arabidopsis. Total RNAs were extracted from 6-week-old inflorescence stems of mutants, transgenic lines, and wild type by Universal Plant Total RNA Extraction Kit (BioTeke). First strand complementary DNAs (cDNAs) were synthesized from 1 mg of total RNAs with the ImProm-II Reverse Transcription System (Promega) and oligo(dT) 15 . ACTIN2 was used as the internal reference, and primers used in this study were listed in the Supplemental Table S4.

Small RNA Analysis by RT-qPCR
Total small RNA was extracted from 6-week-old inflorescence stems using RNAiso for Small RNA Kit (Takara, Dalian, China) according to the manufacturer's instruction. First strand cDNAs were synthesized from total small RNA (0.5 mg) using micro RNA 1st Strand cDNA Synthesis Kit (Vazyme) with specific primers listed in Supplemental Table S4. RT-qPCR was performed using micro RNA Universal SYBR qPCR Master Mix (Vazyme), and U6 was used as the internal reference (Turner et al., 2013).

RNA Sequencing Analysis
Six-week-old inflorescence stems were collected from Col-0 and cer16-2 and used for total RNA extraction. Total RNA was extracted by TRIzol reagent (Invitrogen Life Technologies) and prepared for library construction as previously described by Chen et al. (2015). The library was further sequenced with an Illumina HiSeq 2500 platform and generated 150 bp single-end reads. The deep sequencing data from RNA sequencing were further analyzed as previously described by Zhang et al. (2015).

SUPPLEMENTAL DATA
The following supplemental materials are available.
Supplemental Figure S3. AT5G44150 is not responsible for the branched chain wax deficiency in cer16-1.
Supplemental Figure S5. CER16 expression pattern in different tissues.
Supplemental Figure S6. Gene Ontology (GO) term enrichment analysis of differentially expressed genes in stems of Col-0 and cer16-2.
Supplemental Table S1. Cuticular wax composition of inflorescence stems and leaves of Arabidopsis.
Supplemental Table S2. Genome annotation data and transcriptomes used for BLAST search in this study.
Supplemental Table S3. List of genes differentially expressed in stems of cer16-2 relative to Col-0.
Supplemental Table S4. Primers used in this study.