Tomato locule number and fruit size controlled by natural alleles of lc and fas

Abstract Improving yield by increasing the size of produce is an important selection criterion during the domestication of fruit and vegetable crops. Genes controlling meristem organization and organ formation work in concert to regulate the size of reproductive organs. In tomato, lc and fas control locule number, which often leads to enlarged fruits compared to the wild progenitors. LC is encoded by the tomato ortholog of WUSCHEL (WUS), whereas FAS is encoded by the tomato ortholog of CLAVATA3 (CLV3). The critical role of the WUS‐CLV3 feedback loop in meristem organization has been demonstrated in several plant species. We show that mutant alleles for both loci in tomato led to an expansion of the SlWUS expression domain in young floral buds 2–3 days after initiation. Single and double mutant alleles of lc and fas maintain higher SlWUS expression during the development of the carpel primordia in the floral bud. This augmentation and altered spatial expression of SlWUS provided a mechanistic basis for the formation of multilocular and large fruits. Our results indicated that lc and fas are gain‐of‐function and partially loss‐of‐function alleles, respectively, while both mutations positively affect the size of tomato floral meristems. In addition, expression profiling showed that lc and fas affected the expression of several genes in biological processes including those involved in meristem/flower development, patterning, microtubule binding activity, and sterol biosynthesis. Several differentially expressed genes co‐expressed with SlWUS have been identified, and they are enriched for functions in meristem regulation. Our results provide new insights into the transcriptional regulation of genes that modulate meristem maintenance and floral organ determinacy in tomato.


| INTRODUC TI ON
To develop varieties with improved characteristics, breeders exploit the genetic variation of important agronomic traits such as fruit weight, grain yield, and overall plant architecture. With the advent of Quantitative Trait Loci (QTL) mapping and cloning methods in past decades, many genes contributing to the increase in fruit weight and crop yield have been identified (Bommert, Nagasawa, & Jackson, 2013;Chakrabarti et al., 2013;Frary et al., 2000;Li et al., 2011;Song et al., 2007). The weight of fruits and vegetables is genetically controlled as early as inflorescence and floral meristem development (van der Knaap & Østergaard, 2017). Specifically, the organization of floral meristems is important in determining the final number of carpels, which collectively become the fruit after fertilization of the ovules (Fletcher, Brand, Running, Simon, & Meyerowitz, 1999;Mayer et al., 1998;Xu et al., 2015). Not only is fruit weight affected, the differential regulation of floral meristem development also accounts for variances in rice grain size and maize kernel number (Bommert et al., 2013;Suzaki et al., 2009). Together, these findings indicate that meristem-regulated processes are essential and conserved features for reproductive organ development in higher plants.
Subsequently, AG suppresses WUS expression by binding to the CArG box located downstream of the gene and the recruitment of polycomb group proteins, known to trigger transcriptional repression by enhancing histone methylation ). An ag loss-of-function mutation completely abolishes carpel development in Arabidopsis (Yanofsky et al., 1990).
LC and FAS are two important genes contributing to enlarged fruits with many locules in tomato. Both mutants were selected among different genetic subgroups during tomato domestication because of their positive effects on fruit weight (Blanca et al., 2015;Rodríguez et al., 2011). However, a mutation in FAS often results in unacceptable fruits that are unevenly shaped and therefore, the allele is not commonly found in conventional and commercially grown tomatoes. lc is a mutation near SlWUS, and is expected to cause increased expression by abolishing the binding site of its suppressor AGAMOUS. The causative mutation is comprised of two SNPs that are downstream of the 3′ UTR of SlWUS . fas on the other hand is caused by a ~294 kb inversion with a breakpoint in the promoter region of SlCLV3 (Huang & van der Knaap, 2011;Xu et al., 2015). Although the role of CLV3-WUS in meristem maintenance has been investigated in Arabidopsis and other plant species, it remains unclear whether novel factors are involved in FM regulation in tomato. To understand the molecular mechanism underpinning LC-FAS mediated developmental processes, we conducted a series of genetic and gene expression analyses using backcross populations ( Figure S1). We show that lc is a gain-of-function mutation of SlWUS, whereas fas is a partial loss-of-function mutation of SlCLV3. Importantly, our RNA-seq results led to the identification of a number of co-expressed gene clusters associated with various developmental processes that might act downstream of LC and FAS.

