Cynodon dactylon cultivar Yangjiang was used for genome sequencing and assembly in this study. The bermudagrass turf were grown in turfgrass plots of Yangzhou University (longitude and latitude: 32°35’N, 119°40′E; average annual temperatures: 22.4°C; average annual precipitation: 1, 106 mm; annual average sunshine hours: 1, 960 h; soil type: 80% river sand; and 20% peat soil) under routine management conditions (irrigation: keep the soil moist as required; fertilization: four times/year; and mowing: one times/month) for 3 years. Healthy leaves were randomly collected from the turf plots. Half of the leaf samples were frozen and used for de novo sequencing, whereas another half of fresh leaf samples were used for Bionano and Hi-C sequencing. Oryza sativa subspecies indica cultivar 93–11 was grown in growth chamber at 24°C under 16 h/8 h light/dark conditions.

The genome size of C. dactylon cultivar Yangjiang was estimated using flow cytometry as previously described (Zhang et al., 2020). Specifically, O. sativa cv. 93–11 with a genome size of 430 Mb was used as an internal standard. Young leaves of C. dactylon and O. sativa were homogenized on ice in Galbraith’s buffer (45 mM MgCl₂, 30 mM sodium citrate, 20 mM MOPS, and 0.1% (v/v) Triton X-100, pH 7.0) with 50 μg mL⁻¹ propidium iodide. After filtration with 40 μm nylon cell strainer (BD Biosciences, Franklin Lakes, United States), samples were analyzed on a FACSCanto^™ II flow cytometer (BD Biosciences). The flow cytometry data were analyzed using BD Spectrum Viewer.

Genomic DNA was isolated from the frozen leave samples using the DNeasy Plant Mini Kit (Qiagen, Hilden, Germany). The DNA quality and concentration were tested by 1% agarose gel electrophoresis and Qubit 2.0 Fluorometer (Life Technologies, Carlsbad, United States). Two paired-end libraries with short insert size of 270 bp and 500 bp were constructed using the NEBNext^® Ultra^™ DNA Library Prep Kit for Illumina^® (New England Biolabs, Ipswich, United States) and sequenced on the Illumina HiSeq X Ten platform (Illumina, San Diego, United States). The raw Illumina sequencing reads were processed with SOAPnuke v2.1.6¹ to remove adapters and low-quality reads (Chen et al., 2018). The obtained 161.1 Gb high-quality sequencing reads were used to generate a k-mer depth distribution curve adopting the Jellyfish v2.3.0.² The obtained peak k-mer number (k = 27) and corresponding peak depth were calculated by GenomeScope v2.0³ to estimate the genome size and heterozygosity (Marçais and Kingsford, 2011).

High-molecular weight (HMW) DNA fragments were separated from the extracted genomic DNA samples using BluePippin Size Selection System (Sage Science, Beverly, United States) through pulse-field gel electrophoresis and eight 20-kb sequencing libraries were constructed using SMRTbell Template Prep Kit (Pacific Biosciences, Menlo Park, United States) following the manufacturer’s instructions. The libraries (16 SMRT cells) were sequenced on the PacBio RSII platform (Pacific Biosciences). Contig sequences were assembled from the 151.99 Gb PacBio sequencing reads using Hifiasm v0.12⁴ and polished by Racon v1.4.3⁵ (Cheng et al., 2021). The Illumina sequencing reads were aligned to the assembled contigs using Bwa-mem v2.2.1⁶ and the draft assembly was corrected by the aligned short sequences using Pilon v1.24⁷ (Walker et al., 2014).

HMW DNA was extracted from the agarose-embedded cell nuclei fractions, which were isolated from fresh leaf samples, using the Bionano Prep^™ Plant DNA Isolation Kit (Bionano Genomics, San Diego, United States) following the manufacturer’s instructions. The DNA was digested by the single-stranded nicking endonuclease Nt.BspQI, fluorescently labeled, loaded into a Saphyr Chip^®, and imaged on a Saphyr Optical Genome Mapping Instrument (Bionano Genomics). The 395.4 Gb image data were filtered using a molecule length cutoff of 100 kb and a label number cutoff of 6, and assembled to 954 genome maps. To assist genome assembly, contigs obtained from the above-mentioned PacBio sequencing were transformed into in silico BspQI-digested reference genome maps and compared with the optical genome maps. The aligned and merged genome maps were further transformed into scaffold sequences using the Bionano Solve^™ v3.6.1.⁸

