LIN Hao-wei,WU Zhen,ZHOU Rong,CHEN Bin,ZHONG Zhao-jiang,JIANG Fang-ling

College of Horticulture,Nanjing Agricultural University/Key Laboratory of Horticultural Plant Biology and Germplasm Innovation in East China,Ministry of Agriculture,Nanjing 210095,P.R.China

Abstract Fruit cracking occurs easily during the late period of fruit development when plants encounter an unsuitable environment,dramatically affecting fruit production and marketing.This study conducted the bulked segregant RNA-Seq (BSR) to identify the key regulatory gene of fruit cracking in tomatoes.BSR-Seq analysis illustrated that two regions associated with irregularly cracking were located on chromosomes 9 and 11,containing 127 candidate genes.Further,through differentially expression analysis and qRT-PCR in cracking-susceptible and cracking-resistant genotypes,the candidate gene SlGH9-15(Solyc09g010210) with significantly differential expression levels was screened.Bioinformatics analysis of the GH9 gene family revealed that 20 SlGH9 genes were divided into three groups.The phylogenetic analysis showed that SlGH9-15 was closely related to cell wall construction-associated genes AtGH9B1,AtGH9B6,OsGH9B1,and OsGH9B3.The cis-acting elements analysis revealed that SlGH9-15 was activated by various hormones (ethylene and ABA) and abiotic stresses.The expression pattern indicated that 13 SlGH9 genes,especially SlGH9-15,were highly expressed in the cracking-susceptible genotype.Its expression level gradually increased during fruit development and achieved maximum value at the red ripe stage.Additionally,the cracking-susceptible tomato showed higher cellulase activity and lower cellulose content than the cracking-resistant tomato,particularly at the red ripe stage.This study identified SlGH9-15 as a key gene associated with fruit cracking in tomatoes for the first time and gives new insights for understanding the molecular mechanism and complex regulatory network of fruit cracking.

Keywords: tomato,irregular cracking,BSR-Seq,GH9-15

1.Introduction

Tomato (Solanum lycopersicum),one of the most important vegetable crops worldwide with extremely high commercial and research value,is recognized as a model plant for investigating fruit growth and development.However,when encountering an unsuitable environment,tomato fruit is prone to cracking during growth and development periods,as well as post-harvest transportation and cleaning.Cracking may influence fruit appearance and quality and shorten shelf life.Also,cracking fruits are susceptible to infection by the pathogenic bacterium.These severely handicap the production and marketing of tomato fruits.Using balanced water,fertilizer,and exogenous hormones in production can alleviate cracking to some extent,but it is difficult to solve the problem essentially (Doraiset al.2004;Davarpanahet al.2018;Sanoet al.2018).Clarifying molecular mechanisms and breeding cracking-resistant cultivars might be the countermeasure to address the root cause of fruit cracking.

Tomato fruit cracking is mainly classified as ring cracking,longitudinal cracking,irregular cracking,cuticle cracking,and mixed cracking (Khadivi-Khub 2015).There has been some progress in fruit cracking.Hudson (1956) found that longitudinal and ring cracking in tomatoes was relatively independently inherited,but there may be common or correlated modifying factors.Li (2016) detected five QTLs associated with tomato cracking mapped in chromosomes 1,2,and 3 from cracking-resistant processing tomato‘14803’.Cuiet al.(2017) adopted fine mapping technology and identified theER4.1gene regulating tomato epidermal reticulation cracking.Cr3a,the dehiscence resistance gene associated with saccharide biosynthesis and hormone signaling transduction,was also detected through fine mapping technology (Zhuet al.2020).By bulked segregant analysis (BSA) sequencing and QTL localization,Chenet al.(2021) hypothesized that theMAPKKKgene on chromosome 2 regulates irregular tomato fruit cracking.

The commonly used traditional QTL analysis is an effective mechanism-analyzing method of many important agronomic traits,but it is limited by demanding mapping populations and suitable DNA markers.BSA and bulked segregant RNA-Seq (BSR) have been developed in recent years (Liuet al.2012;Wanget al.2021).BSR-Seq is a combination of BSA and RNA-Seq.Individuals with extreme characters from segregate populations were selected to construct two mixed pools,and the total RNA of the two mixed pools was extracted for transcriptome sequencing.Compared with BSA,BSR-Seq only requires mRNA sequencing and removes great amounts of repeated and redundant sequences (such as transposon sequences).It could localize target traits and give information about gene patterns in mixed pools so that the differential expression analysis can identify candidate genes within localized regions.Thus,BSR-Seq effectively improves sequencing efficiency,substantially reduces costs,and provides attractive gene localization approaches for large genomic crops.It has been broadly applied to investigate the genetic control of agronomic traits in pepper,lettuce,cabbage,kiwi,maize,wheat,and other plants (Liuet al.2012;Duet al.2017;Huang Zet al.2017;Liet al.2018;Weiet al.2019).Liuet al.(2020) identified a stagnant green geneCaSGR1associated with pepper fruit color formation at a 131-kb interval on chromosome 1 by BSR-Seq and linkage analysis.Linet al.(2020) identified a cold tolerance gene in kiwifruit by BSR-Seq.These offer a new approach to studying irregular fruit cracking in tomatoes.