| The effects of lc and fas on fruit morphology and reproductive traits
The natural mutations at the lc and fas loci underlie the orthologs of the Arabidopsis meristem organization genes WUS and CLV3, and were originally identified due to their strong effects on locule F I G U R E 1 The effect of natural lc and fas mutant alleles on floral organs and inflorescence development. (a) Tomato varieties containing lc and/or fas mutant alleles carry multilocular fruits. The wild type tomato (Solanum pimpinellifolium LA1589) typically contains only two locules. Size bar = 3 cm. (b) Genomic sequence changes in lc and fas. The genomic sequence underlying the lc mutation shares similarity to the CArG box of Arabidopsis AP3. The two SNPs underlying lc are marked in red and the putative CArG box is highlighted in gray. The fas mutation is caused by a ~294 kb inversion with a breakpoint in the promoter region of SlCLV3. (c-e) Inflorescences, flowers and fruits of lc, fas, lc/fas NILs and SlCLV3-RNAi lines. Bar = 1 cm. (f-i) The ratio of branched inflorescences and floral organ number in NILs and various transgenic lines. RNAi-CR4 and RNAi-CR9 represent two independent SlCLV3-RNAi transgenic lines. pHC2-6-2 and pHC2-7-2 represent two independent transgenic lines that were transformed with SlCLV3 genomic sequence driven by a 5.5 kb promoter construct. pHC4 contained a shorter SlCLV3 promoter that served as a negative control for the complementation test. Pairwise comparisons were made between the genotypes using ANOVA and means were separated with Tukey's HSD test with p < 0.05 number ( Figure 1a) (Barrero, Cong, Wu, & Tanksley, 2006;Lippman & Tanksley, 2001;Muños et al., 2011;Rodríguez et al., 2011;Xu et al., 2015). The associated nucleotide polymorphisms are 3′ of SlWUS van der Knaap et al., 2014) and the promoter region of SlCLV3 (Huang & van der Knaap, 2011;Xu et al., 2015), respectively ( Figure 1b). The lc mutations may correspond to the CArG box, which is critical to suppress WUS expression in Arabidopsis . We determined the effects of lc and fas on fruit morphology and reproductive traits by creating near-isogenic lines (NILs) in the wild species Solanum pimpinellifolium accession LA1589 background that differed only for the alleles at these two loci. In addition to more locules, the expectation is that these meristem organization mutants would lead to increased inflorescence branching (Park, Jiang, Schatz, & Lippman, 2012). Wild type tomato typically develop a single-branch inflorescence with a reiterating pattern of an IM terminating into an FM and the formation of a new IM along the flank of the FM (Figure S1b), whereas fas and especially lc/fas NILs show a significant increase in inflorescence branching (Figure 1c,f). For example, the first inflorescence in the fas and lc/fas NILs nearly always forms a branched architecture -as two IMs appear to emerge from a single terminating SAM ( Figure S2). lc/fas also resulted in the highest floral organ number among the NILs especially for locule number ( Figure 1d-e,g-i). Overall, lc alone had no effect on inflorescence branching and floral organ number unless in combination with fas (Table S1). Regardless, the most dramatic floral organ number change was for locule number in the natural lc and fas mutants.
SlCLV3 has been demonstrated to underlie fas (Xu et al., 2015), which implies that the fas inversion has compromised the promoter of this gene. We also downregulated SlCLV3 expression by expressing part of the coding region as an RNAi construct in stably transformed tomato. As expected, the tomato plants transformed with the pHC2 construct, which contained approximately 5 kb of the wild type promoter and the entire coding region of the gene, rescued the bi-locular fruit phenotype in the fas background, whereas the shorter promoter construct pHC4 did not (Figure 1i). On the other hand, downregulation of SlCLV3 led to an increase in all floral organs as well as inflorescence branching (Figure 1c-i). Additional phenotypes associated with the severely reduced expression of SlCLV3 included nearly seedless fruits, indeterminate meristematic activity as evidenced by the development of ovaries within the initial ovary, widened leaflets with smooth margins, decreased number of secondary leaflets, and reduced complexity of the compound leaf (Figure 2a-b). Extreme phenotypes in the flowers were also occasionally observed, such as an inflorescence that was reinitiated inside a flower ( Figure 2c).
To determine whether the selection of lc and fas during domestication might have been due to the associated increases in fruit weight, we evaluated the effect of the natural mutations on this trait. In the wild species background, fruit weight was only significantly increased in the lc/fas double NIL (Table 1a). Fruit area was significantly increased in fas and even more in the lc/fas double NIL. This suggested that the fruits were wider but flatter and thus the increase in locule number did not yield much heavier fruits in the LA1589 background. To determine whether lc and fas exerted a synergistic effect on locule number, we evaluated the genetic effect between these two loci (Table 1b). The epistatic analysis was performed using two-way ANOVA and showed a significant interaction between lc and fas for locule number (p-value < 0.001) but not for other traits (Table 1b). The degree of dominance of each locus showed that lc affected locule number in a mostly additive manner. On the other hand, the fas mutation was nearly completely recessive over the wild type with a d/a value of −0.88 (Table S2).