Fresh leaf samples were fixed in 1% formaldehyde to maintain the 3-D structure of genome. Genomic DNA was extracted and digested with restriction endonuclease MboI. The digested DNA fragments were biotin-labeled at the ends and ligated to each other randomly. The ligated DNA was sheared into 300–600 bp fragments, blunt-end repaired, and purified using streptavidin pull-down. The purified DNA was also sequenced on the Illumina HiSeq X Ten platform, which yielded 231.38 Gb of data with 771 million paired-end reads. The paired-end reads were mapped to the assembled scaffold sequences using Juicer v1.6⁹ to discriminate valid and invalid interaction pairs (Durand et al., 2016). The obtained 185 million valid interaction pairs (55.5 Gb data) were further used to adjust the relative locations of the scaffolds and cluster the scaffolds into pseudochromosomes using 3D-DNA¹⁰ (Kronenberg et al., 2021).

Repetitive sequences were annotated by combining the homology alignment and de novo prediction approaches (Zhang et al., 2021). For the homology alignment approach, the assembled genome sequence was blast searched against the RepBase repeat sequence collection¹¹ using RepeatMasker v4.0.9¹² (Tempel, 2012). For the de novo prediction approach, five softwares, including RepeatModeler,¹³ PILER,¹⁴ RepeatScout,¹⁵ LTR_Finder,¹⁶ and Tandem Repeats Finder,¹⁷ were used to find the possible repeat sequences (Price et al., 2005; Flynn et al., 2020). The identified repetitive sequences were manually checked and classified according to the nomenclature system of transposons. The insertion time of different families of long-terminal repeat retrotransposons (LTR-RTs) were calculated using the formula T = k/2r, where k is the divergence distance between the 5′ LTR and 3′ LTR of intact LTR-RTs and r is the base substitution rate (1.38 × 10⁻⁸ substitutions/site/year for grasses; Ma et al., 2004). The LTR Assembly Index (LAI) scores of assembled pseudo chromosomes and whole genome were calculated using LTR_retriever v2.9.0¹⁸ with default parameters (Ou et al., 2018). Putative centromeric repeat arrays were specifically identified using Tandem Repeats Finder with searching parameters “1 1 2 80 5 200 2000 -d –h” as previously described (VanBuren et al., 2020). The identified centromeric repeat array sequences were used to construct a maximum likelihood phylogenetic tree using MEGA v10.0.5 with a bootstrap of 1,000.

Protein-coding genes were identified by combining the homology alignment prediction, ab initio prediction, and transcriptome-assisted prediction approaches (Zhang et al., 2021). For the homology alignment approach, protein sequences of Arabidopsis thaliana and five grass species, including O. sativa, Brachypodium distachyon, Zea mays, Sorghum bicolor, and Oropetium thomaeum, were downloaded from the Phytozome database¹⁹ and blast searched against the assembled genome sequence to identify the homologous proteins, which were then aligned to the genome by GeneWise²⁰ to annotate gene structures (Birney et al., 2004). Ab initio gene prediction was conducted using five softwares, including Augustus v3.4.0,²¹ geneid v1.4.4,²² FgeneSH,²³ GlimmerHMM v3.0.4,²⁴ and Genscan²⁵ with default parameters (Yao et al., 2005; Nachtweide and Stanke, 2019). For transcriptome-assisted prediction, the PacBio single-molecule transcriptome sequencing data of mixed organ samples (Zhang et al., 2018a) were aligned to the assembled genome using GMAP²⁶ and the gene structures were modeled using PASA,²⁷ whereas Illumina transcriptome sequencing data of six different organs (Chen et al., 2021) were aligned to the genome using TopHat v2.1.1²⁸ and the gene structures were modeled using Cufflinks v2.2.1²⁹ (Wu and Watanabe, 2005; Ghosh and Chan, 2016). A non-redundant reference gene set was generated by merging the predicted genes using EVidenceModeler v1.1.1³⁰ (Haas et al., 2008). Functional annotations of the reference gene set were obtained through orthology assignment of the eggNOG v5.0 database³¹ using eggNOG-mapper v2³² (Cantalapiedra et al., 2021). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation of the reference gene set were obtained through BLAST searching against the GO database³³ and the KEGG pathway database,³⁴ respectively. KEGG enrichment analysis were performed using KEGG-Orthology Based Annotation System (KOBAS; Xie et al., 2011).³⁵ Transcription factors (TFs) were annotated using iTAK³⁶ incorporated with PlantTFDB database³⁷ (Zheng et al., 2016).