However,there is little research on tomato irregular cracking.This study obtained the F2segregation from cracking-susceptible genotype ‘NT189’ (male parent) and cracking-resistant genotype ‘NT91’ (female parent).BSRSeq analyses were utilized to identify the genes associated with fruit cracking in tomatoes.Global gene expression analysis demonstrated that candidate genes in mapping regions enriched in cell wall construction,cuticle biosynthesis,certain biochemical processes,and signal transduction pathways.Combined with differentially expressed analysis and the quantitative real-time PCR (qRT-PCR),GH9-15was identified as the key gene that regulates fruit cracking.The whole genome identification and gene family analysis of 20GH9in tomato were also carried out,including chromosomal localization,physicochemical property,phylogenetic tree,genetic structure,conserved domain,cis-acting element,and expression pattern analyses.In addition,the cellulase activity and cellulose content of parent genotypes at three fruit developmental stages were determined.This research can improve the understanding of the molecular mechanism of tomato fruit cracking.It can also offer breeders genes and plant resources to develop desirable cracking-resistant tomato cultivars and lays a theoretical foundation for further research on tomato cracking.

2.Materials and methods

2.1.Plant materials and growth condition

This study investigated the cracking rate from more than 500 tomato germplasm resources and selected the irregularly cracking-susceptible genotype ‘NT189’ (LA2661,from Tomato Genetics Resource Center,USA) and crackingresistant genotype ‘NT91’ (Tomato Cherry Chocolate,from University of California,Davis,USA).The cracking rates of ‘NT189’ and ‘NT91’ were 71.03 and 1.68% at the red ripe stage,respectively.Both genotypes have been selfpollinated and selected for more than six generations with stable fruit cracking traits.

In the 2019 fall,the cracking-susceptible genotype ‘NT189’as the male parent was crossed with the cracking-resistant genotype ‘NT91’ as the female parent.The F2generation was generated from a self-pollinated F1individual,and plants were cultivated in the plastic greenhouse of Jiangsu Agricultural Expo Garden,China on July 26,2020.This region has distinct characteristics of subtropical monsoon climate with four clear seasons and an average annual temperature of 15.7°C.F2segregation progeny used for BSR sequencing consists of 350 individual plants.The plant spacing was 30 cm×50 cm.It adopted the single stem pruning cultivation method and topped when the plant had five fruit clusters.We employed drip irrigation and kept the irrigation constant.Pest,disease control,and fertilization practices followed the conventional management method.

For qRT-PCR analysis,fruits of cracking-susceptible(‘NT189’,‘NT48’,and ‘NT177’) and cracking-resistant(‘NT91’,‘NT54’,and ‘NT191’) genotypes and the parents at 30 h of saturated irrigation were used.The cultivation method was consistent with the F2population.

2.2.ldentification and statistics of fruit cracking traits

Fruit splitting phenotypes were investigated when fruits of the third cluster reached full maturity.A ruler was used to measure the crack length;fruits with cracks longer than 1 cm were considered cracked fruits.The classifications of mature green,turning,and red ripe stages were referred to Laiet al.(2007).The cracking rate was calculated and tested using the one-sample K-S (Kolmogorov-Smirnov)testviaSPSS Statistics 25.

2.3.BSR-Seq analysis

RNA extraction,mixed pool constructionFor BSR-Seq analysis,four cDNA libraries were built with two parent pools (CR and CS) and two F2generation mixed pools (22 extremely cracking-resistant and 22 extremely crackingsusceptible F2individuals were evenly mixed as CRMP and CSMP).Genomic RNA was extracted from red ripe fruit in the parents and F2mix pools.Library construction and RNA sequencing were carried out at Genepioneer Biotechnology Company (Nanjing,China).

Analysis of BSR-Seq dataClean reads were attained by filtering the raw reads with low quality.Filtered reads were aligned to the tomato reference genome using HISAT2(v2.0.5).The single nucleotide polymorphisms (SNPs) and small InDels were screenedviaa GATK toolkit (v4.1.4.1).Sliding-window analysis on all chromosomes was performed for Δ(wSNP-index) and Δ(InDel-index) calculation with a quantile threshold above 99% as association regions relevant to target traits.Moreover,the transcript expression level was calculatedviathe FPKM (fragments per kilobase of transcript per million fragments mapped) method.The differentially expressed genes (DEGs) between BSR-Seq bulks (CRvs.CS,CRMPvs.CSMP) were assessed using the edgeR (Robinsonet al.2010).The fold change (FC)represented the ratio of expression levels of the same transcript between the two samples.The false discovery rate (FDR) was gained from calibratingP-values (Benjamini and Hochberg 1995).Standards of selecting differentially expressed genes are (fold change)≥2,FDR＜0.05.All DEGs and candidate genes in the association regions were annotated to public databases (GO,COG,KOG,and KEGG).

qRT-PCR verificationSaturated field irrigation (Balbontínet al.2013) was adopted to induce fruit cracking when the third cluster of fruits reached full maturity.Four candidate genes and eight randomly selected genes were selected to conduct the qRT-PCR experiment.Three biological replications were performed in each group,with five mixed samples in each replication.

Total genomic RNA was extracted from tomato pericarp by using Trizol reagent (Biotechnology Co.,Ltd.,Beijing,China).The quality and concentration of RNA were determined by using a Q6000 Nucleic Acid Concentration Assay (Quawell,Los Angeles,USA).CDNA was generated by using a 5× All-In-One RT Master Mix reverse transcription kit (Abm,Vancouver,Canada).The qRT-PCR assay was performed following the instruction manual.

Actinserved as a reference geneThe gene-specific primer sequences were all gained from the qPrimer DB database (Luet al.2018),and the primers were synthesized by Nanjing Springen Biotechnology Company.The 2-ΔΔCtmethod (Livak and Schmittgen 2001) was utilized to obtain the relative levels of gene expression in six genotypes with two cracking traits and two-parent genotypes under 30 h saturated irrigation.Thet-test and SPSS were adopted to analyze the data.

2.4.Analysis of tomato glycoside hydrolase 9 (GH9)gene family