| lc and fas produce fasciated inflorescences and enlarged meristems
To determine whether higher locule number was associated with However, almost all detected transcripts mapped to the 5′ region of SlYABBY2b comprising the first exon and before the breakpoint of the fas inversion ( Figure S3). This indicated that the mutation led to the production of truncated RNAs that are likely non-functional. We wanted to evaluate whether the WUS-dependent AG expression found in Arabidopsis is conserved in tomato. The expression of tomato AG homolog TAG1 was significantly increased in lc, fas, lc/fas floral buds at stage 4 dpi ( Figure 4f). These differences were also observed in mutants at 6 dpi albeit to a lesser extent than at 4 dpi. Notably, the synergistic effect of lc and fas was not observed from TAG1 expression, suggesting that SlWUS solely controls TAG1 expression. It is conceivable that TAG1 leads to downregulation of SlWUS though its binding to the CArG box located 3′ of SlWUS. However, the increased expression of TAG1 in lc, fas and lc/fas NILs was not associated with reduced SlWUS expression suggesting that TAG1 might not directly involved in SlWUS repression.

| Pronounced temporal-spatial expression changes of SlCLV3 and SlWUS in the lc and fas NILs
To determine whether temporal-spatial expression patterns of SlCLV3 and SlWUS were altered in the NILs, we performed RNA in situ hybridi-

| Differentially expressed genes in lc and fas
To decipher the genome-wide gene expression changes resulting from mutations in lc and fas, we applied linear factorial analysis on the RNA-seq results across the five developmental stages.
Differentially expressed genes (DEGs) that were consistently up-or down-regulated across different stages were identified as genes with significant genotype effect, while DEGs responding differently to genotypic and developmental variations (ex: showing dynamic expression patterns across stages) were considered as genes with significant genotype by developmental We also identified genes encoding metabolic enzymes such as glucose-6-phosphate isomerase (Solyc04 g076090), ALCOHOL DEHYDROGENASE 2/SlADH2 (Solyc06 g059740) involved in fatty acid degradation and the production of volatiles during fruit ripening (Speirs et al., 1998), and a close homolog of Arabidopsis F I G U R E 4 RNA-seq analysis of SlCLV3, SlWUS, SlYABBY2 and TAG1 during floral development. Tissues were collected from lc, fas, lc/fas, RNAi-SlCLV3, and the wild type plants at five developmental stages: sympodial shoot apical meristem (SYM), floral meristem with inflorescence meristem (F&IM), 2, 4 and 6 dpi. Shown are the normalized expression of SlCLV3 (a, c) and SlWUS (b, d), SlYABBY2b (e) and TAG1 (f). The expression levels obtained from 3′ Tag RNA-seq method (a, b, e, f) were normalized using reads per million reads (RPM), while data obtained from whole mRNA-seq method (c, d) were normalized using reads per kilobase million reads (RPKM). The p-value was obtained from linear-based likelihood ratio test between mutants and the wild type using DEseq2 in R. Data are shown as means ± SD from three to four biological replicates. Significant differences are represented by asterisks. •p < 0.1, *p < 0.05, **p < 0.001 and ***p < 0.0001 (a, b, e, f). * adjusted p < 0.05 (c, d) which affects anthocyanin accumulation (Kubo, Peeters, Aarts, Pereira, & Koornneef, 1999). DEGs with G × D interaction effects also included a NON-PHOTOTROPIC HYPOCOTYL 3

