# ORGAN MODIFICATION FOR EDIBLE PARTS OF HORTICULTURAL CROPS

EDITED BY : Yuke He, Guusje Bonnema, Han Xiao and Yunde Zhao PUBLISHED IN : Frontiers in Plant Science

#### Frontiers Copyright Statement

© Copyright 2007-2019 Frontiers Media SA. All rights reserved. All content included on this site, such as text, graphics, logos, button icons, images, video/audio clips, downloads, data compilations and software, is the property of or is licensed to Frontiers Media SA ("Frontiers") or its licensees and/or subcontractors. The copyright in the text of individual articles is the property of their respective authors, subject to a license granted to Frontiers.

The compilation of articles constituting this e-book, wherever published, as well as the compilation of all other content on this site, is the exclusive property of Frontiers. For the conditions for downloading and copying of e-books from Frontiers' website, please see the Terms for Website Use. If purchasing Frontiers e-books from other websites or sources, the conditions of the website concerned apply.

Images and graphics not forming part of user-contributed materials may not be downloaded or copied without permission.

Individual articles may be downloaded and reproduced in accordance with the principles of the CC-BY licence subject to any copyright or other notices. They may not be re-sold as an e-book.

As author or other contributor you grant a CC-BY licence to others to reproduce your articles, including any graphics and third-party materials supplied by you, in accordance with the Conditions for Website Use and subject to any copyright notices which you include in connection with your articles and materials.

All copyright, and all rights therein, are protected by national and international copyright laws.

The above represents a summary only. For the full conditions see the Conditions for Authors and the Conditions for Website Use. ISSN 1664-8714 ISBN 978-2-88963-103-2 DOI 10.3389/978-2-88963-103-2

#### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

#### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

#### Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

#### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office: researchtopics@frontiersin.org

# ORGAN MODIFICATION FOR EDIBLE PARTS OF HORTICULTURAL CROPS

Topic Editors:

Yuke He, Shanghai Institute of Plant Physiology and Ecology, Chinese Academy of Sciences, China Guusje Bonnema, Wageningen University and Research, Netherlands Han Xiao, Shanghai Institute of Plant Physiology and Ecology, Chinese Academy

of Sciences, China Yunde Zhao, University of California, San Diego, United States

Citation: He, Y., Bonnema, G., Xiao, H., Zhao, Y., eds. (2019). Organ Modification for Edible Parts of Horticultural Crops. Lausanne: Frontiers Media. doi: 10.3389/978-2-88963-103-2

# Table of Contents


Jing Yu, Liwei Gao, Wusheng Liu, Lixiao Song, Dong Xiao, Tongkun Liu, Xilin Hou and Changwei Zhang


Mark D. Lazzaro, Shan Wu, Ashley Snouffer, Yanping Wang and Esther van der Knaap

*41 Genetic Analysis of Chinese Cabbage Reveals Correlation Between Rosette Leaf and Leafy Head Variation*

XiaoXue Sun, Shuangxia Luo, Lei Luo, Xing Wang, Xueping Chen, Yin Lu, Shuxing Shen, Jianjun Zhao and Guusje Bonnema


Lingjie Wu, Zhendong Tian and Junhong Zhang


# Editorial: Organ Modification for Edible Parts of Horticultural Crops

Yuke He<sup>1</sup> \*, Guusje Bonnema<sup>2</sup> , Han Xiao<sup>1</sup> and Yunde Zhao<sup>3</sup>

*<sup>1</sup> National Key Laboratory of Plant Molecular Genetics, Shanghai Institute of Plant Physiology and Ecology, Chinese Academy of Sciences, Shanghai, China, <sup>2</sup> Plant Breeding, Wageningen University and Research, Wageningen, Netherlands, <sup>3</sup> Section of Cell and Developmental Biology, University of California, San Diego, San Diego, CA, United States*

Keywords: morphological modification, corm, bulb, curd, fleshy fruit, fleshy root, leafy head, tuber

**Editorial on the Research Topic**

#### **Organ Modification for Edible Parts of Horticultural Crops**

Plant development refers to various cellular processes that govern the morphogenesis of root, stem, leaf, flower, fruit, and seed. Plant development is regulated by genetic, epigenetic, and environmental factors. Developmental processes includes the determination, initiation, differentiation, and expression of different types of primordia on shoot/root apical meristems, and involve the establishment of polarity, cell division and differentiation of cell and tissue, phase transition, and hormone responses. Plant organs such as leaf, stem, root, flower, and fruit plays roles in water and mineral absorption, transportation, photosynthesis, pollination, fertilization, and other physiological processes. However, many horticultural crops often display enlarged organs to develop into the edible products: leafy head, bulb, tuber, fleshy stem, corm, fleshy root, root tuber, curds, and fleshy fruits. Unlike the grains of many field crops that mainly provide starch, the modified organs in horticultural crops store vitamins, secondary metabolites, minerals, and dietary fiber that are important for human health. Their commercial quality mainly depends on the diversity in size, shape, surface features, and texture. It is very interesting to know how these plant organs were selected and domesticated to become the edible parts of horticultural crops. Recently, molecular biology technology have been applied to address the key scientific questions and research directions with regard to morphological modifications of leaf, stem, root, flower, and fruit. For example, several miRNAs have been identified to play essential roles in determination of size, shape, and timing of leafy heads (Mao et al., 2014; Wang et al., 2014; Ren et al., 2018).

The chapters address genetic and molecular basis of morphological modifications of plant organs in horticultural crops. The focused topic issues cover dynamic phenotyping of the modified organs: leafy head, bulb, tuber, fleshy stem, corm, fleshy root, root tuber, curd, and fleshy fruit during plant development, genetic variation of modified organs, functional analysis of genes of organ-related traits that control size, shape, weight, texture, color and flavor of the modified organs, regulation of non-coding RNAs controlling morphological modification, effects of abiotic stress molecular mechanisms underlying metamorphosis of plant organs.

Genetic basis of leafy heads and bulbs is a critical topic. Many vegetable crops develop into certain types of leaf curvature required for high yield and quality. Recent research has revealed that miRNAs in Brassica crops regulate the timing and shape of leafy heads (Mao et al., 2014; Wang et al., 2014). Several papers in this issues described the origin, diversity, and development of leafy heads. Three papers ("Genetic analysis of Chinese cabbage reveals correlation between rosette leaf and leafy head variation" by Sun et al.; "Characterization of non-heading mutation

Edited and reviewed by: *Roberto Papa, Marche Polytechnic University, Italy*

> \*Correspondence: *Yuke He heyk@sippe.ac.cn*

#### Specialty section:

*This article was submitted to Plant Breeding, a section of the journal Frontiers in Plant Science*

Received: *14 June 2019* Accepted: *09 July 2019* Published: *23 July 2019*

#### Citation:

*He Y, Bonnema G, Xiao H and Zhao Y (2019) Editorial: Organ Modification for Edible Parts of Horticultural Crops. Front. Plant Sci. 10:961. doi: 10.3389/fpls.2019.00961*

**4**

in heading Chinese cabbage (Brassica rapa L. ssp. pekinensis)" by Li et al.; and "Transcription coactivator ANGUSTIFOLIA3 (AN3) regulates leafy head formation in Chinese Cabbage" by Yu et al.). These papers described the phenotyping, DNA resequencing and gene function identification of Chinese cabbage. Modification of flower organ is important for function transformation of edible organs. Two papers are related to flower organ modifications ("Differential gene expression caused by the F and M loci provides insight into ethylene-mediated female flower differentiation in cucumber" by Pan et al. and "defective apetala2 genes lead to sepal modification in Brassica crops" by Zhang et al.). Shape, size, color, and flavor of fleshy fruits are genetically related to leaf features. One paper "Chemical composition and crystal morphology of epicuticular wax in mature fruits of 35 pear (Pyrus spp.) cultivars" by Wu et al. "discloses crystal morphology of "epicuticular wax in mature fruits;"" and another paper "Functional dissection of Auxin Response Factors in regulating tomato leaf shape development" by Wu et al. provide the evidences that auxin response factors regulate tomato leaf shape and fruit development. The paper "Plant organ shapes are regulated by protein interactions and associations with microtubules" addressed the importance of microtubules to formation of organ shape (Lazzaro et al.). Finally, this book contains eight papers from scientists working in institutions from different countries. Thus, we can claim that this book, with its multiple scientific voices, provides also a contribution to understanding

#### REFERENCES


genetic basis and molecular mechanisms underlying organ modification.

# AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

#### FUNDING

This work was supported by grants from the National Key Research and Development Program of China (Grant No.2016YFD0101900), Natural Science Foundation of China (Grant Nos. 31771442 and 31571261), The National Key Research and Development Program of China (2016YFD0100500), Shanghai Agriculture Applied Technology Development Program, China (Grant No. Z20160109), and Shanghai Agriculture Applied Technology Development Program, China (Grant No.G2015060107), Dutch Royal Academy of Sciences China Exchange Program (Grant No. 530-4CDP08), the National Key R&D of China (Grant No. 2017YFD0101802), the National Natural Science Foundation of China (Grant No. 31672151), International Cooperation project in the Science and Technology Support Program of Hebei (Grant No. 17396315D), Hundred Innovative Talent Program of Hebei (Grant No. SLRC2017040), the Natural Science Foundation of Hebei (C2016204170), and the Science and Technology Support Program of Hebei (16226304D-2).

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 He, Bonnema, Xiao and Zhao. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Transcription Coactivator ANGUSTIFOLIA3 (AN3) Regulates Leafy Head Formation in Chinese Cabbage

Jing Yu<sup>1</sup>† , Liwei Gao<sup>1</sup>† , Wusheng Liu<sup>2</sup> , Lixiao Song<sup>3</sup> , Dong Xiao<sup>1</sup> , Tongkun Liu<sup>1</sup> , Xilin Hou<sup>1</sup> and Changwei Zhang<sup>1</sup> \*

<sup>1</sup> State Key Laboratory of Crop Genetics and Germplasm Enhancement, College of Horticulture, Nanjing Agricultural University, Nanjing, China, <sup>2</sup> Department of Horticultural Science, North Carolina State University, Raleigh, NC, United States, <sup>3</sup> Jiangsu Key Laboratory for Horticultural Crop Genetic Improvement, Nanjing, China

#### Edited by:

Yuke He, Institute of Plant Physiology and Ecology, Shanghai Institutes for Biological Sciences (CAS), China

#### Reviewed by:

Xiaowu Wang, Biotechnology Research Institute (CAAS), China Jianjun Zhao, Agricultural University of Hebei, China

\*Correspondence:

Changwei Zhang changweizh@njau.edu.cn

†These authors have contributed equally to the work

#### Specialty section:

This article was submitted to Plant Physiology, a section of the journal Frontiers in Plant Science

Received: 30 May 2018 Accepted: 04 April 2019 Published: 30 April 2019

#### Citation:

Yu J, Gao L, Liu W, Song L, Xiao D, Liu T, Hou X and Zhang C (2019) Transcription Coactivator ANGUSTIFOLIA3 (AN3) Regulates Leafy Head Formation in Chinese Cabbage. Front. Plant Sci. 10:520. doi: 10.3389/fpls.2019.00520 Leafy head formation in Chinese cabbage (B. rapa ssp. pekinensis cv. Bre) results from leaf curvature, which is under the tight control of genes involved in the adaxial-abaxial patterning during leaf development. The transcriptional coactivator ANGUSTIFOLIA3 (AN3) binds to the SWI/SNF chromatin remodeling complexes formed around ATPases such as BRAHMA (BRM) in order to regulate transcription in various aspects of leaf development such as cell proliferation, leaf primordia expansion, and leaf adaxial/abaxial patterning in Arabidopsis. However, its regulatory function in Chinese cabbage remains poorly understood. Here, we analyzed the expression patterns of the Chinese cabbage AN3 gene (BrAN3) before and after leafy head formation, and produced BrAN3 gene silencing plants by using the turnip yellow mosaic virus (TYMV)-derived vector in order to explore its potential function in leafy head formation in Chinese cabbage. We found that BrAN3 had distinct expression patterns in the leaves of Chinese cabbage at the rosette and heading stages. We also found silencing of BrAN3 stimulated leafy head formation at the early stage. Transcriptome analysis indicated that silencing of BrAN3 modulated the hormone signaling pathways of auxin, ethylene, GA, JA, ABA, BR, CK, and SA in Chinese cabbage. Our study offers unique insights into the function of BrAN3 in leafy head formation in Chinese cabbage.

Keywords: AN3, BRM, leafy head formation, VIGS, transcriptome analysis, Chinese cabbage

# INTRODUCTION

Heading Chinese cabbage (Brassica rapa ssp. pekinensis cv. Bre; 2n = 2× = 20) is one of the most important horticultural crops in China and, to a lesser extent, an oilseed crop (Zhao et al., 2005). Leafy head formation undergoes four developmental stages, i.e., the seedling, rosette, folding, and heading stages (He et al., 2000; Yu et al., 2000; Ke, 2010; Wang et al., 2014b). It forms uniform leafy heads with extremely incurved leaves surrounding the shoot apexes after the rosette stage. The leafy heading trait has been selected for several Brassica species including Chinese cabbage and cabbage (B. oleracea; 2n = 2× = 18) during crop domestication and breeding (Cheng et al., 2014). Wellformed leafy heads function normally for photosynthesis and serve as storage organs for essential

**6**

minimal nutrients, fibers, and vitamins, while poorly-developed heads significantly decrease the crop's commercial values (Yu et al., 2000, 2013; Cheng et al., 2016). Leafy head formation is affected by many factors such as low temperature, low light intensity, uneven auxin distribution, and low carbohydrate/nitrogen ratio (Ito, 1957). Plants have to respond to these factors properly for leaf bending and folding (i.e., leaf curvature), which results from the differential cell division and enlargement in leaf regions, and which are under the tight control of genes involved in the adaxial-abaxial patterning during leaf development (Mao et al., 2014).

Recent efforts had been taken to study the molecular mechanisms controlling leafy head formation in Chinese cabbage by using genetic mapping (Zhao et al., 2005; Ge et al., 2011; Yu et al., 2013; Zhang et al., 2013; Zhang X. et al., 2016), transcriptome profiling (Wang et al., 2012; Zhang C.W. et al., 2016; Gu et al., 2017; Li et al., 2018), proteomic analysis (Zhang C.W. et al., 2016), and miRNA expression profiling (Wang et al., 2013). Moreover, the effects of a few individual genes on leafy head formation had been studied in Chinese cabbage. These genes include the Agrobacterium genes Aux1 and Aux2 (He et al., 2000), the Chinese cabbage genes BcpLH (Yu et al., 2000), TCP (Mao et al., 2014), BrpSPL9 (Wang et al., 2014b), and auxin transport genes BrLAX, BrPIN, and BrPGP (Gao et al., 2017). Most of these genes are involved in the adaxial-abaxial patterning during leaf development in Chinese cabbage.

The Arabidopsis ANGUSTIFOLIA3/GRF-INTERACTIN FACTOR1 (AtAN3/AtGIF1), an important transcriptional activator of the GIF family, is also involved in the determination of leaf adaxial-abaxial patterning and growth (Horiguchi et al., 2011). Loss-of-function mutant of AtAN3 exhibited smaller and narrower leaves due to a decrease in cell number (Kim and Kende, 2004; Horiguchi et al., 2005), while ectopic overexpression of AtAN3 resulted in enlarged leaf size (Horiguchi et al., 2005; Lee et al., 2009). AtAN3 binds to the SWI/SNF chromatin remodeling complex formed around ATPases such as BRAHMA (BRM). By using the energy from ATP hydrolysis, the Arabidopsis AN3-SWI/SNF-BRM complex regulates gene expression by changing the interactions between histone octamers and the DNA for the access of transcription factors (Clapier and Cairns, 2009). The brm mutant had small spiral-shaped leaves with downward curling edges (Hurtado et al., 2006; Vercruyssen et al., 2014). The AN3-SWI/SNF-BRM complex also interacts with GROWTH REGULATING FACTOR (GRF) proteins, a class of plant-specific transcription activators in Arabidopsis (Kim and Kende, 2004; Liu et al., 2009; Debernardi et al., 2014). Ectopic overexpression of the GRFs increased leaf size in Arabidopsis due to enhanced cell proliferation (Kim et al., 2003; Horiguchi et al., 2005; Liu et al., 2009; Debernardi et al., 2014; Wang et al., 2014a). However, the regulatory function of AN3 in Chinese cabbage remains poorly understood.

In the present study, we explored the expression patterns of the Chinese cabbage AN3 gene (i.e., BrAN3) in different leaf locations of Chinese cabbage at the rosette and heading stages. We generated BrAN3-silencing Chinese cabbage plants by using the turnip yellow mosaic virus (TYMV)-derived vector, and examined the effect of BrAN3 silencing on the stimulation of leafy head formation. We also conducted transcriptome analysis of the BrAN3-silencing leaves and identified the differentially expressed genes (DEGs) caused by BrAN3 silencing. All of the results provide insights into the function of the BrAN3 gene in leafy head formation in Chinese cabbage.

# MATERIALS AND METHODS

#### Sequence Alignment and Phylogenetic Analysis

The protein sequences of the Arabidopsis AtAN3 and AtBRM genes were used individually as the query sequences to BlastP against the Brassica database<sup>1</sup> in order to obtain their homologous sequences in B. rapa. The protein sequences of the Arabidopsis AtAN3 and AtBRM genes were also used individually as the query sequences to BlastP against Phytozome<sup>2</sup> in order to obtain their homologous sequences in B. oleracea, Zea mays, and Oryza sativa. The genomic DNA sequences and cDNA sequences of these homologous sequences in these plant species were also obtained from Phytozome. All of the homologous protein sequences of each gene in different species were aligned together using CluastalX 2.1 (Larkin et al., 2007), and subjected to the phylogentic analysis using the MEGA5 program<sup>3</sup> and the neighbor-joining method with 1,000 bootstraps (Tamura et al., 2013).

#### Plant Material and Growth Conditions

Inbred line seeds of Chinese cabbage cv. Bre were germinated in soil in a tray with holes. Three weeks old seedlings were transferred to pots containing a mixture of nutrient soil and vermiculite (3:1) and grown in a growth chamber at 22◦C with a 16/8 h light/dark photocycle and 50% relative humidity.

# Gene Cloning and Plasmid Construction

Total RNA was extracted from the leaves of Chinese cabbage for gene cloning using an RNA extraction kit (TaKaRa; Dalian, Liaoning, China). DNase I treatment and first-strand cDNA synthesis were performed sequentially using the PrimeScript <sup>R</sup> 1st Strand cDNA Synthesis Kit (TaKaRa; Dalian, Liaoning, China). One microgram of total RNA was reverse-transcribed with the oligo(dT) primer according to the manufacturer's protocol. The cDNA fragments of BrAN3 and BrBRM were PCR amplified individually using the first-strand cDNA as the templates and the gene-specific primers (**Supplementary Table S1**). The PCR products were purified using an AxyPrep DNA Gel Extraction Kit (Axygen Biosciences; Union City, CA, United States) and then cloned into the pMD19-T vector (TaKaRa; Dalian, Liaoning, China), followed by Sanger sequencing. A gene-specific fragment of 40 nt in length was selected for each of the two genes to create an in-frame stop codon (TAA, TGA, or TAG) on the second, third or fourth amino acid position on each fragment due to frame shift. The gene-specific fragment was selected to target all

<sup>1</sup>http://brassicadb.org/brad/

<sup>2</sup>https://phytozome.jgi.doe.gov/pz/portal.html#!search?show=BLAST <sup>3</sup>http://www.megasoftware.net/

the homologous sequences of each gene in the Chinese cabbage genome. Each fragment and its reverse complementary sequence (**Supplementary Table S1**) were synthesized from TaKaRa (Dalian, Liaoning, China) and used to form a palindromic oligonucleotide dimer after self-hybridization, which was cloned into the plasmid pTY-S with the help of SnaBI (Pflieger et al., 2008). The resulting plasmids were named as pTY-BrAN3 and pTY-BrBRM, respectively.

### Real-Time RT-PCR

fpls-10-00520 April 26, 2019 Time: 15:38 # 3

Total RNA was extracted from the apical location, the lateral 1∼3 locations, and the base location of individual leaves and the whole leaves of Chinese cabbage (Gao et al., 2017) at both the rosette and heading stages, respectively, for the analysis of expression profiles of BrAN3 on different locations of Chinese cabbage leaves. As for the analysis of the effects of virus-induced gene silencing, total RNA was extracted individually from the Chinese cabbage leaves inoculated with each virus-induced gene silencing vector, the negative control vector, and positive control vector. The purity and integrity of the RNA were analyzed using a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies; Wilmington, DE, United States) and verified by gel electrophoresis. The real-time RT-PCR was conducted using the SYBR Green PrimeScript Plus RT-PCR Kit (TaKaRa; Dalian, Liaoning, China) on an ABI PrismR 7900HT (Applied Biosystems; Carlsbad, CA, United States) according to the manufacturer's instructions. Gene-specific primers were designed outside of the 40 nt long fragments for both BrAN3 and BrBRM genes as well as the BrAN3 downstream genes ARABIDOPSIS RESPONSE REGULATOR4 (ARR4) and CYTOKININ RESPONSE FACTOR2 (CRF2) (Vercruyssen et al., 2014) using the Beacon Designer 7.9 (Premier Biosoft International, Palo Alto, CA, United States) (**Supplementary Table S1**). These primers were designed to target all the homologous sequences of each gene in the Chinese cabbage genome. The amplification procedure was as follows: pre-denaturation at 94◦C for 10 s, followed by 40 cycles of 94◦C for 30 s and 60◦C for 30 s, and finally a melting curve was performed (95◦C for 15 s and at increments of 0.5◦C from 60 to 95◦C). Three biological and three technical replicates were carried out for each sample. The 2−11CT method (Schmittgen and Livak, 2008) was used for relative transcript quantification normalized by the internal control gene BrAct (Bra028615) as described (Dheda et al., 2004). The quality cutoff of the realtime RT-PCR efficiency was set at >95% with R <sup>2</sup> > 0.99 for the standard curve.

# Particle Bombardment-Mediated Delivery of the TYMV-Derived Vectors

The plasmids pTY-BrAN3 and pTY-BrBRM, the empty plasmid pTY-S and the positive control plasmid pTY-BrPDS (Yu et al., 2018) were transformed into the 6 weeks old Chinese cabbage seedlings for virus-induced gene silencing (VIGS) by particle bombardment as described in Yu et al. (2018). Leaves were harvested from three individual plants before and after head formation and kept in liquid nitrogen for RNA extraction (see above).

# Library Preparation for Transcriptome Sequencing

Total RNA was extracted individually from the leaves of two Chinese cabbage plants inoculated with each virus-induced gene silencing vector. The purity and integrity of the RNA was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies; Palo Alto, CA, United States) and the NanoPhotometer <sup>R</sup> spectrophotometer (IMPLEN; Westlake Village, CA, United States), while the concentration was measured using the Qubit <sup>R</sup> RNA Assay Kit in a Qubit <sup>R</sup> 2.0 Flurometer (Life Technologies; Carlsbad, CA, United States). Three µg RNA per sample was used for RNA sample preparations. Sequencing libraries were generated using the NEB Next <sup>R</sup> UltraTM RNA Library Prep Kit for Illumina <sup>R</sup> (NEB; Beverley, MA, United States) following the manufacturer's recommendations and index codes were added to attribute sequences to each sample. The library quality was assessed on the Agilent Bioanalyzer 2100 system (Agilent Technologies; Palo Alto, CA, United States).

# Clustering and Reads Alignment to the Reference Genome

The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3 cBot-HS (Illumina; Hayward, CA, United States) according to the manufacturer's instructions. Then, the libraries were sequenced on an Illumina Hiseq 2500 platform and 150 bp paired-end reads were generated. Reference genome and gene model annotation files were downloaded from Brassica database<sup>4</sup> . Index of the reference genome<sup>4</sup> was built and paired-end clean reads were aligned to the reference genome using HISAT2 (Kim et al., 2015).

The expression level was calculated as fragments per kilobase of exon model per million mapped (FPKM) values and gene expression pattern using Clusterv3.0<sup>5</sup> and the heat map of hierarchical clustering was established using TreeView v.3.0<sup>6</sup> based on the log2-converted FPKM values. Bubble Charts were performed using the OmicShare tools, a free online platform for data analysis<sup>7</sup> .

#### GO Enrichment and KEGG Pathway Analysis of Differentially Expressed Genes

Differentially expressed genes (DEGs) were defined if the False Discovery Rate is smaller than 0.01 and the Fold Change equals to or is larger than 2. Gene Ontology (GO) enrichment analysis of the DEGs was implemented by the GOseq R packages based on the Wallenius non-central hyper-geometric distribution (Young et al., 2010), which allows the adjustment for gene length bias in DEGs. KOBAS software was used to test the statistical enrichment of DEGs in KEGG pathways (Mao et al., 2005).

<sup>4</sup>http://brassicadb.org/brad/

<sup>5</sup>http://bonsai.hgc.jp/~mdehoon/software/cluster/

<sup>6</sup>http://jtreeview.sourceforge.net/

<sup>7</sup>www.omicshare.com/tools

FIGURE 1 | Comparison of the homologs of the AtAN3 protein in Brassica rapa, B. oleracea, Zea mays, and Oryza sativa. (A) Phylogenetic relationship and gene structures of the homologs of the AtAN3 protein in those species. The unrooted phylogenetic tree was generated by the neighbor-joining method using the MEGA5 program with 1,000 bootstraps. Lines represent the introns while solid boxes represent exons. (B) Comparison of the deduced amino acid sequences of the homologs of the AtAN3 protein in those species. Dashes represent gaps, dots represent the identical amino acid residues, and the numbers represent the length of the amino acid residues.

# Statistical Analysis

fpls-10-00520 April 26, 2019 Time: 15:38 # 5

Statistical analyses was performed using analysis of variance (ANOVA; SAS 9.4 for Windows; SAS Institute, Cary, NC; p<0.05).

# RESULTS

# Identification of the BrAN3 and BrBRM Genes in Chinese Cabbage

The BlastP search for the homologous sequences of the AtAN3 gene in the Brassica database and Phytozome returned three homologous genes in both B. rapa and B. oleracea, and one homologous gene in both Z. mays and O. sativa (**Figure 1**). Phylogenetic analysis revealed that the homologous genes in B. rapa and B. oleracea formed three monophyletic groups, indicating an origin of gene duplication (**Figure 1A**). The solo Arabidopsis homolog grouped with one of the three groups. The comparison between the cDNA sequences with the genomic DNA sequences revealed that all of these homologous genes have similar exon/intron structures (**Figure 1B**). Each homologous gene consists of four exons and three introns at relatively conserved positions. When using the protein sequence of the AtBRM gene as the query, we found two homologs in both B. rapa and B. oleracea, and one homolog in both Z. mays and O. sativa (**Supplementary Figure S1**). The homologous genes in B. rapa and B. oleracea formed two monophyletic groups, indicating a single gene duplication event (**Supplementary Figure S1**). All of these homologous genes have similar exon/intron structures

FIGURE 3 | Effects of gene silencing of BrAN3 and BrBRM in Chinese cabbage. (A) Representative images of the silencing plants by using the TYMV-derived vector pTY-S at 60 DPI. A gene-specific fragment of 40 nt in length was selected for each of the two genes to create an in-frame stop codon (TAA, TGA, or TAG) on the second, third or fourth amino acid position on each fragment due to frame shift. The gene-specific fragment was selected to target all the homologous sequences of each gene in the Chinese cabbage genome. Six weeks old Chinese cabbage plants were inoculated with each gene silencing vector. The mock treatment and the empty vector pTY-S were used as the negative control, while the pTY-BrPDS was used as the positive control. Leaf head formation was observed for pTY-BrAN3 and pTY-BrBRM-inoculated plants. (B) Relative expression levels of BrAN3, BrBRM, and two of BrAN3 downstream genes BrARR4 and BrCRF2 in the BrAN3-silencing plants measured by Real-time RT-PCR. RNA was extracted from the BrAN3-silencing leaves with the empty vector pTY-S-inoculated leaves as the control at 60 DPI. Relative expression levels were calculated using the standard curve method with BrAct (Bra028615; Dheda et al., 2004) as the internal control gene. Bars represent the means of three replicates ± standard errors (vertical bars). Asterisks indicate a significant difference between the two stages at p ≤ 0.05 as calculated by t-test.

(**Supplementary Figure S1**). The existence of the same numbers of homologous genes of both AtAN3 and AtBRM in the two Brassica genomes is consistent to the finding that parallel selection of homologous genes for leaf heading in the two species (Cheng et al., 2016).

# Distinct Expression Patterns of BrAN3 in the Leaves of Chinese Cabbage at the Rosette and Heading Stages

To explore the potential roles of the BrAN3 gene in head formation in Chinese cabbage, we conducted real-time RT-PCR analysis of the relative transcription levels of BrAN3 in the

Chinese cabbage leaves at both the rosette and heading stages. We found that BrAN3 had a constantly stable relative expression levels in different locations of the leaves as well as the whole leaf at the rosette stage except in the leaf apical tip where BrAN3 expression levels doubled (**Figure 2**). However, BrAN3 exhibited highly variable but significantly decreased relative expression levels in different locations of the leaves of Chinese cabbage at the heading stage (**Figure 2**). Their relative expression peaked in the leaf apical tip, followed by the leaf base, and the lateral 2 and 3 locations. The least relative expression was observed in the lateral 1 location. The significant variations on the relative expression in different locations of the Chinese cabbage leaf at the heading stage indicate a potential role in heading formation since the lateral parts affect the leaf curvature and leaf head formation than the mid-vein part (i.e., the lateral 1 location) does.

### The Effects of VIGS-Mediated Knockdown of BrAN3 and BrBRM Individually in Head Formation in Chinese Cabbage

To further explore the potential roles of the BrAN3 and its genetically interacting BrBRM gene in head formation in Chinese cabbage, VIGS was employed to knock down the expression of BrAN3 and BrBRM individually in Chinese cabbage. At 15 days post inoculation (DPI), all the TYMV-derived vector-inoculated plants displayed typical viral symptoms when compared with the untreated plants. These include plant stunting, leaf distortion, and mosaic symptoms. At 60 DPI, the untreated and the empty vector pTY-S-inoculated plants showed a normal increase in leaf size, while the positive control vector pTY-BrPDSinoculated plants exhibited a photobleaching phenotype as expected (**Figure 3A**; Yu et al., 2018). However, the pTY-BrAN3 and pTY-BrBRM-inoculated plants produced curled leaves and formed smaller leafy heads at the rosette age (**Figure 3A**), indicating an important role for both genes in head formation at the early stage.

To confirm the VIGS effect of BrAN3 in the BrAN3-silencing Chinese cabbage, real-time RT-PCR was performed to detect the relative expression levels of BrAN3 in pTY-BrAN3-inoculated plants, respectively. When compared to its relative expression levels on the empty vector pTY-S-inoculated plants, the relative expression levels of BrAN3 were significantly decreased on the pTY-BrAN3-inoculated plants (**Figure 3B**). In addition, we found that the relative expression levels of BrBRM significantly decreased in the BrAN3-silencing plants (**Figure 3B**), indicating BrAN3 positively regulates the expression of BrBRM.

### The Involvement of the Cytokinin Signaling Pathway in Head Formation in Chinese Cabbage

Since AtAN3 plays an important role in the regulation of cytokinin signaling to stimulate leaf cell proliferation and adaxialabaxial patterning, it represses and induces the expression of its

downstream genes AtARR4 and AtCRF2, respectively, which are cytokinin response regulators (Rashotte et al., 2006; Werner and Schmuelling, 2009; Holst et al., 2011; Vercruyssen et al., 2014). Thus, real-time RT-PCR was also conducted to detect the relative expression levels of BrARR4 and BrCRF2 in the BrAN3-silencing plants. As expected, the relative expression levels of BrARR4 and BrCRF2 significantly increased and decreased, respectively, in the BrAN3-silencing plants when compared to the empty vector pTY-S-inoculated plants (**Figure 3B**).

# RNA-Seq Analysis and Differentially Expressed Gene (DEG) Identification in the BrAN3-Silencing Chinese Cabbage

To further characterize the molecular mechanisms of BrAN3 in leafy head formation, RNA-Seq was performed on the BrAN3 silencing plants with the empty vector pTY-S-inoculated plants being used as the control. About 20,992,641 and 19,553,742 clean reads obtained from the BrAN3-silencing plants and the control plants, respectively, were mapped against the Brassica reference sequence, indicating equal match ratios on the reference genome for both groups of data. When compared to the control plants, a total of 1692 DEGs were detected in the BrAN3-silencing plants based on the transcript abundances; among which, 800 DEGs were significantly up-regulated and 892 DEGs were significantly down-regulated (**Figure 4**).

### Functional Classification of DEGs by Gene Ontology (GO) Enrichment Analysis and KEGG Pathway Analysis

Functional classification of the identified DEGs using the GO term enrichment analysis assigned all the genes and the DEGs to the three GO categories: biological processes, cell components, and molecular functions (**Figure 4**). Overall, DEGs was more enriched in biological processes category than that in the other two categories as well as that in all genes background. Specifically, the over-represented (p < 0.05) GO terms of biological processes category included cellular process, metabolic process, single-organism process, response to stimulus, biological regulation, developmental process, etc., indicating the silencing of BrAN3 may affect the developmental and biological processes. Furthermore, the KEGG pathway analysis of the identified DEGs showed that the DEGs involving in the pathways of photosynthesis-antenna proteins, nitrogen metabolism, phenylpropanoid biosynthesis, phenylalanine metabolism, and plant hormone signal transduction were significantly enriched in the BrAN3-silencing plants (**Figure 5**).

fpls-10-00520 April 26, 2019 Time: 15:38 # 7

We found there existed 41 DEGs in the plant hormone signal transduction term with 23 and 16 DEGs being significantly up-regulated and down-regulated, respectively (**Figure 6**). The remaining 2 DEGs showed controversial results, i.e., a significant increase in one BrAN3-silencing plant and a significant decrease in the other BrAN3-silencing plant (**Figure 6**). Specifically, the DEGs involved in the signaling pathways of ethylene (ETH; Bra023744 and Bra023746) and jasmonic acid (JA; Bra000421) were significantly increased in the BrAN3-silencing leaves when compared to the empty vector pTY-S-inovulated leaves

(**Figure 6**). The DEGs involved in the signaling pathways of gibberellin (GA; Bra039460 and Bra003520), brassinosteroid (BR; Bra002718, Bra002719, Bra020433, Bra012909, and Bra029992), and salicylic acid (SA; Bra004329 and Bra018634) were significantly decreased in BrAN3-silencing leaves in comparison to the empty vector-inoculated leaves (**Figure 6**). However, the DEGs involved in the signaling pathways of auxin (AUX), abscisic acid (ABA), and cytokinin (CK) showed opposite expression patterns. For example, 13 DEGs involved in the AUX signaling pathway (Bra011332, Bra014411, Bra023958, Bra027232, Bra032520, Bra033581, Bra005136, Bra017136, Bra008836, Bra009636, Bra030560, Bra023403, and Bra039186), 3 DEGs involved in the ABA signaling pathway (Bra019121, Bra031574, and Bra015579), and 4 DEGs involved in the CK signaling pathway (Bra025708, Bra019932, Bra001629, and Bra031714) were significantly increased in the BrAN3-silencing leaves (**Figure 6**). And 6 and 1 DEGs in the signaling pathways of AUX (Bra033890, Bra003044, Bra022663, Bra019060, Bra029023, and Bra032521) and CK [Bra019270 (BrCRF2)] showed a significantly decrease in the BrAN3-silencing plants (**Figure 6**). In addition, 1 DEG involved in each of the auxin (Bra041046) and ABA (Bra011815) signaling pathways exhibited controversial expression patterns when compared to the empty vector-inoculated leaves (**Figure 6**).

To validate the accuracy of the RNA-Seq data, we selected 9 out of the 41 DEGs from the plant hormone signal transduction term for real-time RT-PCR analysis. When compared to the empty vector pTY-S-inoculated leaves, we found that 7 (Bra011332, Bra014411, Bra023958, Bra027232, Bra005136, Bra017136, and Bra025708) out of the 9 DEGs were significantly increased while the other 2 DEGs (Bra003044 and Bra022663) were significantly decreased in the BrAN3-silencing plants (**Figure 7**). We also found the same expression patterns of these 9 DEGs were observed in the RNA-Seq data when calculated as the FPKM (**Figure 7**). Thus, the real-time RT-PCR and RNA-Seq data correlated very closely (**Figure 7**).

### DISCUSSION

Leafy head formation in Chinese cabbage is a tightly controlled process. Many proteins such as the SWI/SNF chromatin remodeling protein complex, ATPases, and transcription activators could play important roles in head formation. This is the first time to functionally confirm that the BrAN3 gene could induce head formation in Chinese cabbage. Its distinct expression patterns in different leaf locations and at different head developmental stages are consistent with its function in head formation. Similar expression patterns had also been observed for the BrTCP gene (Mao et al., 2014) and the BrLAX, BrPIN, BrPGP genes (Gao et al., 2017), which exhibited higher expression levels in the leaf apical tip and at the rosette stage than in the leaf base and at the heading stage in Chinese cabbage. However, the BcpLH (Yu et al., 2000), BrGRF (Wang et al., 2014a), and BrpSPL9 (Wang et al., 2014b) genes showed a gradual increase in transcript abundance from the rosette stage to the heading stage, implying different functions.

