A high-quality genome provides insights into the new taxonomic status and genomic characteristics of Cladopus chinensis (Podostemaceae)

Genome sequencing and quality assembly

We estimated a genome size of 820 Mb for C. chinensis, with heterozygosity of 0.73% and a repeat content of 51.0% (Fig. S1). A total of 95.36 Gb of PacBio long reads (~105X coverage of the genome) and 106.2 Gb of Illumina clean reads (~118X coverage of the genome) were generated, resulting in approximately 224.2-fold coverage of the C. chinensis genome (Table 1). The total size of all reads assembled from the C. chinensis genome was 827.92 Mb, consisting of 5629 contigs. The contig N50 was 1.42 Mb, and the longest contig was 8.89 Mb (Fig. 2 and Table 2). The genome size was close to the results based on flow cytometry (Fig. S2) and genome surveys (Fig. S1). The assembly quality was assessed through BUSCO analysis and the alignment of Illumina short reads to the genome. The NGS short reads from the Illumina sequencing platform were mapped to the contigs by using Bowtie2 software, and approximately 99.78% of the Illumina resequencing reads were mapped to the assembly (Table S1). BUSCO analysis revealed that the assembly completeness was 90.7% in the C. chinensis genome (Table S2). Using the contact matrix and the agglomerative hierarchical clustering method in LACHESIS¹⁸, the 538 contigs were successfully clustered into 29 superscaffolds (Fig. 3; Table S3). The scaffold N50 reached 1.42 and 21.22 Mb (Table 2), providing the first high-quality genome assembly for C. chinensis. Moreover, a total of 836 syntenic blocks were detected in the C. chinensis genome, which involved 15,410 genes. Synteny analysis of the 29 superscaffolds of C. chinensis confirmed that superscaffold no. 4 shared the most syntenic blocks with superscaffold no. 6 (Fig. 2).

Gene prediction and functional annotation

A combination of reference plant protein homology support, transcriptome data, and ab initio gene prediction was used to generate all gene models. All gene models were merged, and redundancy was removed with MAKER, leading to a total of 27,370 protein-coding genes. The average transcript length was 1465 bp, and the average CDS length was 868 bp (Table 3). We functionally annotated 22,564, 14,193, 22,564 and 11,950 genes to the eggNOG, GO, COG and KEGG databases, respectively, leading to 23,163 (84.63% of the total) genes showing at least one hit to the public databases (Table 3). In total, 2156 transcription factors were identified in the C. chinensis genome, and these genes were classified into 28 families, including 705 protein kinase family, 459 PPR, 155 MYB and 123 bHLH superfamily proteins (Table S4).

Annotation of noncoding RNAs (ncRNAs)

We identified snRNA, miRNA and rRNA genes in the C. chinensis genome from the Rfam database using BLASTN software (E-value ≤ 1e-5), and we used tRNAscan-SE and RNAmmer to predict tRNAs and rRNAs, resulting in a C. chinensis genome with 79 miRNAs, 1997 tRNAs, 397 rRNAs, 116 sRNAs, and 128 snRNAs (Table S5).

Repeat element (TE) annotation and burst analysis

The C. chinensis genome contained 490.13 Mb of repetitive sequences, accounting for 59.20% of the genome (Table S6). Long terminal repeat (LTR) retrotransposons accounted for 28.97% of the genome, with 22.30% Ty1/copia and 0.72% Ty3/gypsy sequences. Kimura distance analysis indicated an LTR burst (Fig. S3) involving the Ty1/copia and Ty3/gypsy superfamilies. Tandem Repeats Finder identified over 50,071 tandem repeats, accounting for 7.40% of the C. chinensis genome (Table S6). Telomeric and centromeric repeats were identified by searching the ends of contigs, and 8 telomeric repeat sequences and 572 centromeric repeat sequences were identified (Tables S7, S8).

Evolution of the C. chinensis genome

