Sequence characterization of eccDNA content in multiple blackgrass (Alopecurus myosuroides) populations

1. Introduction

Extra-chromosomal circular DNAs (eccDNA) are ring-shaped amplifications of linear genomic segments that exist separate from the autosomal genome. These structures have been observed in various species, including plant and higher eukaryote species [1,2,3,4,5,6]. Bassel and Hoota first discovered eccDNA in 1964 [7], but limited progress in understanding their function or genesis was made until the past decade. The formation of eccDNA can occur through multiple mechanisms, including non-allelic homologous recombination, replication slippage and DNA repair [8,9], and intramolecular recombination and retrotransposon activity [1,10,11,12]. Historically, the detection of eccDNA has been challenging due to their complex and highly repetitive structures. However, recent advancements in hybrid sequencing technologies (Illumina, PacBio, Oxford Nanopore) and computational algorithms (AmpliconArchitect [13], CiderSeq [14], and Circle-Seq [15]) have improved the ability to detect eccDNA and define tissue-specific "DNA circulome" [13,16].

EccDNAs have recently garnered significant attention in scientific communities as a unique hallmark of somatic mutations. These entities are capable of amplifying gene copies, related transcripts, and regulatory elements. eccDNAs can range in size from a few hundred base pairs [1] to much larger structures near a megabase in size [11,17], with a majority being smaller than 10 kb [1,18]. EccDNAs accumulate in cells through either known mechanisms of autonomous replication mediated by a replication origin [11,12,19,20], rolling circle amplification [21], DNA replication and mitosis [22,23] or some elusive mechanism. Notably, the accumulation of many replicate copies of a large eccDNA in Palmer amaranth has been shown to significantly increase the c-value of cells, seemingly without fitness penalty [24].

EccDNAs have been identified in various human physiological and pathological states, including normal cells (blood and muscle) [25], the central nervous system and neurodegeneration [26], aging and telomere maintenance [27,28], and cancer malignancy and progression [16]. eccDNAs have been observed in Drosophilia [4], and are believed to play a role in stress adaptation in yeast [29,30] and herbicide resistance in weeds [19]. The eccDNAs present in the noxious weed Palmer amaranth have been well characterized. This eccDNA, which is the largest gene-containing eccDNA reported thus far in a plant species, confers genetically transmissible resistance to the herbicide glyphosate [11,19]. This eccDNA replicon has multiple copies of the gene encoding for 5-enolpyruvylshikimate-3-phosphate (EPSP) synthase (EPSPS) and serves as the basis for gene copy amplification and therefore resistance to the EPSPS-inhibiting herbicide glyphosate. A study comparing eccDNA from glyphosate-sensitive and resistant populations revealed a significant content of eccDNA in glyphosate-sensitive populations [11]. The EPSPS gene was only present on eccDNA in glyphosate-resistant populations, while glyphosate-sensitive eccDNA displayed the genomic disposition for recombination events that may lead to the formation of the eccDNA replicon [1]. High-resolution cytological analysis of the eccDNA replicon revealed that its persistence in the germline is facilitated through a chromosomal tethering mechanism, rather than genomic re-integration [19]. This observation provides an explanation for the uneven segregation reported in previous studies [19]. Sequence analysis of the eccDNA replicon in Palmer amaranth revealed that it contains 58 genes in addition to the EPSPS gene, as well as a complex distribution of repeat sequences and gene functions that encompass detoxification, transport, replication, and stress response [11]. A synteny analysis of the other genes aligned to a neighboring species revealed that they originated from multiple chromosomes, suggesting a complex biogenesis event [11]. This study suggests that eccDNA serve as a rapid response mechanism and facilitate rapid adaptation [1].

Blackgrass (Alopecurues myosuroides) is a pernicious weed that has rapidly evolved to infest cereal crops across Northern and Western Europe and is now widely regarded as the top weed problem in several regions [31,32]. In response to the extensive use of herbicides, black-grass has evolved resistance to seven herbicide modes-of-action, exacerbating this species’ impact [33]. Blackgrass has been estimated to cause an annual £0.4 billion lost gross profit in the UK alone, due to wheat yield losses of 0.8 million tonnes [34]. Growers are now increasingly relying upon mixtures and sequences of (particularly pre-emergent) herbicides from a range of modes-of-action to manage this species, alongside physical or cultural practices. Nevertheless, further herbicidal control is complicated by the extensive presence of both gene mutations and multiple mechanisms that inactivate herbicides in this species [35], evolving independently within separate populations [36].

