NILs in Solanum pimpinellifolium accession LA1589 background, of which the interval regions harbor one, two or three of the tomato fruit shape loci (sun, ovate and fs8.1), were investigated in this study. Even though the mutant alleles of the fruit shape genes were not caused by the loss-of-function mutations, the three loci were still represented by lowercase italics. Details about the interval regions and the morphological characteristics of the NILs were described in Wu et al., 2015. All eight NILs, including WT, sun, ovate, fs8.1, sun+ovate (s/o), sun+fs8.1 (s/f), ovate+fs8.1 (o/f), sun+ovate+fs8.1 (s/o/f), were grown in the greenhouse with standard fertilizer application and water supply at China Agricultural University in year 2021 and 2022. In each repeat, at least eight plants of each genotype were grown and all the investigated plants were placed randomly in the greenhouse.

Analysis of variance (ANOVA) and Turkey’s mean separation tests were performed in this paper using SPSS software. Significant effects and interactions were shown by the p-values computed from the F ratio in ANOVA. Visualization of the results were conducted using GraphPad.

Sixty to eighty developing seeds were dissected from at least three fruits of each NIL at 8 dpa and immediately frozen in liquid nitrogen. Total RNA extract, library construction and sequencing were performed by Biomarker Technologies Corporation (Beijing, China) on a GAII platform (Illumina). Preliminary data analysis, including short sequences (reads) filtering, alignment, and FPKM calculation, was performed in BMKCloud (www.biocloud.net). Since very few genes were statistically significant when using DESeq2 method, genes that met the following conditions were considered as the differentially expressed genes (DEGs) and selected for the Venn diagram analysis. (i) the average value of a gene’s FPKM must be ≥2 in any of the eight NILs; (ii) the ratio of FPKM value of a gene (any of the seven mutant NILs/WT) must be ≥1.50 or ≤0.67. Venn diagram analyses were conducted using Venny 2.0.2 software (https://bioinfogp.cnb.csic.es/tools/venny/index2.0.2.html). Heatmaps were drawn via using the TBtools software.

For verifying the RNA-seq results, quantitative real-time PCR (qRT-PCR) was performed to clarify the expression patterns of some of the interested genes, including DWF1 (Solyc02g069490), 2ODD (Solyc02g062490), GRP13 (Solyc06g007920), GA2ox2 (Solyc10g007560), NCED3 (Solyc07g056570), VI2 (Solyc08g079080), SWEET14 (Solyc03g097560) and WLIM1 (Solyc03g114000). Total RNA used for RNA-seq analysis was reverse-transcribed into cDNA using a PrimeScript RT reagent kit (Takara, Dalian, Liaoning, China). Primers were designed using Primer5 software and listed in Supplementary Table S2 . SAND (Solyc03g115810) was selected as the internal control gene according to Expositó-Rodríguez et al. (2008). 2^-ΔΔCt method was used to calculate relative expression values. For each gene, three technical repeats were performed.

Phytohormones were detected by MetWare (http://www.metware.cn/) based on the AB Sciex QTRAP 6500 LC-MS/MS platform. The hormones analyzed were: indole-3-acetic acid (IAA), abscisic acid (ABA), cytokinin (CK), gibberellins (GAs), brassinosteroid (BR) and ethylene (ETH). The significance of the effect of a locus on the accumulation of a certain hormone and interactions between loci were calculated by ANOVA. Paired Student’s t-tests were performed to compare the group means when the effect was significant. The effects of each locus were tested in the WT and mutant backgrounds.

Mature seeds collected from tomato fruit shape NILs were investigated visually ( Figure 2A ). Compared with that in WT, seed size was significantly decreased in the double and triple NILs, meanwhile, decreased trends of seed size were recorded in NILs sun and ovate; in contrast, an increased trend of seed size was discovered in fs8.1 NIL ( Figure 2B ). Consistent with that, similar variation patterns were discovered for thousand seed weight ( Figure 2C ). As to seed length, it was significantly down-regulated in all the mutant NILs, except for fs8.1, in which significant decreases were only recorded in two of the three repeats ( Figure 2D ). In regard to seed width, it was significantly up-regulated in NILs ovate and fs8.1 but was down-regulated in s/f NIL; meanwhile, decreased trends of seed width were observed in NILs s/o and s/o/f; different from what is mentioned above, no significant difference of seed width was found among NILs WT, sun and o/f ( Figure 2E ). In the case of seed shape index (SSI) which is determined by both the seed length and width, it was significantly reduced in all the mutant NILs, except for the s/f NIL, in which significant decreases were only recorded in two of the three repeats ( Figure 2F ).

In order to evaluate the effect of sun, ovate and fs8.1 on the above-mentioned parameters as well as clarify the interactions among the three loci, contrasts and three-way ANOVA were performed. Based on the contrast’s results, ovate showed the greatest effect on seed length, meanwhile, all three fruit shape loci showed significant effects on the seed shape index ( Table 1 ). Epistasis was discovered in NILs s/o, s/f and s/o/f for at least one studied parameter. Particularly, interaction between sun and ovate significantly impacted seed size, thousand weight and length; meanwhile, interaction between sun and fs8.1 significantly affected seed size, thousand weight, width and shape index; moreover, interaction among the three loci influenced all the studied parameters ( Table 2 ).

In order to illustrate the histological bases for the seed size and shape variations among the tomato fruit shape NILs, paraffin sections of mature seeds were made and seed structure, seed cell number and cell size were analyzed ( Figure 1 ; Table 3 ). In terms of seed structure, seed coat area, cavity area, embryo area, embryo length and embryo crimp ratio (ratio of embryo diameter (d) to length) were investigated ( Figure 1A ). Except for the embryo crimp ratio, all the above-mentioned parameters were significantly correlated with seed size ( Supplementary Table S1 ), suggesting seed size variations were caused by the changes of seed coat area, cavity area as well as embryo area and length. Particularly, consistent with what observed in seed size, compared with those in WT, seed coat area, cavity area, embryo area and embryo length were not significantly changed in sun and ovate NILs ( Table 3 ). Whereas, in fs8.1 NIL, of which the seed size showed clearly increased trend, cavity area was significantly increased, indicating the increase of seed size in this NIL was mainly attributed to the size change of cavity area ( Table 3 ). In contrast to the single NILs, seed coat area, cavity area as well as embryo area and length were all significantly decreased in s/o and s/f NILs ( Table 3 ), suggesting those tissues/structures determined seed size of those two NILs; meanwhile, only seed coat area was significantly reduced in o/f NIL ( Table 3 ), which may be the major reason for the smaller seed size in this NIL; in the s/o/f NIL, seed coat area and embryo area were significantly downregulated ( Table 3 ), implying the seed size variation in this NIL was mainly due to the changes of those two tissues.

Organ size variation is generally resulted from the alterations of cell division and/or expansion patterns. In this study, compared with that in WT, seed coat cell number in proximal-distal direction (P-D) was significantly decreased in sun, double and triple NILs; whereas, seed coat cell layer was not significantly changed in any of the mutant NILs; as to the seed coat cell size, it was only significantly increased in the distal and body (DB) areas in s/f NILs, but was not changed in the proximal area ( Table 3 ). In regard to the embryo, cotyledon cell number was significantly up- and down-regulated in ovate and s/f NILs, respectively; meanwhile, hypocotyl cell number was significantly increased in ovate, fs8.1, o/f and s/o/f NILs, but was decreased in s/f NILs; as to the embryo cell size, no significant change was observed in cotyledon and hypocotyl ( Table 3 ).

In the aspect of seed shape, it was not only determined by cell number and size, but also controlled by the embryo crimp pattern. However, since cell size of embryo was not affected by tomato fruit shape loci ( Table 3 ) and the significant change of seed coat cell size in the DB of s/f NIL was opposite to variation trends of seed length and width ( Table 3 ; Figures 2D, E ), it can be presumed that cell number may played a more important role in the regulation of seed shape than cell size. Taking together, in this study, tomato fruit shape loci regulated seed size and shape mainly through impacting the cell number of seed coat and embryo.

For clarifying the effects of fruit shape loci on seed development, developing seeds were separated from young fruits at 8, 10, 14 and 20 dpa ( Figure 3A ). Similar development trends were observed in all the NILs for seed size, length and width, which looked like equilibrium curves ( Figures 3B–D ; Supplementary Figure S1 ). In contrast, N-shaped development curves were recorded in all eight NILs for SSI ( Figure 3E ; Supplementary Figure S1 ). At each stage, the effects of fruit shape loci on seed size and shape were similar to those at the mature stage, although statistical significance was not always found between each mutant NIL and WT NIL in every repeat ( Figures 3A–E ; Supplementary Figure S2 ). In another aspect, at 8 dpa, s/o, s/f and s/o/f significantly decreased seed size in all the repeats, and meanwhile, o/f significantly down-regulated SSI in all the repeats ( Figures 3B–E ; Supplementary Figure S2 ), suggesting fruit shape loci and their interactions played critical roles in the regulation of seed size and shape at no later than 8 dpa.

Compared with WT, IAA was significantly decreased in NILs s/f, o/f and s/o/f; total GA content was significantly up-regulated in fs8.1 NIL, but was down-regulated in NILs ovate, s/o, o/f and s/o/f; except for NILs sun and ovate, levels of trans-zeatin and trans-zeatin riboside were reduced in all mutant NILs; BR content was significantly decreased in sun NIL, but was increased in NILs fs8.1 and s/o/f; the accumulation of ABA was only enhanced in fs8.1 NIL ( Figure 4 ).

For exploring the regulatory mechanisms of fruit shape loci on seed size and shape at the transcriptional level, RNA-seq analysis was performed with seeds collected from 8-dpa fruits. Comparisons were conducted between each mutant NIL and WT NIL, and DEGs were identified. sun, ovate, fs8.1, s/o, s/f, o/f and s/o/f led to 255 (up-/down-regulated (hereafter u/d) =124/131), 174 (u/d=69/105), 559 (u/d=274/285), 311 (u/d=133/178), 707 (u/d=329/378), 603 (u/d=298/305) and 1051 (u/d=617/434) DEGs, respectively ( Figure 5A ). Enrichment and pathway analyses indicated that there were common pathways or function classes among sun, ovate and fs8.1, and there were additive and epistatic effects among the three loci ( Supplementary Figures S3-S6 ). Transcription factor enrichment revealed that B3, AP2 and MYB were the most common families that were discovered in the developing seeds ( Supplementary Figure S7 ).

For further dissecting the transcriptional bases of each fruit shape locus and the interactions between different loci, Venn diagram analyses were performed ( Figures 5 , 6 ). In the Venn diagram of NILs sun, ovate and fs8.1, there were 149, 83 and 468 DEGs unique to NILs sun, ovate and fs8.1, respectively; meanwhile, 37, 37 and 22 DEGs were shared between NILs sun and ovate, NILs sun and fs8.1, and NILs ovate and fs8.1, respectively; additionally, 32 DEGs were shared by NILs sun, ovate and fs8.1 ( Figure 5B ). Based on genes’ annotations, 43 interested DEGs were screened out ( Figure 5B ). Among them, SUN was highly expressed in the NILs harboring sun; CWI1 (encoding a cell wall invertase which is involved in sucrose degradation), GRP13 (encoding a gibberellin-regulated protein), GASA4 (encoding a gibberellin-regulated protein), DCR (encoding a BAHD acyltransferase whose homolog controls rice grain size by modulating brassinosteroid homeostasis (Feng et al., 2016)), BEH4 (encoding BES1/BZR1 homolog protein 4 which takes part in BR signal transduction), CML19 (encoding a calmodulin-like protein which mediates multiple reactions involved in regulation of plant growth), NCED3 (encoding a 9-cis-epoxycarotenoid dioxygenase which plays a part in ABA biosynthesis) and MAP70-5 (encoding a microtubule-associated protein 70-5 which may be involved in the regulation of cell anisotropic growth, division patterns, and overall organ shape) were highly expressed in fs8.1 NIL; meanwhile, IQD10/21, CYP724A1 (encoding a cytochrome P450 enzyme which plays a role in brassinosteroid biosynthesis), IAA19/29 (encoding AUX/IAAs which participates in auxin signal transduction), PIN3 (encoding an auxin polar transporter), NCED1, CYP704A2 (encoding a cytochrome P450 enzyme whose homolog is associated with seed size in rice (Tang et al., 2016)), CEPR2 (encoding a receptor protein-tyrosine kinase which mediates nitrate uptake signaling) and CKX7 (encoding a cytokinin oxidase which is involved in cytokinin degradation) were up-regulated in NILs harboring fs8.1; in contrast, expression of NR2 (encoding a nitrate reductase which takes part in nitrogen assimilation), NRT1.1 (encoding a dual-affinity nitrate transporter), IAA36, CKI7 (encoding a cyclin-dependent kinase inhibitor 7 which is involved in cell cycle control), GIF3 (encoding a GRF1-interacting factor whose homolog controls seed size in Arabidopsis), GA2ox2 (encoding a gibberellin 2-oxidase which plays a role in gibberellin degradation), WLIM1 (encoding a F-actin binding protein containing LIM domain which takes part in actin cytoskeleton remodeling (Thomas et al., 2007)) was reduced in the NILs harboring fs8.1; as to CYCD3 (encoding a cyclin-D3-3 which is involved in the switch from cell proliferation to the final stages of differentiation), CYP85A1 (encoding a BR6ox enzyme which is involved in brassinosteroid biosynthesis), IAA1, CYP707A1 (encoding an abscisic acid 8’-hydroxylase which catalyzes the degradation of ABA), they highly expressed in NILs s/f, o/f and s/o/f; expression of SWEET14 (encoding a bidirectional sugar transporter protein), 2ODD (encoding a 2-oxoglutarate-dependent dioxygenases which is involved in gibberellin biosynthesis), DWF1 (encoding a steroid synthetase which plays a role in brassinosteroid biosynthesis), NRT1.11, GA2ox2, OFP10, WLIM1, 2ODD was obviously higher in WT NIL ( Figure 5B ). In regard to left genes, PIN5 was up-regulated in NIL fs8.1 and s/o; GRP6 was activated in NILs fs8.1 and s/f; VI2 (encoding a vacuolar invertase 2 protein which is involved in sugar metabolism), STP1 (encoding a sugar transporter protein) and NRT1.1 clearly highly expressed in o/f NIL ( Figure 5B )

In the Venn diagram of NILs sun, ovate and s/o, 135, 75 and 192 DEGs were unique to NILs sun, ovate and s/o, respectively; meanwhile, 31, 51 and 30 DEGs were only shared between NILs sun and ovate, NILs sun and s/o, and NILs ovate and s/o, respectively; in addition, 38 DEGs were shared by the three NILs ( Figure 6A ). According to their annotations, 14 interested DEGs were screened out and 10 were unique to s/o NILs. Among the 10 DEGs, CEPD2 (encoding a member of the CC-type glutaredoxin family which is induced by nitrogen starvation) and NRT1.1/2/2.4 were obviously up-regulated in the s/o NIL; meanwhile, an abscisic acid receptor gene PYR1-L was clearly down-regulated in s/o NIL but expressed higher in NILs WT, sun, ovate and s/o/f; as to CKX3, it highly expressed in NILs s/f, o/f and s/o/f; expression of SCL3, a gene involved in GA signaling, was higher in NILs WT and sun; as to GA20ox1, it was obviously up-regulated in fs8.1 NIL, but was clearly down-regulated in s/o NIL; in regard to the expression of IAA1 and 2, the highest levels were recorded in NILs s/o/f and s/o, respectively ( Figure 6A ).

In the Venn diagram of NILs sun, fs8.1 and s/f, 119, 289 and 399 DEGs were unique to NILs sun, fs8.1 and s/f, respectively; meanwhile, 29, 67 and 201 DEGs were only shared between NILs sun and fs8.1, NILs sun and s/f, and NILs fs8.1 and s/f, respectively; in addition, 40 DEGs were shared by the three NILs ( Figure 6B ). Based on the putative functions, 45 interested DEGs were screened out, of which 26 were unique to s/f NIL ( Figure 6B ). Among the 26 DEGs, GH3.1 (encoding a Grethchen Hagen 3 family protein which takes part in inactivation of IAA), GH3-10, GA2ox2, CKX1/3/7, IPT1 (encoding an isopentenyl transferase which is involved in trans-zeatin biosynthesis), SWEET1, SUS4 (encoding a sucrose synthase protein), CYCD3.1/6.1, CYP78A10, IQD6 and OFP13 were up-regulated in NILs s/f, o/f and s/o/f; in contrast, CID1L2 and CAR4 were down-regulated in the above-mentioned NILs; expression of SAUR50, encoding a SAUR-like auxin-responsive protein, was increased in NILs harboring fs8.1, while, expression of WLIM1-L was decreased in the same NILs; expression of GH3.2, GAST1-L and IQD19 was enhanced in sun NIL, but was repressed in s/f NIL; as to CYCB2.4 and ATN (encoding a microtubule-binding protein), their expression levels peaked in s/f NIL, and meanwhile, the expression of SAUR71 peaked in fs8.1 NIL; in regard to 2ODD, its expression was up-regulated in NILs WT and s/o ( Figure 6B ).

In the Venn diagram of NILs ovate, fs8.1 and o/f, 88, 320 and 358 DEGs were unique to NILs sun, fs8.1 and s/f, respectively; meanwhile, 26, 32 and 185 DEGs were only shared between NILs ovate and fs8.1, NILs ovate and o/f, and NILs fs8.1 and o/f, respectively; in addition, 28 DEGs were shared by the three NILs ( Figure 6C ). Thirty-eight interested DEGs were screened out of the above-mentioned DEGs, among which GH3-10, IAA3, GA2ox1, CKX1/3, SWEET1-L, SUS4, CYCD6.1, CYP78A10 and CML25 were up-regulated in NILs s/f, o/f and s/o/f, however, CAR4 was down-regulated in the same NIL; expression IAA2 was increased in NILs s/o, o/f and s/o/f; SAUR50, CKX7 and OFP13 were highly expressed in NILs harboring fs8.1, meanwhile, expression of SAUR70 peaked in fs8.1 NIL; expression of GIF1 and TUB3 (encoding tubulin beta chain 3) was increased in NILs sun, fs8.1, s/f, o/f and s/o/f; as to the CML3/8/15/22/24, they all highly expressed in NILs o/f and s/o/f ( Figure 6C ).

In order to clarify the epistatic interaction among sun, ovate and fs8.1 in the triple NIL, Venn diagram analysis was also performed with NILs s/o, s/f, o/f and s/o/f. The results showed that 400 DEGs were unique to s/o/f NIL; 41, 134 and 71 DEGs were only shared between NILs s/o/f and s/o, NILs s/o/f and s/f, and NILs s/o/f and o/f, respectively; 27, 308 and 15 DEGs were only shared by NILs s/o/f, s/o and s/f, NILs s/o/f, s/f and o/f, and NILs s/o/f, s/o and o/f, respectively; 55 DEGs were shared by all four NILs ( Figure 5C ). Sixty interested DEGs were screened out from the above-mentioned DEGs, of which 21 were unique to s/o/f NIL ( Figure 5C ). Among the 21 DEGs, the highest expression levels of YUCCA6 (encoding an indole-3-pyruvate monooxygenase), GH3.1/3.6, GA2ox8-L, DWF2, CYP85A1, CYCA3.2, GRF3-L, WRKY31, dUTPase (encoding a deoxyuridine 5’-triphosphate nucleotidohydrolase whose homologous positively regulates grain length in rice (Utsunomiya et al., 2011)), NRT1.11/1.2, MAP65-6 and CML35/60D were all in the s/o/f NIL; in contrast, for the left genes, including GA2ox1, GRP, GRP13, SPSB (encoding a sucrose-phosphate synthase B) and OFP10, the lowest expression levels were all discovered in the s/o/f NIL ( Figure 5C ).

In another aspect, when considering the expression patterns of all the interested DEGs, they could be classified into 17 groups. Group I to XVII mainly represented the DEGs that highly expressed in NILs harboring fs8.1; only in fs8.1 NIL; in NILs s/f and s/o/f; in NILs s/f, o/f and s/o/f; only in s/o/f NIL; in NILs WT, ovate and fs8.1; in NILs WT and sun; in NILs WT and fs8.1; in NILs WT, ovate, fs8.1 and s/o; in NILs WT, sun and fs8.1; in NILs WT and ovate; in NILs WT, ovate and s/o; in NILs WT, sun, ovate and s/o; in NILs WT, sun, ovate, fs8.1 and s/o; only in sun; only in s/o; and DEGs that lowly expressed in fs8.1 NIL, respectively ( Figure 6D ).

Besides what is mentioned above, for further understanding the mechanism of fruit shape loci in the regulation of seed morphogenesis, expression of SUN and OVATE was screened out from the RNA-seq data of S. pimpinellifolium ovaries (Pattison and Csukasi, http://ted.bti.cornell.edu/cgi-bin/TFGD/digital/experiment.cgi?ID=D009) as well as our developing seeds. In the WT NIL, SUN barely expressed in anthesis ovules, 4-dpa embryos, 4-dpa endosperm, 4-dpa seed coat and 4-dpa funiculus but highly expressed in 8-dpa seeds ( Figures S8A, B ). While, OVATE highly expressed in anthesis ovules and 4-dpa funiculus but barely expressed in embryos, endosperm and seed coat at 4 dpa ( Figures S8C, D ).

In order to verify the expression patterns of some important genes, qRT-PCR was conducted and the results were consistent with those observed in the RNA-seq ( Supplementary Figure S9 ).

Since sun, ovate and fs8.1 have been fine-mapped in tomato, their functions in the regulation of plant organ morphogenesis have been well studied in Arabidopsis, rice and a number of horticultural crops (Wang et al., 2011; van der Knaap et al., 2014; Burstenbinder et al., 2017b; Pan et al., 2017; Dou et al., 2018; Lazzaro et al., 2018; Yang et al., 2018a). sun and some gain-of-function mutations of IQDs regulate organ shape through simultaneously promoting and inhibiting cell proliferation in the longitudinal direction and cross direction, respectively (Wu et al., 2011; Yang et al., 2018a). ovate and the loss-of-function mutations of OFPs elongate fruits and tubers mainly through increasing cell number in proximal-distal direction but decreasing cell number in medio-lateral direction, especially in the proximal part (Wu, 2015; van der Knaap and Ostergaard, 2018). As to fs8.1, it enlarges fruit shape index and size mainly through enhancing cell division along the longitudinal axis (Sun et al., 2015). However, in this study, totally different regulatory patterns of the fruit shape loci were discovered in the control of seed shape and size. In particular, sun significantly reduced SSI through decreasing seed length, which could be attributed to the inhibition of seed coat cell proliferation in proximal-distal direction ( Figures 2A, D–F ; Tables 1 , 3 ). This pattern was not only different from that in fruit, but also in contrast with that of GSE5 in rice (Duan et al., 2017). GSE5 encodes an IQD which is closely related to Arabidopsis IQD25-27. Up-regulation of GSE5 led to increased SSI through increased cell proliferation in spikelet hulls, whereas, down-regulation of its expression resulted in wider grains (Duan et al., 2017). The difference between sun and rice GSE5 mutations may be attributed to the different seed structures between the dicotyledon and monocotyledon. In regard to ovate, it down-regulates SSI through decreasing seed length but at the same increasing seed width ( Figures 2A, D–F ; Table 1 ). Since seed coat cell number was not significantly changed and the total cell number of embryo was increased, the decreased seed length of ovate NIL could be attributed to the down-regulation of embryo crimp ratio ( Figure 2D ; Table 3 ). As to the seed width, now that no significant difference was found in seed coat cell layer, it could be presumed that the increased seed width of ovate NIL was attributed to additive effect of the increased cavity size and the increased embryo cell number in medio-lateral direction ( Figure 2E ; Table 3 ). In the fs8.1 NIL, the reduced SSI was mainly due to the significantly increased seed width, which was high-likely caused by the additive effect of the increased trends of seed coat cell layer and embryo cell number in medio-lateral direction ( Figures 2D–F ; Table 3 ). Besides the single loci, double loci interaction and triple loci interaction also played critical roles in impacting seed length, width and shape index ( Table 2 ), which was also not fully consistent with what has been observed in the regulation of fruit shape (Wu et al., 2015). In the previous investigation of same NILs, only interactions of sun×ovate and sun×fs8.1 showed significant effects on fruit shape index (Wu et al., 2015), whereas, in this study, interactions of sun×fs8.1 and sun×ovate×fs8.1 significantly impacted SSI ( Table 2 ), which further suggested the different patterns and mechanisms underlying the regulation of seed shape by fruit shape loci.

In another aspect, significant decreases of seed size and weight were found in the double and triple NILs ( Figures 2B, C ) and this phenomenon was never observed in fruits of the same NILs (Wu et al., 2015), which also indicated the mechanisms of fruit shape loci in the regulation of seed development was different from those in fruit. Based on the histological analysis, s/o and s/o/f significantly decreased seed coat cell number and embryo crimp ratio, which led to the down-regulation of seed coat size and embryo size, and finally resulted in smaller seeds ( Figures 2A, B ; Table 3 ). o/f significantly decreased seed size mainly through reducing the seed coat size, which could be attributed to the significant decrease of cell number in seed coat ( Figures 2A, B ; Table 3 ). s/f had the smallest seed size which was caused by the significant reduction of seed coat size and embryo size. The decreased embryo size in s/f NIL was mainly attributed to significant down-regulation of cell number in cotyledon and hypocotyl. As to the seed coat, since cell number was significantly decreased, logically the reduction of seed coat cell number should be the major reason for the decrease of seed coat size ( Figures 2A, B ; Table 3 ).

Taking together, although sun, ovate and fs8.1 still mainly acted as cell number regulators in this study, they played obviously different roles in the regulation of tomato seed shape and size than they did in the control of fruit development.

For further understanding the mechanisms underlying the morphological and histological variations of seeds induced by the fruit shape loci, phytohormones and transcriptomes of the developing seeds were investigated, and phytohormone-, sugar- and cytoskeleton-related pathways were discovered to be significantly impacted by the fruit shape loci.

Phytohormones are always considered as critical regulators of seed development. GA, CK and BR generally play positive roles in the regulation of seed size (Werner et al., 2003; Jiang et al., 2013; Fang et al., 2016; Hu et al., 2017; Zhou et al., 2017; Gan et al., 2022), whereas, ABA has been proved to be negatively associated with seed size (Cheng et al., 2014; Wang et al., 2021). As to auxin, its local accumulation and polar transportation also take part in the control of seed development (Liu et al., 2015b). In this study, sun, ovate and fs8.1 showed distinct effects on adjusting the phytohormone-related pathways. Particularly, sun only significantly reduced the accumulation of BR by down-regulating the expression of the BR biosynthesis gene DWF1 ( Figure 4D ; Figure 5B ). While, ovate significantly decreased the GA level through repressing the expression of 2ODD ( Figure 4D ; Figure 5B ), a gene taking part in GA biosynthesis (Huang et al., 2010). Additionally, expression of PIN5, a gene encoding an auxin polar transporter, was also decreased by ovate ( Figure 5B ), which could impact the auxin polar transportation and local auxin gradients in developing seeds, and may finally change seed shape and size. As to fs8.1, it showed a pleiotropic effect on phytohormones that GA, ABA and BR contents were significantly increased, whereas, the CK level was obviously decreased ( Figures 4B–E ). Consistently, in fs8.1 NIL, expression of GA2ox and CYP707A1 was down-regulated, whereas, expression of NCED1, NCED3, CYP724A1 and CKX was up-regulated ( Figure 5B ). The impacts of the fruit shape loci on the phytohormone level were also discovered in the developing flowers at 9 dba (days before anthesis) (Wu et al., 2015). In young flower buds, sun significantly increased JA (jasmonic acid) and JA-Ile (isoleucine conjugated-JA) levels in the absence of the fs8.1; ovate had an impact on iP (isopentenyladenine, a kind of CK) accumulation; fs8.1 not only tend to reduce IAA and increase ABA in 9-dba flower buds but also inhibited TRA (tryptamine, an auxin precursor) accumulation in anthesis ovaries and 5-dpa fruits (Sun et al., 2015; Wu et al., 2015). Therefore, it can be presumed that phytohormones are high-likely the targets of the three fruit shape loci, however, the regulatory mechanisms of different loci or the same locus in different organs are not exactly the same. Interestingly, the increased ABA level was observed in both the developing seeds ( Figure 4E ) and young flower buds (Wu et al., 2015) harboring fs8.1, indicating this may be a conserved mechanism underlying the regulation of organ morphology by fs8.1.

Additionally, interactions were discovered among different loci at the hormone level. In sun background, in addition to the decreased CK level, fs8.1 also significantly reduced IAA and GA levels but did not impact ABA level ( Figures 4A–C ). Consistent with that, in the s/f NIL, GH3.1, GA2ox and GKX1/3/7 were up-regulated, whereas, 2ODD was down-regulated ( Figure 6B ). However, in ovate background, fs8.1 only reduced IAA and CK levels by up-regulating the expression of GH3-10 and CKX1/3/7 ( Figures 4A, C ; Figure 6C ). Thus, it can be deduced that interactions of sun×fs8.1 and ovate×fs8.1 enhanced the impact of fs8.1 on the IAA accumulation but eliminated the positive effect of fs8.1 on ABA and GA accumulations. In another aspect, interaction between sun and ovate significantly reduced GA and CK levels through simultaneously down-regulating the expression of GA20ox1 but up-regulating the expression of CKX3, respectively ( Figure 4B, C ; Figure 6A ), which has never been observed in NILs harboring sun or ovate alone ( Figure 4 ; Figure 5B ), indicating new regulatory pathways were established by the interaction. As to the interaction of sun×ovate×fs8.1, it inherited most of the effects of the interactions of the double NILs, except for that of ABA, which was significantly increased in s/o/f NIL ( Figures 4A–D ). Interactions among sun, ovate and fs8.1 were also observed in 9-dba young flower buds in the regulation of iP, SA, JA and JA-Ile but not in that of IAA and ABA (Wu et al., 2015), indicating the interaction patterns among the three fruit shape loci may be different in the regulation of different organs.

Sugar metabolism and transport also play an important part in the regulation of seed development (Xu et al., 2019; Liao et al., 2020; Wang et al., 2020). In the developing seeds, sun significantly down-regulated the expression of VI2 and SWEET ( Figure 5B ), which could influence the energy transportation from the source to the sink and may repress cell number by inhibiting cell proliferation. Evidences supported the above-mentioned hypothesis came from studies in rice and soybean that a mutation of OsVIN2, a homolog of VI2, resulted in decreased grain size by altering sugar metabolism (Xu et al., 2019), meanwhile, GmSWEET10a was confirmed to be positively correlated with seed size (Wang et al., 2020). Similar to what happened in sun NIL, expression of SWEET and VI2 was also down-regulated in ovate NIL ( Figure 5B ). Furthermore, expression of STP1, which belongs to the monosaccharide transporter (MST) gene family and plays a key role in maintaining source/sink characteristics (Liu et al., 2021), was also repressed ( Figure 5B ), suggesting ovate may have a stronger impact on the source-sink relationship of seed than sun. Different from what is mentioned above, fs8.1 increased the expression of CWI1 ( Figure 5B ). This gene encodes a cell wall invertase which is involved in the hydrolyzation of the sucrose into glucose and fructose and is positively corelated with seed size (Liao et al., 2020). Thus, it can be presumed that the smaller seeds in NILs sun and ovate could be partially attributed to the impaired source-sink relationship, whereas, the larger seeds in fs8.1 NIL are high-likely partially ascribed to the activated hydrolyzation of the sucrose.

Interactions were also discovered in the sugar-related pathways. Compared with that in sun, the inhibition of sugar transportation and degradation was relieved in s/f NIL, because expression of SWEET14 and VI2 was no longer repressed and the transcriptional levels of SWEET-like and SUS4 were up-regulated ( Figure 6B ). SUS4 encodes a sucrose synthase (SuSy) and is capable of catalyzing the first degradative step in sucrose utilization (Yao et al., 2020). Similar to what happed in s/f NIL, the inhibition of sugar transportation and degradation was also relieved in o/f NIL, which was caused, on the one hand, by the up-regulation of SWEET1 and SUS4, and on the other hand, by relieving the repression of the expression of SWEET14 and VI2 ( Figure 6C ). These phenomena suggested that a common interaction pattern may be shared by sun×fs8.1 and ovate×fs8.1 in the sugar transportation and degradation during seed development. As to the interaction between sun and ovate, it is worth mentioning that although sun and ovate alone inhibited the expression of sugar transportation- and degradation-related genes, those genes were not significantly altered in s/o NIL, indicating that the interaction of sun×ovate relieved this inhibition. In regard to s/o/f NIL, only SPSB, a gene encoding a sucrose-phosphate synthase, was identified in the s/o/f unique DEGs ( Figure 5C ), indicating the interaction of sun×ovate×fs8.1 almost eliminated all the above-mentioned effects of the fruit shape loci on sugar metabolism and transport.

Many experimental evidences support that the biological functions of SUN and OVATE are based on their roles in the rearrangement of the cytoskeleton (Abel et al., 2013; van der Knaap et al., 2014; Burstenbinder et al., 2017b; Lazzaro et al., 2018; Wendrich et al., 2018; Wu et al., 2018). In this study, sun reduced the expression of WLIM1 (Solyc03g114000) ( Figure 5B ), which encodes a F-actin binding protein and is involved in actin cytoskeleton remodeling (Thomas et al., 2007). Interestingly, in the developing flower buds, expression of another WLIM1 gene was also impacted by sun (Wang et al., 2019), suggesting this class of gene may be the target of sun in the control of tomato fruit and seed morphology. Similarly, ovate also down-regulated the expression of WLIM1 ( Figure 5B ), on the one hand, indicating the WLIM1-related pathway was also the target of ovate in the regulation of seed morphology, on the other hand, confirming that sun and ovate share more pathways in this process. It is worth noting that another member of the OFP gene family was also down-regulated in ovate NIL, which may further enhance the effect of ovate on the reduction of seed size and shape index. In regard to fs8.1, it significantly reduced the expression of WLIM1 (Solyc04g077780) but increased the expression of MAP70-5 ( Figure 5B ), which encodes a microtubule-associated protein and may regulate fruit shape through interacting with IQD21 (Bao et al., 2022). Therefore, it can be deduced that members of the WLIM1 gene family are high-likely the targets of the three tomato fruit loci in the regulation of fruit and seed morphogenesis.

Interactions among different loci also impacted the cytoskeleton-related pathways. The impacts of sun and fs8.1 on the cytoskeleton-related genes may be enhanced by their interaction, because besides the down-regulated WLIM1 (Solyc04g077780), expression of another WLIM1 (Solyc12g088320) was also reduced in the s/f NIL, and meanwhile, expression of ATN, a gene involved in the identification of the division plane during mitosis, was up-regulated ( Figure 6B ). Similarly, in the o/f NIL, expression of WLIM1 (Solyc04g077780) was still down-regulated, however, TUB3 was up-regulated by the interaction. In contrast, interaction between sun and ovate eliminated the effect of single loci on cytoskeleton. As to the interaction of sun×ovate×fs8.1, it eliminated the effects of single loci, but at the same, established a new regulatory pathway which consisted of MAP65-6 and ROP9.

Besides what is mentioned above, based on the expression patterns of SUN and OVATE in WT NIL ( Figure 5B ), it can be presumed that variations at the phytohormone and transcriptomic levels in NILs harboring sun could be directly attributed to effects of this locus; whereas, the physiological and transcriptomic variations in NILs carrying ovate were probably ascribed to the after effects of the locus or influence of nearby tissue, such as funiculus, since OVATE still highly expressed in this tissue at 4 dpa ( Figure S8D ).

Taking together, sun impacted tomato seed morphology mainly through reducing BR biosynthesis, inhibiting sugar transport and degradation as well as repressing WLIM1-related cytoskeleton arrangement; meanwhile, ovate regulated seed morphology mainly by inhibiting auxin polar transportation as well as GA biosynthesis, sugar transport and degradation, and WLIM1-related cytoskeleton arrangement; as to fs8.1, although a decreased CK level and an increased ABA level were observed in the developing seeds of this NIL, the increased GA and BR levels as well as the enhanced source-sink relationship were still able to result in an increased trend of seed size. Interaction between sun and fs8.1 changed seed morphology mainly through decreasing the IAA, GA and CK contents as well as further repressing the WLIM-related cytoskeleton pathway; interaction between ovate and fs8.1 impacted seed morphology mainly through down-regulation of the IAA and CK contents, and repressing WLIM1-related cytoskeleton arrangement; interaction between sun and ovate changed seed morphology mainly through reducing GA and CK levels; in regard to the interaction of sun×ovate×fs8.1, the effects of the decreased IAA, GA and CK levels may overcome that of increased BR level and finally down-regulated the seed size and shape index together with the rearrangement of cytoskeleton which was caused by the different expression of MAP65-6 and ROP9.