Genome-wide identification of GH9 gene family and the assessment of its role during fruit abscission zone formation in Vaccinium ashei

The genomic location and stage-specific expression pattern of GH9 genes reveal their critical roles during fruit abscission zone formation in Vaccinium ashei. Glycosyl hydrolase family 9 (GH9) cellulases play a crucial role in both cellulose synthesis and hydrolysis during plant growth and development. Despite this importance, there is currently no study on the involvement of GH9-encoding genes, specifically VaGH9s, in abscission zone formation of rabbiteye blueberries (Vaccinium ashei). In this study, we identified a total of 61 VaGH9s in the genome, which can be classified into 3 subclasses based on conserved motifs and domains, gene structures, and phylogenetic analyses. Our synteny analysis revealed that VaGH9s are more closely related to the GH9s of Populus L. than to those of Arabidopsis, Vitis vinifera, and Citrus sinensis. In silico structural analysis predicted that most of VaGH9s are hydrophilic, and localized in cell membrane and/or cell wall, and the variable sets of cis-acting regulatory elements and functional diversity with four categories of stress response, hormone regulation, growth and development, and transcription factor-related elements are present in the promoter sequence of VaGH9s genes. Transcriptomic analysis showed that there were 22 differentially expressed VaGH9s in fruit abscission zone tissue at the veraison stage, and the expression of VaGH9B2 and VaGH9C10 was continuously increased during fruit maturation, which were in parallel with the increasing levels of cellulase activity and oxidative stress indicators, suggesting that they are involved in the separation stage of fruit abscission in Vaccinium ashei. Our work identified 22 VaGH9s potentially involved in different stages of fruit abscission and would aid further investigation into the molecular regulation of abscission in rabbiteye blueberries fruit.


