Digitaria genome analyses indicate introgression may drive local adaptation and herbicide resistance

Domain: neuroscience
Published: 2026-02-12
Editorial outcome: major_revision

Abstract
Crabgrasses (Digitaria spp.) rank among the most invasive upland weeds globally, exhibiting remarkable ecological plasticity. To uncover the genomic basis underlying their weediness and environmental adaptability, we generate a telomere-to-telomere (T2T) reference genome of Digitaria sanguinalis, along with its diploid progenitor D. radicosa and tetraploid progenitor D. milanjiana. In addition, we re-sequence 579 accessions collected from fields of 24 grain-producing provinces in China. A phenotypic classification system for Digitaria species is developed by integrating genome-scale data with detailed morphological observations. Phylogenomic and comparative genomic analyses resolve the evolutionary history of Digitaria and reveal lineage-specific expansions in herbicide-associated gene families. Extensive sympatric introgression is observed within Digitaria, correlated with environmental selection. Through dose-response assay of 196 accessions, we identify 19 candidate genes associated with herbicide resistance, including DsSOH1. Haplotype analysis indicates that the resistant allele at DsSOH1 locus originated from D. ciliaris, suggesting a possible evolutionary contribution of introgression to the distribution of resistance-associated alleles. Together, our findings provide genomic insights into the evolutionary success of D. sanguinalis and offer resources for studying polyploid evolution and environmental adaptation as well as precision weed management in agroecosystems.

Introduction
Weeds are ecologically resilient organisms that thrive in agricultural ecosystems, exhibiting broad adaptability1–3. Digitaria sanguinalis (L.) Scop., commonly known as large crabgrass, is one of the world’s worst weeds. It is widely distributed from tropical to temperate regions and thrives in both cultivated and no-tillage farming systems4,5. At high densities, this species can cause yield losses above 90% and is frequently found in soybean, maize, and sorghum fields6,7. The genus Digitaria comprises over 220 species, whose identification is complicated by extensive phenotypic convergence8–10. This taxonomic ambiguity has led to inconsistent assessments of species distributions, ultimately hindering the development and implementation of precise weed management strategies.

With widespread and intensifying use of agrochemicals, herbicide resistance has emerged as a global threat to effective weed control, currently affecting 273 weed species worldwide11–13. Herbicide resistance typically arises through two primary mechanisms14. The first is target-site resistance (TSR), which results from mutations or amplification of herbicide target genes. The second is non-target-site resistance (NTSR), which is conferred by enhanced metabolic detoxification or sequestration. This process involves major detoxification enzyme families such as cytochrome P450s, Glutathione transferases (GSTs), Aldo-keto reductases (AKRs) and ATP-binding cassette transporters (ABC transporters)15. For D. sanguinalis, resistance has been documented to acetolactate synthase (ALS), photosystem II, Acetyl CoA carboxylase (ACCase), and 5-enolpyruvylshikimate-3-phosphate synthase (EPSPS) inhibitor herbicides16–19. The ALS-inhibitor herbicide nicosulfuron has been widely used for large crabgrass control since its registration in China in the 1990s18. In recent years, however, widespread resistance to nicosulfuron has emerged across numerous D. sanguinalis populations, involving both TSR and NTSR mechanisms20–22. This escalating resistance evolution is a growing concern, posing a significant threat to the effectiveness and sustainability of current chemical weed management strategies. However, the molecular and evolutionary basis of NTSR remains unclear due to the lack of high-quality genomic resources.

Polyploidization is widespread across the plant kingdom and is particularly prevalent among weedy species, where it often contributes to increased resistance to stress23–25. The presence of multiple homoeologous genomes within a single nucleus introduces regulatory complexity, including coordinated gene expression and extensive epigenomic remodeling26,27. In grasses, polyploidization has repeatedly facilitated the acquisition of stress resilience28,29. A well-known example is hexaploid wheat (Triticum aestivum), which exhibits increased fitness by combining root sodium retention, mediated by HKT1.5 from diploid Aegilops tauschii, with enhanced seed germination from tetraploid emmer wheat29. Patterns of gene retention and loss following polyploidization are also shaped by selection30. Disease-resistance genes have undergone significant contraction in barnyardgrass (Echinochloa crus-galli), in contrast to their expansion in wheat31, suggesting divergent evolutionary trajectories driven by ecological pressures32. D. sanguinalis is hexaploid (2n = 6× = 54). Therefore, investigating the evolutionary consequences of polyploidization is crucial for understanding the genetic basis of its ecological success and remarkable adaptability across agroecosystems.

Introgression is increasingly recognized as a powerful evolutionary force across diverse taxa, particularly in response to rapid environmental changes, including anthropogenic stressors such as herbicide exposure and pollution33. While natural interspecific hybridization has been documented in Digitaria, the functional implications of such hybridization events remain poorly understood34–36. In European aspen (Populus tremula), adaptive introgression has been shown to facilitate adaptation to high latitudes37. Another notable example is the Gulf killifish (Fundulus grandis), which acquired pollution tolerance through recent introgression of aryl hydrocarbon receptor loci from F. heteroclitus38,39. Similarly, adaptive introgression may serve as a key mechanism underpinning the wide geographic distribution and rapid evolutionary response of large crabgrass to increasing herbicide selection pressure.

In this study, to investigate the evolutionary trajectory and adaptive mechanisms of Digitaria, we assemble T2T genomes of hexaploid D. sanguinalis, and its tetraploid progenitor D. milanjiana and diploid progenitor D. radicosa. We also conduct large-scale genomic analyses of 579 Digitaria accessions and perform herbicide dose-response assays on a representative panel of 196 accessions. Our results reveal that introgression from closely related species may facilitate ecological adaptation in D. sanguinalis, including the acquisition of NTSR-associated haplotypes, which may have accelerated the recent evolution of herbicide resistance. These findings advance our understanding of Digitaria adaptive evolution, and will be helpful for development of precision weed management strategies.