Gene knockout/knockdown methods such as the CRISPR-Cas9 system, RNA interference, T-DNA insertion, transposons, and some chemical reagent-induced mutations suffer from the limitations of off-target effects, functional redundancy, embryonic lethality, and multi-insertions (Pflieger et al., 2008). To the contrary, VIGS is a rapid and cost-effective RNA-mediated reverse genetics technology, and does not rely on the acquisition of mutants or transgenic plants. It can be easily applied in many species to study gene function either individually or on a large scale. pTY-S is a newly developed VIGS vector, which contains the entire genome of TYMV, and the 35S CAMV promoter and terminator (Pflieger et al., 2008; Huang et al., 2018a,b; Muntha et al., 2018). It possesses a strong infection ability, produces liable and robust VIGS effects that lasts throughout the plant life (Pflieger et al., 2008), and has been successfully applied in Chinese cabbage (Yu et al., 2018). In the present study, the expression of BrAN3 and BrBRM genes had been effectively knocked down by each silencing vector, and the formation of leafy heads had been observed in the inoculated Chinese cabbage plants at the early stage (**Figure 3A**). The results demonstrated the function of the BrAN3 and BrBRM genes in leafy head formation in Chinese

cabbage. Since the efficiency of VIGS is highly dependent on the virus-host affinity, and the infected TYMV eventually took over and killed the host plants at the later stage of head formation, it should be noted that the effect of knockdown of the BrAN3 and BrBRM genes on leafy head formation could be better examined in stable transgenic China cabbage in the near future.

Hormones such as auxin and gibberellins play important roles in the adaxial-abaxial patterning during leaf development in Chinese cabbage (He et al., 2000; Yu et al., 2000; Wang et al., 2014a; Gao et al., 2017). The identification of 41 DEGs from the plant hormone signal transduction term in the BrAN3-silencing leaf by RNA-Seq further demonstrated that plant hormones are key players in leafy head formation. In the present study, we found that cytokinin plays an important role in leafy head formation in Chinese cabbage. The significant increase and decrease in the expression of the BrARR4 and BrCRF2 genes in the BrAN3-silencing plants (**Figure 3B**) indicate that BrAN3 is involved in the cytokinin signaling pathway since both BrARR4 and BrCRF2 genes are cytokinin response regulators (Rashotte et al., 2006; Werner and Schmuelling, 2009; Holst et al., 2011; Vercruyssen et al., 2014).

Moreover, our RNA-Seq results revealed that the BrAN3 protein functions to induce the GA, BR, and SA signaling pathways, which inhibit leafy head formation in Chinese cabbage (**Figure 6**). The BrAN3 also inhibit the ethylene and JA signaling pathways, which promote leafy head formation in Chinese cabbage (**Figure 6**). The BrAN3 protein increases and decreases DEGs involved in the auxin, ABA, and CK signaling pathways, indicating the three pathways play both a positive and a negative roles in leafy head formation in Chinese cabbage (**Figure 6**). The accumulation and uneven distribution of auxin in the head leaves play an important role in Chinese cabbage leafy head formation (He et al., 2000; Gao et al., 2017). This is in accordance with our findings that a total of 20 auxin DEGs were significantly expressed in the BrAN3 silencing plants, including 13 up-regulated and 6 down-regulated DEGs. Among them, 3 DEGs (Bra017136, Bra005136, and Bra 023958, which correspond to BrLAX5, BrLAX6, and BrIAA29) had been demonstrated to be involved in leafy head formation in Chinese cabbage (Wang et al., 2012; Gao et al., 2017). Auxin mediates plant growth by affecting DELLA proteins (Fu and Harberd, 2003), while SWI3C, a core component of Arabidopsis SWI/SNF chromatin-remodeling complexes, also interacts with several DELLA proteins (Sarnowska et al., 2013). DELLA was reported as a repressor of ethylene signaling (Sarnowska et al., 2013), and interacts with JAZ in response of JA signaling (Lor and Olszewski, 2015; Vleesschauwer et al., 2016). Ethylene has been shown to induce petiole unequal growth between the adaxial and abaxial axes in Arabidopsis (Polko et al., 2012), and JAZ represses the JA signaling (Lor and Olszewski, 2015; Vleesschauwer et al., 2016). Accordingly, we found that the ethylene and JA signaling pathways were significantly induced in the BrAN3-silencing plants. The only DEG (Bra000421) which corresponds to the BrJAR1 gene was identified to be a DEG in leafy head formation in Chinese cabbage by transcriptomics and proteomics analyses (Zhang C.W. et al., 2016). In addition, GA signaling was proved to be a positive regulator of the SWI/SNF chromatin remodeling complexes (Rafal et al., 2013; Jose, 2014). BRM regulates primary root elongation by repressing the activity of the ABA signaling pathway (Han et al., 2012). Our findings that the BrAN3-silencing plants inhibited and promoted GA and ABA genes expression, respectively, are consistent with those reports (Han et al., 2012; Rafal et al., 2013; Jose, 2014). Moreover, BRM binds to the promoter of ARR16 for boundary

establishment in Arabidopsis cotyledon, which is an inhibitor of CK responses (Efroni et al., 2013). Similarly, we found that 4 and 1 DEGs involved in the CK signaling pathway were up-regulated and down-regulated in the BrAN3-silencing plants, respectively (**Figure 6**). The potential roles of these hormone signaling pathways in leafy head formation in Chinese cabbage were summarized in **Figure 8**. Further studies should be conducted to investigate the involvement of each hormone signaling pathway in leafy head formation in Chinese cabbage.

#### CONCLUSION

We identified and characterized the function of the BrAN3 gene in leafy head formation in Chinese cabbage. We also found that DEGs involved in auxin, ethylene, GA, JA, ABA, BR, CK, and SA signaling pathways play either positive or negative roles in leafy head formation in Chinese cabbage. This information could be useful for genetic engineering, gene editing, and molecular breeding of Chinese cabbage to improve heading and thus yield.

#### AUTHOR CONTRIBUTIONS

CZ and XH conceived the project. JY and LG conducted the experiments. JY, LG, WL, CZ, DX, XH, JY, and

#### REFERENCES


TL performed the data analysis. JY, LG, WL, CZ, XH, and LS wrote the manuscript. All authors read and approved the manuscript.

### FUNDING

This work was supported by the grants from the National Natural Science Foundation of China (31272172), the Fundamental Research Funds for the Central Universities (KYZ201826), and the Nature Science Foundation of Jiangsu Province (BK20141364).

#### ACKNOWLEDGMENTS

We thank Antoine Bouteilly (Centre National dela Recherche Scientifique) for providing the plasmid pTY-S.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2019.00520/ full#supplementary-material

with heading type in Chinese cabbage. Front. Genet. 8:176. doi: 10.3389/fgene. 2017.00176


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Yu, Gao, Liu, Song, Xiao, Liu, Hou and Zhang. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

fpls-10-00520 April 26, 2019 Time: 15:38 # 12

fpls-10-00112 February 9, 2019 Time: 17:6 # 1

# Characterization of Non-heading Mutation in Heading Chinese Cabbage (Brassica rapa L. ssp. pekinensis)

Jingrui Li<sup>1</sup>† , Xiaomeng Zhang<sup>1</sup>† , Yin Lu<sup>1</sup> , Dongxiao Feng<sup>1</sup> , Aixia Gu<sup>1</sup> , Shan Wang<sup>1</sup> , Fang Wu<sup>1</sup> , Xiangjie Su<sup>1</sup> , Xueping Chen<sup>1</sup> , Xing Li<sup>1</sup> , Mengyang Liu<sup>1</sup> , Shuangxi Fan<sup>2</sup> , Daling Feng<sup>1</sup> , Shuangxia Luo<sup>1</sup> , Shuxin Xuan<sup>1</sup> , Yanhua Wang<sup>1</sup> \*, Shuxing Shen<sup>1</sup> \* and Jianjun Zhao<sup>1</sup> \*

<sup>1</sup> Key Laboratory of Vegetable Germplasm Innovation and Utilization of Hebei, Collaborative Innovation Center of Vegetable Industry in Hebei, College of Horticulture, Hebei Agricultural University, Baoding, China, <sup>2</sup> Plant Science and Technology College, Beijing University of Agriculture, Beijing, China

#### Edited by:

Yuke He, Shanghai Institutes for Biological Sciences (CAS), China

#### Reviewed by:

Jun Yan, East China Normal University, China Lugang Zhang, Northwest A&F University, China

#### \*Correspondence:

Yanhua Wang yywyh@hebau.edu.cn Shuxing Shen shensx@hebau.edu.cn Jianjun Zhao jjz1971@aliyun.com †These authors have contributed

Specialty section:

equally to this work

This article was submitted to Plant Breeding, a section of the journal Frontiers in Plant Science

Received: 29 June 2018 Accepted: 23 January 2019 Published: 12 February 2019

#### Citation:

Li J, Zhang X, Lu Y, Feng D, Gu A, Wang S, Wu F, Su X, Chen X, Li X, Liu M, Fan S, Feng D, Luo S, Xuan S, Wang Y, Shen S and Zhao J (2019) Characterization of Non-heading Mutation in Heading Chinese Cabbage (Brassica rapa L. ssp. pekinensis). Front. Plant Sci. 10:112. doi: 10.3389/fpls.2019.00112 Heading is a key agronomic trait of Chinese cabbage. A non-heading mutant with flat growth of heading leaves (fg-1) was isolated from an EMS-induced mutant population of the heading Chinese cabbage inbred line A03. In fg-1 mutant plants, the heading leaves are flat similar to rosette leaves. The epidermal cells on the adaxial surface of these leaves are significantly smaller, while those on the abaxial surface are much larger than in A03 plants. The segregation of the heading phenotype in the F<sup>2</sup> and BC<sup>1</sup> population suggests that the mutant trait is controlled by a pair of recessive alleles. Phytohormone analysis at the early heading stage showed significant decreases in IAA, ABA, JA and SA, with increases in methyl IAA and trans-Zeatin levels, suggesting they may coordinate leaf adaxial-abaxial polarity, development and morphology in fg-1. RNA-sequencing analysis at the early heading stage showed a decrease in expression levels of several auxin transport (BrAUX1, BrLAXs, and BrPINs) and responsive genes. Transcript levels of important ABA responsive genes, including BrABF3, were up-regulated in mid-leaf sections suggesting that both auxin and ABA signaling pathways play important roles in regulating leaf heading. In addition, a significant reduction in BrIAMT1 transcripts in fg-1 might contribute to leaf epinastic growth. The expression profiles of 19 genes with known roles in leaf polarity were significantly different in fg-1 leaves compared to wild type, suggesting that these genes might also regulate leaf heading in Chinese cabbage. In conclusion, leaf heading in Chinese cabbage is controlled through a complex network of hormone signaling and abaxial-adaxial patterning pathways. These findings increase our understanding of the molecular basis of head formation in Chinese cabbage.

Keywords: Brassica rapa, mutant, genetic analysis, epidermis cell, RNA-Seq, phytohormones

# INTRODUCTION

Chinese cabbage (Brassica rapa ssp. pekinensis) is an important vegetable crop of the Brassica genus containing several species that are of agricultural and horticultural importance. Breeding has transformed the head morphology of this crop from a loose heading to semi-heading and finally a heading type. As the edible organ, the head of Chinese cabbage is the basis for its economic

**18**

fpls-10-00112 February 9, 2019 Time: 17:6 # 2

value. Curling, crinkling and folding of leaves are typical characteristics of heading in Chinese cabbage. The timing and compactness of head formation are affected by the time and degree of inward curling of the leaves. Leaf polarity and phytohormones (especially auxin) are critical for leaf architecture (Liu et al., 2011), but the exact mechanism of leaf folding in Chinese cabbage is still unclear.

Leaf polarity is composed of centro-lateral axis, proximaldistal axis and abaxial-adaxial polarity (Kim and Cho, 2006). The imbalance of abaxial-adaxial polarity is important for head formation (Mao et al., 2014). Many genes involved in abaxial-adaxial polarity have been cloned in Arabidopsis, providing useful insight for exploring head development in Chinese cabbage. The HD-ZIPIII family genes PHABULOSA (PHB), REVOLUTA (REV) (Otsuga et al., 2001), PHAVOLUTA (PHV) (McConnell et al., 2001), and HOMEOBOX8 (HB8), the transcription factors ASYMMETRICLEAVES 1 (AS1) and ASYMMETRICLEAVES 2 (AS2), and ta-siRNAs contribute to adaxial polarity. Auxin response factors (ARF3/ETT, ARF4), the KANAD gene family (KAN1, KAN2, KAN3, YABBY gene family) (Eshed et al., 2001; Kerstetter et al., 2001), and miRNA165/166 contribute to abaxial polarity (Palatnik et al., 2003; Hunter et al., 2006; Kidner and Timmermans, 2010; Townsley and Sinha, 2012). Although there is no heading in Arabidopsis, many genes related to abaxial-adaxial polarity in Arabidopsis also contribute to head formation in B. rapa (Liang et al., 2016). The re-sequencing data of different B. rapa and B. oleracea morphotypes were analyzed to detect signals of artificial selection that have shaped the complex heading trait by comparing genomic variation between heading and non-heading groups (Cheng et al., 2016). Many selection signals, or selective sweeps, including 15 loci that are under selection at syntenic positions in heading Chinese cabbages and cabbages, were detected in these two species. Several genes involved in the abaxial-adaxial patterning and leaf curvature were selected, such as BrARF3.1, BrARF4.1, BrKAN2.1, and BrKAN2.3 in B. rapa, and BoATHB15.2 (belonging to the HD-ZIPIIIs) and BoKAN2.2 in B. oleracea. However, the inheritance of head formation is complicated and the synergistic regulation of key genes in head development is still unclear.

Studies have shown that the synthesis, transport and signaling of phytohormones, especially auxin, play an important role in head formation in Chinese cabbage. He et al. (2000) reported that auxin participates in regulating head formation. Combined with the genome analysis of the convergence of B. rapa and B. oleracea (Cheng et al., 2016), gene enrichment analysis identified gibberellic acid (GA) biosynthesis and auxin- , cytokinin (CK)- and jasmonic acid (JA)-mediated signaling pathways. These pathways are known to be involved in leaf initiation and morphogenesis. Gao et al. (2017) found that the polar transport and uneven distribution of auxin affects head formation in Chinese cabbage. The auxin transport genes BrLAX (LIKE AUXIN RESISTANT), BrPIN (PIN-FORMED) and BrPGP (P-GLYCOPROTEIN) may also regulate the head development (Gao et al., 2017). In our previous study, the candidate genes BrGH3.12 and BrABF1 were identified using a Chinese cabbage-cabbage monosomic alien addition line AC<sup>4</sup> by RNA-seq analysis (Gu et al., 2017). Although these phytohormone-related genes have been associated with head formation, how they communicate together to regulate this process is largely unknown. In Arabidopsis, methyl IAA ester (MeIAA) contributes to leaf curvature (Pérez-Pérez et al., 2010), while there are limited reports about how MeIAA affects the head morphology in B. rapa.

The reference genome of Chinese cabbage was successfully completed in 2011 (Wang et al., 2011). As a result, it revealed a whole genome triplication (WGT) event since diverging from Arabidopsis that likely facilitated the generation of extensive diversity in morphotypes (Wang et al., 2011; Cheng et al., 2014, 2016). The completion of this work has greatly promoted the study of related traits in Chinese cabbage and laid the foundation for accelerating the molecular breeding of Chinese cabbage vegetables. The genome has greatly improved our abilities to characterize mutants for gene discovery and functional research. Ethyl methanesulfonate (EMS) is the most widely used reagent for mutagenesis that provides a high mutation frequency and relatively few chromosomal aberrations. EMS mutagenesis can be used to improve specific traits and is widely used in crop germplasm resource innovation. EMS has proven to be very successful in uncovering key regulatory genes contributing to a wide range of traits in Arabidopsis (Gabrielson et al., 2006; Chiu et al., 2007). In B. napus and other Brassica crops, mutant libraries in various cultivars have been constructed by EMS mutagenesis in order to study a range of variant trait-related genes (Stephenson et al., 2010; Wang N. et al., 2010). However, in Chinese cabbage, EMS mutants are rarely used as a genetic analysis for candidate genes.

A mutant library containing 4253 M<sup>1</sup> lines and the resulting M<sup>2</sup> population was constructed by artificial EMS mutagenesis of the Chinese cabbage inbred line A03 (Lu et al., 2016). One flat growth non-heading mutant, fg-1, was obtained from the EMS-induced mutagenesis population that has flat leaves prior to the heading stage with a wrinkled leaf surface compared to the wild type. Using the mutant fg-1 and its wild type A03, we revealed the genetic structure of the mutant heading trait in Chinese cabbage by creating segregating populations. Combining the RNA-seq and phytohormone quantifications, the molecular regulatory mechanism of head development was investigated by assessing transcript level changes and characterizing leaf growth, phytohormone levels and leaf epidermal cell morphology. In addition, a possible regulatory model is proposed. The purpose of this study was to identify new genes regulating head development in Chinese cabbage and generating new genetic resources for future Chinese cabbage crop improvement studies.

# MATERIALS AND METHODS

#### Plant Materials

A mutant library of Chinese cabbage was developed by EMS treatment of seeds from the inbred line A03 (Lu et al., 2016), from which a non-heading mutant of the M<sup>5</sup> generation with flat growth of heading leaves (fg-1) was isolated. Wild type A03 has an outward-curling heading pattern on the top. The mutant fg-1 has flat leaves during growth before the heading stage and trends to heading at the heading stage. The populations of F1, F<sup>2</sup> and two BC<sup>1</sup> (F<sup>1</sup> × fg-1 and F<sup>1</sup> × wild type) were developed from the cross between A03 and fg-1, which were used as the experimental materials for genetic analysis of the mutant trait. The plants were grown in a plastic tunnel on the experimental farm at Hebei Agricultural University in Baoding (115.47 E, 38.87 N), China, in 2016 and 2017.

In August 2017, 60 plants each of A03 and fg-1 of the M<sup>6</sup> generation were grown in the same plastic tunnel at Hebei Agricultural University. At the early heading stage (80 days after sowing), the 16th leaf from the exterior of the developing head was sampled at four positions: apical, middle, bottom of the soft leaf and basal of the whole leaf, which were named a, b, c and d, respectively (**Supplementary Figure S1**). Three biological replicates were used for further analysis. All leaf samples were flash frozen in liquid nitrogen and stored at −80◦C for RNA or phytohormone analysis.

#### Inheritance of the Mutant Trait

The mutant fg-1, wild type A03, five F<sup>1</sup> lines, 163 F<sup>2</sup> lines, 14 F<sup>1</sup> × fg-1 BC<sup>1</sup> lines and 10 F<sup>1</sup> × wild type BC<sup>1</sup> lines were planted in a plastic tunnel in July 2016. To confirm the results, 145 F<sup>2</sup> lines and their parents were planted again in a plastic tunnel in July 2017. The number of plants with flat growth and normal heading were counted and a Chi square test was performed. The height and expansion of plants were investigated for A03 and fg-1.

#### Morphological Characteristics of fg-1 at Different Developmental Stages

During the late rosette (50 days after sowing), early heading (80 days after sowing) and heading stages (90 days after sowing), the angle was measured in different layers. During the late rosette and heading stages, the leaf length, width and area of the 13th and 19th leaf from the exterior was measured via ImageJ. At the heading stage, the expansion degree (the maximum distance of the outside leaves) and the plant height were measured for fg-1 and the wild type.

#### Analysis of Leaf Abaxial and Adaxial Epidermal Cell Structure and Area

During the rosette (40 days after sowing) and heading stages, 2 mm × 2 mm from the top and the middle edge of leaves were prepared for abaxial epidermal cell structure and area imaging by scanning electron microscopy. Leaves were fixed in 2.5% amyl glycol at room temperature for 24 h then rinsed completely with 0.1 M phosphate buffer and fixed with 1% osmium acid for 2 h. The samples were dehydrated with an alcohol gradient (30, 50, 70, 80, 90, 95, and 100%) followed by an overnight treatment with isoamyl acetate. Samples were dried using a LEICA EM CPD 030 critical point drying. An Eiko IB5 Ion Coater instrument was used to spray gold after the sample was glued to the table. A Hitachi SU-8010 scanning electron microscope was used to observe and photograph the samples.

During the late heading stage (100 days after sowing), the adaxial epidermis cell structure at the top and the middle edge of the leaves were observed using the glue-marking method. The glue was evenly coated on the leaves. After air-drying, the cells were removed by tweezers and spread on the slide for imaging. Cell structure was examined using the optical microscope (Olympus CX41, Tokyo, Japan). The area of the cells was measured using ImageJ software (Collins, 2007).

#### RNA Extraction, cDNA Library Preparation and Sequencing for RNA-Seq

Total RNA was isolated using an RNA extraction kit (TIANGEN, China). RNA purity and concentration were assessed as previously described (Zhao et al., 2014). The cDNA library was prepared and sequenced as previously described (Zhao et al., 2014; Gu et al., 2017).

#### RNA-Seq Data Analysis

Raw reads were filtered by removing adaptors, poly(N) and low quality reads (Li et al., 2015). The reference genome and gene model annotation files were downloaded from the B. rapa database (<sup>1</sup> v1.5). Filtered reads were mapped to the genome using HISAT2 v2.0.5 (Kim et al., 2015). HTSeq v0.6.0 (Anders, 2010) was used to generate counts per gene. The FPKM values of each gene were calculated based on the gene length and read counts. Differential expression analysis of each set of two groups (three biological replicates per group) was performed using the DESeq R package (1.16.1.) (Anders and Huber, 2010, 2012; Wang L. et al., 2010). The P-values were adjusted using the Benjamini and Hochberg (1995). method Genes with a corrected P < 0.05 were considered differentially expressed.

# Validation of RNA-Seq by qRT-PCR

Quantitative real-time PCR (qRT-PCR) was carried out as previously described (Zhao et al., 2014; Gu et al., 2017; Zhang et al., 2017). To verify the RNA-seq data, we analyzed the expression of ten genes by qRT-PCR. A Chinese cabbage actin gene (Bra009081) was used as an internal reference. The 2 −[Ct−Ct(actin)] × 1000 method was used to calculate gene relative expression levels between wild type and fg-1. Primers are listed in **Supplementary Table S1**. Each experiment was repeated three times.

#### Concentration Analysis of Different Phytohormones in fg-1 and Wild Type

The b sections of (**Supplementary Figure S1**) the 16th leaf from the exterior at the early heading stage were harvested, weighed and immediately frozen in liquid nitrogen. Samples (120 mg) were ground into a fine powder and extracted with methanol/water (v/v = 8/2) at 4◦C. The concentrations of six classes of phytohormones [auxin, abscisic acid (ABA), GA, salicylic acid (SA), JA and CK], including 25 forms were analyzed using a high-pressure liquid chromatography-electrospray tandem mass spectrometry (LC-ESI-MS/MS) system as

fpls-10-00112 February 9, 2019 Time: 17:6 # 3

<sup>1</sup>http://brassicadb.org/brad/index.php

previously described (Wang et al., 2018). Three technical replicates and three biological replicates were conducted.

#### RESULTS

#### Inheritance of Mutant Trait

fpls-10-00112 February 9, 2019 Time: 17:6 # 4

Genetic analysis revealed that the phenotypic traits of all F<sup>1</sup> plants were similar to the wild type. In the F<sup>2</sup> populations of years 2016 and 2017, the proportions of wild type and mutant plants were 3.17:1 and 3.68:1, respectively, conforming to the 3:1 ratio (Chi square test: χ <sup>2</sup> = 0.691 and 0.311). In the backcross population F<sup>1</sup> × fg-1, the proportion of wild type and mutant plants was 1.33:1, conforming to the 1:1 ratio (Chi square test: χ <sup>2</sup> = 0.593), and all plants in the backcross group F<sup>1</sup> × wild type were similar to wild type (**Table 1**). This phenotypic analysis suggests that the flat growth trait is controlled by a pair of recessive alleles.

#### Morphological Characteristics of fg-1 During Leaf Development

During whole leaf development, all fg-1 expanded leaves show flat growth with a reduced leaf angle to the ground (**Supplementary Table S2**). All fg-1 leaves are visibly wrinkled and lacking the outward-curling features seen in wild type. In the wild type, the angles from rosette leaves to the ground are larger (**Supplementary Table S2**). Additionally, the rosette leaves are slightly wrinkled; the heading leaves are outward-curling at the top edge. From the early heading stage, the newly expanded leaves curl transversely and longitudinally to form the head (**Figure 1**).

During the late rosette and heading stages, the leaf length and area in fg-1 are significantly smaller than that in wild type, while there is no significant difference in leaf width (**Supplementary Table S3**). The leaf shape in fg-1 tends to be obovate, but that in wild type is close to rectangle.

At the heading stage, the expansion degree of the mutant plants (the maximum distance between the outside leaves) was similar to that of the wild type, while plant height was significantly less (**Supplementary Figure S2**).

#### Leaf Abaxial and Adaxial Epidermal Cell Structure and Area in fg-1

To further characterize the difference between the leaf morphology of wild type and fg-1, abaxial epidermal cells

TABLE 1 | Segregation ratios of F1, F<sup>2</sup> and two BC1F1s between fg-1 and wild type.


of leaf sections were observed by scanning electron microscopy at the rosette and heading stage (**Figure 2**). Adaxial epidermal cells of leaf sections were also observed at the late heading stage (**Supplementary Figure S3**). At both stages, more large and slender cells in the abaxial epidermis were observed at the top and central edge of fg-1 mutant leaves compared to wild type. The abaxial epidermis cell areas of the top and the central leaf edges in fg-1 were significantly greater than that in the wild type at the rosette and heading stages (**Table 2**). At the rosette stage, abaxial epidermis cell area of both the top and the central edge of the leaf in fg-1 was 57.7 and 28.9% larger than that in the wild type, respectively. At the heading stage, the abaxial epidermis cell area of both the top and central edge was 35.8 and 39.1% larger than that of the wild type, respectively (**Table 2**). At the late heading stage, the adaxial epidermis cell area of both the top and central edge was 26.1 and 41.1% lower than that of the wild type, respectively (**Table 3**).

#### Identification of Differentially Expressed Genes via RNA-seq Between fg-1 and A03 Leaves

Following read processing and quality filtering, roughly 6.34 Gb remained for each sample with over 85% of reads uniquely mapped (**Supplementary Table S4**). We identified 7669 DEGs in the top leaf section, with 3966 down-regulated and 3703 up-regulated; 5022 DEGs were found in the middle section, with 2625 down-regulated and 2397 up-regulated; 3907 DEGs were found in the base of the soft leaf section, with 1900 down-regulated and 2007 up-regulated; and 4823 DEGs were found in the basal section, with 2090 down-regulated and 2733 up-regulated (**Supplementary Table S5**). To validate the RNA-seq data, 10 representative DEGs were assessed by qRT-PCR analysis. The results indicated that the expression of the

TABLE 2 | Area of leaf abaxial epidermal cells.


Each value is the mean ± SE for 60 cells from three independent replicates. Different letters indicate significant differences between treatments (P < 0.05) according to Duncan's multiple range test.

TABLE 3 | Area of leafy adaxial epidermal cells in the late heading stage.


Each value is the mean ± SE for 60 cells from three biological replicates. Different letters indicate significant differences between treatments (P < 0.05) according to Duncan's multiple range test.

fpls-10-00112 February 9, 2019 Time: 17:6 # 5

FIGURE 2 | Scanning electron microscopy for abaxial epidermis cells of leaf sections at different developmental stages for the wild type and fg-1. (A,B) represent abaxial epidermis cells of the top and central leaf edge at the rosette stage in wild type. (C,D) represent abaxial epidermis cells of the top and central leaf edge at the heading stage in wild type. (E,F) represent abaxial epidermis cells of the top and central leaf edge at the rosette stage in fg-1. (G,H) represent abaxial epidermis cells of the top and central leaf edge at the heading stage in fg-1. Slender cells are indicated by arrows. Length of black bars is 10 µm.

genes between the two methods were consistent (R <sup>2</sup> = 0.8441) (**Supplementary Figure S4**).

#### Pathway Analysis Involved in Head Formation in Chinese Cabbage

Pathway definitions were derived from the KEGG (Kyoto Encyclopedia of Genes and Genomes) database (**Figure 3**). Among the main pathways of DEGs, plant hormone signal transduction has a known role in regulating leaf development. Except for genes BrLAX2, BrPIN1, and BrPIN6 in leaf section d (**Figure 4A**, subcluster c), the expression levels of the auxin influx and efflux carrier genes, three BrAUX1 genes, four BrLAXs and six BrPINs were reduced (**Figure 4A** and **Supplementary Tables S5, S6**). Additionally, the main expression profiles of genes (including BrARF7, 19 BrAUX/IAAs, six BrGH3 and 50 BrSAURs) involved in IAA signaling were also down-regulated (**Supplementary Figure S6A**, subcluster b and c; **Supplementary Table S7**), but the expression levels of three syntenic genes of AtARF5 in leaf section d were elevated in fg-1 (**Supplementary Figure S6A**, subcluster a; **Supplementary Table S7**). The expression patterns of 30 genes (including 10 BrPYR/PYLs, nine BrPP2Cs, six BrSnRK2s and five BrABFs) in abscisic acid (ABA) signaling (**Supplementary Figure S6B**, subcluster b; **Supplementary Table S7**), 19 genes (including three BrHK3 genes, three AHPs, two BrB-ARRs and 11 BrA-ARRs) in CK signaling (**Supplementary Figure S7A**), 16 genes (including three BrJAR1 syntenic genes, two BrCOI1 syntenic genes, 10 BrJAZs and one BrMYC2) in JA signal transduction (**Supplementary Figure S7B**) and 12 genes (nine BrTGAs and three BrNRRs) in SA signal transduction (**Supplementary Figure S7C**) were altered.

However, the expression levels of BrIAMT1 (Bra002919) encoding the IAA carboxyl methyltransferase 1, known to alter leaf curvature phenotypes through an auxin-regulated developmental process (Qin et al., 2005), decreased in three sections (a, b, c) of fg-1 leaves (**Figure 4B**, subcluster d; **Supplementary Figure S5**). In addition, except for three BrT in subcluster h and BrTCP22 (Bra008001) in subcluster d fpls-10-00112 February 9, 2019 Time: 17:6 # 6

pathways of DEGs in sections a–d of the fg-1 leaf compared with the wild type, respectively.

(**Figure 4B** and **Supplementary Figure S5**), increased expression levels of the remaining nine BrTCP genes that are known to regulate leaf curvature were observed in the fg-1 a, b and c leaf sections (**Figure 4B** and **Supplementary Table S6**). Additionally, 19 genes regulating leaf abaxial-adaxial patterning, including BrKAN1, BrBOP2 (BrNPR5), six BrHD-ZIPIIIs (BrREV, BrHB8.1, BrHB8.2, BrHB9, BrHB14.1, BrHB14.2), four BrYABs (BrYAB1.1, BrYAB1.2, BrYAB2, BrYAB3) and seven BrKNOX class II genes (BrKNAT1, BrKNAT2, BrKNAT3, BrKNAT4.1, BrKNAT4.2, BrKNAT5, BrKNAT6), were differentially expressed in several sections of fg-1 leaves (**Supplementary Figure S8** and **Supplementary Table S9**).

# Phytohormone Quantification in fg-1 and Wild Type

Concentrations of auxin, including IAA and MeIAA, ABA, salicylic acid (SA), JA and CK, including trans-Zeatin (tZ), in fg-1 were significantly different from wild type. Specifically, in fg-1, levels of IAA, ABA, JA and SA were reduced 25.2, 51.3, 55.7, and 28.7% compared with A03, respectively (**Figures 5A,C–E**). However, levels of MeIAA and tZ were increased in fg-1 to 3.2 and 1.7 fold that of the wild type, respectively (**Figures 5B,F**). The levels of other phytohormones such as GA showed no difference between fg-1 and wild type.

# DISCUSSION

To explore the heading mechanism in Chinese cabbage, a comprehensive regulatory network was constructed based on the above results (**Figure 6**).

# Epidermal Cell Shape Variation Between fg-1 and Wild Type

Leaf shape, such as curvature and wrinkling, is closely related to cell structure (Liu et al., 2011). Leaves also have unique abaxialadaxial characteristics, including the distribution of hair and stoma on the surface of the leaf and the type and arrangement of the cells on the abaxial side of the leaf (Tsukaya, 2002). The Arabidopsis hyl1 mutant lacks the long and narrow cells in the abaxial epidermis, resulting in dramatically inward curling leaves (Liu et al., 2011). Interestingly, similar to the inward curling phenotype in Arabidopsis hyl1, leaves of wild-type Chinese cabbage also curl inward and form a head due to the lack of these long and narrow cell types in the abaxial epidermis. The presence of elongated and narrow cells in the abaxial epidermis leads to the non-heading flat growth of fg-1 leaves (**Figure 2** and **Table 1**). The single adaxial cell area in wild type is significantly larger than in fg-1 (**Table 3** and **Supplementary Figure S3**), suggesting that cell shape and size play an important role in

fpls-10-00112 February 9, 2019 Time: 17:6 # 7

maintaining leaf polarity and regulating heading in Chinese cabbage (**Figure 6**).

# The Role of Phytohormones in Leaf Development

Phytohormones play essential roles in regulating leaf shape and polarity (Gazzarrini and Mccourt, 2003; Santner and Estelle, 2009; Stamm and Kumar, 2010; Cheng et al., 2016). Four phytohormones in particular (CK, auxin, GA and JA) have been implicated in heading morphotypes of both Chinese cabbage and cabbage (Cheng et al., 2016). Auxin directly influences cell growth (Teale et al., 2006) but also contributes to axial polarity and patterning of leaves (Barkoulas et al., 2007). IAA is the predominant form of auxin in plants (Teale et al., 2006), while MeIAA, a non-polar and more potent auxin based on hypocotyl elongation assays, plays an important role in regulating auxin homeostasis and plant development (Qin et al., 2005). In fg-1, IAA homeostasis was altered compared to wild type with lower levels of IAA and elevated levels of MeIAA (**Figures 5A,B**). These differences might contribute to the disrupted axial polarity (**Figure 1**) and cell enlargement observed in fg-1 leaves (**Figure 2**). IAA transport and signaling play essential roles in controlling plant growth and development (Teale et al., 2006). Consistent with the reduced accumulation of IAA, the expression levels of auxin influx and efflux carrier genes were reduced, similar to the majority of auxin responsive genes (**Figure 4A** and **Supplementary Figure S6A**), providing further support for the importance of auxin in head formation.

Endogenous ABA is involved in plant responses to environmental stresses in addition to fine-tuning plant development through a regulatory circuit of primary metabolism, cell growth and cell division (Fedoroff, 2002; Himmelbach et al., 2003). Cross-talk between ABA and auxin signaling occurs via AXR2/IAA7 (Timpte et al., 1994; Nagpal et al., 2000; Fedoroff, 2002). The homeostasis of ABA is also disturbed in fg-1 (**Figure 5C**), along with changes in ABA signaling (**Supplementary Figure S6B**). The expression of two AtIAA7 syntenic genes in subclusters b and c were significantly down-regulated in fg-1 (**Supplementary Figure S6A**) suggesting that the cross-talk between ABA and IAA signaling impacts leaf heading (**Figure 6**).

CKs have profound roles in plant growth regulation, such as release of lateral buds from apical dominance and delay of senescence (Kasahara et al., 2004). As a form of CKs, tZ is active and regulates plant development (Inoue et al., 2001; Yamada et al., 2001). In this study, tZ accumulates (**Figure 5F**) in fg-1 suggesting that higher concentrations of tZ might delay leaf senescence via CK signaling resulting in the non-heading phenotype. At the late heading stage of fg-1 plants, leaves begin heading, but under suboptimal conditions (e.g., low temperature and short day length), the leaves do not form a complete head (**Figure 1**).

SA has crucial roles in regulating physiological processes via a complex SA signaling network (Vanacker et al., 2001; Vicente and Plasencia, 2011). Altered SA levels disturb normal growth phenotypes (e.g., plant stature, leaf morphology, and cell size) (Greenberg et al., 2000; Song et al., 2004; Brodersen et al., 2005). JA regulates plant stress response, growth and development (e.g., leaf movement and senescence) through a regulatory network that integrates other plant hormones including SA, IAA and ABA (Nakamura et al., 2006; Wasternack, 2007; Wasternack and Strnad, 2016). In fg-1, the accumulation of both JA and SA was reduced, with

fpls-10-00112 February 9, 2019 Time: 17:6 # 8

fpls-10-00112 February 9, 2019 Time: 17:6 # 9

up-regulated expression patterns of several genes involved in JA and SA signaling (**Figures 5D,E** and **Supplementary Table S8**). The abnormal leaf growth and development seen in fg-1 could be due to disrupted coordination of multiple hormone signals.

In leaf section b of fg-1, the expression patterns of four BrPYR/PYLs, four BrPP2Cs, three BrSnRK2s and BrABF3 involved in ABA signaling were up-regulated, which was contrary to the change of ABA content. This same trend of gene expression patterns was observed in JA and SA signaling pathways in the fg-1 leaf section b. One explanation for this discrepancy is that these genes may be involved in multifunctional roles in other signaling pathways important for leaf development (Wasternack and Strnad, 2016).