rRNA and tRNA genes were predicted using the programs Barrnap³⁸ and tRNAscan-SE-2.0,³⁹ respectively (Chan et al., 2021). miRNA, snoRNA, and snRNA genes were all identified by searching against the Rfam database via Infernal v1.1.4⁴⁰ with default parameters (Nawrocki and Eddy, 2013).

The completeness and accuracy of the assembled genome and predicted reference gene set were both assessed using the embryophyta_odb10 core gene collect (1,375 genes) of the Benchmarking Universal Single-Copy Orthologs (BUSCO) v5.2.2 database⁴¹ (Simão et al., 2015). The number of single-copy and duplicated genes with complete coverage, genes with fragment coverage, and missing genes were all counted.

The protein sequences of A. thaliana, O. sativa, B. distachyon, Z. mays, S. bicolor, O. thomaeum, Panicum hallii, Setaria viridis, Hordeum vulgare, and Tritcum urartu were downloaded from the Phytozome database. Orthologous gene families were clustered using OrthoFinder v2.5.4⁴² through all-against-all blast alignment of these protein sequences and predicted protein sequences of C. dactylon (Emms and Kelly, 2015). The identified 112 single-copy orthologous gene families were aligned using MUSCLE⁴³ and the alignments of each gene family were concatenated to a super-alignment matrix. A phylogenetic tree was then constructed using OrthoFinder with A. thaliana as the outgroup. PAML v4.9⁴⁴ was used to estimate the divergence time of C. dactylon using recorded divergence times of other 10 species in the TimeTree database⁴⁵ as calibrations (Yang, 2007).

Homologous pairs of C. dactylon proteins were identified using the all-to-all search in BLASTP v2.12.0⁴⁶ with an E-value cutoff of 10⁻⁵. Syntenic blocks with at least 50 collinear gene pairs were then identified using MCScanX⁴⁷ with default parameters (Wang et al., 2012). The same method was used to identify the collinear blocks between C. dactylon and O. thomaeum/B. distachyon. Synonymous substitutions per site values (Ks) of syntenic gene pairs were calculated using PAML v4.9 and the distribution of Ks values was plotted to infer the time for speciation or whole-genome duplication (WGD) events using the formula T = Ks/2λ, where Ks is peak Ks value and λ is the average substitution rate (6.5 × 10⁻⁹ substitutions/site/year for grasses; Gaut et al., 1996).

The Illumina transcriptome sequencing data of six different organs of bermudagrass cultivar Yangjiang were aligned to the assembled genome using HISAT v2.1.1⁴⁸ with default parameters (Kim et al., 2015). The numbers of mapped reads for each genes were converted to RPKM (reads per kilobase of transcript per million mapped fragments) values. The log₂ transformed RPKM values were applied to perform Hierarchical clustering using Pearson’s correlation distance (Chen et al., 2021). The significantly expressed genes were defined as RPKM value > 1, the organ-enhanced genes were defined as RPKM value is 5-fold above the average RPKM values of other organs, whereas organ-enriched genes were defined as RPKM value is 5-fold above the RPKM values of any other organs (Uhlén et al., 2016; Nautiyal et al., 2020).