We collected the genome sequences of representative plant species and performed comparative genomic analysis with C. chinensis to reveal the genome evolution and divergence time of C. chinensis. The results suggested that gene family contractions outnumbered expansions in C. chinensis, R. communis, O. sativa, P. alba, B. oleracea, B. rapa, J. curcas, M. esculenta, A. thaliana, C. pepo and F. vesca, in contrast to the other seven species (Fig. 4a). The results support the view that Malpighiales and Rosales share a common Rosidae ancestor, which is the basal taxon of dicotyledons. The phylogenetic tree showed that C. chinensis phylogenetically diverged into an independent branch approximately 106 million years ago (Mya), without clustering into Malpighiales after the divergence of the A. thaliana (Brassicales) lineage 108 Mya. The data further confirmed that the divergence time of C. chinensis (106 Mya) was earlier than those of Malpighiales (82 Mya) and Rosales (102 Mya) and later than that of Brassicales (108 Mya). We also found that the internal branches separating Brassicales, Podestemales, Rosales, and Malpighiales were very short, as shown in Fig. 4a, suggesting the rapid diversification of these major clades. To obtain a highly reliable phylogenetic tree, we identified 91 single-copy homologous genes and performed coalescent analysis with RAxML and ASTRAL to estimate the species tree from gene trees. This result showed that C. chinensis phylogenetically clusters into Malpighiales, which is consistent with the results from NCBI Taxonomy (Fig. 4b). On the basis of the number of transversions at four-fold degenerate sites, we calculated an age distribution for all duplicate gene pairs. Using 8,074 paralogous gene pairs of similar ages and excluding tandem or local duplications, a large peak centered on a synonymous substitution rate (K_s) of approximately 0.14 was observed in the C. chinensis genome (Fig. 5), indicating only one recent whole-genome duplication (WGD) event and one ancient WGD that occurred 4.43 million years ago.

Gene family analysis

We performed comparative genomic analyses among 18 plant species and detected 81,585 families of homologous genes. A total of 17,964 gene families were identified in C. chinensis, among which 2164 and 14,356 gene families showed expansion and contraction, respectively (Fig. 4). A total of 2164 expanded gene families were annotated to KEGG pathways (Table S9) and GO terms (Table S10). KEGG analysis showed that most of the expanded genes were clustered in the categories of energy metabolism, signal transduction and aging. GO analysis showed that the expanded orthogroups were related to metabolic process and cellular component categories and the stimulus response, cysteine synthase, phototropism and sugar carrier terms. A total of 14,356 of the contracted gene families were involved in carbohydrate metabolism, signal transduction, secondary metabolism biosynthesis, amino acid metabolism, the nervous system, the metabolism of terpenoids and polyketides, the immune system, environmental adaptation and aging (Table S11). The GO terms of the contracted genes were related to the cellular component, flavonoid biosynthetic process, developmental process, stimulus response, catalytic activity, biological regulation, transcription regulator activity, immune system process, signaling, growth and reproduction process categories (Table S12). When we sought to investigate the environmental adaptation of C. chinensis, we found that a series of photosynthesis and energy metabolism-associated gene families exhibit significant expansion, including the NADPH oxidase (NDHF, NDHB1, GPSA), NADH dehydrogenase (NAD1, NDHD, NAD4L, NAD7, NDHG), chlorophyll a-b binding protein (LHCB) and cytochrome P450 (CYP86) families (Table S13). Notably, CYP450 genes are involved in the regulation of basic developmental processes such as cell differentiation and growth, indicating effects on the expression of developmental genes¹⁹. Additionally, approximately half of the encoded proteins catalyze specific biochemical reactions in various metabolic pathways. To obtain insights into the evolutionary relationships among C. chinensis CYP450 family proteins, we constructed a phylogenetic trees based on the sequences of full-length CYP450 proteins from 216 amino acid sequences of C. chinensis (121) and A. thaliana (95) using the neighbor-joining method. The numbers of most CYP450 types, including CYP71, CYP86 and CYP90, were higher than in A. thaliana indicating a close relationship of the development and metabolism of C. chinensis plants (Fig. S4).

The comparison of C. chinensis, A. thaliana, O. sativa, P. alba, M. esculenta, P. euphratica, J. curcas, R. communis, C. pepo and H. brasiliensis revealed that 4858 (31.37%) of the 15,485 C. chinensis gene families were shared by the other nine species, whereas 5,636 gene families were unique to C. chinensis (Fig. 6; Fig. S5). These 5,636 unique families have been specific to C. chinensis during its long history of evolution. Some of them may have been lost in other species, although we believe that some gene families originated de novo in Podostemaceae. Functional analysis performed via GO analysis revealed that these 5,636 unique families were enriched in suberin biosynthesis, cationic antimicrobial peptides, folate biosynthesis and ABC transporters (Table S14). KEGG analysis showed that most of the 5636 unique families were clustered in the signal transduction, lipid metabolism, energy metabolism and environmental adaptation categories (Table S15).

Phylogenetic analysis of the WRKY, NAC and resistance gene families

