Genetic analysis and gene mapping of a dwarf and liguleless mutation in barley

Baojian Guo, Jiang Qi, Dongfang Li, Hongwei Sun, Chao Lyu, Feifei Wang, Juan Zhu,Ganggang Guo, Rugen Xu,

a Jiangsu Key Laboratory of Crop Genomics and Molecular Breeding/Key Laboratory of Plant Functional Genomics of the Ministry of Education/Jiangsu Key Laboratory of Crop Genetics and Physiology, Jiangsu Co-Innovation Center for Modern Production Technology of Grain Crops, Yangzhou University, Yangzhou 225009, Jiangsu, China

b Joint International Research Laboratory of Agriculture and Agri-Product Safety, The Ministry of Education of China, Yangzhou University, Yangzhou 225009, Jiangsu, China

c Key Laboratory of Crop Germplasm Resources and Utilization (MOA), The National Key Facility for Crop Gene Resources and Genetic Improvement, Institute of Crop Sciences,Chinese Academy of Agricultural Sciences, Beijing 100081, China

Keywords:Barley EMS Liguleless ALOG Shoot architecture

A B S T R A C T

1. Introduction

Leaf development is a key factor determining crop growth and productivity and has been a major target of crop domestication and improvement [1]. In cereal, each leaf consists of a leaf sheath,a leaf blade,a lamina joint,and boundary organs that contain a pair of ligules and auricles [2,3]. The lamina joint is mechanical tissue that bends the leaf blade away from the leaf sheath to form the leaf angle, which can directly influence canopy structure and increase the leaf area index, even affecting yield [4,5]. The ligule is a thin and tongue-like white membrane that keeps the leaf sheath closed tightly against the auricles to prevent rainwater or pathogens from entering [1]. Boundary organ formation and determination are independently controlled aspects of leaf sheath and blade.Although most cereals contain ligules and auricles,some have been lost [2]. Understanding the genetic basis of leaf boundary organ formation may contribute not only to elucidating shoot architecture but to improving crop grain yield.

Development of leaf boundary organs has been characterized in maize, rice, wheat, and barley at the genetic and molecular levels using mutants [2,6-9]. In maize, loss of function of the recessivelg1gene prevented the proper formation of ligules and auricles[6,7]. Further analysis suggested thatLG1, encodingSQUAMOSApromoter-binding proteins (SBP), was similar to anAntirhinum majusgene. Mutations in theLGgenes homologous to maizeLG1,which results in a loss of the auricle, ligule, and lamina joint in cereals [2,8-10]. Maizelg2mutants showed an age-dependent phenotype differing from that oflg1. TheLG2gene encodes a bZIP(basic-leucine zipper) transcription factor involved in the establishment of the leaf blade-sheath boundary and is expressed in the meristem and developing-ligule region [7]. Double-mutant genetic analysis [11] suggested thatLG1andLG2interact and act in the same pathway.BOPgenes function in ligule development,tillering, and flower identity in cereals. The barley recessiveuniculme4(cul4) gene results in ligulelessness and low tillering.Positional cloning [12] suggested that theCul4gene, encoding a BROAD-COMPLEX, TRAMTRACK, BRIC-à-BRAC (BTB)-ankyrin protein, is related toArabidopsisBLADE-ONPETIOLE-Like Protein 1(BOP1) and BOP2. In rice, three homologs of theArabidopsis BOP1gene determine the leaf sheath:blade ratio, as well as ligule and auricle differentiation [3]. BarleyLaxatum-a(Lax-a), a homolog ofArabidopsis BOP1/2, controls internode length and homeotic changes of the barley inflorescence.The functions ofLax-aare similar to those of theBdUNICULME4(CUL4)andBdLAXATUM-A(LAXA)genes inBrachypodium distachyon[13]. A characteristic of the barleyELI-Amutant is dwarfing and ligulelessness, with weak culms that break at the nodes.Comparative genetic studies predicted thatELI-Aencodes an unannotated protein containing an RNaseH-like domain [14].

To date, only three mutants,liguleless(li),cul4, andeli-aaffecting the leaf blade-sheath boundary have been reported in barley[8,12,14]. In this study, an EMS-induced dwarf and liguleless (dl)mutant with dwarf and liguleless phenotype was characterized.We cloned thedlgene using a map-based approach. These results lay the foundation for analyzing the molecular mechanisms underlyingDLgene function and provide insight into the regulation of boundary organs in barley.