#### Impacts of BrIAMT1 and BrTCPs on Leaf Polarity in Chinese Cabbage

The role of IAMT1 in Arabidopsis leaf development was implicated in an auxin-regulated developmental process (Qin et al., 2005). Decreased expression levels of IAMT1 causes dramatic epinastic leaf phenotypes (including smaller leaves and dwarfism) in IAMT1-RNAi Arabidopsis plants (Qin et al., 2005). Interestingly, the expression level of BrIAMT1 (Bra002919) was also down-regulated in fg-1 (**Figure 4A**, subcluster d) with a similar axial patterning of leaves and plant dwarfism as shown in IAMT1-RNAi Arabidopsis plants.

A subset of TCP genes (e.g., TCP2, TCP3, TCP4, and TCP10) have been shown to play important roles in regulating leaf morphology (Nath et al., 2003; Palatnik et al., 2003; Mao et al., 2014). Higher IAMT1 expression levels were correlated with lower expression of the TCP genes in Arabidopsis iamt1-D mutants that display dramatic hyponastic leaf phenotypes (Qin et al., 2005). In contrast, decreased mRNA levels of BrIAMT1 and increased expression of nine BrTCP genes (**Figure 4B**, subcluster e, f and g) may contribute to the epinastic growth of leaves in fg-1 (**Figure 6**).

It has been reported that the enzyme IAMT1 converts IAA to MeIAA in vitro (Qin et al., 2005; Zhao et al., 2008). Nevertheless, MeIAA is hard to detect (Qin et al., 2005), and it is unclear whether plants make MeIAA or whether IAMT1 could catalyze IAA to MeIAA in vivo (Zimmerman and Hitchcock, 1937; Yang et al., 2008). In this study, the expression levels of IAMT1 were decreased in leaf sections a, b and c of fg-1, while the concentration of MeIAA, detected via the HPLC-ESI-MS/MS system, was elevated in fg-1 (**Figure 5**). There are two possible explanations for these results. First, the IAMT1 enzyme function may be distinct

#### REFERENCES


between in vitro and in vivo conditions. Second, the gene underlying the fg-1 mutant may disrupt the BrIAMT1 catalysis in methylating IAA to MeIAA. The altered hormone levels and expression patterns of hormone responsive genes suggests a complex interplay of hormone signaling is necessary for proper head formation. Deciphering this complex network through additional genetic and genomic studies is needed to inform future breeding efforts for head morphology in Chinese cabbage.

#### AUTHOR CONTRIBUTIONS

JL and XZ performed the research and wrote the manuscript. DxF and XS surveyed the morphological characteristics. ML, AG, and SX analyzed the data of RNA-Seq. SL, SW, and YL performed the genetic analysis. DlF, XL, and FW performed the morphology analysis. XC and SF reviewed the manuscript. JZ, SS, and YW designed the research and reviewed the manuscript. All authors declare no competing financial interests.

#### FUNDING

The financial support for this work was provided by the National Key R&D Program of China (2017YFD0101802 and 2016YFD0100204-17), the National Natural Science Foundation of China (31672151 and 31801857), International Cooperation project in the Science and Technology Support Program of Hebei (17396315D), Hundred Innovative Talent Program of Hebei (SLRC2017040), the Natural Science Foundation of Hebei (C2016204170), the Science and Technology Support Program of Hebei (16226304D-2), the Science and Technology Research Project of Hebei Colleges and Universities (QN2017084).

#### ACKNOWLEDGMENTS

We thank Dr. Kathleen Greenham at Dartmouth College in Hanover, NH, for critical reading of the manuscript.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2019.00112/ full#supplementary-material


fpls-10-00112 February 9, 2019 Time: 17:6 # 10


family of Arabidopsis. Plant Physiol. 147, 1034–1045. doi: 10.1104/pp.108. 118224


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Li, Zhang, Lu, Feng, Gu, Wang, Wu, Su, Chen, Li, Liu, Fan, Feng, Luo, Xuan, Wang, Shen and Zhao. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

fpls-10-00112 February 9, 2019 Time: 17:6 # 11

# Plant Organ Shapes Are Regulated by Protein Interactions and Associations With Microtubules

Mark D. Lazzaro1,2, Shan Wu<sup>3</sup> , Ashley Snouffer<sup>2</sup> , Yanping Wang<sup>4</sup> and Esther van der Knaap2,5,6 \*

<sup>1</sup> Department of Biology, College of Charleston, Charleston, SC, United States, <sup>2</sup> Center for Applied Genetic Technologies, University of Georgia, Athens, GA, United States, <sup>3</sup> Boyce Thompson Institute, Cornell University, Ithaca, NY, United States, <sup>4</sup> National Engineering Research Center for Vegetables, Beijing Academy of Agriculture and Forestry Sciences, Beijing, China, <sup>5</sup> Institute for Plant Breeding, Genetics and Genomics, University of Georgia, Athens, GA, United States, <sup>6</sup> Department of Horticulture, University of Georgia, Athens, GA, United States

Plant organ shape is determined by the spatial-temporal expression of genes that control the direction and rate of cell division and expansion, as well as the mechanical constraints provided by the rigid cell walls and surrounding cells. Despite the importance of organ morphology during the plant life cycle, the interplay of patterning genes with these mechanical constraints and the cytoskeleton is poorly understood. Shapes of harvestable plant organs such as fruits, leaves, seeds and tubers vary dramatically among, and within crop plants. Years of selection have led to the accumulation of mutations in genes regulating organ shapes, allowing us to identify new genetic and molecular components controlling morphology as well as the interactions among the proteins. Using tomato as a model, we discuss the interaction of Ovate Family Proteins (OFPs) with a subset of TONNEAU1-recruiting motif family of proteins (TRMs) as a part of the protein network that appears to be required for interactions with the microtubules leading to coordinated multicellular growth in plants. In addition, SUN and other members of the IQD family also exert their effects on organ shape by interacting with microtubules. In this review, we aim to illuminate the probable mechanistic aspects of organ growth mediated by OFP-TRM and SUN/IQD via their interactions with the cytoskeleton.

#### Keywords: OFP, TRM, SUN, IQD, microtubules, organ shape

# INTRODUCTION

Plant organs display remarkable phenotypic diversity among and within species. Especially for cultivated crops, selection for the harvestable organs has led to greatly increased size and variable shapes of the produce. This diversity is critical for the successful marketing of a wide array of foods such as fruits, vegetables, seeds, leaves, and tubers. Recent studies have revealed many genes

#### Edited by:

Han Xiao, Shanghai Institutes for Biological Sciences (CAS), China

#### Reviewed by:

Li Zhengguo, Chongqing University, China Pengwei Wang, Huazhong Agricultural University, China

> \*Correspondence: Esther van der Knaap vanderkn@uga.edu

#### Specialty section:

This article was submitted to Plant Physiology, a section of the journal Frontiers in Plant Science

Received: 05 August 2018 Accepted: 14 November 2018 Published: 04 December 2018

#### Citation:

Lazzaro MD, Wu S, Snouffer A, Wang Y and van der Knaap E (2018) Plant Organ Shapes Are Regulated by Protein Interactions and Associations With Microtubules. Front. Plant Sci. 9:1766. doi: 10.3389/fpls.2018.01766

**29**

**Abbreviations:** CAP350, centrosome associated protein 350; cMT, cortical microtubules; CMU, cellulose microtubule uncoupling; CSC, cellulose synthase complex; FOP, fibroblast growth factor receptor 1 oncogene partner; IQD, isoleucine glutamine domain; KLCR, kinesin light chain related; MTOC, microtubule organizing center; NIL, nearly isogenic line; OFD1, orofaciodigital syndrome 1; OFP, ovate family protein; POK, phragmoplast-orienting kinesin; PPB, preprophase band; ROP, rho-GTPase of plants; SPR, spiral; TRM, tonneau1 recruiting motif; TTP, tonneau1-TRM-phosphatase 2C; γ-TuRC, tubulin ring complex.

that control the growth form of agriculturally important organs (Zuo and Li, 2014; van der Knaap and Ostergaard, 2017). This includes a newly discovered genetic pathway, which through protein interactions and associations with microtubules is proposed to lead to changes in cell division patterns that accompany the different growth forms (Wu et al., 2018). Mechanistically, how the growth forms are controlled by these genes is largely unknown.

The classification of varieties of the same crop based on morphological descriptors is critical. Organizations such as the Union for the Protection of New Varieties of Plants<sup>1</sup> and the International Plant Genetic Resources Institute (IPGRI)<sup>2</sup> developed descriptors of the shape of many vegetables and fruits such that varieties can be distinguished from one another. These descriptors have become the framework for the identification of genes underlying the morphological variation in crops like tomato (Brewer et al., 2006; Rodriguez et al., 2011a). Consumers recognize the shape and size of vegetables and fruits for the different culinary purposes and/or cultural significances (Pickersgill, 2007; Daunay et al., 2008; De Haan, 2009; Monforte et al., 2014). Similarly important for grains, the slender rice grain shape is associated with improved transparent appearance and reduced undesirable grain quality and is, therefore, highly sought-after in certain cuisines (Calingacion et al., 2014; Harberd, 2015).

#### PROCESSES THAT CONTROL ORGAN MORPHOLOGY

Lateral plant organs such as leaves and fruits typically initiate in the flanks of apical meristems. Together with the hormone auxin, AGAMOUS (for ovaries/fruits) and CUP SHAPED COTYLEDON/NO APICAL MERISTEM initiate organ primordia by setting up organ identity and boundaries (Maugarny-Cales and Laufs, 2018). To change from a meristematic cell fate to an organ fate, the down regulation of KNOXI transcription factors by ASYMMETRIC LEAVES1 and LATERAL ORGAN BOUNDARIES DOMAIN proteins is required (Maugarny-Cales and Laufs, 2018). Many hormones play important roles in the growth of organs, in particular gibberellins and brassinosteroids.

The patterns of further outgrowth occur along different axes: the proximal-distal, the medial-lateral and the abaxial-adaxial axis (Van der Knaap et al., 2014). Simply stated, isotropic growth along all three axes tends to lead to larger and round shapes whereas anisotropic growth leads to alternate shapes. For multidimensional organs such as the fruit, the different tissue types grow in an anisotropic way and together form an overall spherical or elongated shape (van der Knaap and Ostergaard, 2017). At the cellular level, the growth patterns are manifested by a combination of cell proliferation (growth and division) and cell enlargement (growth without cell division) driven by turgor pressure. The directions of cell enlargement are guided and restricted by cellulose microfibrils, which are glucose polymers bundled together by hydrogen bonds and Van der Waals forces. These polymers are deposited into the cell wall by CSCs guided by cMTs, In cells undergoing anisotropic expansion, cellulose microfibrils are deposited perpendicular to the axis of expansion and are coaligned with cMTs (Szymanski and Cosgrove, 2009; McFarlane et al., 2014). During cell proliferation, the plane of cell division is determined by the positioning of the preprophase band (PPB) (Rasmussen et al., 2011). The duration and rates of cell proliferation also affect the pattern of growth. Since plant cells are bound to surrounding cells by cell walls, once division has taken place, including formation of the phragmoplast, plant cells are positioned in the same relative location as when they were formed. Therefore, the orientation of cell division has a dramatic effect on the final shape of plant organs (**Figure 1**; Meyerowitz, 1997; Jenik and Irish, 2000; Van Damme et al., 2007; Schaefer et al., 2017). Mechanical forces provide direct signals leading to coordinated growth toward the final organ shape and size (von Wangenheim et al., 2016). During lateral organ initiation, a highly organized supracellular alignment of microtubule arrays forms along the maximal stress in the region between the meristematic dome and lateral primordia (Hamant et al., 2008). The microtubules guide the directional deposition of cellulose microfibrils, which reinforces the cell wall strength along the appropriate axes to separate the new organs and the undifferentiated cells. During growth, microtubule array dynamics are regulated to respond to the mechanical forces (Uyttewaal et al., 2010). The reorientation of microtubule arrays along the maximal tensile stress can control the directions of cell division and cell expansion leading to heterogeneous growth.

Understanding the molecular basis of shape of harvestable organs comes mostly from studies conducted in tomato and rice. The increase in rice grain size is often accompanied with altered shape, and found to be under the control of proteins involved in diverse pathways such as G-protein signaling, the ubiquitin–proteasome pathway, phytohormone signaling including brassinolides, auxin and cytokinin, as well as transcriptional regulation (Zuo and Li, 2014; Zheng et al., 2015). In the case of tomato, the identified proteins appear functionally less diverse as they seem to interact with the cytoskeleton. Specifically, a mutation in OVATE, the founding member of the OFP class, and another member named SlOFP20 result in a distinct pear shaped tomato fruit (Wu et al., 2018). OVATE and SlOFP20 interact with several members of the Tonneau1 Recruitment Motif (TRM) proteins, which are often found associated with microtubules (Hamant et al., 2008; Wu et al., 2018). SUN, a member of the IQ Domain (IQD) family, also impacts tomato fruit shape (Xiao et al., 2008). Members of the IQD family have been found to interact with calmodulin (CaM) as well as microtubule binding proteins Kinesin Light Chain-Related protein-1 (KLCR1) and SPR2 to regulate microtubule structure based on external auxin and calcium inputs (Burstenbinder et al., 2013, 2017a,b; Wendrich et al., 2018; Yang et al., 2018).

<sup>1</sup>http://www.upov.int/upovlex/en/upov\_convention.html

<sup>2</sup>https://www.bioversityinternational.org/fileadmin/bioversity/publications/Web\_ version/572/ch01.htm

FIGURE 2 | Effect of the fruit shape genes in the wild type tomato background. The loci were introgressed (sun, ovate, sov1) or edited (trm5) in the Solanum pimpinellifolium accession LA1589 background to create near isogenic lines (NILs). WT, wild type; sov1, suppressor of ovate corresponding to SlOFP20. The single natural NILs are shown on the top of the figure, whereas the double and triple NILs are shown below the single NILs.

### THE ROLE OF OFPs, SUNs AND TRMs ON ORGAN SHAPE

#### OFP and SUN

OVATE and SUN are two important genes controlling tomato fruit shape. The shape of many oval shaped varieties, including grape tomatoes, is controlled by OVATE. SUN can be found in very elongated, tapered or oxheart shaped heirloom as well as commercially grown plum tomatoes (Ku et al., 1999; Liu et al., 2003; Rodriguez et al., 2011b; Van der Knaap et al., 2014). OVATE is the founding member of the OFP class. Recently a new fruit shape gene was identified as a suppressor of ovate (sov1). This fruit shape gene is another member of the same family, SlOFP20 (Huang et al., 2009; Rodriguez et al., 2013; Wu et al., 2018). Whereas ovate is a null, the Slofp20 allele shows reduced expression and the effects of both mutations on ovary shape are already apparent at anthesis (Van der Knaap and Tanksley, 2001; Van der Knaap et al., 2014; Wu et al., 2018). This finding implies that the patterning mediated by these OFPs occurs early in the ontogeny of the ovary, perhaps immediately after organ initiation. SUN also affects ovary shape before anthesis and continues to promote fruit elongation immediately after fertilization (Van der Knaap and Tanksley, 2001; Xiao et al., 2009; Wu et al., 2011). sun is due to a transposon-mediated duplication event leading to high expression of the transposed gene during reproductive development (Xiao et al., 2008). Overexpression of SUN in both wild and cultivated tomatoes leads to evenly elongated fruit shape (Xiao et al., 2008). Interestingly, sun synergistically interacts with ovate and together the two promote growth at the proximal end to form a pear-shaped and pointed tomato fruit (Wu et al., 2015; **Figure 2**). ovate and sov1 also synergistically interact to form a pear-shaped tomato but with a round bottom shape (Wu et al., 2018; **Figure 2**). This suggests that obovoid organ shapes may be achieved by alleles from different sets of proteins or that the pathways intersect.

The expression of wild type OVATE in tomato is high in the IM/FM, and its expression reduces 8 days after floral initiation (dpi) (**Figure 3**). In contrast, the expression of wild type SlOFP20 is relatively low in the IM/FM and increases in 6 dpi buds, with a dramatic increase at 16 dpi (**Figure 3**). For SUN, wild type gene expression is very low (first two time points, LA1589) whereas in the NIL with the retrotransposon-mediated duplication (sun introgressed in the LA1589 background), SUN is highly expressed

during floral development (**Figure 3**). The initiation of the gynoecium primordia occurs at 6 dpi (Xiao et al., 2009), which is when OVATE, SlOFP20 and SUN are well expressed. At 8 and 13 dpi, the expression of OVATE and SUN respectively, is much reduced from expression levels at the earlier developmental stages coinciding with when these genes may function in development.

OVATE, SlOFP20 and SUN affect neither floral organ identity, nor the organization or number of floral organs (Wu, 2015). Instead, SUN, OVATE and SlOFP20 regulate organ elongation by altering the directions of cell division along the proximal-distal axis (Wu et al., 2011, 2018). Whereas SUN affects cell division along the entire proximal-distal axis, OVATE and SlOFP20 appear to have a specific role in anisotropic growth primarily at the proximal end of the ovary. In NILs that carry the ovate and sov1 mutant alleles, there is an increased number of cells in the proximal-distal direction and a reduced number of cells in the medial-lateral direction compared to wild type. Cell size is also enlarged but cell shape appears to change little in the mutant background (Wu et al., 2018). Thus, the effect of cell size and shape in fruit elongation is not clear and therefore, cell division patterns are thought to drive the shape of the ovate/sov1 fruits.

Certain OFPs and SUNs are likely to be involved in conserved mechanisms of morphology regulation across plant species. Genetic evidence indicates that the same subclade of OFPs, represented by Arabidopsis OFP1 and tomato OFP20, controls tomato fruit shape as well as aerial organ shapes in Arabidopsis (Wang et al., 2011), tuber shape in potato, and fruit shape in melon (Wu et al., 2018). Specifically, the potato (Solanum tuberosum L.) tuber shape QTL Ro has been fine-mapped in an outcrossing diploid F<sup>1</sup> population to a region on chromosome 10 that contains the potato ortholog of SlOFP20. There is also strong association between tuber shape and StOFP20 in a separate inbred diploid F<sup>2</sup> population. Very elongated tubers lack the StOFP20 gene, consistent with its role in the regulation of organ shape as found in tomato. In melon (Cucumis melo), fine mapping within the fruit shape QTL fsqs8.1 has identified CmOFP13 in a cross of Piel de Sapo and PI124112 (Wu et al., 2018). For SUN and other members of the IQD family, natural mutations affecting organ shape have been found in rice and species in the Cucurbitaceae family. Specifically, the rice gene GSE5 at the GW5/qSW5 locus encodes a SUN member closely related to Arabidopsis IQD25- 27 (Duan et al., 2017). The change in grain shape is due to increased cell proliferation in spikelet hulls. Interestingly, in cucumber and watermelon, a SUN member that is also most similar to Arabidopsis IQD25-27 likely controls fruit shape in these two species (Pan et al., 2017; Dou et al., 2018). Another rice SUN-like gene, OsIQD14 has been shown to affect rice grain shape and this member is most closely related to another subclade of the SUN/IQD family (Yang et al., 2018). Arabidopsis, overexpression of several IQD members leads to altered organ shapes. The overexpression of microtubular localized AtIQD16 and AtIQD11 resulted in elongated aerial organs with left-handed helical growth abnormalities similar to the phenotype of the tomato SUN overexpressors (Wu et al., 2011; Burstenbinder et al., 2017b). Overexpression of AtIQD14 results in organ twisting but not cell elongation as observed in AtIQD11 and AtIQD16 (Burstenbinder et al., 2017b), a phenotype that resembles that of tortifolia/spiral mutants (Furutani et al., 2000; Buschmann et al., 2004; Shoji et al., 2004). Furthermore, overexpression of plasma membrane localized IQD25 resulted in rounder leaves and larger cells, the opposite phenotype from that observed in overexpression of microtubule localized IQDs suggesting that IQD proteins can have diverse functions in regulating the cytoskeleton and cell elongation (Burstenbinder et al., 2017b). Thus, in addition to SUN in tomato, several members of this family have been associated with changing plant organ shape.

#### TRMs

A knockout mutation in tomato's TONNEAU1 Recruiting Motif 5 (SlTRM5) results in a slightly flatter fruit yet its effect is most strongly noticeable in the ovate/sov1 mutant background (**Figure 2**). The expression of wild type SlTRM5 is high in

IM/FM throughout floral development until 13 dpi (**Figure 3**), following similar expression dynamics as OVATE and SUN. SlTRM5 is a member of the Arabidopsis TRM1-5 subclade in which tomato carries only two TRM paralogs. At the cellular level, SlTRM5 controls the number of cells in the proximal-distal and medial-lateral direction such that the mutant allele Sltrm5 rescues the tomato fruit shape phenotype of ovate/sov1. TRM5 orthologs and close paralogs also appear to regulate organ shape in other crops. For example, the cucumber ortholog of TRM5 underlies the fs2.1 QTL controlling fruit shape (Wu et al., 2018). In rice, a major QTL for grain length encodes a TRM member in the TRM1-5 subclade. The discovery was made in three independent studies as GRAIN LENGTH ON CHROMOSOME 7 (GL7)/GRAIN WIDTH 7 (GW7)/SLENDER GRAIN ON CHROMOSOME 7 (SLG7) loci (Wang Y. et al., 2015; Wang S. et al., 2015; Zhou et al., 2015). Copy number variants at the GL7 locus contribute to grain size diversity (Wang Y. et al., 2015) and the increased expression of GW7/SLG7 increases grain length (Wang S. et al., 2015; Zhou et al., 2015). However, these studies show contrasting effects on the cellular mechanisms of grain shape changes. Higher expression of GW7 increased cell division in the proximal-distal direction and decreased cell division in the medial-lateral direction (Wang S. et al., 2015), which is similar to the effect of SlTRM5 on tomato fruit shape. On the other hand, increased expression of SLG7 increased cell length and decreased cell width with no changes in cell division (Zhou et al., 2015). In Arabidopsis, certain members of the TRM1-5 subclade control the elongation of various aerial organs. Overexpression of AtTRM1 (LONGIFOLIA2) or AtTRM2 (LONGIFOLIA1) leads to extremely long cotyledons, leaves, floral organs and siliques (Lee et al., 2006). On the other hand, loss-of-function mutations in AtTRM1 or AtTRM2 cause shortened siliques and cotyledons (Lee et al., 2006; Drevensek et al., 2012), which intriguingly mimic the phenotypes of AtOFP overexpressors (Wang et al., 2007). The more elongated leaf blades seen in the AtTRM1 and AtTRM2 overexpressors are due to increased cell expansion along the proximal-distal axis rather than an altered cell proliferation pattern (Lee et al., 2006).

#### MECHANISTIC INSIGHTS INTO THE REGULATION OF ORGAN SHAPE

#### Interaction Between OFPs and TRMs

As mentioned in the previous section and based on several studies, TRMs play a critical role in regulating organ shape. In tomato, TRMs were first discovered in a Yeast 2-Hybrid (Y2H) experiment using OVATE as bait. The goal of the experiment was to identify molecular interactants of OVATE to learn about cell division patterning mediated by OFP family members. A total of 11 out of 26 members of the TRM superfamily were identified in the screen. What set these OVATE-interactants apart from the other members of the TRM family was the conserved M8 motif (Van der Knaap et al., 2014; Wu et al., 2018). These findings suggest that the genetic interaction of TRM5 is through protein-protein interactions with OVATE via the TRM M8 motif. To validate the findings from Y2H, the interaction motifs were mapped in OVATE, SlOFP20, and several Ovate-interacting TRMs (Wu et al., 2018). OVATE and SlOFP20 both interact through highly conserved negatively charged amino acids in the OFP domain with TRMs via the highly conserved basic residue (K or R) in the TRM M8 motif. It is reasonable to conclude that the electrostatic interactions in the OFP domain and the M8 motif enable the interactions between these proteins (Wu et al., 2018).

The Y2H protein interactions have also been validated in a plant system. OVATE, SlOFP20 and several Ovateinteracting TRMs were expressed as fusion proteins in Nicotiana benthamiana leaf epidermal cells (Wu et al., 2018). When expressed alone, OVATE localizes in the cytoplasm and SlOFP20 is in the nucleus and cytoplasm. When SlTRM3/4 or SlTRM5 (members of the AtTRM1-5 subclade) are expressed alone, they localize to microtubules. Co-expression of OVATE and SlTRM5 dissociates SlTRM5 from microtubules and both proteins are found in the cytoplasm. On the other hand, coexpression of SlOFP20 and SlTRM5 causes the localization of SlOFP20 to microtubules coincident with SlTRM5. Coexpressions of OVATE or SlOFP20 with SlTRM3/4 both lead to a nearly complete dissociation of SlTRM3/4 from microtubules to the cytoplasm (Wu et al., 2018). These re-localizations are much reduced when mutants of OVATE, SlOFP20, SlTRM5,

and SlTRM3/4 lacking the interacting charged amino acid residues are co-expressed in N. benthamiana cells. These findings imply that relocalization occurs through physical protein interactions. Bifluorescence complementation assays further demonstrate that the charged amino acid residues of OVATE and SlTRM5, or OVATE and SlTRM3/4, are responsible for their interactions as well as relocalizations (Wu et al., 2018). The relocalization of OFPs and TRMs to different subcellular compartments upon interaction suggests that a dynamic balance between cytoplasmic- and microtubular-localized OFP-TRM protein complexes regulates cell division and organ growth. A mechanistic model describing the function of OFPs and TRMs to control organ shape is shown in **Figure 4**.

The subcellular relocalization when co-expressing OVATE or SlOFP20 on the one hand and SlTRM5 or SlTRM3/4 on the other hand suggests that certain OFPs play a critical role in localizing protein complexes. The discovery of the OFP-TRM module may provide an explanation for how other OFPs would serve as regulators in various distinct developmental processes by their effect on subcellular localization. AtOFP5 negatively affects the function of BLH1-KNAT3 complex in early embryo sac development (Pagnussat et al., 2007). The authors show that this is due to abnormal migration and positioning of embryo sac nuclei during megagametogenesis. They propose that this could be due to a change in the behavior of microtubules (Pagnussat et al., 2007) that serve as tracks for nuclear movement in plant cells (Vogelmann et al., 1981; Meindl, 1983; Mineyuki and Furuya, 1985). However, no evidence was available at that time to link the function of AtOFP5 with microtubule dynamics. Another example is offered by AtOFP4 that interacts with KNAT7 to regulate secondary cell wall formation (Li et al., 2011). cMTs participate in secondary cell wall development by directing the deposition of cell wall matrix components (Oda et al., 2005). The defects in secondary cell wall formation could be caused by an abnormal microtubule behavior or due to mislocalization of the CSCs due to the loss-of-function of AtOFP4. Again, this suggests that subcellular localization of CSCs may be disrupted by an OFP. Interestingly, SUN-like protein, AtIQD13 is also associated with positioning the CSCs by influencing the organization of the cMT arrays that guide them (Sugiyama et al., 2017). Other research has shown that the BEL1-like homeodomain 1 (BLH1) protein is shuttled from the nucleus to the cytoplasm when interacting with AtOFP1 or AtOFP5. It was therefore proposed that AtOFPs affect the activities of TALE transcription factors by altering their subcellular localization (Hackbusch et al., 2005). AtOFP1 was also identified as a protein partner of AtKu70, which plays a role in non-homologous end-joining DNA repair (Wang et al., 2010). Interestingly, the centromeric function of Ku70 depends on the presence of microtubules (Cabrero et al., 2013). Therefore, it is reasonable to propose that AtOFP1 may function in DNA-repair by affecting the anchoring of AtKu70 to microtubules. Thus, the OFP-TRM module could explain these seemingly unrelated pathways where OFP controls the subcellular localization of protein complexes. This idea is in stark contrast with the notion in the literature that AtOFPs are transcriptional repressors (Wang et al., 2007, 2011). However, these conclusions were primarily made based on protoplast expression assays as well as expression correlation when overexpressing AtOFP1 (Hackbusch et al., 2005; Wang et al., 2007) and thus transcriptional repression was not validated in intact plants.

#### TRMs, TTP and the Cytoskeleton

A broader function of the TRMs is their role in assembling the TTP (TON1-TRM-PP2A) complex. Before members of this family were associated with the TTP complex, two members namely LONGIFOLIA (LNG) 1 and LNG2, were identified to control organ shape in Arabidopsis (Lee et al., 2006). The entire TRM family, however, was identified in Arabidopsis in a Y2H study using TON1 as bait (Drevensek et al., 2012). The Arabidopsis TRMs consists of 34 members and all contain the TON1-interacting M2 motif at the C terminus. A typical TRM in Arabidopsis is AtTRM1, a microtubule-associated protein that localizes to cMTs in vivo and binds microtubules in vitro. AtTRM1 recruits TON1 to microtubule arrays in N. benthamiana leaf cells. A subset of TRMs target the TTP complex to microtubules (Drevensek et al., 2012). The TTP complex has been proposed to regulate the organization of microtubule arrays and PPB formation, and thus cell division patterns and cell growth (Camilleri et al., 2002; Azimzadeh et al., 2008; Drevensek et al., 2012; Spinner et al., 2013; Schaefer et al., 2017). Throughout interphase in plant cells, microtubules are found just beneath the plasma membrane in the cell cortex. These cMTs determine cell shape as they form patterns in the absence of focused nucleation centers like centrosomes in animal and fungal cells. The nucleation of new cMTs is geometrically constrained (Fishel and Dixit, 2013). Most initiate from a nucleation site on the side of a parent microtubule. These nucleation sites contain γ-tubulin and associated γ-tubulin complex proteins (Nakamura et al., 2010; Murata and Hasebe, 2011) to form the γ-TuRC. After nucleation, new microtubules elongate at about 40<sup>o</sup> from the parent microtubule (Chan et al., 2009; Nakamura et al., 2010; Murata and Hasebe, 2011). New microtubules also grow parallel to the parent microtubule and move alongside existing microtubules by polymer treadmilling (Chan et al., 2009; Nakamura et al., 2010). The branched form of nucleation is the dominant pattern, while parallel nucleation occurs about half as frequently. About 1–2% of nucleation events occur de novo, where a new microtubule elongates independently of an existing microtubule (Shaw et al., 2003; Chan et al., 2009; Nakamura et al., 2010). As a component of the TTP complex, the Arabidopsis thaliana B" subunit of protein phosphatase 2A is encoded by the TONNEAU2/FASS (TON2) gene (Camilleri et al., 2002), and microtubule branching nucleation is specifically promoted by this regulatory subunit (Kirik et al., 2012). In ton2- 15 mutants, the frequency of microtubule branching nucleation is reduced 4-fold while the frequency of parallel nucleation is increased 2.4-fold. The branching angle of new microtubules is unchanged. In hypocotyl cells, loss of TON2 function also results in the inability of microtubule arrays to reorient in response to light, suggesting an essential role for TON2 and microtubule branching nucleation in the reorganization of microtubule arrays (Kirik et al., 2012). It has been postulated that TON2 may influence the orientation of initial polymerization through a direct interaction with or phosphorylation of a component

of the γ-TuRC (Fishel and Dixit, 2013). Microtubule assembly at the γ-TuRC is also partly regulated by the microtubule severing protein katanin (Nakamura et al., 2010) and katanin1 mediated microtubule rearrangement is proposed to play a role in regulating rice grain shape controlled by a SUN-like gene, OsIQD14 (Yang et al., 2018). This finding supports the notion that the regulation of organ shape might be functionally linked by the OFP-TRM and SUN/IQD pathways.

The formation and function of different microtubule arrays are regulated by microtubule nucleation, dynamics and stability. Microtubule assembly is a polarized process starting from one or several MTOCs (Murata and Hasebe, 2011). The centrosome is the major MTOC in animal cells to recruit and modify cell cycle proteins (Wu and Akhmanova, 2017). Even though vascular plant cells lack centrosomes, TTP components have sequence similarity to animal centrosomal proteins. For example, the N-terminus of TON1 has sequence similarity to FOP and OFD1, proteins required for microtubule anchoring and stability within the centrosome, respectively (Yan et al., 2006; Azimzadeh et al., 2008; Singla et al., 2010). In addition, three TRM motifs (i.e., M3-M4-M2) are found in the human centrosomal protein CAP350, which interacts with FOP. The C-terminal M2 motif in CAP350 is responsible for FOP recruitment to the human centrosome and facilitates microtubule anchoring within the centrosome (Yan et al., 2006). Similarly, TON1 and the TTP complex bind to the PPB, which marks the future division plane by promoting spindle bipolarity and limiting spindle rotation to ensure properly patterned cell division. The consequence of the sequence similarity and overlapping motifs between TTP complex proteins and animal centrosomal proteins may be the functional similarities among the complexes in plant and animal cell division (Schaefer et al., 2017). The TTP complex is required for proper PPB assembly and division plane establishment (Spinner et al., 2013). The PPB is an array of microtubules and actin filaments that forms a ring at the cell periphery during G2 and persists throughout prophase. Although the PPB is disassembled as the nuclear envelope breaks down and the mitotic spindle forms, its position precisely correlates with the position of the future division plane. The spatial information of the PPB is preserved by selective recruitment and depletion of proteins that lead to the generation of the cortical division zone and the precise positioning of the cell plate during cytokinesis (Rasmussen et al., 2011; Rasmussen and Bellinger, 2018). TON1A, TON1B and the PP2A subunit FASS/TON2 (in Arabidopsis) or DISCORDIA1/ALTERNATIVE DISCORDIA1 (FASS/TON2 orthologs in monocots) are required for PPB formation. Knockout mutants lack PPBs and have incorrectly positioned division planes (Camilleri et al., 2002; Azimzadeh et al., 2008; Wright et al., 2009; Spinner et al., 2010). Thus, a potential mechanistic link between OFPs, OVATE-interacting TRMs and cell patterning is established through interactions with the TTP complex thereby regulating organ shape.

Whether the PPB is absolutely required for division plane patterning is not clear as cells in certain tissues appear to divide without TON1a and PPB formation in other tissues (Zhang et al., 2016a; Costa, 2017). Further insights about the function of the PPB and certain TRMs show that TRM7 is a specific PPB marker whereas TRM6 and TRM8 are constitutively expressed throughout the cell cycle (Schaefer et al., 2017). The frequency of normal PPB formation is reduced in trm7 mutants and no PPBs are found in the trm678 triple mutants. Cells with disrupted PPB formation retain the capacity to define a cortical division zone but lose precision in the orientation of this division zone. Intriguingly, mutant trm678 plants are fertile with normal organs that do not exhibit aberrant cell division patterns. The phragmoplast-orienting kinesin 1 (POK1) is a factor controlling the timing and efficiency of the cortical division zone (Lipka et al., 2014). In trm678 mutant cells, even though POK1 targeting to the cell cortex is altered in the absence of the PPB, POK1 still forms a cortical ring that corresponds with the division zone (Schaefer et al., 2017). These results suggest that the PPB may be less of a causal determinant of the cell division plane and more of a regulator to ensure the fidelity of a division plane defined by another mechanism. Regardless, the position of the cell division plane has profound impact on the shape of plant organs, and therefore, much remains to be discovered of how plane positioning is regulated.

### SUN/IQD AND THE CYTOSKELETON

The tomato SUN/AtIQD12 is a member of the IQ67 domain (IQD) protein family (Xiao et al., 2008; Huang et al., 2013). The IQ67 domain of IQD proteins is a conserved region of 67 amino acids and contains up to three regularly spaced IQ motifs which promote calmodulin (CaM) binding in the presence of Ca2<sup>+</sup> (Rhoads and Friedberg, 1997; Abel et al., 2005, 2013). Ca2<sup>+</sup> is a common secondary messenger in all eukaryotes and is used to regulate many cellular processes in response to both cellular and environmental stimuli, including cell division and shape (Cardenas, 2009; Steinhorst and Kudla, 2013; Burstenbinder et al., 2017b). The IQ67 domain of several IQD proteins interacts with CaM demonstrating that this family of proteins may serve as a large class of CaM binding proteins in plants (Burstenbinder et al., 2013, 2017b; Yang et al., 2018). The founding member of the IQD family, IQD1, localizes to microtubules with CaM2 in Arabidopsis and both IQD1 and IQD20 were found to interact with CaMs by Y2H, suggesting that IQDs may integrate Ca2<sup>+</sup> sensing in regulation of the cytoskeleton (Levy et al., 2005; Burstenbinder et al., 2013).

Expression analyses in N. benthamiana showed that the N-terminus of most IQD proteins localizes to microtubules and that half of the IQDs localize to the plasma membrane (Burstenbinder et al., 2017b). There is also evidence that certain IQDs have differential subcellular localization dependent on the cell cycle stage (Wendrich et al., 2018). Cells in plants overexpressing AtIQD16 had altered orientation of cMTs with more oblique aligned microtubules and significantly elongated cells. Colocalization of IQDs with CaM also suggests that IQD proteins are capable of sequestering or recruiting CaM to specific subcellular domains. Subcellular localization of several Arabidopsis IQDs (IQD12, IQD22, IQD24, and IQD25) showed punctate structures that are reminiscent of regions within the plasma membrane, which may act as signaling centers in the cell.

IQD14 is the rice ortholog to the Arabidopsis IQD15-18 subclade and loss of function alleles result in shorter and wider grains than wild type rice (Yang et al., 2018).