| Cluster analysis and DEGs co-expressed with SlWUS
The linear factorial modeling identified a large collection of DEGs that were consistently up-or down-regulated in all developmental stages as well as those showing genotype-and developmental stage-dependent differential expressions. To identify groups of DEGs with similar expression dynamics among the 675 DEGs dataset, we conducted fuzzy C-means clustering using corresponding genes with normalized expression levels in the WT. This led to the identification of eight co-expressed clusters with cluster 1 representing SlWUS and SlCLV3 (Figure 7a). Clusters with DEGs that F I G U R E 6 Differentially expressed genes with significant genotype × development (G × D) interaction effects. Y-axis represents the RPM value. X-axis represents five different developmental stages showed higher expression at early developmental stages (cluster 1, 2, 6) were highly enriched in GO terms for "stem cell population maintenance" and "microtubule motor activity" (Table 2). DEGs with higher expression at the later developmental stages (cluster 5 and 8) were specifically enriched in "steroid biosynthetic processes".
Cluster 1, 5 and 6 were selected as they contained the GO terms that were most significantly overrepresented in any of the clusters (adj. p-value < 1e-05).
Cluster 5 was enriched with genes encoding enzymes in sterol biosynthesis. They included putative orthologs of DWARF1 Figure 7d). In contrast to genes co-expressed with SlWUS, the expression of the sterol biosynthesis pathway genes increased with the floral development stages. In addition, these genes were expressed at lower levels in fas and lc/fas compared to WT across different developmental stages, implying the negative roles of these genes in meristem maintenance.
As sterols are precursors of brassinosteroids (BRs), membrane components, and signaling molecules during plant development (Vriet, Russinova, & Reuzeau, 2013), these results have provided new insights into the potential roles of sterols/BRs in FM regulation. These genes were expressed at higher levels at 2 dpi, indicating that their expression might be positively regulated by genes mediating outer whorl initiation. Phylogenetic analysis revealed that these kinesin genes belonged to different subfamilies ( Figure S4), and the co-expression signature of these genes suggested that they might act cooperatively in cell division and cell growth during floral development. In addition to the genes encoding kinesin proteins, the genes closely related to Cyclin B2;3 (Solyc03 g032190) and Cellulose Synthase-Like D5 (Solyc09 g075550) were also found in cluster 6.
Although kinesins might also play a role in organelle movement (Zhu & Dixit, 2012), the co-expressed pattern of these kinesin genes with putative tomato Cyclin B2;3 and Cellulose Synthase-Like gene suggested that they were more likely involved in cytokinesis (Hunter et al., 2012;Tank & Thaker, 2011).

| Conserved regulatory mechanism between LC-FAS and WUS-CLV3
In the recent years, a number of genes associated with the control of tomato fruit shape and size have been identified (Mu et al., 2017;van der Knaap & Østergaard, 2017;van der Knaap et al., 2014). The key genes that contribute to the flat shape and multi-loculed tomato fruits are SlWUS and SlCLV3 Xu et al., 2015). The origin and distribution of the mutant alleles of these genes, including their effect on tomato fruit morphology was investigated at the population level (Blanca et al., 2015;Rodríguez et al., 2011). While lc appears to have arisen in wild relatives, fas arose later during domestication (Blanca et al., 2015). Although fas appears to be more effec- Unlike Arabidopsis, in which CLV3 is expressed in the L1, L2 and L3 layer (Fletcher et al., 1999), we found that SlCLV3 was absent from L1 and its expression was above and partially overlapped with the SlWUS expression domain in tomato. Similar results were also observed in soybean, in which GmCLV3 was absent in L1-L3 layer and its expression domain overlapped with GmWUS below the L5 layer in the SAM (Wong, Singh, & Bhalla, 2013). Together, these findings raise the possibility that the CLV-WUS meristem regulatory mechanism has somewhat diverged across different species.
In summary, our results imply that lc and fas mutations cause an expansion of the SlWUS expression domain and delay the termination of SlWUS, resulting in the production of larger FMs and more floral primordia.

| Genes responding to LC and FAS expression dynamics during early floral development
Compared to a previous report using a  Figure S6). Together, these results indicated that SlCLV3, in addition to its involvement in maintaining the meristem cell population with SlWUS, might be involved in other developmental processes as well.
We also observed a trend of decrease in the number of overlapping genes as floral buds developed further ( Figure S6). Because whole flower buds were used in this study ( Figure S1b), the growing mass of sepals, petals and stamens might dilute the meristemspecific transcripts and therefore hinder the detection of DEGs between genotypes.