2. Materials and methods

2.1. Plant materials and growth conditions

A two-rowed malting barley cultivar Yangnongpi 5 (wild type,WT),which is an elite cultivar released by the barley breeding program in Jiangsu province, China, was treated with 30 mmol L-1doses of ethylmethane sulfonate (EMS) [15] and a dwarf and liguleless mutant was identified among 1070 M2plants. Thedlmutant was crossed with both Yangnongpi 5 and Bowman. The parents and the F1and F2progenies (sown in 2019 and 2020)and F2:3lines(sown in 2021)were planted in the field at Yangzhou University experimental farm. The phenotypes (mutant type or wild type) were recorded in plants of the segregating F2population,and the homozygosity of individual F2plants was determined by observing segregation patterns in the F2:3generation. Twelve seeds of each line were planted 10 cm apart with 20 cm between rows. All field trials for fertilization, irrigation, pest control, and weed management were conducted in a manner similar to that used in local fields [16].

2.2. Phenotypic evaluation and statistical analysis

The main culms were collected for measuring plant height,internode length, and spikelet length using rulers at maturity.The number of internodes, productive tiller number, and kernels per spike were also recorded. Seeds were harvested and dried to constant weight. Seed length and width and 1000-kernel weight were measured using a SC-G grain appearance quality image analysis system (Hangzhou WSeen Detection Technology Co., Ltd,Hangzhou, Zhejiang, China). Ten plants were harvested from each of three replications. Comparison of phenotype means of WT anddlmutants was performed using Student’st-test.

2.3. Histology observation