A phylogenetic tree was constructed based on the WRKY gene family members from A. thaliana (72), P. trichocarpa (102) and C. chinensis (68) to examine the evolutionary relationships among them. According to our results, a total of 68 WRKY genes were identified in the C. chinensis genome and classified into 7 groups: I, IIa, IIb, IIc, IId, IIe and III. The subfamily III group was one of the most ancient WRKY types, while subfamily I had fewer members, with no more than six members in each species. The identified WRKY genes of C. chinensis were only classified into the IIc subfamily (Fig. 7a).

To obtain insights into the evolutionary relationships among C. chinensis NAC family proteins, phylogenetic trees were constructed based on the sequences of full-length NAC-box proteins from 350 amino acid sequences of C. chinensis (81), P. trichocarpa (165) and A. thaliana (104) using the neighbor-joining method. The 81 NAC genes of C. chinensis were classified into 18 categories, including NAC1, NAP, NAM and TIP subfamily genes. The number of NAM and NAP subfamily genes from the NAC family was higher than in A. thaliana and P. trichocarpa. We speculate that the members of the NAM and NAP subfamilies could play roles in the development and defense responses of C. chinensis (Fig. 7b).

Furthermore, we analyzed the evolutionary relationships among resistance (R) gene proteins from C. chinensis (26) and A. thaliana (423) via phylogenetic tree analysis. The results showed that the R gene proteins could be classified into three groups, including TNL, CNL and RLP/RLK subfamily genes. Interestingly, the number of R gene proteins from A. thaliana (423) was 16.26-fold higher than that in C. chinensis (26), and the RLP/RLK subfamily genes of C. chinensis were identified and clustered. It is speculated that the presence of most R gene proteins reduces resistance and environmental adaptability, resulting in the survival of C. chinensis only in good-quality water (Fig. S6).

Positively selected genes in C. chinensis

The Ka/Ks ratios of the single-copy genes were evaluated, and positively selected genes were identified in C. chinensis. A total of 490 candidate genes in C. chinensis underwent positive selection (P < 0.05). Most of which were enriched in GO terms related to the biological regulation, cellular process, developmental process, growth, metabolic process, stimulus response, catalytic activity, binding, positive regulation of biological process and cell part categories (Table S16; Fig. S7). KEGG enrichment analyses showed that the positively selected genes were enriched in the categories of environmental adaptation, endocrine system, metabolism of terpenoids and polyketides, lipid metabolism, energy metabolism, carbohydrate metabolism, biosynthesis of other secondary metabolites, amino acid metabolism, translation, signal transduction and cell growth (Table S17; Fig. S8).

Comparative analysis of gene expression

To obtain an overview of the transcriptome profile, principal component analysis (PCA) was performed according to normalized log10 (FPKM + 1) values. The first principal component (PC1), which explained 99.66% of the total variance, showed clearly different gene expression profiles between the root group and the shoot group (Fig. S9). We identified the specific genes responsible for shoot development and defense responses in C. chinensis. In comparison with the root group, we identified 3376 DEGs (log2 FC ≥ 2 and p-value ≤ 0.05) in the shoot group, including 2002 downregulated and 1374 upregulated genes (Table S18). Some of these DEGs were associated with shoot development, including WUS, ASL and STM. Among these genes, CcWUS10, CcASL4, CcASL9, CcSTM3 and CcSTM4 were expressed at low levels in the shoots, whereas CcWUS1, CcASL1, CcASL2 and CcASL7 were highly expressed in the shoots (Fig. 8a, Tables S18–S21). Podostemaceae species have extraordinary body plans showing unusual shoot organogenesis involving the modification of the SAM, making C. chinensis an ideal system for studying the development of roots and shoots in plants. To identify the potential key genes or transcription factors involved in shoot formation in C. chinensis, we identified 609 genes that were highly expressed in shoots (log2 FC ≥ 1.5 and p-value ≤ 0.01) (Table S22; Fig. S10), and these highly expressed genes were significantly enriched in the metabolic process, stimulus response, developmental process, environmental adaptation, signal transduction, energy metabolism, carbohydrate metabolism and secondary metabolite synthesis were identified according to KEGG and GO enrichment analyses (Figs. S11, S12; Tables S23, S24). Additionally, we found that the 609 highly expressed genes in the shoots included hormone biosynthesis and signaling genes, including 11 auxin, 1 cytokinin, 2 abscisic acid, 5 gibberellin and 1 brassinosteroid gene. The development of the plant shoot is dependent on the SAM, MYB and bHLH transcription factors, which could regulate meristem initiation, meristem function and leaf patterning. Notably, 107 transcription factors were identified from the 609 highly expressed genes in the shoots, including 11 MYB, 11 bHLH and HD-ZIP family proteins (Table S25). In addition, we observed that the CcASL1, CcASL2, CcASL7, CcASL8 and CcSTM1 genes were involved in shoot development and were significantly expressed in specific shoot tissues (Fig. S13). Furthermore, we identified 132 and 146 MYB proteins in C. chinensis and A. thaliana. These MYB genes are divided into 13 subfamilies, and the numbers of genes in subfamilies I, III, and X are greater than those in the other subfamilies of A. thaliana (Fig. 8b). MYB transcription factors are widely involved in the regulation of various physiological responses and phenylpropanoid metabolic pathways, indicating a close relationship with the growth of C. chinensis on submerged rock surfaces in rapids and waterfalls during the rainy season.