| Novel mechanisms underlying LC-FAS mediated control of meristem development revealed by co-expressed gene clusters
The transcriptional network controlling tomato meristem development is not fully unexplored. To reveal potential players participating in LC-and FAS-mediated meristem and floral development programs, a time-course RNA-seq gene expression profiling was conducted. From the cluster analysis, we uncovered a cluster enriched with genes related to microtubule motor activity (Figure 7, Data S1). This group of genes only showed higher expression in lc compared to the wild type and fas, especially in the SYM ( Figure S7, Data S1). It is possible that this group of genes is activated by CLV-mediated MAPK signaling cascade (Betsuyaku, Takahashi et al., 2011), as CLV activity is elevated in lc. In Arabidopsis, MAPK cascade targets various TFs involved in a plethora of developmental, defense, and stress responses (Popescu et al., 2009). The tobacco MAPK cascade, positively regulates cytokinesis by phosphorylating NtMAP65-1, a microtubuleassociated protein (Sasabe et al., 2006). In our results, we also found an elevated expression of a tomato microtubule-associated protein gene MAP65-1a (Solyc11 g072280) in lc.

The identification of a cluster of genes involved in phytosterols
and BRs synthesis indicates their putative roles in meristem function during tomato floral development. Genes in this cluster were expressed at lower levels in fas and lc/fas, indicating that the BR level might also be lower in the meristem of fas and lc/fas (Figure 7, Data S1). Phytosterols are the precursors of BRs and the homeostasis of BRs is controlled by a feedback regulatory loops (Vriet et al., 2013). A previous study demonstrates that phytosterols have a BRindependent role in controlling and activating signals for plant development (He et al., 2003). The sterol biosynthesis activities appear to be highly localized to the meristematic region. For example, the tobacco squalene synthase (SQS) enzyme activity is predominantly detected in SAM (Devarenne, Ghosh, & Chappell, 2002). Arabidopsis 3-hydroxy-3-methylglutaryl coenzyme A reductase 2 (HMG2) gene is expressed in SAM and floral tissues (Enjuto, Lumbreras, Marín, & Boronat, 1995). In addition, Arabidopsis FACKEL gene is set involved in embryonic patterning and meristem programming (Jang et al., 2000).A Although we propose that sterol biosynthesis in meristematic regions is affected in fas, we cannot rule out the possibility that the linked genes located within the introgressed fas inversion are the cause of the differential expression.
Brassinosteroids are growth promoting hormones in general and the BR contents are maintained at a low level in the meristem, particularly in the organ boundary (Hepworth & Pautot, 2015). The KNOX genes, such as STM, maintain the identity of meristem and boundary in the SAM by suppressing the BR levels and directly activating genes involved in boundary formation (Bolduc et al., 2012;Johnston et al., 2014;Spinelli, Martin, Viola, Gonzalez, & Palatnik, 2011;Tsuda, Kurata, Ohyanagi, & Hake, 2014). In Arabidopsis, BR biosynthesis mutant, det2, and BR membrane-bound receptor mutant, bri1, cause extra carpel formation (Gendron et al., 2012). to LA1589 and identified close recombination breakpoints, followed by two more self-pollinated generations to generate final BC 9 F 2 population (family 13S133) ( Figure S8). During this selection, lc and fas loci were maintained in the heterozygous state, while the surrounding loci were selected to be homozygous wild type and selected for recombinants around the genes. Three NILs, lc, fas, lc/fas and the wild type (WT), were created from the BC 9 F 2 population. The primers used to select recombinants are listed in Table S3.
The size of the introgressed segment varied for each locus