Results
Genomic landscape of Digitaria reference genomes
A representative D. sanguinalis accession (#YJ2023) was collected from agricultural fields in Shandong Province, China, and its genome was sequenced, as well as its tetraploid progenitor, D. milanjiana (accession #DZ2), and diploid progenitor, D. radicosa (accession #YZGJ2). Cytological analysis confirmed its hexaploid status (2n = 6× = 54) (Supplementary Fig. 1), and k-mer analysis (k = 21, peak depth = 34) of short-read data estimated a genome size of 1.24 Gb (Supplementary Fig. 2), consistent with flow cytometry results (1.35 pg/1C) for D. sanguinalis (Supplementary Fig. 3). The heterozygosity was estimated at 0.05%.

For de novo assembly of D. sanguinalis, we employed complementary long-read technologies: 81× coverage PacBio HiFi reads (N50 = 15.9 kb) and 88× coverage Nanopore ultra-long reads (N50 = 100.2 kb). A chromosome-level assembly was generated using Hi-C data (79× coverage) to scaffold the initial contigs. In total, 419 contigs were anchored into 27 scaffolds, yielding a final assembly size of 1.35 Gb. To delineate subgenomes, we identified subgenome-specific k-mers and clustered homeolog-differentiating scaffolds, enabling consistent partitioning into three distinct subgenomes (Supplementary Fig. 4). Three subgenomes were designated as C, D, and E, based on the markedly low mapping rate (19.72%) and genome coverage (17.53%) with D. exilis, supporting its distinct genomic origin40. Following error correction and scaffold ordering, a chromosome-level high-quality assembly was generated with subgenome sizes of 453.06 Mb (CH, subgenome C in the hexaploid), 419.87 Mb (DH), and 474.80 Mb (EH), respectively. The assembly contains only two unresolved gaps (Supplementary Fig. 5). Assessment with BUSCO (v5.6.1, Poales lineage dataset) showed 99.30% completeness of conserved genes41. Genome annotation was performed using a combination of homology-based, transcript-based, and ab initio prediction approaches, and 114,996 gene models were identified after filtering 709.61 Mb (51.80%) of repetitive sequences. Conserved centromeric regions containing tandem repeats (detailed below) were identified on all 27 chromosomes. Telomeric arrays (TTTAGGG)n were resolved at 50 of 54 chromosomes, comprising 23 fully terminal chromosomes (telomeres at both ends) and 4 partially terminal chromosomes (single telomere detected) (Supplementary Data 1). rDNA sequences were identified on seven chromosomes (Supplementary Table 1). Through nation-wide sampling in China, we identified the diploid (2n = 2× = 18) D. radicosa, and tetraploid (2n = 4× = 36) D. milanjiana, as the progenitors of D. sanguinalis (details provided in the next section) (Fig. 1a). In parallel, two high-quality genomes for D. radicosa (489.39 Mb, CC genome) and tetraploid D. milanjiana (908.85 Mb, with subgenome sizes of 421.50 Mb (DT, subgenome D in the tetraploid) and 469.45 Mb (ET)) were also generated in this study (Supplementary Table 2; Supplementary Note 1).

To assess the assembly quality of three Digitaria reference genomes, we first mapped Illumina paired-end reads to their respective assemblies. Sequencing reads from diploid, tetraploid, and hexaploid Digitaria accessions were uniformly mapped to the D. sanguinalis reference genome, with high mapping rates of 99.51%, 99.48%, and 99.58%, respectively (Supplementary Table 3). RNA-seq reads showed normal alignment ratios to their respective genomes (94.12% for D. radicosa, 89.81% for D. milanjiana, and 93.32% for D. sanguinalis) (Supplementary Data 2). The genome assembly index LAI scores were 14.41, 17.30, and 15.62 for the three genomes, respectively, comparable to those of Arabidopsis thaliana (TAIR10) and Vitis vinifera42,43. We also estimated base-level accuracy and completeness of these assemblies and high assembly consensus quality values, 58.70 (99.74%), 54.51 (99.74%), and 50.08 (99.60%), were achieved for the three genomes, respectively (Supplementary Table 3). For continuity, we detected potential assembly gaps with low-confidence read supports using CRAQ, and high scores were also obtained for the three assemblies44 (Supplementary Tables 3 and 4). Taken together, these results suggest that the three assemblies are of high quality in terms of continuity, completeness, and accuracy.

Whole-genome alignments revealed extensive synteny between D. sanguinalis and its diploid/tetraploid progenitors, with conserved macro-collinearity across 95.0% of the genome (Fig. 1b; and Supplementary Fig. 5). The results also provide independent validation of our chromosome-scale assembly quality.

We calculated the synonymous substitutions per synonymous site (Ks) for orthologous gene pairs among the Digitaria species and other grass species to estimate their divergence times (Fig. 1c). The Ks for orthologs between D. radicosa and Oryza sativa was 0.57, corresponding to an estimated divergence time of ~45 million years ago (mya) (Fig. 1c), while Digitaria diverged from Eleusine ~35 mya, and subsequently from Setaria around 17 mya. This phylogenetic branching timeline was consistent with previous estimates45,46 (Supplementary Fig. 6). The two lineages with diagnostic spikelet arrangement phenotypes, D. exilis and D. sanguinalis, diverged at 12 mya. A large-scale amplification of Gypsy-type transposable elements was observed in the D. sanguinalis genome compared to D. exilis (Supplementary Fig. 7a). Further analysis of subgenome differentiation in D. sanguinalis resolved three ancestral lineages: CH diverged around 6.9 mya, followed by the bifurcation of DH and EH at ~6.4 mya. Using insertion polymorphisms of long terminal repeat retrotransposons, calibrated against synonymous mutation rates, we estimated that the tetraploidization event leading to D. milanjiana occurred around 0.9 mya, while hexaploidization in D. sanguinalis followed at ~0.4 mya. Interestingly, DH and EH in D. sanguinalis exhibit similar patterns to those in D. milanjiana (Supplementary Fig. 8). In addition, a maximum-likelihood (ML) phylogeny based on a concatenated matrix of 2030 single-copy orthologs among the 14 (sub)genomes was reconstructed, and a coalescent-based phylogenetic analysis was also performed, integrating individual gene trees (Fig. 1c; and Supplementary Fig. 9). These divergence estimates were concordant with our Ks-based molecular dating analyses, with robust statistical support at all key nodes.

Using domain-based gene family quantification across 18 genomes, we observed significant lineage-specific contraction in biotic stress-responsive gene families, highlighting differences between polyploid upland and paddy weeds as well as crops (Fig. 1d). Consistent with evolutionary patterns observed in other weeds31, Digitaria exhibited pronounced contraction of NB-ARC domain-containing genes (median 252 ± 48.7 copies vs. 447 ± 94.0 in crops; two-tailed t-test, p = 0.002) and D-mannose-binding lectin genes (61 ± 9.1 vs. 116 ± 22.8 in crops; two-tailed t-test, p = 0.001) compared to crops. The decay in defense-related genes likely reflected ecological trade-offs that favor growth in ruderal habitats, where biotic stress responses may be reduced in importance31,47. Interestingly, comparative analyses revealed a significant expansion of UDP-glucuronosyl/glucosyltransferase (UDP/GT) genes in Digitaria (162 ± 28.5) relative to Echinochloa species (64 ± 13.1; two-tailed t-test, p = 1.63e-5). Furthermore, gene families such as GAI-RGA-SCR (GRAS, 67.0 ± 5.6 vs. 43.5 ± 7.3; two-tailed t-test, p = 0.008) also exhibited expansion (Supplementary Data 3). These expansions may have provided a genomic foundation for adaptation under modern weed management systems. UDP/GTs were significantly enriched on Chr3 in Digitaria (Fisher’s test, p < 0.01) (Supplementary Fig. 10), which also harbors the highest transposon density across the genome (Supplementary Fig. 7b). Several flood-adapted gene families, such as the anaerobic-response Apetala2, showed marked expansion in Echinochloa (copy number: 182.0 ± 12.3 vs. 159.0 ± 9.9 in D. sanguinalis; two-tailed t-test, p = 0.008), likely reflecting adaptation to prolonged submergence (Supplementary Data 3). Compared with ancestral genomes, D. sanguinalis exhibited a lineage-specific amplification of far-red impaired response 1 DNA-binding genes (avg. 71.3 copies vs. 8.5 copies in D. exilis, 59.5 copies in D. milanjiana, and 54.0 copies in D. radicosa) likely associated with its mat-forming growth habit and the shaded microhabitats it occupies (Fig. 1d). Concurrently, D. exilis showed a distinctive expansion in Heat shock protein 70, consistent with its adaptation to extreme thermal environments in west African provenances (Fig. 1d).

Population structure and demography of Digitaria species
To investigate the genetic diversity and population structure of Digitaria, we re-sequenced the genomes of 579 Digitaria accessions collected over the past decade from a wide range of habitats across China (Fig. 2a; and Supplementary Fig. 11; Supplementary Data 4). Based on read mapping rates and genome coverage to the D. sanguinalis reference genome (Fig. 2b; and Supplementary Data 4), combined with morphological traits and genome size estimates (flow cytometry and k-mer analysis) (Supplementary Fig. 12; and Supplementary Data 4), the 579 accessions were classified into two major clades defined by distinct spikelet arrangements (Supplementary Note 2). The first clade, the ternate-type group, includes D. exilis, D. ischaemum, and D. violascens, species with primarily European distributions. The second clade, the binate-type group, comprises D. ciliaris, D. bicornis, and D. sanguinalis, which collectively represent the most globally invasive agricultural weeds48,49. On average, D. bicornis showed a 76.82% mapping rate with average coverage of 65.44%, 83.89%, and 85.12% across the CH, DH, and EH, respectively (Supplementary Fig. 13). D. ciliaris achieved a 94.63% mapping rate, with CH, DH, and EH coverage of 95.36%, 84.20%, and 85.94%, respectively (Supplementary Fig. 13). D. sanguinalis accessions showed the highest mapping rates, averaging 97.46%, with CH, DH, and EH coverage of 98.11%, 98.04%, and 97.03%. In particular, D. milanjiana and D. radicosa had differential subgenome coverage patterns (Fig. 2b). For example, D. radicosa exhibited a mapping rate of 95.05%, with subgenome compartment coverage of CH: 95.63%, DH: 2.07%, and EH: 1.80%, while D. milanjiana showed a mapping rate of 93.17%, with CH, DH, and EH coverage of 7.96%, 86.57%, and 87.75%, respectively (Supplementary Data 4). This asymmetric subgenome coverage provides genomic evidence supporting their ancestral contributions (i.e., tetraploid and diploid progenitor) to the allohexaploid genome of D. sanguinalis. In contrast, ternate-type Digitaria species consistently exhibited low mapping rates (~18.2%) and subgenome coverage (~13.6%), consistent with Ks-based molecular divergence estimates (Figs. 1c and 2b), indicating a more distant evolutionary relationship with the binate-type clade.

Our nationwide sampling effort covered 25 provinces across China, spanning from Hainan (18.3°N) in the tropical south to Heilongjiang (47.3°N) in the subarctic north (Fig. 2a). D. sanguinalis was the dominant species (78.90%) in the sampling population, with notable prevalence in major agricultural regions such as Henan, Shandong, and Shanxi (Supplementary Data 4). In contrast, the majority (68.75%) of D. violascens accessions were collected from northeastern provinces, Heilongjiang and Jilin. Regional differences in diversity were evident. Southern coastal provinces exhibited higher diversity, with Hainan showing the highest Simpson diversity index and Shannon index, 0.60 and 0.99, respectively, followed by Guangxi and Guangdong. In contrast, lower diversity levels were observed in Anhui and Jiangsu, with Simpson index values of 0.07 and 0.11, respectively (Supplementary Table 5).

We called single nucleotide polymorphisms (SNPs) across the CH, DH, and EH, yielding 10.76 million SNPs and 23.76 SNPs/kb for CH, 9.07 million and 21.60 SNPs/kb for DH, and 7.45 million and 15.68 SNPs/kb for EH. These SNPs were used for population structure and phylogenetic analyses of binate-type group (ternate-type accessions were excluded from subsequent analyses due to low mapping rates) (Supplementary Figs. 14 and 15). Within the retained binate-type specimens, each species formed a monophyletic clade in maximum-likelihood phylogenies (Fig. 2c). Meanwhile, we identified phylogenetically intermediate outliers that displayed three convergent signatures of potential hybridization: firstly, they occupied intermediate positions in species phylogeny between species clades (Fig. 2d); secondly, principal component analysis (PCA) positioned them centrally between primary species clusters (e.g., between D. ciliaris and D. sanguinalis) (Supplementary Fig. 15); thirdly, ancestry composition analysis revealed admixed genomic profiles (Supplementary Fig. 14). Together, these findings point to possible historical introgression among Digitaria species, although further analyses are needed to robustly resolve the extent and directionality of these events. Within D. sanguinalis, four distinct varieties and three admixed populations were identified, each showing marked biogeographic specialization across China (Fig. 2b, d; Supplementary Fig. 16). The phylogenetically basal var. glabra was predominantly distributed in northeastern regions and is characterized by the narrowest leaf blades, most compact tiller angles, and the heaviest 1000-seed weight among all varieties (Fig. 2d, e). In contrast, var. parvispicula is dominant in southern provinces and is distinguished by elongated, slender leaves and the lightest grain weight (Fig. 2d, e; and Supplementary Figs. 17 and 18). In the middle-lower Yellow River basin, the remaining two varieties exhibit north-south geographical separation relative to the river (Fig. 2b; and Supplementary Fig. 16). Although morphologically similar, var. pubescens, localized on the northern slopes, is distinguished by heavier seeds and the presence of golden bristles on the lemma surface (Fig. 2b, e; and Supplementary Figs. 17 and 18). Field surveys revealed that admixed populations were found in transitional ecotones between discrete varietal ranges (Fig. 2b; and Supplementary Fig. 16), further supporting their role as dynamic genetic intermediates. This classification into four distinct varieties lay the groundwork for exploring whether historical introgression has contributed to the observed genetic and phenotypic diversity.

By comparing the unfolded joint site frequency spectrum (SFS) of unlinked SNPs, polarized using known ancestral alleles, we inferred that archaic introgression from D. ciliaris into D. sanguinalis started at ~431,653 years before present (yr BP, 95% confidence interval (CI) 350,044−2,420,993 yr BP), following the speciation of D. sanguinalis, which itself was estimated to have occurred ~1.072 mya (95% CI: 712,489-2,921,965 yr BP) (Fig. 2f). Var. glabra diverged from the ancestral D. sanguinalis lineage at ~53,926 yr BP (95% CI: 37,101-1,423,489 yr BP). This was followed by the separation of the common ancestor of var. sanguinalis and var. parvispicula from var. pubescens at ~46,719 yr BP (95% CI: 18,199-1,156,488 yr BP). During this diversification process, introgression was detected between the ancestral var. pubescens and the common ancestor of the other two varieties. The best-fitting demographic model also predicted ongoing, bidirectional introgression between D. ciliaris and multiple D. sanguinalis varieties since the emergence of var. parvispicula around 39,561 yr BP (95% CI: 2,919-438,092 yr BP) (Fig. 2f; and Supplementary Fig. 19; Supplementary Table 6 and 7). Given the weed’s prevalence in upland agricultural systems, these patterns likely reflect human-mediated dispersal. Agricultural practices, especially mechanized harvesting and commercial seed exchange, may have accelerated the nationwide spread of D. sanguinalis and facilitated genetic exchange among historically isolated subspecies. To examine more recent demographic trends over the past 10,000 years, both folded and unfolded SFS were analyzed, and a sharp population bottleneck was detected in var. glabra around 30,000 years ago (Supplementary Fig. 20). For D. ciliaris, a severe and continuous decline in the effective population size began ~1000 years ago, following an earlier bottleneck. In contrast, after recovering from a bottleneck around 5000 years ago, the population sizes of var. sanguinalis and var. parvispicula had remained relatively stable over the past 2000 years.

Overall, nucleotide diversity (π) displayed an asymmetric distribution across the Digitaria genomes, with DH in both D. ciliaris and D. sanguinalis exhibiting the lowest levels of diversity (Supplementary Fig. 21). Except in var. pubescens, EH showed higher diversity than CH, suggesting that var. pubescens underwent distinct genomic alterations leading to elevated diversity. Moreover, the divergence between DH and EH was more pronounced in D. sanguinalis than in other Digitaria species (Supplementary Fig. 9), implying that D. sanguinalis experienced stronger post-origin introgression than other species, with a directional bias favoring CH and EH genomic compartments.

Comparative analyses of population structure across spatially clustered and temporally distinct populations of D. sanguinalis from 2013 to 2023 in the Huang-Huai-Hai agroecosystems revealed clear evidence of population genetic restructuring over time (Supplementary Fig. 22a; and Supplementary Table 8). Accessions collected in 2013 and 2015 were designated as historical groups, whereas geographically proximate accessions collected in 2023 were defined as modern groups. Modern groups exhibited markedly higher admixture diversity compared to historical collections (Supplementary Fig. 22b). While core geographic clusters have retained notable spatial structure, temporal analyses revealed genetic homogenization over time. For example, in GR7, the proportion of the k4 ancestral component declined sharply from 89.47% in 2013 to 31.95% in 2023. Concurrently, the admixture index increased from 0.29 to 0.55, reflecting a substantial rise in genomic intermixing (Supplementary Fig. 23; and Supplementary Table 8).

Sympatric introgression driving local environmental adaptation
Introgression is known to have played a primary role in increasing genetic diversity. We further constructed a ML phylogeny using SNPs from chloroplast genomes to infer the maternal origins of Digitaria topologies (Fig. 3a). Species-level divergence patterns in the tree were congruent with those inferred from the nuclear SNPs, notably with D. violascens and D. ciliaris forming a monophyletic clade (Figs. 2b and 3a). Interestingly, a subset of D. sanguinalis accessions (n=9) was scattered within the clade containing D. ischaemum and D. violascens, suggesting persistent chloroplast introgression between lineages that diverged more than 10 mya (Fig. 3a, b). Although D. bicornis formed a monophyletic group in the nuclear phylogeny, its chloroplast genomes clustered within Clade 4, which is dominated by D. sanguinalis accessions (95.68%) (Figs. 2d and 3b). This implied a D. sanguinalis origin of the D. bicornis chloroplast genome, likely resulting from historical chloroplast capture through hybridization or horizontal transfer. The chloroplast phylogeny further supports the hypothesis of at least two distinct maternal donors in the origin of D. sanguinalis. One lineage corresponded to the monophyletic Clade 4, while the other branches as a sister group to the D. ciliaris monophyly, collectively forming Clade 3, a topology that was discordant with the nuclear genome (Fig. 3a, b). Both diploid and tetraploid ancestors of D. sanguinalis were positioned within Clade 3, with D. radicosa nested inside the D. ciliaris subclade (Fig. 3a). This placement aligned with the high CH coverage observed in D. ciliaris (95.36%), suggesting shared ancestry of CH between D. ciliaris and D. sanguinalis (Fig. 3a; and Supplementary Fig. 13). Overall, the chloroplast phylogenies indicate extensive interspecific introgression across Digitaria, most likely resulting from historical hybridization events50,51.

Given sympatric distribution observed for D. ciliaris and D. sanguinalis, we applied the genome-wide D-statistic and fd value to test for introgression between 7 sympatric population pairs of D. ciliaris and D. sanguinalis (Supplementary Fig. 24a, b). All 7 populations exhibited significant signals of introgression (D > 0.11 and Z > 3 across all 7 populations; Supplementary Table 9). The strongest signal was detected in the SX population (D = 0.16), whereas the JN exhibited the weakest evidence of introgression (D = 0.11) (Supplementary Table 9). This pattern was further supported by the same analysis using sliding windows, revealing that SX harbored the greatest proportion of introgression across the genome (PIG) across the genome (Fig. 3c; and Supplementary Table 10). In contrast, eastern populations, including DZ, JN, and LF, displayed fewer introgressed regions, indicating minimal impact of introgression in these lineages (PIG in Fig. 3c).

To assess the potential role of environmental factors in shaping sympatric introgression, we conducted matrix regression analyses (Mantel/partial Mantel tests) on pairwise shared introgressed genomic regions (PSIG), defined as the Jaccard similarity index of introgressed haplotypes between “D. sanguinalis - D. ciliaris” population pairs. Environmental and geographic distance matrices were constructed using data from 7 population pairs (Supplementary Data 5). Both Mantel tests revealed significant negative correlations between PSIG and geographic distance (Pearson’s r = −0.70, p = 0.02), as well as between PSIG and environmental divergence (Pearson’s r = −0.64, p = 0.01) (Supplementary Table 11). However, partial Mantel tests did not detect significant associations when controlling for the alternative variable, likely due to high collinearity between geographic and environmental distances (Pearson’s r = 0.61, p = 1.48e-5) (Supplementary Table 11). Furthermore, correlation analysis between PSIG and specific environmental variables identified a significant association with precipitation (Supplementary Fig. 24d; and Supplementary Table 11).

To investigate the functional relevance of introgressed genomic regions, we performed Pfam domain enrichment analyses across populations (Fig. 3d). The number of introgressed genes varied from 2181 to 3472 among populations, and enrichment profiles were largely population-specific, with 13–25 significantly enriched domains per population and only 14 domains shared across more than three populations (Supplementary Fig. 25). Several functionally relevant domains were recurrently enriched across populations, spanning categories such as detoxification (e.g., ABC transporters, UDP/GTs), signaling (e.g., protein kinases, WD40 repeats), and transcriptional regulation (e.g., Myb-like DNA-binding domains) (Fig. 3d; Supplementary Data 6). The Hsp20 domain, associated with heat stress responses, was uniquely enriched in the SX population, which experiences the highest temperatures among all sampled regions (Fig. 3d; Supplementary Fig. 24a; Supplementary Data 5), suggesting a possible case of environmentally mediated introgression linked to thermal adaptation.

To identify genetic variants associated with environmental adaptation, we conducted genotype-environment association (GEA) analyses using the latent factor mixed model (LFMM), which accounts for background population structure while testing for associations between genotypes and environmental variables52 (Fig. 3e; Supplementary Figs. 26-28). The analysis included 19 environmental variables, comprising 10 temperature-related and 9 precipitation-related factors (Supplementary Data 7). In total, we identified 9437 SNPs significantly associated with one or more environmental variables, corresponding to 2117 genes (Supplementary Figs. 26-28; Supplementary Data 8). These environment-associated variants exhibited a genome-wide distribution, without clustering in specific genomic regions, suggesting that adaptation in Digitaria is governed by polygenic architectures rather than by hard selective sweeps at a handful of loci (Supplementary Figs. 26-28).

Several genes previously implicated in climate adaptation were identified in our analysis as harboring variants significantly associated with environmental variables53 (Supplementary Figs. 26-28; Supplementary Data 8). For example, the gene DsRZ2, which is homologous to OsRZ2, exhibited strong associations with the minimum temperature of the coldest month (Bio6) (Fig. 3e). This gene encodes a protein containing both a zinc knuckle domain (PF00098) and an RNA recognition motif (PF00076), and is known to play a critical role in plant protection against cold and freezing stress54,55. Expression analysis, based on RNA-seq data from root, stem, old leaf, and young leaf tissues, confirmed ubiquitous expression of DsRZ2 with RPKM values ranging from 4.40 to 7.24 (Supplementary Fig. 29a). We identified three major DsRZ2 haplotypes, exhibiting distinct geographic distributions (Supplementary Table 12). Hap1 was fixed (> 90%) in accessions located in the northeast that experience coldest winter temperatures −7.65°C, whereas Hap3 dominated (87% frequency) in southern populations subjected to −4.77°C (Fig. 3f; and Supplementary Fig. 29b, c). The observed haplotype frequencies also differed across D. sanguinalis varieties, indicating a possible association with subspecies identity (Fig. 3g). Furthermore, extended haplotype homozygosity (EHH) analysis at the DsRZ2 locus did not exhibit significant differences between haplotypes carrying the T or the C allele, which may result from polygenic additive effects (Supplementary Fig. 30; standardized |iHS| score = 0.72).

Given the recurrent introgression observed in Digitaria, we tested whether local adaptation in D. sanguinalis may have been facilitated by introgression. In ABBA-BABA statistics, significantly elevated fd values were detected at the DsRZ2 locus in sympatric populations (Fig. 3h; and Supplementary Fig. 29d). This indicates that in these D. sanguinalis accessions, genetic variation within an environmentally associated region on Chr4 more closely resembled sympatric D. ciliaris haplotypes than those of allopatric D. sanguinalis populations, supporting a hybrid origin for this ~500 kb genomic segment. This introgressed segment contains a functionally coordinated cluster of six stress-adaptive loci, including RZ2 and four genes previously characterized in rice (OsRALyase, OsBIHD1, OsUBC26, BK-PP2A)56–59 (Fig. 3h, i). The presence of these genes in an introgressed block suggests that introgression from D. ciliaris acted as a beneficial reservoir of adaptive alleles.

Temporal escalation and geographic spread of ALS-inhibitor resistance
Nicosulfuron has been widely used for weed management since its registration in the 1990s18. To quantify the temporal dynamics of ALS-inhibitor resistance in Digitaria, we performed nicosulfuron single-dose bioassays on accessions collected from the Huang-Huai-Hai agroecosystems in 2013, 2015, and 2023, respectively (Supplementary Table 13). As expected, resistance levels among surveyed populations exhibited a marked increase over a decade. Survival rates (number alive after treatment) following herbicide treatment rose from 36.5% in 2013 to 74.2% in 2023 (Supplementary Table 13), while the mean herbicide control efficacy concurrently declined from 86.0% to 70.7%. To further characterize resistance variation, we conducted dose-response assays on 196 representative accessions, determining GR50 (herbicide dose causing 50% plant growth reduction) and GR90 values (Supplementary Data 9). The GR50 estimates varied by more than 2000-fold, ranging from 0.083 g a.i. ha-1 (observed in accession #15-9 from Anhui) to 168 g a.i. ha-1 (accession #W-21 from Shandong), reflecting substantial inter-population differences in resistance levels.

Despite being predominantly collected from ecologically comparable environments, different Digitaria species displayed substantial variation in herbicide resistance (Fig. 4a; Supplementary Data 9). Accessions of D. bicornis, primarily sourced from the Hainan Island, exhibited the lowest GR values, indicating high susceptibility to nicosulfuron (Fig. 4a). In contrast, D. ischaemum, D. ciliaris, and D. sanguinalis consistently showed higher GR values (Fig. 4a; and Supplementary Data 9). The median GR90 value of the three weeds reached 84 g a.i. ha-1, with 100 g a.i. ha-1, 80 g a.i. ha-1, and 84 g a.i. ha-1 of D. ischaemum, D. ciliaris, and D. sanguinalis accessions, respectively, above the recommended field application dose (60 g a.i. ha-1). These variations in GR50/GR90 values in surveyed accessions have already rendered nicosulfuron weed control ineffective. To assess the temporal dynamics of resistance development, D. sanguinalis accessions collected from the lower-middle Yellow River region (32°−40°N, 108°−120°E) were used in the analysis. Both GR50 (one-tailed t-test, p = 4.26e-2) and GR90 (one-tailed t-test, p = 4.16e-3) values showed positive correlations with collection years, indicating a progressive increase in resistance over time (Fig. 4b). This trend highlights the rapid adaptation of D. sanguinalis populations to escalating herbicide selection pressure. Accessions exhibiting resistance (GR90 > 60 g a.i. ha-1) were found to be geographically widespread across multiple provinces, highlighting the increased challenges such a distribution poses for effective weed management (Supplementary Fig. 31).

To assess potential TSR, we first examined non-synonymous mutations in ALS, the herbicide target gene of nicosulfuron. Whole-genome sequencing of the Digitaria accessions (n = 554) identified 25 non-synonymous mutations in three DsALS gene copies (Supplementary Fig. 32a; and Supplementary Table 14; Supplementary Data 10). The DsALS-C copy, located on the CH, harbored the highest number of mutations (12 on the CH, 7 on the DH, and 6 on the EH). Most mutations (15/25) were present in a heterozygous state across accessions, in contrast to the predominantly homozygous resistance alleles observed in Echinochloa31 (Supplementary Table 14). These alleles also remained at low frequencies (less than 0.16) (Supplementary Table 14; Supplementary Data 10). No shared mutation sites were observed among three DsALS copies, and no accession carried variants in all three copies, either in homozygous or heterozygous state (Supplementary Data 10). In addition, no accessions exhibited abnormal coverage at any of the three DsALS gene copies (avg. coverage: 0.96 for DsALS-C, Shapiro-Wilk test, p > 0.01; 0.95 for DsALS-D, Shapiro-Wilk test, p > 0.01 and 0.92 for DsALS-E, Shapiro-Wilk test, p > 0.01). While no fixed mutations were identified in these copies, the complexity of the hexaploid genome does not preclude the presence of a dominant mutation. Further studies at both the genetic and functional levels will be necessary to establish whether such a mutation underlies TSR in D. sanguinalis.

Gene amplification-mediated dosage effects represent a mechanism for rapid adaptation to intense selection pressures, as illustrated in blackgrass (Alopecurus myosuroides), and bamboo rats60,61. To determine whether NTSR in D. sanguinalis involves gene amplification, we assessed variation in coverage for six key resistance-related gene families (cytochrome P450s, GSTs, UDP/GTs, AKRs, NB-ARC, and ABC transporters). Most accessions showed coverage levels for resistance-related gene families similar to those of single-copy genes, although 6 to 73 accessions exhibited high-coverage genes across different gene families (Supplementary Fig. 32b). However, correlation analysis between resistance levels and NTSR gene coverage revealed no association with resistance (Supplementary Fig. 32b; and Supplementary Data 11). Overall, no evidence was found to support metabolic resistance mediated by gene amplification in D. sanguinalis, though further investigation is needed to fully exclude this possibility.

To identify genes either constitutively expressed or induced by herbicide, we performed time-course RNA sequencing on leaf tissue from a resistant (R; accession #21-17) and a susceptible (S; accession #15-2) D. sanguinalis accession treated with nicosulfuron (Supplementary Data 12). In the R accession, 152 genes exhibited induced upregulation following herbicide treatment, including AKR1, UGT706D1, CYP81A6, and ABCG43, which are involved in reactive oxygen species (ROS) scavenging and stress response processes62–65. We identified DsSOH1, encoding a short-chain dehydrogenase/reductase (SDR), as a candidate resistance gene. This gene family is widely conserved across plants and implicated in detoxification and abiotic stress responses66,67. It exhibited nicosulfuron-inducible upregulation specifically in the R but not in the S accession (Supplementary Fig. 33). A total of 149 genes exhibited constitutive differential expression, comprising two genes (CYP75B3 and CYP92C21) with characterized homologs in rice involved in biotic stress responses68,69. Additional analyses were described in the Supplementary Note 3 (see Supplementary Figs. 34–36).

Adaptive introgression enables evolutionary escape from herbicide
To identify genomic loci underpinning NTSR, we conducted a genome-wide association study (GWAS) using 141 accessions based on their GR50 values under nicosulfuron treatment70. A total of 40 resistance-associated SNPs (rSNPs) were identified, corresponding to 19 genes, including three homologs (DsCYP81A6, DsAKR4C16 and DsABCC1) functionally validated to confer non-target-site herbicide resistance in other grass species64,71,72 (Fig. 4c). rSNPs were identified 1898 bp upstream of DsSOH1 and 1891 bp downstream of DsCYP81A6, both of which exhibited differential expression in response to nicosulfuron treatment (Fig. 4c; Supplementary Fig. 33).

Herbicide resistance levels were quantitatively assessed using the Resistance Index (RI). Based on RI values, accessions were categorized into four distinct resistance groups: Susceptible (S), Low Resistance (LR), Moderate Resistance (MR), and High Resistance (HR) (Supplementary Data 9). The rSNP located near DsSOH1 was present in 40% of HR accessions but was absent in both LR and S groups (Fig. 4d). Phenotypic assessments further demonstrated that accessions carrying the alternate allele (A) exhibited significantly higher herbicide resistance compared to those with the reference allele (C) (Fig. 4e). To trace the origin of the haplotype at DsSOH1, we conducted ABBA-BABA analysis using the topology (((P1, P2), P3), O), with D. bicornis, susceptible D. sanguinalis, resistant D. sanguinalis, and D. ciliaris were assigned to O, P1, P2, and P3, respectively (Supplementary Fig. 37a). A prominent fd peak was detected specifically at the DsSOH1 locus in the HR group, suggesting introgression from D. ciliaris, whereas the MR and LR groups exhibited progressively weaker signals of admixture (Fig. 4f). The DsSOH1 variant was not fixed in the HR group, with an observed allele frequency of 0.40, supporting a scenario of recent adaptive introgression (Fig. 4d). Concurrently, elevated genome-wide Fst values revealed population subdivision within D. sanguinalis (Fig. 4f). The regions with the highest Fst differentiation between resistant and susceptible populations coincided precisely with the fd peak interval, strongly implicating introgression from D. ciliaris at the DsSOH1 locus (Fig. 4f). Although EHH analysis revealed difference between the A and C alleles, the |iHS| value (1.39, within the 90th percentile of the genome-wide distribution) did not reach the conventional threshold for detecting strong recent positive selection (Supplementary Fig. 37b).

To clarify the ancestry of the introgressed genomic region, we employed Bayesian Phylogenetics and Phasing (BPP v4.0)73 under an explicit introgression model (Fig. 4g). Parameter estimates supported a topology in which the DsSOH1 haplotypes in resistant D. sanguinalis accessions originated from recent introgression with D. ciliaris (φ = 0.30, 2.5% HPD = 0.24, 97.5% HPD = 0.36), while alternative scenarios were statistically rejected (Supplementary Data 13). Moreover, divergence time estimates indicate that speciation between D. sanguinalis and D. ciliaris substantially predated the inferred introgression event (divergence τ = 0.002010 vs. introgression divergence τ = 0.000056), providing further confirmation that resistance alleles originated via admixture. Consistent with the wide geographic distribution of resistant D. sanguinalis accessions, maximum-likelihood phylogenies constructed from introgressed SNPs consistently grouped resistant accessions within clades of sympatric D. ciliaris (Fig. 4h). This recurrent phylogenetic pattern suggests that the NTSR-associated haplotype at DsSOH1 likely evolved through parallel, recent introgression from neighboring D. ciliaris populations.

To investigate the correlation between introgression frequency from D. ciliaris and herbicide resistance levels in D. sanguinalis, we compared the fd values at rSNPs to the genome-wide background. The HR group exhibited significantly higher introgression frequencies at rSNPs compared to the genomic background (Mann-Whitney U test, p = 0.02), a pattern not detected in the MR or LR groups (Supplementary Fig. 37c). In addition, the HR group exhibited the greatest cumulative length of introgressed windows across the genome (HR vs. MR, Mann-Whitney U test, p = 1.60e-2; HR vs. LR, Mann-Whitney U test, p = 3.90e-5) (Supplementary Fig. 37d). These findings suggest that introgression from D. ciliaris may be associated with the emergence and spread of NTSR in D. sanguinalis. Interestingly, we observed a higher frequency of rSNPs in accessions collected in 2023 compared to those collected in 2013 (Supplementary Fig. 37e), which may correlate with the elevated herbicide resistance observed in recent populations.

Methods
Genome sequencing and assembly
Voucher specimens of the three sequenced species were deposited in the Herbarium of Zhejiang University (HZU), with accession numbers HZU60147516 (D. sanguinalis, #YJ2023), HZU60147511 (D. milanjiana, #DZ2), and HZU60147514 (D. radicosa, #YZGJ2).

For D. sanguinalis, total genomic DNA was extracted from young leaves using the cetyltrimethylammonium bromide (CTAB) method. High-molecular-weight DNA was prepared using a nuclei method104 for Nanopore ultralong sequencing and library preparation, followed by sequencing on a PromethION platform. ~116.98 Gb of ultralong reads (N50 > 100 kb) were de novo assembled with NextDenovo v2.5.2 using default parameters105. In parallel, single-molecule real-time (SMRT) sequencing libraries were prepared following PacBio protocols and sequenced on the Sequel II platform using the circular consensus sequencing approach. HiFi reads and ultralong reads were co-assembled using hifiasm v0.19.8-r603, and the contiguity of initial assembly was improved by merging with the initial NextDenovo assembly using Quickmerge106,107, with validation based on long-read mapping depth. The assembly was further polished by PacBio HiFi and Illumina data via NextPolish2 (v0.2.1) with default parameters108. Hi-C data were processed with YaHS v1.1 to generate a chromosome-scale contact map (Supplementary Fig. 4a), which was manually curated in Juicebox v1.11109–111. Subgenome phasing was performed using a k-mer based approach implemented in SubPhaser112.

For D. milanjiana, genomic DNA was extracted from young leaves using the CTAB method. Illumina paired-end libraries were prepared according to the manufacturer’s protocol (Illumina, USA). PacBio long reads were generated, error-corrected, and assembled into contigs using hifiasm. The high-contiguity D. sanguinalis assembly was used as a reference to anchor and order the contigs, and the contigs were further scaffolded into pseudochromosomes corresponding to the DH and EH using Ragtag v2.1.0113.

For D. radicosa (#YZGJ2), DNA extraction, library preparation, and genome assembly followed the same procedures described for D. milanjiana.

Genome annotation
Repeat families were identified de novo and initially classified using RepeatModeler v1.0.10114 followed by genome-wide repeat annotation performed with RepeatMasker v4.0.7114.

Protein-coding genes were predicted using an integrative approach that combined ab initio prediction, homology-based inference, and transcriptome-supported annotation. Ab initio predictions were generated with Fgenesh and AUGUSTUS v3.2.2, whereas homology-based evidence was obtained using GMAP115–117. All evidence types were integrated using EVidenceModeler v1.1.116,118, resulting in a non-redundant consensus gene set. Gene models were retained only if supported by homology to known proteins, transcript alignments, or by at least two independent ab initio predictions. Low-confidence models, defined as those encoding peptides shorter than 50 amino acids or showing significant similarity to repetitive elements in Repbase (E-value < 1e-5, identity > 30%, and coverage > 25%), were filtered out to improve annotation quality.

Functional annotation of predicted protein-coding performed was conducted using InterProScan v5.24-63.0 across Digitaria spp., Echinochloa spp., O. sativa, S. italica and Panicum hallii40,119–122. Homologous genes previously validated in rice and maize were annotated using BLAST, retaining hits with an E value < 1e-5 and sequence identity > 50%53,123.

Repetitive elements annotation
Centromeric satellite repeats were predicted using the Tandem Repeat Annotation and Structural Hierarchy (TRASH) pipeline124 in two iterative rounds. In the first round, genomic sequences were partitioned into 1-kb windows, and local k-mer frequencies were calculated to detect repeat-enriched regions using default parameters. The most abundant repeat templates were clustered and extracted using CD-HIT125. In the second round, TRASH was executed with the parameters “--simpleplot --frep 10 --N.max.div 5 --par 5 --seqt,” with the extracted templates provided via the “--seqt” option. Windows were scored based on the proportion of repeated k-mers, and regions exceeding the threshold were classified as repeat-rich. Tandem repeats within these windows were further characterized based on the spacing between identical k-mers, enabling the identification of consensus repeat units, including CEN113, CEN159, and CEN178, hereafter collectively referred to as cenSat.

Ribosomal DNA (rDNA) regions were annotated via BLAST searches against maize rDNA references sequences, including 5S rDNA (DQ351339), 5.8S rDNA (AF019817), and the intergenic spacer (AF013103). Telomeric repeats were identified by scanning chromosome termini for high-copy-number tandem arrays of the canonical monomer “TTTAGGG”.

Genome quality assessment
Assembly quality for each genome was evaluated based on completeness, correctness, and continuity. For completeness, NGS short reads and PacBio HiFi long reads were mapped to their respective assemblies using BWA mem v0.7.17-r1188126 and Minimap2 v2.03127 with default parameters, respectively. Mapping statistics were summarized with Sambamba Flagstat v1.0.1128. Additional completeness evaluations were performed using the LTR Assembly Index (LAI)129 and BUSCO v5.5.0 with the poales_odb10 database41. Short reads were split into 21-mers using Meryl v1.4 with low frequency k-mers removed (meryl gt 1) prior to generating a k-mer reference library. Then, k-mer based completeness was estimated with Merqury v1.3130. For correctness, base-level accuracy was evaluated by estimating quality values with Merqury, based on the same k-mer library. For continuity, potential assembly gaps and structural inconsistencies were identified using CRAQ44. Long and short reads were aligned against assemblies, and AQI scores were calculated by analyzing clipped alignments to identify local and large-scale errors.

Phylogenetic analysis
Homoeologous exchanges (HEs) in polyploids can distort phylogenetic inference; therefore, candidate HE regions in the D. sanguinalis genome were identified and excluded prior to downstream phylogenetic analyses. ~99× short read coverage from D. radicosa and 50× from D. milanjiana were mapped to the D. sanguinalis reference genome using Bowtie2 with default settings131. Mapping depths were calculated in 100-kb sliding windows across the genome. For the DH and EH, windows with coverage exceeding 20× and showing higher mapping depth for D. radicosa than D. milanjiana were designated as candidate HE regions; the same approach was applied to the CH. Adjacent candidate HE windows were merged, resulting in four distinct HE regions. Genes located within these regions were excluded from subsequent phylogenetic reconstruction.

Phylogenetic relationships among Digitaria and related species were inferred using three approaches. First, genetic divergence was estimated based on synonymous substitution (Ks) values, and pairwise Ks peak distributions were used to infer relative divergence times. Second, a concatenated alignment of 2,030 single-copy orthologs, identified by OrthoFinder132 and excluding HE-associated genes, was used to construct a ML tree. Amino acid sequences were aligned using MAFFT v7.310133 and trimmed with TrimAl under default parameters134. The ML tree was generated with IQ-TREE v1.6.12135 with 1000 bootstrap replicates and the best-fit substitution model (GAMMA+JTT+F4) as determined by ModelFinder. Third, 2,030 single gene trees were inferred individually using RAxML under the best-fitting amino acid substitution models, and a coalescent-based species tree was subsequently estimated with ASTRAL v5.7.8136,137.

Collection and phenotyping of Digitaria accessions
A total of 579 Digitaria accessions were collected from 24 provinces across China from 2013 to 2023 (Supplementary Data 4). Seeds collected prior to 2023 were preserved at the Plant Protection Research Institute, Shandong Academy of Agricultural Sciences. In 2023, 496 accessions were grown under uniform field conditions at the Jiyang Research Station of the Shandong Academy of Agricultural Sciences (Jinan), with five representative individuals planted per accession for phenotypic evaluation. The geographic distribution of collected accessions was visualized using R v4.3.1. Based on detailed collection-site metadata, 272 accessions were further classified into three ecological habitat types: natural habitats (n = 40; riverbanks or wilderness), agricultural habitats (n = 221; active croplands), and disturbed habitats (n = 11; parks or roadsides).

Definition of spatial groups and comparisons
Accessions collected in 2013 or 2015 within a of 100-km radius were clustered into historical groups using the “great_circle” method in the geopy package, with the centroid of each cluster designated as the groups’ location. Historical groups with at least four accessions were retained for temporal comparison. To create modern groups, we selected modern accessions in close proximity (< 100 km) to the historical groups. Finally, 11 groups comprising 167 accessions were used in the analysis (Supplementary Fig. 22). For each spatial group, the normalized Shannon entropy was calculated as a diversity score using R package entropy based on k=6 from the admixture analysis. Nucleotide diversity and Tajima’s D were estimated with VCFtools v0.1.17138 in non-overlapping 100-kb windows. Heterozygosity was calculated for each accession.

Resequencing and variant calling
Genomic DNA was extracted from fresh leaves using a standard CTAB-based protocol. Paired-end resequencing libraries were prepared and sequenced on the DNBSEQ T7 platform. Raw reads were quality-filtered using fastp v0.24.2139, and high-quality paired-end reads were mapped to D. sanguinalis reference genome (#YJ2023) using Bowtie2 with default parameters131. Whole-genome variant calling and filtering were performed using an integrated pipeline140. To minimize false positives due to high sequence similarity between subgenomes, resequencing reads of YJ2023 were realigned to its reference, and variants detected in this process were excluded from the final dataset. Variants were further filtered by the minor allele frequency (MAF) greater than 0.01 and missing rate less than 30%. Functional annotation of all high-confidence variants was conducted using SnpEff v3.652141.

Population structure and genetic diversity
Species identity of each accession was initially assessed by calculating read mapping rates to the YJ2023 reference genome and estimating genome coverage using sambamba v1.5 and bamdst128,142. Population structure was inferred with FastStructure143 based on genome wide SNPs and synonymous SNPs from CH, DH and EH, with k values ranging from 2 to 9. Phylogenetic relationships were reconstructed with FastTreeMP144 based on SNPs from 499 accessions (excluding D. violascens and D. ischaemum due to low genome coverage), using 51 D. ciliaris accessions as the outgroup; branch support was assessed with 1000 bootstrap replicates. Trees were visualized using iTOL v7 (http://itol.embl.de)145. PCA was conducted in PLINK v1.90b6.20146 using a linkage-disequilibrium pruned SNP set (10 SNPs per 50-kb sliding window, r² < 0.5). Nucleotide diversity (π) was calculated with VCFtools v0.1.17138 in non-overlapping 20-kb, 50-kb and 100-kb windows, and genome-wide or subgenome-level diversity was reported as the mean π across all windows within each population.

Demographic history inference based on site frequency spectrum (SFS)
To reconstruct the demographic history of Digitaria and evaluate the role of introgression during species and lineage divergence, we performed composite ML inference based on the site frequency spectrum (SFS). Joint folded two-dimensional SFSs (2D-SFSs) were generated from four-fold synonymous SNPs using easySFS.py (https://github.com/isaacovercast/easySFS). Eight demographic scenarios were designed to test for the occurrence and timing of introgression between lineages and species (Supplementary Fig. 18), with divergence time priors informed by fossil calibrations in TimeTree5147. Likelihood estimation for each scenario was conducted using fastsimcoal2148, with 100,000 coalescent simulations per likelihood estimation (-n 100,000) and 40 expectation-conditional maximization (ECM) cycles (-L 40). Model selection was performed using the Akaike information criterion (AIC), calculated as AIC = 2k + 2ln(L) (MaxEstLhood), where k represents the number of estimated parameters and MaxEstLhood is the ML function value for the model. To avoid convergence to local optima, each analysis was repeated at least twice, and the best-supported model based on AIC was re-run 100 times to obtain refined parameter estimates. Subsequently, 100 independent DNA polymorphism datasets were simulated as joint SFSs conditional on estimated demographic parameters. ML analysis was then applied to each joint SFS over 40 ECM cycles to obtain confidence intervals (CIs) for final estimates.

Changes in effective population size over time were further inferred using Stairway Plot v2.0149, based on both folded and unfolded SFSs. Population-specific SFSs were constructed and folded from the same SNP dataset used for fastsimcoal2.

Phylogeny of chloroplast genomes
The reference chloroplast genome of YJ2023 was assembled de novo using GetOrganelle v1.7.7.1150. For other accessions, chloroplast-derived reads were extracted from paired-end resequencing data with GetOrganelle and mapped to YJ2023 chloroplast reference genome. Variants were called using the GATK pipeline, resulting in 4,742 high-quality chloroplast SNPs after filtering for MAF 0.05 and missing rate 0.2 across 579 Digitaria accessions. A ML tree was constructed in IQtree with 1500 bootstrap replicates135, and the tree was rooted with S. italica as the outgroup.

Detection of sympatric introgression between D. ciliaris and D. sanguinalis
Sympatric population pairs were defined based on pairwise geographic distances calculated with the “great_circle” method in the geopy package, with a threshold of 150-km. Populations pairs comprising at least two D. ciliaris and eight D. sanguinalis accessions were retained for downstream analyses. Introgression in sympatric population pairs was assessed using D-statistics151,152, and the extent of introgressed genomic regions was quantified using the modified fd statistic153. The phylogenetic topology (((P1, P2), P3), O) was applied, with D. bicornis designated as the outgroup (O), sympatric D. sanguinalis as P2, sympatric D. ciliaris as P3, and allopatric D. sanguinalis accessions as P1. Under the null hypothesis, ABBA and BABA site patterns are expected to occur at equal frequencies due to incomplete lineage sorting, whereas an excess of ABBA sites indicates introgression between P2 and P3. To minimize sample bias, eight individuals were randomly selected for P1 and P2, and three for P3 and the outgroup.

Introgression was evaluated in two steps. First, D-statistics were calculated genome-wide. To account for linkage disequilibrium, significance was assessed using a block jackknife approach151, in which the genome was divided into 1-Mb blocks, and one block was sequentially removed at a time to estimate the mean and variance of D. Second, introgressed genomic regions were identified using sliding windows of 10-, 50-, and 100-kb, with corresponding step sizes of 5-, 25-, and 50-kb, respectively. Windows containing fewer than 200, 100, or 20 SNPs (for 100-, 50-, and 10-kb windows, respectively) were excluded. Both D and fd statistics were computed with “ABBABABAwindows.py” and Fst values for the same windows were calculated with “popgenWindows.py”153. According to the filtering criteria used by Zhou et al154., fd values outside the 0-1 range were set to zero. The overall introgression level was quantified as the Proportion of Introgression across the Genome (PIG) and Pairwise Shared Introgressed Genomic regions (PSIG) was used to quantify shared introgression among sympatric populations. Regions with fd > 0.6 were defined as significantly introgressed, and genes within these regions were subjected to functional enrichment analysis.

Functional enrichment analysis of introgressed genes
Gene families were annotated based on Pfam domains identified in the YJ2023 genome. For enrichment analysis, all Pfam domains in YJ2023 were used as the background set. To test whether the Pfam domains found within introgressed regions represented a non-random sampling, we applied Fisher’s exact test to compare their frequencies against the genome-wide background. Multiple testing correction was performed using the Benjamini-Hochberg and domains with adjusted p-values < 0.05 were considered significantly enriched. Fold enrichment was calculated as the ratio of each Pfam domain’s frequency in introgressed regions to its frequency within the background.

Environmental factors and correlation with introgression
To investigate the relationship between introgression and environmental variation, 91 environmental variables were retrieved from the WorldClim database (http://www.worldclim.org/) (Supplementary Data 5). These variables encompassed four primary categories, temperature, precipitation, wind speed, and solar radiation, as well as 19 additional bioclimatic variables. Pairwise environmental distance matrices for the primary categories were calculated using Euclidean distances with the R package “ecodist”155, and a separate matrix was generated from the 19 bioclimatic variables. Longitude and latitude were also used to generate distance matrices. Geographic distances among populations were computed using the “great_circle” method in the “geopy” package, and statistical significance in downstream analyses was assessed using 10,000 permutations.

Mantel and partial Mantel tests were conducted using “ecodist” to assess correlations between PSIG and pairwise geographic, environmental and Fst distances.

Identification of environment-associated genetic variants
Variants with MAF > 0.01 and missing rate > 0.2 were retained, yielding a dataset of 3,981,523 SNPs for downstream analyses. Associations between allele frequencies and 19 environmental variables were initially assessed using a univariate latent factor mixed model (LFMM) implemented in the R package LEA v3.14.0156. Latent factors were determined based on the ancestry clusters inferred from FastStructure, consistent with results from the demographic history analysis. For each environmental variable, five independent runs were conducted, with 5000 burn-in iterations followed by 10,000 sampling iterations. Median p values from five independent runs were calculated and adjusted for multiple testing using both the false discovery rate (FDR, 5%) and the Bonferroni method.

Single-dose bioassays on Digitaria
A single-dose herbicide bioassay was conducted on Digitaria accessions (n = 53 in 2013, n = 56 in 2015, n = 52 in 2023) at the 2-3 leaf stage. Plants were grown in 9-cm diameter pots filled with sterilized potting soil and maintained under controlled environmental conditions (30/25°C day/night, 16h of photoperiod) until herbicide application. A commercial formulation of nicosulfuron (field-recommended dose: 60 g a.i. ha-1) was applied using an enclosed cabinet sprayer calibrated to deliver 450 L ha-1 at a pressure of 0.4 MPa (ASS-5, Information Technology Research Center, Beijing, China). Untreated control (0 g a.i. ha-1) was also included for comparison. Each treatment, including the control, was replicated four times with one pot per replicate.

Survival rate (%) was calculated as:
Survival rate (%) = (Ns / NT) × 100% (1)
Ns represents the number of surviving plants, and NT represents the total number of plants per pot. Plants were classified as surviving if they exhibited active regrowth or retained green tissue.

Aboveground biomass from each pot was harvested, oven-dried at 70°C for 72 h, and weighed. Control efficacy (%) was calculated as the percentage of biomass reduction relative to the untreated control:
Control efficacy (%) = (1 - Bt / Bc) × 100% (2)
Bt is the biomass in the treatment, and Bc is the biomass in the control.

Whole-plant resistance bioassay
Seeds of 196 Digitaria accessions were collected, air-dried, and stored at 4°C prior to use. Seeds were sown in moist loam soil in 9-cm diameter pots, covered with 1 cm of soil, and grown under greenhouse conditions (30/25°C day/night, 16h of photoperiod) with sub-irrigation.

Herbicide resistance to nicosulfuron was evaluated using whole-plant bioassays. At the three-leaf stage, ten seedlings were foliar treated with the commercial formulation of nicosulfuron at seven doses, 0, 30, 60, 120, 240, 480, and 960 g a.i. ha-1 (field-recommended dose: 60 g a.i. ha-1). Applications were performed using an enclosed cabinet sprayer (ASS-5, Information Technology Research Center, Beijing, China) calibrated to deliver 450 L ha-1 at 0.4 MPa. Plant survival and aboveground dry biomass were assessed 21 days after treatment. Each treatment included four biological replicate pots (10 seedlings per pot) and was independently repeated.

Dose-response relationships for biomass reduction (expressed as a percentage of untreated controls) were analyzed by nonlinear regression in SigmaPlot v12.5 (Systat Software Inc.), fitting a four-parameter log-logistic model:
Y = C + (D - C) / (1 + (X / GR50)^b) (3)
where Y is the response (percentage of control), X is the herbicide dose, C and D are the lower and upper asymptotes, GR50 is the herbicide dose required to plant growth reduction by 50%, and b is the slope of the curve. A mixed-model ANOVA was conducted to assess differences in percent control and biomass reduction across treatments. Resistance indices (RIs) for individual accessions were calculated as the ratio of their GR50 values to the average GR50 of susceptible accessions, which were defined as those controlled at the recommended field dose (GR90 < 60 g a.i. ha⁻¹). Resistance levels were categorized into four groups based on RI values: S (susceptible, RI < 2), LR (low resistant, RI > 2 and RI < 4), MR (moderate resistant, RI > 4 and RI < 10), and HR (high resistant, RI > 10).

Herbicide resistance analyses
ALS genes were identified in the D. sanguinalis genome using BLASTP. The predicted ALS protein sequences were aligned with orthologous sequences from O. sativa and A. thaliana using MAFFT133, and known resistance-associated sites were annotated across the three DsALS copies. A comprehensive catalog of causative mutations was then established for all accessions (Supplementary Data 10). Genealogical relationships among haplotypes were reconstructed using HapNetworkView based on variant profiles across the three ALS copies157.

NTSR-related genes were identified based on Pfam domain annotations. Sequencing depth for each annotated gene was calculated, and gene family abundance was estimated as the sum of average depths across all loci, using custom scripts. Coverage was normalized by dividing the average read depth of each gene family by the mean sequencing depth of 2030 single-copy genes. Accessions with sequencing depth < 5 were excluded to minimize bias in coverage estimation.

Plant material and nicosulfuron treatment for RNA-seq
The nicosulfuron-resistant (#21-17, GR50 = 122.53 g a.i. ha-1) and susceptible (#15-2, GR50 = 0.59 g a.i. ha-1) D. sanguinalis accessions were selected for RNA-seq analysis. Plants were grown in 9-cm diameter pots filled with moist loam soil, covered with a 1-cm soil layer, and maintained under greenhouse conditions (30/25°C day/night, 16h of photoperiod) with sub-irrigation at the Jiyang Research Station of the Shandong Academy of Agricultural Sciences (Jiyang, Shandong, China).

At the three-leaf stage, seedlings were treated with 2 g a.i. ha-1 nicosulfuron based on the GR50 value distribution of susceptible populations from dose-response assays, to assess herbicide response. Each treatment included three biological replicates, with three individual plants pooled to constitute one replicate. Leaf tissues were collected at three time points: 0h (untreated control), 6h (early response), and 24h (late response) post-treatment. Samples were immediately frozen in liquid nitrogen and stored at −80°C until RNA extraction and sequencing.

RNA sequencing and transcriptome analysis
Total RNA was extracted from 18 samples, comprising three time points (0h, 6h, and 24h) for each accession, with three biological replicates per time point. Messenger RNA was purified using poly-T oligo-attached magnetic beads, and sequencing libraries were constructed according to the DNBSEQ standard protocol. Libraries were pooled based on effective concentration and target data yield. The 5′ ends of the libraries were phosphorylated and circularized, followed by rolling circle amplification to generate DNA nanoballs, which were subsequently loaded onto a flow cell for sequencing on the DNBSEQ-T7 platform.

After trimming adapter sequences and filtering out low-quality reads, clean reads were mapped to the YJ2023 reference genome using Hisat2 v2.1.0158 and gene expression levels were quantified by StringTie v2.2.1159 with default parameters. DEGs were identified with the pyDESeq2 package v0.5.2160, with significance thresholds set at |log2(fold change)| ≥ 1 and FDR-adjusted p-value < 0.05. GO enrichment analysis was performed using the clusterProfiler package161.

GWAS of nicosulfuron resistance
GWAS were conducted on 141 D. sanguinalis accessions using GR50 values using EMMAX70. SNPs were filtered with MAF > 0.05 and missing genotype rate < 0.1. A pairwise genetic distance matrix, derived from simple matching coefficients of SNPs, was used to define the variance-covariance structure of random effects in the linear mixed model. Median p-values from five independent runs were calculated and adjusted for multiple testing using both the false discovery rate (FDR, 5%) and the Bonferroni correction. Manhattan plots were generated using the R package “qqman”162 for visualization of association signals.

To examine selection pressures acting on adaptive variants associated with herbicide resistance and climatic adaptation, extended haplotype homozygosity (EHH) was assessed for strongly associated loci, and the integrated haplotype score (iHS) was calculated for common variants using hapbin v1.3.0163.

Introgression analysis at resistance associated loci
To assess potential introgression from D. ciliaris into resistant D. sanguinalis accessions, we applied the ABBA-BABA test under a four-taxon topology of the form ((P1, P2), P3, O), where P1 represented the susceptible D. sanguinalis population, P2 comprised resistant accessions (including HR, MR, and LR), P3 corresponded to D. ciliaris, and D. bicornis served as the outgroup (Supplementary Fig. 37a). We estimated the fd values across the genome with sliding windows of 10-, 50-, and 100-kb and respective step sizes of 5-, 25-, and 50-kb. Windows containing fewer than 200, 100, or 20 SNPs (for 100-, 50-, and 10-kb windows, respectively) were excluded. For each window, we also calculated Fst between susceptible and resistant D. sanguinalis accessions. To identify windows with significantly elevated introgression, we performed two-tailed Z-tests and defined outlier windows as those with |Z-score| > 4. To establish significance thresholds for Fst, we conducted randomization tests by permuting the group labels, generating a null distribution from which one-tailed p-values were derived (e.g., one-tailed p = 0.001 corresponds to the 0.1% or 99.9% quantile range).

To assess the phylogenetic structure and potential introgression events, we applied coalescent-based modeling using BPP v4.6.273. Analyses were based on a 100-kb genomic region flanking the rSNP (Chr20: 34,032,644). First, to test for the presence of introgression under the MSC-I model, we implemented the A00 method, which estimates the posterior probabilities of divergence times and introgression probabilities (Supplementary Data 13). This approach allowed us to quantify the likelihood of introgression between D. ciliaris and resistant D. sanguinalis accessions within this focal genomic region. In parallel, we used the A11 method to evaluate whether the hypothesized species tree topology receives the highest posterior probability. All analyses were conducted using three independent MCMC runs with 10,000 burn-in iterations followed by 100,000 samples drawn every 10 iterations to ensure convergence and robustness.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Things that are good:\nThis manuscript covers an enormous amount of data – with its 4 multipaneled figures and 34 supplementary figures, 24\nsupplementary tables, 3 supplementary notes and 4 supplementary datasets that cover genome assembly and annotation of\nD. sanguinalis and progenitors, population genomics of 579 accessions (structure, diversity, demography), evidence and\nfunctional analysis of introgression, environmental adaptation (cold tolerance, genotype-environment associations),\nherbicide resistance phenotyping with dose-response curves and temporal and geographic trends, analysis of target-site vs.\nnon-target-site resistance in the populations, GWAS for NTSR loci (40 SNPs and associated candidate genes),\ntranscriptomics showing DEGs in resistant vs. susceptible, correlation of introgression with resistance, and full discussion of\nevolutionary and management implications – and yet it reads clearly. It is a long and detail-rich paper, which could be split\ninto different papers which focus on each aspect, but I am not sure this is the only or best approach. Overall, it presents a\nsubstantial body of well-connected research that will be useful for and interesting to many researchers and with its wide-\nranging well-covered topics, will be cited widely.\nThe introduction clearly and accurately states the rationales behind doing this research. The authors make it clear why the\ngenome of Digitaria is needed (e.g. so that taxonomic ambiguities can be resolved and that the molecular nature of its NTSR\ncan be investigated) alongside the bigger pictures of using it as a model for polyploidization-driven acquisition of\necologically important traits. However, I would recommend a small rearrangement of these paragraphs to improve the flow of\nthe introduction. See below.\nThe description of the data and abbreviations as well as the signposting between sections are all good and helpful to guide\nthe reader through the massive data content. One suggestion is to explain the “CD”, “DT” and “ET” in the text (see below). I\nwould also ask that the authors check that full names of the species are used in the first instance of usage as well as in the\nfigure legends.\nAs considerable text is dedicated to the physical descriptions of these Digitaria, Supplementary Figures 16 and 17 are very\nuseful – thank you for including them. There are several places where much text is devoted to supplementary figures, so\nplease ensure that the supplementary data is readily available.\nThroughout the analysis and results are presented clearly and discussed well – occasionally the cross referencing needs to\nbe double checked to ensure that the right figure number and panel is being referenced. The numerical and statistical data\npresented make it clear that this is a high-quality body of work that will be widely useful for current and future studies in\nDigitaria. These data are put into context and made useful for a generalist audience with the details in the supplementary\nfigures / text / data / notes.\nThe Discussion, albeit brief compared to the weighty results, is helpful to bring the different datasets together. It is largely\nfocused on the herbicide resistance phenotypes rather than the other data presented, but this is understandable as Digitaria\nis largely “a problem” because of its herbicide resistance.\nThings to consider:\nThe abstract rather undersells the depth and breadth of the data – the authors may want to rewrite this to ensure that they\nsell all of the results appropriately to the scientific audience who inevitably be drawn in or not by a skim-read of the abstract.\nI would suggest a small rearrangement of the introduction to improve its flow. Now the paragraphs read: 1) weeds\n(especially Digitaria) are bad, 2) polyploidization may drive adaptive plasticity, 3) TSR and NTSR are present in Digitaria,\nyet not fully understood, then 4) interspecific hybridisation may underpin the wide geographic distribution and rapid\nevolutionary responses of Digitaria. I suggest that an order of 1) weeds (especially Digitaria) are bad, 3) TSR and NTSR are\npresent in Digitaria, yet not fully understood, 2) polyploidization may drive adaptive plasticity, and 4) interspecific\nhybridisation may underpin the wide geographic distribution and rapid evolutionary responses of Digitaria would be more\nlogical.\nThe map in Figure 2b is incredibly interesting, but a bit difficult to read as there is considerable overlap between the shapes.\nPerhaps it could be cropped to focus in on the area that was sampled rather than trying to include all of China.\nOn line 381, the authors state “especially considering the rarity of natural grafting in grasses” which seems out of context.\nPerhaps the authors could explain.\nThings that need to be corrected:\nTheres is consistent usage of “hexaploidy” when perhaps “hexaploid” is probably more correct. As in “T2T genomes of\nhexaploidy D. sanguinalis….”\nIt would help the reader if the text in lines 188-191 explained that “CD” refers to the C genome from the diploid and “DT” and\n“ET” as line 158 explains “H” denotes hexaploid (not “hexaploidy”) origin.\nWithin the text, figures and their legends, please ensure that the full species names are given in the first instance usage in\nthe text (e.g. line 196 Oryza sativa rather than O. sativa) and again in the figure legends. This is important not only to follow\n“publishing standards”, but also “E. indica” in Figure 1c as Eleusine indica and Eragrostis indica are both applicable and\nprobable. Perhaps some explanation on why these other grass species were chosen might improve the text. It was\nsurprising that “E. indica” was used in Figure 1c and then not in Figure 1d – perhaps this could be explained.\nI believe that the text in line 249 should refer to Figure 2b not 2a and then on line 250, it is referring to 2a not 2b. If that is the\ncase, perhaps consider restructuring the figure so the map is on the left and the genome coverage vs read maps and\nunrooted maximum-likelihood phylogenetic tree reversed. Similarly, there are several places where Figure 2f is referred to\n(e.g. lines 300, 302,) where I suspect the authors were intending to highlight Figure 2e. Because of the huge number of\nsubpanels and Supplementary Figures with subpanels, can the authors please double check that each and every cross\nreference is correct. Also, I suggest that the panel layout reflects the order in the text - Figure 3e and 3f are only referenced\non line 461 while 3g on 429 and 3h on 448… I understand that some of this is just making them visually appealing and to\nmake the figures fit in the space, but with so many essential Supplementary Figures and subpanels, it makes it slightly more\ndifficult jumping back and forth between the text and the supplements.\nAs far as I can tell, Figure 2c is not referred to anywhere in the text. Nor is Figure 4f. Nor is Figure 4g. I may have missed\nthese, but I did not find them despite trying. Moreover, I did not check each and every cross reference; please ensure that\nEACH of the figures and presented data are referenced, discussed, and analysed. If not… there are several instances where\nthe Supplementary Figures are covered by considerable amounts of text or even get their own sections (e.g. Supplementary\nFigures 21, 23). I’d ask the authors to consider if there are places to rethink the figures or rearrange them to make space for\ndata that is essential to your story.\nPlease make sure that in the results section discussing the candidate genes associated with local adaptation, that all\nreferences and data are cross referenced. Line 444 talks about expression analysis of DsRZ2 across “different tissues” but\nas far as I can tell you haven’t discussed sampling different tissue types.\nPlease move the label the Y axes for Figure 4b from above to the more conventional next to the axis location.\nThere is a strange symbol on line 567 in the text describing allele frequency, please correct.", "reviewer_2": "The manuscript by Huang et al., entitled “The Digitaria genomes reveal local adaption and herbicide resistance mediated by\nintrogression” presents an impressive dataset and array of analyses to study the genomic diversity, evolutionary history, and\nadaptation to herbicides in a pressing agricultural weed, Digitaria sanguinalis, commonly known as crab grass. In producing\nvery high quality reference genomes for not only the focal species but its diploid progenitors, the authors provide a valuable\nresource for systematically classifying Digitaria accessions in agricultural landscapes. The expansive analyses of an\nequally impressive population genomic dataset include model based inference of demographic history, population structure,\nintrogression gene-environment association, bio-assays for herbicide tolerance and accompanying GWA— the analyses on\nthe genetic architecture of non-target site resistance are especially impressive and convincing. However, the motivation for\nseveral of the analyses presented are difficult to follow, and the conclusions given throughout the manuscript seem to\noverstep the results presented. Generally, some work to improve the flow of writing, including explicitly stating the\nhypotheses being tested for each analysis, would help readability of this elaborate work.\nLine 229-230, 233-234: It seems very tentative to suggest gene family expansions relate to herbicide resistance when this\nanalysis was at the genus level, and there is plenty of variation in herbicide use / resistance levels within and across\nspecies. To claim this, one would want to compare lineages contrasting in a history of herbicide use to see if this result\nstands up. Similarly, that UDPs/GTs are enriched on Chromosome 3, which happens to also have the highest TE density,\ndoes not convincingly show that TEs have driven an adaptive expansion of this gene family.\nLines 348-349: It would be useful to provide more context to help the reader interpret these results in the main text - what\nwas the temporal scale of longitudinal surveillance? How many years were sampled? Were the same localities revisited?\nEither way, I would want to make sure that geographic coordinates were controlled for when conclusions are made about\nstructure changing through time to be sure it's not just driven by differences in space.\nLines 384: The authors claim early on that there was evidence for pervasive introgression between sympatric populations,\nbut more context in how these tests were performed would really help with interpreting these results. Were these tests\nbetween the diploid, tetraploid, and hexaploid populations? Was there more introgression in some directions compared to\nothers? The authors note there is evidence for introgression in 7 tests, but how many total populations were tested?\nLine 389: the acronym PIG is used before it is defined\nLine 393-398: It is hard to interpret these results in the absence of a well defined null expectation. I would like to see a test of\nwhether correlations between PSIG and geographic distance / environmental divergence are greater than expected from the\ngenome-wide background levels of divergence /population structure. The fact that environmental divergence did not remain\na significant predictor after accounting for geographic distance to me, suggests that this pattern is in part due to population\nstructure, as we know it to be strongly structured by geography. Therefore, I do not believe this is evidence for\nadaptive/environment mediated introgression.\nLine 417-419: As there has been work focused on these gene families in NTSR, it seems a bit like an over\nsimplification/cherry-picking to mention only these gene families when there are supposedly 3500 genes introgressed. I also\nwas not able to find details in the methods about how these tests were performed. Generally, I find these results to provide\nweak evidence of environment-driven selection shaping introgression patterns, which is concluded in the final sentence. It\nseems like overlaying this work with the GEA results in the next section would first be necessary to say this. This is done for\nonly one gene – DsRZ2, but I would like to see a more systematic comparison of introgression regions with respect to\nintrogressed regions to support this more general conclusion the authors make.\nLine 436: Genomics islands typically refer to a debate centered around speciation and reproductive isolation, not adaptation.\nA more standard contrast would be polygenic adaptation vs. selective sweeps.\nLine 450: Can you provide some context for interpreting the lack of significant difference between the T and C haplotype?\nLine 473-474: How do we interpret survival rate vs control efficacy? I’m not sure I understand the difference.\nLine 497: There is no data presented to support the conclusion that resistance has evolved independently multiple times.\nThat resistance accessions were widespread could be the result of gene flow and recombination onto local backgrounds.\nLine 505: Are any of the mutations non-synonymous? This sentence could also use grammatical edits, and the entire\nmanuscript should be carefully edited for similar grammatical errors.\nLines 516-518: The authors ignore dominance in this discussion - a single dominant mutation present on only one genome-\ncopy could confer resistance, even in a hexaploid, so this does not conclusively rule out target-site resistance mechanisms.\nLine 523: I think in this case it would be more accurate to say gene amplification rather than expansions? I think that the\nGWAS is done based on coverage, so it would be more descriptive/informative to you “assessed variation in coverage”\nrather than copy number. To help with readability more generally, I think this section could be shortened given the lack of\nevidence for resistance associated with copy number variation.\nLine 553: What is a GWAS prioritized gene?\nLine 565: Again, what was the source population of the fd peak at DsSOH1? Where do we think this haplotype came from?\nHow many tests gave such a result?\nLine 572/589: DsSOH1 is one of the numerous (~40) NTSR-associated loci you have detected. Therefore it feels like an\noverstatement to conclude that introgression is a driving force behind NTSR. Do any of the other loci show such a signal? Is\nit beyond what you would expect compared to the frequency of introgression across the rest of the genome?\nLine 595-598: Unclear how these two analyses relate to one another\nLine 681-684: I do not follow this sentence for both conceptual and grammatical reasons.\nLine 687: I do not think this novelty statement is necessary or true. Adaptive introgression of herbicide resistance has been\nshown in Teosinte from Maize (https://pubmed.ncbi.nlm.nih.gov/32989136/) and in numerous other systems in response to\npesticides more broadly."}, "author_response": "A Point-to-point response to the reviewers’ comments\nReviewers’ comments are shown as black text. Response to comments is shown as blue text.\nReviewer #1 (Remarks to the Author):\nThings that are good:\nThis manuscript covers an enormous amount of data – with its 4 multipaneled figures\nand 34 supplementary figures, 24 supplementary tables, 3 supplementary notes and 4\nsupplementary datasets that cover genome assembly and annotation of D. sanguinalis\nand progenitors, population genomics of 579 accessions (structure, diversity,\ndemography), evidence and functional analysis of introgression, environmental\nadaptation (cold tolerance, genotype-environment associations), herbicide resistance\nphenotyping with dose-response curves and temporal and geographic trends, analysis\nof target-site vs. non-target-site resistance in the populations, GWAS for NTSR loci\n(40 SNPs and associated candidate genes), transcriptomics showing DEGs in resistant\nvs. susceptible, correlation of introgression with resistance, and full discussion of\nevolutionary and management implications – and yet it reads clearly. It is a long and\ndetail-rich paper, which could be split into different papers which focus on each\naspect, but I am not sure this is the only or best approach. Overall, it presents a\nsubstantial body of well-connected research that will be useful for and interesting to\nmany researchers and with its wide-ranging well-covered topics, will be cited widely.\nThe introduction clearly and accurately states the rationales behind doing this\nresearch. The authors make it clear why the genome of Digitaria is needed (e.g. so\nthat taxonomic ambiguities can be resolved and that the molecular nature of its NTSR\ncan be investigated) alongside the bigger pictures of using it as a model for\npolyploidization-driven acquisition of ecologically important traits. However, I would\nrecommend a small rearrangement of these paragraphs to improve the flow of the\nintroduction. See below.\nResponse:\nWe appreciate your time and constructive comments. We find your comments very\nhelpful for improving our manuscript and have tried our best to address them below.\nThe description of the data and abbreviations as well as the signposting between\nsections are all good and helpful to guide the reader through the massive data content.\nOne suggestion is to explain the “CD”, “DT” and “ET” in the text (see below). I\nwould also ask that the authors check that full names of the species are used in the\nfirst instance of usage as well as in the figure legends.\nResponse:\nWe are sorry for the confusion caused by the incomplete information and the improper\nabbreviation of species names in the manuscript and figure legends. These issues have\nbeen thoroughly revised and corrected in the updated version. Abbreviations used for\ndifferent subgenomes (“CD”, “DT”, and “ET”) have been revised to “CC”, “DT”, and\n“ET”, respectively, to improve clarity and consistency and an explanatory note has\nbeen added at their first mention in the revised manuscript, as follows:\n“In the same way, two high-quality genomes for D. radicosa (489.39 Mb, CC\ngenome) and tetraploid D. milanjiana (908.85 Mb, with subgenome sizes of 421.50\nMb (DT, subgenome D in the tetraploid) and 469.45 Mb (ET) were also generated in\nthis study (see Methods and Supplementary Note 1 for details).”\nAs considerable text is dedicated to the physical descriptions of these Digitaria,\nSupplementary Figures 16 and 17 are very useful – thank you for including them.\nThere are several places where much text is devoted to supplementary figures, so\nplease ensure that the supplementary data is readily available.\nResponse:\nThanks for your valuable suggestion. Accordingly, we have carefully reviewed and\nupdated the supplementary files in the revised version to ensure they are fully\naccessible and easily viewable.\nThroughout the analysis and results are presented clearly and discussed well –\noccasionally the cross referencing needs to be double checked to ensure that the right\nfigure number and panel is being referenced. The numerical and statistical data\npresented make it clear that this is a high-quality body of work that will be widely\nuseful for current and future studies in Digitaria. These data are put into context and\nmade useful for a generalist audience with the details in the supplementary figures /\ntext / data / notes.\nResponse:\nWe sincerely thank you for your positive comments. Following your suggestion, we\nhave carefully re-checked all cross-references to figure numbers and panel labels\nthroughout the manuscript, including the main text and supplementary materials.\nCorrections have been made where necessary to ensure that each reference accurately\ncorresponds to the correct figure or panel. For example, the Supplementary Data\nfiles have been reordered to reflect the correct sequence. Specifically, the sentence:\n“These tests aimed to quantify the relative contributions of geographic distance or\nenvironmental divergence to observed introgression patterns (Supplementary Fig.\n23c; Supplementary Data 3).” now correctly refers to the revised Supplementary\nData 2.\nThe Discussion, albeit brief compared to the weighty results, is helpful to bring the\ndifferent datasets together. It is largely focused on the herbicide resistance phenotypes\nrather than the other data presented, but this is understandable as Digitaria is largely\n“a problem” because of its herbicide resistance.\nResponse:\nWe thank you for this insightful comment. As noted, the Discussion section is\nrelatively concise and primarily centers on herbicide resistance, which reflects its\ncritical importance in the context of Digitaria as a major agricultural weed. In\nresponse, we have slightly expanded the Discussion on polyploidization to provide a\nmore balanced synthesis of our findings.\nThings to consider:\nThe abstract rather undersells the depth and breadth of the data – the authors may\nwant to rewrite this to ensure that they sell all of the results appropriately to the\nscientific audience who inevitably be drawn in or not by a skim-read of the abstract.\nResponse:\nFollowing your suggestion, we have substantially revised the abstract to better reflect\nthe depth and breadth of our datasets and the significance of our findings. The revised\nabstract now more clearly highlights the multiple dimensions of our study, including\ngenome assembly quality, structural and evolutionary analyses, and the origins of\nherbicide resistance.\nI would suggest a small rearrangement of the introduction to improve its flow. Now\nthe paragraphs read: 1) weeds (especially Digitaria) are bad, 2) polyploidization may\ndrive adaptive plasticity, 3) TSR and NTSR are present in Digitaria, yet not fully\nunderstood, then 4) interspecific hybridisation may underpin the wide geographic\ndistribution and rapid evolutionary responses of Digitaria. I suggest that an order of 1)\nweeds (especially Digitaria) are bad, 3) TSR and NTSR are present in Digitaria, yet\nnot fully understood, 2) polyploidization may drive adaptive plasticity, and 4)\ninterspecific hybridisation may underpin the wide geographic distribution and rapid\nevolutionary responses of Digitaria would be more logical.\nResponse:\nThank you for your constructive suggestion. We agree that reordering the paragraphs\nwill improve the logical flow and enhance the clarity of the narrative.\nIn the revised version, we have followed the suggested order: “1) Weeds are\necologically resilient organisms that thrive in agricultural ecosystems, exhibiting\nbroad adaptability and robust population dynamics…2) The escalating use of\nagrochemicals has consequently driven the evolution of herbicide resistance in\npopulations of 273 weed species worldwide…3) Polyploidy is common across the\nplant kingdom and is particularly prevalent among weedy species, where it often\ncontributes to enhanced adaptive plasticity…4) Hybridization is common in nature”\nWe believe this adjustment strengthens the coherence of the introduction and provides\na clearer rationale for the aims of our study.\nThe map in Figure 2b is incredibly interesting, but a bit difficult to read as there is\nconsiderable overlap between the shapes. Perhaps it could be cropped to focus in on\nthe area that was sampled rather than trying to include all of China.\nResponse:\nWe are sorry for the information loss caused by overlapping elements. We agree that a\ncropped figure is needed to better focus on the regions where sampling occurred. This\nadjustment reduces overlap among symbols and improves the clarity of the\ngeographic distribution patterns.\nFig. 2a | Geographic distribution of Digitaria accessions used in this study.\nGeographic distribution of Digitaria accessions used in this study. The main\ncollection region is displayed in detail, with an inset in the upper-left corner showing\nthe full geographic extent of sampling across China. Different species are denoted by\ndistinct symbols, while D. sanguinalis varieties are indicated by colored circles.\nRange boundaries of major taxa are shaded.\nOur dataset includes multi-year samples from Huang-Huai-Hai region in China from\n2013 to 2023. We have further supplemented this information with additional\nsupplementary figures to detail the collection specifics across different years:\nFig. R1 | Distribution of Digitaria accessions collected from 2013 to 2023 across\nHuang-Huai-Hai region.\nOn line 381, the authors state “especially considering the rarity of natural grafting in\ngrasses” which seems out of context. Perhaps the authors could explain.\nResponse:\nSorry for the confusion caused by the lack of clarity. In this section, our intention was\nto integrate nuclear and plastid phylogenies to propose preliminary hypotheses\nregarding the reticulate evolutionary history of Digitaria. Either sexual hybridization\nor asexual grafting can lead to chloroplast capture or horizontal transfer (Bock, 2010).\nHowever, monocots are generally considered graft-incompatible due to the absence of\na vascular cambium and the dispersed arrangement of their vascular bundles (Feng et\nal., 2024; Gaut et al., 2019; Melnyk and Meyerowitz, 2015). Thus, the observed\ncytonuclear discordance is more plausibly explained by post-speciation hybridization.\nWe have revised the relevant text to better highlight the inference of widespread\nhybridization events in Digitaria and to improve semantic clarity.\n“Overall, the chloroplast phylogenies indicate extensive interspecific introgression\nacross Digitaria, likely driven by historical hybridization events, given the\nimprobability of natural grafting in grasses (Feng et al., 2024; Gaut et al., 2019).”\nReferences for this response:\nBock, R, 2010. The give-and-take of DNA: horizontal gene transfer in plants. Trends\nPlant Sci. 15, 11–22.\nFeng, M., Augstein, F., Kareem, A., Melnyk, C.W., 2024. Plant grafting: molecular\nmechanisms and applications. Mol. Plant, Focus Issue on “Unique Aspects of Plant\nSciences” 17, 75–91.\nGaut, B.S., Miller, A.J., Seymour, D.K., 2019. Living with two genomes: grafting and\nits implications for plant genome-to-genome interactions, phenotypic variation, and\nevolution. Annu. Rev. Genet. 53, 195–215.\nMelnyk, C.W., Meyerowitz, E.M., 2015. Plant grafting. Curr. Biol. 25, R183–R188.\nThings that need to be corrected:\nThere is consistent usage of “hexaploidy” when perhaps “hexaploid” is probably more\ncorrect. As in “T2T genomes of hexaploidy D. sanguinalis….”\nResponse:\nSorry for the mistakes and we have made corrections in the main text.\nIt would help the reader if the text in lines 188-191 explained that “CD” refers to the\nC genome from the diploid and “DT” and “ET” as line 158 explains “H” denotes\nhexaploid (not “hexaploidy”) origin.\nResponse:\nThank you for your feedback. Abbreviations used for different subgenomes (“CD”,\n“DT”, and “ET”) have been revised to “CC”, “DT”, and “ET”, respectively, to improve\nclarity and consistency, and an explanatory note has been added at their first mention\nin the revised manuscript, as follows:\n“In the same way, two high-quality genomes for D. radicosa (489.39 Mb, CC\ngenome) and tetraploid D. milanjiana (908.85 Mb, with subgenome sizes of 421.50\nMb (DT, subgenome D in the tetraploid) and 469.45 Mb (ET) were also generated in\nthis study (details see Methods, Supplementary Note 1).”\n“D. radicosa CD” in Figure 1 has been revised to “D. radicosa” to improve clarity.\nWithin the text, figures and their legends, please ensure that the full species names are\ngiven in the first instance usage in the text (e.g. line 196 Oryza sativa rather than O.\nsativa) and again in the figure legends. This is important not only to follow\n“publishing standards”, but also “E. indica” in Figure 1c as Eleusine indica and\nEragrostis indica are both applicable and probable. Perhaps some explanation on why\nthese other grass species were chosen might improve the text. It was surprising that\n“E. indica” was used in Figure 1c and then not in Figure 1d – perhaps this could be\nexplained.\nResponse:\nWe apologize for the oversight and confusion about species names. In response, we\nhave carefully revised the manuscript and figure legends to ensure that the full species\nnames are provided upon first mention in the main text and consistently included in\nall figure legends.\nSpecifically, we have clarified the identity of “E. indica” in Figure 1c as Eleusine\nindica, to avoid confusion with Eragrostis indica, which is also a plausible\ninterpretation. We apologize for the ambiguity in the original figure legend.\nThe presence of Eleusine indica in Figure 1c but not in Figure 1d reflects differences\nin sampling. In Figure 1c, our goal was to present a broad phylogenetic context\nacross Poaceae by including representative grass species; therefore, Eleusine indica\nwas included. However, Figure 1d focuses specifically on tetraploid paddy and\nupland weeds, along with cultivated crops, and for this reason Eleusine indica was not\nincluded in the comparison. We have now clarified as follows:\n“Using domain-based gene family quantification across 18 genomes, we observed\nsignificant lineage-specific contraction in biotic stress-responsive gene families,\nhighlighting differences between polyploid upland and paddy weeds as well as crops\n(Fig. 1d).”\nI believe that the text in line 249 should refer to Figure 2b not 2a and then on line 250,\nit is referring to 2a not 2b. If that is the case, perhaps consider restructuring the figure\nso the map is on the left and the genome coverage vs read maps and unrooted\nmaximum-likelihood phylogenetic tree reversed. Similarly, there are several places\nwhere Figure 2f is referred to (e.g. lines 300, 302,) where I suspect the authors were\nintending to highlight Figure 2e. Because of the huge number of subpanels and\nSupplementary Figures with subpanels, can the authors please double check that each\nand every cross reference is correct. Also, I suggest that the panel layout reflects the\norder in the text - Figure 3e and 3f are only referenced on line 461 while 3g on 429\nand 3h on 448… I understand that some of this is just making them visually appealing\nand to make the figures fit in the space, but with so many essential Supplementary\nFigures and subpanels, it makes it slightly more difficult jumping back and forth\nbetween the text and the supplements.\nResponse:\nWe sincerely apologize for these oversights, which indeed represent basic errors, and\nwe truly appreciate your patience in identifying them. We have carefully re-examined\nand verified all figure and supplementary figure references throughout the manuscript\nto ensure they correctly correspond to the intended content. Specifically, after\nreorganizing the layout of Figure 2, the references around Lines 249–250 have been\ncorrected to refer to the corresponding figures. Similarly, the incorrect references to\nFigure 2f in Lines 300 and 302 have been amended to correctly point to Figure 2e.\nFigure 2 | Population structure and demographic history of Digitaria species.\na, Geographic distribution of Digitaria accessions used in this study.\nGeographic distribution of Digitaria accessions used in this study. The main\ncollection region is displayed in detail, with an inset in the upper-left corner showing\nthe full geographic extent of sampling across China. Different species are denoted by\ndistinct symbols, while D. sanguinalis varieties are indicated by colored circles.\nRange boundaries of major taxa are shaded. b, Species identification based on\nmapping rates and coverage of resequencing reads aligned to D. sanguinalis genome\n(#YJ2023). c, An unrooted maximum-likelihood tree of D. sanguinalis, D. ciliaris, and\nD. bicornis. d, A rooted maximum-likelihood tree of D. sanguinalis based on genome-\nwide SNPs, using D. ciliaris as an outgroup. Varieties and putative admixed\npopulations are color-coded. e, Phenotypic variation in seed morphology (aspect\nratio and 1000-seed weight) across Digitaria species (sample sizes: D. ischaemum [n\n= 9], D. violascens [n = 16], D. bicornis [n = 5], D. ciliaris [n = 51], D. sanguinalis\n[n = 485]) and D. sanguinalis varieties (pubescens [n = 111], sanguinalis [n = 105],\nparvispicula [n = 78], glabra [n = 35]). Central dots represent medians, and\nwhiskers denote the interquartile range. Statistical significance was assessed by\ntwo-tailed t-tests: *, p < 0.05; **, p < 0.01; ***, p < 0.001. f, Demographic\ninference of D. sanguinalis varieties. Time points T1-T4 denote estimated divergence\nevents. Nₑ, effective population size; mya, million years ago.\nAdditionally, we have reorganized the layout of the main figures, including Figure 2\nand Figure 3, to better reflect the logical order in which they were discussed in the\ntext. For example, the map in Figure 2 is now placed on the left, and the genome\ncoverage plots and phylogenetic tree panels have been reordered accordingly. We\nhave also adjusted the order of subpanels in Figure 3 to follow the sequence in which\nthey were introduced in the main text.\nFigure 3 | Sympatric introgression from D. ciliaris to D. sanguinalis contributed to\ncold adaptation.\na, Maximum-likelihood phylogeny based on chloroplast genome SNPs. D. sanguinalis\nand its putative progenitors (#YJ2023, #DZ2, and #YZGJ2) are marked by asterisks.\nb, Bar plots summarizing species and varieties composition across major clades. c,\nDistribution of introgression intensity between sympatric D. ciliaris and D.\nsanguinalis population pairs. Red dots represent mean genome-wide proportions of\nintrogression (PIG, right Y-axis). The violin plots show the genome-wide f value\nd\ndistribution across 8,517, 9,002, 12,958, 14,703, 9,673, 9,640 and 11,143 sliding\nwindows for populations LF to ZJ, respectively. d, Enriched Pfam domains in\nintrogressed regions across sympatric populations. Detoxification-related domains\nare highlighted in green and heat-shock-related is shown in red. ABC, ABC\ntransporter; Myb-like, Myb-like DNA-binding domain; UDP/GTs, UDP-glucoronosyl\nand UDP-glucosyl transferase; P450, Cytochrome P450; C2, C2 domain; AP2,\nAP2 domain; Hsp20, Hsp20/alpha crystallin family; D-mannose, D-mannose binding\nlectin. e, Manhattan plot showing SNP associations with minimum temperature of the\ncoldest month (Bio6) using LFMM. Dashed horizontal lines denote significance\nthresholds (blue represents the FDR correction, adjusted p = 0.05; red represents the\nBonferroni correction, adjusted p = 0.05). Genes previously validated in Oryza sativa\nare labeled. f, Geographic distribution of candidate haplotypes of DsRZ2. The\nbackground map is shaded according to Bio6 variation. g, Frequency of DsRZ2\nhaplotypes in D. ciliaris and four D. sanguinalis varieties. Accessions with missing\ngenotypes in the focal region were excluded. GLA, D. sanguinalis var. glabra; PUB,\nD. sanguinalis var. pubescens; SAN, D. sanguinalis var. sanguinalis; PAR, D.\nsanguinalis var. parvispicula. h, Proportion of introgression f distributions on Chr4\nd\nin DZ, SJZ and SX (see details in Supplementary Fig. 24). i, Stress-responsive genes\nidentified within the 54.2-54.6 Mb region of Chr4.\nAs far as I can tell, Figure 2c is not referred to anywhere in the text. Nor is Figure 4f.\nNor is Figure 4g. I may have missed these, but I did not find them despite trying.\nMoreover, I did not check each and every cross reference; please ensure that\nEACH of the figures and presented data are referenced, discussed, and analysed. If not… there\nare several instances where the Supplementary Figures are covered by considerable\namounts of text or even get their own sections (e.g. Supplementary Figures 21, 23).\nI’d ask the authors to consider if there are places to rethink the figures or rearrange\nthem to make space for data that is essential to your story.\nResponse:\nWe are sorry for the oversight regarding the missing references to several figures in\nthe text. We fully acknowledge that it is both basic and avoidable error, and we deeply\nregret any confusion it may have caused.\nIn response, we have carefully reconsidered the arrangement of the supplementary\nfigures as well as the organization of Figure 4, to ensure that the presentation of\nvisual data aligns with the logical flow of the analysis.\nFigure 4 | Genomic basis of nicosulfuron resistance in D. sanguinalis\na, Differences in nicosulfuron resistance of five Digitaria species. Circles indicate\nmean GR₅₀ values (g a.i. ha-1), with error bars representing interquartile ranges. b,\nTemporal trends of nicosulfuron resistance in D. sanguinalis accessions sampled from\n2013 to 2023. Linear regression lines are shown with shaded 95% confidence\nintervals; p-values from one-tailed tests are indicated. c, Manhattan plot of GWAS for\nnicosulfuron resistance in D. sanguinalis. Dashed horizontal lines denote significance\nthresholds (yellow represents the FDR correction, adjusted p = 0.05; red represents the\nBonferroni correction, adjusted p = 0.05). SNPs near known NTSR-related gene\nfamilies are highlighted in green. Previously validated gens are highlighted in yellow.\nd, Allele frequency distribution of the SNP at Chr20: 34,032,644 among accessions\nwith high (HR), moderate (MR), low (LR), and susceptible (S) resistance. e, GR₅₀\nvalues of accessions with different genotypes at Chr20: 34,032,644 (CC: n = 126;\nCA: n = 3; AA: n = 6). Statistical significance was assessed by two-tailed t-tests (***,\np < 0.001; ****, p < 0.0001). f, Zoomed-in view of GWAS signals in Chr20: 34.0-\n34.1 Mb. Bottom panel shows local genomic signatures (left: f ; right: F ) across HR,\nd st\nMR, and LR groups. Gray rectangles represent 99% bootstrap quantiles, dark green\npoints denote windows with Z-score > 4 and yellow points denote windows with one-\ntailed p < 0.001 in randomization tests. g, Divergence time estimates of the\nintrogression scenario for D. sanguinalis under MSC-I model within a 100-kb window\nsurrounding Chr20: 34,032,644. The orange dotted arrow represents introgression\nfrom D. ciliaris to resistant D. sanguinalis population along with probability. h, An\nunrooted maximum-likelihood tree inferred from SNPs within a 100 kb window\nflanking Chr20: 34,032,644.\nIn addition, we have conducted a thorough review of the entire manuscript to ensure\nthat all figures, tables, and datasets have been correctly referenced, discussed, and\nanalyzed in the main text. Specifically, references to Figure 2c, Figure 4f, and Figure\n4g have now been added and integrated appropriately into the narrative.\nPlease make sure that in the results section discussing the candidate genes associated\nwith local adaptation, that all references and data are cross referenced. Line 444 talks\nabout expression analysis of DsRZ2 across “different tissues” but as far as I can tell\nyou haven’t discussed sampling different tissue types.\nResponse:\nThank you very much for your valuable comment, and we sincerely apologize for the\nconfusion caused by the phrasing in Line 444. The expression analysis of DsRZ2 was\nbased on RNA-seq data obtained from root, stem, old leaf and young leaf tissues of D.\nsanguinalis. These tissues were selected for the analysis to provide a comprehensive\nview of gene expression across different plant organs (Fig. R2). To avoid\nmisunderstanding, we have revised the text accordingly in the Results section.\nSpecifically, we have now clearly stated that the expression data were obtained from\nroot, stem, old leaf and young leaf tissues. The sentence has been updated as follows:\n“Expression analysis, based on RNA-seq data from root, stem, old leaf and young leaf\ntissues, confirmed ubiquitous expression of DsRZ2 with RPKM values ranging from\n4.40 to 7.24 (Supplementary Fig. 27a).”\nIn addition, we have carefully checked and cross-referenced all relevant data and\nreferences in this section to ensure clarity and consistency.\nFig. R2 | Expression pattern of DsRZ2 across old leaf, young leaf, root and stem.\nOleaf, old leaf; Yleaf, young leaf.\nPlease move the label the Y axes for Figure 4b from above to the more conventional\nnext to the axis location.\nResponse:\nSorry for the unconventional placement of the Y-axis label. In response, we have\nrevised the figure to reposition the Y-axis label to the standard location adjacent to the\naxis as follows. We appreciate the reviewer’s suggestion, which has improved the\nclarity and consistency of the figure presentation.\nFig. 4b | Temporal trends of nicosulfuron resistance in D. sanguinalis accessions\nsampled from 2013 to 2023. Linear regression lines are shown with shaded 95% confidence\nintervals; p-values from one-tailed tests are indicated.\nThere is a strange symbol on line 567 in the text describing allele frequency, please\ncorrect.\nResponse:\nSorry for the mistakes. The formatting error was caused by an unintended symbol on\nline 567 in the allele frequency description. This has now been corrected in the\nrevised manuscript. The sentence has been updated as follows:\n“Notably, the DsSOH1 variant was not fixed in HR accessions (allele frequency\naround 0.4), collectively suggesting recent adaptive introgression through\nintrogression (Fig. 4e).”\nReviewer #2 (Remarks to the Author):\nThe manuscript by Huang et al., entitled “The Digitaria genomes reveal local adaption\nand herbicide resistance mediated by introgression” presents an impressive dataset\nand array of analyses to study the genomic diversity, evolutionary history, and\nadaptation to herbicides in a pressing agricultural weed, Digitaria sanguinalis,\ncommonly known as crab grass. In producing very high quality reference genomes for\nnot only the focal species but its diploid progenitors, the authors provide a valuable\nresource for systematically classifying Digitaria accessions in agricultural landscapes.\nThe expansive analyses of an equally impressive population genomic dataset include\nmodel based inference of demographic history, population structure, introgression\ngene-environment association, bio-assays for herbicide tolerance and accompanying\nGWA— the analyses on the genetic architecture of non-target site resistance are\nespecially impressive and convincing. However, the motivation for several of the\nanalyses presented are difficult to follow, and the conclusions given throughout the\nmanuscript seem to overstep the results presented. Generally, some work to improve\nthe flow of writing, including explicitly stating the hypotheses being tested for each\nanalysis, would help readability of this elaborate work.\nResponse:\nWe appreciate your positive comments and constructive suggestions, which are very\nhelpful for improving our manuscript.\nLine 229-230, 233-234: It seems very tentative to suggest gene family expansions\nrelate to herbicide resistance when this analysis was at the genus level, and there is\nplenty of variation in herbicide use / resistance levels within and across species. To\nclaim this, one would want to compare lineages contrasting in a history of herbicide\nuse to see if this result stands up. Similarly, that UDPs/GTs are enriched on\nChromosome 3, which happens to also have the highest TE density, does not\nconvincingly show that TEs have driven an adaptive expansion of this gene family.\nResponse:\nThank you for this insightful comment, and we sorry for the confusion caused by the\ninitial phrasing of this section. We have realized that the causal relationship expressed\nhere requires more robust evidence to verify. While the observed expansion of\nherbicide-related gene families in Digitaria may have provided a genomic foundation\nfor adaptation under herbicide exposure, there is currently no direct evidence linking\nthis expansion to herbicide selection. To avoid overinterpretation, we have now\nremoved the statement in Results section “contributed to adaptation to herbicide\nselection” and revised the manuscript as:\n“These expansions may have provided a genomic foundation for adaptation under\nmodern weed management systems.”\nTransposable elements (TEs) have been implicated in rapid plant adaptation through\ntheir role in promoting gene copy number variation (Lisch, 2013). In our case, we\nobserved that Chr3 exhibits both the highest TE density and an enrichment of UDP-\nglucuronosyl/glucosyl transferase (UDP/GT). While this co-localization may suggest\na possible association, we concur that the current evidence is insufficient to establish a\ncorrelation between TE activity and the expansion of this gene family.\nTo avoid misunderstanding, we have now removed the statement in the Results\nsection “transposable elements may have contributed to the dynamic expansion” and\nrevised the manuscript as:\n“Notably, UDP/GTs were significantly enriched on Chr3 in Digitaria (Fisher’s test, p\n< 0.01) (Supplementary Fig. 10), which also harbors the highest transposon density\nacross the genome (Supplementary Fig. 7b).”\nReferences for this response:\nLisch, D., 2013. How important are transposons for plant evolution? Nat. Rev. Genet.\n14, 49–61.\nLines 348-349: It would be useful to provide more context to help the reader interpret\nthese results in the main text - what was the temporal scale of longitudinal\nsurveillance? How many years were sampled? Were the same localities revisited?\nEither way, I would want to make sure that geographic coordinates were controlled\nfor when conclusions are made about structure changing through time to be sure it's\nnot just driven by differences in space.\nResponse:\nWe thank you for this important point and sincerely apologize for the lack of clarity in\nthe original text. The longitudinal sampling was conducted across the Huang-Huai-\nHai region from 2013 to 2023. The majority of samples were collected in 2013 (72\naccessions), 2015 (83 accessions), and 2023 (285 accessions), with additional samples\ncollected in 2016, 2018, 2019, 2020, and 2021 in smaller numbers (Supplementary\nTable 7). In addition, these accessions were collected through field surveys conducted\nalong fixed routes that were revisited across different years. We have supplemented\nthe manuscript with information on the spatial distribution of samples across sampling\nyears in the revised manuscript (Fig. R1).\nFig. R1 | Distribution of Digitaria accessions collected from 2013 to 2023 across\nHuang-Huai-Hai region.\nTo control for geographic effects in the temporal comparison of population structure,\naccessions collected within a 100 km physical distance were clustered into the same\nlocality group. Groups with fewer than four samples sampled from 2013, 2015 or\n2023 were excluded from subsequent analyses to ensure robustness. Finally, 11\ngroups were selected for downstream analysis (Supplementary Table 11).\nWe have included this information in the revised manuscript, both in the Results and\nMethods sections, to provide clearer context for interpretation.\nLines 384: The authors claim early on that there was evidence for pervasive\nintrogression between sympatric populations, but more context in how these tests\nwere performed would really help with interpreting these results. Were these tests\nbetween the diploid, tetraploid, and hexaploid populations? Was there more\nintrogression in some directions compared to others? The authors note there is\nevidence for introgression in 7 tests, but how many total populations were tested?\nResponse:\nWe thank you for these helpful comments and agree that additional methodological\ncontext was necessary. In brief, we applied a genome-wide test for introgression\nbetween seven sympatric population pairs of hexaploid D. ciliaris and hexaploid D.\nsanguinalis. To define sympatry, we restricted population pairs to those collected\nwithin a 150 km geographic radius. Each pair included at least two D. ciliaris and\neight D. sanguinalis individuals to minimize sample size bias for statistical analysis\n(Supplementary Table 12). Among them, the strongest signal was detected in the SX\npopulation (D = 0.16), while the weakest was observed in the JN population (D =\n0.11). These clarifications and methodological details have been added to the revised\nmanuscript.\nLine 389: the acronym PIG is used before it is defined\nResponse:\nSorry for the mistakes. We have added the description in the main text:\n“This pattern was further supported by the same analysis using sliding-windows,\nrevealing that SX harbored the greatest proportion of introgression across the genome\n(PIG) across the genome (Fig. 3c; Supplementary Table 13).”\nLine 393-398: It is hard to interpret these results in the absence of a well defined null\nexpectation. I would like to see a test of whether correlations between PSIG and\ngeographic distance / environmental divergence are greater than expected from the\ngenome-wide background levels of divergence /population structure. The fact that\nenvironmental divergence did not remain a significant predictor after accounting for\ngeographic distance to me, suggests that this pattern is in part due to population\nstructure, as we know it to be strongly structured by geography. Therefore, I do not\nbelieve this is evidence for adaptive/environment mediated introgression.\nResponse:\nThank you for your thoughtful comment. We recognize that our original analysis\nlacked a clearly defined null expectation. In the revised manuscript, we have clarified\nthat the null hypothesis tested in the Mantel and partial Mantel analyses was that\nshared genomic regions of introgression between D. sanguinalis and D. ciliaris,\nmeasured by PSIG, were unrelated to either geographic distance or environmental\ndivergence.\nWe also agree that it is essential to account for background population structure,\nwhich may confound observed correlations. To evaluate whether environmental\ndivergence may be indirectly associated with PSIG via its correlation with genetic\ndivergence, we first computed pairwise F values among groups. Mantel tests\nst\nrevealed a significant association between F and environmental distance matrices\nst\n(Pearson's r = 0.5501, p = 0.0240). We then performed a partial Mantel test to assess\nthe relationship between PSIG and environmental distance while controlling for F .\nst\nThis analysis revealed a significant negative correlation, suggesting that\nenvironmental divergence explains variation in PSIG beyond what can be attributed to\ngenetic divergence (Table R1). Additionally, we tested the correlation between PSIG\nand specific environmental variables and found a significant association with\nprecipitation, further supporting a possible role of environmental factors in shaping\nintrogression patterns.\nPartial Mantel test control for\nenvironmental distance Partial Mantel test\nMantel test\n(Geographic distance) or control for Fst\nGeographic distance (Others)\ncorrelation p-value correlation p-value correlation p-value\nGeographic distance -0.700 0.024 -0.358 0.378 -0.286 0.255\nEnvironmental distance -0.644 0.007 -0.249 0.246 -0.549 0.019\nLatitude -0.380 0.133 -0.207 0.230 -0.264 0.209\nLongitude -0.641 0.077 0.020 0.392 -0.153 0.465\nPrecipitation -0.754 0.009 -0.148 0.345 -0.551 0.018\nSolar radiation -0.408 0.131 -0.167 0.324 -0.101 0.412\nWind speed -0.536 0.026 -0.108 0.378 -0.367 0.100\nTemperature -0.245 0.236 -0.083 0.397 -0.143 0.342\nTable R1 | Mantel regression and partial Mantel regression for the correlation\nbetween environmental distances and geographical distances among 7 groups of\nD. sanguinalis.\nWe have carefully reorganized this section of the manuscript, clarified the statistical\nframework, and explicitly accounted for genetic divergence in our analysis of the\ncorrelation between environmental distance and PSIG.\nLine 417-419: As there has been work focused on these gene families in NTSR, it\nseems a bit like an over simplification/cherry-picking to mention only these gene\nfamilies when there are supposedly 3500 genes introgressed. I also was not able to\nfind details in the methods about how these tests were performed. Generally, I find\nthese results to provide weak evidence of environment-driven selection shaping\nintrogression patterns, which is concluded in the final sentence. It seems like\noverlaying this work with the GEA results in the next section would first be necessary\nto say this. This is done for only one gene – DsRZ2, but I would like to see a more\nsystematic comparison of introgression regions with respect to introgressed regions to\nsupport this more general conclusion the authors make.\nResponse:\nThank you for your comments. Firstly, regarding the potential oversimplification in\nhighlighting only a subset of gene families, we acknowledge the limitation and have\nreorganized the corresponding section in the Results to provide a balanced summary\nof functional domains observed across introgressed regions: “Several functionally\nrelevant domains were recurrently enriched across populations, spanning categories\nsuch as detoxification (e. g., ABC transporters, UDP/GTs), signaling (e. g., protein\nkinases, WD40 repeats), and transcriptional regulation (e. g., Myb-like DNA-binding\ndomains) (Fig. 3d; Supplementary Table 15).”\nSecondly, we agree that our original conclusion may have overstated the evidence for\nenvironment-driven introgression. In response, we conducted a systematic comparison\nby overlaying the significant loci identified in the GEA analysis with introgressed\ngenomic regions. This revealed that environment-associated loci were not\nsignificantly enriched within introgressed regions in any population (Fisher’s exact\ntest, p = 0.999 for merged introgressed region; Table R2). This suggests that\nintrogression patterns may be influenced by factors beyond environmental selection,\nsuch as heterogeneous recombination landscapes or genetic divergence. Consequently,\nwe have revised the conclusion of this section to reflect this finding and removed\nspeculative language that could imply stronger evidence than supported by the data.\nRegion length Environment associated SNPs p-value\nJN 26250000 53 0.999994\nLF 32300000 120 0.500959\nHZ 30050000 55 1\nSX 26200000 69 0.99889\nDZ 26550000 77 0.988996\nSJZ 35750000 137 0.357647\nZJ 36100000 113 0.970588\nTable R2 | One-sided Fisher’s exact test results for the enrichment of\nenvironment-associated SNPs within introgressed genomic regions across\npopulations.\nFinally, in response to the request for methodological clarification, we have now\nadded a detailed description of the Pfam enrichment analysis in the Methods section.\nLine 436: Genomics islands typically refer to a debate centered around speciation and\nreproductive isolation, not adaptation. A more standard contrast would be polygenic\nadaptation vs. selective sweeps.\nResponse:\nThank you very much for pointing out this imprecise terminology. We have revised\n“genomic islands” in the context of local adaptation and corrected as follows:\n“These environment-associated variants exhibited a genome-wide distribution,\nwithout clustering in specific genomic regions, suggesting that adaptation in Digitaria\nis governed by polygenic architectures rather than by hard selective sweeps at a\nhandful of loci (Supplementary Fig. 25).”\nLine 450: Can you provide some context for interpreting the lack of significant\ndifference between the T and C haplotype?\nResponse:\nThank you for your valuable suggestion. We fully understand your concern about the\ninsufficient clarification regarding the lack of significant differences in the EHH\nresults.\nTo detect genes underlying heritable phenotypic variation, both top-down approaches\n(e. g., genome-wide association studies and genotype-environment association\nanalyses) and bottom-up approaches (e. g., population genomics) are commonly\nemployed. In this study, we adopted a bottom-up approach, extended haplotype\nhomozygosity (EHH), to detect adaptive signatures in DsRZ2. However, consistent\nwith previous findings in African indigenous chickens and Populus koreana (Gheyas\net al., 2021; Sang et al., 2022), we observed no significant haplotype differentiation in\nthe regions potentially under local adaptation. As noted by Rellstab et al. (2015),\nEHH-based methods have limited power in detecting subtle allele frequency shifts\nassociated with polygenic adaptation. Therefore, it is not unexpected that loci strongly\nassociated with phenotypes may lack clear haplotype differentiation signals. We have\nrevised the relevant sections in the manuscript to clarify this point more effectively.\n“Furthermore, extended haplotype homozygosity (EHH) analysis at the DsRZ2 locus\ndidn’t exhibit significant differences between haplotypes carrying the T or the C\nallele, which may result from polygenic additive effects (Supplementary Fig. 27;\nstandardized |iHS| score = 0.72).”\nReferences for this response:\nGheyas, A.A., Vallejo-Trujillo, A., Kebede, A., Lozano-Jaramillo, M., Dessie, T.,\nSmith, J., Hanotte, O., 2021. Integrated environmental and genomic analysis reveals\nthe drivers of local adaptation in African indigenous chickens. Mol. Biol. Evol. 38,\n4268–4285.\nSang, Y., Long, Z., Dan, X., Feng, J., Shi, T., Jia, C., Zhang, X., Lai, Q., Yang, G.,\nZhang, H., Xu, X., Liu, H., Jiang, Y., Ingvarsson, P.K., Liu, J., Mao, K., Wang, J.,\n2022. Genomic insights into local adaptation and future climate-induced vulnerability\nof a keystone forest tree in east Asia. Nat. Commun. 13, 6541.\nRellstab, C., Gugerli, F., Eckert, A.J., Hancock, A.M., Holderegger, R., 2015. A\npractical guide to environmental association analysis in landscape genomics. Mol.\nEcol. 24, 4348–4370.\nLine 473-474: How do we interpret survival rate vs control efficacy? I’m not sure I\nunderstand the difference.\nResponse:\nSurvival rate (%) measures the proportion of individual plants that physically\n“survived” the herbicide insult (binary dead / alive scoring). A survival rate of 100 %\nmeans every plant stayed alive; 0 % means all plants were completely necrotic.\nControl efficacy (%) quantifies how much total biomass was reduced by the treatment\nrelative to the non-treated control, irrespective of how many individuals lived or died.\nAn efficacy of 100 % indicates total biomass suppression (no living tissue), whereas\n0 % means biomass was identical to the control.\nThus, the survival rate reflects the number of surviving plants, whereas control\nefficacy reflects the amount of living tissue. A treatment could kill only half of the\nplants (50 % survival) yet strongly stunt the survivors, giving a high control efficacy\n(e. g., 85 %). We now explicitly state this distinction in the revised manuscript and\nhave added the explanatory sentences above to the Results section as follows:\n“Survival rates (number alive after treatment) following herbicide treatment rose\nfrom 36.5% in 2013 to 74.2% in 2023 (Supplementary Table 18), while the mean\nherbicide control efficacy (proportion of tissue successfully controlled) concurrently\ndeclined from 86.0% to 70.7%.”\nLine 497: There is no data presented to support the conclusion that resistance has\nevolved independently multiple times. That resistance accessions were widespread\ncould be the result of gene flow and recombination onto local backgrounds.\nResponse:\nThank you for this insightful comment. We agree that the original interpretation\nsuggesting multiple independent origins of resistance was not sufficiently supported.\nWe have removed the original claim and revised the text to instead emphasize the\npractical implications of the observed resistance pattern. The revised sentence now\nreads:\n“Notably, accessions exhibiting resistance (GR₉₀ > 60 g a.i. ha-1) were found to be\ngeographically widespread across multiple provinces, highlighting the increased\nchallenges such a distribution poses for effective weed management (Supplementary\nFig. 28).”\nLine 505: Are any of the mutations non-synonymous? This sentence could also use\ngrammatical edits, and the entire manuscript should be carefully edited for similar\ngrammatical errors.\nResponse:\nWe are sorry for the confusion caused by the original phrasing. Among the three\nDsALS gene copies, we identified a total of 25 non-synonymous mutations. We have\nrevised the sentence accordingly to clearly state this information. We also appreciate\nthe reviewer’s note regarding grammatical issues. The manuscript has been carefully\nedited throughout to improve grammar, clarity, and overall readability. This section\nhas been rewritten in the revised manuscript.\nLines 516-518: The authors ignore dominance in this discussion - a single dominant\nmutation present on only one genome-copy could confer resistance, even in a\nhexaploid, so this does not conclusively rule out target-site resistance mechanisms.\nResponse:\nWe appreciate the reviewer’s insightful comment and fully agree that the potential\nrole of dominant mutations was not adequately considered in our original phrasing.\nAlthough no fixed mutations in the three DsALS copies were detected, we recognize\nthat, given the hexaploid nature of D. sanguinalis, a dominant TSR mutation present\nin a single copy could still confer resistance. In response, we have revised the text to\nprovide a more cautious interpretation:\n“While no fixed mutations were identified in the three DsALS copies, the complexity\nof the hexaploid genome does not preclude the presence of a dominant mutation.\nFurther studies at both the genetic and functional levels will be necessary to establish\nwhether such a mutation underlies TSR in D. sanguinalis.”\nLine 523: I think in this case it would be more accurate to say gene amplification\nrather than expansions? I think that the GWAS is done based on coverage, so it would\nbe more descriptive/informative to you “assessed variation in coverage”\nrather than copy number. To help with readability more generally, I think this section\ncould be shortened given the lack of evidence for resistance associated with copy\nnumber variation.\nResponse:\nWe fully agree that \"gene amplification\" is a more accurate term than \"expansion\" in\nthis context. We also agree that “variation in coverage” better reflects the nature of the\ndata derived from the GWAS. Accordingly, we have revised the terminology in this\nsection and reworded the relevant sentences for clarity. In line with the reviewer’s\nrecommendation, we have also shortened the section to improve readability as\nfollows:\n“Gene amplification-mediated dosage effects represent a mechanism for rapid\nadaptation to intense selection pressures, as illustrated in Alopecurus myosuroides\n(blackgrass), and bamboo rats (Cai et al., 2023; Klure et al., 2025). To determine\nwhether NTSR in D. sanguinalis involves gene amplification, we assessed variation in\ncoverage for six key resistance-related gene families (cytochrome P450s, GSTs,\nUDP/GTs, AKRs, NB-ARC and ABC transporters). Most accessions showed coverage\nlevels for resistance-related gene families similar to those of single-copy genes,\nalthough 6 to 73 accessions exhibited high-coverage genes across different gene\nfamilies (Supplementary Fig. 29b). However, correlation analysis between resistance\nlevels and NTSR gene coverage revealed no association with resistance\n(Supplementary Fig. 29c; Supplementary Table 22). These findings do not support a\nmodel of metabolic resistance mediated by gene amplification in D. sanguinalis.”\nLine 553: What is a GWAS prioritized gene?\nResponse:\nWe are sorry for the confusion. In this context, we intended to refer to genes\nsignificantly associated with herbicide resistance identified by GWAS. We have\nrevised the corresponding statement in the manuscript for improved clarity, as\nfollows:\n“A total of 40 resistance-associated SNPs (rSNPs) were identified, corresponding to\n19 genes, including three homologs of genes functionally validated to confer non-\ntarget-site herbicide resistance in other grass species (Fig. 4c) (Chen et al., 2025;\nGoldberg-Cavalleri et al., 2023; Pan et al., 2022)”\nLine 565: Again, what was the source population of the fd peak at DsSOH1? Where\ndo we think this haplotype came from? How many tests gave such a result?\nResponse:\nWe apologize for the lack of a clearly defined null expectation in our analysis of the\norigin of DsSOH1, which may have contributed to ambiguity in the text.\nTo infer the potential origin of DsSOH1, we employed the ABBA-BABA test under\nthe topology (((P1, P2), P3), O), with D. bicornis as the outgroup (O), susceptible D.\nsanguinalis as P1, resistant D. sanguinalis as P2, and D. ciliaris as P3\n(Supplementary Fig. 31a). Regions introgressed from D. ciliaris into the resistant\ngroup were expected to display elevated f values. The prominent f peak at the\nd d\nDsSOH1 locus thus suggests that D. ciliaris may have contributed alleles to this\nregion. This inference is further supported by the phylogenetic tree constructed from\nthe 100-kb regions surrounding DsSOH1, where the topology is consistent with\nintrogression from D. ciliaris (Fig. 4h).\nHowever, since resistant accessions did not form a monophyletic clade, we\nacknowledge that this pattern could alternatively be explained by independent\nevolution or more complex gene flow events. We currently lack direct evidence to\ndistinguish these scenarios.\nTo assess whether different levels of resistance are associated with distinct patterns of\nintrogression, we conducted separate analyses for the Low Resistance (LR), Moderate\nResistance (MR), and High Resistance (HR) groups. Specifically, we included 10\naccessions for HR, 43 for MR, 27 for LR, and 55 for the susceptible group in our\ntests.\nTo clarify this issue, we have added a description of the null expectation in the revised\nmanuscript and provide a detailed representation of the null model in Supplementary\nFig. 32a.\n“To trace the origin of the haplotype at DsSOH1, we applied ABBA-BABA analysis\nusing the topology (((P1, P2), P3), O), with D. bicornis, susceptible D. sanguinalis,\nresistant D. sanguinalis, and D. ciliaris assigned to O, P1, P2, and P3, respectively\n(Supplementary Fig. 32a). A prominent f peak was detected specifically within the\nd\nDsSOH1 locus in the HR group, implying potential introgression from D. ciliaris,\nwhereas the MR and LR groups exhibited progressively weaker admixture signals\n(Fig. 4f).”\nLine 572/589: DsSOH1 is one of the numerous (~40) NTSR-associated loci you have\ndetected. Therefore it feels like an overstatement to conclude that introgression is a\ndriving force behind NTSR. Do any of the other loci show such a signal? Is it beyond\nwhat you would expect compared to the frequency of introgression across the rest of\nthe genome?\nResponse:\nWe agree that the introgression origin of DsSOH1, as one of the NTSR-associated\nSNPs (rSNPs), is not sufficient to generalize introgression as the primary driving\nforce of NTSR.\nTo assess whether introgression broadly contributes to NTSR, we compared f values,\nd\nused here to quantify the strength of introgression from D. ciliaris, within 50-kb\nwindows surrounding rSNPs with those from genome-wide background regions. As\nshown in Fig. R3, the HR group exhibited significantly elevated introgression\nfrequencies at rSNPs compared to the genomic background (Mann-Whitney U test, p\n= 0.02). This pattern was not observed in the MR or LR groups.\nNotably, the bimodal distribution of f values in the HR group suggests that while\nd\nintrogression may contribute to herbicide resistance, it is unlikely to be the only\nunderlying mechanism.\nFig. R3 | Distribution of f values in LR, MR, and HR groups across 50-kb\nd\nwindows surrounding rSNPs and genome-wide background regions. The dashed\nline indicates the 90th percentile f value derived from the genomic background.\nd\nStatistical significance was assessed using a one-tailed Mann-Whitney U test.\nTaken together, these results suggest a positive correlation between the frequency of\nintrogression from D. ciliaris and herbicide resistance levels in D. sanguinalis,\nparticularly in the most resistant populations. To avoid further misunderstanding, we\nrevised the statements at Line 572 as “Remarkably, the regions with the highest F\nST\ndifferentiation between resistant and susceptible populations coincided precisely with\nthe f peak interval, strongly implicating introgression from D. ciliaris at the DsSOH1\nd\nlocus (Fig. 4f).” and at Line 589 as “This recurrent phylogenetic pattern suggests\nthat the NTSR-associated resistant haplotype at DsSOH1 likely evolved through\nparallel and recent introgression from neighboring D. ciliaris populations.”\nWe further present the association between introgression frequency from D. ciliaris\nand herbicide resistance levels in D. sanguinalis in the next section.\nLine 595-598: Unclear how these two analyses relate to one another\nResponse:\nWe appreciate your comment and agree that the relationship between the two analyses\nwas not clearly articulated in the original version. To improve clarity, we have revised\nthis section as follows:\n“To investigate the correlation between introgression frequency from D. ciliaris and\nherbicide resistance levels in D. sanguinalis, we compared the f values at rSNPs to\nd\nthe genome-wide background. The HR group exhibited significantly higher\nintrogression frequencies at rSNPs compared to the genomic background (Mann-\nWhitney U test, p = 0.02), a pattern not observed in the MR or LR groups\n(Supplementary Fig. 31c). In addition, the HR group exhibited the greatest\ncumulative length of introgressed windows across the genome (HR vs. MR, Mann-\nWhitney U test, p = 1.6e-2; HR vs. LR, Mann-Whitney U test, p = 1.18e-4)\n(Supplementary Fig. 31d). These findings suggest that introgression from D. ciliaris\nmay be associated with the development of NTSR in D. sanguinalis. Interestingly, we\nobserved an increase in the frequency of rSNPs in accessions collected in 2023\ncompared to those collected in 2013 (Supplementary Fig. 31e), which may be linked\nto the higher herbicide resistance observed in 2023.”\nLine 681-684: I do not follow this sentence for both conceptual and grammatical\nreasons.\nResponse:\nWe agree that the original phrasing was grammatically problematic and conceptually\ncondensed. Our intention was to convey that introgression within the Digitaria\ncomplex may contribute to local adaptation across diverse environments. We have\nrevised the sentence for clarity and precision as follows:\n\"Although introgression is increasingly recognized as a potential mechanism of\nadaptive evolution, instances of beneficial introgression remain limited (Edelman et\nal., 2019; Taylor and Larson, 2019). In our study, we identified recurrent post-\ndivergence introgression between D. ciliaris and D. sanguinalis, two species with\noverlapping distributions across ecologically heterogeneous regions. This frequent\ngene flow appears to have contributed to environmental adaptation in D. sanguinalis,\nparticularly with respect to cold tolerance (Fig. 3h). Similar patterns of adaptive\nintrogression have also been reported in oaks (Quercus spp.), where introgression\nbetween sympatric species may contribute to climatic adaptation (Fu et al., 2022;\nRen et al., 2025). Given their widespread geographic overlap, such processes may be\nprevalent within globally distributed weedy genera, such as Echinochloa\n(barnyardgrass) and Setaria (green foxtail), where recurrent hybridization and\nselection may drive local adaptation in response to heterogeneous environments\nworldwide.\"\nLine 687: I do not think this novelty statement is necessary or true. Adaptive\nintrogression of herbicide resistance has been shown in Teosinte from Maize\n(https://pubmed.ncbi.nlm.nih.gov/32989136/) and in numerous other systems in\nresponse to pesticides more broadly.\nResponse:\nWe sincerely apologize for the oversight in our literature citation and phrasing. We\nfully acknowledge that adaptive introgression of herbicide resistance has been\npreviously reported, including the well-documented case of introgression from maize\nto teosinte (Le Corre et al., 2020), as well as in other systems such as weedy rice\n(Wedger et al., 2024, 2022) and Brassica (Pandolfo et al., 2016).\nThese studies have primarily focused on crop-to-weed gene flow. Our study highlights\na different scenario in which adaptive introgression occurs between two weedy\nspecies within the Digitaria complex. This weed-to-weed gene flow may represent an\nalternative and underappreciated route through which herbicide resistance can evolve\nand spread. Furthermore, our study integrates extensive phenotypic screening of 196\nrepresentative accessions using dose-response assays, together with genome-wide\nresequencing, allowing us to identify resistance-associated SNPs and provide genomic\nevidence supporting their introgressive origin. We believe this combined approach\nadds depth to the current understanding of how herbicide resistance can emerge via\ninterspecific introgression in weedy taxa.\nWe have revised the novelty statement to more accurately reflect the contribution of\nthis work within the context of existing literature, as follows:\n“While introgression events in weeds have been increasingly documented (Ribeiro et\nal., 2025; Wedger et al., 2024), our findings suggest that adaptive introgression\nbetween closely related weedy species may contribute to the rapid evolution of\nherbicide resistance.”\nReferences for this response:\nLe Corre, V., Siol, M., Vigouroux, Y., Tenaillon, M.I., Délye, C., 2020. Adaptive\nintrogression from maize has facilitated the establishment of teosinte as a noxious\nweed in Europe. Proc. Natl. Acad. Sci. 117, 25618–25627.\nWedger, M.J., Roma-Burgos, N., Olsen, K.M., 2022. Genomic revolution of US\nweedy rice in response to 21st century agricultural technologies. Commun. Biol. 5,\n885.\nWedger, M.J., Xiao, E., Butts, T.R., Chlapecka, J.L., Webster, L.C., Olsen, K.M.,\n2024. Recent crop-to-weed adaptive introgression has reshaped the genomic\ncomposition and geographical structure of US weedy rice (Oryza spp.). Mol. Ecol.\ne17604.\nPandolfo, C.E., Presotto, A., Carbonell, F.T., Ureta, S., Poverene, M., Cantamutto, M.,\n2016. Transgenic glyphosate-resistant oilseed rape (Brassica napus) as an invasive\nweed in Argentina: detection, characterization, and control alternatives. Environ. Sci.\nPollut. Res. 23, 24081–24091."}, {"round": 2, "reviewer_comments": {"reviewer_1": "This version of the manuscript still presents a very large and diverse dataset, including\nphenotyping and then genome assembly and annotation of D. sanguinalis and its\nprogenitors, population genomics of 579 accessions (covering structure, diversity, and\ndemographic analyses), evidence of introgression and environmental adaptation,\nherbicide resistance phenotyping with dose–response curves and spatiotemporal\ntrends, comparative analysis of target-site versus non-target-site resistance, GWAS for\nNTSR loci, transcriptomic profiling of resistant versus susceptible biotypes, and\ncorrelations between introgression and resistance. The manuscript’s discussion is\ncomparatively brief, focusing on evolutionary and management implications.\nThe manuscript is supported by extensive resources: 4 multipaneled figures, 38(?)\nsupplementary figures, 24 supplementary tables, 3 supplementary notes, and 4\nsupplementary datasets. Despite its complexity, the text is generally well-structured\nand readable, though challenging for any single reviewer to evaluate in full given the\nbreadth of the science that is presented. The integrated approach provides a\ncomprehensive resource; I assume the authors have made a considered decision to\nmaintain this format rather than divide the work into separate papers.\nI thank the authors for addressing many of my previous suggestions and for providing\na clear rebuttal outlining the changes made. Improvements include a more logical\nintroduction, smoother transitions, corrected terminology, revised abbreviations, and\nclearer figure legends. These changes enhance readability and scientific clarity.\nAreas for further attention:\nFigure referencing: I pointed out previously that the figure references were incorrect\nand sometimes the narrative did not match panel order. Because of the complexity of\nthe datasets and the edits that were made in the revision and time allocated to do the\nrevision, I have not systematically checked that the associations between text and\nfigures are all correct. I hope instead that authors have been more diligent in this\nsecond revision and expect that this will be double checked in the production process\nif the paper progresses.\nAbstract: Now more comprehensive but still omits phenotyping, which according to\nthe text is novel and central to Figures 1a, 1c, and 2b.\nClarity: A few sentences remain unclear or incomplete. Please revisit these for\nprecision.\n• 165 = “sequencing reads could be successfully mapped, with properly paired ____\n(Supplementary Table 3).”\n• 222 – 225 = “Flood-adapted genes, Echinochloa displayed marked expansion of …”\n• 554 = “…collectively suggesting recent adaptive introgression through introgression…”\nClarity: Sentence at 553-554 “allele frequency around 0.4…” lacks the strange symbol that was in the previous version, but\nwhat do you mean “around 0.4”? Is “around” equivalent to “approximately” and if so, why not simply use the appropriate\nlevel of precision instead of stating “around”?\nInterpretation: The text at lines 367-369 was rewritten to emphasize that interspecific\nintrogression events are more probable than hybridization and it now cites monocot\ngrafting incompatibility literature. Those additions, although useful do not explain the\nrationale for introgression based on chloroplast phylogenies – perhaps because this is\nnot my speciality. Is introgression the only explanation for the observation? If so,\nstate this explicitly.\nInterpretation: Although your data do not support a model of metabolic resistance mediated by gene amplification (lines 514-\n515) do they exclude it? If not, please refine this sentence – lack of evidence is not evidence of lack. Particularly as you then\ngo on to explore NTSR mediated by e.g. CYP, AKR and ABC genes in the sections following line 532 and emphasize that\n“NTSR has become the dominant resistance strategy” (line 644) in the Discussion – which remains quite brief compared to\nthe body of data presented.\nFigures: Newly added figures (e.g., R1, R2, R3) appear strikingly similar to existing\nsupplementary figures. Please clarify whether these replace previous versions or\nrepresent distinct data.\n• How does the added figure Fig. R1 with the year by year collections differ from\nwhat was already sent as Supplementary Figure 11? They seem to be the same except\nthat R1 is slightly more spaced out than Sup 11?\n• The text at 430–431 now specifies root, stem, old leaf, young leaf; expression range\nreported; Supp. Fig. 27a & Fig. R2 show profiles. However, how does Fig. R2 differ\nfrom the data shown in Supplementary Fig. 27a except that now there are the three\ndatapoints indicated in R2? Are these the same data? If so, it is inappropriate to show\nthem twice. Will R2 replace Sup Fig 27a?\n• Similarly, will R3 replace Supplementary Figure 32c? For all three of the\nSupplementary Figures added in the revision, it is not clear how they differ from what\nwas already presented except perhaps minor cosmetic changes. Please explain.\nLanguage: Some informal phrasing (e.g., “didn’t” at line 438) remains and should be\nrevised for scientific tone.\nMinor Issues:\n• “Aldo/keto reductases” should be “Aldo-keto reductases”\n• The text at line 247 has CH, DH, and EH without subscripting.\n• The authors alternate between using common names followed by bracketed Latin\nnames (e.g. “oaks (Quercus spp.”) on line 651 and then “Echinochloa (barnyardgrass) and Seteria (green foxtail)” on line 654). Choosing one way to state the\nnomenclature or the other would make the manuscript more appropriate for the\nscientific literature.", "reviewer_2": "I appreciate the effort the authors put into reviewing this manuscript, updating the\nconclusions to better reflect the analyses run and evidence provided. I maintain this is\nan interesting manuscript full of useful data and insights to the field. In terms of\npresentation, I do find this version easier to follow--however, there is still quite a lot\nof room for improvement throughout. While I will give a few examples of recurring\ncommunication issues, I encourage the authors to take a very careful look at the\nspelling, grammar, and narrative arc throughout.\nCommunication Issues:\n- Connection between subsections: When reading the manuscript, sections feel very\ndisjunct from each other, and only start to come together near the end and with the\nsynthesis in the discussion. If the authors could incorporate some more\ntransitional sentences throughout, and weave in why it was first important to understand\n(for e.g.) demographic history & pop structure before studying introgression etc, it\nwould strengthen the cohesiveness of the paper. The same is true for paragraphs of the\nintroduction -- as written, they seem almost unrelated, although clearly that is not true.\n- Motivation for analyses: not completely unrelated to the above point, but I still feel\nthere is opportunity to help the reader better understand the reasoning and progression\nfrom analysis to analysis - given how many analyses the authors present, this is really\nimportant. For e.g. in the Section on Resistance, the authors have subsections on the\ntemporal/spatial distribution of phenotypic resistance, TSR and NTSR variation, and\nthen Transcriptomic analysis. While experts in the field would intuit why all of these\nwere done, motivation should be given at the beginning of each section why and how\nthe authors expect these approaches to be informative about Resistance evolution.\nGrammatical and spelling issues throughout: There remain numerous grammatical\nerrors that affect readability and professionalism. This is not an exhaustive list, but\nexamples include:\n--> Title: \"local adaption\" should be \"local adaptation\"\n-->Line 17 (Highlights): \"Introgression from sympatric relatives have contributed\" - subject-verb disagreement (should be\n\"has contributed\")\n-->Introduction: \"Polyploidy is across the plant kingdom\" - missing word (should be \"Polyploidy is common across\" or\n\"Polyploidy is found across\")\n-->Results, line ~50: \"we significant lineage-specific contraction\" - missing verb (should be \"we observed significant\")\n-->Not a gramatical error but the section header “Population structure restructuration in ten years” reads very awkwardly.\n-->Multiple instances of missing spaces after punctuation, incomplete sentences, and\nwords missing letters.\n-->Some terminology is used inconsistently (e.g., \"NTSR-associated SNPs\" vs \"resistance-associated SNPs\" vs \"rSNPs\")\n-->Please look to see where run-on sentences can be split for clarity.\n- Figure legends: Several figure legends could be more informative and self-contained:\n-->All abbreviations should be defined in each legend\n-->Statistical tests used should be mentioned in the legend when significance is\nindicated\n-->The legends should allow readers to understand the main takeaway without having\nto refer back to the main text\n- Data availability: While the authors mention data will be made available, specific\ndetails are still missing:\n-->Accession numbers for genome assemblies\n-->Repository information for resequencing data\n-->Code availability for custom analyses should be included in a Data Availability\nstatement."}, "author_response": "A Point-to-point response to the reviewers’ comments\nReviewers’ comments are shown as black text. Response to comments is shown as\nblue text.\nReviewer #1 (Remarks to the Author):\nThis version of the manuscript still presents a very large and diverse dataset, including\nphenotyping and then genome assembly and annotation of D. sanguinalis and its\nprogenitors, population genomics of 579 accessions (covering structure, diversity, and\ndemographic analyses), evidence of introgression and environmental adaptation,\nherbicide resistance phenotyping with dose–response curves and spatiotemporal\ntrends, comparative analysis of target-site versus non-target-site resistance, GWAS for\nNTSR loci, transcriptomic profiling of resistant versus susceptible biotypes, and\ncorrelations between introgression and resistance. The manuscript’s discussion is\ncomparatively brief, focusing on evolutionary and management implications.\nThe manuscript is supported by extensive resources: 4 multipaneled figures, 38(?)\nsupplementary figures, 24 supplementary tables, 3 supplementary notes, and 4\nsupplementary datasets. Despite its complexity, the text is generally well-structured\nand readable, though challenging for any single reviewer to evaluate in full given the\nbreadth of the science that is presented. The integrated approach provides a\ncomprehensive resource; I assume the authors have made a considered decision to\nmaintain this format rather than divide the work into separate papers.\nI thank the authors for addressing many of my previous suggestions and for providing\na clear rebuttal outlining the changes made. Improvements include a more logical\nintroduction, smoother transitions, corrected terminology, revised abbreviations, and\nclearer figure legends. These changes enhance readability and scientific clarity.\nResponse:\nWe are grateful for your constructive comments/suggestions, which help us\nsubstantially improve our manuscript.\nAreas for further attention:\nFigure referencing: I pointed out previously that the figure references were incorrect\nand sometimes the narrative did not match panel order. Because of the complexity of\nthe datasets and the edits that were made in the revision and time allocated to do the\nrevision, I have not systematically checked that the associations between text and\nfigures are all correct. I hope instead that authors have been more diligent in this\nsecond revision and expect that this will be double checked in the production process\nif the paper progresses.\nResponse:\nWe thank the reviewer's reminder regarding figure referencing and consistency\nbetween the narrative and figure panels. In this revised version, we have conducted a\ncareful and thorough review of all figure references to ensure accuracy and\nconsistency.\nAbstract: Now more comprehensive but still omits phenotyping, which according to\nthe text is novel and central to Figures 1a, 1c, and 2b.\nResponse:\nIn the revised Abstract, we have now explicitly included a summary of classification,\nin which we provided comprehensive and accurate phenotypic categorization of\nDigitaria accessions. This classification was based on integrated morphological\nassessments and genome-scale evidence, and it underpins the genotype-phenotype\nassociations presented in the main text. We agree this addition enhances the\ncompleteness of the Abstract and better reflects the scope and novelty of our work.\nClarity: A few sentences remain unclear or incomplete. Please revisit these for\nprecision.\n• 165 = “sequencing reads could be successfully mapped, with properly paired ____\n(Supplementary Table 3).”\nResponse:\nWe thank the reviewer for pointing this out. The sentence has now been revised for\nclarity. It now reads:\n“Sequencing reads from diploid, tetraploid, and hexaploid Digitaria accessions were\nuniformly mapped to the D. sanguinalis reference genome, with high mapping rates of\n99.51%, 99.48%, and 99.58%, respectively (Supplementary Table 3).”\n• 222 – 225 = “Flood-adapted genes, Echinochloa displayed marked expansion of …”\nResponse:\nThis sentence has been restructured for grammatical correctness and improved\nreadability. It now reads:\n“Several flood-adapted gene families, such as the anaerobic-response Apetala2,\nshowed marked expansion in Echinochloa (copy number: 182.0 ± 12.3 vs. 159.0 ± 9.9\nin D. sanguinalis; two-tailed t-test, p = 0.008), likely reflecting adaptation to\nprolonged submergence (Supplementary Data 1).”\n• 554 = “…collectively suggesting recent adaptive introgression through\nintrogression…”\nResponse:\nWe have revised it for clarity and to avoid repetition. The sentence now reads:\n“Notably, the DsSOH1 variant was not fixed in HR group, with an observed allele\nfrequency of 0.40, supporting a scenario of recent adaptive introgression (Fig. 4d).”\nClarity: Sentence at 553-554 “allele frequency around 0.4…” lacks the strange\nsymbol that was in the previous version, but what do you mean “around 0.4”? Is\n“around” equivalent to “approximately” and if so, why not simply use the appropriate\nlevel of precision instead of stating “around”?\nResponse:\nWe thank the reviewer for the helpful suggestion and have now revised the sentence\nto improve clarity and precision. Specifically, we replaced “around 0.4” with the exact\nvalue “0.40” (as calculated from the genotype calls). The revised sentence reads:\n“Notably, the DsSOH1 variant was not fixed in HR group, with an observed allele\nfrequency of 0.40, supporting a scenario of recent adaptive introgression (Fig. 4d).”\nInterpretation: The text at lines 367-369 was rewritten to emphasize that interspecific\nintrogression events are more probable than hybridization and it now cites monocot\ngrafting incompatibility literature. Those additions, although useful do not explain the\nrationale for introgression based on chloroplast phylogenies – perhaps because this is\nnot my speciality. Is introgression the only explanation for the observation? If so,\nstate this explicitly.\nResponse:\nWe apologize for the confusion. In light of the chloroplast phylogenies, interspecific\nintrogression remains the most plausible and well-supported explanation for the\nobserved patterns. While alternative possibilities such as incomplete lineage sorting or\nsample misidentification were considered, these are less consistent with the combined\nevidence from bootstrap support, nuclear-chloroplast discordance, and morphological\ncharacters. We also acknowledge the theoretical possibility of natural grafting as a\nsource of chloroplast capture; however, this phenomenon is extremely rare in\nmonocots, especially grasses. Therefore, interspecific introgression-likely mediated\nby historical hybridization events-is the most likely explanation, as has been similarly\ndemonstrated in genera such as Secale and Malus (Liu et al., 2022; Sun et al., 2022).\nWe have revised the text in the manuscript accordingly to state this interpretation\nmore clearly.\n“Overall, the chloroplast phylogenies indicate extensive interspecific introgression\nacross Digitaria, most likely resulting from historical hybridization events (Feng et\nal., 2024; Gaut et al., 2019).”\nReferences for this response:\nLiu, B.-B., Ren, C., Kwak, M., Hodel, R.G.J., Xu, C., He, J., Zhou, W.-B., Huang, C.-\nH., Ma, H., Qian, G.-Z., Hong, D.-Y., Wen, J., 2022. Phylogenomic conflict\nanalyses in the Apple genus malus s.l. reveal widespread hybridization and\nallopolyploidy driving diversification, with insights into the complex\nbiogeographic history in the northern hemisphere. J. Integr. Plant Biol. 64,\n1020–1043.\nSun, Y., Shen, E., Hu, Y., Wu, D., Feng, Y., Lao, S., Dong, C., Du, T., Hua, W., Ye, C.-\nY., Zhu, J., Zhu, Q.-H., Cai, D., Skuza, L., Qiu, J., Fan, L., 2022. Population\ngenomic analysis reveals domestication of cultivated rye from weedy rye. Mol.\nPlant 15, 552–561.\nInterpretation: Although your data do not support a model of metabolic resistance\nmediated by gene amplification (lines 514-515) do they exclude it? If not, please\nrefine this sentence – lack of evidence is not evidence of lack. Particularly as you then\ngo on to explore NTSR mediated by e.g. CYP, AKR and ABC genes in the sections\nfollowing line 532 and emphasize that “NTSR has become the dominant resistance\nstrategy” (line 644) in the Discussion – which remains quite brief compared to the\nbody of data presented.\nResponse:\nThank you for this important clarification. We agree that the original sentence may\nhave been overly definitive. We have revised the sentence to reflect that while our\ndata do not provide evidence supporting metabolic resistance mediated by gene\namplification, and, however, we cannot exclude this possibility entirely. The revised\nsentence now reads as follow:\n“Overall, no evidence was found to support metabolic resistance mediated by gene\namplification in D. sanguinalis, though further investigation is needed to fully exclude\nthis possibility.”\nFigures: Newly added figures (e.g., R1, R2, R3) appear strikingly similar to existing\nsupplementary figures. Please clarify whether these replace previous versions or\nrepresent distinct data.\nResponse:\nFigures R1, R2, and R3 were initially included in the point-by-point response to\nprovide clearer, high-resolution versions of existing analyses for ease of review. Upon\nfurther consideration, we found that these figures provide valuable visual\nenhancements and clarification of the corresponding results. Therefore, we have\nintegrated them into the revised supplementary materials as additional panels to\ncomplement the original figures: R1 is now Supplementary Figure 11, R2 has been\nadded as Supplementary Figure 27a, and R3 as Supplementary Figure 32c. These\ndo not replace prior figures but serve as refined and informative extensions of the\nexisting content.\n• How does the added figure Fig. R1 with the year by year collections differ from\nwhat was already sent as Supplementary Figure 11? They seem to be the same except\nthat R1 is slightly more spaced out than Sup 11?\nResponse:\nSorry for the confusion. Figure R1 was originally provided in the point-by-point\nresponse as a high-resolution version of the revised collection timeline. This same\nfigure has now been formally included in the revised Supplementary Materials as\nSupplementary Figure 11. The current Supplementary Figures (including\nSupplementary Figure 11) represent the finalized versions for submission. Figures\nR1, R2, and R3 were included solely to facilitate reviewer evaluation and will not be\nsubmitted as separate figures in the final manuscript.\n• The text at 430–431 now specifies root, stem, old leaf, young leaf; expression range\nreported; Supp. Fig. 27a & Fig. R2 show profiles. However, how does Fig. R2 differ\nfrom the data shown in Supplementary Fig. 27a except that now there are the three\ndatapoints indicated in R2? Are these the same data? If so, it is inappropriate to show\nthem twice. Will R2 replace Sup Fig 27a?\nResponse:\nWe appreciate your careful check. To avoid redundancy, Figure R2 will not be\nsubmitted as a separate figure in the final manuscript. Instead, it has been fully\nintegrated as the Supplementary Figure 27a.\n• Similarly, will R3 replace Supplementary Figure 32c? For all three of the\nSupplementary Figures added in the revision, it is not clear how they differ from what\nwas already presented except perhaps minor cosmetic changes. Please explain.\nResponse:\nWe confirm that Figure R3 is not an additional figure but has been incorporated into\nthe revised Supplementary Materials as Supplementary Figure 32c, replacing the\nprevious version. No duplicate figures will be submitted.\nLanguage: Some informal phrasing (e.g., “didn’t” at line 438) remains and should be\nrevised for scientific tone.\nResponse:\nThank you for careful attention to language. We have now revised all instances of\ninformal or colloquial phrasing (e.g., “didn’t”) throughout the manuscript to ensure a\nconsistent and professional scientific tone.\nMinor Issues:\n• “Aldo/keto reductases” should be “Aldo-keto reductases”\nResponse:\nSorry for the mistakes and we have made corrections in the main text.\n• The text at line 247 has CH, DH, and EH without subscripting.\nResponse:\nSorry for the mistakes. The formatting has been corrected and properly subscripted in\nthe revised manuscript.\n• The authors alternate between using common names followed by bracketed Latin\nnames (e.g. “oaks (Quercus spp.”) on line 651 and then “Echinochloa (barnyardgrass)\nand Seteria (green foxtail)” on line 654). Choosing one way to state the nomenclature\nor the other would make the manuscript more appropriate for the scientific literature.\nResponse:\nThank you for pointing this out. We have revised the manuscript to ensure consistency\nin formatting, adopting the common name (Genus name) style throughout the text\n(e.g., line 506).\nIn line 654, to improve clarity and avoid the mixed use of common and scientific\nnames, we have modified the phrasing to:“Given the broad geographic overlap of\nmany weedy species, such processes may also be widespread in globally distributed\ngenera, such as Echinochloa (including E. crus-galli) and Setaria (including S.\nitalica), where recurrent hybridization and local selection pressures likely drive\nadaptation in response to heterogeneous environments worldwide.” This revision was\nmade to explicitly refer to representative weedy species within each genus, while\nmaintaining formal scientific nomenclature.\nReviewer #2 (Remarks to the Author):\nI appreciate the effort the authors put into reviewing this manuscript, updating the\nconclusions to better reflect the analyses run and evidence provided. I maintain this is\nan interesting manuscript full of useful data and insights to the field. In terms of\npresentation, I do find this version easier to follow--however, there is still quite a lot\nof room for improvement throughout. While I will give a few examples of recurring\ncommunication issues, I encourage the authors to take a very careful look at the\nspelling, grammar, and narrative arc throughout.\nResponse:\nWe appreciate your time and constructive comments and have tried our best to address\nthem below.\nCommunication Issues:\n- Connection between subsections: When reading the manuscript, sections feel very\ndisjunct from each other, and only start to come together near the end and with the\nsynthesis in the discussion. If the authors could incorporate some more\ntransitional sentences throughout, and weave in why it was first important to understand\n(for e.g.) demographic history & pop structure before studying introgression etc, it\nwould strengthen the cohesiveness of the paper. The same is true for paragraphs of the\nintroduction -- as written, they seem almost unrelated, although clearly that is not true.\nResponse:\nThank your insightful and constructive suggestion. We fully agree that improving the\nlogical flow and strengthening the connection between sections can significantly\nenhance the readability and cohesiveness of the manuscript.\nIn response, we have revised the Introduction by adding transitional sentences to\nclarify the conceptual progression from agronomic significance of weediness to its\nunderlying evolutionary mechanisms. For example, in the section introducing\nherbicide resistance, the opening sentence has been rephrased to emphasize its central\nrole in modern weed management: “With widespread and intensifying use of\nagrochemicals, herbicide resistance has emerged as a global threat to effective weed\ncontrol, currently affecting 273 weed species worldwide (Délye, 2013; Heap, 2025;\nPowles and Yu, 2010).” We also refined the Results sections by including brief\nforward-looking statements or transitional phrases at the beginning and end of\nrelevant subsections.\n- Motivation for analyses: not completely unrelated to the above point, but I still feel\nthere is opportunity to help the reader better understand the reasoning and progression\nfrom analysis to analysis - given how many analyses the authors present, this is really\nimportant. For e.g. in the Section on Resistance, the authors have subsections on the\ntemporal/spatial distribution of phenotypic resistance, TSR and NTSR variation, and\nthen Transcriptomic analysis. While experts in the field would intuit why all of these\nwere done, motivation should be given at the beginning of each section why and how\nthe authors expect these approaches to be informative about Resistance evolution.\nResponse:\nThank you for this valuable suggestion regarding analytical motivation and narrative\nclarity. We agree that it is important to clearly articulate the reasoning behind each\nanalytical step to guide the reader through the logic of the study.\nIn response, we have revised the manuscript to include brief motivating sentences at\nthe beginning of major analytical sections, particularly within the Resistance section.\nFor example, in the section on TSR and NTSR variation, we added the following\nsentence to clarify the purpose of the analysis: “To assess potential TSR, we first\nexamined non-synonymous mutations in ALS, the herbicide target gene of\nnicosulfuron.” We believe these revisions help make the analytical flow more\ntransparent and strengthen the overall readability of the manuscript for both\nspecialists and broader audiences.\nGrammatical and spelling issues throughout: There remain numerous grammatical\nerrors that affect readability and professionalism. This is not an exhaustive list, but\nexamples include:\nResponse:\nWe apologize for the remaining errors in the previous version. The manuscript has\nbeen carefully revised to correct all identified grammatical and spelling errors, and we\nhave made every effort to ensure that the revised version meets the standards of\nclarity and professionalism.\n--> Title: \"local adaption\" should be \"local adaptation\"\nResponse:\nSorry for the mistakes and we have corrected it.\n-->Line 17 (Highlights): \"Introgression from sympatric relatives have contributed\" -\nsubject-verb disagreement (should be \"has contributed\")\nResponse:\nWe have double-checked the current version of the sentence, which now reads (Line\n33): “Introgression from sympatric relatives has contributed to ecological adaptation\nand rapid evolution of non-target-site resistance in D. sanguinalis.” The subject-verb\nagreement is correct in this revised sentence.\n-->Introduction: \"Polyploidy is across the plant kingdom\" - missing word (should be\n\"Polyploidy is common across\" or \"Polyploidy is found across\")\nResponse:\nThanks for your suggestion. We have checked the current version of the manuscript,\nand the sentence in the Introduction now reads: “Polyploidization is widespread\nacross the plant kingdom and particularly prevalent among weedy species, where it\noften contributes to increased resistance to stress (Jiao et al., 2011; Rutland et al.,\n2021; Soltis et al., 2015).” We appreciate your attention to clarity and language\naccuracy.\n-->Results, line ~50: \"we significant lineage-specific contraction\" - missing verb\n(should be \"we observed significant\")\nResponse:\nThe sentence has been revised for clarity and now reads: “we observed significant\nlineage-specific contraction in biotic stress-responsive gene families”. We appreciate\nyour careful reading and helpful suggestion.\n-->Not a gramatical error but the section header “Population structure restructuration\nin ten years” reads very awkwardly.\nResponse:\nWe agree that the original subheading was awkwardly phrased. To improve clarity and\nreadability, we have revised the section title to: “Spatiotemporal population\nstructure.” This new title more accurately and professionally reflects the content of\nthe section.\n-->Multiple instances of missing spaces after punctuation, incomplete sentences, and\nwords missing letters.\nResponse:\nIn the revised version, we have thoroughly reviewed and corrected instances of\nmissing spaces, incomplete sentences, and truncated words throughout the\nmanuscript.\n-->Some terminology is used inconsistently (e.g., \"NTSR-associated SNPs\" vs\n\"resistance-associated SNPs\" vs \"rSNPs\")\nResponse:\nSorry for the inconsistency in terminology. In the revised manuscript, we have\nstandardized the usage throughout the text by consistently referring to these variants\nas “rSNPs” (resistance-associated SNPs). This change improves clarity and\nconsistency across the manuscript.\n-->Please look to see where run-on sentences can be split for clarity.\nResponse:\nThank you for this valuable suggestion. In the revised manuscript, we have carefully\nreviewed the text and revised several run-on or overly long sentences to improve\nclarity and readability.\n- Figure legends: Several figure legends could be more informative and self-\ncontained:\n-->All abbreviations should be defined in each legend\nResponse:\nSorry for the mistakes. We have added abbreviations in the legends.\n-->Statistical tests used should be mentioned in the legend when significance is\nindicated\nResponse:\nWe have added all used statistical tests in the legends accordingly.\n-->The legends should allow readers to understand the main takeaway without having\nto refer back to the main text\nResponse:\nWe have revised the figure legends to clearly convey the main findings, so that they\ncan be understood independently without needing to refer back to the main text.\n- Data availability: While the authors mention data will be made available, specific\ndetails are still missing:\n-->Accession numbers for genome assemblies\n-->Repository information for resequencing data\n-->Code availability for custom analyses should be included in a Data Availability\nstatement.\nResponse:\nThank you for this helpful comment regarding data availability. In the revised\nmanuscript, we have ensured that all data and code availability details are clearly\nstated in the Data and code availability section.\nSpecifically:\n– Accession numbers for all newly generated Digitaria genome assemblies are\nprovided and deposited at the National Genomics Data Center (NGDC) and described\nin the manuscript as: “Genome assemblies and annotations for the three Digitaria\nspecies generated in this study have been deposited in National Genomics Data\nCenter (NGDC) database (https://bigd.big.ac.cn) under the following accession\nnumbers: GWHGGEX00000000 (D. sanguinalis, #YJ2023), GWHGGEY00000000\n(D. radicosa, #YZGJ2) and GWHGGEZ00000000 (D. milanjiana, DZ2).”\n– Resequencing data is deposited at NCBI under accession number PRJNA1296823\nand it reads as: “Raw resequencing data for all accessions have been deposited in the\nNational Center for Biotechnology Information (NCBI) database\n(https:/www.ncbi.nlm.nih.gov) under the project accession number PRJNA1296823.”\n– Code availability has now been explicitly included in the Data and code\navailability section.\nWe appreciate the reviewer’s attention to this important aspect of reproducibility, and\nwe have revised the section heading and wording to ensure clarity and compliance\nwith journal standards."}, {"round": 3, "reviewer_comments": {"reviewer_1": "Overall, the manuscript's text has changed substantially in this last revision, and these\nchanges have clearly improved it. I have not checked all technical details (for\nexample, cross-referencing of figure numbering between the main text and figures /\nsupplementary materials), but the authors have addressed the points in the text that\npreviously raised scientific concerns.\nI appreciate the authors’ efforts to use succinct and precise language; however, in\nsome places a simpler sentence using more common language may better serve the\nmanuscript’s aims. If further changes are made, focusing on the abstract and\ndiscussion may have the greatest impact.\n- In the abstract, the sentence “Phenotypic classification of Digitaria was established\nbased on genome-scale data and extensive morphological profiling” might be\nclarified, for example: “This genomic data was combined with detailed observations\nof physical characteristics to generate a system for classifying Digitaria species.”\n- In the Discussion, the sentence “Unlike their large genome sizes (17 Gb for wheat\nand 11 Gb for oats), the relatively compact genomes of D. sanguinalis and its\nancestral lineages provide a system for studying the evolutionary consequences of\npolyploidization, particularly the structural dynamics of repetitive elements” could be\nsimplified, for example: “While wheat and oats have very large genomes (17 Gb and\nuseful model for studying how polyploidization influences genome evolution,\nparticularly the structural dynamics of repetitive elements.”\nI agree with the other reviewer that additional or improved transition and explanatory\nsentences in the Introduction would help make this substantial body of work more\napproachable. This may partly reflect my earlier suggestion to rearrange this section,\nand I apologise if so. If there is another revision, I suggest the authors to consider this\nagain. The added transitions in the Results section improve readability, a few similar\nadditions in the Introduction may be beneficial.\nThere remain a few places where the language could be tightened (for example,\n“Polyploidization is widespread across the plant kingdom and [is] particularly\nprevalent among weedy species, where it often contributes to increased resistance to\nstress…”), but these are largely cosmetic and could be addressed during final editing.\nThe changes shown in the tracked version have clearly improved the manuscript, and\nI recognise that any manuscript can always be further refined.\nIn summary, this represents a considerable and substantial amount of work, and the\nauthors have done a very good job presenting it in a compact and readable form so\nthat users can use it.", "reviewer_2": "Thank you for your revisions."}, "author_response": "A Point-to-point response to the reviewers’ comments\nReviewers’ comments are shown as black text. Response to comments is shown as\nblue text.\nReviewer #1 (Remarks to the Author):\nOverall, the manuscript's text has changed substantially in this last revision, and these\nchanges have clearly improved it. I have not checked all technical details (for\nexample, cross-referencing of figure numbering between the main text and figures /\nsupplementary materials), but the authors have addressed the points in the text that\npreviously raised scientific concerns.\nI appreciate the authors’ efforts to use succinct and precise language; however, in\nsome places a simpler sentence using more common language may better serve the\nmanuscript’s aims. If further changes are made, focusing on the abstract and\ndiscussion may have the greatest impact.\nResponse:\nWe are grateful for your constructive suggestions, which help us\nsubstantially improve our manuscript.\n- In the abstract, the sentence “Phenotypic classification of Digitaria was established\nbased on genome-scale data and extensive morphological profiling” might be\nclarified, for example: “This genomic data was combined with detailed observations\nof physical characteristics to generate a system for classifying Digitaria species.”\nResponse:\nThank you for the helpful suggestion. We have revised the sentence in the Abstract to\nimprove clarity as follows: “A phenotypic classification system for Digitaria species\nis developed by integrating genome-scale data with detailed morphological\nobservations.”\n- In the Discussion, the sentence “Unlike their large genome sizes (17 Gb for wheat\nand 11 Gb for oats), the relatively compact genomes of D. sanguinalis and its\nancestral lineages provide a system for studying the evolutionary consequences of\npolyploidization, particularly the structural dynamics of repetitive elements” could be\nsimplified, for example: “While wheat and oats have very large genomes (17 Gb and\nuseful model for studying how polyploidization influences genome evolution,\nparticularly the structural dynamics of repetitive elements.”\nResponse:\nThank you for your helpful suggestion regarding this sentence in the Discussion.\nFollowing your recommendation, we have revised the sentence to improve clarity,\nwhile maintaining a formal academic tone, as follows:\n“In contrast to the large genomes of wheat (17 Gb) and oats (11 Gb), the compact\ngenomes of D. sanguinalis and its ancestors serve as an effective model system for\nstudying how polyploidization influences genome evolution, particularly the\nstructural dynamics of repetitive elements.”\nI agree with the other reviewer that additional or improved transition and explanatory\nsentences in the Introduction would help make this substantial body of work more\napproachable. This may partly reflect my earlier suggestion to rearrange this section,\nand I apologise if so. If there is another revision, I suggest the authors to consider this\nagain. The added transitions in the Results section improve readability, a few similar\nadditions in the Introduction may be beneficial.\nResponse:\nThank you very much for this valuable suggestion. We have taken prior comments\ninto account to enhance the transitions in the Introduction, and we believe the current\nversion reads more smoothly.\nThere remain a few places where the language could be tightened (for example,\n“Polyploidization is widespread across the plant kingdom and [is] particularly\nprevalent among weedy species, where it often contributes to increased resistance to\nstress…”), but these are largely cosmetic and could be addressed during final editing.\nThe changes shown in the tracked version have clearly improved the manuscript, and\nI recognise that any manuscript can always be further refined.\nResponse:\nWe sincerely thank you for your thoughtful comments regarding language refinement.\nWe have revised the sentence according to the example you provided. In addition, we\nhave made broader edits throughout the manuscript to further enhance clarity,\nprecision, and overall professionalism.\nIn summary, this represents a considerable and substantial amount of work, and the\nauthors have done a very good job presenting it in a compact and readable form so\nthat users can use it.\nResponse:\nThank you very much for the thoughtful and encouraging feedback. We truly\nappreciate your recognition of the clarity and usability of our work.\nReviewer #2 (Remarks to the Author):\nThank you for your revisions."}]

Decision letter