To obtain ZMM (acronym for Zip1-4, Msh4-5, and Mer3) protein-coding genes in C. dactylon and other 10 plant species, ZMM genes from A. thaliana were used as baits to search against the assembled genome of C. dactylon and other plant species recorded in Phytozome database or Ensembl Plants database⁴⁹ using BLASTP v2.12.0 with an E-value cutoff of 10⁻⁵. The gene copy numbers and chromosome locations of different genes were manually summarized based on their identities. For PROG1, LA1, and TAC1 genes, PROG1, LA1, and TAC1 proteins from O. sativa were used as baits to search against the assembled genome of C. dactylon and other seven species of Oryza genus recorded in Ensembl Plants database as described above. The amino acid sequences of proteins encoded by each gene families were searched against the Pfam database⁵⁰ for domain comparisons (Mistry et al., 2021).

The C. dactylon cultivar Yangjiang was used for genome sequencing. As a national authorized C. dactylon cultivar., cultivar Yangjiang is a typical turf-type common bermudagrass and is widely used for turf planting in China (Zhang et al., 2018a; Supplementary Figure S1). Based on the K-mer genome survey result, the estimated genome size of C. dactylon cultivar Yangjiang is approximately 984 Mb, which is in line with the flow cytometry genome size estimation result of 1.02 Gb (Supplementary Figure S2). K-mer analysis also revealed that the genome of C. dactylon cultivar Yangjiang has a high heterozygosity (1.92%) with a repeat frequency of 56.91%.

To overcome the impact of high heterozygosity on the genome assembly, we adopted an integrated assembly strategy incorporating PacBio sequencing, Illumina sequencing, and Bionano and Hi-C techniques with the haplotype-resolving Hifiasm algorithm (Cheng et al., 2021; Supplementary Figure S3). Firstly, 151.99 Gb PacBio long reads (about 150× coverage of the genome) were de novo assembled into contigs, which were polished by 161.1 Gb Illumina paired-end reads (about 160× coverage of the genome; Supplementary Tables S1 and S2). Totally, 3,703 contigs with a N50 contig length of 2.65 Mb and a sum contig length of 1.295 Gb were obtained (Table 1). Secondly, 395.4 Gb Bionano optical maps (about 390× coverage of the genome) were used to integrate the contigs into scaffolds (Supplementary Figure S4; Supplementary Table S3). This procedure generated 241 scaffolds with a N50 scaffold length of 9.38 Mb and a sum scaffold length of 1.26 Gb (Table 1). Lastly, 231.3 Gb Hi-C data (24% useful information, about 55× coverage of the genome) were used to further cluster the scaffolds into pseudo chromosomes (Supplementary Figure S5; Supplementary Table S4). The finally obtained genome assembly (1.01 Gb) contained 36 chromosome-level superscaffolds, among which the longest and the shortest are 52.77 Mb and 14.32 Mb, respectively (Figure 1; Table 1). The assembly size was consistent with the estimated genome size. Furthermore, BUSCO analysis against the 1,375 Embryophyta gene sets indicated that 96.2% complete genes were successfully identified in the genome assembly, among which 88.1% were duplicated genes (Supplementary Table S5). These results collectively suggested that the assembled C. dactylon genome is high quality and complete.

A total of 76,879 protein-coding genes with an average gene length of 3,535 bp and an average transcript number per gene of 1.9 were successfully predicted from the assembled genome (Table 2). The predicted gene model was also evaluated by BUSCO analysis. The result indicated that 1,324 (96.3%) complete core Embryophyta genes were identified and the majority (1,272, 96.07%) was duplicated genes (Supplementary Table S5). Among the predicted 146,743 transcripts, 87.89% (128, 966) were annotated by various functional database (Supplementary Figure S6; Supplementary Table S6). Functional classification further indicated that signal transduction mechanism, post-translation modification/protein turnover/chaperones, and transcription are the top three categories containing the largest number of transcripts (Supplementary Figure S7). Specifically, 4,888 transcription factors (TFs) belonging to 65 classes were successfully identified. Compared with other grass species, gene numbers of HSF, WRKY, NF-X1, NF-YA, NF-YC, CPP, GARP-G2-like, and DDT TF families were greatly increased in common bermudagrass (Supplementary Table S7). In addition, 6,265 non-protein-coding genes were also identified, including 1349 rRNAs, 2047 tRNAs, 1025 miRNAs, 1441 snoRNAs, and 403 snRNAs (Table 2).