| Cloning of the RNAi-SlCLV3 constructs
To reduce the expression of SlCLV3 in wild type tomatoes, a hairpin RNAi construct was created using the pKYLX80 vector similar to the method described in Siminszky, Gavilano, Bowen, and Dewey (2005). A gene-specific fragment of 355 bp was amplified from SlCLV3 with the following primer pairs: CRF1: 5′-AATTCTAGAAGCTTTCAATCTCT TTGTCTTGCTGA-3′ and CRR1: 5′-ATGGAGCTCTCGAGATGAA ACCATATACTACCCT-3′. The amplified product was digested with HindIII/XhoI and SacI/XbaI to construct sense and antisense fragments flanking the 151 bp region of soybean ω-3 fatty acid desaturase (FAD3) intron ( Figure S10).
Next, both digested fragments were inserted into vector pKYLX80.
The resulting EcoRI-XbaI fragment from pKYLX80 containing the CaMV35S2 promoter, SlCLV3 sense hairpin-stem, FAD3 intron and SlCLV3 antisense hairpin-stem was subcloned into binary vector pKYLX71 between TL border and the RBCS subunit terminator to produce the final RNAi-SlCLV3 construct, named pRNAi-CR. The pRNAi-CR was stably transformed into S. pimpinellifolium accession LA1589. We selected two independent T0, pRNAi-CR4 and pRNAi-CR9, which were shown to contain six and one copy of the transgene respectively, and were further evaluated by phenotypic and expression analysis in the T1 generation. Target specificity of this RNAi experiment was examined through BLASTN using Tomato genome cDNA database (SL3.20) with the designed hairpin-stem sequence ( Figure S11). Comparisons of the expression level of all the SlCLE gene families between the wild type and RNAi-SlCLV3 lines in tomato FM and 2dpi floral buds are shown in Figure S12.