ldentification ofGH9gene family and chromosomal mappingThe protein sequences of the GH9 family inArabidopsis thaliana,Oryza sativaandS.lycopersicumwere downloaded from theArabidopsisinformation resources (TAIR,https://www.arabidopsis.org/),the rice database (http://rice.uga.edu/index.shtml),and the tomato genome website (Solanaceae Genomics Network,https://solgenomics.net),respectively,which were used for BLASTP homology search.In addition,the NCBI CD-Search Tool(https://www.ncbi.nlm.nih.gov/) and Pfam (http://pfam.xfam.org/) databases were adopted to identify sequences of tomato candidateGH9.Genes without complete GH9 domains were discarded.The chromosome annotation information was extracted from genome structure annotation files,and the specific chromosomal distributions were generated with TBtools.

Physicochemical property analysisThe isoelectric point(pI),molecular weights (mWs),number of amino acids (aa),instability index,and average coefficient of hydrophilicity of the 20 SlGH9 amino acid sequences were predicted with ExPASy (https://web.expasy.org/).

SignalP 5.0 (http://www.cbs.dtu.dk/services/SignalP/)and TMHMM-2.0 (http://www.cbs.dtu.dk/services/TMHMM/)were used to predict transmembrane domains and signal peptides.

Phylogenetic analysisMultiple sequence alignment ofA.thaliana,O.sativa,and the SlGH9 proteins were performed using Clustal W.The results were imported into MEGA7.0 to construct the phylogenetic tree by using maximum likelihood (ML) method with 1 000 bootstrap replicates and LG+G model.EvolView online website(https://evolgenius.info//evolview-v2/) was utilized to beautify and edit the evolutionary tree.

Structural characterizationThe conserved structural domain information was obtainedviaNCBI-CDD.The structural patterns of exons and introns were determined using the GSDS online website (http://gsds.cbi.pku.edu.cn/).The conserved motifs of protein sequences were analyzedviaMEME (http://meme-suite.org/).

Analysis ofcis-elements of promotersThe promoter sequences in theSlGH9genes (the 2 000 bp upstream) were predictedviaPlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/).All results were visualized using TBtools.

Expression pattern analysisThe expression levels of 20 members of the tomatoGH9family in four groups(CR,CS,CRMP,and CSMP) were analyzed using BSR sequencing data,and the heat map was drawn using TBtools.Gene expression levels were also visualized using the Tomato Expression Atlas (http://tea.solgenomics.net)and Tomexpress (http://tomexpress.toulouse.inra.fr/).The total pericarp of cultivated tomato ‘M82’ and ‘MicroTom’ fruits was the plant material at different fruit developmental stages,including anthesis,5,10,20,and 30 days after anthesis,respectively,mature green,breaker,pink fruit,light red fruit,and red ripe stages.

2.5.Measurements of cellulase activity and cellulose content

This study essayed the cellulase activity and cellulose content between the cracking-resistant genotype ‘NT91’ and the cracking-susceptible genotype ‘NT189’ at the mature green,turning,red ripe stages according to the methods of Caoet al.(2007) and Vicenteet al.(2007) with minor modifications.

For cellulase activity,the total pericarp was homogenized in 10 mL of pre-chilled 95% ethanol and centrifuged at 4°C for 20 min at 3 000 r min-1.The supernatant was repeated,and the supernatant was combined to form the enzyme extract solution.The 0.5 mL of enzyme extract solution was mixed with 1.5 mL of 10 g L-1carboxymethyl cellulose (CMC)solution at 37°C for 60 min,then 1.5 mL of 3,5-dinitrosalicylic acid (DNS) reagent was added to terminate the reaction.The mixture was incubated in boiling water for 5 min,cooled to room temperature,and the absorbance was determined at 540 nm.The boiled enzyme extract solution was used as a control.

For cellulose content,the total pericarp was placed in 100 mL of 95% ethanol and boiled for 1 h.The material was filtered through a fiberglass filter and sequentially washed with 100 mL of 95% ethanol,chloroform,methanol,and acetone.Then,the residue was dried at 37°C.Finally,cell wall extraction was yielded.Moreover,the residue was washed with acetone,dried at 37°C to a constant weight,and weighed.

3.Results

3.1.Analysis of fruit cracking traits

According to the field observation,the cracking-susceptible tomato fruits were easily cracked during fruit developmental stages,while the cracking-resistant tomato fruits hardly cracked even at the red ripe stage (Fig.1-A).By investigating the cracking rates of two parental genotypes(CR and CS) and two mixed pools selected from 22 crackingresistant plants and 22 cracking-susceptible plants (CRMP and CSMP),it was found that the cracking rates increased with fruit development (Fig.1-A).The CR genotype and CRMP had 0 and 0% cracking rates at the turning stage and reached 1.68 and 3.22% at the red ripe stage,respectively.The fruit cracking rates in the CS genotype and CSMP were 11.25 and 14.78% at the turning stage and surged to 71.03 and 82.15% at the red ripe stage,respectively.Great amounts of fruit dehiscence occurred in the crackingsusceptible genotype from the turning to red ripe stages,which indicated that tomato fruits were extremely prone to dehiscence during these periods.

Frequency statistics were conducted on red ripening fruits in the parents ‘NT189’ and ‘NT91’,as well as 350 F2individuals (Fig.1-B).The result showed that the cracking rate ranges from 0 to 100%,with a mean value of 44.76%and a standard deviation of 21.73%.The asymptotic significance (two-tailed) of the one-sample K-S test was 0.2,which was greater than the significance level of 0.05.This supported the original hypothesis: the sample obeyed a normal distribution,and the irregular cracking trait in tomato was in accordance with the quantitative traits.Moreover,the cracking trait has a pattern of transgressive inheritance.The average cracking rates of crackingresistant and cracking-susceptible plants were 3.22 and 82.15%,respectively.

Fig.1 Fruit cracking traits analysis.A,cracking-resistant (CR,‘NT91’) and cracking-susceptible (CS,‘NT189’) tomato fruits at the mature green (MG),turning (TS),and red ripe stages (RR).B,fruit cracking rates of CR,CS,CRMP,and CSMP.CRMP,CR mixed pools consist of 22 cracking-resistant F2 individuals;CSMP,CS mixed pools are composed of 22 cracking-susceptible F2 individuals.C,frequency distribution histogram of fruit cracking rate in F2 segregation progeny.

3.2.BSR-Seq analysis to discover the crackingrelated regions and candidate genes

Sequencing data analysisRNA sequencing was performed on four groups,including two-parent pools(cracking-resistant pool (CR) and cracking-susceptible pool (CS)) and two F2generation mixed pools (CRMP and CSMP).A total of 57.19 Gb data with Q30 above 92.71%was obtained.The four groups were aligned to the reference genome using tomato SL4.1 as the reference genome.The alignment efficiency was over 93.65%,the average coverage depth was from 161.37× to 213.65×,and the genome coverage was from 14.78× to 20.35× (Appendix A).A total of 116 160 SNPs and 41 591 small InDels were detected among the groups,then 10 533 SNPs and 10 599 small InDels with different genotypes between the CRMP and CSMP groups were chosen for association analysis.

Association analysis of BSR-Seq dataThe distribution of ΔSNP-and ΔInDel-index values between the mixed pools on the whole genome was analyzed by loess regression fitting (Fig.2).The regions above the 99% threshold line (red) were selected as the association regions that were responsible for fruit cracking.As shown in Fig.2-A,on the whole,the distribution of Δ(SNP-index) value on chromosome 9 fluctuated around the 99% threshold line,and its corresponding fitted value (black curve) almost coincided with the 99% threshold line.However,the fitted value of Δ(SNP-index) was noticeably below the 99% threshold line in other loci,e.g.,chromosomes 4 and 8.Similarly,the fitted values of Δ(InDel-index) on chromosomes 9 and 11 were closer to the 99% threshold line than on other chromosomes(Fig.2-B).In sum,three candidate regions were determined with a distinguished peak,in which one region was identified by association analysis with specific SNPs,including 122 genes in chromosome 9,while other regions were identified with specific small InDels,including 22 genes in chromosomes 9 and 11.A total of 127 candidate genes related to fruit cracking traits were obtained by merging the genes in three candidate regions (Table 1).

Table 1 Candidate regions identified by association analysis of BSR-Seq

Fig.2 Calculation of Δ(SNP-index) and Δ(InDel-index) values in the whole genome to identify the candidate genes regulating tomato irregular fruit cracking.The scattered points represent the Δ(SNP-index) or Δ(InDel-index) value.The black curve represents the corresponding fitted value,and the red dash line represents the 99% threshold line set by Loess regression.A,candidate regions were located by the association analysis of specific SNPs between the two segregate bulks of BSR-Seq(CRMP and CSMP).B,candidate regions were located by the association analysis of specific small InDels between CRMP and CSMP.CRMP,cracking-resistant (CR) mixed pools consist of 22 CR F2 individuals;CSMP,cracking-susceptible (CS) mixed pools consist of 22 CS F2 individuals.

A total of 127 candidate genes within association regions were annotated to GO and KEGG public databases to investigate the key genes that regulate fruit cracking (Fig.3).The GO enrichment analysis showed that 127 candidate genes were involved in biological processes such as nitric oxide-mediated signal transduction,iron ion homeostasis,positive regulation of fatty acid biosynthesis,and cellulose microfibril construction.KEGG pathway revealed that the candidate genes were remarkably enriched in the metabolic processes (starch and sucrose;galactose) and biosynthetic processes (cutin,suberine,and wax;glycan;pantothenic acid and coenzyme A).The results suggested that irregular fruit cracking in tomatoes may be associated with cell wall construction,cuticle biosynthesis,certain biochemical processes (glycan,acid biosynthesis,and metabolism),and signal transduction pathways.

Fig.3 GO enrichment and KEGG pathway analysis of candidate genes identified by the association analysis of specific SNPs (A and B) and small InDels (C and D).

Screening differentially expressed genes,GO and KEGG analysisThe differential expression analysis revealed that there were 7 738 DEGs in the CR and CS genotypes,of which 5 715 were up-regulated,and 2 023 were downregulated (Appendix B).The number of up-regulated genes was higher than that of down-regulated genes in the CS compared to the CR.There were 4 111 DEGs in the CRMP and CSMP,of which 3 230 genes were up-regulated,and 881 were down-regulated.The number of up-regulated genes was higher than the number of down-regulated genes in the CRMP compared to the CSMP.

The GO and KEGG annotations were performed on DEGs from two sets of sequencing bulks (CRvs.CS,CRMPvs.CSMP) to enhance the understanding of DEGs’ trends in gene functions.The GO enrichment analyses of the screened DEGs in CR and CS genotypes were performed(Fig.4-A).The results showed that in the biological process category,2 923 and 2 670 DEGs were dramatically enriched in cellular and metabolic processes,and 2 525 and 2 295 DEGs were remarkably enriched in binding and catalytic activities in the molecular function category,respectively.Moreover,DEGs in CR were mainly enriched in fatty acid metabolism,protein phosphorylation,and DNA repair processes;those in CS were mainly enriched in cell wall composition and modification,biogenesis,xyloglucan metabolism,and lipid storage processes.The KEGG analysis of DEGs indicated that they were remarkably enriched in metabolic (nitrogen;glycan),biosynthetic processes (fatty acid;cutin,suberine,and wax;phenyl propane;amino acid);and other pathways (Fig.4-B).In general,they were mainly involved in processes relevant to cuticle biosynthesis and cell wall metabolism.

The GO analysis of the screened DEGs in the CRMP and CSMP showed that 1 665 and 1 612 DEGs were largely enriched in cellular and metabolic processes,and 1 413 and 1 411 DEGs were enriched in catalytic and binding activities,respectively.Moreover,DEGs in CRMP were mainly enriched in processes such as respiratory chain and iron ion homeostasis;those in CSMP were mainly enriched in processes such as hydrogen peroxide metabolism,sterol biosynthesis,negative regulation of apoptosis,and cell wall construction (Fig.4-C).The KEGG enrichment analysis revealed that DEGs were significantly enriched in biosynthetic pathways (phenol propane;cutin,suberine,and wax),metabolic pathways (sugar;amino acid;glutathione),and phytohormone signaling pathways(Fig.4-D).Therefore,DEGs were probably involved in the processes of cell wall construction,redox reaction,and phytohormone responsiveness.

Fig.4 Graph of GO enrichment and KEGG pathway analysis of differentially expressed genes (DEGs) between two set pairs of RNA sequencing bulks.A and B,CR vs.CS.C and D,CRMP vs. CSMP.CR,cracking-resistant;CS,cracking-susceptible;CRMP,CR mixed pools consist of 22 CR F2 individuals;CSMP,CS mixed pools consist of 22 CS F2 individuals.

qRT-PCR validationA total of four crossed genes were screened by combining the DEGs in two RNA sequencing bulks (CRvs.CS,CRMPvs.CSMP) with candidate genes in association regions.They wereSolyc09g010060(kinesinlike protein),Solyc09g010080(beta-fructofuranosidase),Solyc09g010210(EG,β-1,4-endoglucanase),andSolyc09g010230(B3 domain-containing protein).The qRT-PCR was performed using the red ripening fruits of parent genotypes at 30 h of saturated irrigation.The results exhibited that only the expression level ofSolyc09g010210was distinctly higher in CS at 30 h of saturation irrigation than in CR (Fig.5-A).Solyc09g010210,a member of theglycoside hydrolase 9gene family,participated in the metabolic process of cellulose (a major component of the cell wall),which exerted its function on degrading the cell wall.This study also investigated the expression levels of eight randomly selected genes of association region in the three cracking-susceptible genotypes and three crackingresistant genotypes (Fig.5-B).The results in qRT-PCR were accorded with the sequencing data,which confirmed the reliability of BSR-Seq.

Fig.5 qRT-PCR analysis of differentially expressed genes (DEGs).A,qRT-PCR verification of four crossed genes.Relative transcript levels are expressed relative to the tomato ACTIN gene internal control,expressed as 2-ΔΔCt.Solyc09g010060,kinesin-like protein;Solyc09g010080,beta-fructofuranosidase;Solyc09g010210, EG,β-1,4-endoglucanase;Solyc09g010230,B3 domain-containing protein.B,qRT-PCR verification of candidate gene expression in cracking-sensitivity (CS) genotypes‘NT189’,‘NT48’,and ‘NT177’ and cracking-resistant (CR) genotypes ‘NT91’,‘NT54’,and ‘NT191’.Solyc09g009990,auxin-induced protein 15A-like;Solyc09g010010,ACO1-like,1-aminocyclopropane-1-carboxylate oxidase 1 like;Solyc09g010020,ACO1-like;Solyc09g010000,ACO1;Solyc09g009420,unknown;Solyc09g009940,signal recognition particle protein;Solyc11g071530,unknown;Solyc11g071510,unknown.Data were the mean of three biological repeats±standard error (SE).Different letters represented significant differences (P＜0.05).

3.3.Analysis of the tomato glycoside hydrolase 9(GH9) gene family

Genome-wide identification and chromosomal localization of GH9 gene familyBased on the tomato SL4.0 reference genome,a total of 26 tomatoGH9genes with the structural domain of glycoside hydrolase 9 were confirmed with SMART and CDD tools.Genes without complete GH9 domains were discarded.Finally,20 tomatoGH9genes were screened.The distribution ofSlGH9son chromosomes was mapped using TBtools according to the gene location information (Appendix C).The 20SlGH9swere found to be distributed on 11 chromosomes except for chromosome 10.Based on their locations on the chromosome,SlGH9swere designated asSlGH9-1toSlGH9-20.The candidate geneSolyc09g010210was identified asSlGH9-15.

Physicochemical properties analysis ofSlGH9 genesAs presented in Table 2,the number ofGH9amino acids ranged from 450 to 633 in tomato.The amino acids and molecular weight ofSlGH9-15were 490 aa and 54 118.36 KDa,ranging in the lower middle of the 20 genes.The isoelectric point of SlGH9 protein ranged from 5.04 (SlGH9-20) to 9.05(SlGH9-16),and the isoelectric point of SlGH9-15 protein was 8.47,which was greater than seven and belonged to the basic protein.The instability indexes were between 25.28(SlGH9-13) and 45.59 (SlGH9-2),indicating that most GH9 proteins are stable.The grand average of hydropathicity values indicated that allSlGH9genes encoded hydrophilic proteins.Meanwhile,10 and 11 SlGH9 proteins contained transmembrane structures and signal peptides,respectively.The SlGH9-15 protein contained zero transmembrane domain and one signal peptide.

Table 2 Properties and locations of predicted GH9 proteins in tomato1)

Phylogenetic analysisTo further explore the evolutionary relationships and homology ofGH9family members,we constructed a phylogenetic tree with 26A.thaliana,24O.sativa,and 20 tomatoGH9gene-encoded proteins using the ML method in MEGA 7 (Fig.6).Twenty tomatoGH9genes were composed of three branches (clades A,B,and C),including three clade A genes (SlGH9-1,SlGH9-7,SlGH9-17),13 clade B genes (SlGH9-2,SlGH9-4,SlGH9-5,SlGH9-6,SlGH9-8,SlGH9-9,SlGH9-10,SlGH9-11,SlGH9-12,SlGH9-14,SlGH9-15,SlGH9-16,andSlGH9-18),and four clade C genes (SlGH9-3,SlGH9-13,SlGH9-19,andSlGH9-20),according to the structural differences,clustering results,and available taxonomic information.Clade A contained threeA.thalianagenes and threeO.sativagenes,clade B contained 21A.thalianagenes and 18O.sativagenes,and clade C included twoA.thalianagenes and threeO.sativagenes.These results indicated that the three clades were formed before the differentiation of monocotyledons and dicots,and theGH9gene family is highly conserved in plants.SlGH9-15belongs to clade B and has high homology withAtGH9B1,AtGH9B6,OsGH9B1,andOsGH9B3,indicating that they might have similar biological functions.

Fig.6 Phylogenetic tree of GH9 proteins from tomato,Arabidopsis thaliana,and Oryza sativa.An unrooted tree was constructed using the maximum like (ML) method with 1 000 bootstrap replicates in MEGA 7.The red solid stars represent the SlGH9 genes of tomato,the green solid squares represent the GH9 genes of A.thaliana,and the yellow triangle represents that of O. sativa.Different clades are shown in different colors;clades A,B,and C were highlighted in pink,salmon,and blue,respectively.

Structural characterization analysisThe gene structure of tomatoGH9Aswas relatively conserved,with the glycosyl hydrolase structural domain distributed on four or six exons,and that in clade C was five or ten.The gene structure of members of clade B was more variable,distributed from four to eight exons,which may be related to their complex and variable functions.However,SlGH9-20contains only one exon,which is presumed to be originally derived from the retrotransposon.Meanwhile,the structural distribution of introns and exons ofSlGH9showed that clade A genes had four introns,clade C genes had four or six introns,while clade B genes had 0,two,or four introns (Fig.7-C).Further analysis of the conserved motifs revealed that all 20 identified conserved motifs were within the glycoside hydrolase structural domain,indicating that this structural domain is relatively conserved in tomato (Fig.7-B).Most of theSlGH9genes in each clade contained 14 conserved motifs,withSlGH9-13andSlGH9-20containing the largest number of 17 motifs.In addition toSlGH9-17,genes of clade A contained specific motif 18,genes of clade C contained specific motif 14 and motif 19,those in clade B contained approximately the same kinds of motifs,and some members contained specific motif 17.The results indicated that the members of the same clade have similar motifs,and unique motifs in different clades also enhance the support of the evolutionary tree (Fig.7-A).

Fig.7 The structure characterization of GH9 genes in tomato.A,phylogenetic tree of GH9 in tomato.Green round,pink triangle,and blue diamonds represent SlGH9As,SlGH9Bs,and SlGH9Cs,respectively.B,distribution of 20 conserved motifs.C,distribution of exons and introns.Coding sequence (CDS) and glyco hydro 9,colored in green and yellow,represent exons,the black line represents intron,and the numbers represent the numbers of introns.UTR,untranslated region.

Analysis of cis-acting elements in the promoter regions of SlGH9 genesTo further understand the regulatory mechanisms ofSlGH9sin response to various hormones and stresses,the 2 000 bp promoter regions ofSlGH9sgenes were analyzed forcis-acting elements prediction using PlantCARE.Two kinds ofcis-elements,hormone-responsive and biotic and abiotic stressresponsive elements,were visualized and demonstrated in Fig.8.Initially,10 hormone-responsive elements were substantially found in promoter regions,which were composed of AuxRR-core,TGA element (IAA,auxin responses),ABRE (ABA,abscisic acid responses),ERE(ethylene responses),GARE-motif,P-box,TATC-box(GA,gibberellin responses),TGACG-motif,CGTCA-motif(MeJA,methyl jasmonate responses),and TCA-element(SA,salicylic acid responses).A total of 101 of these elements were associated with ABA,ethylene,and IAA responses.Then,four stress-responsive elements were analyzed,namely MBS (drought-induced),LTR(low temperature-induced),ARE (anaerobic-induced),and TC-rich repeats (disease and stress-induced).In addition,most of theSlGH9sgenes also have elements that act in relation to MYB binding sites,light response,circadian control,seed regulation,meristem expression,endosperm expression,flavonoid,and zein metabolism.Among them,theSlGH9-15promoter hascis-elements that are relevant to phytohormone (IAA,ethylene,ABA,MeJA,and SA) responses,anaerobic induction,drought induction,and flavonoid biosynthesis gene regulation.This implied thatSlGH9-15gene may play a vital role in some stress and hormone response processes by participating in distinctive regulatory pathways,which in turn affects fruit cracking (Table 3).

Table 3 cis-Acting elements of SlGH9-15 promoter regions

Fig.8 cis-Acting element analysis in the promoters of SlGH9 genes.

Expression pattern analysisUsing BSR data,a heat map of the mRNA expression levels of 20SlGH9genes in CR,CS,CRMP,and CSMP was created (Fig.9-A).ThirteenSlGH9genes were found to be more highly expressed in cracking-susceptible tomatoes than in cracking-resistant tomatoes.Among them,SlGH9-1andSlGH9-15had higher expression (FPKM＞30).SlGH9-1was highly expressed in cracking-resistant tomatoes,whileSlGH9-15was highly expressed in crackingsusceptible tomatoes.The expression levels ofSlGH9genes at different developmental stages of fruit were visualized using Tomato Expression Atlas (http://tea.solgenomics.net) and Tomexpress (http://tomexpress.toulouse.inra.fr/) (Fig.9-B and C).It was found that the expression levels ofSlGH9-6,SlGH9-7,andSlGH9-13were significantly high at the immature green to mature green stage and decreased dramatically during fruit ripening.While the expression level ofSlGH9-15was quite low at the mature green stage,it gradually increased from the breaker stage and reached the maximum in red ripe,which is consistent with the increase of fruit cracking rate in cracking-susceptible tomatoes.The differences in expression level between BSR-Seq and website data might be due to the differences in varieties,treatments,and screening standards.

Fig.9 Expression patterns of SlGH9 gene in tomato.A,expression patterns of SlGH9 gene from BSR-Seq data.B,expression cubes of SlGH9 gene in cultivated tomato ‘M82’ at different fruit developmental stages.C,expression patterns of SlGH9 gene in‘Micro Tom’ at different fruit developmental stages.CR,cracking-resistant;CRMP,CR mixed pools consist of 22 CR F2 individuals;CS,cracking-susceptible;CSMP,CS mixed pools consist of 22 CS F2 individuals.DPA,days post anthesis.

3.4.The cellulase (Cx) activity and the cellulose content between CR and CS genotypes during fruit development

By measuring Cx activity and cellulose content in mature green,turning stages,and red ripening fruits of CR and CS genotypes (Fig.10),we found that Cx activity tended to increase,and cellulose content tended to decrease during fruit development.In all three fruit developmental stages,Cx activity in CS genotype ‘NT189’ was higher than that of CR genotype ‘NT91’,while the cellulose content was lower than that of CR tomato,with significant differences.The activity of Cx in irregularly cracked tomatoes was 213,433,and 911 μmol h-1g-1at the mature green,turning,and red ripe stages,respectively,which were 1.58,1.17,and 2.14 times higher than that in CR tomato.Moreover,cellulose content in CR fruits at the mature green,turning,and red ripe stages were 1.41,1.57,and 2.07 times higher than that of CS fruits.The massive differences between CR and CS revealed that high cellulase activity might decrease cellulose content and then accelerate the disintegration of the fruit cell wall in cracking-susceptible tomato,subsequently triggering fruit cracking.

Fig.10 The cellulase activity (A) and cellulose content (B) between cracking-resistant (CR,‘NT91’) and cracking-susceptible (CS,‘NT189’) tomato fruits at the mature green (MG),turning (TS),and red ripe (RR) stages.Data were the mean of three biological repeats±standard error (SE).Different letters represented significant differences (P＜0.05).

4.Discussion

The tendency of fruits to dehiscence has been correlated with heredity,fruit characteristics (e.g.,fruit size,shape,hardness;exocarp mechanical strength,water absorption capacity;number and distribution of stomata),external environmental conditions (e.g.,humidity,temperature,and light),and cultivation management conditions (e.g.,irrigation and fertilization).It is not a sole factor but various factors working collectively to influence fruit dehiscence (Hahn 2011;Balbontín 2013;Khadivi-Khub 2015).Previous studies also illustrated that fruit cracking is heritable (Yamaguchiet al.2002;Correiaet al.2018).However,it seems that no single gene but many genes may contribute to the inheritance of fruit cracking (Vaidyanathan and Harrigan 2005).Moreover,plant hormone,cell wall,and possibly the cuticle layer were involved in fruit cracking (Cosgrove 2000;Liet al.2003;Bargel and Neinhuis 2004;Demirsoy and Demirsoy 2004;Gine-Bordonabaet al.2017;Liaoet al.2020;Xueet al.2020).It has been identified that genes such asSlTBG6,SlPG,SlEXP1,SlER4.1,ClERF4,andSlMAPKKKare associated with fruit cracking,using QTL,BSA-Seq,gene function identification,and GWAS (Moctezumaet al.2003;Cuiet al.2017;Jianget al.2019;Liaoet al.2020;Chenet al.2021).Research findings varied regarding key genes that regulate fruit cracking due to the different plant materials,environmental conditions,and experimental approaches.The molecular mechanisms of fruit cracking are still poorly understood and require further investigation.

In this study,the BSR-Seq was constructed with the F2segregation individuals to investigate the molecular mechanism of tomato fruit cracking at the genetic and transcriptional levels.Two regions associated with irregularly fruit cracking traits were located on chromosomes 9 and 11,which varied from previous locations,and should be new sites for fruit cracking.In previous research,Li (2016)detected five QTLs on chromosomes 1,2,and 3 associated with crack resistance property index from cracking-tolerant processed tomato ‘14803’.Capelet al.(2017) identified QTL related to cracking in tomatoes on chromosomes 1,3,8,10,and 12.Zhuet al.(2020) localized the crack tolerance geneCr3ain a 349-kb interval on chromosome 3,containing the cell wall-related genesα-1,6-xylosyltransferaseandgalactosyltransferase.Chenet al.(2021) combined BSA sequencing and genetic linkage analysis,locating the tomato irregular cracking trait in a 387-kb interval on chromosome 2.

The localization of the cracking trait at different chromosomes might be due to the different plant materials and treatments.In this study,a total of 127 genes were obtained in the two candidate regions,and four genes were screened by differential expression analysis.Subsequently,qRT-PCR analysis was conducted in the red ripening fruits of parent genotypes at 30 h of saturated irrigation.Finally,Solyc09g010210with significantly differential expression was screened out,which encodesβ-1,4-endoglucanasegene,hydrolyzing polysaccharides with a glucan backbone.It was involved in cellulose degradation and exerted a great influence on cell wall remodeling.Pericarp cell wall composition and biochemical modifications affect fruit firmness and mechanical properties of the pericarp (e.g.,elasticity,plasticity,and toughness),which in turn impact the occurrence of fruit cracking (Considine and Brown 1981;Chapmanet al.2012;Balbontínet al.2013;Gaoet al.2021).The fruit cell wall consists of cellulose,hemicellulose,and pectin.Cellulose is an important component of the cell wall and constitutes the skeletal structure of the cell wall.Therefore,we hypothesized thatβ-1,4-endoglucanasemight be a candidate gene for the regulation of irregular fruit cracking in tomatoes.The evolutionary tree analysis revealed that the tomatoSlGH9-15gene is a close affinity withArabidopsis AtGH9B1andAtGH9B6and riceOsGH9B1andOsGH9B3.Previous studies found that silencingAtGH9B1resulted in cell wall ruffling with reduced cellulose and lignin content (Tsabaryet al.2003).OsGH9B1andOsGH9B3were characterized to have cellulase activity for decreasing cellulose crystallinity (Huang J Fet al.2017).Overexpression ofOsGH9B1andOsGH9B3slightly changed cell wall content but remarkably increased cellulase activity and reduced cellulose polymerization and crystallinity index,which in turn modified cellulose microfibrils (Huanget al.2019).Therefore,it should suggest that tomatoSlGH9-15also has an enzymatic activity affecting cellulose content,lignin content,cellulose polymerization,and crystallization to modify the cell wall.The expression pattern analysis showed thatSlGH9-15was more highly expressed in the cracking-susceptible genotype than in the cracking-resistant genotype.The expression level ofSlGH9-15was quite low at the mature green stage,and it was dramatically upregulated from the breaker stage to the red ripe stage.This might be due to the fact that gene expression precedes the phenotype.The differences in fruit cracking rate between CR and CS during fruit development revealed that the cracking-susceptible genotype is highly prone to dehiscence from the turning to red ripe stages.The fruit cracking rates of CS under the mature green,turning,and red ripe stages were 1.25,11.25,and 71.03%,while those of CR were 0,0,and 1.68%,respectively.Meanwhile,compared with CR fruits,CS fruits had higher cellulase activity and lower cellulose content at the mature green,turning,and red ripe stages,especially at the red ripe stage.Compared with CR,the expression level ofSlGH9-15was dramatically upregulated in the CS genotype.Its homologousOsGH9B1andOsGH9B3remarkably increased cellulase activity.In this experiment,we also observed increased cellulase activity and decreased cellulose content,which can accelerate the disassembly of the cell wall,in turn facilitating fruit cracking.Thus,we suppose thatSlGH9-15plays a critical role in fruit cracking.

Previous studies have also shown that ethylene signaling can affect fruit development and ripening by regulating the expression ofPGandEXPgenes related to cell wall degradation (Giovannoni 2004;Gaoet al.2020).Chenet al.(2021) identifiedMAPKKK17in a 387-kb interval on chromosome 2 as a candidate gene for irregular fruit splitting in tomato,which is associated with the signaling pathway of ethylene-responsive protein kinaseCTR1.ClERF4was detected on chromosome 10 and confirmed to be a key gene regulating watermelon rind hardness and its crack resistance.Additionally,ethylene and auxin signals are cross-linked.Auxin activates the production of ROS and H2O2crumbled cell wall polymers by producing OH-,while cell wall degradation decreases rind elasticity and toughness,which may lead to fruit cracking (Liszkayet al.2004;Liet al.2006;Xueet al.2020;Shiet al.2021).The present study found that the promoter region ofSlGH9-15contained severalcis-acting elements related to phytohormones (e.g.,ethylene,auxin,and ABA) and stress responses (e.g.,drought and anaerobic induction),indicating thatSlGH9-15may be induced by hormones such as ethylene and auxin and certain abiotic stresses.We can speculate that the cell wall-associatedSlGH9-15and phytohormone (ethylene,auxin) signaling pathways may jointly regulate tomato fruit dehiscence.

5.Conclusion

This study used BSR-Seq analysis to locate two regions associated with irregularly cracking on chromosomes 9 and 11,containing 127 candidate genes.Global gene expression analysis demonstrated that candidate genes enriched in cell wall construction,cuticle biosynthesis,some biochemical processes,and signal transduction pathways.Further,through differentially expression analysis and qRTPCR in cracking-susceptible and -resistant genotypes,theSlGH9-15(Solyc09g010210) was identified.Bioinformatics analysis revealed thatSlGH9-15belongs to clade B of the GH9 gene family and is closely related to the cell wall construction-associated genesAtGH9B1,AtGH9B6,OsGH9B1,andOsGH9B3.Thecis-acting elements showed thatSlGH9-15was activated by various phytohormones (ethylene and ABA) and abiotic stresses.The expression pattern demonstrated that the expression level ofSlGH9-15was relatively low at the mature green stage,gradually increased during fruit ripening,and reached its maximum at the red ripe stage.Simultaneously,the fruit cracking rates of cracking-susceptible genotype were remarkably increased from the turning to red ripe stages.During these periods,cracking-susceptible fruits had enhanced cellulase activity and decreased cellulose content.This study identifiedSlGH9-15as the key gene involved in fruit cracking in tomatoes for the first time.Our findings can provide the foundation for understanding the molecular mechanism and complex regulatory network of fruit cracking in tomatoes.

Acknowledgements

This work was supported by the National Key Research and Development Program of China (2019YFD100190200),the Jiangsu Agricultural Science and Technology Innovation Fund,China (CX(20)3101),and a grant from the Fundamental Research Funds for the Central Universities,China (KYZZ2022004).

Declaration of competing interest

The authors declare that they have no competing interests.

Appendicesassociated with this paper are available on http://www.ChinaAgriSci.com/V2/En/appendix.htm