OsIQD14 was found to localize to the nucleus and cytoplasm and also in punctate locations along the microtubules, suggesting that the protein may function at specific points in microtubule regulation or on particular microtubule structures. Interestingly, expression of rice IQD14-GFP N- or C-terminal regions in N. benthamiana showed that the C-terminal region localized to the microtubules while the N-terminal region localized to the nucleus. This result is similar to the localization observed in the Arabidopsis IQD15-18 clade where the full length IQD protein was found on both microtubules and in the nucleus (Burstenbinder et al., 2017b; Yang et al., 2018). Rice IQD14 was also found to interact with Arabidopsis SPR2 by Y2H, and the orthologous Arabidopsis IQD15-18 subclade members were found to interact with SPR2 and CaM as well (Wendrich et al., 2018; Yang et al., 2018). SPR2 is a microtubule binding protein involved in protecting the minus end of microtubules and promoting severing and reorientation of the cMT arrays. SPR2 generally localizes to microtubules and does not distinguish dynamic from stable microtubules. However, IQD proteins may serve to direct the location of SPR2 function to specific regions to regulate reorganization the cytoskeleton in response to a certain signal (Buschmann et al., 2004; Shoji et al., 2004; Yao et al., 2008; Nakamura et al., 2018).

Auxin has been suggested to influence MT dynamics, but the mechanism is unclear. However, recent studies suggest that auxin-mediated cytoskeletal changes may involve IQD proteins (Chen et al., 2014; Wendrich et al., 2018) and IQDs are likely downstream targets of AUXIN RESPONSE FACTOR5/ MONOPTEROS (Boer et al., 2014; Moller et al., 2017). It has been proposed that IQDs in low auxin/ Ca2<sup>+</sup> environments do not bind to CaM and instead bind to SPR2, inhibiting its function. This results in stabilized microtubules and a less dynamic cytoskeleton (Leong et al., 2018; Nakamura et al., 2018). Auxin leads to an increase in Ca2+, which promotes CaM binding to the IQD, and CaM binding then prevents IQDs from binding to SPR2. Unbound SPR2 can then bind the minus end of microtubules and promote microtubule branching and changes to cytoskeletal architecture in response to auxin (Wendrich et al., 2018). IQD proteins may function to resolve these signals and changes in Ca2<sup>+</sup> levels within developing tissues in response to the environment, thereby directing cell elongation and expansion to ultimately drive organ shape.

As organ shape is controlled by both cell division and directed cell expansion, the regulation of microtubule dynamics is important in both of these processes to determine morphology. Some IQD proteins regulate cytoskeletal architecture by guiding the formation of ROP domains in the plasma membrane. ROPs are plant specific Rho GTPases with diverse functions (Yalovsky, 2015). One function of ROPs is in organizing the microtubule and actin cytoskeleton to determine a cell's final shape. Some ROPs have been shown to promote aggregation of fine actin filaments and inhibit the assembly of organized microtubule arrays (ROP2 and ROP4), while ROP6 has been shown to have the opposite role and promotes accumulation of organized microtubule arrays (Fu et al., 2005, 2009; Ivakov and Persson, 2013). Furthermore, downstream of ROP6, the microtubule severing protein katanin is activated and promotes microtubule reorganization (Hamant, 2013; Lin et al., 2013). Formation of distinct subcellular domains of ROPs with these opposing functions can alter cytoskeleton composition and fine-tune the overall cell shape (Fu et al., 2002, 2005, 2009; Lin et al., 2013). The plasma membrane localization of Arabidopsis IQD13 has been shown to regulate ROP function in the xylem by promoting cMT growth and interaction with the membrane surface, thus restricting the formation of ROP11 domains (Oda and Fukuda, 2012; Sugiyama et al., 2017). In the presence of IQD13, the active ROP11 is restricted within narrow plasma membrane domains where it can recruit additional proteins to ultimately deplete the region of cMTs and form narrow pits in the secondary cell wall. In the absence of IQD13 or in the presence of truncated IQD13 lacking the plasma membrane associated domain, ROP11 forms circular domains that are independent of the cMTs and ultimately forms round pits in the secondary cell wall (Sugiyama et al., 2017). Secondary cell wall pit architecture is further refined by an interplay between the restriction of ROP11 to narrow domains by IQD13 and the impairment of ROP11 restriction and the resulting delineation of cMTs by CORTICAL MICROTUBULE DISORDERING1 protein (CORD) (Sasaki et al., 2017). AtIQD5 has recently been shown to regulate microtubule dynamics that affect cMT organization and subsequent cell shape formation in leaf pavement cells. IQD5 is enriched in lobed regions of pavement cells where cMTs are organized in parallel arrays and cell expansion is restricted (Liang et al., 2018; Mitra et al., 2018). In iqd5 mutant cells, regions lacking IQD5 expression no longer form lobes and have altered cellulose deposition (Mitra et al., 2018). In summary, these results suggest that IQDs may regulate microtubule organization in distinct subcellular regions through interactions with ROPs in order to impact a cells final shape and ultimately organ shape.

Since organ shape is also influenced by anisotropic cell expansion (Maugarny-Cales and Laufs, 2018), the regulation of this expansion is an important factor controlling morphology. While cell expansion is driven by isotropic turgor pressure, the direction of expansion is controlled by the pattern of cellulose microfibrils, which are generally deposited perpendicular to the axis of expansion. This pattern resists turgor driven expansion in the direction parallel with the cellulose microfibrils, thus promoting expansion in the perpendicular direction (Szymanski and Cosgrove, 2009; McFarlane et al., 2014). cMTs regulate the deposition pattern of cellulose microfibrils within the cell wall by interacting with cellulose synthase complexes in the plasma membrane (Paredez et al., 2006; McFarlane et al., 2014). Therefore, the organization of microtubules, in part controlled by IQDs also impacts the arrangement of cellulose microfibrils which will impact how the cell wall expands and the ultimate shape of the cell. While randomly aligned cMTs are located in the cytoplasm away from the plasma membrane, cMTs closely anchored to the plasma membrane are organized in parallel bundles to preferentially serve as tracks for cellulose synthase movement (Barton et al., 2008). These cMTs rarely display lateral displacement from their parallel organization due to their tight

association with the plasma membrane (Shaw et al., 2003). Cellulose microfibrils are synthesized by large CSCs, composed of 18–36 cellulose synthase subunits and their accessory proteins (McFarlane et al., 2014). Cellulose synthases are assembled and matured in the Golgi and sorted by the Trans Golgi network into small cellulose synthase compartments/ microtubule-associated cellulose synthase compartment vesicles for its secretion to the plasma membrane (Crowell et al., 2009; Gutierrez et al., 2009; Sampathkumar et al., 2013; Zhang et al., 2016b). CSCs are tethered to cMTs through cellulose synthase interactive protein 1 (CSI1) (Gu et al., 2010; Bringmann et al., 2012; Li et al., 2012; Mei et al., 2012), which determines their trajectory along the cMTs (Paredez et al., 2006). Delivery to the plasma membrane is also mediated by CSI1 and fusion with the plasma membrane is mediated by the plant specific protein PATROL1 and the exocyst complex (Zhu et al., 2018). AtIQD1 interacts with KLCR1 in Y2H and in planta where its recruitment to microtubules is dependent on IQD1 (Abel et al., 2013; Burstenbinder et al., 2013). KLCR and cellulose-microtubule uncoupling (CMU) are in the same protein family (Burstenbinder et al., 2017a) and the same proteins in Arabidopsis where At4g10840 encodes KLCR1 and CMU1 and At3g27960 encodes KLCR2 and CMU2 (Burstenbinder et al., 2013; Liu et al., 2016). CMU1 and CMU2 are localized as static puncta along microtubules. Disruption of CMU function causes lateral microtubule displacement. This compromises microtubule-based guidance of CSCs leading to cell twisting and altered growth (Liu et al., 2016). Together these results indicate that IQD proteins regulate how microtubules direct CSCs and the pattern of cellulose microfibril deposition, which influences cell expansion and organ shape. It is intriguing that AtOFP4 also regulates cell wall formation through its interaction with KNAT7 (Li et al., 2011). Thus, both IQDs and OFPs may influence organ shape at the cellular level through their regulation of cell wall formation.

With respect to the regulation of organ shape by SUN/IQDs and the role in cell division patterning, these proteins may directly influence this process. AtIQD5 is localized to the PPB, spindle, and phragmoplast of dividing root cells (Liang et al., 2018). ROPs may also be involved in PPB formation (Oda, 2018) because two ROP GTPase activating proteins interact with POK1 and are required for

#### REFERENCES


accurate orientation of the PPB, phragmoplast, and cell plate (Stockle et al., 2016). POK1 is required for division plane maintenance (Lipka et al., 2014; Stockle et al., 2016) and its function is also influenced by TRMs (Schaefer et al., 2017). It is plausible that SUN/IQD and TRM proteins coordinate cell division planes, which contribute to organ shape.

#### CONCLUSION

A model to describe organ shape in the context of interactions of SUN/IQD, OFP and TRM, and associations with microtubules is shown in **Figure 4**. Assuming round shape as the default, the down regulation of OFP would lead to the association of more TRMs to the microtubules and hence elongated shape. This organ shape might also be attained via a similar mechanism when up-regulating TRM (**Figure 4**). Conversely, up-regulation of OFP and down regulation of TRM would result in rounder or even flat shapes due to less microtubular association of TRMs. SUN/IQD proteins are often found at the microtubules where their interaction with SPR2 and CaM might lead to altered cytoskeleton activities. Higher expression of SUN/IQD would lead to more association with the microtubules and hence elongated shape. Together with findings in Arabidopsis and crop plants, further information has shown that OFPs, TRMs and SUN/IQDs impact microtubular activities, offering mechanistic insights into how the different shapes of plant organs are realized.

#### AUTHOR CONTRIBUTIONS

ML and EvdK wrote and edited the review with significant contributions and further edits from SW, AS, and YW.

# FUNDING

This research was supported by National Science Foundation IOS 0922661, NSF IOS 1564366, and AFRI 2017-67013-26199 of the USDA National Institute of Food and Agriculture.





**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Lazzaro, Wu, Snouffer, Wang and van der Knaap. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Genetic Analysis of Chinese Cabbage Reveals Correlation Between Rosette Leaf and Leafy Head Variation

XiaoXue Sun1,2† , Shuangxia Luo<sup>1</sup>† , Lei Luo<sup>1</sup> , Xing Wang<sup>1</sup> , Xueping Chen<sup>1</sup> , Yin Lu<sup>1</sup> , Shuxing Shen<sup>1</sup> \*, Jianjun Zhao<sup>1</sup> \* and Guusje Bonnema1,2 \*

#### Edited by:

Jacqueline Batley, The University of Western Australia, Australia

#### Reviewed by:

Ryo Fujimoto, Kobe University, Japan Yang Zhang, St. Jude Children's Research Hospital, United States

#### \*Correspondence:

Shuxing Shen shensx@hebau.edu.cn Jianjun Zhao jjz1971@aliyun.com Guusje Bonnema guusje.bonnema@wur.nl †These authors have contributed equally to this work

#### Specialty section:

This article was submitted to Plant Breeding, a section of the journal Frontiers in Plant Science

Received: 29 June 2018 Accepted: 12 September 2018 Published: 04 October 2018

#### Citation:

Sun X, Luo S, Luo L, Wang X, Chen X, Lu Y, Shen S, Zhao J and Bonnema G (2018) Genetic Analysis of Chinese Cabbage Reveals Correlation Between Rosette Leaf and Leafy Head Variation. Front. Plant Sci. 9:1455. doi: 10.3389/fpls.2018.01455 <sup>1</sup> Key Laboratory of Vegetable Germplasm Innovation and Utilization of Hebei, Collaborative Innovation Center of Vegetable Industry in Hebei, Department of Horticulture, Hebei Agricultural University, Baoding, China, <sup>2</sup> Plant Breeding, Wageningen University & Research, Wageningen, Netherlands

To understand the genetic regulation of the domestication trait leafy-head formation of Chinese cabbages, we exploit the diversity within Brassica rapa. To improve our understanding of the relationship between variation in rosette-leaves and leafy heads, we phenotyped a diversity set of 152 Chinese cabbages. This showed correlation between rosette-leaf traits and both head traits and heading capacity. Interestingly, the leaf number of the mature head is not correlated to heading degree nor head shape. We then chose a non-heading pak choi genotype to cross to a Chinese cabbage to generate populations segregating for the leafy head traits. Both a large F2 (485 plants) and a smaller Doubled Haploid (88 lines) mapping population were generated. A high density DH-88 genetic map using the Brassica SNP array and an F2 map with a subset of these SNPs and InDel markers was used for quantitative trait locus (QTL) analysis. Thirty-one quantitative trait loci (QTLs) were identified for phenotypes of rosette-leaves in time and both heading degree and several heading traits. On chromosome A06 in both DH-88 and F2-485 QTLs for rosette leaf length and petiole length at different developmental days and an F2 QTL for head height co-located. Variation in head height, width and weight all correlate with variation in heading degree with co-locating QTLs, respectively, on chromosome A03, A05, and A08 in F2-485. The correlation between rosette-leaf and heading traits provides not only insight in the leaf requirements to form a head, but also can be used for selection by Chinese cabbage breeders.

Keywords: Chinese cabbage, rosette leaf, leafy head, heading degree, correlation, co-location of QTLs

# INTRODUCTION

Brassica rapa is an economical important species that includes many morphotypes that are cultivated as vegetable, oil and fodder crops (Prakash and Hinata, 1980; Gómez-Campo, 1999; Hanelt, 2001). It contributes many important leafy vegetable crops grown worldwide, with main production and consumption area's in Asia and Europe (Dixon, 2007; Ge et al., 2011). The two

**41**

most important B. rapa leafy morphotypes are Chinese cabbage with leafy heads and pak choi's with flat smooth leaves and fleshy petioles that do not form a head. The general hypothesis is that Chinese cabbages have been domesticated from non-heading leafy pak choi's in China (Li, 1984). Compared to the non-heading accessions, the heading types can be transported easily, stored long time and have better cold tolerance and higher yield (Daly and Tomkins, 1997).

Chinese cabbages go through four stages of vegetative growth: seedling, rosette, folding and heading stage (Ge et al., 2011). During seedling stage the leaves are round and have long petioles. Rosette leaves (RL) differentiate at the rosette stage and become large and round, with short petioles. At the end of the rosette stage, the folding leaves (FL) are also large and round with short petioles, but curve more inward and thereby form a mold for further head development. Unlike species that do not form a leafy head, the rosette and FL in heading plants remain metabolically active during the entire head development in order to supply the head with the products by photosynthesis (Yu et al., 2000; Wang et al., 2012). Leafy heads of Chinese cabbage are formed following the rosette and folding stages, when heading leaves develop an extreme upward and inward curvature giving rise to the leafy head that in effect becomes a nutrient storage organ. However, details of leaf development are rarely described and the correlation between variation of leaves in rosette/folding stage with timing of head formation and head shape is still unclear. As far as we know, only few genetic studies have been conducted on variation in rosette leaf traits and the relationship with head formation. Mao et al. (2014) phenotyped the RL and head shape at the heading stage and found a correlation between their morphology. Generally round heads had flat RL, cylindrical heads had RL with wavy margins and cone-like heads had RL that curved inward. The same team also recently published about the effects of microRNAs on leaf curvature, which is also highly relevant to understand leafy head formation in Chinese cabbage (Wang et al., 2014). In addition, several genomic studies are conducted that describe gene expression profiles during Chinese cabbage development (Wang et al., 2012, 2013; Gau et al., 2017). When these are combined with phenotyping of both the complete plant and the anatomy of leaves and meristems, our understanding of leafy head formation will increase.

Head formation is a quantitative trait, controlled by multiple genes, however, genetic mechanism of head formation is still unclear. Currently, few studies have identified quantitative trait loci (QTLs) for heading traits in Chinese cabbage, such as head weight, head diameter, head height and number of leaves. Kubo et al. (2010) constructed a map based on a mapping population consisting of 188 F2 plants derived from a cross between a Chinese cabbage and a vegetable turnip. For head formation, two QTLs were detected in both 2005 and 2007, with no common head formation QTLs detected. Ge et al. (2011) phenotyped 154 F3 families and identified twenty-seven QTL distributed over all nine linkage groups for gross weight, head weight, head length, head width, numbers of wrapper leaves and head forming leaves. Inoue et al. (2015) evaluated eight traits over 2 years (2010 and 2011) for a F2 population of 188 plants from a cross between two Chinese cabbage cultivars. They reported QTL for four heading traits (degree of head-top-leaf overlap, head height, head weight and the ratio of head height to head diameter) on linkage groups A1, A5, A7, and A8 in both years. Xiao et al. (2013) used a genetical genomics approach to study variation in leaf shape and size in B. rapa. For this purpose they used a DH population derived from a cross between non-heading pak choi and an oil type, yellow sarson. They identified several candidate genes for leaf variation and also revealed correlations between leaf traits and both plant architecture and flowering time.

In this paper we aimed to understand further the complexity of head formation during Chinese cabbage development. To achieve this goal, we performed the first global analysis of the correlation between leaf traits during development and head formation at heading stage using 152 Chinese cabbage cultivars. The data will substantially increase knowledge of Chinese cabbage heading related traits during leafy head development and provide research directions for QTL studies of heading related traits. We also generated an F1 population derived from a cross between heading Chinese cabbage (CC-48) and non-heading pak choi (PC-101). A doubled haploid population (DH-88) and an F2 population (F2-485) were developed from F1 by microspore culture and self-pollination, respectively. DH-88 and F2-485 linkage maps were constructed using 66 DH-88 and 485 F2- 485 individuals for QTL analysis to detect loci related to RL related traits at rosette stage and leafy head formation at heading stage.

### MATERIALS AND METHODS

#### Plant Materials

#### Chinese Cabbage Genotypes for Morphological Traits Analysis

A collection of 152 B. rapa L. ssp. pekinensis (Chinese cabbage) accessions from different locations of China was used in this study (**Supplementary Table S1**). Among them, 69 accessions were inbred lines and 83 accessions were hybrid varieties (F1). Plants were grown in a random block design in the field in Hebei province in China, September 2013. Three blocks with fifteen plants of each accession were used for collection of morphological trait data.

#### DH-88 and F2-485 Mapping Populations

The mapping populations used in this study consisted of a doubled haploid population (DH-88) and an F2 population (F2- 485) (**Figure 1A**). The DH-88 and F2-485 have the same father line, heading Chinese cabbage CC-48 (B. rapa ssp. pekinensis: CGN06867) and mother line, non-heading pak choi PC-101 (B. rapa ssp. chinensis: CGN13926). The F1 was obtained from a cross between CC-48 and PC-101. The 66 lines of DH-88 were obtained by F1 plant microspore culture. From this F1, 485 F2 plants were produced after selfing.

For QTL analyses, the experiments were conducted in Hebei Agriculture University, Hebei, China. Sixty-six DH-88 lines, plus their parental lines and their F1's were sown in seeding soil in

FIGURE 1 | The strategy for constructing the Chinese cabbage (DH-88 and F2-485) populations (A) and description of the measurements of leaf and heading traits (B). Head degree (HDe), head diameter (HD), head height (HH), core length (CL) and core width (CW) were recorded for the leafy heads. Leaf length (LL), leaf width (LW), leaf mid vein length (LVL), leaf mid vein width (LVW) and leaf mid vein thickness (LVT) were recorded for leaves.

the greenhouse, transplanted September 2014 to an open field at the Hebei Agriculture University, and grown under short day conditions (11 h photoperiod) until December 2014. The 485 F2 seeds were sown on the same experimental field September 2015, and grown under natural conditions until November 2015.

# Trait Evaluations and Phenotypic Data Analyses

The collection of 152 Chinese cabbage accessions were evaluated for 17 agronomic traits during heading stage and the traits scored in each individual plant are listed in **Table 1**. Five plants as biological replicates for each accession were evaluated for their agronomic traits.

For QTL analysis, four traits of DH-88 and eight traits of the F2-485 population were evaluated at rosette and heading stages. In rosette stage, the Sixth rosette leaf was phenotyped for leaf width (LW), leaf length (LL), leaf petal length (PL) and leaf trichomes (LT) at 36, 41, 44, 47, 50, and 53 days after sowing for the DH-88 and 30, 34, 37, 41, 44, and 48 days after sowing for the F2-485 population. At the heading stage, the degree of heading (HDe) was classified on a 1–4 scale (1) nonheading, (2) leaf up curling, (3) loose heading, and (4) tight heading) by visual observation before maturity of the heads for both the DH-88 and F2-485 population (**Figure 1B**). At maturity, head diameter (HD), head height (HH), and head fresh weight (HW) were measured after discarding the loose outer leaves of all F2-485 plants at the same time. In DH-88 population, five biological replicates per DH line were used for traits observation.

The mean, standard deviation and variance were calculated for all traits by SPSS (Statistical Package for the Social Sciences, IBM) software. R package ("corrplot" package<sup>1</sup> ) with Pearson correlation coefficient was used to calculate the correlation coefficients between traits and between populations.

<sup>1</sup>https://www.r-project.org

#### Genetic Map Construction and QTL Mapping DNA Extraction

Genomic DNA was extracted from young leaves of DH-88 and F2-485 plants by the CTAB method (Murray and Thompson, 1980; Rogers and Bendich, 1985). Subsequently, the DNA was checked on gel, diluted to the same concentration with water and used to screen for SNP markers.

#### DH-88 and F2-485 Population Genetic Maps

TraitGenetics (Germany) developed a 60 K B. napus array. The parental genotypes CC-48 and PC-101 were screened for this array, which resulted in 5795 polymorphic SNPs, which were used to construct the genetic map. The DH-88 individuals were genotyped for a set of 5795 SNPs covering the 10 linkage groups. A subset of 899 specific SNPs was used for constructing the DH88 genetic linkage map (**Supplementary Table S2**) in this study. The F2-485 genetic map was based on 99 InDel markers and 36 SNP markers (**Supplementary Table S3**). The InDel markers were designed using the B. rapa data base<sup>2</sup> and SNP markers were selected from the TraitGenetics set. SNP markers were used to anchor DH-88 and F2-485 populations. InDel markers were used for F2-485 genotyping by Native-PAGE method. SNP genotyping was carried out by high resolution melting analysis of small amplicons and HRMA was performed on 96-Well LightScanner System. Genetic maps for both populations were constructed by using JoinMap4.1 (Van Ooijen, 2006). After creating population nodes, LOD scores 8.0 to 10.0 were used to assign the markers into linkage groups (LGs) and Kosambi's mapping function (Kosambi, 2016) was used to construct genetic maps. Consensus LGs of the two populations were constructed with common markers in both populations and drawn using Mapchart 2.3 (Voorrips, 2002).

<sup>2</sup>http://brassicadb.org/brad/


TABLE 1 | List of 17 evaluated phenotypic traits of 152 Chinese cabbage accessions, their ranges and descriptive statistics.

<sup>1</sup>Std. Deviation; <sup>2</sup>Coefficient of Variation (CV) = (Std. Deviation / Mean) <sup>∗</sup> 100%.

#### DH-88 and F2-485 Population QTL Mapping of Heading Related Traits

Based on the genetic map, QTL mapping analysis was performed using MapQTL6 (Van Ooijen, 2009). QTL regions for testing traits were identified by interval mapping (IM) analysis. Based on the results of IM, multiple- QTL model mapping (MQM) was performed using the closest markers as cofactors. A QTL was considered significant when LOD thresholds were estimated with a genome-wide confidence level p < 0.005 (Doerge and Churchill, 1996).

# RESULTS

### Diversity and Correlation Analyses of 17 Heading Related Traits Scored Over the 152 Chinese Cabbage Accessions

Phenotypic data for seventeen traits scored for the diversity set of 152 Chinese cabbages are shown in **Table 1**. All the quantitative traits showed normal distribution (**Supplementary Figure S1**). The coefficient of variation for "head-top-leaf overlap level" (HtO, 0.65) was higher than that of other traits, while coefficient of variation for "heading degree" (HDe, 0.13) was the lowest among the seventeen investigated traits.

The 17 traits were divided into three clusters based on the outcome of principal component analysis (**Figure 2A**). HDe, CL, HD were classified as "heading capacity traits" with regard to the Chinese cabbage head formation. Thirteen traits (LVL, LL, PL, HtO, HI, HH, CW, PW, HW, LVW, LVT, PD, and LW) were classified as "head type traits." Based on the Eigen Value, CW, PW, HW, LVW, LVT, PD, LW traits mainly belong to the Head Trait component, but are also partly related to "Head capacity trait" component. This is also supported by Pearson correlation analysis (**Figure 2B**). Significant phenotypic correlations were found within the head type trait and within the heading capacity groups. CL and HDe clustered into one group and LVL, LL, PL, HtO, HI, HH formed another group, with significant correlation among these traits. The traits shared between Head Trait components and Heading capacity component grouped together based on their correlation rate. Only head leaf number (HLN) was negatively correlated with other traits in both component and correlation analysis.

# QTL Analysis in DH-88 and F2-485

# Linkage Map Construction

#### **Construction of DH-88 linkage maps**

The 5795 SNPs were profiled over the DH-88 population, which resulted in 899 bins distributed over ten LGs spanning a total genetic distance of 1469.5 cM and average distance between markers is 0.92 cM. The LGs correspond to the Chinese cabbage reference map. We only found one inversion of adjacent SNP bins located at the bottom of A10 and few inconsistent bins located in chromosomes A02 and A07, respectively, compared to the Chiifu Chinese cabbage reference genome (**Supplementary Figure S2**).

#### **Construction of F2-485 linkage map**

The DH-88 population consists of only 66 DH lines, limiting mapping precision and detection of significant QTLs. In order to identify significant QTLs, a large F2-485 population with 485 F2 plants was generated and a linkage map was constructed. A variety of different types of molecular markers were used for the construction of this linkage map. These included 36 SNPs and 99 InDel markers distributed over the ten LGs spanning a genetic distance of 1145.7 cM with an average distance of 8.93 cM between markers (**Supplementary Figure S3**). In general, marker order was in agreement with the Chinese cabbage genome sequence v1.5<sup>2</sup> with few exceptions. There were 36 common SNPs markers between DH-88 and F2-485 linkage

maps (**Supplementary Table S4**). Common markers assist in comparison of the QTL positions on the two different linkage maps.

#### Correlation Analysis of Rosette Leaf and Heading Traits in Both DH-88 and F2-485

#### **Development of RL in DH-88 and F2-485**

Rosette leaves were phenotyped six times during their development; rosette leaf length and width continuously increase throughout the development (**Figure 3A**). In DH-88, LL increased significantly between 36DAS (days after sowing) and 47DAS, but growth ceased 50DAS as there was no significantly increase between 50DAS and 53DAS (P = 0.146). For LW, growth ceased after 47DAS (P = 0.167). In F2-485, LL increased significantly during all time points measured, while for LW no significant difference was seen between 44DAS and 48DAS (P = 0.763), and we can conclude that leaf growth ceased at 44DAS.

Correlations between LL, LW, and PL at different DAS with the trait "Final heading degree" are shown in **Figure 3B**. In both DH-88 and the F2-485 population, LL and PL are not significantly correlated with heading degree. LW however, is positively correlated with "Final heading degree" in both DH-88 and F2-485. In DH-88, the correlation between LW and heading degree was clear from 36DAS to 53DAS, while this correlation was only significant at 37DAS in F2-485.

#### **Head trait analysis in F2 population**

At the heading stage, the heading degree (HDe) was evaluated for 66 DH-88 lines and 485 F2 plants, while HH, HW, and HWe were only evaluated in the F2-485 population (**Supplementary Table S5**). In the F2 population, the variance was higher for HH than for HW and HWe. The correlation coefficient analyses showed that HH, HW, HWe, and HD were significantly correlated with each other. Correlation between HH and HD was lower, and not significant at the 0.01 level.

#### Detection of QTLs

Quantitative trait locus analyses for both leaf traits at the rosette stage and heading traits were performed for the DH-88 and F2- 485 populations (**Table 2**).

In the rosette stage, QTLs for rosette leaf area length (LL), LW and petiole length (PL) were detected for both the DH-88 and F2-485 population. For LL one QTL was detected in DH-88 (R\_DHLL1), while four QTLs (R\_F2LL2, 3, 4, 5) were detected in the F2-485. R\_F2LL4 and R\_F2LL5 had high LOD scores (9.93 and 9.65), and explained around 18% of the phenotypic variation. There was one LW QTL detected on A03 in DH-88 that explained 33% of the phenotypic variation and one on A04 in the F2 that explained 6.8%. For PL, 6 QTLs were detected in DH-88 on A06. Four QTLs on A06 for PL in F2-485 were identified for different developmental days. The QTL region of A06 showed multiple QTL for LL and PL.

In the heading stage, two QTLs for head height (H\_F2HH1 and HH2), one for head width (H\_F2HW1) and two for head weight (H\_F2HWe1 and H\_F2HWe2) were detected in the F2- 485. For the heading degree trait, a total of five QTL in four LGs in F2-485 and two on two LGs in DH-88 were detected. Of the QTL detected for the HDe trait, H\_F2De7 showed a higher LOD (8.54) value on A08 in F2-485. We found that the HDe QTL region of A03 and A08 were detected in both DH-88 and F2-485, while the HDe QTL on A04 and A05 were only detected in the F2-485 population.

In both DH-88 and F2-485, a major QTL for LT was detected on A06, explaining 64 and 57.8% of the variation.

FIGURE 3 | Sixth rosette leaf traits analysis in DH-88 and F2-485 populations through development. (A) Box plots of the rosette leaf length (LL), leaf width (LW), petiole length (PL), grouped by growth days (the days after sowing: DAS). (B) Correlation analysis between rosette leaf traits and "Final heading degree."

#### TABLE 2 | Quantitative trait loci found in DH-88 and F2-485 population.


#### Co-localization of QTLs Between Rosette Leaf and Head Traits

For the F2-485 population both rosette leaf and head traits were recorded, while for DH-88 population we measured rosette leaf traits and heading degree only. A total of 9 QTLs in rosette stage and 11 QTLs in heading stage were detected in the F2- 485. Several QTLs were found to co-localize (**Figure 4A**). The QTL region of A03, A04, A05, A06, and A08 showed multiple QTLs for two till five traits. Head trait QTL H\_F2HH2 and H\_F2HDe3 were located at the top position of A03 (0.000– 16.071 cM). Heading degree and rosette LW (H\_F2HDe4 and R\_F2LW2) were mapped closely linked on the top of LG A04. The bottom part of A05 (60.708–77.975 cM) showed overlapping QTL for head weight and heading degree traits: H\_F2HW1 and H\_F2HDe5, H\_F2HDe6. QTL for both rosette leaf and head traits R\_F2LL (2, 3, 4, 5), R\_F2PL (9, 10, 11, 12), R\_F2LT2, H\_F2HH1 and H\_F2HWe2 were clustered at the interval of 32.085–62.563 cM region on F2-485 LG A06. QTL for head traits H\_F2HDe7 and H\_F2We1 were detected at the top of A08 around 10.978 cM.

The positions of the rosette leaf RLL, RPL, LT, and heading degree (HDe) QTLs in this study were compared between DH-88 and F2-485 maps. A QTL for rosette leaf length (RLL: 46.8– 49.5 cM) was mapped on DH88 LG A06. QTLs for RLL in F2 (46.3–57.2 cM) were detected on LG A06 below the common marker Bn-A06-p4452778. The SNP marker Bn-A06-p4452778 was linked with the QTL for RPL in both DH-88 (41.5–53.5 cM) and F2-485 (33.4–52.2 cM) in this study. QTLs for LT were mapped on LG A06 and a comparison of map positions could be made using the common Bn-A06-p22910339 marker in the maps. QTL for heading degree (HDe) in DH-88 (5.4–9.5 cM) and F2-485 (0.0–11.0 cM) on the top of LG A08 overlap, as they were linked to the common marker Bn-A08-p4528943 (**Figure 4B**).

### DISCUSSION AND CONCLUSION

Heading is an important breeding trait, as size and shape of the Chinese cabbage heads determine their market share. In addition, the heading trait is an important domestication trait that was domesticated around 1600 years ago (Tan, 1979; Li, 1981). Identification of QTL and the underlying causal genes can both shed light on the mutations and pathways that were selected for during domestication, but will also generate tools for breeders to select for optimal size and shape of Chinese cabbages.

#### Head Formation Related Traits

So far, only few studies are published that study the genetics of leafy head formation (He et al., 2000; Wang et al., 2012; Yu et al., 2013), but only very few combined morphological variation at the rosette stage with variation in heading degree and shape. Mao et al. (2014) demonstrated that variation in expression of TCP4 miRNA (miRNA319) explained head shape and the size, shape and curvature of RL. They showed clear correlation between morphology of leaves in rosette stage and final leafy head shape. In this study, we also extensively phenotyped leaf development during rosette stage and phenotyped the leafy head in a diverse collection of Chinese cabbage hybrids and inbred lines. This allowed us to evaluate variation and correlation between rosette and heading traits. We found that rosette leaf petiole- and leaf

linkage group. Common markers between DH-88 and F2-485 are indicated by "Sxx" and marker information, which is also shown in Supplementary Table S4. (B) Zoom in of QTL hotspots of the rosette traits RLL (green), RPL (yellow), leaf trichomes LT (red) and heading degree HDe (black) on DH-88 and F2-485 on A06 and A08.

blade length were not correlated with final heading degree (HDe). However, rosette leaf width measured at several developmental stages was correlated with HDe during plant development in both F2-485 and DH88 populations. This correlation was weaker at early rosette stage and increased following the head formation process. The observation that rosette leaf width, which largely determines size and heading degree, are correlated fits with the physiology of Chinese cabbage. The special feature of heading Chinese cabbages is that younger leaves that initiate and grow inside the head are not photosynthetically active. As a consequence the outer leaves should remain active until harvest to generate photosynthetates and stimulate N uptake (Wang et al., 2012). Also, the redistribution from older outer leaves to the young inner leaves forming the head is of importance (Lammerts van Bueren and Struik, 2017). At later rosette stages when the head growth increases, the RL are thus very important and their width may affect the growth of the head.

A number of leafy head related traits were studied in the heading stage of the collection of Chinese cabbage accessions. By phenotype variation analysis, we found that in Chinese cabbage the variation in the extent to which the outer heading leaves overlap (HtO) was the largest, while the trait Heading Degree (HDe) showed little variation in the collection. We showed clear correlations among traits, with a cluster of traits related to heading capacity and a cluster of traits related to the shape of the head. The core of a Chinese cabbage head is the compact stem from which the leaves branch off. The shorter the core is, the denser the head. We detected strong correlation between core length (CL) and heading degree (HDe). Surprisingly, head leaf number (HLN) was not significantly correlated to the core length or to other head related traits. In this study the head weight was measured, however, head density (weight divided by volume) was not determined, a trait that is likely related to both leaf number and core length. Currently we developed tools to analyze 3D images that make it possible to calculate head shape, volume and density.

#### Co-localized QTLs

fpls-09-01455 October 1, 2018 Time: 16:45 # 9

By phenotyping a large F2 population derived from a cross between a PC and a CC parent for both rosette and heading traits and in addition phenotyping rosette traits in a DH population generated from the same cross, we found that QTLs for leaf size traits LL, LW, PL, the leaf trichome trait (LT) and leafy head traits HH, HW, HWe, HDe collocated in many cases. This was expected based on the correlations found between these traits.

For rosette leaves traits, common QTLs for leaf length (RLL), petiole length (RPL) and leaf trichome (LT) traits were identified on LG A06 in both the DH-88 and the F2-485 populations. However, several QTLs are found in a single population only. The two populations have very different sizes (485 F2 plants and 65 DH lines), which affects statistical power to detect QTL. In addition, F2 phenotypic data are based on single plants while DH data are based on several biological repeats of the same DH line. This, and the fact that all loci are homozygous in the DH, while in F2 plants homo- and heterozygotes both occur, and dominance relationships of the traits are generally unknown, can all affect QTL detection. QTLs for rosette leaf width (RLW) were identified on chromosome A03 in DH-88 and on chromosome A04 in F2-485. Very interestingly, the rosette leaf width (RLW) QTLs all co-located with heading degree (HDe) QTL. This as discussed based on their correlation, matches the physiological explanation that RL are the source of nutrients for the growing leafy head that serves as a sink. Colocation of QTL for RLL and RPL on LG A06 in both F2-485 and DH-88 population illustrates that total leaf length is correlated to petiole length.

For head related traits in the F2-485, we detected three QTLs for head length (HH) on LG A03, A06, A08, one QTL for head width (HW) on LG A05 and two QTLs for head weight (HWe) on LG A06, A08. Another QTL study in F3 families from a cross between small-and large heading Chinese cabbage lines identified QTLs for head related traits on the same LGs: HH on LG A03 and for HWe on LG A06 and A08 (Ge et al., 2011). In yet another study a head width QTL on LG A05 was detected (Yu et al., 2013). In the study of Yu et al. (2013) HW was phenotyped using 150 recombinant inbred lines from a cross between heading and nonheading Chinese cabbage in both 2011 and 2012 and one HW QTL on LG A05 was identified. In our study, we found that the head traits height, width and weight (HH, HW, and HWe) were

#### REFERENCES


highly correlated with heading degree (HDe) and their QTLs colocated on LG A03, LG A05 and LG A08. Very interestingly, colocation of QTL for RL traits (RLL, RPL) and head related traits (HH, HWe) were identified on LG A06 in F2-485. This again shows the relation between RL and leafy heads, where RL provide the nutrients for the growing leafy head.