| Fruit weight and dimension analysis
For fruit weight analysis, 20 ripe representative fruits per plant were selected and weighted. For fruit dimension analysis, eight to ten fully mature fruits from each of the genotypes were cut horizontally and scanned. Tomato analyzer 3.0 (Rodríguez et al., 2010) was used to analyze the scanned images for fruit perimeter and area following the instructions (http://vande rknaa plab.uga.edu/tomato_analy zer. html).

| Morphological analysis of inflorescence structure and meristem size measurement
The first young inflorescences of lc, fas, lc/fas NILs and the wild type were collected in the greenhouse and immediately placed in ice-cold RNAlater (QIAGEN) to preserve the tissue structure. The inflorescences were imaged using an Olympus SZH10 stereomicroscope and Olympus DP-10 digital camera. For the meristem size measurement, paraffin slide sections were made the same way as described in in situ hybridization procedures. Paraffin slide sections were rehydrated through ethanol series and stained with 1% Toluidine Blue (Sigma-Aldrich). Stained tissues were further dehydrated through ethanol series and finally mounted with Cytoseal 60 (Thermo Scientific). Images were taken under a fluorescence microscope and meristem size was measured using ImageJ software (NIH). The width of floral meristems was measured along a line between two sepal and petal primordia in floral buds at stage 3 and 4 days post floral initiation (dpi), respectively. For each genotype and time point, at least five meristems were measured. A two-tailed ttest was performed for statistical analysis.

| Statistical analysis
Analysis of variance (ANOVA) and Tukey's mean separation tests (HSD) were performed using the average of 20-40 measurements from each plant except for fruit size dimension, in which an average of eight fruits were used. Comparisons were made using the average per plant and 3-7 plants per genotype. Epistasis between the two loci was determined using two-way ANOVA with the following model:

| In situ hybridization to determine LC and FAS expression in floral meristems
RNA in situ hybridization was performed with digoxignin-labeled RNA using by Wu, Xiao, Cabrera, Meulia, and van der Knaap (2011) with minor modifications. To generate the RNA probes, full-length SlCLV3 (Solyc11 g071380) and SlWUS (Solyc02 g083950) cDNA was amplified from M82 cDNA using Phusion Taq (Invitrogen) and ligated into the pSC-A-amp/kan vector containing T7 and T3 promoter (provided by Lippman's lab, CSHL). Clones pSL-CLV3-1 and pSL-CLV3-4 were generated to make SlCLV3 antisense and sense probes, respectively. Clone pSL-WUS-4 was used to make both SlWUS sense and antisense probes. Depending on the orientation of the insert, T7 or T3 RNA polymerase was used to transcribe sense or antisense RNA.
Young inflorescences were fixed with ice-cold 4% (w/v) paraformaldehyde and vacuum infiltrated with a pressure of 25-28 in Hg for 20-30 min until the samples had sunken. Samples were dehydrated through ethanol series, followed by histoclear replacement.
Paraffin wax (Polyscience) was used for sample embedding with at least six fresh exchanges for 3 days. Microtome sections were taken to obtain 10 μm thick ribbons. Slides were rehydrated in an ethanol series followed by Proteinase K digestion and acetylation.
The hybridization reactions were conducted at 55°C overnight with the gene-specific DIG-labeled probes. Excess probes were washed off with saline-sodium citrate buffer (SSC buffer) and slides were blocked with blocking reagent (Roche). To detect the signal, slides were incubated with alkaline phosphatase-conjugated antibody (anti-DIG-AP Fab fragments, Roche) at room temperature for 2 hr.
Non-specific binding of antibody was washed three times with BSA buffer for 1.5 hr each. Finally, Western Blue (Promega) was applied to each slide and incubated overnight in dark at room temperature for the color reaction. Images were taken under the fluorescence microscope (Leica) equipped with a digital camera in Molecular and Cellular Image Center in OARDC.

| Gene expression analysis
Raw reads were checked for their quality in FASTQC. Total raw reads for each sample ranged between 5 to 10 million (Table S4).
Filtering steps were performed using fqtrim program (http://ccb. jhu.edu/softw are/fqtri m/index.shtm) to remove up to 50% low quality reads with QC score < 20 as well as poly-A tail contaminations. About 94% of clean reads mapped to ITAG2.4 Released Tomato Genome through Rsubread 3.4 (Liao, Smyth, & Shi, 2013) and 81% of clean reads mapped to annotated genes. Most of the genes were either not expressed or expressed at low levels ( Figure S13). To obtain an overview of enriched Gene Ontology (GO) terms for the 675 DEGs, Arabidopsis homologs of these DEGs were used as inputs in the Cytoscape plug-in, GlueGO v2.1.6 (Bindea et al., 2009). To identify co-expressed genes, samples were clustered based on the normalized expression (Z-score) of the wild type across four developmental stages (F&IM, 2, 4, and 6 dpi). By using the Mfuzz package (Kumar & Futschik, 2007) with fuzzy c-mean algorithm in R, DEGs were grouped into eight clusters. Furthermore, to identify a core of genes showing similar expression dynamics in mutants within each cluster, expression values from lc, fas and lc/fas NILs were Z-score normalized with WT.
The normalized expression values were used to calculate PCCs between gene pairs within each cluster. To select the core co-expressed genes, the PCC matrix was hierarchically clustered using Ward's method and visualized through heatmap3 in R (Zhao, Guo, Sheng, & Shyr, 2014). DEGs clustered within the hierarchical sub-group were objectively selected as a sub-cluster as they were tightly co-expressed among WT and mutants through different developmental stages.

| Whole mRNA-seq sample preparation, library construction, and data analysis
For the whole mRNA-seq, tissues were collected from the wild type

| Phylogenetic analysis
Phylogenetic tree analysis was performed to assigned nine differentially expressed kinesin genes to one of 10 plant kinesin families.
To retrieve the kinesin protein sequences in tomato, full-length protein sequences were downloaded from the International Tomato Maximization for Motif Elicitation (MEME) tool (Bailey & Elkan, 1994) was used to define the conserved motifs in tomato kinesin genes. The multiple sequence comparison by log-expectation (MUSCHEL) algorithm implemented in the Molecular Evolutionary Genetics Analysis program (MEGA, version7.0) was used to perform multiple sequence alignment with default settings (Kumar, Stecher, & Tamura, 2016).
Phylogeny tree was constructed using the neighbor-joining method with nucleotide p-distance and 1,000 bootstrap replicates.

ACK N OWLED G M ENTS
This research was funded by NSF-IOS 0922661. We thank Dr. Tea Meulia (The Ohio State University) for microscopy assistance, Dr.
Thomas Juenger (University of Texas at Austin) for help with the 3′ Tag RNA-Seq and the OSU-Wooster and OSU-Columbus greenhouse facilities for greenhouse facilities and plant care.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest associated with the work described in this manuscript.

AUTH O R CO NTR I B UTI O N S
YHC and ZH performed the research; YHC contributed to new analyses tools and analyzed the data; JCJ and EvdK supervised the research and writing; EvdK designed the research; YHC wrote the manuscript with revisions from JCJ and EvdK. All authors approved the manuscript.