To elucidate the detailed function of thedlmutation in barley dwarfism,cytological investigations of the WT anddlmutant were performed.For paraffin section study,approximately 2 cm lengths of culm were collected and fixed in FAA solution(90:5:5 v/v/v,70%ethanol: 100% formaldehyde: 100% acetic acid). The samples were dehydrated in a graded ethanol series (70%, 85%, and 95%) and embedded in paraffin. Sections approximately 10 μm thick were cut with a Leica EM UC7(Leica Microsystems,Bensheim,Germany)and stained with 1%(w/v)cresyl blue, scanned with a Pannoramic MIDI slice scanner (3D HISTECH, Budapest, Hungary), and photographed. The cell size was measured with ImageJ version 1.32 software (http://rsb.info.nih.gov/ij/). Ten samples per genotype were analyzed.

2.4. Bulked-segregant analysis (BSA) allele frequency visualization

Thirty F2plants displaying the mutant or wild-type phenotype were collected for DNA extraction. The modified CTAB method[17] was used for genomic DNA extraction. An Illumina HiSeq 2500 sequencer (Illumina, San Diego, CA, USA) was used for 2×100 bp paired-end sequencing. Raw reads were trimmed with Cutadapt 1.9.1 [18] and aligned to the barley cv. Morex reference genome[19].SNP calling was performed with Samtools and Bcftool v1.3 [20], and a custom Perl script, using a threshold of Q40 and 15-fold read depth. Δ(SNP-index) was defined as the difference between the mutant or wild-type SNP indices [21]. The Δ(SNPindex) value was plotted along the barley physical map [19] with a 100 kb sliding window size and 10 kb shift.

2.5. InDel and dCAPS marker development

The genomic DNA sequence of Bowman was retrieved from the IPK database (https://galaxy-web.ipk-gatersleben.de/). Genomic variation information was obtained using the SnpHub database model for Bowman and Yangnongpi 5 (http://wheat.cau.edu.cn/Wheat_SnpHub_Portal/collaboration_GBJ_191007/) [22]. Primers were designed to span the InDel or SNP with flanking region length of 400 bp. Amplicons ranged from 80 to 300 bp. SNPs were converted to dCAPS markers for SNP validation [23]. Primer design was performed using Primer 3 [24].

2.6. Mapping of the putative dl mutation

The F2population and F2:3lines ofdland Bowman were used for gene mapping. InDel and dCAPS markers showing polymorphism betweendland Bowman were selected for linkage analysis based on dwarf and liguleless segregants from the F2population(Table S1). PCR was performed in a final reaction volume of 10 μL containing 1 μL genomic DNA (100 ng μL-1), 1 μL solution of forward and reverse primers (10 pmol L-1), 3 μL sterile ddH2O,and 5 μL 2× Taq Master Mix (Vazyme Biotech Co., Ltd. Nanjing,Jiangsu, China). PCR products were separated on 8% nondenaturing polyacrylamide gels and visualized by silver staining.

2.7. Sequence analysis of candidate gene

DNA sequences of the candidate gene were amplified from both Yangnongpi 5 anddlusing the primer pairs shown in Table S1.PCR products were separated by 1% agarose gel electrophoresis. DNA fragments were cut from the gel and purified with the GeneJET Gel Extraction kit (ThermoFisher Scientific, Waltham, MA, USA).The fragments were connected to the pEASY-T1 cloning vector and sequenced. Sequence analysis was performed using DNAMAN software version 10 (https://www.lynnon.com/dnaman.html).

Domain searches were performed using HMMSCAN (https://www.ebi.ac.uk/Tools/hmmer/search/hmmscan) [25] with the default cutoff parameters.Multiple sequence alignments were performed using the CLUSTAL X program [26]. The following ALOG proteins were used as queries: rice (Os02g0166800and

Os06g0672400), maize (GRMZM5G816289andGRMZM2G087267),Brachypodium distachyon(Bradi3g04960andBradi1g32445), andSorghum bicolor(Sobic.003G421300,Sobic.010G225100andSobic.004G052800)[27].A phylogenetic tree was constructed based on this alignment result using the neighbor-joining(NJ)method in MEGA version 6 [28] with the following parameters: Poisson correction, pairwise deletion, uniform rates, and bootstrap (1000 replicates).

2.8. Isolation of total RNA and quantitative real-time PCR

Total RNA of each sample was isolated using an RNA extraction kit (TRIzol reagent, Invitrogen, Waltham, MA, USA) and incubated with RNase-free DNase I (TaKaRa, Kyoto, Japan) to remove DNA contamination. Quality and yield of RNA were characterized by agarose gel electrophoresis and NanoDrop 1000 Spectrophotometer (ThermoFisher Scientific). First -strand cDNA was generated from 2 μg total RNA with M-MLV reverse transcriptase (TaKaRa)using random primers. Primers used for quantitative real-time PCR are listed in Table S1. The reaction was performed in a 20-μL volume containing 10 mmol L-1Tris-HCl (pH 8.5), 50 mmol L-1KCl, 2 mmol L-1MgCl2, 0.4 μL DMSO, 200 mmol L-1dNTPs,10 pmol specific PCR primers, 1 U Taq DNA polymerase, and 0.5 μL SYBR GREEN I fluorescence dye. Quantitative real-time PCR was performed with a ViiA 7 Real-Time PCR System (Applied Biosystems, Foster, USA). The running protocol was as follows:94°C for 3 min, followed by 40 cycles at 94°C for 30 s, 58°C for 30 s, 72°C for 30 s, and a final extension of 72°C for 5 min. The amplification ofHvActin(accession numberHORVU1Hr1G002840)was employed as an internal standard. All reactions were run in triplicate.Ctvalues were determined by ViiA 7 software with default settings (Applied Biosystems). The relative expression levels of target genes were determined with the 2-ΔΔCTmethod[29]. For each sample, PCR was performed with three biological replicates.

3. Results

3.1. Phenotype of the dl mutant

At seedling stage, erect leaves lacking ligules were observed indlmutants compared with WT plants(Fig.S1).Thedlmutant consistently produced stable phenotypes with severe dwarfing and without ligules, with leaf sheaths directly connected to blades without a lamina joint region (Fig. 1A-D). At maturity, WT plants averaged 78 cm in height anddl28 cm (Table 1). Each internode of the mutant but the fifth showed a relatively consistent reduction in length. Productive tiller and internode numbers showed no significant difference between WT anddlmutants. The spike length and kernel number per spike of the mutant were respectively 68.4% and 71.4% those of WT plants (Table 1). Grain length, grain width, and 1000-kernel weight were decreased significantly indlmutants (Fig. 1E,F; Table 1).

Cell width varied significantly in the same regions of the third internode (Fig. 2A-D). The smaller cells at the internode of thedlmutant appeared to contribute to the dwarf phenotype(Fig.2E,F).

3.2. Genetic analysis

To determine whether the dwarf and liguleless phenotype was controlled by a single gene, thedlmutant was crossed to Yangnongpi 5 and to Bowman. F1plants of both crosses were normal in plant height and ligule indicating recessiveness of the mutant gene. In the two F2populations (sown in 2019), the mutant type was easy to distinguish by the dwarf and liguleless phenotype. In thedl×Yangnongpi 5 cross,201 F2plants were similar to the wild type and 58 were dwarf and liguleless, indicating segregation at a single locus (χ2=0.33 ＜χ2(0.05,1)=3.84). The F2population fromdl×Bowman showed a similar segregation with 178 wild-type and 56 dwarf and liguleless plants, fitting a 3:1 ratio(χ2=0.71 ＜χ2(0.05,1)=3.84).

3.3. Preliminary mapping of the dl mutation

After read-depth and quality filtration,only 22,817 high quality SNPs remained in each plant for subsequent SNP-index and Δ(SNPindex) calculation. In the visualization of Δ(SNP index), one sharp peak symmetrically distributed with respect to thex-axis was observed on the long arm of chromosome 7H (600-650 Mb)(Fig. 3A). This region was accordingly assigned as the candidate gene interval.

3.4. Fine map-based cloning of dl

Based on the preliminary mapping ofdl, 30 InDel markers evenly distributed in the preliminary mapping interval were designed to detect polymorphism between Bowman anddl.Twelve were polymorphic between Bowman anddland allowed the genotyping of 56 F2plants homozygous(from the 2019 trial)for thed1mutant phenotype (Table S1). Thedlmutant was mapped to a 4.14-Mb interval between markers 7H-59 and 7H-50 (Fig. 3B).Three polymorphic markers were used to genotype the 177 F2plants withdlmutant phenotype (2020 trial) (Table S1) and thedllocus was further fine-mapped to an interval of 1.21 Mb flanked by markers 7H-22 and 7H-28 (Fig. 3C).

A set of 1668 F2:3homozygous dwarfed liguleless mutants identified in the 2021 trial were genotyped.Markers 7H-22 and 7H-28 as well as five newly developed polymorphic markers(7H-73,7H-90, 7H-79, 7H-88, and 7H-76) were used to narrow the genomic interval of thedllocus (Table S1). Fifty recombinants between 7H and 22 and 7H-28 and respectively, 12, 7 and 2 recombinants between 7H and 73/32/90 and 7H-88 were identified (Fig. 3D).Thedllocus was localized downstream of marker 7H-90 and upstream of marker 7H-88 (Fig. 3D).

3.5. Candidate gene analysis

Based on the reference genome of Morex V1 barley[19], markers 7H-90 and 7H-88 were separated by 56.58 kb. This genomic region contained one annotated gene,HORVU7Hr1G106960(Fig. 3E). To identify the genetic variation indl, primers were designed using the reference genome for Morex V1 to amplify exons,and the upstream promoter sequence of the candidate gene,1.6 kb in length. Comparison of the DNA sequences ofHORVU7Hr1G106960in the WT anddlrevealed a C-to-T singlenucleotide substitution at exon position 790,which is a functional mutation resulting in a proline-to-serine substitution at th amino acid 264. The functional mutation was observed in the conserved C-terminus peptide (PLSVF) (Fig. 4A,B).

HORVU7Hr1G106960was annotated as Protein LIGHTDEPENDENT SHORT HYPOCOTYLS 4 (LSH4) [19] belonging to theALOGfamily,which is a regulator of the development of reproductive and lateral organs. Sequence analysis revealed that theHORVU7Hr1G106960gene consisted of 2126 bp, including a coding sequence (CDS) of 810 bp, 727 bp of the 5′untranslated region(UTR), and 589 bp of the 3′UTR (Fig. 4A). The CDS encodes a 269 amino-acid protein with an ALOG domain (DUF640) according to prediction by HMMER (Fig. 4B).HORVU7Hr1G106960shared 73.74%, 54.78%, and 69.79% sequence similarity with riceOsG1L2,OsG1L1and sorghumAWN1/AWN1-10, respectively (Fig. 4B).

Fig.1. Morphological comparison of the wide type and dl mutant.(A)Whole-plant phenotype of wild type(left)and dl mutant(right)at the filling stage.Scale bar,10 cm.(B)Internode of wild type cv dl mutant. Scale bar, 1 cm. (C) Ligular regions of wild type and dl mutant (indicated by red arrows). (D) Lamina joint of wild type and dl mutant(indicated by red arrows). (E) Grain length of wild type and dl mutant. Scale bar, 1 cm. (F) Grain width of wild type and dl mutant. Scale bar, 1 cm.

The relative expression levels of theHORVU7Hr1G106960gene are characterized in the barley genome database (https://galaxyweb.ipk-gatersleben.de/).HORVU7Hr1G106960was highly expressed in 16 selected tissues not including young developing inflorescences and roots (28 days after pollination, DAP) and senescing leaves (56 DAP) (Fig. S2). Real-time PCR showed that theDLallelic gene showed no expression changes in root, leaf and ligule region (Fig. 5), suggesting that the single-nucleotide substitution in the coding sequence does not affectDLallelic gene expression in the selected tissues of both WT anddlmutants.These findings indicatedHORVU7Hr1G106960as a candidate gene for the dwarf and liguleless character in barley.

4. Discussion

4.1. A novel shoot-architecture locus, dl, resides on chromosome 7HL

Fig.2. Paraffin section of the third internode of the wild type(Yangnongpi 5)and dl mutant.(A and B)Longitudinal sections of wild type and dl mutant stems.(C and D)Cross sections of wild type and dl mutant stems. (E) Cell length of stem. (F) Cell width of stem. Scale bars, 500 μm. ** indicates significant difference at P ＜0.01.

Morphological features of the ligule, auricles, and lamina joint,as well as their presence or absence, are used as key identifying features of agricultural and horticultural crop cultivars. For example,theliguleless(li)mutant of barley lacks ligule and auricle structures on all leaves and is easily recognized at every stage of plant growth from the seedling onwards[30].Recessive mutations at theCul4gene resulted in reduced tiller number and led to a liguleless phenotype but with intact auricles [12]. Barleyeli-amutants produce shorter plants with fewer tillers and liguleless leaves with weak culms [14]. In the present study,dlshowed a dwarf and liguleless phenotype with normal productive tiller number,in contrast to other reported mutants.In the F2and F2:3population, individual plants with a liguleless phenotype were accompanied by a dwarfism phenotype, suggesting that thedlmutant affects many aspects of barley shoot architecture. The reduced spike length, kernel number of spikes, and seed size relative to Yangnongpi 5 invites further analysis.

Several loci controlling dwarf and/or liguleless have been mapped to chromosomes in barley. Thelilocus has been located on the long arm of chromosome 2 of barley by linkage to other morphological mutants and bulk segregant analysis[30,31].Previous mapping localizedCul4,controlling tiller and leaf pattern, to the distal end of chromosome arm 3HL [12].ELIGULUM-Aon 2HS regulates lateral branch and leaf development in barley [14]. We anticipate that the analysis ofdlwill provide a new genetic resource for improvement of plant architecture in barley.

4.2. HORVU7Hr1G106960 as a candidate gene for DL

Sequence analysis of the cloneddlrevealed that it is distinguished by one SNP from theDLgene derived from Yangnongpi 5, with one amino acid change at the protein level (Figs. 3E,4A).DLwas annotated asLIGHT-DEPENDENT SHORT HYPOCOTYLS,which belongs to theALOGfamily.DLshared high amino acid sequence identity with both riceOsG1L1/2and sorghumAWN1/-AWN1-10(Fig. 4B). Variable sites on the conserved C-terminus peptide were observed in rice, maize, sorghum, andBrachypodium distachyon.

There were 10 predicted ALOG proteins in barley, and none of them have been described(Fig.S2).Phylogenetic analysis of ALOG proteins was performed to elucidate their evolutionary relationships in rice,Arabidopsisand barley (Fig. S2). The barleyALOGgenes shared higher sequence identity with rice than withArabidopsis.DLis an ortholog of riceOsG1L2,which showed an expression profile similar to that ofTAWAWA1(TAW1) [32]. Loss of function of theosg1l2mutant produced shorter panicles similar to those of thetaw1mutant [32]. Thedlmutant displayed dramatic phenotypic changes not only in vegetative but in reproductive organs.

Fig.3. BSA-seq-based cloning of the dl gene.(A),Allele frequency analysis.(B and C)Physical positions of DNA markers(according to the Morex V1 gene model[19])used for rough mapping of dl using 56 and 177 F2 plants homozygous for the d1 mutant phenotype.(D)physical positions of DNA markers used for fine mapping of dl using 1668 F2:3 plants homozygous with dwarf accompanied liguleless phenotype for the d1 mutant phenotype.(E)annotated genes in the fine-mapped region of the Morex V1 gene model.Arrows indicate the transcriptional orientation of genes.

4.3. ALOG gene family involvement in plant architecture

The ALOG protein, a family of plant-specific transcription factors that regulate reproductive growth in angiosperms, emerged before the evolution of land plants and has shown functional conservation and diversification during the evolution of land plants[33]. The ALOG protein contains a DNA-binding domain and has weak transcriptional activity [34,35].LSH1is the first functionally characterized gene and is involved in phytochrome-dependent light signaling inArabidopsis[36]. Overexpression of theArabidopsis ALOGfamily genes,LSH4andLSH3, induced extra flower differentiation within a flower [37]. Mutation ofLSH4also modifies cell orientation in the rib zone in thereplumlessmutant background,although neitherlsh4norlsh3lsh4show obvious phenotypes [38].TERMINATING FLOWER(TMF) encodes a member of theALOGfamily. Mutation inTMFinduces simplification of primary inflorescences into single flowers [34].TAW1is the closest homolog ofLSH3/OBO1inArabidopsisand a unique regulator of meristem activity in rice, regulating inflorescence architecture via suppression of the phase change to spikelet meristem identity[35].Mutations ofLONG STERILE LEMMA1(G1) andTRIANGULAR HULL(TH1),as well as the tomato TMF family member (tfam) mutant, display abnormal spikelet organs and floral organs, respectively [39,40].

Marchantia LATERAL ORGAN SUPPRESSOR 1(MpLOS1) belonging to theALOGprotein family,regulates meristem maintenance and lateral organ development inMarchantia[41]. Sorghumawn1, which originated in a gene duplication on chromosome 10,is responsible for awn loss during sorghum domestication or improvement [42].In the present study, phenotypic variation in spikes was observed between WT and thedlmutant,implying that theDLgene not only is involved in regulating shoot architecture but also functions in barley inflorescence development.

Fig.4. Gene structure of DL and its homologs in plants.(A)Schematic of candidate gene,showing a nonsynonymous mutation(P264S)in the dl mutant.White and gray box represent untranslated region and coding sequence, respectively. The start codon (ATG) and the stop codon (TGA) are indicated. (B) Amino acid sequence alignment of DL,HORVU.MOREX.r2.6HG0469780 (Morex V2 gene model) and its homolog in Oryza sativa (Os02g0166800 and Os06g0672400), Zea mays (GRMZM5G816289 and GRMZM2G087267), Brachypodium distachyon (Bradi3g04960 and Bradi1g32445), and Sorghum bicolor (Sobic.003G421300, Sobic.010G225100 and Sobic.004G052800 [27].

Fig.5. Expression of the DL gene in the wild type and the dl mutant in root,leaf,and ligule-region tissues.Expression was determined by qRT-PCR using HvActin gene as endogenous control and normalization relative to Yangnongpi 5(Yangnongpi 5=1).

We propose thatHORVU7Hr1G106960is a strong candidate gene forDL, a gene affecting shoot architecture. Functional characterization of the gene may be achieved by targeted mutagenesis,for example by genome editing.HORVU7Hr1G106960appears to be a valuable gene for studying shoot architecture in barley.

CRediT authorship contribution statement

Baojian Guo:Investigation, Visualization, Data curation, Writing - original draft.Jiang Qi:Investigation, Visualization, Data curation,Writing-original draft.Dongfang Li:Investigation,Visualization, Visualization.Hongwei Sun:Investigation, Data curation.Chao Lyu:Data curation.Feifei Wang:Software, Validation.Juan Zhu:Validation.Ganggang Guo:Funding acquisition,Writing-review&editing.Rugen Xu:Conceptualization,Funding acquisition, Writing - review & editing.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

We thank Prof.Weilong Guo(China Agricultural University)for bioinformatic analyses in his laboratory. We thank Meixue Zhou for helpful comments and suggestions on the manuscript. This work was supported by the Open Project Program of Joint International Research Laboratory of Agriculture and Agri-Product Safety,the Ministry of Education of China, Yangzhou University (JILARKF202002),Natural Science Foundation of the Jiangsu Higher Education Institutions of China (19KJA560005), China Agriculture Research System of MOF and MARA (CARS-05), and a Project Funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions.

Appendix A. Supplementary data

Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2022.01.006.