In conclusion, we provide information on the relationship between leaf characteristics at rosette stage and leafy head formation. Colocation of QTL for heading related traits and for rosette leaf and heading traits suggests pleiotropic gene regulation of these complex traits. Knowledge about the relationship between the size and shape of RL that affect heading-degree, -size and -shape, is relevant for breeders that select for optimal quality of Chinese cabbage leafy heads.

#### AUTHOR CONTRIBUTIONS

XS and SL designed the experiment. XS performed the data analysis and contributed to manuscript writing. SL arranged the field experiments and collected phenotypic data. LL and XW did F2-485 genotyping. XC and YL helped in the discussion. SS and JZ provided funds for this study. JZ was involved in discussions. GB supervised the whole experiment and wrote the manuscript. All authors approved the final manuscript.

#### FUNDING

Financial support for this work was provided by the Dutch Royal Academy of Sciences China Exchange Program (Grant No. 530-4CDP08), the National Key R&D of China (Grant No. 2017YFD0101802), the National Natural Science Foundation of China (Grant No. 31672151), International Cooperation project in the Science and Technology Support Program of Hebei (Grant No. 17396315D), Hundred Innovative Talent Program of Hebei (Grant No. SLRC2017040), the Natural Science Foundation of Hebei (C2016204170), and the Science and Technology Support Program of Hebei (16226304D-2).

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2018.01455/ full#supplementary-material

responsive patterns associated with leafy head formation in Chinese cabbage. Sci. Rep. 7:42229. doi: 10.1038/srep42229



**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Sun, Luo, Luo, Wang, Chen, Lu, Shen, Zhao and Bonnema. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Differential Gene Expression Caused by the F and M Loci Provides Insight Into Ethylene-Mediated Female Flower Differentiation in Cucumber

Jian Pan, Gang Wang, Haifan Wen, Hui Du, Hongli Lian, Huanle He, Junsong Pan\* and Run Cai\*

School of Agriculture and Biology, Shanghai Jiao Tong University, Shanghai, China

#### Edited by:

Yuke He, Institute of Plant Physiology and Ecology (SIBS-CAS), China

#### Reviewed by:

Xuehao Chen, Yangzhou University, China Nobutaka Mitsuda, National Institute of Advanced Industrial Science and Technology (AIST), Japan

#### \*Correspondence:

Junsong Pan jspan71@sjtu.edu.cn Run Cai cairun@sjtu.edu.cn

#### Specialty section:

This article was submitted to Plant Breeding, a section of the journal Frontiers in Plant Science

Received: 05 April 2018 Accepted: 05 July 2018 Published: 14 August 2018

#### Citation:

Pan J, Wang G, Wen H, Du H, Lian H, He H, Pan J and Cai R (2018) Differential Gene Expression Caused by the F and M Loci Provides Insight Into Ethylene-Mediated Female Flower Differentiation in Cucumber. Front. Plant Sci. 9:1091. doi: 10.3389/fpls.2018.01091 In cucumber (Cucumis sativus L.), the differentiation and development of female flowers are important processes that directly affect the fruit yield and quality. Sex differentiation is mainly controlled by three ethylene synthase genes, F (CsACS1G), M (CsACS2), and A (CsACS11). Thus, ethylene plays a key role in the sex differentiation in cucumber. The "one-hormone hypothesis" posits that F and M regulate the ethylene levels and initiate female flower development in cucumber. Nonetheless, the precise molecular mechanism of this process remains elusive. To investigate the mechanism by which F and M regulate the sex phenotype, three cucumber near-isogenic lines, namely H34 (FFmmAA, hermaphroditic), G12 (FFMMAA, gynoecious), and M12 (ffMMAA, monoecious), with different F and M loci were generated. The transcriptomic analysis of the apical shoots revealed that the expression of the B-class floral homeotic genes, CsPI (Csa4G358770) and CsAP3 (Csa3G865440), was immensely suppressed in G12 (100% female flowers) but highly expressed in M12 (∼90% male flowers). In contrast, CAG2 (Csa1G467100), which is an AG-like C-class floral homeotic gene, was specifically highly expressed in G12. Thus, the initiation of female flowers is likely to be caused by the downregulation of B-class and upregulation of C-class genes by ethylene production in the floral primordium. Additionally, CsERF31, which was highly expressed in G12, showed temporal and spatial expression patterns similar to those of M and responded to the ethylene-related chemical treatments. The biochemical experiments further demonstrated that CsERF31 could directly bind the promoter of M and promote its expression. Thus, CsERF31 responded to the ethylene signal derived from F and mediated the positive feedback regulation of ethylene by activating M expression, which offers an extended "one-hormone hypothesis" of sex differentiation in cucumber.

Keywords: cucumber, ethylene response, sex differentiation, unisexual flower, floral development

# INTRODUCTION

Cucumber (Cucumis sativus L.), a horticultural crop consumed worldwide, has the third highest production quantity (61 million tons in 2016) after tomato and onion, and China leads in the production of cucumber with 76% of all production (Un Food and Agriculture Organization Corporate Statistical Database [Faostat], 2017). The cucumber is more than an agronomically

important vegetable, and due to the diversity of its flower sexual types and plant sexual systems (flower sexual types that express on a single plant), the cucumber is an ideal model for investigating the mechanism of sex differentiation in unisex flowers (Malepszy and Niemirowicz-Szczytt, 1991). The sexual types of the cucumber flower include the female flower, male flower, and bisexual flower. According to these sexual types, the following five cucumber sexual systems have been categorized: monoecy (presence of female and male flowers), gynoecy (only female flowers), androecy (only male flowers), andromonoecy (male and bisexual flowers), and hermaphrodite (only bisexual flowers) (Galun, 1962; Kubicki, 1969a,b,c; Fruchard and Marais, 2017). Notably, the variety seen in the cucumber sexual systems is rare in the plant kingdom, covering ∼95% of sexual systems among angiosperms (Käfer et al., 2017). Treatment with exogenous ethylene or ethylene-related chemicals can affect the sex of the flowers in cucumber plants (MacMurray and Miller, 1968; Iwahori et al., 1969). The ethylene-mediated differentiation of unisexual flowers in cucumber has been explained by the "one-hormone hypothesis" (Yin and Quinn, 1995). This hypothesis posits that ethylene in cucumber inhibits maleness and induces femaleness by modulating the expression levels of F and M (Yin and Quinn, 1995; Trebitsh et al., 1997). According to this hypothesis, the F locus encodes a gene for ethylene production, and the M locus encodes an ethylene-response factor that might be an ethylene receptor. This hypothesis has been modified and improved based on the map-based cloning of M, and an M-mediated positive feedback regulatory mechanism of ethylene has been proposed (Li et al., 2008, 2009, 2012). In our previous study, sex expression was associated with the fruit shape (Tan et al., 2015). In most M mutant lines (mm genotype), the gynoecium presents shorter ovaries and oval/round fruit, whereas the MM genotype presents elongated ovaries and normal long fruit, suggesting that M promotes an ovary of better quality via a positive feedback mechanism. The expression analysis also demonstrated that the M transcript began to accumulate beneath the pistil primordia of the flower buds during the bisexual stage (Saito et al., 2007).

Many studies have attempted to identify the "sex genes" in cucumber, and three Mendelian loci – F/f, M/m, and A/a, have been identified to be responsible for the different sexual systems (Pierce and Wehner, 1990). Using map-based cloning, the "sex genes" (F, M, and A) have been cloned from the aforementioned loci (Mibus and Tatlioglu, 2004; Li et al., 2009; Boualem et al., 2015). Unsurprisingly, F(CsACS1G), M(CsACS2), and A(CsACS11) all encode 1-aminocyclopropane-1-carboxylate synthase (ACS), which catalyzes the rate-limiting step in ethylene biosynthesis; the association between ACS and ethylene has been known for decades.

The melon (C. melo) serves as another well-known model in sex differentiation studies. In the melon, the sex-determination pathway model integrating CmACS11 (i.e., A in both cucumber and melon), CmWIP1 (i.e., G in melon and CsWIP1 in cucumber), and CmACS7 (i.e., M in both melon and cucumber) explains the primary mechanism of male, female, and bisexual flower generation (Boualem et al., 2015). This model can be similarly applied to the cucumber according to the current study. Briefly, female flowers are initiated due to the expression of A, which represses the expression of CsWIP1. Thus, the nonexpression of CsWIP1 promotes the expression of M, which inhibits stamen development. If a non-functional M is expressed, hermaphroditic flowers develop instead of female flowers. The loss of function of A generates androecious plants, leading to the expression of CsWIP1 at the whole plant level. Notably, F, as a gain-of-function structural variation (Zhang et al., 2015), was not mentioned in this model. However, these findings have triggered further questions regarding the mechanism by which ethylene (produced by F, M, and A) controls sex expression; furthermore, the mechanism by which the ACS genes were recruited to form a regulatory mechanism in cucumber female flower evolution remains unclear.

In addition to F, M, and A, other ethylene biosynthetic genes, such as CsACO2 and CsACO3, which encode ACC oxidases, are involved in sex expression. However, the transcript levels of CsACO2 and CsACO3 in the shoot apices are inversely correlated with femaleness, indicating the potential presence of feedback inhibition that controls ethylene production (Kahana et al., 1999). In addition, CsACO2 can influence floral organ development in both cucumbers and Arabidopsis (Duan et al., 2008; Chen et al., 2016). CsWIP1 can directly bind the promoter of CsACO2 to repress its expression, suggesting that CsWIP1 plays the same key role as G (CmWIP1) in sex differentiation (Chen et al., 2016). An ethylene receptor gene, CsETR1, has been shown to participate in stamen arrest via the induction of DNA damage (Wang et al., 2010). Further studies suggest that a cucumber nuclease-encoding gene, CsCaN, responds to the ethylene signal and allows the DNA of the stamen primordia to damage the developing female flowers (Gu et al., 2011).

Previous studies have provided information and theories suggesting that hormones, particularly ethylene, promote femaleness in cucumbers. However, molecular and biochemical data are scarce, limiting our understanding of the sex determination mechanism in cucumber and preventing our ability to artificially control the quantity and quality of female flowers. To investigate the ethylene-mediated sex expression process and the interplay mechanism between F and M, three cucumber near-isogenic lines (NILs), namely H34, G12, and M12, with different F and M loci were produced by backcrossing. We found that several novel B- and C-class floral homeotic genes were differentially expressed between the NILs by comparing the transcriptomes of the shoot apices. ERF transcription factors, which act as transacting regulators in the final step of the ethylene signaling pathway, bind to conserved elements, such as the GCC-box or ERE-box, in the promoters of many ethylene-inducible genes (Ohme-Takagi and Shinshi, 1995; Solano et al., 1998; Han et al., 2016). An ERF family gene, CsERF31, was identified among the differentially expressed genes (DEGs) because its spatiotemporal expression correlated with M. Furthermore, we demonstrated that CsERF31 could directly bind the promoter of M to promote its expression. In this study, we propose an extended "onehormone hypothesis" of sex differentiation and development of cucumber.

# MATERIALS AND METHODS

fpls-09-01091 August 11, 2018 Time: 19:32 # 3

# Plant Materials and Growth Conditions

To explore the genes and gene networks that control the sex expression in cucumbers, we performed RNA-Seq analyses of shoot apices from three NILs – G12 (FFMMAA), M12 (ffMMAA), and H34 (FFmmAA). The G12 was substituted for M loci from H34 via seven generations of backcrosses by our marker-assisted selection program (Li et al., 2009). The M12 was a spontaneous mutant from G12 with a loss of F (CsACS1G), and has been stabilized via five generations of self-crosses. The NILs belong to the European greenhouse cultivar and have short gloss fruits covered with small spines without warts (**Figures 1A4–C4**). No visible difference was observed among the NILs in terms of the plant morphology, except that G12 and H34 express female and bisexual flowers, respectively, at the whole plant level, while M12 bears only male flowers before the 25th ± 3 node in the mainstem (**Figure 1B4**); however, the female flowers can be observed in the lateral branches. As the F and M loci cause significant differences prior to the morphological changes, we chose shoot apices that contained 1- to 5-mm buds for the transcriptome analyses.

The seeds of the NILs were germinated on a wet filter paper in a Petri dish at 28◦C in the dark overnight. This was followed by growing the resulting seedlings in an artificial climate room for 16 h/8 h at 25◦C/18◦C during the day/night. The samples used for the gene expression analysis were collected from the plants in the artificial climate room. The chemically treated and control plants were transferred to a greenhouse at the Shanghai Jiao Tong University in April 2017.

#### Preparation of RNA Samples for RNA-Seq and Quantitative Real-Time RT-PCR (qRT-PCR)

The shoot apices from the aforementioned cucumber plants were collected at the four-leaf stage. The samples were immediately frozen in liquid nitrogen and stored at −80◦C. The total RNA was isolated using an RNA extraction kit (Epigenetics, United States). The RNA was checked for contamination using RNase-free agarose gel electrophoresis, and then, the RNA purity was examined using a NanoDrop Spectrophotometer (THERMO, United States). The RNA integrity was measured and assessed using an Agilent 2100 Bioanalyzer system (Agilent, United States).

For the qRT-PCR experiment, G12 flower buds with different lengths of 1, 2, 5, and 10 mm were collected. A day before flowering, the petal, stamen, stigma, mesocarp, and exocarp were isolated from the G12 female flowers. Shoot apices treated with chemicals were isolated from the chemically treated and control plants as described below.

# Transcriptome Library Construction and Sequencing

The RNA-Seq for comparative transcriptomic analyses of three genotypes were performed with three biological replicates. The library construction and sequencing were performed using a

Bars = 4 cm. P, pistil; St, stamen.

BGISEQ-500 by Beijing Genomic Institution (BGI, China). The genomic DNA was removed with two digestions using Amplification grade DNase I (Epigenetics, United States). The RNA was sheared and reverse transcribed using random primers to obtain the cDNA, which was used for the library construction. The library quality was determined using a Bioanalyzer 2100 (Agilent), and then the library was used for sequencing using the sequencing platform BGISEQ-500 (BGI, China) (Fehlmann et al., 2016). All the generated raw sequencing reads were filtered to remove reads with adaptors and reads in which unknown bases were greater than 10% of low-quality reads. The clean

reads were obtained and saved in the FASTQ format. The clean data have been uploaded to the European Nucleotide Archive (ENA) (Project ID: PRJEB25272, Sample ID: SAMEA104656171, SAMEA104656172, SAMEA104656173, SAMEA104656174, SAMEA104656175, SAMEA104656176, SAMEA3469513, SAMEA3469514, SAMEA3469515).

#### Bioinformatics Analysis of RNA-Seq Data

We used Bowtie2 (Langmead et al., 2009) to map the clean reads to the reference genome, (Cucumber\_ChineseLong\_v2<sup>1</sup> , Huang et al., 2009). The read counts were summarized, and the Fragments Per Kilobase of exon per Million fragments mapped (FPKM) was calculated for each annotation on the reference sequence. The NOISeq method (Tarazona et al., 2011) was used to screen for DEGs between the groups. In this study, an absolute value of log2foldchange ≥ 0.6 and diverge probability ≥ 0.6 were used as cut-offs to screen for gradual expression differences among the three NILs (**Supplementary Tables S2, S3**).

A Venn diagram analysis was performed using Venny tools (Oliveros, 2007–2015). An expression trend analysis was performed using OmicShare<sup>2</sup> tools. To identify the homolog in Arabidopsis, we used the cucumber protein IDs (Huang et al., 2009) to batch query the Arabidopsis proteins (TAIR10) using BLASTP on the Cucurbit Genomics Database<sup>1</sup> with an e-value cut-off of 1e−<sup>1</sup> . Based on the cucumber homolog in Arabidopsis, the protein interaction network was predicted using STRING<sup>3</sup> , which is a database that aims to provide a critical assessment and integration of protein–protein interactions, including direct (physical) and indirect (functional) associations. The prediction of protein–protein interaction was measured by the combined score, which was computed by combining the probabilities from the different evidence channels and corrected for the probability of randomly observing an interaction.

#### Chemical Treatments

The ethylene-related chemicals used in this study included ethephon (an ethylene-releasing agent), aminoethoxy vinyl glycine (AVG, ethylene biosynthesis inhibitor), and AgNO<sup>3</sup> (ethylene action inhibitor). These three chemicals were purchased from Sigma (China). All cucumber plant treatments were performed according to our previous studies (Li et al., 2012). Each treatment employed ten plants, and five plant shoot apices were excised for RNA extraction after the treatment. To observe the effects of the chemical application on the sexual expression of cucumber plants, the other five plants of each treated line were allowed to grow to maturity. The sexual types were recorded up to the 25th node on the main stem.

#### Gene Expression Analyses by qRT-PCR

The cDNA used for the qRT-PCR was reverse transcribed from 1 mg of total RNA using a QuantScript RT Kit (TIANGEN, China). The qRT-PCR analyses were carried out using FastStart Essential DNA Green Master (Roche, United States) on a

<sup>1</sup>http://cucurbitgenomics.org

<sup>2</sup>http://www.omicshare.com/tools

Rotor-Gene Q Real-Time System (QIAGEN, United States). The gene-specific primers used for the qRT-PCR were designed using Primer 3 software. The cucumber actin was used as an internal control to normalize the expression data (Li et al., 2012). Five biological replicates for chemical treatments and three biological replicates for other samples were used per gene. Each qRT-PCR experiment was performed with three technical replicates. The gene-specific primers are listed in **Supplementary Table S5**.

All data were expressed as the mean values ± standard deviation (SD) of biological replicates. Differences were analyzed with one-way ANOVA using SAS 9.1.3 software. The P-values < 0.05 were considered to be significant.

#### Yeast One-Hybrid Assays

For the yeast one-hybrid assay, the open reading frames (ORFs) of CsERF31 were amplified using PrimerStar GXL DNA polymerase (TAKARA, Japan) and then cloned into pB42AD (EcoRI/XhoI) using a ClonExpress II One Step Cloning Kit (Vazyme, China). The two ERE-boxes in the M promoter (18 bp length containing the ERE-box element, with 3 tandem repeats) were inserted into the pLacZ vector (**Figure 5C**). The analyses were performed using the yeast EGY48a strain as previously described, and empty vectors were used as negative controls. The transformants were cultivated on SD/-Leu/-Ura medium and tested on the SD/-Leu/-Ura medium with X-gal (5-Bromo-4 chloro-3-indolyl β-D-galactopyranoside) (Yan et al., 2017).

#### Dual-Luciferase (Dual-LUC) Assay

This assay was performed as previously described (Luo et al., 2014). Briefly, the ORFs of CsERF31 were cloned into the pHB vector to be driven by a 35S promoter for overexpression. This was followed by the recombinant vectors being separately transformed into Agrobacterium tumefaciens GV3101 as effectors. The 2 and 1-kb promoters of M were cloned into the pGREEN0800 vector to drive the firefly luciferase (LUC) reporter gene. Each M promoter vector was then co-transformed with the helper plasmid pSoup19 into Agrobacterium GV3101 as reporters. The reporter and effector were mixed at a 1:2 volume ratio and injected into tobacco (Nicotiana benthamiana) leaves. An empty pHB vector was used as a negative control. The constitutive 35S promoter driving Renilla luciferase (REN) was used as an internal reference. The leaf samples used in the Dual-LUC assay were collected after 2 days using commercial Dual-Luciferase reaction reagents according to the manufacturer's instructions (Promega, United States). Four biological replicates of each sample were measured using a GloMax 20/20 Luminometer (Promega, United States). The primers used in the yeast one-hybrid and Dual-LUC assay are listed in **Supplementary Table S5**.

# RESULTS

#### Comparison of Gene Expression Among Three NILs With Different Genotypes

The RNA-Seq sequencing generated ∼24 million single-end reads for each sample, and three biological replicates were performed

<sup>3</sup>https://string-db.org

for each line (**Supplementary Table S1**). By performing a bioinformatics analysis of the RNA-Seq data, we identified 631 DEGs (**Supplementary Table S2**), including 281 genes that were upregulated and 350 genes that were downregulated in G12 compared to those in H34; in addition, we identified 534 DEGs (**Supplementary Table S3**), including 327 genes that were upregulated and 207 genes that were downregulated in G12 compared to those in M12.

To verify the DEGs identified by RNA-Seq, we performed qRT-PCR assays using independently collected samples that were at the same developmental stage as those used for the RNA-Seq analysis. A total of 20 predicted sex differentiation-related genes were selected from the DEGs. All the 20 genes showed the same expression patterns in the qRT-PCR assays as the RNA-Seq data (**Figure 3** and **Supplementary Figure S1**). The Pearson's correlation coefficients between the qRT-PCR and RNA-Seq data were H34 = 0.92904 and M12 = 0.99016 (P < 0.001) when compared with G12, respectively, indicating that the RNA-Seq data were highly reliable.

#### Phytohormone and Floral Regulators Co-regulate With M in Ethylene-Mediated Female Development

We performed a series of analyses to screen for M-related genes among the DEGs. Since M was activated by the ethylene produced by F (Li et al., 2012), the DEGs between G12 and M12 contained genes which activated M and were involved in the M positive-feedback. However, the DEGs between G12 and H34 only contained genes which were involved in the M positive-feedback, because M (in G12) and m (a loss-offunction M, in H34) are both activated by F but m could not initiate the positive-feedback. A total of 264 genes filtered by the Venn and Expression trend analyses were selected (**Figures 2A– C** and **Supplementary Table S4**), and the genes that were most homologous to Arabidopsis genes were identified. We then performed a STRING analysis (Szklarczyk et al., 2015) to predict the association networks of the homologs. The network showed that 11 homologs that were associated with ACS7 (homolog of M) (**Figure 2D**) and 7 homologs which had the opposite expression trends with M were related with floral organs development (**Figure 2E**).

In the predicted protein interaction network, several phytohormone and floral regulators were linked to ACS7. For example, JAZ1 and JAZ8 have a protein-protein interaction (Chini et al., 2009) and regulate ERF1 via EIN3/EIL1 in Arabidopsis (Song et al., 2014) (**Figure 2D**). Accordingly, SUP, which is a well-known negative transcriptional regulator of B-class floral genes, and SHP2 (i.e., CAG2, AG-like C-class floral gene) were indirectly linked to ACS7. Notably, ERF1, which is a homolog of CsERF31, belongs to the EREBP (ethylene-response element binding protein) subfamily, which mediates the response to ethylene, suggesting that CsERF1 may play a key role in the ethylene response pathway of F and M (**Figure 2D**). Among the homologs which had the opposite expression trends with M, several well-study genes, such as AP3 (CsAP3), WUS (CsWUS), and SPL (CsSPL) were screened out, suggesting that these genes had interactions with B-class genes during floral developments in cucumber (**Figure 2E**).

### CsERF31 and M Had Similar Expression Patterns and Responses to the Ethylene-Related Chemical Treatments

To further explore the role of CsERFs in sex expression, we verified the expression patterns of seven ERFs of DEGs by qRT-PCR. Only CsERF31 showed an expression pattern similar to that of M (**Figure 3**). Furthermore, we tested the spatiotemporal expression of CsERF31 and M in G12 flowers. The expression levels of CsERF31 and M were very low in different organs of female flowers that were ready to blossom. However, CsERF31 and M were expressed at a high level in the 1- and 2-mm buds (i.e., at ∼stages 6 and 10, respectively), while the expression levels were reduced in both the 5- and 10-mm buds/ovaries (**Figure 4A**).

The ethylene-related chemical treatments modulated the sexual types of the cucumber flower buds. The treatment of the shoot apices with ethephon generated female flowers at low nodes in M12, whereas no effect was observed in G12. The treatments with AVG and AgNO<sup>3</sup> generated male flowers in G12, whereas no effect was observed in M12 (**Figure 4B** and **Supplementary Figure S2**). Accordingly, the expression levels of CsERF31 and M were significantly upregulated in ethephon-treated G12 and M12. However, the expression levels were downregulated only in G12 treated with AVG or AgNO<sup>3</sup> (**Figure 4C**). In conclusion, CsERF31 and M showed a highly synchronous expression pattern among a variety of samples.

### CsERF31 Bind the ERE-Box in the M Promoter and Activate Its Expression

To obtain additional evidence of the activation effect of CsERF31 on the M expression, we conducted heterologous transient expression experiments in tobacco. The −2 and −1 kb sequences upstream of the translation start site of M were fused to sequence encoding LUC to construct the reporters proM2kb:LUC and proM1kb:LUC, respectively (**Figure 5A**). In addition, a construct containing the coding sequence of CsERF31 downstream of the constitutive 35S promoter was generated. The two constructs were simultaneously introduced into tobacco leaves. Compared to the control (empty vector + proM2kb:LUC), the transcriptional activation activity (LUC/REN) of proM2kb:LUC, but not proM1kb:LUC, was significantly promoted by CsERF31 (**Figure 5B**). We further investigated whether CsERF31 directly activated the M expression. The yeast one-hybrid assay revealed that CsERF31 was bound to the distant ERE-box (D EREbox, −1334) but not the near ERE-box (N ERE-box, −154) of the M promoter (**Figure 5C**). In conclusion, CsERF31 could bind to the ERE-box in −1334 upstream of M to activate its transcription but could not bind to the near ERE-box −154 upstream of M (**Figure 5D**). In addition, the ORFs sequence of CsERF31 and the 2-kb promotor sequence of M are consistent among the NILs (unpublished data).

genotypes. Expression trend analysis of DEGs was performed to screen genes that have downregulated expression trends (Trend 1, 2, and 3) similar to M (C, left) and upregulated trends (Trend 4 and 5) that are opposite. A total of 264 genes were classified into 5 trends and the best homolog was identified in Arabidopsis (C) (Supplementary Table S4). The predicted protein interaction network (D) showed the homologs that were associated with ACS7, homolog of M (D, left red box). The interaction network (E) showed the homologs that were associated with B-class floral protein (AP3 and PI) and shared opposite expression trend with M. In the network, the links between proteins signify the various interaction data supporting the network, colored by evidence type (E, right inset). Information of genes presented in (D) and (E) are listed in Table 1.

# DISCUSSION

In general, the wild and semi-wild cucumber varieties (e.g., hardwickii and xishuangbanna) are monoecious (ffMMAA), indicating that the sexual types of flower buds in the shoot

differences (P < 0.05) of 7 ERF family genes and M in G12, H34, and M12.

apices are indeterminate. However, a different result is observed if the genotype is FFMMAA or FFmmAA because the sex of the buds could be predetermined at the whole plant level. Previous study indicates that the sex expression is directly related to ethylene release rate in different genotype

(P < 0.05) of CsERF31 and M in different treatments.

cucumbers (e.g., FFMMAA, ffMMAA, and FFmmAA) (Cui et al., 2016). Besides, in our previous study, the ethylene release rate and M expression both show the gradually decreasing trend in FFMMAA, FFmmAA, and ffMMAA cucumber lines (Li et al., 2012). Thus, we propose that NILs with three genotypes could be ideal materials to investigate the process of female flower differentiation, considering both F and M. The main purpose of this study was to provide valuable data and clues for exploring the molecular mechanism of cucumber unisexual development, particularly female flower differentiation. In this study, we identified that CsERF31 acts as a regulator of M. Although the function of CsERF31 requires further investigation, our data should be further explored to identify the relevant genes and their networks in the process of female flower differentiation in the cucumber.

### The Unisexual Flower in Cucumber May Be Characterized by the Special Expression Patterns of B- and C-Class Genes

The identification of ABC genes that determine the floral organ identities (Coen and Meyerowitz, 1991) was considered

FIGURE 5 | Effects of CsERF31 on the transcriptional activity of M promoter. (A) Schematic diagram of the constructs using Dual-LUC analysis in a tobacco transient expression system. Effector constructs contained the CaMV 35S promoter fused to the transcription factors CsERF31. Reporter constructs contained 2 or 1 kb of M promoter, upstream of the translation initiation sites, fused to the LUC reporter gene. REN was used as the internal control. (B) Effects of CsERF31 on the activity of 2 kb (two ERE-box in it) or 1 kb (only the near ERE-box in it) of M promoter. Error bars represent the SD from four biological replicates, and different letters (a, b) indicate significant differences (P < 0.05) of the transcriptional activation activity with different promoter regions. (C) Yeast one-hybrid assay of interaction between CsERF31 and the distant ERE-box (D ERE-box) or the near ERE-box (N ERE-box). (D) Schematic diagram of CsERF31 binding the D ERE-box in M promoter to activate M transcription.


TABLE 1 | Information of the genes present in the predicted protein interaction network.

a breakthrough in the understanding of sex differentiation (Bai and Xu, 2012). In this study, as **Figure 2** and **Table 1** depict, the expression of the B-class genes CsPI (putative ortholog of PISTILLATA in Arabidopsis) and CsAP3 was markedly suppressed in G12, which was consistent with our predictions, suggesting that the B-class genes are functionally conservative in cucumber. Moreover, the CsSUP (putative ortholog of SUPERMAN/FLO10 in Arabidopsis) transcripts were only increased in G12, showing an expression pattern opposite to that of CsPI and CsAP3. The SUP/FLO10 have been well studied and are considered cadastral genes that act indirectly to prevent B-class genes from acting in the gynoecial whorl in Arabidopsis (Schultz et al., 1991; Bowman et al., 1992; Gaiser et al., 1995; Jacobsen and Meyerowitz, 1997). Thus, the reduced transcripts of CsPI and CsAP3 might be caused by the upregulation of CsSUP in G12. In addition, CsWUS was recently reported to play a key role in promoting flower development by activating CsAP3 and CUM1 (Zhao et al., 2017), which showed 10-fold

higher expression in M12 and H34 when compared to that in G12. Therefore, we speculated that CsWUS and CsSUP act as regulators that stimulate and suppress the transcription of CsAP3, respectively, in different genotypes to control sex expression. Furthermore, previous studies reveal that CsAP3 also performs unique functions in regulating cucumber flower development by activating the CsETR1 expression (Sun et al., 2016), thus establishing a direct connection between ethylene and B-class genes.

Interestingly, the transcripts of a C-class AG homolog, CAG2, which did not appear to be modulated by gibberellin or ethylene in the previous studies (Kahana et al., 1999), presented a gradient reduction in G12, H34, and M12. Previous studies have also demonstrated that CAG2 was specifically expressed in the ovary (Kahana et al., 1999), which is an organ that develops relatively late in the female flower. Hence, we speculated that CAG2 may not directly respond to ethylene but is likely to play an essential role in ovary development after the female flower is initiated. The insufficient expression level of CAG2 in H34 helps explain the reason for its short fruit. Altogether, we suggest that both B- and C-class genes are involved in sex expression and are selectively suppressed or promoted by endogenous ethylene to form male or female flowers, respectively. However, evidence supporting the causality between ethylene and B- and C-class genes must be provided in further studies.

#### Investigation of the Positive Feedback Regulation of M Increases Our Understanding of Cucumber Unisexual Flower Development and Evolution

The feedback regulation of M renders the "one-hormone hypothesis" more powerful in explaining the sex expression of cucumbers. Meanwhile, the amplified ethylene signal triggered further questions regarding why and how this process evolved. To study the positive feedback of M, we performed a co-expression trend analysis to screen for M-related genes among the DEGs. Consequently, CsERF31 was identified with the characteristics of co-expression (**Figure 4A**), co-responding (**Figure 4C**), and directly interacting (**Figure 5**) with M. Notably, CsERF31 and M were specifically expressed in 1- to 2-mm flower buds and faded in 5- to 10-mm flower buds in G12 (**Figure 4A**). Thus, we speculate that the expression of CsERF31 and M may be activated at the initial stages of flower buds and then immediately suppressed via a cucumber unique mechanism. The highly synchronous expression pattern implied a strong correlation between CsERF31 and M, which led us to verify the direct transcriptional activation of M by CsERF31. Interestingly, CsERF31 could activate the M expression only by binding to the distant, not the near, ERE-box from the M ORF (**Figure 5**). Thus, other ERF family genes are likely to bind to the near EREbox to regulate M. In this study, CsERF31 and M presented different expression patterns with six other ERF genes (**Figure 3**). Thus, it is likely that the other six ERFs are not involved in the positive feedback regulation or serve as negative regulators in ethylene signaling (Han et al., 2016). However, the regulation

mechanisms of M remain obscure and should be explored in further studies.

Combining the understanding from previous data regarding F and M, the following extended model of sex differentiation in cucumber is proposed (**Figure 6**):


# Transcriptome Analysis Suggests That the Female Reproductive Organs May Be Regulated via Ethylene-Auxin Crosstalk

Successful fertilization in plants requires the properly coordinated development of the female reproductive organs. In this study, as the DEGs depict in **Supplementary Tables S2, S3**, the predicted cucumber HECATE (HEC) (Csa2G285890) showed 45-fold higher expression in G12 than in M12 and H34. The HEC encodes a basic-helix-loop-helix (bHLH) protein that can dimerize with SPTTULA (SPT), which is involved in the auxin-mediated control of gynoecium patterning in Arabidopsis (Gremski et al., 2007). In addition, previous studies have demonstrated that SPT coordinates with additional bHLH members, such as CRABS CLAW (CRC), to regulate carpel and fruit development by regulating auxin distribution in Arabidopsis (Alvarez and Smyth, 1999; Heisler et al., 2001; Girin et al., 2011). Unsurprisingly, the putative cucumber CRC showed a 3.8-fold increase in G12 compared to that in H34 (**Figure 2D**). These similar expression trends suggest that the regulatory mechanism of floral organogenesis may be conserved between cucumber and Arabidopsis. Thus, we speculate that the crosstalk network between ethylene and auxin may be involved in modulating the floral organogenesis in cucumbers.

The Cucumber fruit is a type of fleshy fruit that is harvested while still immature. The ovaries and fruits of most cucumber varieties with an M locus are elongated, whereas those of hermaphroditic plants with m alleles are oval/round shaped (Tan et al., 2015). A recent study has shown that the AUX/IAA genes SlIAA29 can respond to ethylene signaling, and this study revealed a strong correlation between the auxinand ethylene-related genes, suggesting a significant crosstalk between auxin and ethylene during tomato ripening (Zhang, 2017). According to our data, the transcripts of the indole-3-acetic acid (IAA)-inducible gene CsIAA29 (Csa3G877650) were upregulated in G12. Furthermore, the putative cucumber histidine phosphotransfer protein (AHP) (Csa2G285890), which is a positive regulator of cytokinin signaling in Arabidopsis (Hutchison et al., 2006), displayed a 3.3-fold higher expression in G12. Besides, a previous study has shown that AHP can interact with the ethylene receptor ETR1 (Scharein et al., 2008), suggesting that the ethylene-cytokinin interaction is likely to exist in cucumber fruit development. Similarly, the putative cucumber FRUITFULL (FUL) (Csa1G039910), which is a MADS-box gene regulating cell differentiation during fruit development in Arabidopsis (Gu et al., 1998) and fruit ripening in tomato (Fujisawa et al., 2014), showed a twofold reduction in H34. Notably, the aforementioned genes, CsIAA29, CsCRC, and CsFUL, displayed similar expression trends in a previous study investigating the fruit length in cucumbers (Jiang et al., 2015), thus suggesting that these genes may play pivotal roles in early fruit development. Interestingly, CsSPL, which encodes a protein that directly interacts with CsWUS (WUSCHEL) and positively regulates CsWUS and CsARF3 (AUXIN RESPONSE FACTOR3) expression (Liu et al., 2018), showed a 4.4-fold higher expression in M12 and H34. Thus, CsSPL not only regulates sex organ development by interacting with auxin signaling, as recently reported, but may also participate in the process of sex differentiation in cucumbers.

#### Exploring Gene Function and Regulatory Networks of Cucumber Sex Differentiation in Model Plants, Are Equally Important

Previous studies have demonstrated that the organ-specific downregulation of the ethylene receptor gene (CsETR1) or upregulation of the ethylene synthesis gene (CsACO2) in transgenic Arabidopsis plants can generate female flowers (Duan et al., 2008; Wang et al., 2010). Subsequent studies have demonstrated that CsACO2 is essential for female flower initiation and is regulated by CsWIP1 in cucumbers (Chen et al., 2016). Previous studies have also shown that At3G23240 (AtERF1) is upregulated in CsETR1 downregulated transgenic Arabidopsis, and "cucumber homolog of At3G23240" has a highly specific expression pattern at stage 6 in female flowers (Wang et al., 2010). In our study, CsERF31 was also a "cucumber homolog of At3G23240," and its expression was consistently increased significantly at stage 6 in the female flowers (∼1 mm flower buds of G12) (**Figure 4A**) to promote the transcription of M (**Figure 6**). Additionally, previous studies investigating CsAP3 have shown a conservative function by rescuing the Arabidopsis ap3 mutant, but this is a novel function in cucumber by activating CsETR1 transcription (Sun et al., 2016). Besides, CUM1 was identified in petunia (Petunia hybrida) with the function of inducing reproductive organ fate (Kater et al., 1998), and the function is supported by in situ hybridization in a recent study (Ran et al., 2018). Overall, these findings in cucumber and Arabidopsis or other species complement and verify each other and have inspired the understanding of sex differentiation in plants in a broader perspective.