Orthologous clustering of protein-coding genes of C. dactylon with other ten plant species totally identified 32,695 orthologous gene families, including 7,792 commonly shared gene families and 3,173 bermudagrass-specific gene families consisting of 9,152 genes (Supplementary Tables S8 and S9). KOBAS enrichment analysis indicated that these bermudagrass-specific genes were enriched in glutathione metabolism, zeatin biosynthesis, ubiquitin mediated proteolysis, and other eight pathways (q value <0.05; Supplementary Table S10). In agreement with the BUSCO analysis result, orthologous gene clustering further revealed that as many as 91.2% (70117) of C. dactylon genes are members of 17,632 multiple-copy gene families, which is much higher than that of other ten species (Figure 2A; Supplementary Table S8). A phylogenetic tree was constructed based on the 112 shared single-copy orthologous genes (Figure 2A). The result indicated that O. thomaeum was closest to C. dactylon and the estimated divergence time of the two species was between 17.85 to 29.19 (midvalue of 23.52) million years ago (MYA). In line with phylogenic relationships, C. dactylon shared more orthologous gene families with members of the PACMAD (acronym for Panicoideae, Aristidoideae, Chloridoideae, Micrairoideae, Arundinoideae, and Danthonioideae) clade of grasses, including O. thomaeum, S. bicolor, and S. viridis, compared with O. sativa belonging to the BEP (acronym for Bambusoideae, Ehrhartoideae, and Pooideae) clade of grasses (Supplementary Figure S8).

A total of 381.3 Mb of repetitive sequences were also annotated in the assembled C. dactylon genome (Table 2 and Supplementary Table S11). The most abundant repetitive sequences are retrotransposons (70.95% of repetitive sequences and 26.9% of genome assembly), among which LTR-RTs and non-LTR-RT represent 84.73 and 15.27%, respectively. DNA transposons make up 11.29% of the repetitive sequences (4.28% of genome assembly), whereas tandem repeats and unclassified repeat sequences account for 2.79 and 3.75% of the assembled genome, respectively. Interestingly, the total repetitive sequence content and retrotransposon content in C. dactylon (37.91 and 26.9%, respectively) were similar to those of closely related species, including O. thomaeum (45 and 26%, respectively), P. hallii (36 and 23%, respectively), and S. viridis (46 and 29%, respectively), but much lower than those of distantly related species, including Z. mays (82 and 76%, respectively), T. urartu (81 and 72%, respectively), and H. vulgare (80 and 75%, respectively; Figure 2B). It is also noteworthy that genes are unevenly distributed in different chromosomes (39.87 to 104.80 Mb⁻¹ in density), whereas similar distributions of repetitive sequences were found on all chromosomes of C. dactylon (482.23 to 559.93 Mb⁻¹ in density; Supplementary Table S12). The annotated LTR-RTs were further used to calculate the LAI of the assembled genome. The total LAI score of C. dactylon genome is 13.63, implying that the current assembly of C. dactylon genome reached the reference genome level (Supplementary Table S13; Ou et al., 2018).

Intra-genomic syntenic analysis totally detected 845 syntenic blocks containing 84,649 pairs of homoeologous genes in the C. dactylon genome, whereas 643 syntenic blocks containing 52,590 pairs of homoeologous genes were found between C. dactylon and O. thomaeum through inter-genomics syntenic analysis (Figure 3A and Supplementary Figure S9). Interestingly, the syntenic depth ratios of C. dactylon versus O. thomaeum and C. dactylon itself were 4:1 and 4:4, respectively, implying that C. dactylon genome is composed of four haplotypes containing the same number of chromosomes (Supplementary Figure S9). To distinguish homoeologous chromosomes from the four haplotypes of C. dactylon, putative centromeric array tandem repeat sequences were identified from the 36 chromosomes and were used to construct a maximum likelihood phylogenetic tree as previously described (Supplementary Table S14; VanBuren et al., 2020). The result indicated that the 36 centromeric array sequences showed distinguishing polymorphisms and could be clustered in four clades (Supplementary Figure S10). Based on this classification result and chromosome length variance, the four haplotypes, which were named as A1, A2, B1 and B2, respectively, were successfully resolved in the C. dactylon genome (Supplementary Table S15). In addition, syntenic analysis also revealed that chromosome 2, 3, and 10 of O. thomaeum are split and merged into chromosome 2 and 7 in four haplotypes of C. dactylon, whereas other chromosomes all have one-to-one correspondence (Figure 3A and Supplementary Figure S9).