Advancements in comprehending the genetic basis of herbicide resistance have been limited. Nevertheless, recent studies have uncovered the genome of blackgrass, revealing independent and semi-overlapping genetic mechanisms underlying Non-Target Site Resistance (NTSR) in blackgrass at both DNA and RNA levels [36]. We used the same genetic material in this study to perform sequence characterization of extra-chromosomal circular DNA (eccDNA) content in multiple blackgrass populations with well-characterized levels of heritable NTSR-based herbicide resistance. We compare the genomic similarities and differences among individuals from these distinct populations. Additionally, we determine the gene content of the eccDNA and identify genomic hotspots for eccDNA biogenesis, evaluating their role in the rapid evolution of herbicide resistance.

2. Materials and Methods

2.1. Blackgrass populations, plant material, and genomic DNA extraction

Seven Blackgrass seed populations were grown for extraction of eccDNA, representing both herbicide sensitive (HS) and herbicide resistant (HR) lines. The Roth (HS) population originates from the Rothamsted ‘Broadbalk’ long-term field experiment and represents a wild-type population which has never experienced herbicide. The ‘Peldon’ and ‘Lola91’ are field collected HR populations confirmed to have non-target-site resistance to a range of herbicide chemistries. For these seed populations, we have selected against any known target-site mutations [37]. Further HR and HS samples were derived from segregating families created through biparental crossing of Roth and Peldon (CC2) and Roth and Lola 91 (CC5). The CC2R and CC2S are F₃ generation lines, caused by identifying and bulk-crossing the most and least resistant CC2 individuals within the F₂ generation. The CC5R and CC5S were produced in an identical way, but from the CC5 seed family. In total this provides three herbicide sensitive (Roth, CC2S, CC5S) populations, and four herbicide resistant (Peldon, Lola91, CC2R, CC5R) populations with which to examine eccDNA.