Our main purpose for studying cucumber sexual expression was to understand the role played by ethylene in this process, particularly the characteristics that distinguish the ethylene regulatory mechanism in cucumbers from that in other plants. Therefore, we consider that the functional studies investigating cucumber sex differentiation-related genes in well-studied plants, such as Arabidopsis, are efficient and provide practical approaches to discovering the regulatory mechanism from a multispecies evolutionary perspective in angiosperm. In conclusion, cucumber sex differentiation should be studied and understood as an evolutionary event more than traditional floral development.

# AUTHOR CONTRIBUTIONS

JSP and RC designed the experiments. JP and HW performed the experiments. GW, HL, and JP analyzed the data. JP wrote the paper along with HD. HH planted the cucumber plants.

# FUNDING

This work was supported by the National Natural Science Foundation of China (Grant Nos. 31772308, 31272185, and 31471879).

#### ACKNOWLEDGMENTS

fpls-09-01091 August 11, 2018 Time: 19:32 # 11

We thank Dr. Hongquan Yang (Fudan University) for providing the pHB vector, and Kexuan Tang (Shanghai Jiao Tong University) for providing pB42AD and pLacZ vector.

#### REFERENCES


#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2018.01091/ full#supplementary-material



**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Pan, Wang, Wen, Du, Lian, He, Pan and Cai. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

fpls-09-01091 August 11, 2018 Time: 19:32 # 12

# Functional Dissection of Auxin Response Factors in Regulating Tomato Leaf Shape Development

Lingjie Wu, Zhendong Tian and Junhong Zhang\*

Key Laboratory of Horticultural Plant Biology, Ministry of Education, Huazhong Agricultural University, Wuhan, China

The phytohormone auxin is involved in many aspects of plant growth and developmental processes. The tomato Aux/IAA transcription factor SlIAA9/ENTIRE/E plays an important role in leaf morphogenesis and fruit development, and the E gene encodes a protein from the Aux/IAA family of auxin response repressors. Both SlIAA9-RNAi transgenic and entire (e) mutant plants reduce the leaf complexity in tomato, but the underlying mechanism is not yet completely resolved. Auxin signaling is known to regulate target genes expression via Aux/IAA and ARFs (auxin response factors) transcriptional regulators. ARFs mediate a wide range of developmental processes. Through an Y2H (yeast two-hybrid) assay coupled with expression profiling of the SlARF genes family, we identified a group of ARFs: SlARF6A, SlARF8A, SlARF8B, and SlARF24. Pull-down and BiFC (Bimolecular Fluorescence Complementation) results demonstrated that these SlARFs interact with SlIAA9 in vitro and in vivo, and the e mutation altered the expression patterns of multiple SlARFs. The simple leaves of the e mutant were partially converted to wild-type compound leaves by VIGS (virus-induced gene silencing) of these four SlARFs. Furthermore, IAA content in these samples was significantly increased compared to the e mutant. In addition, SlARF6A and SlARF24 bound to the SlPIN1 promoter and act as transcriptional activators to regulate genes expression involved in leaflet initiation. It may also suggest that SlARFs regulate leaf morphology through direct binding to auxin-responsive genes in the absence of SlIAA9, providing an insight for the role of SlARFs in leaf shape development.

#### Keywords: auxin, development, leaf shape, SlARFs, SlIAA9, tomato

# HIGHLIGHTS

We firstly found that SlARF6A, SlARF8A, SlARF8B, and SlARF24 could regulate tomato leaf development in a redundant manner; Furthermore, SlARF6A and SlARF24 bound to the SlPIN1 promoter to regulate genes expression involved in leaflet initiation.

# INTRODUCTION

Leaves are one of the main organs of flowering plants, and exhibit a tremendous diversity in shape and size. The shape of leaves varies enormously within the same species and individual plants, and can be ascribed to ranging from simple to compound. Variation is one of the most conspicuous aspects of plant diversity in leaf shape. This diversity is often achieved by the adjustment of leaf

#### Edited by:

Han Xiao, Shanghai Institutes for Biological Sciences (CAS), China

#### Reviewed by:

Shiv Tiwari, Intrexon, United States Barbara Molesini, University of Verona, Italy

\*Correspondence: Junhong Zhang zhangjunhng@mail.hzau.edu.cn

#### Specialty section:

This article was submitted to Plant Breeding, a section of the journal Frontiers in Plant Science

Received: 24 November 2017 Accepted: 14 June 2018 Published: 04 July 2018

#### Citation:

Wu L, Tian Z and Zhang J (2018) Functional Dissection of Auxin Response Factors in Regulating Tomato Leaf Shape Development. Front. Plant Sci. 9:957. doi: 10.3389/fpls.2018.00957

**Abbreviations:** ARF, auxin response factor; IAA, indole-3-acetic acid; LUC, firefly luciferase; qRT-PCR, quantitative realtime PCR; REN, renilla luciferase; VIGS, virus-induced gene silencing; YFP, yellow fluorescent protein.

blade dissection to form lobes or leaflets (Ben-Gera et al., 2012). Simple leaves comprise of a single continuous blade, whereas compound leaves are composed of multiple discontinuous blade units termed as leaflets (Koenig et al., 2009). Leaves are formed at the flanks of the shoot apical meristem (SAM). Following the initiation of new shoot morphogenesis, leaves establish the basic framework for shape and size. Subsequently, organogenesis of lateral appendages occurs through differentiation and expansion of leaf tissue (Shwartz et al., 2016). Furthermore, the wild type leaves consist of primary, secondary, and intercalary leaflets with lobed margins in tomato (Berger et al., 2009).

After the formation and differentiation of the leaf primordia in the SAM, the development of the leaf primordium occurs (Merelo et al., 2016). Previous studies have discovered that two mechanisms are involved in the development of the leaf primordia. The first occurs through mutual repression between KNOX proteins and ARP (AS1/RS2/PHAN) MYBdomain proteins (Waites et al., 1998; Tsiantis et al., 1999; Byrne et al., 2000; Ori et al., 2000). The second mechanism demonstrates that the formation of leaf delimitation is controlled by PINOID (PID) and auxin efflux carrier PIN-FORMED1 (PIN1) which mediates local auxin accumulation (Furutani et al., 2004). In leaf development, AS1 represses the expression of KNOX gene BP (BREVIPEDICELLUS), while the function of the local auxin maxima alongside AS1 remains partly dependent on BP regulation. The SHOOT MERISTEMLESS (STM) gene is needed for SAM formation and maintenance, which prevents AS1 gene expression in the meristem (Byrne et al., 2000). In addition, auxin activities and KNOX proteins might form a feedback loop to facilitate leaf meristem delimitation (Zgurski et al., 2005; Hay et al., 2006).

Studies have indicated that leaf development is coordinated by a cross-talk between different hormones (Shwartz et al., 2016) with auxin playing a crucial role. The precise distribution and location of auxin signaling regulates proper leaf development in a specific spatiotemporal developmental context (Bilsborough et al., 2011). The auxin maxima is concentrated on the leaflet initiation area in the development of tomato leaves, resulting in lamina growth patterns (Ben-Gera et al., 2016). Auxin acts as an inducer of organogenesis. There are postulated inhibitory fields around existing primordia which are thought to result from low concentrations of auxin (Reinhardt et al., 2003; de Reuille et al., 2006; Jonsson et al., 2006). Endogenous auxin levels and localization are altered in developing leaves leads to leaf simplification phenotypes (Shwartz et al., 2016). SlIAA9 is an auxin regulator belonging to the Aux/IAA transcription factors gene family. Not only do SlIAA9-RNAi plants display simple leaves and parthenocarpy instead of compound leaves and seeded fruit characteristic typically seen from wild type (AC), but also the silenced plants have auxin-related growth alterations (Wang et al., 2005). The adjustment of Aux/IAA and SlARF genes and the downregulation of MADS box genes mediate fruiting, and early fruit development is regulated by the regulatory and metabolic events both in the absence and presence of pollination/fertilization (Wang et al., 2009). Meanwhile, tomato e mutant is reported as a single-base deletion in the coding region of the SlIAA9 gene and exhibits single based lamina with primary leaves partially fused (Zhang et al., 2007), E mRNA is discovered throughout the leaf margin (Koenig et al., 2009). Thus, SlIAA9 plays a role in limiting lamina growth between developing leaflets by locally inhibiting auxin responses. GOB (GOBLET) encodes a NAC-domain transcription factor and its expression is intact in the simplified leaves of entire (e) mutants in tomato. Leaves of single gob or e mutants formed only primary leaflets, and downregulation of both GOB and E (SlIAA9) contributed to the complete abolishment of leaflet initiation. This indicates those auxin response and leaflet morphogenesis are modulated by GOB and E via partly redundant pathways (Blein et al., 2008; Ben-Gera et al., 2012). The tomato clau (clausa) mutant exhibits elaborate compound leaves. CLAU might negatively regulate the expression of GOB, and GOB expression is up-regulated in the compound leaf mutant lyr (lyrate). However, the enhancement of the clau phenotype by lyr indicates that clau and lyr affect GOB and leaf development in different pathways (Bar et al., 2015). Higher expression of LA (LANCEOLATE) during the early stages of leaf development result in a simpler leaf shape, likely regulated in part by gibberellic acid (GA) levels (Yanai et al., 2011). The expression of TKn1 in the leaf primordium is needed for compound structure formation (Hareven et al., 1996). miR164 negatively regulates GOB-like genes, and leafspecific overexpression of miR164 induces a loss of secondary leaflet initiation and smooth leaflet margins (Berger et al., 2009). The miR160 targets a group of ARFs which antagonize lamina growth and auxin response in conjunction with E plants. Leaflet separation is assured by different type of auxin signal antagonists (Ben-Gera et al., 2016). Auxin, E, GOB, LYR, and mir160 targeted ARFs collaborate to specify leaflet initiation and promote leaflet separation (Kimura et al., 2008). However, the underlying molecular mechanism of how functional redundancies among SlARF proteins regulate leaf shape development in tomato remains an open question.

Aux/IAA protein can repress ARF transcription factors via protein to protein interaction (Tiwari et al., 2001), and the degradation of Aux/IAA proteins can relieve ARF proteins for auxin-responsive gene transcription (Tan et al., 2007). SlIAA9 protein mediates leaf morphogenesis by participating in auxin signal transduction (Wang et al., 2005; Goetz et al., 2007). Furthermore, ARFs represent essential factors in the transduction of auxin signaling, and multiple ARFs were previously shown to interact with IAA9 (Korasick et al., 2014; Piya et al., 2014). In this work, we showed that tomato plants with a silenced SlIAA9 complex change from simple to complex leaf morphology, providing an insight for the significant role of functional redundancies among SlARF proteins in leaf morphogenesis. This indicates that SlARFs may mediate phenotypic plasticity in foliar organogenesis, and that further studies on SlARFs may reveal insights into the evolution of plant leaves.

#### MATERIALS AND METHODS

#### Plant Materials and Growth Conditions

Tomato plants (Solanum lycopersicum cv. Ailsa Craig) and e mutants in the Ailsa Craig background were grown under

standard greenhouse conditions (14 h day/10 h night cycle, 25/20◦C day/night temperature, 60–75% relative humidity). The e mutation plants prepared for VIGS assay were kept in a growth chamber (16 h day/8 h night cycle, 20–22◦C, 50% relative humidity).

#### qRT-PCR

Total RNA from all samples was isolated using the TRIzol reagent (Invitrogen, United States). The RNA was treated with DNase I at 37◦C for 30 min to remove residual genomic DNA. Using the HiScript II 1st cDNA Synthesis Kit (Vazyme, China) to synthesis the first-strand cDNA according to the manuscript's protocol. The cDNA concentrations were normalized according to actin expression levels for qRT-PCR analysis. qRT-PCRs were performed using the power SYBR Premix Ex Taq kit and the TaKaRa two-step method (TaKaRa, Japan). PCR products were quantified using the Roche Light Cycler 480 Real-Time PCR Detection System and the SYBR Green I Master Kit (Roche, Switzerland). The PCR program was as follows: 95◦C for 45 s; 40 cycles of 95◦C for 10 s, 58◦C for 25 s, and 72◦C for 20 s. For all qRT-PCR experiments, at least three biological replicates were performed, and each reaction was run in triplicate.

# Yeast Two-Hybrid Assay (Y2H)

For yeast two-hybrid assay, the full-length coding sequence of each SlARF (Supplementary Table S1) and SlIAA9 (Solyc04g076850) were cloned from various tissues of Ailsa Craig. Recombinant plasmid pGBKT7-SlIAA9 and pGADT7- SlARF were constructed by inserting SlIAA9 and SlARF into pGBKT7 and pGADT7 vectors separately. The two plasmids were co-transformed into the yeast strain AH109 by small-scale yeast transformation method. The transformants grew on the SD/-Trp/-Leu drop-out medium. After colony formation, transformants were transferred to SD/-Leu/-Trp/-His/-Ade drop-out medium with 40 µg ml−<sup>1</sup> X-gal.

# Bimolecular Fluorescence Complementation (BiFC) Assays

The SlIAA9 ORF without the stop codon was constructed using the pUC-SPYNE/pSPYNE-35S vector to produce SlIAA9-YFP<sup>N</sup> fusions, and the SlARF ORFs without the stop codon were cloned into the pUC-SPYCE/pSPYCE-35S vector to generate SlARFs-YFP<sup>C</sup> fusion proteins. Protoplasts were extracted from 1-week-old Arabidopsis Col-0 suspension cell culture and the corresponding constructs were co-transformed into them. The transfected protoplasts were assayed for fluorescence after 12– 18 h of expression. All primer sequences used in this analysis are listed in Supplementary Table S2.

#### In Vitro Pull-Down Assay

The recombinant plasmids were transformed into BL21(DE3)pLysS chemically competent cells. SlIAA9-GST was purified with Glutathione Agarose (Thermo Fisher Scientific, United States) according to the company manual instruction. MBP, MBP-SlARF6A, MBP-SlARF8A, MBP-SlARF8B, and MBP-SlARF24 were purified as fusion proteins immobilized with amylose resin (New England Biolabs, United States) following standard protocols. Five micrograms of GST-SlIAA9 protein were pre-incubated with 10 µL pre-washed amylose resin in 150 µL incubation buffer (1 mM NaCl, 20 mM MgCl2, 0.2% Triton X-100, and 0.1 M HEPES at pH7.2) for 1 h at 4◦C. The resin was collected by centrifugation and washed five times with washing buffer (20 mM Tris-HCl at pH 7.5, 300 mM NaCl, 0.1 mM EDTA, 0.5% Triton-X100). The pull-down proteins were detected by western blot with an α-GST antibody (Thermo Fisher Scientific, United States). All primer sequences used in this analysis are listed in Supplementary Table S2.

# VIGS Assays

To generate the VIGS constructs, 208 bp, 157 bp, and 238 bp fragments of the gene SlARF6A, SlARF8A, and SlARF24 was amplified by sequence-specific primers (Supplementary Table S2), respectively. Since SlARF8A and SlARF8B are highly homologous, the fragment from SlARF8A was used to silence both SlARF8A and SlARF8B. All of the pTRV1, pTRV2, pTRV2- PDS, and pTRV2-host target genes were transformed into the Agrobacterium tumefactions strain GV3101 by electroporation. Cultures containing the pTRV1 and pTRV2 vectors were mixed in a 1:1 ratio, either individually or simultaneously (pTRV2-PDS, pTRV2, pTRV2-SlARF6A, pTRV2-SlARF8A, pTRV2-SlARF24, pTRV2-SlARF6A/8A, pTRV2-SlARF6A/24, pTRV2-SlARF8A/24, pTRV2-SlARF6A/8A/24). These were used to infect cotyledon of tomato e mutant plants before the emergence of true leaves. The infected plants were transferred to a growth chamber at 16 h day/8 h night cycle, 20–22◦C and 50% RH. The phenotypes were analyzed 6–7 weeks after inoculation. The VIGS method was following the published protocol (Velasquez et al., 2009).

#### Quantification of the Free IAA Using UFLC-ESI-MS/MS

The sample leaves were frozen in liquid nitrogen. Three replicates were prepared for each leaf sample. The biomass for each replicate was 0.1 g. Subsequently, IAA extraction was performed by ESI-MS/MS following the published protocol (Liu et al., 2012).

# Yeast One-Hybrid Assay (Y1H)

For yeast one-hybrid assay, the −1479 bp fragment (upstream from the start codon) from the SlPIN1 (Solyc03g118740) promoter was amplified from Ailsa Craig genomic DNA and cloned into the pAbAi vector (Clontech). Recombinant plasmid pAbAi-SlPIN1 and pGADT7-SlARF were co-transformed into the yeast strain Y1HGold (Clontech) by small-scale yeast transformation method respectively. The transformants were plated on the SD/-Ura drop-out medium. Colonies were picked and diluted in sterile ddH2O to an OD600 of 0.5, and 3 µl of suspension was spotted on SD/-Ura/-Leu drop-out medium with or without AbA antibiotic at 30◦C. Both pGAD-p53+p53-AbAi (positive control) and pGADT7+P1-AbAi (negative control) were included.

### Transient Expression in Tobacco Leaves

The full-length SlARFs ORF were amplified and cloned into the effector vector, pGreen II 62-SK. A −1479 bp fragment (upstream from the start codon) from the SlPIN1 promoter was amplified and cloned into the reporter vector, pGreen II 0800-LUC. Both the effector and reporter vector were respectively co-transformed into the Agrobacterium tumefactions strain GV3101 cells with the pSoup vector, then infiltrated into N. benthamiana young leaves and incubated 72 h in the dark. LUC and REN were analyzed using the dual luciferase assay reagents (Promega) with an Infinite M200 (Tecan). All primers used in this analysis are listed in Supplementary Table S2.

#### RESULTS

#### SlIAA9 Interacts With Multiple SlARF Proteins

In order to dissect the mechanism of SlIAA9 regulating leaf shape development in tomato, a yeast two-hybrid (Y2H) screening was performed to identify the SlIAA9 interacting proteins from tomato cDNA Y2H library. After several screens, SlARF24 was screened out (Supplementary Table S3). To identify more SlARFs that may participate in this pathway, the full-length coding sequences of 15 SlARFs, including the candidate SlARF24 previously identified, were isolated and inserted into the yeast two hybrid vector pGADT7, including SlARF1, SlARF2B, SlARF3, SlARF4, SlARF5, SlARF6A, SlARF6B, SlARF8A, SlARF8B, SlARF9A, SlARF9B, SlARF10A, SlARF10B, and SlARF16A. The recombinant plasmid pGBKT7- SlIAA9 and pGADT7-SlARFs were co-transformed into the yeast, respectively. Specifically, we observed that yeast cells containing pGBKT7-SlIAA9 mated with pGADT7-SlARF6A, pGADT7-SlARF8A, pGADT7-SlARF8B, and pGADT7-SlARF24 respectively to grow under selection conditions (**Figure 1A**). This result suggests that SlIAA9 may interact with SlARF6A, SlARF8A, SlARF8B, and SlARF24 in yeast cells.

Subsequently, the interactions between SlIAA9 and SlARF6A, SlARF8A, SlARF8B, and SlARF24 were examined using BiFC in Arabidopsis mesophyll protoplasts. We transformed the Arabidopsis Col-0 protoplasts with SlIAA9-YFPN/SlARF6A-YFP<sup>C</sup> , SlIAA9-YFPN/SlARF8A-YFPC, SlIAA9-YFPN/SlARF8B-YFPC, SlIAA9-YFPN/SlARF24-YFPC, and SlIAA9-YFPN/YFPC. Strong YFP fluorescence signal was detected throughout the nucleus when SlIAA9-YFP<sup>N</sup> was co-expressed with SlARF6A-YFPC, SlARF8A-YFPC, SlARF8B-YFP<sup>C</sup> , and SlARF24-YFPC, whereas no fluorescence was detected in the control cells (**Figure 1B**). These results indicate that SlIAA9 might interact with multiple SlARF proteins in vivo.

In order to further confirm that the SlIAA9 protein could directly interact with SlARF6A, SlARF8A, SlARF8B, and SlARF24 proteins, a pull-down assay was performed. In this experiment, SlARF6A, SlARF8A, SlARF8B, or SlARF24 fused to maltose binding protein (MBP) immobilized on amylose-agarose beads were used as bait against GST-SlIAA9 fusion proteins. As shown in **Figure 1C**, GST-SlIAA9 could be pulled down by MBP-SlARF6A, MBP-SlARF8A, MBP-SlARF8B, as well as MBP-SlARF24, but not by MBP alone, demonstrating that SlARF6A, SlARF8A, SlARF8B, and SlARF24 proteins physically interact with SlIAA9 in vitro.

# The e Mutation Alters the Expression Patterns of Multiple SlARFs

A previous study showed that the basal gene expression of SlIAA9 was high in roots, leaves, flowers, and fruits (Wang et al., 2005). For comparative purposes, qRT-PCR was performed to investigate the expression levels of SlIAA9 in wild type root, stem, leaf, flower, and fruit tissue (Supplementary Figure S1). These results indicated that SlIAA9 expressed in all tissues, but exhibits lower expression at the mature green (MG) stage in fruit tissue.

We analyzed the expression levels of the SlARF genes family in wild type (compound leaves) and e mutant (simple leaves) leaves. The expression levels of SlARF1, SlARF2A, SlARF2B, SlARF5, and SlARF18 in leaves of the e mutant were induced 2.14-, 5.79-, 5.90-, 3.19-, and 2.73-fold more than that of those in wild type leaves. While the SlARF3, SlARF6A, SlARF6B, SlARF7B, SlARF10A, SlARF10B, SlARF16A, SlARF16B, SlARF19, and SlARF24 showed decreased expression corresponding to 0.05-, 0.48-, 0.45-, 0.47-, 0.15-, 0.15-, 0.44-, 0.01-, 0.05-, and 0.48-fold, respectively. Other SlARFs showed no significant difference (**Figure 2**). These results establish that multiple SlARFs are regulated in the absence of SlIAA9 in the auxin signaling model. This is particularly indicated for the SlARF16B gene, which might have a unique and significant function in this pathway. In addition, the expression levels of SlARF6A, SlARF8A, SlARF8B, and SlARF24 displayed varying degrees of attenuation in e mutant plants compared to wild type plants, supporting the hypothesis that the SlARF6A, SlARF8A, SlARF8B, and SlARF24 interaction with SlIAA9 has compromised biological function in e mutant plants.

#### Diminished Expression of Multiple SlARFs Can Rescue the Leaf Phenotype of Tomato e Mutation

The genetic mechanism by which SlIAA9 regulates leaf shape development in tomato was examined by silencing SlARFs in the e mutant background. Using TRV-mediated VIGS we carried out a functional characterization assay of the four candidate SlARF genes identified through the previous expressional analysis. SlARF8A and SlARF8B share 82% amino acid identity (Supplementary Figure S2), therefore both genes were silenced in one VIGS construct. There was no detectable change in leaf shape after individually silencing pTRV2-SlARF6A, pTRV2- SlARF8A, and pTRV2-SlARF24 (**Figure 3A**). Gene expression analysis showed that the mRNA levels of SlARF6A, SlARF8A, SlARF8B, and SlARF24 were reduced to approximately 43%, 48%, 49%, and 51% respectively compared to the empty vector control (**Figure 3B**). Functional redundancies among ARF proteins have been described in Cucumis sativus, A. thaliana, and S. lycopersicum (Okushima et al., 2005; Liu and Hu, 2013; Hao et al., 2015). In accordance with this redundancy, we inoculated pTRV2-SlARF6A, pTRV2-SlARF8A, and pTRV2- SlARF24 constructs into e mutant plants using double or triple

co-cultures of Agrobacterium. Interestingly, 80%, 73%, and 77% of the e mutant plants infiltrated with Agrobacterium triple co-cultures expressing pTRV2-SlARF6A/8A/24 were partially converted to compound leaves in three replicate experiments (**Table 1**). Meanwhile, the mRNA levels of SlARF6A, SlARF8A, SlARF8B, and SlARF24 in the e mutant plants inoculated with the Agrobacterium co-cultures of pTRV2-SlARF6A/8A/24 were reduced to approximately 60%, 57%, 69%, and 64% compared with the empty vector control, respectively (**Figure 3C**). Furthermore, the gene expression of other members of the SlARF family were not significantly changed after silencing the four candidate SlARFs (Supplementary Figure S3). Moreover, the fifth leaves of silencing of candidate SlARF genes in tomato e mutation were observed, the results presented that e mutant plants inoculated with pTRV2-SlARF6A/8A/24 triple combination cultures were partially converted to wild-type compound leaves, which generating more leaflets (**Figure 4A**). The total leaflets on the mature first five leaves of the pTRV2-SlARF6A/8A/24 inoculated plants were significantly increased compared to the e mutant plants (**Figure 4B**). The double co-cultures pTRV2- SlARF6A/8A, pTRV2-SlARF6A/24, and pTRV2-SlARF8A/24 could not restore the development of compound leaves (**Figure 4A** and **Table 1**) despite that the expression analysis revealed the target genes were down-regulated (**Figure 3B**). These results illustrated that simultaneously silencing of these four genes could restore the compound leaf shape in e mutant plants, and suggested functional redundancies among SlARF proteins in regulating tomato leaf shape development.

To analyze whether these phenotypic changes are regulated by auxin, the free IAA levels were quantified using UFLC-ESI-MS/MS. IAA levels of young leaves that had all four SlARFs silenced reached 20.14 ng/g, which were significantly higher than the e mutant plants (13.14 ng/g). Interestingly, the concentration of IAA in e mutant was observed to have a higher basal level than

the wild type control plants (9.74 ng/g) (**Figure 4C**). This result suggested that the development of compound leaves inoculated with pTRV2-SlARF6A/8A/24 were induced by auxin through promoting auxin response.

#### Downregulation of Multiple SlARFs Can Rescue the Expression of Leaf Shape Determining Genes

The leaf growth and development in tomato is likely driven by the leaf shape related genes (Hareven et al., 1996; David-Schwartz et al., 2009; Illing et al., 2009; Yanai et al., 2011; Pinon et al., 2013). Subsequently, the transcript levels of these leaf shape related genes in tomato were evaluated through qRT-PCR. The SlPIN1 and phan expression in pTRV2-SlARF6A/8AB/24 silenced plants was induced 1.8- and 2.5-fold compared to e mutant plants, respectively. In contrast, the expression level of LYR was reduced by 4.6-fold compared to e mutant plants (**Figure 5A**). Thus, after silencing the four SlARFs in e mutant, the expression of those genes were reverted back to wild type level. We hypothesize that SlARF proteins may regulate leaf shape development by regulating the expression of SlPIN1.

#### SlARF6A and SlARF24 Bind to the SlPIN1 Promoter

To investigate the relationship between the candidate SlARFs proteins and SlPIN1 promoter, a fragment of the SlPIN1 promoter 1,479 bp upstream from the start codon was used in an Y1H assay (**Figure 5B**). Y1H results demonstrated that this fragment could interact with SlARF6A and SlARF24 protein, confirming that SlARF6A and SlARF24 protein recognize the cis-element in the SlPIN1 promoter in yeast.

To determine whether SlARF6A and SlARF24 function as activator or repressor, we used a LUC (dual luciferase) assay to test how SlARF6A and SlARF24 interact with the SlPIN1 promoter. The same fragment of SlPIN1 used in the Y1H assay was introduced into the pGreen II 0800-LUC vector to generate the reporter construct (**Figure 5C**). The effector and reporter construct were transiently expressed in tobacco leaves and the relative LUC activity was determined. This result revealed that LUC activity was 2.29- and 1.69-fold higher in the presence of the SlARF6A and SlARF24 effector and reporter construct than in the negative control (**Figure 5C**), implying that both SlARF6A and SlARF24 may function as a transcriptional activator. This result revealed that the ciselement from the SlPIN1 promoter was bound by SlARF6A and SlARF24.

# DISCUSSION

### SlARF Proteins Regulate Tomato Leaf Shape in a Functionally Redundant Manner

Over the past 10 years, SlIAA9 has been shown to be involved in fruit development, leaf morphogenesis, and fruit parthenocarpy in A. thaliana and S. lycopersicum (Wang et al., 2005, 2009; Goetz et al., 2007). In this study, we aimed to identify interacting partners of SlIAA9. The initial Y2H assay identified several candidates, including Aux/IAA proteins, ubiquitin related proteins, gibberellin beta-hydroxylase protein, MADs box interactor-like protein, and SlARF24 (Supplementary Table S3). More SlARFs (SlARF6A, SlARF8A, SlARF8B, and

TABLE 1 | Statistical information describing the TRV-mediated virus-induced gene silencing (VIGS) of genes in e mutant plants.


a Indicates that the percentage of albino plants injected with pTRV2-PDS cultures reached to 100%.

SlARF24) were further found to interact with SlIAA9 in vivo and in vitro (**Figure 1**), indicating that other SlARFs may play redundant roles with SlARF24 in regulating leaf development. SlIAA9 and miR160-targeted ARFs SlARF10A, SlARF10B, or SlARF17, appear to partially act in a functionally redundant manner, but remain necessary for local inhibition of lamina growth between initiating leaflets (Ben-Gera et al., 2016). However, the elucidation of the molecular mechanism on how these functional redundancies among SlARF proteins regulate tomato leaf shape has been hampered by complexity of the protein in planta. Understanding the mechanisms involved in leaf shape development in tomato can provide new insights into understanding these same mechanisms in other species such as A. thaliana, Glycine max, and C. hirsuta.

SlARFs may regulate leaf morphology through binding to the promoter of SlPIN1 or other auxin-responsive genes in the absence of SlIAA9 (**Figure 5**). SlARF8 and SlIAA9 proteins, together with another unknown protein, may form a regulatory complex to control fruiting and growth, offering a possible explanation for the role of SlIAA9 in parthenocarpy (Goetz et al., 2007). SlARF6 and SlARF8 also play conserved roles in regulating development and growth of flower and vegetable organs in dicots (Liu et al., 2014). The Osarf24-1 mutant presents reduced sensitivity to aberrant auxin signaling and auxin-deficient phenotypes (Sakamoto, 2013). Here, we firstly found that SlARF6A, SlARF8A, SlARF8B, and SlARF24 could regulate tomato leaf development in a redundant manner.

#### SlARF Genes Play Distinct and Vital Roles in the Auxin Signaling Model

It has been illustrated that several SlARF genes might serve unique functions in tomato development. SlARF2A functions in the regulation of tomato fruit ripening as a recognized auxin signaling component (Breitel et al., 2016). Down-regulation of SlARF4 results in a dark-green fruit phenotype with increased chloroplasts densities (Jones et al., 2002). Furthermore, SlARF4 involves in the control of sugar metabolism during fruit development in tomato (Sagar et al., 2013). Both auxin and gibberellin responses are modulated by SlARF7 during fruit formation and development in tomato (De Jong et al., 2011). Compared with wild type fruits, the fruits of SlARF7-RNAi transgenic lines presented seedless, heart-shaped, and thick pericarp phenotypes in tomato (De Jong et al., 2009). Cell division is negatively regulated by SlARF9 during early fruit development in tomato (De Jong et al., 2015). Primexine

formation is modulated by ARF17, which is crucial for pollen wall patterning, partially through regulation of CalS5 gene expression in Arabidopsis (Yang et al., 2013). A 165-bp deletion in ARF18 gene simultaneously affects silique length and seed weight in polyploid rapeseed (Liu et al., 2015).

Moreover, there are functional redundancies among ARF proteins in Cucumis sativus, A. thaliana, and S. lycopersicum (Okushima et al., 2005; Liu and Hu, 2013; Hao et al., 2015). A constitutive expression pattern was exhibited in almost all of the ARF genes in cucumber (Liu and Hu, 2013). In A. thaliana, arf7 arf19 double mutant presented an obvious auxin-related phenotype that were not detectable in the single mutant, suggesting that there are functional redundancies between ARF7 and ARF19 proteins (Okushima et al., 2005). Simultaneous silencing of SlARF2A and SlARF2B genes leaded to severe ripening inhibition, clarifying a functional redundancies between SlARF2A and SlARF2B proteins (Hao et al., 2015). Our VIGS results provided evidence that functional redundancies among SlARF proteins resulted in the change from a simple leaf to a complex one in tomato e mutant plants (**Figure 3**). Due to the far evolutionary relationship among candidate SlARFs, clarifying there are functional compensation among the candidate SlARFs. It should be noted that this study has only examined the function of SlARF6A, SlARF8A, SlARF8B, and SlARF24, because only 15 tomato SlARF genes were isolated from the full-length cDNA sequences out of the 22 SlARF genes family. It has been reported that SlARF17 and SlIAA9 do not interact in yeast (Ben-Gera et al., 2016). There are still another six SlARF genes that have not yet been characterized by us. Our result showed that the expression levels of SlARF16B were 0.01-fold less in leaves of the e mutant compared to those in wild type. A previous study describes that the Pto-ARF16 was affected by Pto-miR160a associated with tree growth and wood property traits in Populus tomentosa (Tian et al., 2016). In A. thaliana, ARF10 and ARF16 were targeted by miR160 to control the formation of root cap cell, and miR160-uncoupled production of ARF16 reflected pleiotropic effects (Wang et al., 2005). Thus, we hypothesize that SlARF16B is regulated in the absence of SlIAA9 in the auxin signaling model, which needs to further evaluated with additional experiments. Consequently, the e mutation likely alters the expression patterns of other SlARF genes through this mechanism. The hetero-dimerization between Aux/IAA and ARF proteins likely able to play unique cellular functions (Piya et al., 2014). However, how the SlIAA9-SlARFs complex functions during tomato leaf development is still not yet completely resolved.

# A Proposed Model of SlIAA9 Complex in the Control of Tomato Leaf Forms

We chose several leaf shape related genes in A. thaliana and S. lycopersicum to evaluate whether these genes could

fpls-09-00957 July 2, 2018 Time: 15:32 # 8

indirectly by preventing SlARFs from functioning as transcriptional activators. In the absence of SlIAA9, SlARFs may regulate leaf growth and development through direct binding to the promoter of SlPIN1 or other auxin-responsive genes, which promoted auxin response.

regulate the development of leaf shape. As a result, the expression levels of SlPIN1, phan, and LYR in leaves having SlARF6A, SlARF8A, SlARF8B, and SlARF24 simultaneously silenced through VIGS were restored to levels similar to those in wild-type plants (**Figure 5A**). The IAA levels in leaves having those same four SlARFs silenced were significantly increased compared to the control (**Figure 4C**). Auxin acts as a positional cue during leaf organogenesis, and auxin efflux carrier PIN1 is one of the main contributors to auxin localization (Reinhardt et al., 2003; Cheng et al., 2007). PIN1 localizes on the periphery of apical meristems directing auxin to convergence points, where auxin maxima is formed, subsequently auxin becomes directed subepidermally at the leaf initiation site to regulate leaf development (Heisler et al., 2005; Martinez et al., 2016). Genetic analyses have also demonstrated that PIN1 is required for leaflet initiation in compound leaves (Barkoulas et al., 2008). Accordingly, the ciselement of SlPIN1 was used for deep analysis. The cis-elements from SlPIN1 promoter was recognized and bound specifically by SlARF6A and SlARF24 in yeast and plants (**Figure 5**). Here, we propose that SlARF6A and SlARF24 may regulate leaf growth and development through direct binding to the SlPIN1 promoter. However, the effects of enhanced SlPIN1 transcription still needs to be further evaluated. This enhanced transcription may result to increased expression of SlPIN1 protein, changed SlPIN1 protein modification, or shift the location of SlPIN1.

Our data provide an insight to suggest that SlARF proteins work with SlIAA9 in a functionally redundant manner to dictate leaf shape. We propose that SlIAA9 interacts with multiple SlARF proteins to promote the formation of a regulatory complex which can directly block leaflet initiation genes. We also assume that this complex may act indirectly by preventing SlARFs from

functioning as transcription activators. In the absence of SlIAA9, SlARFs may regulate leaf growth and development through direct binding to the promoter of SlPIN1 or unknown X genes induced by auxin (**Figure 6**). Future studies will be directed to dissect the relationship between SlIAA9 and SlARF6A, SlARF8A, SlARF8B, and SlARF24. We also intend to clone other SlARF genes to ascertain whether SlIAA9 has additional interactors with unique biological functions.

#### CONCLUSION

In conclusion, this study posits a proposed molecular mechanism of SlIAA9 complex in the control of tomato leaf forms (**Figure 6**). Our results firstly demonstrate that SlARF6A, SlARF8A, SlARF8B, and SlARF24 directly interact with SlIAA9, and the simple leaves of the e mutant are partially converted to wild-type compound leaves by silencing of all the four SlARFs. Meanwhile, SlARF6A and SlARF24 bind to the SlPIN1 promoter to regulate genes expression involved in leaflet initiation. Further studies are still needed to explore the underlying mechanism of SlARF proteins in modulating tomato leaf shape.

#### REFERENCES


### AUTHOR CONTRIBUTIONS

LW designed and performed the experiments, data analysis, and drafted this manuscript. JZ supervised the experiments and revised the manuscript. JZ and ZT designed all the experiments.

# FUNDING

This work was supported by grants from the National Natural Science Foundation of China (31471888 and 31772317) and the Fundamental Research Funds for the Central Universities (Program No. 2662015PY224).

#### ACKNOWLEDGMENTS

We thank Kevin L. Cox for critically revising of the manuscript.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2018.00957/ full#supplementary-material


promote leaf development in Arabidopsis. Development 133, 3955–3961. doi: 10.1242/dev.02545


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Wu, Tian and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

fpls-09-00957 July 2, 2018 Time: 15:32 # 11