Calculation of Ks of homologous gene pairs in the inter-genomic and intra-genomic synteny blocks not only confirmed the phylogenic analysis result that C. dactylon and O. thomaeum diverged at approximately 21.54 MYA (Ks = 0.28), but also indicated that two rounds of WGD events occurred in the evolutionary history of C. dactylon (Supplementary Figure S11). Specifically, the first WGD event occurred at approximately 5.38 MYA (Ks = 0.07), whereas the second WGD event occurred lately at about 0.77 MYA (Ks = 0.01; Supplementary Figure S11). Interestingly, the two WGD time points equivalent exactly to the divergence time of haplotypes A1/A2 with haplotypes B1/B2 and haplotype A1 with haplotype A2 (the same as haplotype B1 with haplotype B2), respectively (Figure 3B).

In combination with the orthologous gene clustering result, syntenic analysis totally identified 20,849 alleles in C. dactylon (Supplementary Table S16). Among these alleles, 11,614 have four allelic copies in all haplotypes, 2,711 have three allelic copies in three of four haplotypes, and 6,524 have two allelic copies in two of four haplotypes (Figure 3C and Supplementary Figure S12; Supplementary Table S17). Meanwhile, 3559, 4954, 1716, and 3530 orphan genes that exist as single-copy genes were also identified from haplotype A1, A2, B1, and B2, respectively (Figure 3C and Supplementary Figure S12; Supplementary Table S17). KOBAS enrichment analyses indicated that alleles were enriched in valine/leucine/isoleucine degradation, proteasome, brassinosteroid biosynthesis, and other eight pathways, whereas orphan genes were enriched in plant-pathogen interaction, base excision repair, DNA replication, and other nine pathways (q value <0.05; Supplementary Table S18). The expression abundance of alleles and orphan genes were further analyzed using the organ-specific transcriptome sequencing data of C. dactylon cultivar Yangjiang (Chen et al., 2021; Supplementary Table S19). The result indicated that similar portions of alleles and orphan genes in the four haplotypes were significantly expressed in six organs of bermudagrass; however, gene numbers of alleles and orphan genes enhance- or enrich-expressed in different organs, especially the three types of stems, varied greatly in the four haplotypes (Supplementary Figure S13; Supplementary Tables S20 and S21). Accordingly, the 11,614 four-copy alleles of the four haplotypes showed similar total expression abundance in the six organs (Figure 3D).

The distribution of repeat sequences in C. dactylon was also analyzed at the haplotype level. Among the four haplotypes, haplotype A2 and B1 has the minimum and maximum number of RTs, respectively (Supplementary Figure S14). By contrast, maximum number of four types of DNA transposons, including Tcl/mariner, EnSpm/CACTA, hAT, and muDR, was observed in haplotype B2, while haplotype A1 has the fewest muDR- and Helitron-type of DNA transposons (Supplementary Figure S14). Notably, total sequence length of Ty3-Gypsy LTR-RTs in haplotype B1 was 2.2 Mb larger than that of haplotype B2, which contributed approximately 40% of size variance between the two haplotypes, whereas another type of LTR-RTs, Ty1-Copia, showed similar sequence length in the two haplotypes (Supplementary Figure S14; Supplementary Tables S22 and S23). Moreover, 5,066 intact LTR-RTs were further used to estimate the insertion time of different families of LTR-RTs in C. dactylon genome (Supplementary Table S24). The results indicated that four families of LTR-RTs, including Athila, SIRE, TAR, and Tork, inserted into the four haplotypes of C. dactylon genome at different time, whereas other nine families showed similar insertion time range in the four haplotypes (Figure 3E). Interestingly, among the 244 active LTR-RTs with an insertion time of zero, 153 were located in three chromosomes of haplotype B1, 47 were located in two chromosomes of haplotype A1, whereas only 29 and 16 were located in single chromosome of haplotype A2 and B1, respectively (Supplementary Table S24).