Introduction
Abscission is a natural process by which organs are shed through cell separation at specific sites called abscission zones (AZs).AZs are narrow bands of functionally specialized cell layers that respond to a range of cues such as developmental, hormonal, and environmental signals.Abscission occurs via at least four stages: (i) pre-formation stage in 1 3 which cells are differentiated to form AZs at the future sites of organ detachment (Estornell et al. 2013;Ma et al. 2021); (ii) induction stage in which AZs cells become competent in response to abscission signals and activated to begin elongation (Patharkar and Walker 2018); (iii) separation stage in which cell wall degradation enzymes are produced to weaken cell adhesion (Roberts et al. 2002); and (iv) protective layer formation stage in which protective compounds are secreted to cover the wound surface of the separated cells and shield from desiccation (Estornell et al. 2013).The regulation of the rate and degree of abscission is heavily influenced by the delicate balance and interplay of phytohormones.Ethylene and abscisic acid are known to promote abscission, while gibberellins and auxins inhibit it (Botton and Ruperti 2019;Colle et al. 2019).Interestingly, jasmonic acid has been found to have a distinctive effect on different plants; it stimulated fruit abscission in grapes, while it delayed ethylene-induced petal abscission in roses (Fidelibus et al. 2022;Priya et al. 2022;Sawicki et al. 2015).The abscission process involves multiple cell wall remodeling enzymes including cellulase, an endo-β-1,4-glucanase that is widely found in both plants and cellulolytic microorganisms (Kalaitzis et al. 1999).Glycoside hydrolase family 9 (GH9) genes (http:// www.cazy.org/), which primarily code for endo-β-1,4-glucanases, have been identified in various plant species, including, but not limited to, Arabidopsis thaliana, Oryza sativa, Populus L., and these genes were shown to be associated with diverse aspects of plant growth and development (Cantarel et al. 2009;Dheilly et al. 2016;Xie et al. 2013;Yu et al. 2014).The crystal structure showed a high degree of conservation in both secondary and tertiary structures, as well as catalytic and carbohydrate-binding amino acid residues in certain species (Gray et al. 2018).Plant-derived GH9 enzymes typically consist of 441-694 amino acid residues, with a calculated and observed molecular mass ranging from 42 to 70 kDa.Their structures feature 12 alpha helices that are comprised of 6 inner and outer alpha helices forming an α/α barrel shape with 6 parallel alpha helices (Du et al. 2015;Yu et al. 2014).The GH9 family proteins can be divided into three distinct structural subclasses, with GH9A having a single NH 2 -terminal transmembrane helix and a catalytic domain, GH9B having a signal peptide and a catalytic domain, and GH9C having a COOH-terminal carbohydrate-binding module (CBM) in addition to a signal peptide and catalytic domain (Del et al. 2012;Urbanowicz et al. 2007).
The GH9 genes in higher plants have evolved specific functions such as cell division and elongation, cytokinesis, anther dehiscence, pollen tube growth, organ abscission, and fruit ripening (Campillo et al. 2012;Szyjanowicz et al. 2004;Yu et al. 2014;Zhou et al. 2006).Although many GH9 proteins have been functionally studied in multiple plants, including Arabidopsis thaliana (Zuo et al. 2000), Oryza sativa L. (Xie et al. 2013), Populus L. (Yu et al. 2014), Triticum (Szyjanowicz et al. 2004), and Fragaria ananassa (Luo et al. 2022), their roles in cellulose synthesis and assembly remain poorly understood.All 74 Triticum GH9 genes (TaGH9s) identified encode hydrophilic proteins, and subcellular localization predictions revealed that TaGH9 proteins act in the cytoplasm, Golgi complex, mitochondria, endoplasmic reticulum, and nucleus (Szyjanowicz et al. 2004).Previous studies have found that individual genes in GH9 family are involved in the regulation of cellulose content and cellulase activity.For example, OsGH9B1, 3, and 16, are associated with cellulase activity according to the catalytic domain analysis (Xie et al. 2013).Cellulase activity was promoted by the increased expression of CEL1 and CEL2 genes during fruit development in Fragaria ananassa (Mercado et al. 2010).However, mutation of OsGLU1 in Oryza sativa resulted in the decreased cellulose content, suggesting that OsGLU1 affects the cellulose biosynthesis of Oryza sativa (Campillo et al. 2012).The KOR1 in GH9 family, which is located primarily in the plasma membrane, plays an important role in cellulose metabolism of both primary and secondary cell walls in Arabidopsis and Populus L. (Yu et al. 2014;Zuo et al. 2000).
Rabbiteye blueberry (Vaccinium ashei) has been increasingly cultivated in China due to their good adaption to high temperatures and acidic (pH 4.0-5.3)soil conditions (He et al. 2009).This fruit is known for its potent benefits for consumers (Davidson et al. 2018;Mercado et al. 2010) and has unique fruit development and ripening, with all berries in the same inflorescence taking at least 1 month to fully abscise from the main plant body.In recent years, abscission has become a focus of research due to its significant negative impact on crop yield and fruit quality when it occurs prematurely.Additionally, in some crops like cotton and citrus, accelerated abscission at a proper time is necessary to facilitate improved harvest and labor savings.As a result, many abscission studies have been carried out in fruit crops including, but not limited to, citrus (Sabbione et al. 2019), apple (Arseneault et al. 2016), litchi (Wang et al. 2019), sweet cherry (Qiu et al. 2020), blueberry (Vashisth et al. 2013;Vashisth et al. 2015) orange (Merelo et al. 2017) and tomato (Brummell et al. 1999).It has been pointed out that fruit abscission mainly results from the disassembly of cell layers at fruit AZs which is driven by cell wall-degrading enzymatic activities (Liu et al. 2019;Wang et al. 2019).As the fruit matures, there is an increasing tendency of cellulase activity in fruit AZs (Li et al. 2019a;Merelo et al. 2017).In recent years, 25 GH9 genes in Arabidopsis, 24 GH9 genes in rice, 22 GH9 genes in barley, 29 GH9 genes in corn, 25 GH9 genes in poplar, and 7 GH9 genes in wheat have been identified (Buchanan et al. 2012;Du et al. 2015;Xie et al. 2013;Yu et al. 2014).To date, there is no systematic study of the GH9 gene family in Vaccinium ashei despite the crucial role of cellulase in abscission.The recent completion of Vaccinium ashei genomic release allowed us to study the GH9 gene family at the genome level (Colle et al. 2019;Main et al. 2018).In this article, we investigated the gene structure, chromosomal location, conserved motifs and domains of Vaccinium ashei GH9 genes (VaGH9s), conducted phylogenetic and syntenic analyses between GH9 genes in Arabidopsis thaliana, Populus L., Vitis vinifera, Citrus sinensis, and Vaccinium ashei, and analyzed the expression of VaGH9s in abscission zone during fruit development in the presence and absence of ethephon and gibberellins.This study presents the first report on genome-wide VaGH9s profile analysis and the differential expression of VaGH9s in fruit AZs during fruit development of Vaccinium ashei.

Plant material and in vitro treatments
Ten-year-old rabbiteye blueberry (Vaccinium ashei) cv.'Brightwell', grown in Zhejiang A&F University was used in this study.During the veraison stage, the trees were treated with either gibberellin (200 mg/L), ethephon (500 mg/L), or water (served as control).The berries at relatively consistent maturity and without mechanical damage were picked at 7-day intervals.The abscission zone (AZ) tissues at fruitpedicel junction were immediately isolated and brought to the Key Laboratory of Quality and Safety Control for Subtropical Fruit and Vegetable, Ministry of Agriculture and Rural Affairs for the analysis of microstructure and the measurement of cellulase activity, reactive oxygen species (ROS), malondialdehyde (MDA), and relative conductivity.The AZ tissues from the berries at the veraison stage were isolated 0 h and 24 h after gibberellin (200 mg/L) or ethephon (500 mg/L) treatment and frozen into liquid nitrogen for transcriptome analysis.

Collection and identification of VaGH9 genes in Vaccinium ashei
The reference genome and gene model annotation files were retrieved from the blueberry genome website (https:// www.vacci nium.org/ crop/ blueb erry) (Colle et al. 2019).To obtain all the GH9 sequences in Vaccinium ashei, we performed BLASTn homology searches against the Vaccinium ashei genome datasets using the full-length coding sequences of Arabidopsis GH9s as queries.Candidate GH9 genes were selected if it was satisfied with E < 10 −10 .Subsequently, the deduced amino acid sequences of candidate genes were analyzed to detect glycosyl hydrolase family 9 (IPR001701) domains from the InterPro protein signature database (https:// www.ebi.ac.uk/ inter pro/) (Mistry et al. 2021), and the related and aligned sequences were retrieved from the Vaccinium ashei genome datasets using HMMER3.0software (https:// www.ebi.ac.uk/ Tools/ hmmer/) (Potter et al. 2018).We provisionally named the candidate sequences VaGH9A1 to VaGH9C12 following the standardized nomenclature of the GH9 family in Arabidopsis to distinguish them (Urbanowicz et al. 2007) (Table S1).

Exon/intron structure, conserved motifs, domains and promoter analysis
The exon-intron arrangement was visualized with the gene structure display server (GSDS) using both genomic DNA sequences and the corresponding coding sequences (http:// gsds.gao-lab.org/) (Hu et al. 2015).Conserved motifs of GH9 proteins were analyzed using the MEME Version 4.11.4 program (http:// meme.nbcr.net/ meme/ tools/ meme).The following parameters was used: motif width set to 6-200 and the maximum number of motifs set to 10 ( Bailey et al. 2009).The MEME motifs were annotated using the Inter-Pro database (https:// www.ebi.ac.uk/ inter pro/), and the conserved domain database was used to predict the conserved domains of the GH9 family members of Vaccinium ashei.The 2.0-kb DNA sequences upstream of ATG were retrieved from the genome database to predict the cis-acting regulatory elements using PlantCARE (http:// bioin forma tics.psb.ugent.be/ webto ols/ plant care/ html/) (Lescot et al. 2002).

Chromosomal location and gene synteny analysis
The positions of GH9-encoding genes on the chromosomes of Vaccinium ashei and aforementioned plants were obtained from the NCBI (https:// www.ncbi.nlm.nih.gov).The synteny analysis among members of Vaccinium ashei GH9 family and between Vaccinium ashei and Arabidopsis, Populus L., Vitis vinifera, and Citrus sinensis were performed using MCScanX and TBtools (Chen et al. 2020;Wang et al. 2012).

RNA isolation and transcriptome sequencing analysis
The total RNAs from AZs samples were extracted using #RC401 test kit according to the manufacturer's instructions (Vazyme Biotech Co., Ltd), and RNA quality was verified with an Agilent 2100 Bioanalyzer.Transcriptome libraries were constructed following the standard Illumina RNA-seq protocol (Illumina, Inc., San Diego, CA, catalog no.RS-100-0801).RNA and DNA libraries were sequenced using a 150-bp paired-end Illumina sequencing, with libraries having inserts of 350 bp.For RNAseq, the raw reads were filtered using fastp (v 0.12.2) with default parameter to eliminate those reads with poly-N, adapter, and low qualities (Chen et al. 2018).The filtered high-quality reads were mapped on the V_corymbosum_ genome_v1.0genome (Kulkarni et al. 2020) using Hisat2 (v 2.1.0)with parameters (-dta, -x) followed by gene annotation in GTF format (Trapnell et al. 2012).Genes with a minimum FPKM (Fragments Per Kilobase of transcript per Million mapped reads) value > 1 in any sample were defined as being expressed.Differentially expressed genes were identified based on P = 0.05 and a two-fold difference cut-off using DEseq2 (Liu et al. 2021a).The heat map analysis of gene expression was performed with TBtools software (Wang et al. 2012).

qRT-PCR analysis
Total RNA was extracted from the AZs tissues at four phenological stages of young fruitlet (S1), swelling (S2), veraison (S3), and full maturity (S4).The qRT-PCR expression analysis of selected genes was conducted using the primers designed by Primer 3 (Table S4).The POLYUBIQUITIN 3 (UBQ3b) gene was used as an internal control since it was found to be one of the most stably expressed genes across multiple organs in rabbiteye blueberry (Vashisth et al. 2011).The qRT-PCR with SYBR Premix Ex Taq II (TaKaRa) was repeated at least three times, and the data were analyzed using the 2 −(ΔΔCt) method (Livak et al. 2001).

Observation of AZs structure
The AZs tissues were fixed in FAA (50% ethanol, 5% glacial acetic acid, 10% formalin, and 35% water by volume) for a minimum of 24 h and stored at 4 °C.They were then dehydrated with ethanol, embedded in molten paraplast, and sectioned into 10 μm thick ribbons using a Microm HM 355S microtome (ThermoFisher, Waltham, MA, USA).The slides were viewed using an Olympus BX60 microscope, and images were captured by an Olympus DP70 camera.

Cellulase activity, and other characteristics assays
Cellulase activity was measured according to the carboxymethyl cellulose (CMC) method (Ghose et al. 2013).The production rate of O 2 − was measured by monitoring the formation of nitrite (NO 2 − ) from hydroxylamine, and the absorbance was measured at 530 nm (Jambunathan et al. 2010).The determination of H 2 O 2 content was based on the absorbance of the titanium-peroxide complex at 415 nm as described by Qian et al. (2013).Relative conductivity and malondialdehyde (MDA) were measured as previously described (Li et al. 2016).

Identification and in silico analysis of VaGH9 genes in Vaccinium ashei
It has been documented that GH9 is involved in cell wall modification by controlling cellulose synthesis and hydrolysis, thereby playing the significant roles in various processes of plant growth and development (Bhandari et al. 2006;Campillo et al. 2012;Shani et al. 2006;Szyjanowicz et al. 2004;Zhou et al. 2006;Zuo et al. 2000).Thus, we implemented in silico analysis of GH9 gene family of Vaccinium ashei in terms of gene structures, physicochemical properties, structural motifs, chromosomal locations, and phylogenetic relationship before VaGH9s were functionally addressed in the AZs formation.A total of 61 VaGH9s were genome-widely identified using HMMER3.0method, and the redundant forms of the same gene were removed simultaneously (see "Materials and methods").GH9 family was divided into three distinct structural subclasses (A, B, and C) and were named VaGH9A1 to VaGH9C12 according to the standardized nomenclature of the GH9 family in Arabidopsis (Urbanowicz et al. 2007) (Table S1).
We analyzed and compared the physicochemical characteristics of VaGH9s proteins, such as amino acids (aa) composition, molecular weight (MW), isoelectric point (pI), hydrophobicity, and subcellular location (Table S1).Analysis of deduced protein sequence revealed that all 61 candidate VaGH9s proteins comprised the Glyco_hydro_9 domains (PF00759) with similar physicochemical properties.The predicted proteins encoded by VaGH9s had a range of length from 156 to 1368 aa, with a molecular mass ranging from 17.5 kDa (for VaGH9C9) to 152.6 kDa (for VaGH9A1).There were 35 acidic and 26 basic proteins, with the pI varying between 4.68 (for VaGH9C8) and 10.16 (for VaGH9B36).Most of the VaGH9 proteins were hydrophilic, with a grand average of hydropathicity (GRAVY) score less than 0. However, VaGH9A1 was a hydrophobic protein with a GRAVY score of 0.088.In silico analysis for subcellular localization predicted that 33 VaGH9s reside in both the cell membrane and cell wall, with 23 VaGH9s in the cell membrane and 5 VaGH9s in the cell wall.Taken together, the physicochemical properties of the proteins encoded by members of VaGH9 gene family appear to be distinct and divergent.

Phylogenetic analysis and classification of VaGH9 genes
To explore the evolutionary relationships between VaGH9s and GH9 proteins from other plants, we constructed a neighbor-joining phylogenetic tree using MEGA 7.0 (Kumar et al. 2016).The analysis revealed that VaGH9 proteins formed three Clades (A, B, and C) and clustered together with GH9 enzymes from other plant species (Fig. 1).Clade B, the largest subclass, consisted of 38 VaGH9 proteins (designated as VaGH9 B1-38); Clade C had 12 VaGH9 proteins (designated as VaGH9 C1-12); and Clade A, the smallest subclass, included 11 VaGH9 proteins (designated as VaGH9 A1-11).
It is notable that certain VaGH9s belong to the same clade as GH9s with known function in the phylogenetic tree.For example, VaGHA1-3 were grouped with AtGH9A1 and PtrKOR, both of which are membrane-bound endo-1,4-ßd-glucanases involved in the cellulose biosynthesis of both primary and secondary cell walls, and they play a role in cell wall loosening during root hair formation and growth (Boraston et al. 2004;Shani et al. 2006).The VaGH9s family may have undergone evolutionary diversification similar to that of its counterparts in four other plant species.The degree of sequence similarities between GH9 proteins in Vaccinium ashei and Populus L. (13.04%) were found to be higher compared to the similarity between Vaccinium ashei and Arabidopsis (12.83%),Vitis vinifera (12.24%), and Citrus sinensis (12.78%).This suggests that GH9 proteins in Vaccinium ashei are more closely related to Populus L. than to the other plants.

Gene structure analysis of VaGH9s in Vaccinium ashei
The VaGH9 gene family was divided into three subgroups according to the phylogenetic tree constructed from the amino acid sequences (Fig. 2A).To examine the organization of exons/introns, we compared the genomic and cDNA sequences of all VaGH9 members.The coding sequences of the entire VaGH9 family were interrupted by introns, with the number of exons ranging from 2 to 17 (Fig. 2B).Comparison of intron-exon structure of VaGH9A showed that most genes in Clade A have 5-7 exons, while VaGH9A1 contains 17 exons.Additionally, VaGH9A8-11 have one more exon than VaGH9A2-4.Clade B, represented by 38 members of VaGH9B group, is the largest and most diverse subclass in the VaGH9 family, and the number of introns in this group ranges from 1 to 13, with lengths varying from 52 to 9798 bp within the transcriptional regions (Fig. 2B).The majority of VaGH9s in Clade B have 5-8 exons, but there are some exceptions: VaGH9B4 and VaGH9B30 both contain 11 exons, while VaGH9B3, VaGH9B5, VaGH9B6, VaGH9B32, and VaGH9B33 all have 14 exons.The number of introns ranges from 1 (VaGH9B36) to 16 (VaGH9A1).
VaGH9C1-4 has one large intron (2965-3332 bp) and some small introns.Meanwhile, VaGH9C5-7 have seven small introns with length varying from 89 to 906 bp (Fig. 2B).Therefore, individual exons and introns vary considerably in number and size within and between the clades.In particular, 29 VaGH9s in clades A, B, and C have both 5′ and 3′ UTR, and the remaining VaGH9s lack either 5′-UTR or 3-UTR, suggesting that expression of VaGH9 genes is subject to at least three different post-transcriptional regulation processes.These findings imply that VaGH9s within the same clade having similar gene structures characterized by comparable numbers and phases of introns have evolved to the gene family expansion with limited functional diversification of the proteins, but with distinct regulation, stability, and localization of mRNA.

Conserved motifs and domains analyses of VaGH9s
The conserved motifs and domains were analyzed to further understand the functions of VaGH9s.In this study, we identified ten motifs and seven domains (Fig. 3A), and their annotations are provided in Table S2.Many closely related members in the phylogenetic tree were found to exhibit common motifs at the same alignment and position, suggesting that members of GH9 that are clustered in the same subgroup may have similar biological functions.A total of 44 VaGH9 proteins were found to have the motifs 1-10 in order, except for VaGH9A11,.The VaGH9A1-7 proteins within clade A lack only motif 5, while all members of clade C contain motif 5.All VaGH9B proteins in clade B have motifs 9-10, except for VaGH9B35 and VaGH9B36, which only have motifs 1 and 6.Particularly, VaGH9C12 carries motif 5 only.In general, the protein sequences of the members of the VaGH9 family are relatively conserved, and all the VaGH9 proteins possess Glyco_hydro_9 domain at comparable positions where all ten motifs are located (Fig. 3B).It is notable that VaGH9B30 protein has both PME and PMEI domains that are involved in pectin degradation (Wormit and Usadel 2018), whereas VaGH9B7 had PME domain only.The CBM49 domain that promotes the interaction  B30, B32-33, B35-36, C9-12) are very similar.This suggests the presence of limited functional diversity within a high degree of functional redundancy in the paralogous genes in the same clades.

Cis-acting regulatory elements in the promoter regions of VaGH9s
Genes with similar expression patterns are likely to have shared transcriptional regulation and common trans-acting binding sites.The cis-acting regulatory DNA sequences play a critical role in transcriptional control of temporal and spatial gene expression in response to developmental and environmental cues (Schmitz et al. 2022).In addition, expression divergence resulting from distinct cis-acting elements of the members of gene family may enhance understanding of the biological significance of different mechanisms of gene duplication.Thus, to further understand the regulatory mode of VaGH9 gene expression, we analyzed the promoter region of each VaGH9 for cis-acting regulatory elements (CAREs) by scanning a 2-kb upstream genomic sequence (Fig. 4A).We identified four categories of CAREs: stress response, plant growth and development, hormone regulation and transcription factor (Fig. 4B).Fifteen types of cis-elements were found to be associated with stress response, such as myeloblastosis (MYB) and myelocytometosis (MYC) elements, anaerobic induction element (ARE), low-temperature responsiveness (LTR), defense and stress responsiveness (TC-rich repeats) and wound-responsive element (WUN-motif) among others.The MYB is the most prevalent element in this category, accounting for 23.78%, followed by MYC at 20.48%.Notably, all of 61 VaGH9s possess MYB in common.In the category of plant growth and development, 16 types of ciselements were detected, including activating sequence-1 (as-1), meristem expression element (CAT-box), zein metabolism regulation element (O2-site), circadian control element, endosperm expression element (GCN4_motif) among others.The as-1 was the most common element in this category, accounting for 26.95%, followed by AAGAA-motif at 19.04%.In the hormone-responsive category, 16 types of cis-element were found, including ABRE related to ABA (abscisic acid) responsiveness (39.1%),CGTCA-motif and TGACG-motif responsive to jasmonic acid methyl ester (MeJA) (29.3%), four types of cis-elements related to auxin responsiveness (AuxRR-core, AuxRE, TGA-element, and TGA-box) (10.2%),P-box and GARE-motif related to gibberellin (GA) responsiveness (7.9%), TCA-element, TCA, SARE responsive to salicylic acid (SA)(7.8%),and ERE responsive to ethylene (ETH) (5.7%).In the transcription factor category, four types of cis-elements were determined: MYB-binding site responsive to drought inducibility (MBS) (54.35%),MYB-binding site involved in light responsiveness (MRE) (8.71%), and MYB-binding site involved in flavonoids biosynthetic genes regulation (MBS1) (6.51%).In addition, CCAAT-box was found to be the binding site of transcription factor MYBHv1, representing 30.43% of the transcription factor category (Fig. 4C).

Chromosomal localization and synteny analysis of VaGH9s
Gene duplication is a key mechanism that contributes to functional diversity.To identify the duplication events in VaGH9 genes, a collinearity analysis was performed using MCScanX software (Fig. 5).Among 61 VaGH9s, 50 were unevenly distributed in 34 Vaccinium ashei scaffold chromosomes, and the exact locations of the remaining 11 were not determined (Fig. 5).Twenty one chromosomes harbor 1 VaGH9, while two VaGH9s are separately distributed on each of 10 chromosomes and 3 VaGH9s in each of 3 chromosomes 12, 27 and 28.A total of 86 pairs of fragment duplication were found in 48 VaGH9 genes, with the highest number of duplications observed in chromosome 12 and no duplicated gene pairs in chromosomes 9 and 36.VaGH9B30 on chromosome 15 was found to have collinearity to VaGH9B7, VaGH9B25,VaGH9B26,VaGH9B27,VaGH9B28,VaGH9B29,and VaGH9B31.This expansion of VaGH9 gene family across the chromosomes is indicative of the occurrence of whole genome duplication (WGD) resulting in a balanced gene drive (Freeling et al. 2009).In addition, the number of gene pairs between Vaccinium ashei and other species was not correlated with the size of their genomes.The evolutionary rates and selective pressure of the VaGH9 gene family was estimated by calculating non-synonymous/synonymous substitution ratio (Ka/Ks) (Table S5).The results showed that the ratio is less than 1, with values ranging from 0 to 0.92, suggesting that VaGH9s in the Vaccinium ashei genome primarily evolved under the influence of a purifying selection process that removes deleterious mutations from a population (Del et al. 2021).

Transcriptomic profiles and differential expression of VaGH9s in fruit AZs at the version stage of Vaccinium ashei in response to ethylene and gibberellins
To identify the differentially expressed genes (DEGs) associated with AZs formation, we analyzed RNA-seq data from fruit AZs at the version stage in the absence (Ctrl) and presence of ethylene and gibberellins with three biological replicates.DEGs were determined by pairwise comparison of the expression data.As compared with Ctrl, ethephon treatment significantly up-regulated 2958 DEGs and down-regulated 2732 DEGs, while gibberellin treatment significantly up-regulated 121 DEGs and downregulated 277 DEGs.The Venn diagram showed that the number of genes identified in Ctrl, ethephon, and gibberellins were 398, 5690, and 2623, respectively (Fig. 7).Among them, 107, 4833, and 1838 DEGs were specifically expressed in Ctrl, ethephon, and gibberellin, respectively.To identify the function of sets of genes impacted, we conducted a GO-term enrichment analysis of DEGs in ethephon and gibberellins relative to Ctrl.DEGs genes were enriched in cell wall degradation, secondary metabolism, and growth-promoting hormone synthesis.
In total, 22 DEGs (marked with stars in Fig. 1) were identified in fruit AZs tissues.To identify group of DEGs that are co-regulated, their expression profiles were displayed as a heatmap.As shown in Fig. 8A, they can be categorized into five groups based on their expression patterns.The first group of VaGH9s (B11, B12, B13, C12) showed significant down-regulation after 24 h without treatment (C0 versus C24C), and they were not considerably affected by ethephon treatment (C24C versus C24E).However, their expression levels were significantly increased after gibberellin treatment (C24C versus C24G), although they remained much lower than those in control group (C0 versus C24G).The second group of VaGH9s (B9, B10) was strongly upregulated by gibberellin (C24C versus C24G) and downregulated by ethephon (C24C versus C24E).The third group of VGH9s (B14, B15, B16, B17, C5, C6, C7, C8, C10, C11) was strongly up-regulated by ethephon and either weakly up-regulated or unaffected by gibberellin.The fourth group of VaGH9s (A1, A2, A3, B1, B2) was fairly up-regulated by both ethephon and gibberellin.Finally, VaGH9B24 expression was not significantly affected by ethephon and gibberellin treatment but markedly reduced after 24 h.

VaGH9s expression profiles in AZs at four phenological stages
Eight VaGH9s (A1, A2, B2, B9, B12, B14, C7, C10) were chosen to investigate their qRT-PCR expression profiles in AZs tissues at four fruit developmental stages (Fig. 8B).As shown in Fig. 8C, the expression levels of VaGH9A1 and VaGH9C7 continually decreased, whereas those of VaGH9B2 and VaGH9C10 continually increased as the fruit development advanced.Interestingly, expression of VaGH9A2, VaGH9B12, and VaGH9B14 was culminated at the S2 stage (swelling stage), after which their expression gradually declined during fruit maturation.The expression of the VaGH9B9 also reached its peak at the S2 stage, after which their expression decreased in the veraison (S3) stage and increased during fruit ripening.VaGH9B24 was not expressed in S3 and S4 stages but expressed in S1 and S2 stages.

AZs structure and cellulase activity of Vaccinium ashei
The cellular morphology of AZs at which fruit detachment takes place was examined by scanning electron microscopy (Fig. 9).As the fruit matured, a group of small, closely interconnected cells arranged in a rectangular shape in the AZs appeared to transform into a more rounded and elongated form with cell walls being less rigid (Fig. 9B).It was observed that cell division in AZs occurred adjacent to the phloem in the longitudinal profile of ethephon-treated fruits, as indicated an arrow, while a cleft only formed in control and ethephon-treated fruits (Fig. 9C).Notably, an indentation between the pedicel and fruit was deepened in ethephon-treated fruits, whereas it remarkably diminished in gibberellins-treated fruits due to the expansion of the surrounding cells.
Cellulase plays a crucial role in abscission process in many fruit crops such as apples, grapes, orange, and blackberries (Dal et al. 2020;Li et al. 2020;Merelo et al. 2017;Zhang et al. 2019).It was reported that abscission was accompanied by burst of reactive oxygen species (ROS), and that ethephon induced oxidative stress in abscission zone (Goldental-Cohen et al. 2017;Yang et al. 2015).Thus, activities of cellulose, ROS and malondialdehyde (MDA) levels, and relative conductivity in the AZs at the version stage were measured every seven days over 28 days during the fruit's maturation.As shown in Fig. 10A, cellulase activities in all three samples rapidly increased as the fruit ripened, with the highest − and H 2 O 2 in the fruit AZs for 14-28 days.However, the application of gibberellins resulted in the suppression of the cellulase activity, O 2 − , H 2 O 2 , MDA, and relative conductivity in the AZs (Fig. 10A-E).

Discussion
It has been well documented that GH9s are involved in the breakdown of cell walls during plant growth and development including, but not limited to, fruit and leaf abscission, grain germination, pollen and seed dehiscence, cytokinesis, fruit softening, and senescence (Brummell et al. 1999;Campillo et al. 2012;Szyjanowicz et al. 2004;Xie et al. 2013;Yu et al. 2014;Zhou et al. 2006;Zuo et al. 2000).Due to their essential role in the organism's survival and fitness, GH9 gene family is ubiquitously present in plants.This study identified 61 VaGH9s genes that can be grouped  S1).Most VaGH9 proteins have a length that falls within the range of 417-694 aa, which is comparable to the range of GH9 proteins observed in other plants, such as 471-622 aa in Populus L., 364-675 aa in Arabidopsis, and 441-694 aa in Oryza sativa L. (Du et al. 2015;Nunan et al. 2001;Xie et al. 2013).Proteins encoded by genes with a broad range of sequence similarity in the same family often have similar secondary and tertiary structures (Hayashi et al. 2005).Therefore, the similarity in structure and length across different species suggests that those GH9 proteins likely share similar functions as a result of divergent evolution from a common ancestor.
Our analyses of gene structure, along with conserved protein motifs and domains, revealed that most VaGH9s encode the highly conserved amino acid sequences with analogous number and distribution of motifs and domains, despite considerable variation in the number and length of introns and UTR among the 61 VaGH9s (Figs. 2 and 3).Similar observations were previously made in poplar and rice (Du et al. 2015;Xie et al. 2013).The UTRs are known to play an important role in post-transcriptional regulation of gene expression, such as controlling mRNA stability, localization, and translation (Mayr et al. 2019;Sangar et al. 2007).The presence or absence of specific regulatory elements within introns can also influence gene expression patterns (Jia et al. 2020).Gene with more and longer introns may have a slower transcription rate since splicing process takes longer to complete.This may result in lower mRNA levels and potentially lower protein levels, as well as a longer halflife of the mRNA in certain tissues or developmental stages.Moreover, longer introns may also contain regulatory elements that can affect gene expression, such as enhancers and silencers (Swinburne and Silver 2008).These elements can affect the transcription rate and alter the spatiotemporal gene expression patterns.On the other hand, genes with fewer and shorter introns may have a faster transcription rate since the splicing process is simpler and faster.This can result in higher mRNA levels and potentially higher protein levels, as well as a shorter half-life of the mRNA in other tissues or developmental stages.However, shorter introns may also have fewer regulatory elements, which can limit the ability of the gene to be regulated in response to environmental and developmental cues.Accordingly, the highly divergent structures of exon-intron and UTR present in VaGH9s may give rise to functionally distinct paralogs that differ in gene expression level and/or spatiotemporal patterns.This notion is further substantiated by the presence of the four categories of CARES involved in stress response, growth and development, hormone regulation, and transcription factors in the promoter regions of VaGH9s.Among many CAREs, the promoters of all 61 VaGH9s were found to have MYB in common.The MYB transcription factors are known to be involved in almost all aspects of plant development and metabolism by controlling the morphology and patterns of the cells (Li et al. 2019b;Wang et al. 2021a).
Besides the common glyco-hydro_9 domain, six other discrete domains were found in VaGH9 proteins (Fig. 3B).In this regard, three groups of proteins attracted our attention.First, VaGH9A1, whose gene was shown to be expressed in fruit AZs and its expression decreased as fruit ripened (Fig. 8), has three unique features in that its gene has the largest number of exon and introns and encodes hydrophobic proteins harboring three distinct domains, ferric reductase transmembrane component-like, FAD-and NAD-binding domains.Ferric reductase is known to be an integral membrane protein containing intra-membrane binding site for heme and cytoplasmic-binding sites for nucleotide cofactors FAD and NAD, and this enzyme is essential for plants to acquire soluble ferrous iron (Fe 2+ ) from soils (Riethoven et al. 2010).Ferric reductase 3 gene was shown to be necessary for correct iron localization in both the root and shoot of Arabidopsis plants (Robinson et al. 1999).Hence, VaGH9A1 could also play a role in plant root and shoot growth and development.Secondly, VaGH9B30 has two separate PME and PMEI domains, whereas VaGHB7 has only PME domain involved in degradation of pectin that facilitates cell adhesion as the main component of the middle lamella of cell walls (Green et al. 2004).The VaGH9B7 and VaGHB30 genes, both of which lack ethylene-responsive CAREs (Fig. 4A), were not expressed in fruit AZs.Along this line, it was reported that the lack of a membrane-bound endo-1,4-β-glucanase may result in significant changes in pectin composition in the primary cell wall (Haas et al. 2021).The interplay between PME and PMEI was assumed to play a critical role in changing cell wall properties by regulating the extent of cell adhesion, cell wall porosity and elasticity in numerous biological processes (His et al. 2001).Thirdly, all VaGH9C5-8 proteins, which harbor non-catalytic C-terminal CBM49 domain, are very similar with respect to the size of proteins as well as the number and size of motifs and domains (Fig. 3).Except for VaGH9C8 gene, VaGH9C5-7 genes all have similar length and number of exon-introns and UTRs, as well as CARES (Figs. 2, 4).They were all expressed and strongly induced by ethephon in fruit AZs, and VaGH9C7 gene expression decreased as fruit ripened (Fig. 8).This result is consistent with the fact that CBMs are known to promote enzymatic hydrolysis of plant structural or storage polysaccharides by glycoside hydrolase through enhanced binding affinity (Wormit et al. 2018).
As for motif, it is notable that all the motifs are localized in the glycol_hydro_9 domain, irrespective of the clades.However, there are differences in the number and type of motifs between and within the clades (Fig. 3).For example, all VaGH9Cs in clade C harbor conserved motif 5, while  2).Due to the significant impact of motifs on protein function (Boraston et al. 2004), the presence or absence of a set of motifs within the glycoside hydrolase 9 domain of VaGH9s can lead to variations not only in the rate of catalytic activity and/or substrate affinity and specificity but also in the regulation of protein function.These unique sets of motifs may allow proteins to target specific substrate or respond to various signaling molecules.As previously reported (Cheng et al. 2019;Corbi-Verge and Kim 2016;Miller et al. 2018), some motifs could be subject to post-translational modifications or interact with regulatory proteins, leading to dissimilarities in the activation or inhibition of GH9s activity.Furthermore, different sets of motifs may also reflect changes in the protein's structural environment; for instance, if a GH is localized in a different subcellular compartment, it may require unique motifs to carry out its function in that environment.Overall, the presence of different sets of motifs in the catalytic domains in the members of GH9 family suggests that these proteins have diversified their functions over time, potentially to respond to dissimilar cellular contexts or to perform specialized tasks.Further studies to understand the differences in motif sets among GH9 family proteins would be useful for predicting their functional divergence and designing experiments to study their roles in biological processes.The chromosomal localization and synteny analysis showed that VaGH9s, which are unevenly distributed across the chromosomes of Vaccinium ashei, share conserved collinear blocks.This suggests that VaGH9 gene family has arisen from whole genome duplication event, which can lead to the presence of ortholog with similar functions between different genomes, and the paralogs within the same genome (Figs. 5, 6).The Ka/Ks values of VaGH9 gene pairs were found to be less than 1 (Table S5), which implies that the duplicated genes could be evolved under negative selection pressure to retain their original function as a result of purifying selection.Similar observations were previously made in the GH9 gene family of Populus L. and Zea mays (Buchanan et al. 2012;Du et al. 2015), as well as other gene families of COMT, ARF, and EXO70 in blueberry and grape (Li et al. 2022;Liu et al. 2021b;Wang et al. 2021b).
Considerable evidence suggests that GH9s have crucial functions in cell wall metabolism, particularly in relation to fruit ripening and abscission.For instance, the involvement of CEL1, the homolog of GH9B1, in regulating fruit maturation in strawberry has been demonstrated in several studies (Flors et al. 2007;Jara et al. 2019).Likewise, FaEG1, a member of GH9A, has been found to be involved in fruit ripening and organ abscission (Jara et al. 2019).In the present study, we identified 22 DEGs in fruit AZs tissues that could be assorted into five groups based on expression profiles in response to ethephon and gibberellins.Expression of VaGH9B24 lacking ETH and GA CAREs (Fig. 4A) was not influenced by ethephon and gibberellins treatment but rather strongly declined after 24 h, while the other 21 genes were either strongly or weakly up-and/or down-regulated by ethephon and gibberellins treatment (Fig. 8A).Thus, it appears that VaGH9B24 expression is entirely dependent on the progression of fruit development and its role is diminished during abscission process.Our qRT-PCR analysis of selected nine VaGH9s at four phenological stages of fruits (Fig. 8B) suggested their involvement in different stages of cell division process during abscission.Expression of VaGH9A1 and VaGH9C7 continuously declined, whereas that of VaGH9B2 and VaGH9C10 continuously increased in the fruit AZs during fruit maturation.This could be explained by the fact that VaGH9A1 and VaGH9C7 are involved in the pre-formation stage where AZs cell layers are differentiated, whereas VaGH9B2 and VaGH9C10 are involved in the separation stage where cell wall disassembly actively occurs.Along this line, our phylogenetic analysis (Fig. 1) showed that VaGH9A1 was grouped with AtGH9A1 and PtrKOR1 (Fig. 1), both of which were found to be involved in the cellulose biosynthesis of both primary and secondary cell walls (Bhandari et al. 2006;Shani et al. 2006).VaGH9B2 was grouped with AtGH9B1 (Arabidopsis cellulase 1) that is thought to be involved in cell wall loosening that drives anisotropic cell elongation (Lipchinsky et al. 2013).Both VaGH9C7 and VaGH9C10 were grouped with AtGH9C2 that was assumed to be engaged in modifying cell wall crystallinity (Glass et al. 2015).There appears to be no direct relationship between the number of motifs and the pattern of expression since VaGH9B2 and VaGH9C10 have 10 and 5 motifs, respectively (Fig. 3A).VaGH9B24, which has ABAresponsive element in its promoter region (Fig. 4A), was not expressed at S3 and S4 stages but expressed at S1 and S2 stages (Fig. 8C).VaBH9B24 belongs to the same phylogenetic group as AtGH9B8, which is assumed to have a cell wall-organizing function and interact with ABA signaling terminator (Turco et al. 2017;Wang et al. 2020).VaGH9A2, VaGH9B12, and VaGH9B14 were expressed at the highest level at the swelling (S2) stage, after which their expression gradually decreased during fruit maturation.VaGH9A2 is classified in the same phylogenetic group as VaGH9A1, while VaGH9B12 and VaGH9B14 are part of the group that includes AtGH9B5, which was believed to be involved in preserving the mechanical durability of the secondary cell wall (Glass et al. 2015).Thus, like VaGH9A1 and VaGH9C7, the four VaGH9s (B24, A2, B12, and B14) also likely take part in the pre-formation stage of abscission, after which their expression decreases as the fruit matures.
It has been postulated that AZs are generated by the differentiation of the cells at a specific site into functionally competent cells that are sensitive to abscission signals (Estornell et al. 2013;Ma et al. 2021).Thus, we microscopically examined the fruit AZs tissues of Vaccinium ashei (Fig. 9).The result showed that epidermal cells in AZs are elongated, and that cell separation occurred in the site of fruit-pedicel junction, which was in agreement with the research on citrus and litchi (Li et al. 2019a;Merelo et al. 2017).However, when the fruits were treated with gibberellins, fruit abscission was inhibited and cell expansion occurred (Fig. 9), as previously reported on strawberry by Gu et al. (2019).It was observed that abscission was accompanied by the increased cell wall metabolism-related activities such as cellulose, ROS, relative conductivity, and MDA (Fig. 10).Similar observations were previously noted in mango, sweet cherries, and peach (Brahem et al. 2017;Chen et al. 2017;Gu et al. 2019;Huan et al. 2016;Ji et al. 2020;Luo et al. 2020).

Conclusions
In this study, we conducted a genome-wide analysis to identifyVaGH9 genes in Vaccinium ashei and grouped them into three distinct clades based on sequence homology.The study has also revealed that most VaGH9s encode highly conserved protein sequences with analogous motifs and domains, despite considerable variation in the number and length of introns and UTRs among 61 VaGH9s.This divergent structure may give rise to functionally distinct paralogs that differ in gene expression level and spatiotemporal patterns.The presence of four categories of CAREs associated with stress response, growth and development, hormone regulation, and transcription factors in the promoter regions of VaGH9s suggests that these genes may play a crucial role not only in plant growth and development but also stress and defense response.Moreover, this study has identified 22 VaGH9s that could play a role in the different stages of fruit abscission.The levels of cellulase activity and oxidative stress indicators rapidly increased in the fruit AZs tissues from S3 to S4 stage, with further increments observed after ethephon treatment and decrements observed after gibberellins treatment.These findings have significant implications for the breeding and improvement of fruit quality in Vaccinium ashei and other plants.However, further studies are needed to elucidate the precise roles and functions of VaGH9s in different tissues and under different environmental conditions.

Fig. 1
Fig. 1 Phylogenetic trees of the GH9s proteins from rabbiteye blueberry, Arabidopsis, poplar, grape, and orange.The complete amino acid sequences of 61 Vaccinium ashei, 25 Arabidopsis, 25 Populus L., 24 Vitis vinifera, and 39 Citrus sinensis GH9 proteins were retrieved from NCBI (www.ncbi.nlm.nih.gov) database, aligned and compared with ClustalX, and the phylogenetic trees were constructed using MEGA7.0 according to the neighbor-joining method with 1000

Fig. 2 Fig. 3
Fig. 2 The polygenetic relationship and gene structure of the VaGH9s in Vaccinium ashei.A The phylogenetic tree was constructed using MEGA7.0 according to the NL method.B Exon-intron structures of the 61 putative VaGH9 genes.The exons are shown by blue rectangle,

Fig. 4
Fig.4Analysis of cis-regulatory element distribution in the VaGH9 promoters.A The cis-elements in the VaGH9s were marked by different colors.B The number of cis-elements numbers in four categories; dark blue histograms (type1), stress responses; blue histograms (type 2), plant hormone regulation; cyan histograms (type 3), growth and development; and orange histograms (type 4), transcription factor.C The proportions of distinct cis-elements from four categories ◂

Fig. 5 A
Fig. 5 A Circos plot showing distribution of VaGH9 genes among 34 scaffold chromosomes and synteny analysis in Vaccinium ashei genome.The annotations on the fragments represent different chromosomes, and the numbers in the outermost circle denote the positions on the corresponding chromosomes.The approximate position

Fig. 6
Fig. 6 Synteny analysis of GH9 genes between Vaccinium ashei and other four species.Gray lines in the background indicate the collinear blocks within Vaccinium ashei and other plant genomes, whereas the red lines highlight the syntenic GH9 gene pairs

Fig. 8
Fig. 8 Expression analysis of VaGH9 genes during fruit development.A Heatmap of differentially expressed VaGH9 genes in AZs of rabbiteye blueberries treated with ethephon and gibberellins at the veraison stage.The color scales with Z scores beside the heat map indicate gene expression levels.The low and high transcript abundance are indicated by blue and red color, respectively.C0, 0-h untreated control; C24C, 24-h untreated control; C24E, 24-h treated ethephon; C24G, 24-h treated gibberellins.B Photograph of fruits at four devel-

Fig. 9
Fig. 9 Anatomical and morphological AZs structure of Vaccinium ashei.A Fruit-spike and abscission zone tissue.The rectangular box refers to the fruit-pedicel junction.B Epidermal cells in S3 and S4.C Morphological structure of AZ cell layers in veraison (S3) and fully ripe (S4) stage, and those treated with ethephon and gibberellins.Phloem (Po), pith (Pt), pedicel (PL), and fruit (FT)

Fig. 10
Fig. 10 Cellulase activity, and the contents of ROS, MDA, and relative conductivity in fruit AZs of Vaccinium ashei.The data represent mean of three replicates with three biological repeats.Error bars indicate SE