Seeds from these seven populations were broadcast sown onto a standard seed tray (37.3 x 24.5 x 5.5 cm) filled with Rothamsted Prescription Mix (75% Medium grade (L&P) peat, 12% Screened sterilized loam, 3% Medium grade vermiculite, 10% Grit (5mm screened, lime free, with 3.5kg Osmocote Exact 3/4 month per m3 - Supplier: Scotts UK Professional, Ipswich, Suffolk.). Plants were grown in standard glasshouse (GH43 101 t) at 10°C/5°C 16-hour day length without supplementary lighting for 40 days (between 15 December 2021 to 25 January 2022) until they were at the 2-3 leaf stage. From each population, two bulk harvests of aerial tissue from these plants totaling 0.5 g were cut, flash frozen in liquid nitrogen, and ground completely to a fine powder in a mortar and pestle. 100 mg of this ground material was then used for DNA extraction.

DNA was extracted from 100 mg of finely ground flash frozen bulked leaf material from 40-day old seedlings using a modified CTAB method. More specifically, 20μl of beta-mercaptoethanol was added to the samples along with 800μl of 2X CTAB buffer (2% Cetyltrimethylammonium Bromide (Merck 219374 – Calbiochem), 1% PVPP, 1.4M NaCl, 20mM EDTA pH 8.0, 100mM Tris pH 8.0) that had been heated to 80-90˚C. Samples were vortexed to mix. 800μl of chloroform was added, vortexed until completely mixed, and tubes were incubated at 65˚C for 10-15min allowing the phases to separate. Samples were centrifuged at max speed for 30min. The supernatant was transferred to new, labeled 2ml tube. A second aliquot of 800μl chloroform was added, tubes were inverted to mix, and centrifuged at maximum speed for 20min. 700μl of the supernatant was transferred to new labeled 1.5ml microfuge tubes. 450μl of cold (4˚C) isopropanol was added and tubes inverted gently. Samples were left to incubate at room temperature for exactly 5min, then centrifuged at maximum speed for 15min. Supernatant was removed, and the pellet was washed with 1ml of 70% ethanol at -20˚C. 70% Ethanol was removed and pellets dried using a speedvac for 7 minutes with medium heat. Pellet was resuspended in 100μl TE buffer with a 1:1000 dilution of RNase A Solution at 10mg/ml (Merck 70856-3) and incubated at 37˚C for 45min to resuspend the DNA into solution. Samples were quantified by nanodrop analysis then frozen and stored at -20˚C until being shipped on dry ice for eccDNA enrichment and sequencing.

2.2. EccDNA sequence enrichment, sequencing, and analysis

Circular DNA enrichment and sequencing was performed according to the CIDER-seq methods [14]. Circular DNA sequencing was performed on a Pacific Biosciences Sequel II instrument (Pacific Biosciences) with a 30-hour movie time and HiFi/CCS read generation. Sequence reads were deconcatenated into circular DNA using the CIDER-seq software [14]. The circular DNA sequences were then aligned to the reference genome of Alopecurus myosuroides [14] to predict the extra-chromosomal circular DNA (eccDNA) content. To reduce redundant sequences, the eccDNA sequences were clustered and collapsed into non-redundant representatives using the CD-hit program [38] with a sequence identity threshold of 0.9. The representative eccDNA sequences were masked for repetitive elements using the Repeat-Masker [39] with default settings and annotated using the MAKER annotation pipeline [40].

GO enrichment analysis was conducted on the predicted coding elements of the eccDNA using the GOseq 1.42.0 R package [41]. GO terms were considered significantly enriched at a False Discovery Rate (FDR) of less than 0.05. The enriched GO terms were visualized using the ggplot2 package [42]. tRNAs were exclusively predicted for the eccDNA sequences without coding sequences (CDS) using the tRNAscan-SE program [43] with default settings. To determine the genomic origin of the eccDNA, the predicted sequences were first aligned to the reference genome using the Minimap2 program [44]. The number of eccDNAs that overlapped with the reference genome in each 500 kb interval was counted using the BEDtools software [45]. The number of genes and transposable elements (TEs) were also counted within 200 kb downstream and upstream of each gene using BEDtools [45]. A permutation test was performed to compare the mean number of genes and TEs between the eccDNA-present and absent regions.

3. Results

3.1. EccDNA content and coding structure in multiple blackgrass populations

The CIDER-seq pipeline was used to identify a high abundance of extra-chromosomal circular DNA (eccDNA) in all blackgrass samples (Table 1). The number of non-redundant representative eccDNAs detected in the blackgrass samples ranged from 4233 in sample CC2S to 5663 in sample Roth, with an average of 4886 (Table 1). The size of the detected eccDNAs showed high variation within the samples, as evidenced by the length distributions that ranged from 31 bp to 29,081 bp and had mean lengths of around 6,900 bp (Table 1 and Figure 1). The length distributions of the eccDNAs were similar among all sequenced samples (Table 1 and Figure 1).

We further analyzed the coding content of the predicted eccDNA sequences. The number of eccDNAs with predicted genes ranged from 958 in sample CC5S to 1291 in sample Roth (average of 1079 per sample) and accounted for up to approximately 22% of the total predicted eccDNAs (Table 1). The number of predicted genes per eccDNA ranged from 1 to 15, with an average of 2 genes per eccDNA across all analyzed samples (Table S1 and S2). Transfer RNA (tRNA) prediction was performed for the eccDNAs and revealed a relatively low number of predicted tRNAs in all samples, ranging from 45 in sample CC5S to 77 in sample CC2S.

3.2. Coding content of eccDNAs in herbicide resistant and sensitive blackgrass populations

To identify common and unique eccDNA-encoded genes with predicted functional protein domains among different blackgrass populations, we compared the predicted gene content from both herbicide resistant (HR) (CC2R, CC5R, Lola91 and Peldon) and sensitive (HS) (CC2S, CC5S and Roth) blackgrass samples. A comparison of the predicted protein domains (PFAM id) identified a total of 80 protein coding domains that were shared by all HR samples (Figure 2 and Table S3). Among these functional protein domains, the number of eccDNAs annotated to them ranged from 6 to 416. There was a range of 7 to 25 functional annotations that were shared by at least two HR samples. The different HR samples had varying numbers of unique functional protein domains, ranging from 68 (Peldon) to 111 (CC5R) (Table S3). Overall, the three HS samples shared a total of 91 functional protein domains (Figure 3 and Table S4). Pairwise comparisons between the HS samples revealed 18 to 38 protein domains were shared by each of the two samples. In addition, 96 (CC2S) to 115 (Roth) unique functional domains were identified in different HS samples. By comparing the 91 functional protein domains shared by the three HS samples to the 80 protein coding domains shared by all four HR samples, a list of 69 functional protein domains that are common to all measured samples was obtained. Additionally, there are 11 functional protein domains specific to all the HR samples and 22 functional protein domains specific to all the HS samples (Table S5). Counts of these core set of functional protein domains found in all samples range from 416 to 6 in the HR samples and 5 in the HS samples and include putative gypsy type transposons, Proton-conducting membrane transporters, ribosomal protein domains, and domains associated with the respiratory chain, such as NADH-ubiquinone/plastoquinone oxidoreductase or NADH dehydrogenase, and the photosynthetic machinery, such as Photosystem I psaA/psaB protein and Photosystem II protein (Table S5). Counts for the functional protein domains that were unique to either HR (14 to 6) or HS (15 to 3) were comparatively low and are discussed below.

In the HR samples, the functional domains with the highest abundance comprised proteins associated with photosystems, gypsy type transposons, NADH dehydrogenase, PPR repeat, and Cytochrome C assembly protein (Table 2). Moreover, we observed the presence of eccDNA genes in the analyzed HR samples, which were previously reported to be associated with herbicide detoxification. Cytochrome P450 was detected in all HR samples, whereas ATP-binding cassette transporters (ABC transporters) and glutathione transferases were exclusively detected in the CC2R, CC5R, and Peldon populations (Table S3). Furthermore, within the HR samples, various stress response domains were identified. Specifically, ribosomal proteins and leucine-rich repeats were present in all HR samples. The WRKY DNA-binding domain was shared among the CC2R, CC5R, and Lola91 samples. Moreover, the Myb-like DNA-binding domain and WD domain were shared by the CC5R, Lola91, and Peldon samples. Finally, peroxidase was detected in the CC2R, Lola91, and Peldon samples.

The functional domains detected in herbicide sensitive (HS) samples were found to be abundant and largely overlapping with those detected in herbicide-resistant (HR) samples, as demonstrated in Tables 2, 3 and S3. Three functional domains, namely the Jacalin-like lectin domain, the GDSL/SGNH-like Acyl-Esterase family, and the PMR5 N-terminal domain, were present exclusively in the three HS samples. Additionally, annotations related to herbicide detoxification were identified in the HS samples, including the presence of cytochrome P450 and glutathione transferases in all HS samples and an ABC transporter annotated only to predicted eccDNA genes in CC5S and Roth (Table S4).

To better understand the inheritance of eccDNA in blackgrass, we further compared the functional domains separately in the parents (Figure 4), CC2 (Figure 5) and CC5 (Figure 6) populations. In total, 99 functional domains were shared by the three parents, and 26 functional domains were shared by the two HR parents (Peldon and Lola91), while 123 functional domains were only detected in the HS parent (Roth) (Figure 4 and Table S6). Some defense related functional domains such as DnaJ domain, glycosyl hydrolase and response regulator receiver domains were uniquely identified in the HR parents.

Upon comparing the two distinct biotypes within the CC2 population, we found that functional annotations associated with ABC transporters were exclusively identified in the CC2R (HR) samples (Table S7). In contrast, within the CC5 population, functional annotations linked to ABC transporters were detected in both the CC5R and CC5S samples.

To gain insights into the functions of the predicted extra-chromosomal circular DNA (eccDNA) genes, we conducted a homology-based analysis with previously identified blackgrass herbicide-resistant candidate genes [46]. Our results showed that eccDNAs homologous to the crucial herbicide-resistant gene, GSTF1 (ALOMY3G11302 [36]), were detected exclusively in the herbicide-resistant (HR) samples CC2R and CC5R (Table 4). Of particular interest, this sequence is 39 amino acids shorter than the other four gene sequences previously identified [47] but otherwise most closely aligns with the GST2c sequence. ALOMY3G11302 was differentially expressed in both CC2R and CC5R populations, as reported by [36]. EccDNAs homologous to GSTU2 were present in both HR and HS samples. The presence of eccDNAs homologous to OPR1 and GSTF2 was restricted to Peldon and CC5S, respectively. We also found sequences that associated with PF00662 (NADH-Ubiquinone oxidoreductase) common all the samples (Supplementary Table S5) as well as those associated with PF01370 the NAD-dependent epimerase/dehydratase in the CC5R, Peldon and Roth samples (Supplementary Table S1).

3.3. Gene ontology enrichment of blackgrass eccDNA

Gene ontology (GO) enrichment analysis of predicted coding contents in eccDNA samples showed a list of GO terms associated with biological processes, cellular components, and molecular function were enriched in HR and HS samples. In HR samples, the enriched biological processes include translation, transcription, protein neddylation, photosynthesis, cytochrome complex assembly, ATP-related processes, and galactose metabolic process (Figure 7a and Supplementary Table S9). For the cellular component category, GO terms that related to photosystem, thylakoid, and ribosome were enriched in HR samples (Figure 7b and Supplementary Table S9). Representative molecular functions for HR eccDNAs include zinc ion binding, UDP-glucose, structure constituent of ribosome and chlorophyll binding (Figure 7c and Supplementary Table S9).

Compared to the HR samples, similar but slightly different GO terms were enriched in HS samples. Biological processes such as chlorophyll catabolic process, glucose catabolic process, RNA metabolic process, heme oxidation, phosphate ion transport and proteolysis were only enriched in HS samples (Figures 7a, 8a and 9 and Supplementary Table S9). While biological process such as ATP synthesis, protein neddylation, transcription, tryptophan metabolic process and cellular manganese ion homeostasis were uniquely enriched in HR samples. All GO terms related to cellular components that were enriched in HR samples were also enriched in HS samples except for DNA replication factor, exosome and H4/H2A histone acetyltransferase complex (Figures 7b, 8b and 9 and Supplementary Table S9). Molecular functions such as chlorophyllase activity, enzyme inhibitor activity, phosphoglycerate mutase activity and double-stranded DNA binding were only identified in HS samples (Figures 7c, 8c and 9 and Supplementary Table S9). While GO terms related oxidoreductase activity, ADP binding, exonuclease activity and manganese ion transmembrane transporter activity were uniquely enriched in HR samples.

Figure 9. Circos plot showing the count of eccDNA identified in herbicide resistant (blue) and sensitive (orange) blackgrass population. Outer ring makes the number of genes (purple) and transposable elements (green) in the reference genome.

Figure 9. Circos plot showing the count of eccDNA identified in herbicide resistant (blue) and sensitive (orange) blackgrass population. Outer ring makes the number of genes (purple) and transposable elements (green) in the reference genome.

3.4. Genomic origins of blackgrass eccDNA

To determine the genomic origins of eccDNA, the predicted eccDNAs were mapped to the Alopecurus myosuroides reference genome [48]. Numbers of mapped eccDNAs were counted within the non-overlapped windows of 500kb through the whole genome (Supplementary Table S11 and Figure 6). Over 90% of the detected were successfully mapped to the established chromosomes, and around 8.5% of the predicted eccDNAs were mapped to the unanchored sequences. Several genomic segments were identified with a high frequency of eccDNAs being mapped to these regions. These regions include the distal end of chromosome 5 and the near center of chromosome 1 and 7 (Figure 6 and Table 5). The 500kb window localized at the near center of chromosome 1 contained 83 eccDNAs from HR samples and 48 eccDNA eccDNAs from HS samples. The two 500kb window located at the distal end of chromosome 5 contained 237 (HR 133, HS 104) and 137 (HR 82, HS 55) predicted eccDNAs, respectively. And the 500kb window on chromosome 7 mapped 130 and 80 eccDNAs from HS and HR samples, respectively. Stress-related genes, such as ribosomal family, lysine decarboxylase and HPPK, were located in these regions. However, no eccDNA coding contents were mapped to these regions.

In order to gain further insight into the genomic features of the regions enriched with eccDNA-associated genes, we analyzed whether the eccDNA-associated genes originated from regions with a higher gene density than expected by chance (Figure 6). We observed that eccDNA present regions tend to have higher number of genes in comparison to eccDNA absent region in both HR and HS samples. In both HR and HS sample, we observed significant differences (permutation test for difference of means, n = 10,000 replicates) between the number of genes located at eccDNA present region and eccDNA absent region. The mean number of observed genes in eccDNA present region was increased by 0.24 (p = 0.0019) and 0.26 (p = 0.0044) in comparison to eccDNA absent region in HR and HS samples, respectively. Comparing to eccDNA absent region, eccDNA present region tend to have less TE (HR = 3.16, HS = 1.37), however, the difference of number of TEs in eccDNA present region and absent region was only significant in HR samples (p = 0.026). Notably, no significant differences were observed between the number of observed genes or TE in eccDNA present region and eccDNA absent region when all samples were analyzed together.

A prior study identified a list of quantitative trait loci (QTL) regions that are associated with herbicide resistance in various blackgrass populations [48]. To investigate the mapping of eccDNAs to these regions, we found that 496 HR eccDNAs and 352 HS eccDNAs were mapped to the blackgrass herbicide resistant QTL regions (Supplementary Table S11). Our analysis revealed that within the QTL regions, two non-overlapping windows of 500kb located on chromosome 3 (230.0–230.5 Mb and 312.5–313.0 Mb) were uniquely mapped to at least three of the HR samples (Supplementary Table S11). EccDNA coding sequences were mapped to 12 QTL regions (Supplementary Table S13), where the qtl-cc2-2-1, qtl-cc2-2-3, qtl-cc2-5-1, and qtl-cc2-5-2 regions only contained HS eccDNA coding sequences (CDS), while the qtl-cc5-3-3 region only contained HR eccDNA CDS. Additionally, we observed that among the eccDNA CDS mapped to the QTL regions, only HR eccDNA genes were successfully annotated with functional predictions. Most of these genes were predicted to be related to photosystem and gypsy-type transposon (Supplementary Table S13).

4. Discussion

Extra-chromosomal circular DNA (eccDNA) are DNA molecules that are separate from the main chromatin body, or chromosomes, within a cell and are gaining significant attention with their role in trait biology [1,5,9,11,16,20,27]. This form of DNA has been demonstrated to play a crucial role in various physiological processes. One of the most prominent functions of eccDNA is its capacity to harbor and amplify segments of transcriptionally active chromatin, resulting in gene focal amplifications. Furthermore, eccDNA acts as a reservoir of genetic diversity that the cell can draw upon during exposure to stress, whether it be biotic or abiotic, facilitating rapid responses and evolutionary adaptations [1,49]. The extensive use of herbicides over the years has led to the evolution of blackgrass populations that are resistant to multiple herbicides, including those from different chemical classes [35]. Aside from the rapid evolution of herbicide resistance, adaptive traits such as a lifecycle that occurs within a standard winter-cropping time, high fecundity, rapid growth, and ability to compete with the crop for nutrients, light and water, have enabled blackgrass to thrive in agricultural landscapes across Northern and Western Europe and become a major challenge for these farmers.

In this research, we have uncovered a wide array of eccDNAs within several black-grass populations, showcasing varying levels of herbicide resistance. Our findings in black-grass, concerning size, prevalence, and coding content, align closely with Amaranthus palmeri [1]. In general, we identified a suite of encoded genes and transposable elements associated with signaling, transport, and DNA mobility in both herbicide resistant and sensitive populations indicating a possible role in adaptability. For example, a number of gypsy retroelements were found on eccDNA which have been previously implicated with a role in creating genetic diversity [50]. We identified a significant number of proton conducting membrane transporter genes (Table S3) which function in the regulation of proton movement across cell membranes and have been associated with enabling plants to modulate their responses to various environmental stressors such as salinity, drought, and metal toxicity, which enhance their adaptability [51,52]. We also found an abundance of NADH dehydrogenase functional domains generally in the eccDNA coding dataset. This gene class is part of the mitochondrial electron transport chain and is involved in energy production, redox balance, and respiratory metabolism. Presence on an eccDNA may influence expression amplitudes of these genes and contribute with the plant’s ability to cope with stress through maintenance of cellular energy levels and redox homeostasis during adverse environmental conditions [53,54]. Ribosomal proteins were also highly abundant indicating a need for the ability for protein turnover, which has also been reported in [1]. Transcription factors such as WRKY, MYB, and other classes (Table S3) were also generally found in the blackgrass populations on eccDNA. These classes of transcription factors are pivotal regulators in plants that enable the activation of key processes such as stress responses, defense mechanisms, and various developmental process which are critical in allowing plants to integrate these signals and fine-tune the adaptative potential [55,56,57]. Presence of these transcription factors on eccDNA indicate their relevance as rapid response elements.

A principal goal of our study was to examine signatures of adaptation in the eccDNA content between blackgrass populations that have rapidly evolved metabolic herbicide resistance with those that have not. Interestingly, gene ontology enrichment analysis between the contrasting populations did not directly identify detoxification as being enriched exclusively in HR populations (Figure 7). Many of the same processes were enriched in both biotypes that include chlorophyll and glucose catabolism, RNA metabolism, heme oxidation, phosphate transport, and proteolysis. However, these processes are fundamental to plant growth and serve as an antioxidant in some capacities.

Investigation at the individual gene level did identify genes previously implicated with herbicide resistance [46]. A key gene found on an eccDNA in our study exclusive to the two resistant populations (CC2R and CC5R) was the GSTF1 (Glutathione S-Transferase Phi 1 - ALOMY3G11302 [36]). This gene belongs to the GST gene family which encodes enzymes involved in detoxifying herbicides by catalyzing the conjugation of glutathione to herbicide molecules making them less toxic and more water-soluble for excretion or transport out of plant cells [47]. This gene has been associated with resistance to multiple classes of herbicides, including triazine and atrazine, indicating a role in multiple herbicide resistance [47]. Alternate alleles for GSTF1 have been discovered where some alleles may confer higher herbicide resistance than other. Notably, the GSTF1 in this study is 39 amino acids shorter than previously reported isoforms [47]. The presence of this gene on an eccDNA with novel genetic variation in our study suggests a possible role in herbicide resistance and that eccDNA based gene amplification and/or ultra-expression may also be critical factors to consider. Our study also identified GSTU2, OPR1, and GSTF2 as present of eccDNA, but not exclusively to HR populations (Table 4). These genes have also been implicated in herbicide resistance; however these results may also suggest that perhaps detoxification may not fully support the explanation for the resistant phenotype observed in our HR populations.

The genomic mechanisms that contribute to the biogenesis of eccDNA are still not well understood. Some explanations include illegitimate replication (replication slippage), non-homologous end joining, homologous recombination, microhomology-mediated break induced replication [58,59,60], and likely some organized or directed mechanisms, such as the EPSPS gene amplification that confers glyphosate resistance [11]. It is speculated that eccDNA can exist in multimeric forms and undergo recombination as a mechanism of compounded building into larger complex structures [1,19]. Analysis of eccDNA origins revealed eccDNA biogenesis on every chromosome with hotspots (high frequency of eccDNA formation) on chromosomes 1 and 7 (Figure 6). The identification of eccDNA originating from all over the genome is consistent with Amaranthus palmeri [1], as well as recent studies in mammalian cells [61,62]. These data provide further evidence that eccDNA biogenesis contributes to evolutionary innovation by contributing global genetic heterogeneity, plasticity, and a critical element to a plant’s trajectory to rapid adaptation. To further expand on this thought, we considered the concept of a two-speed genome recently proposed by eccDNA analysis of the rice blast fungal pathogen [63]. This concept considers that segments or regions within a genome evolve at different rates, and that perhaps, this concept can correlate with genes found as genomic focal amplifications. The study considered genes found on eccDNA and their origin as within a gene rich/poor region of the genome, in addition to proximity to repetitive elements and found that genes on eccDNAs were in fact under faster evolutionary rates when compared to genes not found as eccDNA [63]. We examined this concept with eccDNA in this study and only found a statistically significant result when comparing the HR and HS populations together, suggesting alternate evolutionary pressures between HR and HS populations presented here. Perhaps the selective pressure and adaptative signatures that correspond to herbicide resistance have broader implications genome-wide and further support mechanisms that provide resistance beyond metabolic genes. It has been proposed that plants that have rapidly evolved to a xenobiotic pressure also have a new disposition to withstand other abiotic stresses [64]. The compartmentalized genome architectures of the pathogens could also be organized in a way that these evolutionary signatures are more pronounced when compared to large plant genomes like blackgrass (3.2Gb).

5. Conclusions

Our analysis provides evidence of evolutionary innovation and useful insights into the abundance, coding content, functional domains, and biogenesis patterns of eccDNA in various blackgrass populations, with a focus on differences between herbicide-resistant and herbicide-sensitive populations. These findings demonstrate that the pernicious weed blackgrass (Alopecurus myosuroides) carries eccDNA with functional domains associated with genes previously reported to be associated with herbicide detoxification. As the eccDNA content is not identical between herbicide resistant plants descending from different families exhibiting equivalent NTSR phenotypes, these differences may help to explain the polygenic enhanced metabolic resistance seen in these plants. It is interesting to note that unlike EPSP synthase in Palmer amaranth, as blackgrass plants exhibiting sensitivity to herbicides also carry eccDNA with protein domains associated with herbicide detoxification, differences in eccDNA content cannot be fully responsible for the different phenotypes observed. Our data add blackgrass to the expanding list of organisms including human cells, yeast, Drosophilia, and Palmer amaranth and add to the mechanistic knowledge of how they confer herbicide resistance, alongside acting in stress adaptation and other physiological and pathological states.