# Chemical Composition and Crystal Morphology of Epicuticular Wax in Mature Fruits of 35 Pear (Pyrus spp.) Cultivars

Xiao Wu<sup>1</sup>† , Hao Yin<sup>1</sup> \* † , Zebin Shi<sup>2</sup> , Yangyang Chen<sup>1</sup> , Kaijie Qi<sup>1</sup> , Xin Qiao<sup>1</sup> , Guoming Wang<sup>1</sup> , Peng Cao<sup>1</sup> and Shaoling Zhang<sup>1</sup> \*

<sup>1</sup> Center of Pear Engineering Technology Research, State Key Laboratory of Crop Genetics and Germplasm Enhancement, Nanjing Agricultural University, Nanjing, China, <sup>2</sup> Institute of Horticulture, Zhejiang Academy of Agricultural Sciences, Hangzhou, China

#### Edited by:

Guusje Bonnema, Wageningen University & Research, Netherlands

#### Reviewed by:

Pinarosa Avato, Università degli Studi di Bari Aldo Moro, Italy Brian Farneti, Fondazione Edmund Mach, Italy

#### \*Correspondence:

Hao Yin yinhao85@gmail.com Shaoling Zhang slzhang@njau.edu.cn †These authors have contributed equally to this work.

#### Specialty section:

This article was submitted to Plant Physiology, a section of the journal Frontiers in Plant Science

Received: 04 January 2018 Accepted: 03 May 2018 Published: 23 May 2018

#### Citation:

Wu X, Yin H, Shi Z, Chen Y, Qi K, Qiao X, Wang G, Cao P and Zhang S (2018) Chemical Composition and Crystal Morphology of Epicuticular Wax in Mature Fruits of 35 Pear (Pyrus spp.) Cultivars. Front. Plant Sci. 9:679. doi: 10.3389/fpls.2018.00679 An evaluation of fruit wax components will provide us with valuable information for pear breeding and enhancing fruit quality. Here, we dissected the epicuticular wax concentration, composition and structure of mature fruits from 35 pear cultivars belonging to five different species and hybrid interspecies. A total of 146 epicuticular wax compounds were detected, and the wax composition and concentration varied dramatically among species, with the highest level of 1.53 mg/cm<sup>2</sup> in Pyrus communis and the lowest level of 0.62 mg/cm<sup>2</sup> in Pyrus pyrifolia. Field emission scanning electron microscopy (FESEM) analysis showed amorphous structures of the epicuticular wax crystals of different pear cultivars. Cluster analysis revealed that the Pyrus bretschneideri cultivars were grouped much closer to Pyrus pyrifolia and Pyrus ussuriensis, and the Pyrus sinkiangensis cultivars were clustered into a distant group. Based on the principal component analysis (PCA), the cultivars could be divided into three groups and five groups according to seven main classes of epicuticular wax compounds and 146 wax compounds, respectively.

Keywords: wax, GC-MS, crystal morphology, PCA, pears

#### INTRODUCTION

Cuticular wax is the product of a complex mixture of very-long-chain (VLC) aliphatic compounds and their oxygenated derivatives, including fatty acids, alkanes, alcohols, esters, aldehydes, ketones, and triterpenes (Kolattukudy, 1996). It has been widely reported that wax plays important roles in moderating gas exchange, limiting non-stomatal water loss, protecting plants against ultraviolet (UV) radiation and extreme temperature damage, self-cleaning behavior and providing mechanical support to maintain the integrity of plant organs (Wang et al., 2016). More and more studies have suggested that the fruit cuticular wax layer acts as the first protective barrier against fruit splitting and plays pivotal roles in the reduction of pathogenic and insect attacks, protection against mechanical damage (Chu et al., 2017). Studies also showed that cuticular wax play important roles in maintaining the postharvest quality and delaying fruit senescence. For example, the removal of the natural wax of blueberry fruits can accelerate the postharvest water loss and decay, reduce the sensory and nutritional qualities, and then shorten the fruit shelf-life (Chu et al., 2018). In addition, as one of the most important quality factors determining consumer demand, the apple's appearance during the postharvest storage were also determined by the

physical and chemical properties of wax composition (Glenn et al., 1990).

The chemical composition and morphology of the cuticular wax layer of several fruit species have been detected. For example, triterpenoids and β-diketones were the most prominent compounds in blueberry fruits, and a large amount of tubular wax was deposited on their surfaces (Chu et al., 2017). In the mature fruits of olives, triterpenes were the main compounds, and only epicuticular waxes were observed, arranged in various crystalloid structures, including granules, platelets, plates and rodlets (Lanza and Di Serio, 2015). Aldehydes, triterpenes and secondary alcohols were the most predominant components of cuticular wax from citrus and plum fruits (Ismail et al., 1977; Wang et al., 2014). The wax concentration, chemical composition and morphology were also reported to be positively correlated with disease resistance capabilities (Wu et al., 2017). Therefore, a clear understanding of the components and amounts of fruit wax is important for obtaining better fruit quality, improving disease resistance and developing postharvest treatment strategies. More importantly, the structure, composition and concentration of fruit wax vary among cultivars of the same species; thus, knowledge of the cuticular wax traits from different germplasms will also be helpful for the selection of breeding parents. To the best of our knowledge, the cuticular wax profiles of different cultivars from several fruit and vegetable species have been detected, such as various cultivars of apple (Belding et al., 1998), grape (Pensec et al., 2014), blueberry (Chu et al., 2017), persimmon (Tsubaki et al., 2012), tomato (Bauer et al., 2004), and pepper (Parsons et al., 2012).

As the third most important temperate fruit species, the pear belongs to the Rosaceae family. The Pyrus genus is genetically diverse with 1000s of cultivars, which can be divided into two major groups, Occidental (European) and Oriental (Asiatic) pears. The Oriental pears comprise Pyrus bretschneideri, Pyrus ussuriensis, Pyrus pyrifolia, and Pyrus sinkiangensis and are mainly grown in China, Korea, Japan, and other Asian countries, whereas the Occidental pear (Pyrus communis Linn.) is mainly produced in America, Italy, Spain, and Germany. Previous studies have reported that the amount of wax obtained from four pear varieties ('Pingguoli' and 'Xuehua', P. bretschneideri; 'Kuerle', P. sinkiangensis; 'Yuluxiang', hybrid cultivars) varied dramatically, with wax concentrations ranging from 0.32 mg/cm<sup>2</sup> ('Pingguoli') to 1.43 mg/cm<sup>2</sup> ('Kuerle') (Chen et al., 2014; Wu et al., 2017); however, in the mature fruit of pears belonging to the other three species, information about the epicuticular wax profile is very limited. More generally, the genetic relationships of several important fruit-related traits of pear have been reported, such as sugar, acid, and flesh color (Liu et al., 2016; Wu et al., 2014), but genetic relationship of pear cuticular wax was less reported. Thus, the investigation of the pear cuticular wax at the germplasm level is of significant value for both the genetic study and the higher wax pear breeding program.

In this study, 35 pear cultivars belonging to five different species and hybrid interspecies were selected for determination of the concentration, composition, structure and genetic relationship of epicuticular wax through GC–MS, FESEM and cluster analysis. The results in the present study will not only help further studies of the potential role of wax components on postharvest quality but also provide a foundation for plant breeding aimed at improving fruit epicuticular wax.

# MATERIALS AND METHODS

# Plant Material

The composition, concentration and structure of epicuticular wax compounds of the fruits of 32 pear cultivars (**Table 1**) were evaluated in 2016, and the data of other three pear cultivars ('Kuerle', 'Xuehua' and 'Yuluxiang') was originated from our previous research (Wu et al., 2017). The fruits at commercial maturity, without disease, infection or physical injuries, were selected for the following experiments. Thirty fruits of each cultivar were packed individually in a plastic net bag and delivered immediately to the laboratory at Nanjing Agricultural University. All the samples were examined immediately.

# Extraction of the Epicuticular Wax

The epicuticular wax was extracted according to the method of Li et al. (2014). After being washed with tap water and air dried, three groups of five pear fruits were fully immersed and agitated twice for 1 min in 600 mL chloroform under a fume hood. The solvent containing the waxes was transferred into pre-weighed vials and evaporated by a nitrogen-blowing instrument (JHD-001S, Shanghai Jiheng Industries Company, Ltd.) at 40◦C, and the wax weights were recorded.

# Determination of the Epicuticular Wax Concentration

The surface area of the pear fruits was calculated according to the method of Yin et al. (2011). The wax concentration (µg/cm<sup>2</sup> ) was calculated with the following formula:

$$\text{Wax concentration} = (\text{W}\_1 - \text{W}\_0) / \text{Sa}$$

where W<sup>1</sup> is the final weight of the vials (µg); W<sup>0</sup> is the initial weight of the vials (µg); and Sa is the total surface area of 5 pear fruits (cm−<sup>2</sup> ).

# Chemical Analysis by GC-MS

The components present within each wax extraction prepared were analyzed using the method of Li et al. (2014). The wax extraction (1 mg) dissolved with 1.2 mL chloroform was carried out on a Bruker 450-GC, coupled with a Bruker 320-MS and a BR-5MS capillary column (30 m FS, 0.25 µm ID, 0.25 µm df). Helium was used as a carrier gas at a flow rate of 1.2 mL min−<sup>1</sup> . The following parameters were employed: inlet temperature, 280◦C; MS transfer line temperature, 280◦C; ion source temperature, 250◦C; quadrupole temperature, 150◦C; electron impact (EI), 70 eV; and m z−<sup>1</sup> range, 50–650.

GC was carried out at the following temperature settings. First, the temperature was set to 50◦C for 2 min. Next, it was increased to 200◦C at a rate of 40◦C min−<sup>1</sup> and held at this temperature for 2 min. Finally, it was increased to 320◦C at a rate of 3◦C min−<sup>1</sup> and held at this temperature for 30 min.

#### TABLE 1 | Cultivars of pears used in the study.

fpls-09-00679 May 18, 2018 Time: 16:55 # 3


(Continued)

TABLE 1 | Continued


The data of 'Kuerle', 'Xuehua' and 'Yuluxiang' was originated from our previous research (Wu et al., 2017).

#### Electron Microscope Observations

The epicuticular wax was observed by FESEM according to the method of Wu et al. (2017). Pericarp pieces (3 × 3 × 1 mm) from the equatorial sections of three fruits for each cultivar were excised using a blade, fixed in 2.5% glutaraldehyde. The dehydrated samples were attached to a sample stage with conductive tape and coated with gold particles using a Hitachi E-1045. The coated samples were examined using a Hitachi 4800 field emission scanning electron microscope.

#### Statistical Analysis

The chemical composition of epicuticular wax was analyzed using the NIST 2013 Library. Box-plot, principal component analysis (PCA) and heatmap analysis of the dataset for wax concentration and composition were used to detect clustering and to establish relationships; these analyses were carried out using the R programming language. For the heatmap and clustering analysis, the scaled data file with Log transformed was calculated, and analyzed by the 'complete' algorithm in the heatmap.2 package (Murtagh and Legendre, 2014). We performed the principal component analysis (PCA) by the algorithm of 'Multivariate Analysis' in the 'prcomp' module of the R package, and 3D plot of PCA was drawn by the scatterplot3d packages (Mardia et al., 1979). The SPSS20 statistical software package (IBM Software Group) was used for all statistical analysis. Data were compared by Student's t-test, and significant differences are marked with different letters (P < 0.05). Data are shown as the mean ± SD (n = 3).

#### RESULTS AND DISCUSSION

#### Wax Amounts and VLC Aliphatic Compounds

The total wax of mature fruits from 35 pear cultivars was collected via chloroform extraction (**Table 1**). The amount of wax obtained from the 35 pear varieties varied dramatically (**Figure 1** and Supplementary Tables S1, S2), and the wax concentrations ranged from 0.46 ± 0.03 mg/cm<sup>2</sup> ('Huagai', P. ussuriensis) to 2.44 ± 0.41 mg/cm<sup>2</sup> ('Docteur Jules Guyot', P. communis) (**Figure 1A**). On the species level, the epicuticular wax concentration was the highest in P. communis (1.53 ± 0.77 mg/cm<sup>2</sup> ) and the lowest in P. pyrifolia (0.62 ± 0.12 mg/cm<sup>2</sup> ) (Supplementary Table S2). Similar results have been identified among other species, such as pepper (2.15–7.52 µg/cm<sup>2</sup> ) (Parsons et al., 2012), persimmon (337–770 µg/cm<sup>2</sup> ) (Tsubaki et al., 2012), tomato (27–79 µg/cm<sup>2</sup> ) (Bauer et al., 2004), blueberry (48–172 µg/cm<sup>2</sup> ) (Chu et al., 2017), grape (3.5–5.5 mg/g) (Pensec et al., 2014) and apple (366–1038 µg/cm<sup>2</sup> ) cultivars (Belding et al., 1998). Compared with these species, the average wax concentration of pears (1.05 mg/cm<sup>2</sup> ) was the second highest, lower than only that of grapes. Various methods of wax extraction and surface area calculation were employed in these species; thus, the considerable variation of epicuticular wax obtained from different species was reasonable. In total, 146 epicuticular wax compounds were detected in the fruits of 35 cultivated pear species (Supplementary Table S1 and Supplementary Figure S1), mainly primary alcohols (26 compounds), alkanes (25 compounds), fatty acids (25 compounds), terpenoids (23 compounds), esters (16 compounds) and aldehydes (8 compounds), and some other small amounts of compounds (23 compounds, Supplementary Table S1); the content of these wax compounds accounted for 24.47, 40.72, 7.09, 11.80, 3.59, 3.45, and 8.88% of the total content, respectively (**Figure 1A**). Although the primary alcohols comprised the highest number of wax compounds, the alkanes were the dominant compounds in the epicuticular wax of pear fruits in terms of amounts. It has also been reported that alkanes are the major wax component in many fruits, including tomato (Leide et al., 2007), mandarin (Sala, 2000) and grapefruit (McDonald et al., 1993). However, limited studies reported that primary alcohols are essential

cultivars. The horizontal lines in the interior of the box are the median values. The height of the box is equal to the interquartile distance, indicating the distribution for 50% of the data. Approximately 99% of cultivars fall inside the whiskers, and the cultivars outside these whiskers are considered outliers or extreme values, indicated by horizontal lines and circles.

components in fruits. Similar to the total wax content, the amount of each component also varied among the five different species and hybrid interspecies (Supplementary Figure S2). For example, the concentrations of alkanes, primary alcohols, fatty acids, esters and other epicuticular wax compounds were much higher in P. communis than in the other four species and hybrid interspecies (Supplementary Figures S2B,D,F,L,N), whereas the hybrid cultivars and P. ussuriensis had the highest concentrations of terpenoids and aldehydes, respectively (Supplementary Figures S2H,J). Similar results were also reported in blueberry fruits from nine blueberry cultivars belonging to three different species (Chu et al., 2017).

Very-long-chain aliphatic compounds were the major epicuticular wax components found in the fruit surfaces of mature pears. Abundant VLC aliphatic compounds were detected in the 35 pear cultivars, and their concentrations varied greatly, ranging from 336.85 µg/cm<sup>2</sup> ('Huagai', P. ussuriensis) to 769.87 µg/cm<sup>2</sup> ('Kuikeamute', P. sinkiangensis) (Supplementary Table S2). Our results also showed that the VLC aliphatic compounds of pear fruits are mainly composed of alkanes (C16–C44), primary alcohols (C17–C41), fatty acids (C16–C26), and aldehydes (C16–C18) (**Figure 2**).

Of all the detected 35 pear cultivars, 20 (10 hybrid combinations) were genetically related (Supplementary

Figure S3). This was designed to further understand the regular relationship of the epicuticular wax between the hybrids and their parents. Our results showed that five hybrid cultivars ('Hongxiangsu', 'Yuluxiang', 'Xinli No. 7', 'Chuxialü' and 'Cuiyu') with higher concentration were consistent with that of their parents (Supplementary Figures S3A–D), other three hybrid cultivars ('Huangguan', 'Xueqing' and 'Cuiguan') were also similar with their parents obtained lower concentration of epicuticular wax (Supplementary Figures S3G,H), and only two hybrid cultivars ('Jinfeng' and 'Zhongli No. 2') with higher concentration were contrary with that of their parents (Supplementary Figures S3E,F). These indicate that the epicuticular wax concentration of hybrid pear cultivars was positively correlated with that of the parents. Take an example, 'Kuerle' (Pyrus sinkiangensis Yü.) with relatively higher wax concentration is the most famous pear cultivar in China, not only for its pleasant appearance, aroma and taste, but also the longer postharvest storage period. Here, three hybrid cultivars ('Hongxiangsu', 'Yuluxinag' and 'Xinli No. 7') have all inherited the higher wax concentration from 'Kuerle', making it an excellent parent for pear breeding. In addition, 'Docteur Jules Guyot' (P. communis), 'Tianjianba' (P. ussuriensis), 'Qiubai' (P. bretschneideri) and 'Shinseiki' (P. pyrifolia) should be the corresponding optimal parent cultivars in each species for the pear breeding with higher wax concentrations.

#### Alkanes

Alkanes are important VLC aliphatic compounds in the epicuticular wax of pear fruits, with concentrations ranging from 114.19 µg/cm<sup>2</sup> ('Xuehua', P. bretschneideri) to 1017.49 µg/cm<sup>2</sup> ('Clapp Favorite', P. communis), respectively (**Figure 1A** and Supplementary Figure S2A). On the species level, the alkane concentrations were the highest in P. communis (601.14 µg/cm<sup>2</sup> ), followed by hybrid cultivars (442.75 µg/cm<sup>2</sup> ), P. ussuriensis (408.77 µg/cm<sup>2</sup> ), P. sinkiangensis (578.50 µg/cm<sup>2</sup> ) and P. pyrifolia (408.77 µg/cm<sup>2</sup> ), but they were the lowest in P. bretschneideri (281.72 µg/cm<sup>2</sup> ) (Supplementary Figure S2B). Many studies have shown that alkanes are predominant in epicuticular wax on many fruits, such as mandarin (148– 234 mg/g) (Sala, 2000), pepper (1.4–47.4 µg/dm<sup>2</sup> ) (Parsons et al., 2012) and tomato (55%, 13.75 µg/cm<sup>2</sup> ) (Bauer et al., 2004); leaves, including those of Schefflera elegantissima (∼80%, 400 µg/cm<sup>2</sup> ), Garcinia spicata (∼72%, 864 µg/cm<sup>2</sup> ) (Jetter and Riederer, 2016), castor bean (64.1%–70.7%, 88.1–139.4 µg/cm<sup>2</sup> ) (de Araújo Silva et al., 2017) and aloe (34.4%, 8.8 µg/cm<sup>2</sup> ) (Racovita et al., 2015); and both fruits and leaves of cucumber

(46.1%, 59.6 µg/cm<sup>2</sup> ; 61.9%, 110.3 µg/cm<sup>2</sup> ) (Wang W. et al., 2015). Compared with these species, the average wax concentration of pears (422.73 µg/cm<sup>2</sup> ) is the second highest, only lower than that of Garcinia spicata. VLC alkanes ranging from hexadecane (C16) to tetratetracontane (C44) were detected in the fruit wax of the 35 pear species, and C<sup>28</sup> to C<sup>31</sup> were the most abundant. For example, the average amount of hentriacontane (C31, 238.94 µg/cm<sup>2</sup> ) was 3.70 times higher than the amount of the least abundant alkane (**Figure 2A**). Similarly, hentriacontane (C31) was also detected as the predominant alkane in cuticle wax of 'Newhall' navel orange (Wang et al., 2014), tomato (Bauer et al., 2004), pepper, eggplant (Bauer et al., 2005) and aloe (Racovita et al., 2015). It has also been reported that alkanes are associated with limiting water loss (Zhang et al., 2007). Furthermore, low temperature induced the alkane-forming pathway resulting in the accumulation of very long chain alkanes (including C22, C27, C29, C31 alkane) and their derivatives on the surface of apple fruit peel, which may be positive response to the environment signals of temperature and extension of the storage period (Hao et al., 2017). Thus, the higher proportion of alkanes in epicuticular wax of pear fruits may play a significant role in limiting water loss and prolonging the shelf-life to keep the fruits juicy.

#### Primary Alcohols

Twenty-six primary alcohols were detected as the second dominant component in the epicuticular wax of pear fruits, accounting for 24.47% of the total wax (**Figure 1A** and Supplementary Table S1). The concentrations of primary alcohols in the epicuticular wax of pear fruits ranged from 50.54 µg/cm<sup>2</sup> ('Cuiguan', hybrid cultivars) to 1242.11 µg/cm<sup>2</sup> ('Docteur Jules Guyot', P. communis), and the average concentration was 271.04 µg/cm<sup>2</sup> (**Figure 2C**). On the species level, the primary alcohols were the highest in P. communis (505.53 µg/cm<sup>2</sup> ), followed by P. sinkiangensis (337.77 µg/cm<sup>2</sup> ), P. ussuriensis (259.60 µg/cm<sup>2</sup> ), hybrid cultivars (228.15 µg/cm<sup>2</sup> ) and P. bretschneideri (174.87 µg/cm<sup>2</sup> ), but the primary alcohol levels were the lowest in P. pyrifolia (65.75 µg/cm<sup>2</sup> ) (Supplementary Figure S2D). Higher concentrations of primary alcohols have been detected only in the epicuticular wax of leaves, such as wheat (42.5–72.9%) (Wang Y. et al., 2015), salix (32– 51%) (Teece et al., 2008) and citrus plants (23.0–38.4%) (Baker et al., 1975). By contrast, the epicuticular wax of fruits with lower concentrations of primary alcohols have been reported in several species, such as 3.2% in blueberries (Chu et al., 2017), 0–9.1% in apples (Belding et al., 1998) and 6.0–15.0% in citrus fruits (Baker et al., 1975), respectively. A wide distribution of VLC primary alcohol homologues ranging from C<sup>17</sup> to C<sup>41</sup> was detected in the wax of pear fruits. As shown in **Figure 2B** and Supplementary Table S1, the C<sup>30</sup> alcohols, including triacontanol and triacontane-1,30-diol (average 115.66 µg/cm<sup>2</sup> ), were the most abundant primary alcohols; a similar result has been reported in the epicuticular wax of blueberry fruits (Chu et al., 2017). In addition, tetracosanol (C24), dotriacontanol (C32), and octacosanol (C28) were the dominant primary alcohols in the fruit wax of lemons (Baker et al., 1975) and tomatoes (Bauer et al., 2004) and the leaf wax of wheat (Wang Y. et al., 2015), respectively, whereas hexacosanol (C26) has been detected as the main primary alcohol in the leaf wax of both orange (Baker et al., 1975) and salix plants (Teece et al., 2008).

#### Fatty Acids

Sixteen VLC fatty acids ranging from C<sup>16</sup> to C<sup>26</sup> were detected in the fruit wax of 35 pear species. The highest concentration of fatty acids was observed in 'Docteur Jules Guyot' (P. communis, 225.27 µg/cm<sup>2</sup> ), whereas the lowest was detected in 'Kucheamute' (P. sinkiangensis, 12.28 µg/cm<sup>2</sup> ), and the average concentration was 74.44 µg/cm<sup>2</sup> (Supplementary Figure S2E). On the species level, the fatty acids were the highest in P. communis (132.95 µg/cm<sup>2</sup> ), followed by P. pyrifolia (92.00 µg/cm<sup>2</sup> ), hybrid cultivars (74.31 µg/cm<sup>2</sup> ), P. ussuriensis (67.48 µg/cm<sup>2</sup> ), and P. bretschneideri (57.99 µg/cm<sup>2</sup> ), but they were the lowest in P. sinkiangensis (32.66 µg/cm<sup>2</sup> ) (Supplementary Figure S2F). The fatty acid profile was predominantly composed of C<sup>16</sup> and C<sup>18</sup> carbons, and a similar result has been reported in the epicuticular wax of peach fruits (**Figure 2C**) (Belge et al., 2014). Fatty acids with moderate concentrations, being a common component, were widespread in the epicuticular wax. Compared with pears (C16–C26), the carbon chain length of fatty acids in other fruit species (C16–C34) is longer, and the prominent components of fatty acids are more than 20 carbon atoms. For example, the fatty acid profiles in fruits of blueberry (C16–C30) (Chu et al., 2017) and lemon (C16–C34) were predominantly composed of C<sup>30</sup> chains, and the dominant fatty acids in fruits of orange (C16–C34), clementine (C16–C32) and mandarin (C16–C28) consisted of 28 carbon atoms, while in the leaves of these three citrus species (all C16–C34), chains with 32 carbon atoms were the major fatty acid component (Baker et al., 1975).

#### Aldehydes

Aldehydes were found in all cultivars but at low concentrations in the epicuticular wax of pear fruits, accounting for 3.45% of the total wax (**Figure 1A**). The highest concentration of aldehydes was observed in 'Hongnanguo' (P. ussuriensis, 128.69 µg/cm<sup>2</sup> ), whereas it was not detected in 'Cuiguan' (hybrid cultivars), and the average concentration was 35.32 µg/cm<sup>2</sup> (Supplementary Figure S2I). On the species level, the aldehydes were the highest in P. ussuriensis (87.17 µg/cm<sup>2</sup> ), followed by P. sinkiangensis (35.04 µg/cm<sup>2</sup> ), P. communis (32.00 µg/cm<sup>2</sup> ), hybrid cultivars (25.35 µg/cm<sup>2</sup> ) and P. bretschneideri (24.40 µg/cm<sup>2</sup> ), and they were the lowest in P. pyrifolia (6.84 µg/cm<sup>2</sup> ) (Supplementary Figure S2J). A series of aldehydes composed of 16, 17, and 18 carbon atoms were detected in pear wax, and octadecanal, 9-octadecenal,(Z)- and 13-octadecenal,(Z)- (C18, average 25.87 µg/cm<sup>2</sup> ) was the most abundant aldehyde (**Figure 2D** and Supplementary Table S1). However, similar to the fatty acids, the carbon chain length of aldehydes was longer in other species than in pears, and the concentration and main type of aldehydes in different species showed considerable variation. The aldehydes were the main components in the total wax of citrus fruits including 'Satsuma' (C16–C32, 45% of the total wax) and 'Newhall' navel orange (C16–C33, 46% of the total wax) fruits (Wang et al., 2014). Aldehydes comprising 24 and 32 carbon atoms were found in olive wax and accounted for the smallest

proportion (3–5%) of total wax (Bianchi et al., 1992). As direct precursors of alkanes, aldehydes determine the alkane formation. Similar to alkanes, aldehydes have been identified as playing an important role in avoiding water loss (Wang et al., 2014).

#### Esters

Similar to aldehydes, the proportion of esters accounted for 3.59% of the total wax. The highest concentration of esters was observed in 'Abbe Fetel' (P. communis, 101.63 µg/cm<sup>2</sup> ), whereas the lowest was detected in 'Huagai' (P. ussuriensis, 9.34 µg/cm<sup>2</sup> ), and the average concentration was 36.77 µg/cm<sup>2</sup> (Supplementary Figure S2K). On the species level, the ester concentrations were the highest in P. communis (75.01 µg/cm<sup>2</sup> ) and the lowest in P. sinkiangensis (21.78 µg/cm<sup>2</sup> ) (Supplementary Figure S2L). It is noteworthy that 1-monopalmitin and 1-monostearin were the only two glycerides compounds simultaneously detected in all 35 pear cultivars in this study (Supplementary Table S1). Glycerides were interesting ester compounds which have been rarely detected in waxes, till now, they have only been detected in the cuticular waxes of potato tuber and olives fruits (Vichi et al., 2016; Guo and Jetter, 2017). In addition, the epicuticular wax of fruits, esters were also detected as being the smallest proportion of the total wax in many species, including plums (8.2%) (Ismail et al., 1977), grapes (<5.1%) (Casado and Heredia, 1999) and apple fruits (1.39–4.85%) (Veraverbeke et al., 2001), while esters were not detected in blueberry (Chu et al., 2017) or peach fruits (Belge et al., 2014).

# Terpenoids

In addition to aliphatic compounds, a total of 23 terpenoid compounds, accounting for 11.80% of the total wax, were detected in the epicuticular waxes of the 35 pear cultivars. The highest concentration of terpenoids was observed in 'Jinfeng' (hybrid cultivars, 514.02 µg/cm<sup>2</sup> ), whereas the lowest was detected in 'Red Clapp Favorite' (P. communis, 14.96 µg/cm<sup>2</sup> ), and the average concentration was 36.77 µg/cm<sup>2</sup> (Supplementary Figure S2G). On the species level, the terpenoids were highest in hybrid cultivars (165.91 µg/cm<sup>2</sup> ), followed by P. ussuriensis (118.62 µg/cm<sup>2</sup> ), P. bretschneideri (115.77 µg/cm<sup>2</sup> ), P. sinkiangensis (76.25 µg/cm<sup>2</sup> ) and P. communis (62.58 µg/cm<sup>2</sup> ) and the lowest in P. pyrifolia (26.89 µg/cm<sup>2</sup> ) (Supplementary Figure S2H).

Terpenoids was the predominant component of epicuticular wax in many species, such as the blueberry (64.2%) (Chu et al., 2017), grape (34–49%) (Pensec et al., 2014), peach (44.05–51.92%) (Belge et al., 2014) and sweet cherry (75.6%) (Peschel et al., 2007). Similar to the fruit wax of the tomato (Bauer et al., 2004), grapefruit (Nordby and McDonald, 1994) and citrus (Wang et al., 2014) and the leaf wax of eucalyptus (Viana et al., 2010), we detected lanosterol, farnesene, amyrin, lupeol, squalene, sitosterol and urs-12-en-28-al in the epicuticular wax of pear fruits (Supplementary Table S1). It has been reported that the accumulation of farnesene in the cuticular waxes of apple was the significant cause of greasiness (Christeller and Roughan, 2016). Several certain pear cultivars, such as 'Yuluxiang', 'Kuerle' and 'Hongxiangsu', developing a greasy surface, have accumulated more fluid wax constituents, which were mainly of farnesene. Although ursolic acids and oleanolic acids are the prominent terpenoid components in the cuticular waxes of apple (32.03–69.8%) (Belding et al., 1998), blueberry (16.4–26.0%) (Chu et al., 2017), peach (17.21–29.19%) (Belge et al., 2014), plum (1.0–5.4%) (Ismail et al., 1977), and sweet cherry (7.5–60.0%) (Peschel et al., 2007), they were not detected in epicuticular wax of pear fruits. In addition to the common terpenoid compounds, the pear fruits also contained specific compounds, such as 20α-dihydro pregnenolone (0–7.77 µg/cm<sup>2</sup> ) (Supplementary Table S1).

Terpenoids are of significant importance for maintaining the mechanical strength of the cuticular membrane in persimmon fruit, as well as for defending against biotic and abiotic stress, delaying fruit senescence and being involved in important biological activities. Here, we detected four types of tocopherols in 30 cultivated pear species (except 'Hongnanguo', 'Huangguan', 'Nanguo', 'Docteur Jules Guyot' and 'Kousui'), including deltatocopherol, beta-tocopherol, gamma-tocopherol and alphatocopherol (Supplementary Table S1). As promising sources of biological activities, all four types of tocopherols were simultaneously detected in the mature fruits of 'Xizilü.' Tocopherol plays important roles in inhibiting the growth of pathogens and protecting fruits against biotic stresses, and squalene may extend the storage time by removing free radicals and enhancing oxygen loading in citrus fruits (Wang et al., 2014; Wu et al., 2017). Tocopherol and squalene concentrations were 0–81.07 µg/cm<sup>2</sup> and 0–7.34 µg/cm<sup>2</sup> in the epicuticular waxes of the 35 pear cultivars (Supplementary Table S1). Prior studies have noted that tocopherol and squalene can lead to enhanced resistance to Alternaria rot in pear fruits (Wu et al., 2017).

# Crystal Morphology

To elucidate the crystal morphology of the 35 pear cultivars, wax structures were detected by FESEM. Interestingly, the wax morphology showed various amorphous structures in different cultivars, including crystalline plates, irregular ovate crystals, platelets and rodlets with wax crystals (**Figure 3** and Supplementary Figures S4, S5). For example, the mature fruits of 'Qiubai', 'Kuerle', 'Clapp Favorite' and 'Jinfeng' were covered by wax with irregular ovate crystals (**Figures 3D1–D4** and Supplementary Figures S4C1–C4,I1–I4, S5I1–I4), whereas the fruit wax of 'Chili', 'Eli', 'Hongnanguo', 'Balixiang', 'Lüamute', 'Abbe Fetel', 'Xueqing' and 'Yuluxiang' were composed of numerous transversely polygonal rodlets and platelet crystals (**Figures 3A1–A4** and Supplementary Figures S4A1– A4,G1–G4,H1–H4,K1–K4, S5C1–C4,L1–L4,O1–O4), and glossy transversely polygonal rodlet wax covered the fruit surfaces of the cultivars including 'Zhongli No. 2' and 'Pingguoli' (**Figures 3G1–G4** and Supplementary Figures S5D1–D4). The mature fruits of 'Kucheamute' and 'Hongxiangsu' were covered by wax crystals with vertically polygonal rodlets with numerous platelets (**Figures 3C1–C4,G1–G4**), and 'Nanguoli,' 'Huagai' and 'Xizilü' were covered by wax crystals with glossy vertically polygonal rodlets (**Figures 3B1–B4** and Supplementary Figures S4E1–E4,S5K1–K4). Remarkably, on the fruit surfaces of 'Zhongli No. 2', 'Eli', 'Abbe Fetel', 'Bartlett Max Red' and 'Yuluxiang', these structures separated from the wax

FIGURE 3 | The morphology of mature fruits and the magnification series of the FESEM images of epicuticular wax in seven pear cultivars. (A1–G1) Are mature fruits of 'Chili', 'Nanguoli', 'Kucheamute', 'Clapp Favorite', 'Shinseiki', 'Hongxiangsu' and 'Zhongli No. 2', respectively. Scale bars represent 100 µm (magnification 500×) in (A2–G2), 10.0 µm (magnification 5000×) in (A3–G3) and 1.00 µm (magnification 50,000×) in (A4–G4). The white arrow denotes the most prominent wax crack.

FIGURE 4 | Heatmap and clustering of 146 wax components in 35 pear cultivars. Colors indicate chemical composition levels. The numbers of each row represent the wax component numbers corresponding to the accession numbers in Supplementary Table S1. Each column represents a pear cultivar.

layer, producing cracks and discontinuity in the outer layer of wax (**Figures 3G1–G4** and Supplementary Figures S4A1–A4, S5C1–C4,D1–D4,O1–O4).

It is generally believed that the various epicuticular wax crystal types, such as rodlets, tubules and platelets, deriving from self-assembly processes, were mainly based on their different chemical compositions (Chu et al., 2017). It has been reported that the tubule-shaped wax crystal is mainly determined by higher proportions of secondary alcohol tubes and β-diketone wax components in blueberry fruits (Chu et al., 2017), and the platelet wax crystals are caused by a higher proportion of aldehydes and alkanes in the epicuticular waxes of 'Newhall' navel orange fruits (Wang et al., 2014). However, it seems that similar wax structures can be formed by different chemical compositions in different species. For example, platelets are the most common wax structures, and wax platelets in Sedum rupestre are characterized by a high amount of triterpenoids, whereas platelets found in the Poaceae (e.g., Triticum, Zea) are generally dominated by primary alcohols (Stevens et al., 1995; Koch et al., 2006). However, many diverse structural and intermediate forms, environmental effects, and erosion of the waxes lead to undetermined shapes (Koch and Ensikat, 2008). Alkanes and primary alcohols were the dominant compounds of the epicuticular wax of pear fruits. The epicuticular wax had an amorphous structure in 35 pear cultivars, which makes classification of the pear fruits unclear. This denotes that wax structures in pear species are not only based on their dominant chemical composition but also influenced by other factors. It has been reported that wax structures on plant surfaces could be influenced by the molecular order (crystallinity) and orientation of the polar groups (polarity) of cutin, and the resulting template effects of the substrates could be used to influence the orientation and preferred direction of crystal growth (Koch et al., 2006). Given the present lack of knowledge about the molecular order and polarity of the cutin network, the factors contributing to the growth of wax structures on the cuticles of pear fruits still need further study.

# The Cluster Analysis of the 35 Pear Cultivars

Prior studies have noted that various chemical compositions of epicuticular wax, such as alkanes (Li et al., 2013), terpenes (Medina et al., 2006), flavonoids (Essokne et al., 2012), and fatty acids (Velasco and Goffman, 1999), could be a useful taxonomic marker for classifying the family, genus or species in higher plants and can reflect both ecological and genetic relationships. To elucidate epicuticular wax relationships among 35 pear cultivars, the concentrations of 146 epicuticular wax compounds found in the mature fruits of 35 pear cultivars were clustered through heatmap analysis (**Figure 4**). The dendrogram generated from the heatmap of wax composition could be divided into two major groups (**Figure 4**). Group one, including P. bretschneideri, P. ussuriensis, P. communis, P. pyrifolia and the hybrid cultivars, could be further divided into two subgroups (I–II). Subgroup (I) contains all five P. ussuriensis, four P. bretschneideri, seven hybrid cultivars and two P. communis cultivars, whereas the remaining three P. communis, one P. bretschneideri, two P. pyrifolia and four hybrid cultivars were clustered into Subgroup (II). The close genetic relationships between P. bretschneideri, P. ussuriensis, and P. pyrifolia based on our data are consistent with previous results using isozymes, SSR markers and RAPD markers (Bao et al., 2007). Group two was composed of P. sinkiangensis and two hybrid cultivars. All five cultivars of P. sinkiangensis were grouped into the same clade (**Figure 4**), indicating that the components and contents of their epicuticular wax are