As a polyploid plant species with four sets of chromosomes, C. dactylon might develop a mechanism to control proper pairing and segregation of chromosomes during meiosis thus maintain its genome stability (Svačina et al., 2020). ZMM proteins, which stabilize the D-loop crossover intermediate of synapsis, are important homologous recombination regulators in all eukaryotes (Pyatnitskaya et al., 2019; Figure 4A). Previous studies have illustrated that gene copy number reduction of a ZMM protein, MSH4, could prevent meiotic crossovers between non-homologous chromosomes and stabilize the genome in allotetraploid Brassica napus (Gonzalo et al., 2019). Similar gene copy number reduction of MSH4 was also observed in other two polyploidy plants, allotetraploid Gossypium hirsutum and hexaploid Triticum aestivum (Figure 4B; Supplementary Table S25). However, MSH4 and other four ZMM genes, including ZYP1, MER3, SHOC1, and MSH5, all existed as four-copy alleles in C. dactylon genome. By contrast, two ZMM genes, PTD and HEL10, existed as two-copy alleles, and another ZMM gene, ZIP4, existed as single-copy orphan gene in haplotype B1 of C. dactylon genome (Figure 4B; Supplementary Tables S25 and S26). Syntenic analysis further indicated that different sizes of chromosomal fragments containing ZIP4 were lost in other three haplotypes (Figure 4C). These results collectively implied that C. dactylon also evolved a ZMM-dependent regulatory mechanism to maintain its genome stability as other polyploidy plants did; however, the key regulator might be ZIP4 rather than MSH4.

As a widely used turfgrass species with two types of specialized stems, stolon and rhizome, C. dactylon exhibits a prostrate plant architecture owing to increased tiller angles of the two specialized stems (Dong and de Kroon, 1994). Previous studies have successfully identified several tiller-angle-regulating genes, including PROG1, TAC1, and LA1, in rice and other plants (Li et al., 2007; Yu et al., 2007; Jin et al., 2008; Tan et al., 2008; Figure 5A). Eight PROG1-like genes, four TAC1-like genes, and two LA1-like genes were also identified in C. dactylon (Figure 5B). Similar to semi-prostrate and prostrate growing Oryza genus plants, the family/genome gene number ratio of two prostrate growth-promoting genes, PROG1-like and TAC1-like, were higher than that of erect-growth-promoting LA1-like gene in C. dactylon (Figure 5B). Syntenic and phylogenic analysis revealed that six of the eight PROG1-like genes existed as three-copy alleles and the remaining two genes existed as two-copy alleles, the four TAC1-like genes existed as four-copy alleles, whereas two LA1-like genes existed as two-copy alleles (Figure 5C and Supplementary Figures S15, S16). Pfam domain analysis further indicated that all eight PROG1-like proteins have the conserved C₂H₂-type zinc finger domain identified in the functional OsPROG1 and OgPROG7 proteins, however, both two LA1-like proteins of C. dactylon lack the functional C-terminal conserved region V (Figure 5D and Supplementary Figure S15). Moreover, the LA1-like protein encoded by the allele of haplotype A2 further lack the functional N-terminal conserved region I and two other conserved regions II and III (Yoshihara and Spalding, 2020; Figure 5D). In combination with the observation that large chromosome fragments containing the LA1-like gene locus were lost in other two haplotypes (Supplementary Figure S16), sequence variation of LA1-like genes in the two residual alleles suggested that LA1 protein activity was inhibited in C. dactylon. In addition, both two LA1-like genes were weakly expressed in stolon and rhizome, whereas three of four TAC1-like genes were preferentially expressed in the two specialized stems (Supplementary Figure S16; Supplementary Table S27). These results collectively implied that different tiller-angle-regulating genes were synergistically evolved to promote a prostrate plant architecture in C. dactylon.