The expression of WUS, STM and ASL is involved in the initiation and maintenance of shoot development in Tristichoideae and Podostemoideae, respectively⁸. To confirm the accuracy of the obtained unigene expression levels, nine target genes (CcWUS1, CcWUS10, CcASL1, CcASL2, CcASL4, CcASL7, CcASL9, CcSTM3, CcSTM4) associated with shoot development were analyzed in the present study (Fig. S13). The relative levels of the amplified mRNAs were evaluated according to the 2^−ΔΔCt method. Compared with the FPKM values, the relative expression levels of all nine unigenes determined by qRT-PCR were consistent with the RNA-seq data.

Genome sequencing and de novo assembly

The root, leaf and shoot samples of C. chinensis (sample number CBH03085, Fig. S16) used for genome sequencing and assembly were collected from the Tingjiang River in Changting, Fujian, China (E116^。25'10", N25^。51'03", Alt. 281.3 m)²⁸. The voucher specimens of C. chinensis were deposited in the Plant Herbaria of the College of Life Science, Fujian Normal University (collection number FNU0039809). To generate enough short and long reads for the genome assembly of C. chinensis, next-generation sequencing (NGS) on the Illumina HiSeq X Ten platform and third-generation sequencing (TGS) on the PacBio SEQUEL platform were applied for genome sequencing. Genomic DNA was produced via the CTAB method, and the integrity of the DNA was checked via agarose gel electrophoresis. The purity and concentration of the DNA were analyzed by using a NanoDrop 2000 spectrophotometer (Thermo Scientific, USA). For Illumina X Ten sequencing, we constructed a paired-end library with 250 base pair (bp) sequences using the method indicated by the manufacturer. As a result, 125.95 Gb (~140X coverage of the estimated genome size, Table 1) of short reads were generated from the Illumina platform, which were further cleaned by using Trimmomatic (version 0.36) with the default parameters, resulting in 106.2 Gb (~118×, Table 1) of cleaned data for the following analysis. For PacBio library construction and sequencing, SMRTbell libraries (~20 Kbp) were obtained according to the PacBio protocol. After the removal of adapters, 95.36 Gb of subreads (~105×, Table 1) were corrected, trimmed and assembled using CANU (version 1.6) with the parameters corOutCoverage = 80 and corMinCoverage = 0. To improve accuracy, the primary contigs were further filtered with the Pilon program²⁹ using 106.2 Gb (118×) of Illumina paired-end reads. We summed the statistics of the assemblies in Table 1. Genome completeness was assessed using the BUSCO databases³⁰ with embryophyta_odb10 models.

Estimation of genome size and heterozygosity

Weused NGS short reads and a K-mer-based method to estimate the genome size, heterozygosity and repeat content of C. chinensis. Approximately 106.2 Gb of reads (250 bp) were generated and used to determine the abundance of 21-K-mers in the Illumina data. The frequency of 21-K-mers was counted with Jellyfish software³¹.

Root, leaf and shoot samples of C. chinensis were prepared according to a previously reported protocol³² and stained with propidium iodide (50 mg/ml). Quantification was performed using flow cytometry in a BD FACSCalibur cytometer (Becton Dickinson, San Jose, CA), and the results were calculated as the ratio of the mean fluorescence of C. chinensis to that of Zea mays B73. The estimated genome size of C. chinensis was 835 ± 5.52 Mb.

High-quality assembly using Hi-C technology

Samples were examined for the integrity of nuclei by DAPI staining to ensure the quality of the Hi-C library. Samples with confirmed high-quality nuclei were subjected to the Hi-C procedure^33,34. Chromatin was digested with the restriction enzyme Mbo I or Hind III and ligated together in situ after biotinylation. DNA fragments were enriched via the interaction of biotin and blunt-end ligation and then subjected to HiSeq sequencing. From Hi-C library sequencing, approximately 96.4 Gb of data were generated (Table 1). The sequencing reads were mapped to the C. chinensis genome with Bowtie software. We applied an iterative alignment method to increase the read mapping ratio. We aligned the two read ends to the genome independently, and only the read pairs in which both ends were uniquely aligned to the reference genome were used for the detection and filtering of valid interaction products by using HiC-Pro (version 2.7.8)³⁵. The order and direction of scaffolds/contigs were clustered into superscaffolds by using LACHESIS³⁶ based on the relationships among valid reads.

RNA extraction and sequencing

Total RNA was prepared from C. chinensis roots and shoots (Fig. S17) using TRIzol reagent (Invitrogen, California, USA). A NanoDrop 2000 spectrophotometer (Waltham, MA, USA) and an Agilent 2100 Bioanalyzer (Agilent Technologies, USA) were applied to check RNA quality; the absorbance at 260 nm/280 nm was 1.8, and the RIN value was 9.1. Equal amounts of RNA from each tissue were used for cDNA library construction. Approximately 77.7 Gb of transcript data were produced for C. chinensis from the Illumina HiSeq X Ten sequencing platform and processed using Trimmomatic (version 0.36) with the default parameters. Reads originating from RNA-seq were aligned to the reference genome using HISAT2³⁷. FPKM values and read counts were estimated using Stringtie³⁸ and Ballgown³⁹. The differential expression of genes was analyzed using edgeR⁴⁰, for which the criteria were a log2-fold change (FC) ≥ 1 and a false discovery rate (FDR) ≤ 0.05.

Genome annotation

In the C. chinensis genome, de novo- and homology-based approaches were combined to search TEs and other repetitive sequences. We identified repeat sequences using Tandem Repeats Finder⁴¹, LTR_FINDER⁴², RepeatProteinMask and RepeatMasker⁴³ and used Tandem Repeats Finder to search for tandem repeats in the genome assembly with the following parameters: Mismatch = 7, Match = 2, Delta = 7, PI = 10, Minscore = 50, PM = 100, and MaxPerid = 2,000. Using LTR_FINDER (version 1.0.6), a de novo repeat library was built. Subsequently, we aligned the genome sequences to Repbase TE (version 3.2.9)⁴⁴ by using RepeatMasker for the searching of known repeat sequences and mapping onto the de novo repeat libraries to identify novel types of repeat sequences. tRNAscan-SE (version 2.0.3)⁴⁵ was used to detect reliable tRNA positions, and noncoding RNA (ncRNA) was predicted by searching the RFAM (version 12.0)⁴⁶ databases by using Infernal software (version 1.0)⁴⁷ with the default parameters. Centromere and telomere repeats were identified with Tandem Repeats Finder (version 4.07b)⁴¹. We transformed the resulting ‘dat_dir file’ into a GFF3 file, which was used to identify centromeric and telomeric repeats.

The C. chinensis genome assembly was annotated via the following approaches: homology-based, transcriptome-based, and ab initio annotation. Thirteen representative species were selected to perform homology annotation, including Manihot esculenta (M. esculenta)⁴⁸, Populus euphratica (P. euphratica)⁴⁹, Populus alba (P. alba)⁵⁰, Jatropha curcas (J. curcas)⁵¹, Ricinus communis (R. communis)⁵², Hevea brasiliensis (H. brasiliensis)⁵³, Arabidopsis thaliana (A. thaliana)⁵⁴, Cucurbita pepo (C. pepo)⁵⁵, Fragaria vesca (F. vesca)⁵⁶, Malus domestica (M. domestica)⁵⁷, Populus trichocarpa (P. trichocarpa)⁵², Pyrus x bretschneideri (P. bretschneideri)⁵⁸ and Rosa chinensis (R. chinensis)⁵⁹. The protein sequences of these species were aligned to C. chinensis genome sequences using TBLASTN software⁶⁰ with an E-value ≤ 1e-5. Genewise (version 2.2.0)⁶¹ was utilized to predict the exact gene structures based on all TBLASTN results. We used Cufflinks (version 2.2.1)⁶² to preliminarily identify gene structures. Augustus⁶³ was used for ab initio annotation with the repeat-masked genome sequences. We integrated all genes predicted from the three annotation procedures with MAKER software⁶⁴.

Functional annotation of the protein-coding genes was carried out by using BLASTP with an E-value ≤ 1e-5 in four integrated protein sequence databases: eggNOG, GO, COG and KEGG. We used InterProScan (version 4.8)⁶⁵ and HMMER (version 3.1)⁶⁶ to annotate protein domains by searching the INTERPRO (version 32.0)⁶⁷ and Pfam (version 27.0)⁶⁸ databases. Gene Ontology (GO) terms were produced from the InterPro or Pfam entry⁶⁹. The pathways were assigned through BLAST searches in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (release 53) with an E-value ≤ 1e-5⁷⁰. The functions of the genes were predicted and classified using the Clusters of Orthologous Groups (COG) database⁷¹.

Phylogenetic analysis and divergence time estimation

We investigated the relationships of C. chinensis with 17 other species and performed a phylogenetic analysis based on protein-coding genes from the C. chinensis genome and 17 other species. We extracted and downloaded the protein sequences of single-copy genes from 18 species from the NCBI database, including M. esculenta, P. euphratica, P. alba, J. curcas, R. communis, H. brasiliensis, A. thaliana, C. pepo, F. vesca, M. domestica, P. trichocarpa, P. bretschneideri, R. chinensis, Brassica oleracea (B. oleracea), Brassica rapa (B. rapa), Raphanus sativus (R. sativus) and Oryza sativa (O. sativa). The similarities among the proteins from all species were searched in an all-to-all manner by using BLASTP software with an E-value ≤ 1e-5. By using OrthoFinder software (version 2.27)⁷², we generated multiple sequence alignments for the protein sequences in each single-copy family with the default parameters as well as for phylogenetic tree construction. We designated O. sativa as the outgroup of the phylogenetic tree. The phylogenetic relationships were constructed through the superalignment of the coding DNA sequences (CDSs) using the maximum likelihood (ML) method, and the divergence time between species was estimated using the MCMCtree program of PAML (http://abacus.gene.ucl.ac.uk/software/paml.html). The TimeTree database (http://www.time.org/) was used for the divergence time recalibration of these plant species. The CDSs were aligned with the guidance of the protein alignments and then concatenated into the superalignment matrix of each family. We compared the cluster size differences between the ancestor and each species and analyzed the expansion and contraction of the gene families by using CAFE software (version 2.1)⁷³. For coalescent analysis, the concatenated alignment of single-copy homologous genes for 18 species was used as the input for ML inference with RAxML (version 7.2.8)⁷⁴. The alignment was partitioned by the genes, and ProtTest⁷⁵ was used to select the appropriate model of amino acid substitution for each partition. We used the -f a option of RAxML to generate 200 rapid bootstrap replicates, followed by a search for the best-scoring ML tree, and then used ASTRAL⁷⁶ to estimate the species tree from the gene trees.

Synteny analysis

The conserved paralogs of the protein sequences of C. chinensis were obtained with BLASTP (E-value ≤ 1E-5). By using MCScanX (http://chibba.pgml.uga.edu/mcscab2), we identified collinearity blocks in the genome. The Circos tool (http://www.circos.ca) was used to map gene density, GC content and repeat content as well as gene synteny on individual pseudochromosomes.

Whole-genome duplication

We took advantage of the high-quality genome of C. chinensis, analyzed WGD events and determined the source of the high number of genes in C. chinensis. First, the protein sequences from C. chinensis were searched to identify syntenic blocks by using BLASTP with an E-value ≤ 1e-5. We identified gene synteny and collinearity by using MCScanX software⁷⁷ and calculated the synonymous substitution rate (Ks) for syntenic gene pairs using KaKs_Calculator software⁷⁸ and the Nei-Gojobori method⁷⁹.

Quantitative real-time PCR (qRT-PCR) validation

Total RNA was extracted from the three replicates using the TransZol Up Plus RNA Kit (Transgen Biotech, Beijing). Primers were designed using Primer Premier 5.0, and the sequences are listed in Table S27. qRT-PCR was performed using the ABI 7300 Real-time PCR System (Framingham, MA, USA) with SYBR Green PCR Master Mix (TaKaRa) following procedures described previously⁸⁰. All PCR assays were performed in triplicate. The reference gene was glyceraldehyde-3-phosphate dehydrogenase (GAPDH). The relative expression level was quantified via the 2^−ΔΔCt method⁸¹.