similar, and their relationships are phylogenetically close. In addition, 'Yuluxiang' and 'Kuerle' were clustered into group two (**Figure 4**). Though 'Yuluxiang' is a hybrid cultivar of 'Kuerle' (P. ussuriensis) and 'Xuehua' (P. bretschneideri), the results of this study indicate that the 'Yuluxiang' pear is much more phylogenetically close to the P. sinkiangensis from the hereditary perspective of epicuticular wax. The chemical composition of two major groups (group one and two) and two subgroups (I–II) were compared with t-test. We found that the statistical significant number have 60 and 62 according to two major groups and two subgroups, respectively (Supplementary Table S3).

#### Principal Component Analysis of Wax Compounds

The seven classes of epicuticular wax compounds (alkanes, primary alcohols, fatty acids, aldehydes, esters, terpenoids, and others) found in 35 pear cultivars were analyzed through principal component analysis (PCA). The scatter plot of the scores of the PCs of 35 pear cultivars and the corresponding loadings plot of epicuticular wax compounds are shown in **Figure 5**. PCA1 (47.12%), PCA2 (16.65%) and PCA3 (16.03%) accounted for 79.80% of the total variance, which was high enough to represent all the variables (**Figure 5**). The 35 pear cultivars could be divided into three groups on the basis of the relationships between cultivars (scores) and their seven classes of wax compounds (loadings). Group one was on the positive side of PCA1 with the longest distance and on the negative side of PCA3 and included cultivars 'Jinfeng', 'Tianjianba' and 'Clapp Favorite', which were characterized by high concentrations of alkanes (**Figure 5AI**). Cultivars on the positive side of PCA2 with the longest distance formed the second group and contained five pear cultivars ('Kousui', 'Docteur Jules Guyot', 'Shinseiki', 'Bartlett MaxRed' and 'Abbe Fetel'); these cultivars were characterized as cultivars with lower concentrations of terpenoids (**Figure 5AII**). The third group contained all the other cultivars, and there was no consistency with respect to the composition of wax compounds (**Figure 5AIII**).

To further understand how pear cultivars were affected by cuticular wax compositions, we performed the PCA analysis using 146 wax compounds. Compared with the former PCA analysis using seven classes of epicuticular wax compounds, different similar result showed that the PCA1 (64.45%), PCA2 (8.12%) and PCA3 (6.40%), explained 78.97% of the variation (Supplementary Figure S6). However, the 35 pear cultivars could be divided into five groups based on the relationships between cultivars (scores) and their 146 wax compounds (loadings). 'Clapp Favorite' formed the first group, which were characterized by high concentrations of octacosane. Group two contained six pear cultivars ('Bartlett Max Red', 'Hongnanguo', 'Nanguoli', 'Huagai', 'Xuehuali' and 'Zhongli No. 2'), which were characterized with high concentrations of triacontane-1,30-diol. Group three including two cultivars ('Eli' and 'Huangguan'), which were mainly characterized by high concentrations of triacontane. 'Abbe Fetel', 'Hongxiangsu' and 'Xueqing' characterized with lack of triacontane and hentriacontane were clustered into group four. The fifth group contained all the other cultivars, and these cultivars were characterized with high concentrations of hentriacontane (Supplementary Figure S6A).

# CONCLUSION

The total wax content of mature fruits was highly variable among 35 pear cultivars, ranging from a high of 2.44 ± 0.41 mg/cm<sup>2</sup> ('Docteur Jules Guyot', P. communis) to a low of 0.46 ± 0.03 mg/cm<sup>2</sup> ('Huagai', P. ussuriensis). A total of 146 epicuticular wax compounds were detected through GC-MS analysis. Although the composition and concentration of wax compounds varied greatly, alkanes (422.73 ± 205.79 µg/cm<sup>2</sup> , 40.72 ± 12.03%), primary alcohols (271.06 ± 210.96 µg/cm<sup>2</sup> , 24.48 ± 10.57%) and terpenoids (116.39 ± 98.72 µg/cm<sup>2</sup> , 11.80 ± 8.76%) were the predominant wax compounds in all 35 pear cultivars. FESEM images showed that epicuticular wax crystals of the 35 pear cultivars were mostly amorphous structures, including crystalline plates, irregular ovate crystals, platelets and rodlets. The cluster analysis of wax composition and concentration revealed that the Pyrus bretschneideri cultivars were grouped much closer to Pyrus pyrifolia and Pyrus ussuriensis, whereas Pyrus sinkiangensis were clustered into a distant separate group. The 35 pear cultivars could be divided into three groups and five groups based on the principal component analysis (PCA) using seven main classes of epicuticular wax compounds and all of the 146 wax compounds, respectively. Taken together, 'Docteur Jules Guyot' (P. communis), 'Tianjianba' (P. ussuriensis), 'Kuerle' (Pyrus sinkiangensis Yü.), 'Qiubai' (P. bretschneideri) and 'Shinseiki' (P. pyrifolia) should be the most representative parent cultivars for the future pear breeding with higher wax concentrations.

#### AUTHOR CONTRIBUTIONS

HY and SZ: conceived and designed the experiments. XW, YC, GW, and PC: performed the experiments. XW and YC: analyzed the data. ZS and KQ: contributed reagents, materials, analysis tools or data. XW, XQ, HY, and SZ: wrote the paper.

# FUNDING

This study was supported by the National Natural Science Foundation of China (31701890), Fundamental Research Funds for the Central Universities (KJQN201818), the Key Project for New Agricultural Cultivar Breeding in Zhejiang Province, China (2016C02052-4), Natural Science Foundation of Jiangsu Province in China (BK20160715), and China Postdoctoral Science Foundation funded project (2015M570456).

#### ACKNOWLEDGMENTS

We would like to thank many contributors of experimental materials, including Qingyu Li (Yantai Academy of Agricultural

Sciences, Yantai, Shandong), Shoufeng Sha (Liaoning Institute of Fruit Sciences, Xiongyue, Liaoning), Yongjie Liu (Research Center of Korla Fragrant Pear, Korla, Xinjiang), Mudan Bai (Pomology Institute, Shanxi Academy of Agricultural Sciences, Taigu, Shanxi), Yingtao Wang (Institute of Fruit Tree Research, Hebei Academy of Agriculture and Forestry Sciences, Shijiazhuang, Hebei), Yufen Cao (National Germplasm Repository of Pear in the Research Institute of Pomology, Chinese Academy of Agricultural Sciences (CAAS), Xingcheng, Liaoning

#### REFERENCES


Province), and Long Wang (Zhengzhou Fruit Research Institute, CAAS, Zhengzhou, Henan, China).

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2018.00679/ full#supplementary-material


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Wu, Yin, Shi, Chen, Qi, Qiao, Wang, Cao and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

fpls-09-00679 May 18, 2018 Time: 16:55 # 14

# Defective APETALA2 Genes Lead to Sepal Modification in Brassica Crops

Yanfeng Zhang1,2† , Shuhua Huang<sup>2</sup>† , Xuefang Wang<sup>1</sup> , Jianwei Liu<sup>2</sup> , Xupeng Guo<sup>1</sup> , Jianxin Mu<sup>1</sup> \*, Jianhua Tian<sup>1</sup> and Xiaofeng Wang<sup>2</sup> \*

<sup>1</sup> Hybrid Rapeseed Research Center of Shaanxi Province, Yangling, China, <sup>2</sup> State Key Laboratory of Crop Stress Biology for Arid Areas, College of Horticulture, Northwest A&F University, Yangling, China

#### Edited by:

Yuke He, Institute of Plant Physiology and Ecology, Shanghai Institutes for Biological Sciences (CAS), China

#### Reviewed by:

Xianzhong Feng, Northeast Institute of Geography and Agroecology (CAS), China Hui Feng, Shenyang Agricultural University, China

#### \*Correspondence:

Xiaofeng Wang wangxff99@nwsuaf.edu.cn Jianxin Mu jxmsxyc@163.com †These authors have contributed equally to this work.

#### Specialty section:

This article was submitted to Plant Physiology, a section of the journal Frontiers in Plant Science

Received: 17 January 2018 Accepted: 06 March 2018 Published: 20 March 2018

#### Citation:

Zhang Y, Huang S, Wang X, Liu J, Guo X, Mu J, Tian J and Wang X (2018) Defective APETALA2 Genes Lead to Sepal Modification in Brassica Crops. Front. Plant Sci. 9:367. doi: 10.3389/fpls.2018.00367 Many vegetable and oilseed crops belong to Brassica species. The seed production of these crops is hampered often by abnormal floral organs, especially under the conditions of abiotic conditions. However, the molecular reasons for these abnormal floral organs remains poorly understood. Here, we report a novel pistil-like flower mutant of B. rapa. In the flower of this mutant, the four sepals are modified to one merged carpel that look like a ring in the sepal positions, enveloping some abnormal stamens and a pistil, and resulting in poor seed production. This novel mutant is named sepalcarpel modification (scm). DNA sequencing showed that the BrAP2a gene, the ortholog of Arabidopsis APETALA2 (AP2) that specifies sepal identity, losses the function of in scm mutant due to a 119-bp repeated sequence insertion that resulted in an early transcription termination. BrAP2b, the paralog of BrAP2a featured two single-nucleotide substitutions that cause a single amino acid substitution in the highly conserved acidic serine-rich transcriptional activation domain. Each of the two BrAP2 genes rescues the sepal defective phenotype of the ap2-5 mutant of Arabidopsis. Furthermore, the knockout mutation of the corresponding BnAP2 genes of oilseed rape (B. napus) by CRISPR/Cas9-mediated genome editing system resulted in scm-like phenotype. These results suggest that BrAP2 gene plays a key role in sepal modification. Our finding provides an insight into molecular mechanism underlying morphological modification of floral organs and is useful for genetic manipulation of flower modification and improvement of seed production of Brassica crops.

Keywords: APETALA2 (AP2), carpeloid sepal, Brassica, flower development, oilseed rape, organ modification

#### INTRODUCTION

Eudicot flowers are composed of four distinct floral organs, sepals, petals, stamens, and carpels, which are arranged in four concentric whorls from the outer to the inner layer of the flower. Genetic analysis of floral homeotic mutants in Arabidopsis thaliana and Antirrhinum majus led to the development of the classic ABC model of floral organ identity (Bowman et al., 1991; Coen and Meyerowitz, 1991). This model proposes three classes of floral homeotic genes that act alone or cooperate with each other to specify the four floral organs: class A specifies the sepals, A and B specify the petals, B and C specify the stamens, and C specifies the carpels. A mutation in the A-function gene generates a homeotic transition of sepals to carpels and petals to stamens. A B-functional gene mutation gives rise to a homeotic transformation of petals to sepals and

**88**

stamens to carpels. A C-functional gene mutation leads to a homeotic translation of stamens to petals and carpels to sepals (Meyerowitz et al., 1991; Weigel and Meyerowitz, 1994).

In Arabidopsis, transcription factor APETALA1 (AP1) and APETALA2 (AP2) are A-functional genes, APETALA3 (AP3) and PISTILLATA (PI) are B-functional genes, while AGAMOUS (AG) is a C-functional gene (Yanofsky et al., 1990; Bowman et al., 1993; Jofuku et al., 1994; Krizek and Meyerowitz, 1996). In addition, it has been shown that SEPALLATA orthologs (SEP1–SEP4) are also required to specify floral organ identity (Pelaz et al., 2000). These genes are defined as E-function genes, and their discovery led to the proposal of the ABCE model (Jack, 2001). All of these well-characterized genes, except for AP2, are members of the MADS-box family encoding MIKC-type MADS-box transcription factors (Saedler et al., 2001; Theissen, 2001; Kaufmann et al., 2005). These encoded proteins constitute four different tetrameric complexes that specify the various corresponding floral organs. Accordingly, the AP1, AP1, SEP, and SEP tetramer determines the sepal development; AP1, AP3, PI, and SEP complex determines the petal development; AP3, PI, AG, and SEP determines the stamen development; while the AG, AG, SEP, and SEP tetramer determines the development of carpel (Yanofsky et al., 1990; Honma and Goto, 2001; Pelaz et al., 2001; Melzer and Theißen, 2009). This protein complexesbased model of floral organ specification is referred to as the quartet model, and is an improved successor of the ABCE theory (Theissen and Saedler, 2001). So far, the ABC model and its improved, more detailed versions have proven to be useful frameworks for interpreting floral development in a variety of species. However, there is little evidence to support the fact that the A-function genes specify the development of sepals and petals in non-Arabidopsis plant species.

In Arabidopsis, both A-functional genes AP1 and AP2 are required for sepal and petal development (Kunst et al., 1989; Irish and Sussex, 1990; Drews et al., 1991; Bowman et al., 1993; Gustafson-Brown et al., 1994). Indeed, Arabidopsis is the only species in which floral organ identity defects have been linked to A-functional genes (Davies et al., 1999; Litt, 2007; Causier et al., 2010). Interestingly, the SQUAMOSA (SQUA) gene, an AP1 ortholog in Antirrhinum majus, plays a key role in determining floral meristem identity, but does not affect the specification of sepals and petals in this species (Huijser et al., 1992). Meanwhile, AP2 was originally identified in Arabidopsis as a floral homeotic gene involved both in the A-functional roles of specifying organ identity and in restricted C-functional activities. Loss of AP2 function has previously been shown to cause a homeotic transition of sepals to carpels, and petals to stamens (Drews et al., 1991; Bowman et al., 1993; Causier et al., 2010). Further investigations revealed that the AP2 gene has numerous additional functions in the regulation of flowering time (Aukerman and Sakai, 2003), maintenance of the shoot apical meristem (Wurschum et al., 2006), and the determination of seed and fruit development (Maes et al., 2001; Jofuku et al., 2005; Ohto et al., 2005; Ripoll et al., 2011; Zhou et al., 2012). However, the exact role of the AP2 gene in flower development, especially in sepal and petal specification, has not been determined. Therefore, exploring of novel A-functional homeotic mutants in non-Arabidopsis species is currently critical to verify the functions of the AP1 and AP2 genes in sepal and petal specification, which could contribute new details to complement the classic ABC model of floral organ identity.

Many vegetable and oilseed crops belong to Brassica species. The seed production of these crops is hampered often by abnormal floral organs, especially under conditions of abiotic stress (Zhou et al., 2008; Huang et al., 2010). Similar to Arabidopsis, both B. rapa and B. napus are members of the Cruciferae family. B. napus is an allodiploid species derived from spontaneous hybridization between B. rapa and B. oleracea, and possesses the complete diploid chromosome sets of the parental genomes (Parkin et al., 1995; Snowdon et al., 2002). In our previous study, we recovered a novel pistil-like flower mutant with sepal-to-carpel modification (scm) in an F<sup>3</sup> selfing line of an interspecific cross between B. napus and B. rapa, which prompted us to further explore the molecular mechanism underlying its phenotype. Here we report that the scm mutant was apetalous resulting in poor seed production. Further studies indicated that the mutations in the two AP2 orthologs in B. rapa (BrAP2 genes) genes of B. rapa (Song et al., 2013) were the fundamental causes of the defects.

# MATERIALS AND METHODS

#### Plant Materials

The scm and scm2 from the same F<sup>3</sup> selfing line Nc116 were derived by an interspecific cross between B. napus and B. rapa. Both wild-type B. rapa (WTr) and B. napus (WTn) were used as controls. The AP2 weak mutant ap2-5 of Arabidopsis thaliana (stock name CS6239, purchased from the website<sup>1</sup> ) with a mild sepal-to-carpel and petal-to-stamen phenotype was used as transgenic receptor material.

# RNA Extraction and Gene Expression Assays

Total RNA was extracted from young inflorescence at the budding stage using RNAiso Reagent (TaKaRa, Dalian, China). The first strand cDNA was synthesized using the Prime Script RT-PCR kit (TaKaRa, Dalian, China). Primers specific for the AP1, AP2, AP3, PI, AG, and SEP3 genes were designed based on sequences retrieved from the Brassica database<sup>2</sup> and the NCBI nucleotide and EST databases (details of primers are given in Supplementary Table 1). The expression of these genes was determined using semi-quantitative RT-PCR. The expression levels of the target genes were normalized to the internal control, the cytosolic 18S rRNA gene. The PCR products were separated on a 1% agarose gel, EtBr stained, purified using the Omega DNA Gel Extraction Kit (Omega Bio-Tek), and cloned into the pGEM-T Easy Vector (Promega, Madison, WI, United States) for sequencing (Beijing Sunbiotech Co., Ltd.).

<sup>1</sup>www.arabidopsis.org

<sup>2</sup>http://brassicadb.org/brad/

# Generation of Transgenic Arabidopsis Lines of the Two B. rapa AP2 Genes

To produce GFP-tagged AP2 genes from B. rapa (Br) and Arabidopsis (At), the coding sequences (CDSs) without the stop codon of the BrAP2a, BrAP2b, and AtAP2 (control) were PCR-amplified using the following primer sets: BAP2-F+Sac I/BAP2-R2a+Xba I, BAP2-F+Sac I/BAP2-R2b+Xba I, and AtAP2-F+Kpn I/AtAP2-R+Xba I, respectively (Supplementary Table 1). The amplicons were checked by sequencing, and the full-length CDSs were then cloned into the Sac I and Xba I sites of Cam-35S-GFP-containing vectors, resulting in the C-terminal fusion to GFP.

The resultant Cam-35S-BrAP2a-GFP, Cam-35S-BrAP2b-GFP, and Cam-35S-AtAP2-GFP fusion constructs were transformed into the Agrobacterium tumefaciens strain GV3101 using electroporation and subsequently infiltrated into Arabidopsis ap2-5 and Col-0 plants are using the floral dip method (Clough and Bent, 1998). T1 seeds from the infiltrated plants were selected and plated on 1/2 MS medium containing 33 mg/L hygromycin B.

# Western Blot Assay

To examine the expression of the AtAP2-GFP, BrAP2a-GFP, and BrAP2b-GFP genes in independent transgenic lines, 50 mg of young leaves were ground in liquid N<sup>2</sup> and treated with 200 µL 1×SDS loading buffer without Bromophenol Blue (50 mM pH 6.8 Tris–HCl, 5% β-mercaptoethanol, 10% Glycerol, and 2% SDS), followed by further grinding. The homogenate was boiled for 5 min and centrifuged for 5 min at 10,000 g. The supernatant was collected, and its protein concentration was determined using a Bradford assay. Samples containing 50 µg of protein were subjected to 10% SDS–PAGE and transferred to a PVDF membrane (Roche, Mannheim, Germany). Immunoblot analysis was performed using an anti-GFP (1:2000) monoclonal antibody (Beijing TransGen Biotech Co., Ltd.). The chemiluminescent HRP substrate (Millipore) was then added, incubated for 5 min, and the blots were imaged using Image Lab Software (Bio-Rad, Hercules, CA, United States).

### Generation of ap2 Quadruple Mutants of B. napus by CRISPR/Cas9-Mediated Genome Editing System

To construct the desired recombinant plasmids using CRISPR/Cas9-based system, the same target sequence of AP2 for the four cDNAs of B. napus was designed using the web tool CRISPR-P<sup>3</sup> . The target dsDNA molecules with sticky ends were generated via direct annealing of two oligonucleotides primers in which the sticky ends sequences match those of the Bsa I ends in the pHSN401 vector (Xing et al., 2014). The resultant pHSN401-AP2 constructs were transformed into the Agrobacterium tumefaciens strain GV3101:pSoup using electroporation, and subsequently transformed into B. napus by the previously described method (Bhalla and Singh, 2008). The genomic DNA extracted from T1 plants was used for PCR

<sup>3</sup>http://crispr.hzau.edu.cn/CRISPR/

amplification with primer sets AP2-F250/AP2-R507, and the PCR products were cloned into pMD19-T Vector (Takara) for sequencing using M13-47 primer (Beijing Sunbiotech Co., Ltd.).

# RESULTS

### Origin and Floral Identities of scm, a Novel Flower Mutant of B. rapa

During our breeding experiments, we recovered a novel pistillike flower mutant in the F<sup>3</sup> selfing line Nc116, which was generated from an interspecific cross between B. napus and B. rapa. Unlike the wild-type Brassica (**Figures 1A,D,G,I,L1,L2**, WT), the flowers of the mutant consisted of a very large pistil and lacked sepals, petals, or stamens (**Figures 1B,E1**).

FIGURE 1 | Characteristics of the sepal carpeloid mutants of Brassica. (A–C): inflorescence; (D,E1,E2,F1–F3): flower; (G,H): floral organs; (I–K): SEM photos of flowers; (L1,L2,M1,M2,N1,N2): SEM photos of floral organs; (A,D,G,I,L1,L2): wild type Brassica (WT). (L1,L2) Are sepal (left) and stamen (right); (B,E1,E2,H,J,M1,M2): the typical sepal carpeloid mutant (scm). M1 is a sepal carpeloid organ, and M2 is an incomplete stamen to carpel organ with one ovule; (C,F1–F3,K,N1,N2): the incomplete sepal carpeloid mutant (scm2). Arrows indicate sepal carpel-like transition in F2 and petal stamen-like transition in F3. N1 is an incomplete petal to stamen organ, and N2 is a stamen of scm2. Bars in E1 and E2 are 5 mm; bars in I–K,L1,M1 are 1 mm; bars in L2,M2,N1,N2 are 500 µm; and bars in the others are 10 mm.

After 2 days of growth, the mutants showed a small pistil stretching out of the large, initial pistil (**Figure 1E2**). When the large pistil was torn, it appeared to be separated into three parts: the outer layer was a large pistil consisting of four fused carpels with green ovules inside; the intermediate layer only consisted of only 1–4 abnormal stamens with chimeric carpel-like homeotic identity; and the inner layer was a normal pistil (**Figures 1H,J,M1,M2**). Surprisingly, it appeared that the four whorls and four types floral organs (sepals, petals, stamens, and carpels) found in normal flowers had transformed into only three whorls containing just two types of organs (carpels, stamens, and carpels, respectively), with the petals being completely abolished. Clearly, the mutant had severe defects in both the sepal and petal whorls; the ring of four merged carpels developed at the four sepal positions representing a typical homeotic transition of sepals to carpels, and the complete loss of the four petals indicating a typical apetalous trait. Evidently, the mutant had typical sepal carpeloid and apetalous characteristics, which we therefore designated as scm. Due to the reduced number of stamens in the scm flower compared to wild-type Brassica flowers, we could not conclude whether scm had undergone a homeotic transition of petals to stamens.

We also obtained another incomplete sepal carpeloid mutant, named scm2, in the same population. Compared to scm, this mutant lost the large pistil shape (**Figures 1C,K**), and only a few sepals developed into carpels (**Figures 1F1,F2**). Moreover, some flowers of scm2 were also apetalous, some of which containing few small abnormal petals, and some showing stamen-like petals (**Figures 1F3,N1,N2**). These results indicated that scm2 was an incomplete mutant with the homeotic transitions of sepals to carpels and petals to stamens.

#### Expression of ABCE Class of Genes in the scm Mutant

According to the classical ABC model, floral organ identity genes, especially those of class A and B, specify sepal and petal development. Mutations in the A class gene AP2 in Arabidopsis lead to homeotic transitions of sepals to carpels and petals to stamens (Kunst et al., 1989; Coen and Meyerowitz, 1991). Consistent with the previous findings, our mutants displayed a similar homeotic transition. This prompted us to examine the AP2 gene expression pattern, as well as the expression of the other ABCE classes of genes, in scm and in wild-type B. rapa (WTr) and B. napus (WTn).

Thus, the expression of AP1, AP2, AP3, PI, AG, and SEP3 genes in these samples was studied by semi-quantitative RT-PCR. When AP1-specific primers were applied, all the three samples formed two identical PCR amplicons (**Figure 2**, top left). By contrast, when BrAP2-specific primers were used, both WTr and WTn each produced one clear band for amplicons of identical size, whereas scm formed two bands (**Figure 2**, top right). The length of one of the bands was similar to both controls, while the other band was slightly larger than its counterpart. This larger fragment was named SCM-a, and the smaller fragment was defined as SCM-b. Using AP3 primers, we detected fragments of the same size in the WTr and scm samples, but an additional

larger band appeared in the gel of the WTn sample (**Figure 2**, second from the top left). Meanwhile, the amplification with PI, AG, and SEP3 primers resulted in the production of identical amplicons in all the three samples, with no significant differences (**Figure 2**). These data indicate that a Brassica BrAP2 mutation, like the sepal carpeloid AP2-gene mutant in Arabidopsis, is probably relevant to the scm mutant. Moreover, the results from the AP3 gene amplification indicated that the genetic background of scm was very similar to that of its B. rapa parent.

### The Two AP2 Orthologs in B. rapa Are Mutated

To explore the AP2 orthologs in B. rapa, a BLAST search was performed in the Brassica database, and two genes, Bra011741 and Bra017809, were identified. Bra011741 is located on chromosome A01 between base pairs 00793728 and 00795868, while Bra017809 is found on the chromosome A03 between base pairs 30612490 and 30614651. We also detected two AP2 orthologs in B. oleracea, Bol028934 and Bol018627, which are located on chromosomes C01 and C06, respectively (Supplementary Figure 1). Our results showed that B. rapa has two BrAP2 genes, whereas B. napus might have four BnAP2 genes, two of which from B. rapa and two from B. oleracea (Parkin et al., 1995; Snowdon et al., 2002).

According to the BrAP2 sequences from B. rapa and B. oleracea, primers spanning the full length cDNA of the two BrAP2 genes were designed: BrAP2-F and BrAP2-R2a target to the Bra011741 and Bol028394 genes; and BrAP2-F and BrAP2- R2b specific for the Bra017809 and Bol018627 genes. PCR products obtained from WTr samples using AP2-F/R2a primers were 1302-bp in length, and were named BrAP2a, sharing 99.4% nucleotide identity and 99.1% amino acid identity with Bra011741. Using the same AP2-F/R2a primer pair, the amplified genes from scm were 100% identical to BrAP2a except for a 119-bp insertion, which was a repeat of the sequence from nucleotides 485–603 (**Figure 3A** and Supplementary Figure 2).

The amplicons obtained from WTr samples using the BrAP2- F/R2b primer set were identified as a 1299-bp fragment named BrAP2b. This gene was 100% identical to Bra017809. The corresponding sequences from scm were 100% identical to BrAP2b except for two single-nucleotide substitutions, among which only one A to C mutation, translating into the alteration of the 18th glutamic acid of the corresponding protein into aspartic acid, was observed (**Figure 3A** and Supplementary Figure 3 E18D). Among all the tested AP2 sequences obtained from scm, we could not detect any AP2 genes derived from B. oleracea, thus confirming that the genetic background of scm was very similar to its B. rapa parent.

primers in the scm and wild-type B. rapa (WTr) sequences are shown (upper part). PCR amplification of genome DNA and cDNA of WTr and the scm using

primer pairs of AP2-F/R, AP2-F/R2, and AP2-F2/R (lower part).

The amplification of genomic DNA using the AP2-F/R2a primer set revealed that the BrAP2a gene of scm had an additional 291-bp fragment, which was a repeat of BrAP2a nucleotides 648–938, inserted after nucleotide 938. The above-mentioned 119-bp cDNA insertion was derived from transcription of this 291-bp DNA repeat sequence (**Figure 3A** and Supplementary Figure 2). Further sequence analysis indicated that this insertion led to an early translation termination of the BrAP2a mRNA.

We then designed additional specific primers to further confirm the above results. Using the AP2-F/R primer set, both DNA and cDNA of scm displayed an additional band compared to WTr (**Figure 3B**). When the BrAP2a specific primer set AP2- F/R2a was applied, both DNA and cDNA of scm produced one band that was larger than that of WTr. Using a mutual primer set (AP2-F2/R) for both of the two AP2 genes, both DNA and cDNA of scm formed one amplicon, the band for which was the same as that of the WTr-derived product. These results showed that the BrAP2a gene of scm underwent an insertion mutation that led to the loss of function of this gene. Altogether, these results suggested that mutations in BrAP2a and BrAP2b are most likely associated with the sepal carpeloid mutant scm. However, it was still unclear whether the loss of function of the two AP2 genes was the key cause of the observed scm phenotype.

#### The Two AP2 Genes of B. rapa Can Both Rescue the ap2-5 Mutant of Arabidopsis

The loss of function of AP2 in Arabidopsis has been previously demonstrated to cause a homeotic transition of sepals to carpels and petals to stamens (Kunst et al., 1989). Moreover, a weak AP2-defective Arabidopsis mutant ap2-5 shows a few sepal-tocarpel and petal-to-stamen homeotic transitions (Kunst et al., 1989). Both AP2 genes in B. rapa, BrAP2a and BrAP2b, shared high similarity with the Arabidopsis AP2 gene in their protein coding sequences (Supplementary Figure 3). To clarify whether BrAP2a and/or BrAP2b have biologically significant roles in floral organ specification, we transformed the two B. rapa AP2 genes and Arabidopsis AP2 (AtAP2) gene (control) separately into the wild-type and ap2-5 mutant Arabidopsis.

All the three AP2 genes were fused to green fluorescent protein (GFP) under the control of constitutive 35S promoter and transferred into Arabidopsis Col-0 to generate p35S::AtAP2/Col, p35S::BrAP2a/Col, and p35S::BrAP2b/Col over-expressing transgenic lines. These three transgenic lines displayed normal flower phenotypes (**Figure 4A**), suggesting that over-expression of the AP2 gene does not affect floral organ development.

When the three AP2 fusion constructs were transformed into the ap2-5 mutant of Arabidopsis, the transgenic lines p35S::AtAP2/ap2-5, p35S::BrAP2a/ap2-5, and p35S::BrAP2b/ap2- 5 displayed normal development of the sepals and petals (**Figure 4B**), and thus the AP2 defect was rescued. Moreover, we recovered some strong sepal carpeloid mutants not only in the three knockdown transgenic groups: p35S::AtAP2/ap2-5 KD, p35S::BrAP2a/ap2-5 KD, and p35S::BrAP2b/ap2-5 KD (**Figure 4C**), but also in the three over-expressing groups caused by a cosuppression effect (data not shown). These mutants had some common characteristics: the sepals were transformed into carpels, the petals were absent, and the number of stamens was reduced to 1–4, but the pistil showed normal development.

These lines displayed a strong AP2-knockdown phenotype, which is most likely also observed in scm. We further carried out RT-PCR and Western blot assays on these transgenic lines. Our RT-PCR results confirmed that the three AP2 rescued lines exhibited high expression whereas the three knockdown phenotypes showed low expression level of the AP2 gene. Interestingly, strong protein bands were present in the Western blot analyses of all the three AP2 rescued lines but could not be detected for the three knock-down lines (**Figure 4D**). Therefore, we inferred that the exogenous AP2 gene could interfere with the expression of endogenous AtAP2 gene through RNAi effect. These results confirmed that the AP2 gene is involved in the regulation of sepal and petal development, and the knockdown of these genes could affect sepal and petal development to generate sepal carpeloid and apetalous mutants.

### The ap2 Quadruple Mutants of B. napus Produced by CRISPR/Cas9 System Resulting in Typical Sepal Carpeloid Phenotype

To further confirm that the AP2 genes are involved in regulating sepal and petal development, a CRISPR/Cas9-mediated genome

editing experiment was performed to test whether the loss of AP2 function in B. napus could lead to a similar sepal carpeloid phenotype as in scm. We designed and constructed CRISPR/Cas9 nuclease system to cut at four identical target sites in the conserved sequences of AP2 (**Figure 5A** and Supplementary Figure 4A).

After that, we performed Agrobacterium-mediated transformation to introduce CRISPR/Cas9 construct into

B. napus to generate 33 transgenic plants (Yan et al., 2012). Among these, some plants showed typical sepal carpeloid phenotypes and some exhibited weak carpeloid phenotypes (**Figure 5B**). Furthermore, the genetic composition of one typical mutant was identified by sequencing, in which five or four nucleotides were deleted in each single chromosome of BnaC-AP2a gene; a "G" nucleotide was inserted into BnaA-AP2a; a "T" nucleotide was inserted or deleted in each single chromosome of BnaA-AP2b gene; and a "A" nucleotide was inserted into BnaC-AP2b (**Figure 5C** and Supplementary Figure 4B). In this way, all of the four AP2 genes in B. napus acquired reading frameshift mutations and early translation termination, reflecting the fact that the four AP2 genes completely lost their functions. The resulting ap2 quadruple mutant exhibited carpels sepals, missing petals, and a reduced number of stamens, but showed a normal pistil, which are identical phenotypes to those displayed by scm (**Figure 1B**). These results suggest that knocking out the AP2 genes could severely affect floral organ development and lead to the generation of a sepal carpeloid and apetalous mutant in Brassica.

#### DISCUSSION

#### The Novel Pistil-Like Flower Mutant Is a Sepal Carpeloid and Apetalous Mutant

The scm displayed some distinct characteristics: the sepals and petals were absent; the four merged carpels that looked like a very large pistil developed in the sepal positions; stamens reduced in numbers to 1–4; and only the pistil developed normally. Clearly, the mutant underwent two substantial changes with respect to its floral organs: one was a typical homeotic transformation of sepals into carpels, and the other was an apetalous phenotype. Thus, the pistil-like flower mutant was not only a typical sepal carpeloid mutant but also an apetalous mutant.

Moreover, we found some petals undergoing transformation into stamens in scm2 indicating that the scm2 underwent an incomplete homeotic transition of sepals to carpels and petals to stamens. It has been well-documented in Arabidopsis that the defective A-function AP2 gene affects sepal and petal development and leads to the homeotic transformation of sepals to carpels and petals to stamens (Kunst et al., 1989). However, our weak and strong mutants presented different phenotypes in which the weak mutant exhibited an incomplete transformation of sepals to carpels and petals to stamens, whereas the strong mutant exhibited a complete sepal carpeloid transition and an apetalous trait. These results suggest that the AP2 strong mutant generates only the sepal-carpeloid transformation and the missing petals trait, but could not produce the petal-to-stamen transformation. Although the phenotypes of the two mutants were not identical, both had severe sepal and petal developmental defects.

#### Defective AP2 Genes Are the Main Cause of the scm Phenotype

Our data revealed that the scm mutant, derived from the F<sup>3</sup> generation of an interspecific cross between B. napus and B. rapa, contains only two, B. rapa-similar AP2 genes, suggesting that its genetic background is close to that of the parent B. rapa. Both the AP2 genes in scm were mutated: one contained a partial repeat sequence that led to an early translation termination signal, and the other included an A to C single-nucleotide substitution that transformed the 18th glutamic acid of the corresponding protein into an aspartic acid. This mutation occurred in a highly conserved acidic serine-rich region (amino acids 14–50) that may serve as a transcriptional activation domain (Jofuku et al., 1994). Overall, the mutation of both AP2 genes in scm most likely resulted in the loss of function of these genes. Correspondingly, Yan et al. (2012) showed that knocking out the BnAP2 gene in B. napus generated a sepal carpeloid and apetalous mutant as well as some incomplete mutants. Similarly, we also generated several AP2-knockdown plants in the AP2-overexpressing transgenic Arabidopsis lines, as well as a number of AP2-knockout B. napus transgenic plants with a sepal carpeloid and apetalous phenotype. These results further confirmed that loss-of-function mutations of the AP2 genes resulted in the sepal carpeloid and apetalous phenotypes.

Our transgenic results demonstrated that the two normal AP2 genes from B. rapa can both rescue the sepal and petal developmental defects of the ap2-5 mutant of Arabidopsis, indicating that both normal AP2 genes have the ability to regulate sepal and carpel development, and reflecting that mutation of neither BrAP2a nor BrAP2b alone can lead to the scm phenotype.

Taken together, our results confirm that the AP2 gene plays a key role in regulating sepal and petal development and that its function is highly conserved in Cruciferae family, many of which are important crop plants. The complete loss of AP2 function severely affects sepal and petal development and generates sepal carpeloid and apetalous mutants. The two defective AP2 genes of B. rapa were the core cause for the sepal carpeloid mutant scm phenotype. However, the molecular mechanism of floral development regulated by AP2 gene is still unclear, and needs to be investigated further.

#### AUTHOR CONTRIBUTIONS

XiW and YZ designed the experiments. YZ, SH, JL, JM, and XG performed the experiments. YZ and SH wrote and revised the article. XuW and JT bred the scm mutant. All the authors gave the final approval for submission of the manuscript.

#### FUNDING

This work was supported by a grant from National Key R&D Program of China (2016YFD0101900), National Science Foundation of Shaanxi Province (2014JQ3087), and the Ph.D. Programs Foundation of China (No. 2011M500441).

#### ACKNOWLEDGMENTS

fpls-09-00367 March 17, 2018 Time: 18:31 # 9

We thank Dr. Weirong Xu at Ningxia University for comments on the manuscript and Dr. Qijun Chen at China Agricultural University for the kind gift of CRISPR/Cas9 vector pHSN401.

#### REFERENCES


#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2018.00367/ full#supplementary-material


gene agamous resembles transcription factors. Nature 346, 35–39. doi: 10.1038/ 346035a0


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Zhang, Huang, Wang, Liu, Guo, Mu, Tian and Wang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.