THE TRANSCRIPTOMIC AND EPIGENOMIC RESPONSE OF KOCHIA SCOPARIA TO SUBLETHAL GLYPHOSATE By Carly Abbegail Claucherty A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Crop and Soil Sciences–Master of Science 2022 ABSTRACT THE TRANSCRIPTOMIC AND EPIGENOMIC RESPONSE OF KOCHIA SCOPARIA TO SUBLETHAL GLYPHOSATE By Carly Abbegail Claucherty Weed populations respond and adapt to herbicide stress by evolving resistance. Glyphosate resistance is primarily caused by the amplification of the target site gene, EPSPS, where multiple copies produce a large enough protein pool so that field rates do not kill the plant. This mechanism has evolved independently in at least nine divergent weed species. It has been demonstrated that EPSPS gene duplication may be transposon mediated in Kochia scoparia. A key regulator of transposable element (TE) activity is DNA methylation. The role of the epigenome and subsequent transcriptome in transient responses to herbicides of their primary target, weeds, is not well understood In this study, we performed RNA-Seq and bisulfite sequencing on leaf tissue from glyphosate-sensitive kochia before and three weeks after treatment with two sublethal doses to determine if glyphosate causes hypomethylation of the genome, allowing for the activation of transposons and upregulation of stress-related genes. Our results shows that overall gene expression was suppressed by glyphosate and increases in CHH methylation through development were also ceased. We did not observe significant global changes in cytosine methylation, and overall responses were stochastic. When combining the two datasets together, there was no direct correlation between changes in methylation and changes in gene expression suggesting that DNA methylation is not the primary cause of differential expression in our study. Our results broaden the knowledge pool of weedy species epigenomics and aid in understanding the contribution of DNA methylation to plant resilience in response to herbicide stress. ACKNOWLEDGEMENTS I would like to thank Dr. Eric Patterson for his continuous support and guidance over the past few years. Thank you for believing in my abilities and taking a chance on me as your first graduate student. Thank you to the rest of the Patterson lab, both past and present, for their support, especially Dr. Jinyi Chen for her patience with me in the lab, and Dr. Nathan Hall whom I got the privilege of sharing an office with and bugging with questions constantly. I am thankful for their mentorship and willingness to share their knowledge and skills with me over the course of my graduate school experience. Additionally, I’d like to thank Dr. Margaret Fleming for asking the right questions and pushing me to become a better scientific writer. I would also like to acknowledge the support of my graduate committee, Dr. Erin Burns, Dr. Emily Josephs, and Dr. Addie Thompson for their contributions to my personal and professional growth. Thank you to Dr. Chad Niederhuth as well, for your invaluable advice and help on epigenomics analyses. This work would not have been possible if not for generous funding received through the Academic Achievement Graduate Assistantship Award, and for the fellowship funding through the NRT-IMPACTS Research Traineeship Program. I am thankful for the IMPACTS coordinators, instructors, and the cohort for giving me a sense of community when it was an especially hard time to be a graduate student. Finally, a thank you to my family and friends that got me through these last two years. iii TABLE OF CONTENTS LIST OF TABLES ........................................................................................................................ vi LIST OF FIGURES ...................................................................................................................... vii KEY TO ABBREVIATIONS ........................................................................................................ ix CHAPTER I .....................................................................................................................................1 A REVIEW OF EVOLUTION AND HERBICIDE STRESS.........................................................1 Summary ...................................................................................................................................1 Introduction ..............................................................................................................................2 Evolutionary theory and herbicide resistance .....................................................................3 How herbicides cause stress and its perception through the plant .....................................5 Herbicide resistance mechanisms and their origins ...........................................................6 Assessing structural variation as an evolutionary tool for adaptation ...............................9 Unraveling the role of the epigenome ................................................................................11 Bringing it all together: The future of studying the origins of resistance ........................15 Conclusion ..............................................................................................................................16 REFERENCES ..............................................................................................................................18 CHAPTER II..................................................................................................................................27 KOCHIA SCOPARIA TRANSCRIPTOMIC RESPONSE TO GLYPHOSATE ..........................27 Summary ................................................................................................................................27 Introduction ............................................................................................................................28 Methods ...................................................................................................................................31 Plant material and sublethal dose determination ..............................................................31 Experimental design and tissue collection .........................................................................32 RNA Extraction and sequencing ........................................................................................33 Data processing and identification of differentially expressed transcripts .......................33 Gene Ontology analysis .....................................................................................................34 Analysis of specific gene families (shikimate pathway, stress response, transposons, epigenome) ...................................................................................................................34 Results .....................................................................................................................................35 Glyphosate injury symptomology.......................................................................................35 RNA sequencing mapping and data visualization .............................................................35 Analysis of differentially expressed transcripts .................................................................36 Characterizing the transcriptome response specific to high dose treatment .....................37 Gene Ontology analysis of genes with dose dependent expression ...................................37 Differential expression of stress related and known herbicide resistance gene families ..39 Transposon activation and suppression in response to glyphosate ...................................40 Genes regulating the epigenome .......................................................................................40 Discussion................................................................................................................................41 Glyphosate represses gene expression...............................................................................41 iv Treatment specific response to glyphosate ........................................................................42 Analysis of glyphosate induced stress-related genes .........................................................43 The role of transposons in glyphosate stress response ......................................................44 Conclusion ..............................................................................................................................44 APPENDIX ....................................................................................................................................46 REFERENCES ..............................................................................................................................72 CHAPTER III ................................................................................................................................80 EPIGENOMIC RESPONSE TO GLYPHOSATE IN KOCHIA SCOPARIA................................80 Summary .................................................................................................................................80 Introduction ............................................................................................................................81 The epigenome - DNA methylation maintenance and regulation of gene expression .......81 Epigenetics and plant stress...............................................................................................83 Transposable elements .......................................................................................................83 Herbicide resistance and epigenomics ..............................................................................84 Methods ...................................................................................................................................86 Plant material and sublethal dose determination for experiment .....................................86 Experimental design and leaf tissue collection ..................................................................87 DNA sample preparation and sequencing .........................................................................87 Data processing, mapping, and methylation extraction ...................................................88 Weighted methylation levels genome-wide and by different genomic features .................88 Pattern of CG, CHG and CHH methylation in genes ........................................................89 Analysis of differentially methylation regions ...................................................................89 High dose specific karyotypes ............................................................................................90 Correlation with previously generated gene expression data ...........................................90 Transposon families ...........................................................................................................91 Results .....................................................................................................................................91 Genome wide DNA methylation patterns in kochia ...........................................................91 DNA methylation patterns of genes ...................................................................................92 Analysis of differentially methylated regions and their associated genes .........................93 Changes in DNA methylation do not directly correlate to a change in expression ...........95 Methylation of transposons ................................................................................................96 Discussion................................................................................................................................97 Conclusion ............................................................................................................................101 APPENDIX ..................................................................................................................................103 REFERENCES ............................................................................................................................123 v LIST OF TABLES Table 2.1: RNA-Seq sample summary statistics. ............................................................................48 Table 2.2: High sublethal glyphosate responsive genes. ...............................................................53 Table 2.3: Stress related gene families. .........................................................................................56 Table 2.4: Upregulated genes shared between low and high dose glyphosate. ............................58 Table 2.5: Downregulated genes shared between low and high dose treatment. ..........................59 Table 2.6: Shikimate pathway. .......................................................................................................69 Table 3.1: Summary statistics from bisulfite sequencing data. ....................................................104 Table 3.2: High dose hypermethylated and hypomethylated CG DMRs. ....................................111 Table 3.3: Top CHG DMGs in high treatment. ...........................................................................114 Table 3.4: Top CHH DMGs in high treatment. ...........................................................................116 Table 3.5: DMRs with transposons. .............................................................................................122 vi LIST OF FIGURES Figure 2.1: Glyphosate dose response curve. ................................................................................47 Figure 2.2: Plant height and phenotypic response.........................................................................49 Figure 2.3: Gene expression MDS plot. .........................................................................................50 Figure 2.4: Volcano plot of DETs by treatment. ............................................................................51 Figure 2.5: Distribution of significant DETs between treatments. ................................................52 Figure 2.6: DET heatmap and cluster summary. ...........................................................................55 Figure 2.7: FHY3/FAR1 related genes heatmap. .......................................................................... 57 Figure 2.8: Biological Process terms from heatmap cluster one. ..................................................60 Figure 2.9: Cellular Component terms from heatmap cluster one. ...............................................61 Figure 2.10: Molecular Function terms from heatmap cluster one. ..............................................62 Figure 2.11: Biological Process terms from heatmap cluster four. ...............................................63 Figure 2.12: Cellular Component terms from heatmap cluster four. ............................................64 Figure 2.13: Molecular Function terms from heatmap cluster four. .............................................65 Figure 2.14: Biological Process terms from heatmap cluster five. ................................................66 Figure 2.15: Cellular Component terms from heatmap cluster five. .............................................67 Figure 2.16: Molecular Function terms from heatmap cluster five. ..............................................68 Figure 3.1: Genome wide methylation. ........................................................................................105 Figure 3.2: Methylation level of genes and their features. ..........................................................106 Figure 3.3: Methylation patterns of genes. ..................................................................................107 Figure 3.4: Volcano plot of DMRs in the CG context ..................................................................108 Figure 3.5: Venn diagrams of differentially methylated genes. ...................................................109 vii Figure 3.6: Volcano plot of DMRs in the CHG context. ..............................................................113 Figure 3.7: Volcano plot of DMRs in the CHH context. ..............................................................115 Figure 3.8: Karyotypes. ................................................................................................................117 Figure 3.9: Correlation of methylation to gene expression. ........................................................118 Figure 3.10: Weighted methylation of transposons......................................................................120 viii KEY TO ABBREVIATIONS ABC ATP-binding cassette ACCase Acetyl-CoA Carboxylase ALS Acetolactate Synthase CMT2 Chromomethylase 2 CMT3 Chromomethylase 3 BP Biological Process CC Cellular Component CNV Copy Number Variation CYP450 Cytochrome P450 DEG Differentially Expressed Gene DET Differentially Expressed Transcript DMG Differentially Methylated Gene DMR Differentially Methylated Region DMS Differentially Methylated Site eccDNA Extrachromosomal Circular DNA epiRIL epigenetic recombinant line EPSPS 5-enolpyruvylshikimate-3-phosphate synthase GO Gene Ontology GR Glutathione Reductase GST Glutathione-S-Transferases Log2FC Log2 Fold Change ix LRR Leucine-Rich-Repeats MET1 DNA Methyltransferase 1 MF Molecular Function NTSR Non-Target-Site Resistance phyA Phytochrome A POX Peroxidase PPO Protoporphyrinogen Oxidase RdDM RNA-Directed DNA Methylation ROS Reactive Oxygen Species ROS1 Repressor of Silencing 1 siRNA Small-Interfering RNA SNP Single Nucleotide Polymorphism SOD Super Oxide Dismutase SV Structural Variation TE Transposable Element TF Transcription Factor TGS Transcriptional Gene Silencing TMV Tobacco Mosaic Virus TSR Target-Site Resistance TSS Transcriptional Start Site TTS Transcription Termination Site UV Ultraviolet WAT Weeks After Treatment x CHAPTER I A REVIEW OF EVOLUTION AND HERBICIDE STRESS Summary Pesticides are vital tools for modern agricultural producers that safeguard crop yield from losses due to insects, pathogens, and weeds. Currently, herbicides account for half of all pesticide uses worldwide with glyphosate being the most common active ingredient. Over 8.6 billion kilograms of glyphosate have been applied across the globe since 1974, a historic magnitude that grew exponentially following the introduction of Roundup Ready crops in the mid-late 1990’s. However, severe reliance on this single synthetic herbicide has resulted in 55 different weed species evolving glyphosate resistance. In terms of evolution, this is rapid and yet, little is known about the origins of novel herbicide resistance mechanisms. Traditionally, novel variants in weed populations are thought to arise randomly, and through selection, increase in frequency. Evidence suggests that beyond natural mutations, induced mutations by various biotic and abiotic stressors are also an important source of novel resistance mechanisms. The influence of the environment on the genome incorporates Lamarckian concepts. We propose here that herbicide stress induces directed variation in the genome, by way of stress detection and epigenome modification, and this variation can then be selected upon. The relative contribution of “natural” or “induced” variation for herbicide resistance should be investigated further as it has implications for plant resilience and evolution. The sheer volume of herbicide application has profound impacts on plants and their ecosystems; therefore, a zoomed-in understanding of how plants respond to herbicides at the genomic and epigenomic level will give us perspective on the evolutionary consequences of weed control in modern agriculture. 1 Introduction Weeds are a major threat to agricultural systems worldwide. Left uncontrolled, it is estimated that yield losses caused in North American corn and soybeans alone would cost $43 billion annually (Soltani et al., 2016, 2017). Herbicides are an effective tool for control, especially when paired with herbicide resistant crops (Swinton & Van Deynze, 2017), and growers have become heavily reliant on them for weed management. However, the continual use of the same chemistries has resulted in herbicide resistance developing in 266 weed species (Heap, 2022). Glyphosate is of particular importance because it has become the most used herbicide in the world (Duke & Powles, 2008). It is favored for its efficacy, low cost, relatively low environmental impact, and ease of use (Swinton & Van Deynze, 2017; Heap & Duke, 2018). Although it was first introduced in the 1970’s, genetically engineered Roundup Ready crops caused its application to increase 15-fold (Benbrook, 2016). Consequently, there are now 55 glyphosate resistant weed species globally (Heap, 2022). The continued selection pressure allows the most resistant plants to survive and produce progeny, worsening the problem year after year (Concenço, 2016). In turn, glyphosate resistance decreases the value of the Roundup Ready technology, forcing the agriculture industry to develop new chemistries, revert to older chemistries, adopt different cultural management strategies, and develop new herbicide resistance traits to other chemistries (Van Deynze et al., 2022). Unfortunately, there have been no new herbicide sites of action (specific process in the plant the herbicide targets) released in over 30 years and the tools available for weed management are now dwindling (Heap & Duke, 2018). Understanding how resistance evolves in the field is essential to maintaining the efficacy of herbicides and the productivity of our agricultural systems. An understanding of weed biology and this rapid 2 evolution occurring in the field will help identify ways to address this problem that has both economic and environmental costs. Evolutionary theory and herbicide resistance Evolution is often thought of as both 1) a gradual process that takes multiple generations of small selection pressures slightly altering fitness each generation, and 2) as periods of microevolution where rare but abrupt adaptations take place within a population due to high selection pressure (Gressel, 2009; Koonin & Wolf, 2009). In the case of herbicide resistance, the latter way of thought is usually applied due to the extreme selection pressure imposed by herbicides and the large population size of most weeds (Busi et al., 2013). Herbicide resistance has been documented to evolve as rapidly as three generations (Busi et al., 2012; Busi & Powles, 2009). The origin of resistance alleles in herbicide resistant populations has not been deeply studied in many species, apart from in Amaranthus tuberculatus (Kreiner et al., 2022), but we believe that modern omics approaches will allow us to untangle the origin of novel resistance mechanisms in weed species. From standard evolutionary thinking, we assume resistance alleles exist because of natural variation in any given weed population (Hawkins et al., 2019). When herbicides are sprayed across a field, individuals with resistance traits are selected; they survive, reproduce, and pass resistance alleles to their offspring (Concenço, 2016). The emerging seed bank contains more resistant individuals than prior to selection and the change in resistance allele frequency defines microevolution from a population genetics viewpoint (Gressel, 2009; Hendry & Kinnison, 2001). If there is a fitness penalty associated with resistance alleles, such as reduced reproductive ability, the alleles will eventually be purged from the population over time in the absence of herbicide selection (Vila-Aiub, 2019). On the other hand, continual stress of the same 3 herbicide keeps the resistance alleles in the population at high frequencies (Hawkins et al., 2019). Modern evolutionary theory assumes novel variation is not biased towards increasing fitness in a given environmental condition (Svensson & Berger, 2019). In plants, this is not always the case. Various forms of genetic variation (e.g., single nucleotide polymorphisms (SNPs) and insertion/deletions (in/dels)) appear at a fixed rate in the genome (Bromham et al., 2015; Lynch, 2010), while other forms of variation (structural variants (SV), transposable element (TE) activity, chromosome cross-over) may be punctuated based on environmental stresses or other stimuli (Wellenreuther et al., 2019). The theory that the environment shapes an organism's phenotype, and that in turn offspring inherit these phenotypes, was first proposed by Jean Baptiste Lamarck. However, Darwin’s theory has overshadowed Lamarck's (Day & Bonduriansky, 2011). With modern understanding of genetics and epigenetics, we now know that Lamarck is, at least in part, correct: the environment can affect a phenotype through epigenetic modifications and genomic SV, and be heritable. The question then becomes, how much phenotypic variability can be accounted for by random mutation versus being induced systematically by environmental stress? When we apply this question to the origin of herbicide resistance alleles, it is still unclear the relative contribution of each to phenotypic variability. We suggest a combination of both scenarios is taking place: 1) Herbicide application selects resistant individuals that exist in the population due to natural variation; and 2) Herbicide application acts as an environmental stress, changing plant genomes, which in turn increases genetic diversity, allowing weeds to adapt faster to herbicide application. We want to make a point to differentiate our focus here from other reviews on related topics, which have hinted at epigenetic involvement in stress response, and highlight the mechanisms by which this can occur more in depth (i.e., transposon activation, 4 gene amplification, structural variants) (Dyer, 2018; Markus et al., 2018). There is a gap in knowledge of the extent to which the epigenome contributes to rapid evolution of weedy species, the mechanisms that drive this, and their consequences under herbicide stress. In this review, we will walk through various discoveries in stress biology, resistance evolution, genome biology and epigenetics in the context of herbicides. When put together, they make a compelling hypothesis for the contribution of the environment to the evolution of novel variation and herbicide resistance. We will describe how a macroscale event, spraying herbicides over a field to control weeds, can have unintended microscale effects, changing the genome and epigenome of an individual plant, profoundly impacting resistance evolution and agriculture. How herbicides cause stress and its perception through the plant Plants are sedentary organisms that survive through fluctuating and often adverse environmental conditions by changing their molecular physiology (Laitinen & Nikoloski, 2019). These environmental stresses include ultraviolet radiation, drought, salinity, heavy metals, and extreme temperatures (J. K. Zhu, 2016). A series of signaling networks and plant hormone responses enable the plant to maintain homeostasis. One important aspect of plant stress response are reactive oxygen species (ROS), which are a natural byproduct of metabolism in the plant (J. K. Zhu, 2016). Under normal conditions, ROS function as signaling molecules and are scavenged and detoxified immediately to prevent harm to cellular functions. However, under stress conditions, partially reduced oxygen intermediates can build up creating an oxidative burst and, potentially, cell death (Saxena et al., 2016). This phenomenon is not exclusive to the environmental stresses listed above, as most herbicides cause overproduction of ROS that either directly or indirectly kill the plant (Caverzan et al., 2019). Enzymes that buffer this process include super oxide dismutase (SOD), peroxidase (POX), glutathione S-transferases (GSTs) and 5 glutathione reductase (GR). Both GSTs and GR participate in detoxifying xenobiotics including herbicides (Katerova & Miteva, 2010; Thom et al., 2002). Cytochrome P450s (CYP450s) are also important enzymes in metabolic pathways and 30 different CYP450s have been linked to herbicide detoxification and metabolism (Dimaano & Iwakami, 2021). Herbicide resistance mechanisms and their origins Herbicide resistance mechanisms can be classified as target-site resistance (TSR) or non- target-site resistance (NTSR), with cases where both occur in the same individual (Gaines et al., 2020). TSR includes mutations that alter the sequence of the target enzyme and impair the herbicide’s ability to bind (Powles & Yu, 2010), and is primarily monogenic (Délye, Jasieniuk, et al., 2013). The acetyl-CoA carboxylase or ACCase inhibitors (WSSA Group 1) are used to control grasses, and multiple weed species have evolved TSR at different point mutations (Murphy & Tranel, 2019). Studies have shown that these point mutations were in populations of Alopecurus myosuroides long before the introduction of ACCase chemistries at high rates (Délye, Deulvot, et al., 2013). Furthermore, the relative frequency of resistance mutations appears to coincide with the selection pressure of the ACCase herbicides used (Kaundun, 2014). Resistance to the acetolactate synthase (ALS) inhibitor (Group 2) herbicides is also commonly due to target-site mutations. Other cases of TSR through point mutations have been reported in the microtubule inhibitors (Group 3), synthetic auxins (Group 4), photosystem II inhibitors (Groups 5-7), 5-enolpyruvyl-shikimate-3-phosphate synthase (EPSPS) inhibitor (Group 9), glutamine synthetase inhibitors (Group 10), phytoene desaturase inhibitors (Group 12) and protoporphyrinogen oxidase (PPO) inhibitors (Group 14), among others (Murphy & Tranel, 2019). 6 Another mechanism of TSR is target site duplication, where there are multiple copies of the target gene causing over expression of the target protein (Gaines et al., 2010). So far, in weeds, this mechanism has only been shown to confer resistance to glyphosate (Sammons & Gaines, 2014). Glyphosate inhibits the EPSPS enzyme, which prevents the catalyzation of the reaction needed to form aromatic amino acids essential to plant growth (Duke & Powles, 2008). Glyphosate resistance in kochia (Bassia scoparia (L.) A.J. Scott syn. Kochia scoparia (L.) Schrad.) is caused by having as little as two to more than ten copies of the EPSPS gene that are aligned in tandem, where the increase in expression means field rates of glyphosate are not sufficient for control (Patterson et al., 2019). Remarkably, there are at least eight other weed species that have shown copy number variation (CNV) of the EPSPS gene conferring glyphosate resistance as well, each with their own seemingly novel structural variant (Patterson et al., 2018). Of these, most remarkably, is the case of Amaranthus palmeri where EPSPS duplication is facilitated by a large molecule of extrachromosomal circular DNA (eccDNA) that contains up to 120 copies of the gene (Koo et al., 2018)! Not only can this gene amplification confer glyphosate resistance, but it also seemed to protect A. palmeri from the oxidative stress associated with herbicide response, due to the incomplete inhibition of EPSPS (Eceiza et al., 2022). Further, it is intriguing that more than eight species have converged on EPSPS CNV as the primary resistance mechanism to one herbicide, though each achieve this goal through novel rearrangements. We believe this indicates something unique about glyphosate and the way plants respond to glyphosate selection (Patterson et al., 2018). On the other hand, NTSR mechanisms include reduced translocation and uptake of the herbicide, sequestration, metabolism, and reduced activation of the herbicide (Jugulam & Shyam, 2019). It is thought to be more common in grasses, but there is still evidence of NTSR 7 mechanisms in broadleaves as well (Scarabel et al., 2015). In comparison with TSR, NTSR can be, and often is, multigenic (Délye, Jasieniuk, et al., 2013). The genes involved are part of general stress response pathways. allowing for cross resistance to multiple chemistries, including CYP450s, GSTs, and translocation due to ATP-binding cassette (ABC) transporters (Jugulam & Shyam, 2019). Interestingly, weeds have adapted a myriad of glyphosate resistance mechanisms, ranging from target site point mutations and copy number variants to NTSR mechanisms such as reduced translocation and metabolism, and even sequestration in the vacuole (Pan et al., 2019; Sammons & Gaines, 2014). In both TSR and NTSR, de novo mutations are more likely to occur than gene flow from other populations, until the allele frequency reaches a sufficient level (Délye, Jasieniuk, et al., 2013). Recurrent selection with sublethal doses selects for minor-effect resistance genes consistent with NTSR mechanisms (Busi & Powles, 2009; Tehranchian et al., 2017). This may occur when there is drift from a nearby field, poor coverage from a sprayer, or the weeds were larger than what is recommended at application time (Tehranchian et al., 2017). For example, one study on the herbicide pyroxasulfone showed resistance developing in a susceptible population in only three generations of sublethal dose application in Lolium rigidum, with a similar finding under glyphosate application (Busi et al., 2012; Busi & Powles, 2009). In contrast, TSR is the result of major-effect resistance genes that undergo selective sweep adaptation, where favorable mutations continue in a population (Délye, Jasieniuk, et al., 2013). This is of course, dependent on the mutation rate, initial allele frequencies, fitness costs, reproduction, and gene flow, which can differ between species, as reviewed here (Vila-Aiub, 2019). 8 Assessing structural variation as an evolutionary tool for adaptation Structural variation is a major source of genetic novelty. Structural variants encompass almost any significant changes in genomic architecture such as large inversions, translocations, insertions, deletions, and CNVs (Wellenreuther et al., 2019). They usually occur due to errors in recombination, DNA double strand breaks, and replication errors during either mitosis (somatic) or meiosis (germline) (Gabur et al., 2019). CNVs of functional genes have the potential to alter phenotypes and change the evolutionary trajectory of a species, from humans to insects, and even plants (Feuk et al., 2006; Gorkovskiy & Verstrepen, 2021; Vollmer, 2008; Wellenreuther et al., 2019). In humans, CNVs have been linked to an array of diseases and neurological disorders such as Alzheimer’s disease, asthma, and cancer (Almal & Padh, 2012). In insects, such as aphids and Drosophila, insecticide resistance has been shown to be caused by CNVs, and in plants, CNVs help prevent insect herbivory (Broekgaarden et al., 2011; ffrench-Constant, 2013). Mosquito control for malaria is also threatened by the development of insecticide resistance through CNVs (Weetman et al., 2018). Research has shown that genomic changes such as single nucleotide polymorphisms (SNPs) as well as CNVs can be environmentally induced (Guerrero-Bosagna, 2020). This has been demonstrated in the model plant Arabidopsis thaliana, with higher rates of CNV’s occurring within generations under temperature and salicylic acid stresses (Debolt, 2010). Structural variation has become an important topic in crop species as a source of new traits (Gabur et al., 2019). The cultivated potato, Solanum tuberosum, displays extensive CNVs in regions of greater than 100 kb concentrated near stress tolerant loci (Hardigan et al., 2015). Uncovering this diversity is a key motivation for the development of crop pangenome projects, in which there are both core genes and dispensable genes that are often lost or gained in response 9 to environmental stresses (Golicz et al., 2020). To this extent, CNVs are being extensively characterized in depth in Oryza sativa, Zea mays, Hordeum vulgare, Triticum aestivum and Glycine max (Cook et al., 2012; Juery et al., 2021; Manching et al., 2017; Muñoz-Amatriaín et al., 2013; Yu et al., 2013). These new CNVs are shown to contribute to phenotypes like disease resistance, defense responses, drought tolerance, and a wide range of other environmental responses that could be useful in the context of global climate change threatening our food supply (Gabur et al., 2019). In addition to CNVs, transposable elements are a significant source of genetic variation, and their epigenetic regulation is of foremost importance (Weil & Martienssen, 2008). Transposons that “copy and paste” can proliferate in high copies when not properly regulated through silencing by the epigenome (Henderson & Jacobsen, 2007). The tomato panSV (structural variants) genome, created from 100 different accessions, demonstrates that structural variants are responsible for major improvements in fruit quality, and frequently have underlying transposons (Alonge et al., 2020). An important aspect of gene duplication is also considering what genes are co-duplicated in CNV events and how they also contribute to phenotypic changes (Patterson et al., 2019). The kochia EPSPS gene duplication event was due to a combination of these, beginning with a transposon insertion upstream and downstream of the EPSPS locus, causing unequal crossing over followed by a double strand break. In the case of glyphosate resistance in kochia, six genes are conserved within the repeat region associated with EPSPS duplication, with only four having a linear increase in expression with the number of repeats (Patterson et al., 2019). More research is needed to decipher genes being co-duplicated and their relationship with plant phenotype. Specifically, epigenetic variation of DNA methylation could serve as a mechanism of 10 silencing or controlling the harmful effects of genome instability induced by stress (Seymour & Becker, 2017). Specifically, the convergent evolution of glyphosate resistance due to EPSPS CNV in multiple weed species has scientists perplexed: why are CNVs the mechanism for glyphosate resistance but no other herbicide, apart from an isolated case of ACCase resistance in large crabgrass in 2017 (Laforest et al., 2017). Unraveling the role of the epigenome Plants have countless mechanisms to respond to stress that are broadly grouped as tolerance, avoidance, or resistance (Boyko & Kovalchuk, 2008). However, plant species also can rapidly acclimate to stress in their environment at short trans-generational timescales, sometimes even within a single generation (Boyko et al., 2010). Even more surprising, it can be reversible. Epigenetic modifications allow for heritable changes in gene expression and therefore environmental stress response (Mirouze & Paszkowski, 2011). Epigenetic modifications generally either increase or decrease transcriptional machinery’s access to genes, modulating their expression. These include changes to histones, proteins responsible for packaging DNA and changes in DNA methylation (Lämke & Bäurle, 2017). DNA methylation is the transfer of a methyl group to a cytosine at the C5 position, forming 5-methylcytosine (Zhang et al., 2018). It has particular importance in stress response and has been shown to have critical roles in development, stress protection, and regulating the activation of transposable elements (Boyko & Kovalchuk, 2008). DNA methylation occurs in three different contexts: CG, CHG, and CHH (where H = A, T, or C), each with their own significance and maintenance system. CG methylation (mCG) is the most common form in both plants and animals and is maintained by DNA Methyltransferase 1 (MET1) after DNA replication on both strands, meaning both cytosines will be methylated to 11 maintain symmetry (Niederhuth & Schmitz, 2017). CHG methylation is maintained through a different pathway, involving chromomethylase 3 (CMT3), and is also symmetrical (Weil & Martienssen, 2008). The final context, CHH methylation, is asymmetrical and therefore must be established de novo after each DNA replication (Niederhuth & Schmitz, 2017). De novo methylation across all three contexts is due to the RNA-directed DNA Methylation (RdDM) pathway, which produces small-interfering RNA (siRNA) that in turn recruit DNA methyltransferases for transcriptional silencing at specific sites (Matzke & Mosher, 2014; Weil & Martienssen, 2008). The frequency of each methylation context throughout the genome varies by species, with 24% of CGs methylated in Arabidopsis compared to 52% in Oryza sativa (Markus et al., 2018). For Beta vulgaris, a close relative of kochia, cytosine methylation levels by sequence context were reported as 66.8% for CG, 41.7% for CHG, and only 9.4% for CHH (Gutschker et al., 2022). Intriguingly, methylated cytosines are subject to higher mutation rates, especially from C to T, and increase the mutability of neighboring nucleotides (Kusmartsev et al., 2020). The three methylation contexts are also specific to their location in the genome. Gene body methylation (gbM) is characterized by only having mCG sites, with depletion around transcriptional start sites (TSS) and transcription termination sites (TTS) (Bewick et al., 2017). This is characteristic of genes, while high mCG near the TSS will prevent transcription (Niederhuth et al., 2016). If there is CHG and CHH methylation in the gene body, it is a target of RdDM and will be transcriptionally silenced. CHH methylation is commonly found around transposons, repetitive elements, and enforces the boundary between euchromatin and heterochromatin (Niederhuth et al., 2016). In Zea mays, it has been shown that “CHH islands,” or places with CHH methylation upstream and downstream of a gene result in an increase in 12 gene expression (Li et al., 2015). Methylation at transcription factor (TF) binding sites prevents their binding in most cases (O’Malley et al., 2016). Demethylation is regulated by REPRESSOR OF SILENCING 1 (ROS1), a DNA glycosylase, especially at TF binding sites. ROS1 also plays a role in environmental responses and development (J. Zhu et al., 2007). In Arabidopsis, ROS1 is proposed to be involved in regulating genes associated with herbicide response, as ros1 mutants were more susceptible to imazethapyr than wild type plants (Markus et al., 2021). Therefore, epigenetic regulations must be tied with herbicide detoxification. Transposon activation is a common result of stress response, in which demethylation allows their excision. When DNA is demethylated under stress, gene expression patterns change and TEs are activated, allowing the plant to acclimate quickly and in novel ways to the stress (Boyko & Kovalchuk, 2008). In Arabidopsis for example, researchers were able to link transcriptomic data with loci affected by transcriptional gene silencing (TGS) following stress treatment (salt, temperature, UV radiation, and water) (Mirouze & Paszkowski, 2011). They also demonstrated the stability of such epigenetic markers, as they were silenced again two days after the response, showing the reversibility of methylation changes. Most epigenome studies have focused on the transient effects, as opposed to heritable, transgenerational changes. This is the case for studies that have shown the effect of herbicides on the epigenome in both model species and crop plants. It has been shown that glyphosate application induces global hypermethylation in T. aestivum and reduces the genomic template stability with increasing doses, meaning there was more variability in the genomic DNA (Nardemir et al., 2015). In contrast, sublethal glyphosate application did not change global methylation levels in Arabidopsis, it showed a dose-dependent response in context-specific 13 methylation (CG/CHG) compared between 5% and 10% field rate applications (Kim et al., 2017). Interestingly, glyphosate application does not result in methylation of EPSPS, which is constitutively unmethylated in Arabidopsis. However, researchers did find differentially methylated regions (DMRs) that were unique to glyphosate response and not shared with phosphate starvation or biotic stress. Hypomethylation of TEs was also observed, which leaves them vulnerable to activation (Kim et al., 2017). In O. sativa, atrazine application changed global methylation, where specific genes were either hyper- or hypomethylated, resulting in changes in expression of detoxification genes. They also identified differential methylation and expression in transposable elements (Lu et al., 2016). Still, little is known about how herbicides effect the epigenome of weedy species. The role of epigenetics was hinted at when researchers discovered Eleusine indica resistance to dinitroanilines was caused by a mutation of a methylated cytosine (Anthony et al., 1998). Though more recently, preliminary data in Conyza canadensis reported increased gene body methylation of EPSPS in plants resistant to glyphosate compared to sensitive biotypes, explaining a possible mechanism for increased resistance (Margaritopoulou et al., 2018). These studies have focused on the possibility of transient changes in response to herbicides and focused on the relationship to NTSR. We also need to identify whether epigenetic modifications have directly contributed to SV and their relationship to directed, heritable evolution of herbicide TSR, namely through copy number variation. Transgenerational epigenetic studies in plants are limited with no information on weed species; but there is evidence that the progeny of stressed plants experience heritable epigenetic modifications. In the case of Nicotiana tabacum infected with tobacco mosaic virus (TMV), progeny was shown to be hypermethylated, apart from loci containing leucine-rich-repeats 14 (LRR), which were hypomethylated (Boyko et al., 2007). These LRR containing loci also had a higher frequency of genomic rearrangements, highlighting the importance of epigenetic regulation in silencing and activating expression. The progeny of Arabidopsis subject to abiotic (drought, flood, heat, salt and UVC) stresses exhibited increased homologous recombination frequency, DNA hypermethylation, as well as increased tolerance to the same stress. In this case, however the response did not persist in future generations (Boyko & Kovalchuk, 2011). Through multi-omics approaches, researchers have been able to link the environmental stress imposed on a plant to differential phenotypes. We can now observe how the epigenome, transcriptome, and genome are affected by the environment, and can prove they are dynamic in their response. Bringing it all together: The future of studying the origins of resistance Environmental stresses (including herbicides), change the epigenome, which in turn changes mutation rates (especially SV) (Pareek et al., 2010). This novel variation may include herbicide resistance traits that can be selected for and proliferate. If this is true, then resistance may be already in the field or be induced by the act of applying herbicides. If resistance evolves via the epigenome, it could also be reversed or silenced. For example, the first example of the role of epigenetics in insecticide resistance was demonstrated through gene amplification that caused resistance but was silenced via the epigenome in the absence of selection to prevent a fitness penalty (Oppold & Müller, 2017). Therefore, we suggest the need for studying the relative contribution of natural genetic variation versus induced variation in herbicide resistance evolution and call for a more unified understanding of resistance evolution which includes both Darwinian and Lamarckian concepts, including epigenomics. 15 Conclusion With the reduction in cost of next generation sequencing technology, more “omics” data can be generated than ever before. It is now possible to examine the relationship of the epigenome with genomic architectural changes under different conditions on non-model or crop species. The field of computational genomics and epigenomics opens a range of opportunities for advancements in weed science. Understanding the biology of weedy genomes and their innate ability to adapt to abiotic stresses will open the door to innovative solutions to combat herbicide resistance. After all, “growers will not be able to spray their way out of a situation they sprayed their way into.” Not only is this research valuable in understanding herbicide resistance evolution, but it will also help unravel the contributions of genomic SV to adaptation across a rapid timescale. These questions are of key interest to evolutionary biologists and can help us interpret comprehensive plant responses to their environment (Neve et al., 2009). Researchers have called for increased priority on epigenetics studies for decades, with the same evolutionary questions in mind: 1) how does the environment influence de novo generation of epialleles and 2) how are they inherited in a population? Which directly applies to what weed biologists are attempting to understand. Just as pan-genomes may become the norm when it comes to genomic resources, it is no longer enough to understand plant genetics without the context of epigenetic mechanisms and how they influence variation in gene expression and phenotypic plasticity. The contribution of the epigenome to the evolution of weedy species needs to be a continued area of research focus. The world’s food, feed, fiber, and fuel supplies depend on having access to effective pest management tools, especially herbicides. However, we must consider the consequences of anthropogenic influences on the world around us. This research will help us to understand and 16 predict the threat of future herbicide resistance evolution, in different areas and to new herbicides, as we work to protect this valuable management tool. We want to preserve the efficacy of both the herbicides and associated genetically modified crops that allow them to be used in a sustainable way. 17 REFERENCES 18 REFERENCES Almal, S., & Padh, H. (2012). Implications of gene copy-number variation in health and diseases. Journal of Human Genetics, 57(1), 6–13. Alonge, M., Wang, X., Benoit, M., Soyk, S., Pereira, L., Zhang, L., Suresh, H., Ramakrishnan, S., Maumus, F., Ciren, D., Levy, Y., Harel, T. H., Shalev-Schlosser, G., Amsellem, Z., Razifard, H., Caicedo, A. L., Tieman, D. M., Klee, H., Kirsche, M., … Lippman, Z. B. (2020). Major Impacts of Widespread Structural Variation on Gene Expression and Crop Improvement in Tomato. Cell, 182(1), 145-161.e23. Anthony, R. G., Waldin, T. R., Ray, J. A., Bright, S. W. J., & Hussey, P. J. (1998). Herbicide resistance caused by spontaneous mutation of the cytoskeletal protein tubulin. Nature 1998 393:6682, 393(6682), 260–263. Benbrook, C. M. (2016). Trends in glyphosate herbicide use in the United States and globally. Environmental Sciences Europe, 28(1), 3. Bewick, A. J., Niederhuth, C. E., Ji, L., Rohr, N. A., Griffin, P. T., Leebens-Mack, J., & Schmitz, R. J. (2017). The evolution of CHROMOMETHYLASES and gene body DNA methylation in plants. Genome Biology, 18(1), 65. Boyko, A., Blevins, T., Yao, Y., Golubov, A., Bilichak, A., Ilnytskyy, Y., Hollander, J., Meins, F., & Kovalchuk, I. (2010). Transgenerational Adaptation of Arabidopsis to Stress Requires DNA Methylation and the Function of Dicer-Like Proteins. PLOS ONE, 5(3), e9514. Boyko, A., Kathiria, P., Zemp, F. J., Yao, Y., Pogribny, I., & Kovalchuk, I. (2007). Transgenerational changes in the genome stability and methylation in pathogen-infected plants (Virus-induced plant genome instability). Nucleic Acids Research, 35(5), 1714– 1725. Boyko, A., & Kovalchuk, I. (2008). Epigenetic control of plant stress response. Environmental and Molecular Mutagenesis, 49(1), 61–72. Boyko, A., & Kovalchuk, I. (2011). Genome instability and epigenetic modification — heritable responses to environmental stress? Current Opinion in Plant Biology, 14(3), 260–266. Broekgaarden, C., Snoeren, T. A. L., Dicke, M., & Vosman, B. (2011). Exploiting natural variation to identify insect-resistance genes. Plant Biotechnology Journal, 9(8), 819–825. Bromham, L., Hua, X., Lanfear, R., & Cowman, P. F. (2015). Exploring the relationships between mutation rates, life history, genome size, environment, and species richness in flowering plants. American Naturalist, 185(4), 508–524. 19 Busi, R., Gaines, T. A., Walsh, M. J., & Powles, S. B. (2012). Understanding the potential for resistance evolution to the new herbicide pyroxasulfone: Field selection at high doses versus recurrent selection at low doses. Weed Research, 52(6), 489–499. Busi, R., & Powles, S. B. (2009). Evolution of glyphosate resistance in a Lolium rigidum population by glyphosate selection at sublethal doses. Heredity 2009 103:4, 103(4), 318– 325. Busi, R., Vila-Aiub, M. M., Beckie, H. J., Gaines, T. A., Goggin, D. E., Kaundun, S. S., Lacoste, M., Neve, P., Nissen, S. J., Norsworthy, J. K., Renton, M., Shaner, D. L., Tranel, P. J., Wright, T., Yu, Q., & Powles, S. B. (2013). Herbicide-resistant weeds: from research and knowledge to future needs. Evolutionary Applications, 6(8), 1218–1221. Caverzan, A., Piasecki, C., Chavarria, G., Stewart, C. N., & Vargas, L. (2019). Defenses Against ROS in Crops and Weeds: The Effects of Interference and Herbicides. International Journal of Molecular Sciences 2019, Vol. 20, Page 1086, 20(5), 1086. Concenço, G. (2016). Evolution, epigenetics and resistance – troublesome weeds. Revista Brasileira de Herbicidas, 15(1), 14. Cook, D. E., Lee, T. G., Guo, X., Melito, S., Wang, K., Bayless, A. M., Wang, J., Hughes, T. J., Willis, D. K., Clemente, T. E., Diers, B. W., Jiang, J., Hudson, M. E., & Bent, A. F. (2012). Copy number variation of multiple genes at Rhg1 mediates nematode resistance in soybean. Science, 338(6111), 1206–1209. Day, T., & Bonduriansky, R. (2011). A unified approach to the evolutionary consequences of genetic and nongenetic inheritance. American Naturalist, 178(2). Debolt, S. (2010). Copy Number Variation Shapes Genome Diversity in Arabidopsis Over Immediate Family Generational Scales. Genome Biology and Evolution, 2(1), 441–453. Délye, C., Deulvot, C., & Chauvel, B. (2013). DNA Analysis of Herbarium Specimens of the Grass Weed Alopecurus myosuroides Reveals Herbicide Resistance Pre-Dated Herbicides. PLoS ONE, 8(10). Délye, C., Jasieniuk, M., & le Corre, V. (2013). Deciphering the evolution of herbicide resistance in weeds. Trends in Genetics, 29(11), 649–658. Dimaano, N. G., & Iwakami, S. (2021). Cytochrome P450-mediated herbicide metabolism in plants: current understanding and prospects. Pest Management Science, 77(1), 22–32. Duke, S. O., & Powles, S. B. (2008). Glyphosate: a once-in-a-century herbicide. Pest Management Science, 64(4), 319–325. Dyer, W. E. (2018). Stress-induced evolution of herbicide resistance and related pleiotropic effects. Pest Management Science, 74(8), 1759–1768. 20 Eceiza, M. V., Gil-Monreal, M., Barco-Antoñanzas, M., Zabalza, A., & Royuela, M. (2022). The moderate oxidative stress induced by glyphosate is not detected in Amaranthus palmeri plants overexpressing EPSPS. Journal of Plant Physiology, 153720. Feuk, L., Carson, A. R., & Scherer, S. W. (2006). Structural variation in the human genome. Nature Reviews Genetics 2006 7:2, 7(2), 85–97. ffrench-Constant, R. H. (2013). The Molecular Genetics of Insecticide Resistance. Genetics, 194(4), 807–815. Gabur, I., Chawla, H. S., Snowdon, R. J., & Parkin, I. A. P. (2019). Connecting genome structural variation with complex traits in crop plants. Theoretical and Applied Genetics, 132(3), 733–750. Gaines, T. A., Duke, S. O., Morran, S., Rigon, C. A. G., Tranel, P. J., Anita Küpper, & Dayan, F. E. (2020). Mechanisms of evolved herbicide resistance. Journal of Biological Chemistry, 295(30), 10307–10330. Gaines, T. A., Zhang, W., Wang, D., Bukun, B., Chisholm, S. T., Shaner, D. L., Nissen, S. J., Patzoldt, W. L., Tranel, P. J., Culpepper, A. S., Grey, T. L., Webster, T. M., Vencill, W. K., Sammons, R. D., Jiang, J., Preston, C., Leach, J. E., & Westra, P. (2010). Gene amplification confers glyphosate resistance in Amaranthus palmeri. Proceedings of the National Academy of Sciences, 107(3), 1029–1034. Golicz, A. A., Bayer, P. E., Bhalla, P. L., Batley, J., & Edwards, D. (2020). Pangenomics Comes of Age: From Bacteria to Plant and Animal Applications. Trends in Genetics, 36(2), 132– 145. Gorkovskiy, A., & Verstrepen, K. J. (2021). The Role of Structural Variation in Adaptation and Evolution of Yeast and Other Fungi. Genes 2021, Vol. 12, Page 699, 12(5), 699. Gressel, J. (2009). Evolving understanding of the evolution of herbicide resistance. Pest Management Science, 65(11), 1164–1173. Guerrero-Bosagna, C. (2020). From epigenotype to new genotypes: Relevance of epigenetic mechanisms in the emergence of genomic evolutionary novelty. In Seminars in Cell and Developmental Biology (97). Gutschker, S., Corral, J. M., Schmiedl, A., Ludewig, F., Koch, W., Fiedler-Wiechers, K., Czarnecki, O., Harms, K., Keller, I., Martins Rodrigues, C., Pommerrenig, B., Neuhaus, H. E., Zierer, W., Sonnewald, U., & Müdsam, C. (2022). Multi-omics data integration reveals link between epigenetic modifications and gene expression in sugar beet (Beta vulgaris subsp. vulgaris) in response to cold. BMC Genomics, 23(1). Hardigan, M. A., Crisovan, E., Hamilton, J. P., Kim, J., Laimbeer, P., Leisner, C. P., Manrique 21 Carpintero, N. C., Newton, L., Pham, G. M., Vaillancourt, B., Yang, X., Zeng, Z., Douches, D. S., Jiang, J., Veilleux, R. E., & Buella, C. R. (2015). Genome reduction uncovers a large dispensable genome and adaptive role for copy number variation in asexually propagated Solanum tuberosum. Plant Cell, 28(2), 388–405. Hawkins, N. J., Bass, C., Dixon, A., & Neve, P. (2019). The evolutionary origins of pesticide resistance. Biological Reviews, 94(1), 135–155. Heap, I. (2022). The International Herbicide-Resistant Weed Database. www.weedscience.org Heap, I., & Duke, S. O. (2018). Overview of glyphosate-resistant weeds worldwide. Pest Management Science, 74(5), 1040–1049. Henderson, I. R., & Jacobsen, S. E. (2007). Epigenetic inheritance in plants. Nature. 447(7143), 418–424. Hendry, A. P., & Kinnison, & M. T. (2001). An introduction to microevolution: rate, pattern, process. Genetica, 112, 1–8. Juery, C., Concia, L., de Oliveira, R., Papon, N., Ramírez-González, R., Benhamed, M., Uauy, C., Choulet, F., & Paux, E. (2021). New insights into homoeologous copy number variations in the hexaploid wheat genome. The Plant Genome, 14(1), e20069. Jugulam, M., & Shyam, C. (2019). Non-Target-Site Resistance to Herbicides: Recent Developments. Plants, 8(10), 417. Katerova, Z. I., & Miteva, L. P. E. (2010). Glutathione and herbicide resistance in plants. Ascorbate-Glutathione Pathway and Stress Tolerance in Plants, 191–207. Kaundun, S. S. (2014). Resistance to acetyl-CoA carboxylase-inhibiting herbicides. Pest Management Science, 70(9), 1405–1417. Kim, G., Clarke, C. R., Larose, H., Tran, H. T., Haak, D. C., Zhang, L., Askew, S., Barney, J., & Westwood, J. H. (2017). Herbicide injury induces DNA methylome alterations in Arabidopsis. PeerJ, 2017(7), 3560. Koo, D. H., Molin, W. T., Saski, C. A., Jiang, J., Putta, K., Jugulam, M., Friebe, B., & Gill, B. S. (2018). Extrachromosomal circular DNA-based amplification and transmission of herbicide resistance in crop weed Amaranthus palmeri. Proceedings of the National Academy of Sciences of the United States of America, 115(13), 3332–3337. Koonin, E. v, & Wolf, Y. I. (2009). Is evolution Darwinian or/and Lamarckian? Biology Direct, 4(1), 42. Kreiner, J. M., Sandler, G., Stern, A. J., Tranel, P. J., Weigel, D., Stinchcombe, J. R., & Wright, 22 S. I. (2022). Repeated origins, widespread gene flow, and allelic interactions of target-site herbicide resistance mutations. ELife, 11. Kusmartsev, V., Drozdz, M., Schuster-Böckler, B., & Warnecke, T. (2020). Cytosine Methylation Affects the Mutability of Neighboring Nucleotides in Germline and Soma. Genetics, 214(4), 809–823. Laforest, M., Soufiane, B., Simard, M. J., Obeid, K., Page, E., & Nurse, R. E. (2017). Acetyl CoA carboxylase overexpression in herbicide-resistant large crabgrass (Digitaria sanguinalis). Pest Management Science, 73(11), 2227–2235. Laitinen, R. A. E., & Nikoloski, Z. (2019). Genetic basis of plasticity in plants. Journal of Experimental Botany, 70(3), 739–745. Lämke, J., & Bäurle, I. (2017). Epigenetic and chromatin-based mechanisms in environmental stress adaptation and stress memory in plants. Genome Biology, 18(1), 124. Li, Q., Gent, J. I., Zynda, G., Song, J., Makarevitch, I., Hirsch, C. D., Hirsch, C. N., Dawe, R. K., Madzima, T. F., McGinnis, K. M., Lisch, D., Schmitz, R. J., Vaughn, M. W., & Springer, N. M. (2015). RNA-directed DNA methylation enforces boundaries between heterochromatin and euchromatin in the maize genome. Proceedings of the National Academy of Sciences of the United States of America, 112(47), 14728–14733. Lu, Y. C., Feng, S. J., Zhang, J. J., Luo, F., Zhang, S., & Yang, H. (2016). Genome-wide identification of DNA methylation provides insights into the association of gene expression in rice exposed to pesticide atrazine. Scientific Reports 2016 6:1, 6(1), 1–15. Lynch, M. (2010). Evolution of the mutation rate. Trends in Genetics, 26(8), 345–352. Manching, H., Sengupta, S., Hopper, K. R., Polson, S. W., Ji, Y., & Wisser, R. J. (2017). Phased genotyping-by-sequencing enhances analysis of genetic diversity and reveals divergent copy number variants in maize. G3: Genes, Genomes, Genetics, 7(7), 2161–2170. Margaritopoulou, T., Tani, E., Chachalis, D., & Travlos, I. (2018). Involvement of Epigenetic Mechanisms in Herbicide Resistance: The Case of Conyza canadensis. Agriculture 2018, Vol. 8, Page 17, 8(1), 17. Markus, C., Pecinka, A., Karan, R., Barney, J. N., & Merotto, A. (2018). Epigenetic regulation – contribution to herbicide resistance in weeds? Pest Management Science, 74(2), 275–281. Markus, C., Pecinka, A., & Merotto, A. (2021). Insights into the Role of Transcriptional Gene Silencing in Response to Herbicide-Treatments in Arabidopsis thaliana. International Journal of Molecular Sciences Article. Matzke, M. A., & Mosher, R. A. (2014). RNA-directed DNA methylation: epigenetic pathway of increasing complexity. Nature Reviews Genetics 2014 15:6, 15(6), 394–408. 23 Mirouze, M., & Paszkowski, J. (2011). Epigenetic contribution to stress adaptation in plants. Current Opinion in Plant Biology, 14(3), 267–274. Muñoz-Amatriaín, M., Eichten, S. R., Wicker, T., Richmond, T. A., Mascher, M., Steuernagel,B., Scholz, U., Ariyadasa, R., Spannagl, M., Nussbaumer, T., Mayer, K. F. X., Taudien, S., Platzer, M., Jeddeloh, J. A., Springer, N. M., Muehlbauer, G. J., & Stein, N. (2013). Distribution, functional impact, and origin mechanisms of copy number variation in the barley genome. Genome Biology, 14(6), 1–17. Murphy, B. P., & Tranel, P. J. (2019). Target-Site Mutations Conferring Herbicide Resistance. Plants, 8(10), 382. Nardemir, G., Agar, G., Arslan, E., & Aygun Erturk, F. (2015). Determination of genetic and epigenetic effects of glyphosate on Triticum aestivum with RAPD and CRED-RA techniques. Theoretical and Experimental Plant Physiology, 27(2), 131–139. Neve, P., Vila‐Aiub, M., & Roux, F. (2009). Evolutionary thinking in agricultural weed management. New Phytologist, 184(4), 783–793. Niederhuth, C. E., Bewick, A. J., Ji, L., Alabady, M. S., Kim, K. do, Li, Q., Rohr, N. A., Rambani, A., Burke, J. M., Udall, J. A., Egesi, C., Schmutz, J., Grimwood, J., Jackson, S. A., Springer, N. M., & Schmitz, R. J. (2016). Widespread natural variation of DNA methylation within angiosperms. Genome Biology 2016 17:1, 17(1), 1–19. Niederhuth, C. E., & Schmitz, R. J. (2017). Putting DNA methylation in context: from genomes to gene expression in plants. Biochimica et Biophysica Acta (BBA) - Gene Regulatory Mechanisms, 1860(1), 149–156. O’Malley, R. C., Huang, S. S. C., Song, L., Lewsey, M. G., Bartlett, A., Nery, J. R., Galli, M., Gallavotti, A., & Ecker, J. R. (2016). Cistrome and Epicistrome Features Shape the Regulatory DNA Landscape. Cell, 165(5), 1280–1292. Oppold, A. M., & Müller, R. (2017). Epigenetics: A Hidden Target of Insecticides. Advances in Insect Physiology, 53, 313–324. Pan, L., Yu, Q., Han, H., Mao, L., Nyporko, A., Fan, L., Bai, L., & Powles, S. (2019). Aldo-keto Reductase Metabolizes Glyphosate and Confers Glyphosate Resistance in Echinochloa colona. Plant Physiology, 181(4), 1519–1534. Pareek, A., Sopory, S. K., Bohnert, H. J., & Govindjee (ed.). (2010). Abiotic Stress Adaptation in Plants. Physiological, Molecular and Genomic Foundation. Photosynthetica, 48(3), 474–474. Patterson, E. L., Pettinga, D. J., Ravet, K., Neve, P., & Gaines, T. A. (2018). Glyphosate Resistance and EPSPS Gene Duplication: Convergent Evolution in Multiple Plant Species. Journal of Heredity, 109(2), 117–125. 24 Patterson, E. L., Saski, C. A., Sloan, D. B., Tranel, P. J., Westra, P., Gaines, T. A., & Gaut, B. (2019). The Draft Genome of Kochia scoparia and the Mechanism of Glyphosate Resistance via Transposon-Mediated EPSPS Tandem Gene Duplication. Genome Biology and Evolution, 11(10). Powles, S. B., & Yu, Q. (2010). Evolution in Action: Plants Resistant to Herbicides. Annual Review of Plant Biology, 61(1), 317–347. Sammons, R., & Gaines, T. (2014). Glyphosate resistance: state of knowledge. Pest Management Science, 70(9), 1367–1377. https://doi.org/10.1002/ps.3743 Saxena, I., Srikanth, S., & Chen, Z. (2016). Cross talk between H2O2 and interacting signal molecules under plant stress response. Frontiers in Plant Science, 7. Scarabel, L., Pernin, F., & Délye, C. (2015). Occurrence, genetic control and evolution of non target-site based resistance to herbicides inhibiting acetolactate synthase (ALS) in the dicot weed Papaver rhoeas. Plant Science, 238, 158–169. Seymour, D. K., & Becker, C. (2017). The causes and consequences of DNA methylome variation in plants. Current Opinion in Plant Biology, 36, 56–63. Skinner, M. K., & Nilsson, E. E. (2021). Role of environmentally induced epigenetic transgenerational inheritance in evolutionary biology: Unified Evolution Theory. Environmental Epigenetics, 7(1), 1–12. Soltani, N., Dille, J. A., Burke, I. C., Everman, W. J., Vangessel, M. J., Davis, V. M., & Sikkema, P. H. (2016). Potential Corn Yield Losses from Weeds in North America. Weed Technology, 30, 979–984. Soltani, N., DIlle, J. A., Burke, I. C., Everman, W. J., Vangessel, M. J., Davis, V. M., & Sikkema, P. H. (2017). Perspectives on Potential Soybean Yield Losses from Weeds in North America. Weed Technology, 31(1), 148–154. Svensson, E. I., & Berger, D. (2019). The Role of Mutation Bias in Adaptive Evolution. Trends in Ecology & Evolution, 34(5), 422–434. Swinton, S. M., & Van Deynze, B. (2017). Hoes to Herbicides: Economics of Evolving Weed Management in the United States. The European Journal of Development Research, 29(3), 560–574. Tehranchian, P., Norsworthy, J. K., Powles, S., Bararpour, M. T., Bagavathiannan, M. v., Barber, T., & Scott, R. C. (2017). Recurrent Sublethal-Dose Selection for Reduced Susceptibility of Palmer Amaranth (Amaranthus palmeri) to Dicamba. Weed Science, 65(2), 206–212. Thom, R., Cummins, I., Dixon, D. P., Edwards, R., Cole, D. J., & Lapthorn, A. J. (2002). 25 Structure of a Tau Class Glutathione S-Transferase from Wheat Active in Herbicide Detoxification. Biochemistry, 41(22), 7008–7020. Van Deynze, B., Swinton, S. M., & Hennessy, D. A. (2022). Are glyphosate‐resistant weeds a threat to conservation agriculture? Evidence from tillage practices in soybeans. American Journal of Agricultural Economics, 104(2), 645–672. Vila-Aiub, M. (2019). Fitness of herbicide-resistant weeds: Current knowledge and implications for management. Plants, 8(11), 469. Vollmer, W. (2008). Structural variation in the glycan strands of bacterial peptidoglycan. FEMS Microbiology Reviews, 32(2), 287–306. Weetman, D., Djogbenou, L. S., & Lucas, E. (2018). Copy number variation (CNV) and insecticide resistance in mosquitoes: evolving knowledge or an evolving problem? Current Opinion in Insect Science, 27, 82–88. Weil, C., & Martienssen, R. (2008). Epigenetic interactions between transposons and genes: lessons from plants. Current Opinion in Genetics & Development, 18(2), 188–192. Wellenreuther, M., Mérot, C., Berdan, E., & Bernatchez, L. (2019). Going beyond SNPs: The role of structural genomic variants in adaptive evolution and species diversification. Molecular Ecology, 28(6), 1203–1209. Yu, P., Wang, C. H., Xu, Q., Feng, Y., Yuan, X. P., Yu, H. Y., Wang, Y. P., Tang, S. X., & Wei, X. H. (2013). Genome-wide copy number variations in Oryza sativa L. BMC Genomics, 14(1), 1–12. Zhang, H., Lang, Z., & Zhu, J. K. (2018). Dynamics and function of DNA methylation in plants. Nature Reviews Molecular Cell Biology 2018 19:8, 19(8), 489–506. Zhu, J. K. (2016). Abiotic Stress Signaling and Responses in Plants. Cell, 167(2), 313–324. Zhu, J., Kapoor, A., Sridhar, V. v., Agius, F., & Zhu, J. K. (2007). The DNA Glycosylase/Lyase ROS1 Functions in Pruning DNA Methylation Patterns in Arabidopsis. Current Biology, 17(1), 54–59. 26 CHAPTER II KOCHIA SCOPARIA TRANSCRIPTOMIC RESPONSE TO GLYPHOSATE Summary Herbicide application is a prime example of human-mediated alteration of the ecosystem in natural, agricultural, and built environments. This is an added abiotic stress that plants must respond and adapt to, most notably by evolving resistance to herbicides. The evolved mechanism of glyphosate resistance in Kochia scoparia is amplification of the EPSPS gene, which is transposon mediated. We hypothesized that glyphosate activates expression of stress-related genes, specifically those involved in xenobiotic detoxification, as well as transposons. To investigate the unique transcriptomic response to glyphosate in kochia, we performed RNA sequencing on leaf tissue from glyphosate-sensitive line ‘7710’ before and three weeks after treatment with two different sublethal doses. We then compared differential expression across time periods and between treatments. Our results show that overall gene expression was suppressed by glyphosate treatment, including downregulation of multiple transcription factors, suggesting that glyphosate impaired the plant’s ability to sense and regulate stress. Gene Ontology annotation of the downregulated gene set following glyphosate application confirmed decreases in nucleic acid binding (GO:0003676), DNA binding (GO:0003677), and protein binding (GO:0005515). We also saw a dose-dependent increase in expression of a FHY3/FAR1 transposon, believed to be responsible for EPSPS amplification in kochia. Our research shows that glyphosate stunts gene expression even at sublethal rates, and the activation of transposons in stress response is supported and should be further studied in non-model species to assess their contributions to stress adaptation. 27 Introduction Kochia (Bassia scoparia (L.) A.J. Scott syn. Kochia scoparia (L.) Schrad.) is an invasive weed that has ravaged grasslands, desert drylands and agricultural fields throughout the U.S. Great Plains and Canada (Friesen et al., 2009). Yield losses reported from kochia can be greater than 90% in sugar beets, corn, sorghum, and sunflower (Geddes and Sharpe, 2022). This C4 tumbleweed can withstand drought and salinity, is a prolific seed producer, and emerges earlier in the spring than other species, giving it a competitive advantage in the climates it invades (Friesen et al., 2009). It has also evolved resistance to multiple herbicides including photosystem II (PSII) inhibitors, acetolactate synthase (ALS) inhibitors, synthetic auxins, and glyphosate (Heap, 2022; Kumar et al., 2019). Glyphosate resistant kochia populations have been reported in 10 states in the U.S. and multiple locations in Canada (Kumar et al., 2019). Glyphosate is the most used herbicide active ingredient in the world and resistance is a costly threat to its continued use (Duke and Powles, 2008). It is a systemic herbicide that enters plant cells through the cuticle via diffusion and is translocated through the phloem, concentrating in “sink” tissue like the meristems (Shaner, 2009). Glyphosate inhibits the shikimate synthesis pathway within the chloroplast by competing with phosphoenolpyruvate (PEP) for the binding site of the enzyme 5-enol-pyruvylshikimate3-phosphate synthase (EPSPS) (Green and Owen, 2010; Shaner, 2009). The shikimate pathway consists of seven enzymatic steps that convert PEP to chorismate, which is needed to form the aromatic amino acids phenylalanine, tyrosine, and tryptophan (Herrmann, 1995; Tzin and Galili, 2010). When this pathway is inhibited by glyphosate, shikimate accumulates in the plant, bursting the chloroplasts, and results in plant death (Shaner, 2009; Singh et al., 2020). Other secondary effects from herbicides include the over-production of reactive oxygen species (ROS), lipid peroxidation, and antioxidant activity 28 which can lead to alterations in photosynthesis, respiration, and overall plant physiology (Eceiza et al., 2022; Gomes and Juneau, 2016). At basal levels, ROS are important signaling molecules essential for regulating plant metabolism and are scavenged by antioxidants such as superoxide dismutase (SOD), peroxidase (POD), glutathione reductase (GR), and glutathione-S-transferase (GST), among others, to maintain homeostasis of the cells (Caverzan et al., 2019). Herbicide non-target-site resistance (NTSR) mechanisms include reduced translocation or uptake of the herbicide, sequestration, reduced activation, or metabolism, which can be mediated by cytochrome P450 monooxygenases (CYP450s) or GSTs (Jugulam and Shyam, 2019). NTSR is mostly thought to be present in standing variation and can evolve in as little as three generations from sublethal herbicide application selecting for multiple minor-effect resistance genes (Busi et al., 2012; Busi & Powles, 2009). Weeds may receive sublethal doses from spray drift, poor application coverage, or when the herbicide is not applied at the recommended development stage (Tehranchian et al., 2017). The mechanism for glyphosate resistance in kochia is target-site resistance through increased copy number of EPSPS, which has also been demonstrated in at least eight other species (Hordeum glaucum - Adu-Yeboah et al., 2019; Eleusine indica - J. Chen et al., 2015; Amaranthus palmeri - Gaines et al., 2010; Amaranthus tuberculatus - Lorentz et al., 2014; Bromus diandrus - Malone et al., 2016; Amaranthus spinosus - Nandula et al., 2014; Chloris truncata - Ngo et al., 2018; Lolium perenne ssp. Multiflorum - Salas et al., 2012; Kochia scoparia - Wiersma et al., 2015). The subsequent increase in expression of EPSPS means that higher rates of glyphosate application are needed to kill the plant (Patterson et al., 2019). In the case of kochia, multiple genes were co-duplicated with EPSPS during the gene amplification event, including a Mutator DNA transposon in the FHY3/FAR1 (FAR-Red Elongated 29 Hypocotyls 3/FAR-Red-Impaired Response 1) transcription factor family (Dupeyron et al., 2019; Patterson et al., 2019; Wang and Wang, 2015). Believed to be “domesticated”, the presence and co-duplication of the FHY3/FAR1 in the kochia EPSPS copy number variant is evidence that this element is still an active MULE transposon (Wang and Wang, 2015). Transposable elements (TEs) are key drivers of evolution and contribute to gene duplication by transposition and by creating homology in the genome for homologous recombination to occur (Hastings et al., 2000; Reams and Roth, 2015). Mutator-like elements, or MULEs, are Class II transposons that exhibit high transposition activity and commonly insert near genes, which is what occurred in kochia (Craig, 2020; Dupeyron et al., 2019; Underwood et al., 2017). TE transcription and transposition is also largely controlled by epigenetic mechanisms such as DNA methylation (Craig, 2020). Previous studies have produced mixed results regarding glyphosate’s effect on gene expression and how quickly the response is taking place, especially in weed species (Josephs et al., 2021; van Etten et al., 2021). In weeds, the primary focus has been on identifying herbicide resistance genes by comparing gene expression in susceptible and resistant biotypes (J. Chen et al., 2017; Duhoux et al., 2015; Gaines et al., 2014). Ipomoea purpurea rapidly responded to glyphosate application by shutting down gene expression eight hours after treatment (Leslie and Baucom, 2014). Additionally, a time-course study of ALS resistant and susceptible Lolium plants displayed considerable variation in gene expression of different gene families between just 0-48 hours after treatment (Duhoux et al., 2015). However, glyphosate acts more slowly in the plant and the most severe visual symptoms are observed more than two weeks after application. This begs the question, how quickly does the transcriptome respond and how long after treatment can this be detected? 30 Our objective is to follow changes in gene expression in glyphosate-susceptible kochia in response to sublethal glyphosate doses three weeks after treatment. We designed this study with the hypothesis that glyphosate acts as an abiotic stress and predicted that signature changes in expression of stress-related genes and transposon activation would be found, with a higher dose resulting in more drastic changes. The plasticity of weedy genomes, the role of transposable elements, and the unique stress-related genes that herbicide application may induce must all be accounted for to explain how weeds sense herbicide stress, how resistance might evolve, and ultimately how it can be prevented and/or delayed. Methods Plant material and sublethal dose determination The previously characterized glyphosate susceptible kochia line ‘7710’ was used in this experiment to ensure homogenous genetic background (Patterson et al., 2019; Pettinga et al., 2018; Preston et al., 2009). Kochia seedlings were grown in a growth chamber at 25°C/23°C (day/night) with a photoperiod of 16 hours, and a light intensity of 600 µmol photon m-2 s-2. To determine appropriate dosage levels, nine glyphosate doses from 0-2x the standard field rate of 840 g ae ha-1 (0, 13, 26, 52.5, 105, 210, 420, 840 and 1,680 g ae ha-1) were assayed with three replicates of 18 seedlings for each dose (Wiersma et al., 2015). Roundup PowerMAX (540 g ae L−1, Bayer Crop Science) with ammonium sulfate (2% w/v) was applied to plants averaging 10 cm in height using the Generation 4 DeVries Research Track Sprayer in the Michigan State University Plant and Pest Diagnostics Pesticide Application Laboratory. The sprayer was calibrated to deliver 187 L ha-1 spray volume at 207 kPa pressure. Plants were then moved to the greenhouse at a temperature of 29.4°C/23.8°C (day/night) with a 16 h photoperiod. All 31 aboveground tissue was collected three weeks after treatment and dried in an Equatherm Oven (Curtin Matheson Scientific, Inc.) at 70°C for 48 hours. Dry weight was the response variable used to generate the dose response curve using a four-parameter log-logistic model in the drc package in R, which was selected using the mselect function (Fig. 2.1) (Ritz et al., 2015). We tested both three- and four-parameter models and chose the one with the smallest information criterion and residual standard error, as well as the largest P value for lack-of-fit test. The median effective dose (ED50) was 92.5 g ae ha-1 and ED25 was 75.7 g ae ha-1. We saw a hormetic effect around the ED25 (Anunciato et al., 2022). The two sublethal doses for this study were chosen as 52.5 and 105 g ae ha-1 (1/16th and 1/8th of a field rate) to ensure we had two “extreme” doses, without impairing the plants’ ability to set seed for future studies (Fig. 2.1). Experimental design and tissue collection Kochia plants from the same line ‘7710’ were grown in individual pots in the greenhouse at 29.4°C/23.8°C with a 16 h photoperiod until plants averaged 10 cm tall. 120 plants per treatment were sprayed as described above with either 0, 52.5 (“low”), or 105 (“high”) g ae ha-1 of glyphosate. Young leaf tissue was collected from four individuals per treatment pre- application (on the date of spraying) and was flash-frozen in liquid nitrogen and stored at -80°C until RNA extraction. After application, plant height was recorded up to two weeks after treatment (Fig. 2.2A). Since visual glyphosate injury symptomology has a long effect, newly expanded leaf tissue three weeks post-application was collected from four individuals per treatment and handled in the same manner as pre-treatment samples. 32 RNA Extraction and sequencing For each plant sampled, total RNA was extracted from leaf tissue using the RNeasy Plant Kit (Qiagen, Germany). RNA concentration, purity and quality were assessed with a NanoDrop One spectrophotometer (Thermo Fisher, USA) and a Qubit 1.0 fluorometer (Invitrogen, USA). RNA sequencing was performed by Novogene Co., Ltd. using Illumina NovaSeq6000 paired end sequencing (PE150) for 24 samples, including four untreated pre (UNT) and post (UNTP), four low dose pre (LOW) and post (LOWP), and four high dose pre (HIGH) and post (HIGHP). Data processing and identification of differentially expressed transcripts Raw FASTQ files of sequence reads were trimmed and quality filtered using fastp with default parameters (S. Chen et al., 2018). Next, a reference transcriptome was created with gffread (Pertea and Pertea, 2020) by extracting transcript sequences from the existing kochia genome annotation file (Patterson et al., 2019). Hisat2 (Kim et al., 2019) was used to align the reads to the reference transcriptome and read counts per gene were obtained with samtools (Danecek et al., 2021). Changes in gene expression resulting from no, low, or high glyphosate treatment were identified by comparison with pretreatment plants, using the package edgeR with quasi-likelihood functionality and a generalized linear model. P values were adjusted to control for the false discovery rate (FDR) by performing multiple testing correction using the Benjamini- Hochberg method (Y. Chen et al., 2016; McCarthy et al., 2012; Robinson et al., 2010). A multidimensional scaling (MDS) plot was generated using the plotMDS function in edgeR to visualize sample variance. Volcano plots were generated with the ggplot2 package in R (Wickham, 2016) to visualize the statistically significant DETs between pre and post timings of each treatment defined as transcripts with both an FDR-adjusted P value < 0.05 and change in expression of absolute Log2Fold Change > 1. 33 Gene Ontology analysis Gene Ontology (GO) terms associated with each DET were retrieved from the kochia genome annotation pipeline, which was used as the “background” set of transcripts. GO term enrichment among DETs associated with response to glyphosate treatment was calculated using a two-proportion z-test, which was confirmed with the AgriGO Single Enrichment Analysis where we saw the same results (Tian et al., 2017). The pheatmap function in R was used to generate a hierarchical clustering heatmap of all DETs and cluster transcripts by their expression pattern across treatments (Kolde, 2019). The number of clusters was chosen by visual interpretation of the response patterns. To test our hypothesis that glyphosate treatment effects gene expression in a dose-dependent manner, all GO terms associated with the transcripts in clusters that followed this pattern were input to Revigo to visualize and interpret their functional categories of Biological Process (BP), Cellular Component (CC), and Molecular Function (MF) based on their annotations (Supek et al., 2011). Analysis of specific gene families (shikimate pathway, stress response, transposons, epigenome) To characterize glyphosate-specific responses in greater detail, we focused on differential expression of pathways with known or strongly suspected roles in glyphosate-mediated injury using the functional annotation of the kochia genome. First, we used the ‘grep’ command to identify genes in the genome annotation file involved in the shikimate pathway (DAHP synthetase, 3-dehydroquinate synthase, 3-dehydroquinate dehydratase, shikimate dehydrogenase, shikimate kinase, EPSPS, chorismate synthase) where glyphosate is acting and compared the significantly regulated genes across treatments. We then found annotated genes involved in abiotic stress response including peroxidases (PODs), superoxide dismutases (SODs), CYP450s, ATP-binding cassette (ABC) transporters, aldo-keto reductases (AKRs), GSTs, and glutathione 34 reductases (GR). These genes also play a role in glyphosate resistance (Jugulam and Shyam, 2019, Gill et al., 2015; Pan et al., 2019). Genes in the FHY3/FAR1 family were also of interest due to their potential to cause glyphosate resistance specifically in kochia. Finally, genes involved in the epigenome through methylation were identified for analysis. Results Glyphosate injury symptomology Glyphosate application had clear effects on growth of susceptible kochia line ‘7710’. There were few noticeable effects after just one week, and mean plant height was not significantly different between any of the treatments. By two weeks after treatment, mean height of plants treated with high sublethal glyphosate was reduced by 9.9% compared to both other doses (P < 0.05) (Fig.2.2A). High dose treated plants displayed chlorotic meristems and stunting (Fig. 2.2B). Three weeks after treatment, even low dose treated plants had chlorotic meristematic tissue leading to some necrosis (Fig. 2.2C). Beyond this time, the treated plants outgrew their injuries and were able to produce seed. RNA sequencing mapping and data visualization Approximately 64 million reads were generated for each of the 24 samples, averaging 9.65 Gb, and representing 45x coverage of the kochia transcriptome. The base error rate for all samples was 0.03%. The average Hisat2 alignment rate to the reference transcriptome was 76%. Summary statistics of sequence data quality are reported in Table 2.1. To visualize the relationships between each of these transcriptomes, an MDS plot was used to graphically represent the pairwise ‘distances’, or dissimilarity measurement. All the pre samples (UNT, LOW, HIGH) cluster together, meaning they have similar gene expression (Fig. 2.3). However, the UNTP and LOWP samples separate from the pre samples along dimension 35 one, and cluster separately from one another. This indicates the development that took place for three weeks between sampling times for UNTP, as well as the treatment effect of glyphosate on LOWP. The HIGHP samples separate along dimension two, but these individuals vary greatly from one another as well. These points are dissimilar from both UNTP and LOWP and could be interpreted to be more like the pre-treatment samples. Analysis of differentially expressed transcripts To identify glyphosate-induced effects independent of developmental changes, we compared differential gene expression between pre and post glyphosate application samples to view changes associated with herbicide treatment. The number of differentially expressed transcripts was similar across treatments, with 24,309 DETs between UNTP and UNT, 23,446 between LOWP and LOW, and 23,740 between HIGHP and HIGH. However, in general, the transcriptome of untreated plants differed more over time than either low or high dosed plants after treatment (Fig. 2.4). The magnitude of transcriptional responses was also greater in untreated, where the range of fold change values of DETs was -15.6 to 14.4, but just -9.2 to 11.1 in low, and -9.3 to 11.4 with high treatment. With significance cutoffs of FDR-adjusted P value < 0.05 and |Log2FC > 1| there were 1,771 upregulated DETs unique to untreated, 170 unique to low, and 92 unique to high treatment (Fig. 2.5A). The two glyphosate treatments shared 13 upregulated transcripts (Fig. 2.5A). Similarly, there were 2,466 downregulated transcripts unique to the untreated plants, 134 for low treatment and 88 for high (Fig. 2.5B). Only 12 transcripts were downregulated in both glyphosate treatments (Fig. 2.5B). Genes shared between untreated and either of the glyphosate treatments are considered as developmental changes and were not analyzed further. 36 Characterizing the transcriptome response specific to high dose treatment The ten genes with the greatest increase in expression unique to high sublethal glyphosate included genes related to RNA processing and degradation, plant hormone signaling, cell wall modification, pathogen defense response, as well as a GST (Table 2.2). Genes with the greatest decrease in expression following treatment have functions in the expansion of the cell wall, telomere regulation, and stress tolerance. Transcription factors involved in stress perception and development as well as DNA binding were also among the top ten downregulated genes that were high-dose-specific (Table 2.2). Gene Ontology analysis of genes with dose dependent expression To characterize functional changes associated with differential gene expression we compared DETs across treatments with a heatmap and then used hierarchical clustering analysis to identify eight clusters of genes with similar expression patterns. Among these clusters, three showed dose-dependent gene expression, with expression decreasing with glyphosate treatment (cluster one, 1,650 DETs) or increasing with glyphosate treatment (clusters four and five, 1,286 and 672 DETs respectively) (Fig. 2.6A). The number of genes with the same GO term annotation was calculated within a given cluster. GO terms associated with the greatest number of genes that were shared between all three clusters are not reported below since their expression could be both increasing and decreasing with glyphosate treatment. The full outputs from Revigo, including all GO IDs, their description, and the number of times each ID is represented in each cluster are listed in Figures 2.8 – 2.16. Genes in cluster one were highly upregulated in UNTP and intermediately upregulated in LOWP, while HIGHP plants had no change in expression (Fig. 2.6B). Based on the number of genes annotated to each GO term, this cluster has high representation of oxidative stress response 37 (GO:0006979), pollen recognition (GO:0048544), and glutathione metabolism (GO:0006749) pathways (Fig. 2.8). Genes in this cluster are generally associated with the cell membrane (GO:0016020, GO:0016021) (Fig. 2.9) and function in heme binding (GO:0020037), DNA- binding transcription factor activity (GO:0003700), and iron ion binding (GO:0005506) (Fig. 2.10). Genes in cluster four had an opposite expression pattern compared to cluster one, with strong downregulation in UNTP and intermediate downregulation in LOWP, with almost no change in HIGHP plants (Fig. 2.6B). BP terms included translation (GO:0006412), photosynthesis (GO:0015979), protein folding and transport (GO:0006457, GO:0015031), and phloem development (GO:0010088) (Fig. 2.11). The top CC terms were ‘membrane’ (GO:0016020) and ‘ribosome’ (GO:0005840) (Fig. 2.12). MF terms included hydrolase activity (GO:0004553), structural constituent of the ribosome (GO:0003735), and zinc and calcium ion binding (GO:0008270, GO:0005509) (Fig. 2.13). Finally, cluster five includes those DETs which were upregulated in the HIGHP plants, downregulated in the LOWP, and further downregulated in UNTP (Fig. 2.6B). BP GO terms in this cluster included genes involved in fatty acid biosynthesis (GO:0006633), DNA repair (GO:0006281), and cell wall modification (GO:0042545) (Fig. 2.14). The CC where these gene products are localized were mostly in the nucleosome (GO:0000786) (Fig. 2.15). MFs included protein kinase activity (GO:0004672), microtubule motor activity (GO:0003777), and protein heterodimerization activity (GO:0046982) (Fig. 2.16). Functional consequences of transcriptome changes due to glyphosate treatment were also assessed by testing for gene ontologies significantly associated with DETs in high dose treatments. We found no significantly enriched GO terms corresponding with all upregulated 38 transcripts in HIGHP (189 DETs) . However, three terms were significantly overrepresented in the downregulated gene set (177 DETs) of HIGHP: nucleic acid binding (GO:0003676, P < 0.01), DNA binding (GO:0003677, P < 0.05), and protein binding (GO:0005515, P < 0.05). Differential expression of stress related and known herbicide resistance gene families One of our hypotheses was that the clear physical stress of sublethal glyphosate doses (slowed or stopped growth, yellowing of meristems) would be reflected in upregulation of stress response genes, including many with known roles in herbicide resistance like ABC transporters, AKRs, CYP450s, GR, GSTs and SOD genes (Jugulam and Shyam, 2019), as well as peroxidases that function in scavenging reactive oxygen species during abiotic stress response, development, and pathogen infection (Bela et al., 2015). Overall, this was not the case as we saw greater changes in expression of these gene families in the untreated plants than those treated with glyphosate. In fact, there was no significant change in expression of ABC transporters, AKRs, GR or SODs after high dose treatment (Table 2.3). There were significant changes in expression of some CYP450s, GST and PODs after high dose treatment. Out of 315 CYP450s, two were upregulated (Bs.00g240410, Bs.00g366530) and two were downregulated (Bs.00g013280, Bs.00g430200) after high dose treatment. One of the upregulated CYPs is a flavonoid 3’-monooxygenase. Additionally, out of 39 GSTs, there was one upregulated (Bs.00g363310) and one downregulated (Bs.00g172490). Out of a total of 135 POD genes, one ascorbate peroxide was upregulated (Bs.00g335350) and two other PODs were significantly downregulated in the high dose treatment (Bs.00g282090, Bs.00g318020). As a proportion of the total annotated genes in these families, those regulated with high sublethal glyphosate application represented 1.5% of CYPs, 5.2% of GSTs and 1.6% of PODs (Table 2.3). 39 We also analyzed 15 genes in the shikimate pathway, including glyphosate’s target, EPSPS, and there were three genes downregulated only in untreated. These genes included an uncharacterized protein suspected to be 3-dehydroquinate synthase, as well as two genes annotated as probable inactive shikimate kinases (Table 2.6). Transposon activation and suppression in response to glyphosate Of the 607 genes in the FHY3/FAR1 family in kochia, 11 were differentially expressed in at least one of the treatments. These genes were mostly changing in expression in the untreated plants, likely functioning in red light sensing. However, one of these is the specific FHY3/FAR1 gene co-duplicated in the EPSPS copy number variant that causes glyphosate resistance in kochia (Gene ID=Bs.00g352450; Hall et al., unpublished), which seems to increase its expression slightly with glyphosate treatment (Fig. 2.7). Whereas it is significantly downregulated in untreated plants over time, low dose and high are not significantly regulating its expression based on the Log2FC values of -1.41, -0.45, and 0.07 in UNTP, LOWP, and HIGHP, respectively. Genes regulating the epigenome The role of the epigenome in regulating gene expression is well documented (Perrone and Martinelli, 2020), therefore we looked at differential expression of 25 genes including DNA methyltransferases, genes in the RNA-directed DNA Methylation pathway, and those involved in demethylation. Of these, there were no differentially expressed genes following high dose glyphosate. However, three methyltransferases responsible for maintaining DNA methylation were downregulated in both untreated and low sublethal glyphosate treated plants, including CMT2 and a MET1 isoform. 40 Discussion The objective of this study was to characterize the transcriptomic response of glyphosate sensitive kochia line ‘7710’ to two different sublethal rates of glyphosate three weeks after treatment. Few studies have characterized the “long-term” effects of herbicide application on gene expression in weedy species. Understanding this can help inform us on how weeds sense herbicide stress and provide evidence for how resistance might evolve in the field. Glyphosate is a slow-acting systemic herbicide and the most severe symptoms can be observed more than two weeks after application, especially at the “sinks” in the plant. We observed significant stunting (10% reduction in plant height) and chlorosis in the meristems of treated plants beginning at two weeks after treatment. This led to necrosis of some axillary buds by three weeks after treatment when we tissue sampled the newly expanded tissue for post- application sequencing. Since it was a sublethal dose, plants were able to recover beyond this time and could be used for future studies on the transgenerational effects of sublethal glyphosate. Glyphosate represses gene expression Our research shows that glyphosate arrested overall gene expression in kochia in a dose- dependent manner three weeks after treatment, beginning at just 1/16th of a field rate (52.5 g ae h- 1 ). Additionally, there was high variability in the transcriptomes of plants treated with the high dose of glyphosate, or 1/8th of a field rate (105 5 g ae h-1) meaning that there was a stochastic response. Untreated plants had 28 times more uniquely significantly upregulated transcripts, and 19 times as many downregulated transcripts compared with the high dose treated plants. The magnitude of the changes in gene expression were also lower with glyphosate treatment. Assuming that changes in untreated are associated with development over this three-week period, this indicates to us that the plants were less able to regulate their expression in response 41 to treatment. However, this seemed to be dose dependent as there were very few shared responses between both the high and low doses. Treatment specific response to glyphosate There were 92 and 88 genes exclusively upregulated and downregulated after high sublethal glyphosate application, respectively. One of the genes with the highest fold change in expression unique to treatment included an RNA exonuclease, functioning in RNA processing and decay. RNA exonucleases are involved in defense response but can contribute to RNA silencing as well (Li and Wang, 2018; Yu et al., 2017), which could be an explanation as to why overall gene expression is reduced with treatment. There were also three genes related to auxin, gibberellin, and jasmonate, confirming further the crosstalk between plant herbicide response and hormone signaling (Alberto et al., 2016). Whether or not the upregulation of these genes is a direct or indirect response should be investigated further. Some of the other highly upregulated genes exclusive to high dose treatment included a chalcone synthase, a putative ripening-related protein 1 (Kiwellin-like), a beta-glucosidase and another antimicrobial protein all involved in defense against pathogens (Dao et al., 2011; Jaswal et al., 2021; Zanatta et al., 2020). Previous research has suggested that sublethal glyphosate makes plants more susceptible to pathogenic bacteria and fungi by suppressing defense response (Sharon et al., 1992). The chalcone synthase is of interest due to its involvement in dicamba resistance in kochia (Pettinga et al., 2018). The genes with the greatest decrease in expression included multiple transcription factors, including a bHLH118-like gene involved in DNA binding and a NAC transcription factor 25- like, which is a key regulator of stress perception and development (Jensen et al., 2010; Toledo- Ortiz et al., 2003). This is yet another explanation for the overall decrease in gene expression with treatment, as the plants’ ability to regulate expression under stress conditions was impaired. 42 Further evidence of this would be the downregulation of a PIF1-like ATP-dependent DNA helicase, which functions in telomere regulation. This DNA helicase prevents telomere formation at DNA double strand breaks (Seraj et al., 2018), which also relates to our GO annotation that revealed DETs involved in DNA repair (GO:0006281) were upregulated under high sublethal doses of glyphosate. Another downregulated gene was a 21 kDa protein-like, involved in stress tolerance to salt and drought, indicating glyphosate induces shared abiotic stress responses (Aarati et al., 2003; Ghodke et al., 2020). Analysis of glyphosate induced stress-related genes We found little evidence of upregulation of various gene families involved in stress response including PODs, SODs, CYP450s, ABC transporters, AKRs, GSTs, and GR, though the oxidative stress effect of glyphosate has been reported in several other plant species (Ahsan et al., 2008; Gomes and Juneau, 2016; Liu et al., 2019; Moldes et al., 2008; Radwan and Fayez, 2016; Sergiev et al., 2006; Soares et al., 2019). High sublethal glyphosate treatment only regulated four CYP450s, one of which is a flavonoid 3’-monooxygenase, two GSTs, and three PODs, including one ascorbate peroxide. These represented a smaller proportion of the total annotated genes in these families compared with both untreated and low dose treated plants. One explanation for this observation is that previous research found that GST expression peaks 48 hours after glyphosate treatment, though RNA was extracted from all aboveground tissue in this case (J. Chen et al., 2016). Tissue and sampling time are key factors in the response we have recorded. Perhaps there were greater responses prior to the three-week mark. Therefore, we reject our hypothesis that glyphosate would activate the expression of these stress related gene families three weeks after treatment, and instead suggest the need for time course studies to determine the peak sampling time for these families to capture the full plant response. 43 The role of transposons in glyphosate stress response Interestingly, we did find that a FHY3/FAR1 gene named Muntjac, believed to be involved in the gene amplification event causing glyphosate resistance in kochia, was constitutively downregulated in untreated plants, but slightly upregulated with glyphosate treatment. As previously mentioned, the epigenome modulates gene expression, and we found no effect of glyphosate on genes related to methylation. So, this could indicate that in untreated plants, Muntjac is being suppressed through methylation, while glyphosate treatment is releasing it to be expressed. Perhaps at even higher doses, we would see a greater increase in expression, possibly leading to the transposition event creating copy number variation of EPSPS for glyphosate resistance in this susceptible line. Future research will continue this study for multiple generations to investigate this hypothesis and potentially witness the transposon mediated evolution of glyphosate resistance in kochia. There are a multitude of other mechanisms involved in stress response modulation, including long-non-coding RNA, micro-RNA (miRNA), and small-interfering RNA (siRNA), as well as other post-transcriptional or post-translational modifications (Giacomini et al., 2018; Jha et al., 2020) shown to be involved in salt, drought, temperature, heavy metal, and oxidative stress (Banerjee, 2020). These could be responsible for some of the differences we are seeing but were not captured with our approach. Conclusion Kochia is a problematic weed and widespread glyphosate resistance is making control efforts difficult for producers and land-managers alike. An understanding of weed biology and this rapid evolution occurring in the field will help identify ways to address this problem that has both economic and environmental costs. Interestingly, the mechanism of resistance in kochia 44 through gene amplification of glyphosate’s target gene, EPSPS, has evolved in eight other species, leading researchers to question whether glyphosate is activating unique herbicide responses. We performed RNA sequencing of glyphosate sensitive kochia to determine the transcriptomic response to sublethal doses of glyphosate three weeks after treatment. Our findings suggest that gene expression is largely suppressed by glyphosate treatment, even at 1/16th of a field rate. We also observed minimal changes in key abiotic stress related genes, again supporting the notion that sublethal doses repress the weed’s ability to perceive and regulate stress. We also found that a gene in the FHY3/FAR1 family of supposed domesticated MULE transposons, thought to be involved in the evolution of glyphosate resistance in kochia, increases expression in a dose-dependent manner under sublethal glyphosate selection. This strengthens the argument for increased studies on the role of transposons in stress adaptation, especially relating to glyphosate resistance due to copy number variation. Overall, this study provides us with valuable insight on how an economically important weed species senses and responds to sublethal herbicides application and helps generate hypotheses on how resistance may evolve in the field. 45 APPENDIX 46 Figure 2.1: Glyphosate dose response curve. Glyphosate dose response of known glyphosate-sensitive Kochia scoparia line ‘7710’ generated by measuring dry weight of treated plants (three replicates, 18 seedlings each) 21 days after treatment with nine glyphosate doses from 0-2x the standard field rate of 840 g ae ha-1 (0, 13, 26, 52.5, 105, 210, 420, 840 and 1,680 g ae ha-1). Plants were grown in a growth chamber and transferred to the greenhouse once sprayed. Solid red line is the ED50 which occurs at 92.5 g ae ha-1 of glyphosate and the red dotted line shows the ED25 at 75.7 g ae ha-1. Gray area shows 95% confidence interval. Orange triangle indicates high sublethal dose (105 g ae ha-1) and blue triangle indicates low dose (52.5 g ae ha-1) used in the rest of the experiment. 47 Sample Total raw reads Raw data (G) Effective (%) Q20 (%) Q30 (%) GC (%) Hisat2 alignment rate (%) HIGH1 65,341,578 9.8 98.48 97.53 93.25 45.36 76.04 HIGH2 64,733,290 9.7 98.36 97.75 93.72 45.53 76.76 HIGH3 61,983,396 9.3 98.34 97.57 93.34 45.22 76.93 HIGH4 64,715,900 9.7 98.00 97.58 93.37 45.36 77.03 HIGHP1 68,656,762 10.3 98.29 97.59 93.43 44.80 75.26 HIGHP2 64,998,950 9.7 97.75 97.57 93.37 45.43 76.06 HIGHP3 70,008,316 10.5 97.92 97.61 93.50 45.58 68.79 HIGHP4 60,146,362 9.0 98.03 97.69 93.59 45.05 75.77 LOW1 59,002,842 8.9 98.31 97.40 92.94 45.15 76.89 LOW2 67,804,816 10.2 98.44 97.52 93.21 45.61 76.80 LOW3 63,111,206 9.5 98.47 97.56 93.32 45.62 76.67 LOW4 60,807,182 9.1 98.24 97.70 93.61 45.53 76.73 LOWP1 61,695,418 9.3 98.42 97.62 93.44 45.25 76.38 LOWP2 65,582,452 9.8 98.17 97.61 93.41 45.26 76.29 LOWP3 66,419,036 10.0 98.03 97.51 93.18 44.86 76.28 LOWP4 70,566,768 10.6 98.35 97.45 93.05 45.09 76.83 UNT1 61,360,226 9.2 98.13 97.69 93.28 45.47 76.36 UNT2 62,489,060 9.4 98.35 97.70 93.32 45.35 76.45 UNT3 68,196,692 10.2 98.31 96.53 91.11 45.63 75.78 UNT4 59,049,900 8.9 98.50 97.69 93.60 45.47 76.63 UNTP1 60,058,460 9.0 98.44 97.59 93.34 44.00 76.04 UNTP2 63,292,456 9.5 98.13 97.36 92.85 44.18 75.75 UNTP3 64,483,540 9.7 98.38 97.54 93.22 43.99 75.98 UNTP4 68,593,614 10.3 98.21 97.55 93.27 44.10 75.20 Table 2.1: RNA-Seq sample summary statistics. Summary of kochia RNA sequencing data adapted from Novogene. RNA was extracted from young leaf tissue from four replicates each in UNT = Untreated Pre, UNTP = Untreated Post, LOW = Low dose pre, LOWP = Low Dose Post, HIGH = High Dose Pre, and HIGHP = High Dose Post treatment. Pre samples were from 10 cm tall plants, which were then sprayed with 0, 52.5, or 105 g ae h-1 of glyphosate. Post samples were collected three weeks after treatment. Effective is the proportion of clean reads out of total reads, and Q20 and Q30 are the proportion of bases with Phred values greater than 20 or 30. GC content is listed as a percentage of total bases. Hisat2 alignment rates to the kochia reference transcriptome are also reported. 48 Figure 2.2: Plant height and phenotypic response. A) Boxplot of kochia height measurements recorded at time of spraying and up to two weeks after treatment (WAT) with 0, 52.5 and 105 g ae ha-1 of glyphosate. Plants were grown in a greenhouse setting. Median is represented by the notches and error bars indicate mean ± standard error. Statistical significance between treatments is shown by a student’s two-tailed T-test with heteroscedastic distribution. Boxes within each week after treatment followed by the same letter are not statistically different at P < 0.05. B) Comparison at 2 WAT of a representative untreated plant (left) and high sublethal dose treated plant (right); the treated plant has chlorotic meristems indicated with white arrows and stunted growth depicted by difference in height showed by the orange tray. C) Consequences of low dose 3 WAT when tissue samples were taken. Severe yellowing is apparent at the growing points, particularly the apical meristem indicated by the white arrow. None of the plants were flowering yet at the time of sampling 3 WAT. 49 Figure 2.3: Gene expression MDS plot. MDS plot used to visualize comparisons between 24 sequenced kochia transcriptomes generated from leaf tissue before and three weeks after treatment with 0, 52.5 and 105 g ae ha-1 of glyphosate. Plants were grown in a greenhouse setting. Each point represents the transcriptome of an individual and the distance between points represents the difference in their transcriptomes where UNT = Untreated Pre, UNTP = Untreated Post, LOW = Low dose pre, LOWP = Low Dose Post, HIGH = High Dose Pre, and HIGHP = High Dose Post treatment. Since UNT and UNTP samples separate along MDS1, this would represent changes in development over the three-week time period. However, there is also clear separation between post samples indicating there is a treatment effect from glyphosate and it is dose dependent since LOWP and HIGHP samples are dissimilar. 50 Figure 2.4: Volcano plot of DETs by treatment. Volcano plot of all differentially expressed transcripts between kochia transcriptomes generated from leaf tissue before and three weeks after treatment with 0 (Untreated), 52.5 (Low) and 105 (High) g ae ha-1 of glyphosate. Plants were grown in a greenhouse setting and included four replicates from each of the treatments both pre and post application timings. Each point represents a transcript, and significant (FDR-adjusted P value < 0.05 and |Log2FC > 1|) differential expression is represented by point color, where red indicates upregulated transcripts, blue indicates downregulated transcripts, and gray indicates no significant regulation. Note that axes are scaled the same across plots. A) 24,309 total DETs between UNTP and UNT, B) 23,446 total DETs between LOWP and LOW, C) 23,740 total DETs between HIGHP and HIGH, where UNT = Untreated Pre, UNTP = Untreated Post, LOW = Low dose pre, LOWP = Low Dose Post, HIGH = High Dose Pre, and HIGHP = High Dose Post treatment. 51 Figure 2.5: Distribution of significant DETs between treatments. Venn diagrams of significantly A) upregulated or B) downregulated differentially expressed transcripts (DETs) between kochia transcriptomes generated from leaf tissue before and three weeks after treatment with 0 (Untreated), 52.5 (Low) and 105 (High) g ae ha-1 of glyphosate. Plants were grown in a greenhouse setting. Overlaps show a shared treatment response. Significant differential expression was defined as an adjusted P value < 0.05 and |Log2FC > 1|. A total of 2,423 DETs were significantly upregulated across all treatments and 3,116 were significantly downregulated. 52 GeneID Annotation High Pvalue Sig. Low Pvalue Sig. Unt Pvalue Sig. Bs.00g286820 PREDICTED: RNA exonuclease 4 5.82 0.03 UP 2.93 0.06 NS 2.29 0.07 NS Bs.00g477580 PREDICTED: probable indole-3-acetic 4.19 0.02 UP -0.31 0.84 NS 0.46 0.77 NS acid-amido synthetase GH3.1 Bs.00g130650 putative ripening-related protein 1 4.16 0.04 UP 0.00 0.00 NS 0.48 0.66 NS (Kiwellin-like) Bs.00g518030 chalcone synthase (Polyketide synthase, 3.93 0.03 UP 1.55 0.15 NS 1.28 0.15 NS type III) Bs.00g122680 23 kDa jasmonate-induced protein-like 3.87 0.03 UP 0.87 0.55 NS 0.04 0.98 NS Bs.00g310490 PREDICTED: gibberellin-regulated 3.73 0.03 UP -0.12 0.96 NS 0.00 0.00 NS protein 1 Bs.00g392340 hypothetical protein BevumaM_p137 3.06 0.03 UP -0.43 0.69 NS 1.20 0.15 NS (mitochondrion) (NADH:quinone oxidoreductase/Mrp antiporter, membrane subunit) Bs.00g180610 beta-glucosidase 11-like isoform X2 3.00 0.03 UP -0.72 0.47 NS 0.16 0.83 NS Bs.00g231030 PREDICTED: protein 2.85 0.02 UP 0.44 0.50 NS 0.14 0.63 NS DETOXIFICATION 35 (Multi antimicrobial extrusion protein) Bs.00g363310 Glutathione-S-transferase U9-like 2.75 0.05 UP 1.00 0.29 NS 1.28 0.06 NS Bs.00g007570 NAC transcription factor 25-like -9.33 0.04 DOWN -1.84 0.15 NS -0.27 0.82 NS Bs.00g393910 putative expansin-B2 -9.02 0.02 DOWN 0.00 0.00 NS 0.00 0.00 NS Bs.00g297780 putative disease resistance protein -7.41 0.04 DOWN 0.00 0.00 NS 0.00 0.00 NS At1g50180 (P-loop containing nucleoside triphosphate hydrolase) Bs.00g228390 transcription factor bHLH118-like -5.95 0.03 DOWN -2.81 0.09 NS 0.00 0.00 NS Bs.00g012430 transducin family protein / WD-40 -5.11 0.05 DOWN -1.27 0.54 NS 0.00 0.00 NS repeat family protein (Major facilitator, sugar transporter-like) Table 2.2: High sublethal glyphosate responsive genes. Top ten differentially expressed genes either upregulated (UP) or downregulated (DOWN) exclusively in kochia leaf tissue three weeks after treatment with 105 (High) g ae ha-1 of glyphosate. The associated expression values in plants treated with 52.5 (Low) or 0 (Unt) g ae h-1 glyphosate are also shown. For each gene, its 53 Table 2.2 (cont’d) GeneID Annotation High Pvalue Sig. Low Pvalue Sig. Unt Pvalue Sig. Bs.00g122100 PREDICTED: uncharacterized protein -4.89 0.02 DOWN -2.28 0.05 NS -1.30 0.14 NS LOC104907727 Bs.00g523690 PREDICTED: ATP-dependent DNA -4.38 0.04 DOWN -1.67 0.21 NS 0.00 0.00 NS helicase PIF1-like Bs.00g492890 PREDICTED: probable mannitol -3.77 0.02 DOWN -0.71 0.20 NS -0.73 0.12 NS dehydrogenase (NAD(P)-binding domain superfamily) Bs.00g439720 PREDICTED: serine carboxypeptidase- -3.49 0.03 DOWN -1.00 0.32 NS 0.63 0.25 NS like 28 Bs.00g534680 probable xyloglucan -3.24 0.02 DOWN -1.84 0.09 NS -1.10 0.48 NS endotransglucosylase/hydrolase protein 33 annotation, log2Fold Change values by treatment, their associated FDR-adjusted P values, and significance (Sig.) are shown. NS = not significant. 54 Figure 2.6: DET heatmap and cluster summary. A) Hierarchical clustering heatmap of all significantly differentially expressed transcripts (FDR-adjusted P value < 0.05 and |Log2FC > 1|) between kochia leaf tissue pre-treatment and three weeks after treatment with 0 (Untreated), 52.5 (Low) and 105 (High) g ae ha-1 of glyphosate. Plants were grown in a greenhouse setting and included four replicates from each of the treatments both pre and post application timings. Clusters were defined based on expression patterns across treatments, using Z-score normalized Log2Fold Change values. The number of genes in each cluster is listed parenthetically. B) Trends of the average Log2 Fold Change of the DETs represented in each heatmap cluster with increased glyphosate treatment, where the transcripts in cluster 1 would be decreasing with an increased dose, and cluster 5 would be increasing their expression with increased glyphosate. 55 Gene Family ABC AKR CYP450 GR GST POD SOD Total Genes Annotated 252 22 315 3 39 135 10 Untreated – UP | DOWN 32 | 19 7|0 49 | 15 0|2 9|4 14 | 20 1|1 Proportion of total 12.7% | 7.5% 31.8% | 0% 15.6% | 4.8% 0% | 66.7% 23% | 10.3% 10.4% | 14.8% 10% | 10% Low – UP | DOWN 3|8 1|0 14 | 7 0|1 0|4 3|3 0|0 Proportion of total 1.2% | 3.2% 4.5% | 0% 4.4% | 2.2% 0% | 33.3% 0% | 10.3% 2.2% | 2.2% 0% | 0% High – UP | DOWN 0|0 0|0 2|2 0 | 0 1|1 1|2 0|0 Proportion of total 0% | 0% 0% | 0% 0.6% | 0.6% 0% | 0% 2.6% | 2.6% 0.07% | 1.5% 0% | 0% Table 2.3: Stress related gene families. Number and proportion of the total genes in each stress related gene family significantly upregulated (UP) and downregulated (DOWN) in kochia leaf tissue three weeks after treatment with 0 (Untreated), 52.5 (Low) or 105 (High) g ae ha-1 of glyphosate. Plants were grown in a greenhouse setting. The total number of genes annotated in the kochia genome in each family is also reported. ABC = ABC transporters, AKR = aldo-keto reductases, CYP450 = cytochrome P450 monooxygenase, GR = glutathione reductase, GST = glutathione-S-transferase, POD = peroxidase, SOD = superoxide dismutase. 56 Figure 2.7: FHY3/FAR1 related genes heatmap. Heatmap of 11 FHY3/FAR1 related genes significantly differentially expressed between kochia leaf tissue pre-treatment and three weeks after treatment with 0 (Untreated), 52.5 (Low) and 105 (High) g ae ha-1 of glyphosate. Plants were grown in a greenhouse setting and included four replicates from each of the treatments both pre and post application timings. A black star indicates the specific FHY3/FAR1 involved in the gene amplification event causing glyphosate resistance in kochia (gene ID Bs.00g352450). 57 Gene ID Annotation Notes Low FDR- High FDR- Log2FC adjusted Log2FC adjusted P value P value Bs.00g237980 bark storage protein A-like 11.12 0.001 9.06 0.046 Bs.00g268300 PREDICTED: early light-induced protein 2, chloroplastic 4.56 0.005 4.93 0.036 isoform X2 Bs.00g253650 60S acidic ribosomal protein P0-like 4.95 0.001 4.82 0.015 Bs.00g055130 oxysterol-binding protein-related protein 4B-like isoform X1 3.24 0.017 4.08 0.022 Bs.00g484920 Unknown Protein 3.32 0.022 2.53 0.030 Bs.00g445390 non-specific lipid-transfer protein 2-like 1.86 0.006 2.34 0.028 Bs.00g480640 PREDICTED: cytosolic sulfotransferase 12 2.95 0.001 2.27 0.022 Bs.00g155130 vacuolar iron transporter homolog 1-like 2.02 0.046 2.13 0.031 Bs.00g301430 PREDICTED: probable sulfate transporter 3.4 1.40 0.003 1.77 0.015 Bs.00g410990 PREDICTED: AAA-ATPase At3g28580-like 1.47 0.005 1.77 0.019 Bs.00g537920 probable carboxylesterase 5 1.53 0.026 1.58 0.030 Bs.00g448670 heat shock cognate 70 kDa protein 2-like 1.14 0.026 1.55 0.041 Bs.00g309340 MADS-box protein AGL24-like 1.18 0.010 1.19 0.034 Table 2.4: Upregulated genes shared between low and high dose glyphosate. The 13 differentially expressed genes that were upregulated in kochia leaf tissue three weeks after treatment with with 52.5 (Low) and 105 (High) g ae h-1 glyphosate but did not change in expression in the untreated plants. Plants were grown in a greenhouse setting. Gene ID, annotation in the kochia genome, and the Log2Fold Change and adjusted P values for both Low and High treatment are reported. 58 Gene ID Annotation Notes Low FDR- High FDR- Log2FC adjusted Log2FC adjusted P value P value Bs.00g193410 xyloglucan endotransglucosylase/hydrolase protein 24-like -2.73 0.002 -2.19 0.036 Bs.00g429720 PREDICTED: LOW QUALITY PROTEIN: systemin receptor -1.03 0.009 -1.77 0.038 SR160-like Bs.00g114750 uncharacterized protein LOC110708747 isoform X2 -1.35 0.003 -1.74 0.015 Bs.00g161180 tryptophan aminotransferase-related protein 4-like -1.30 0.015 -1.65 0.029 Bs.00g318020 peroxidase P7-like -1.86 0.003 -1.63 0.025 Bs.00g034600 PREDICTED: uncharacterized protein LOC104901445 -1.40 0.033 -1.59 0.037 Bs.00g129800 probable serine/threonine-protein kinase At1g01540 -1.21 0.035 -1.58 0.044 Bs.00g070260 PREDICTED: disease resistance protein RPM1 isoform X1 -2.11 0.009 -1.39 0.027 Bs.00g070550 PREDICTED: probable LRR receptor-like serine/threonine- -1.09 0.031 -1.10 0.038 protein kinase At4g26540 Bs.00g462550 uncharacterized protein LOC110735216 -1.04 0.004 -1.08 0.046 Bs.00g035680 PREDICTED: proline synthase co-transcribed bacterial homolog -1.24 0.003 -1.07 0.041 protein isoform X1 Bs.00g352280 PREDICTED: serine carboxypeptidase-like 34 -1.28 0.013 -1.03 0.039 Table 2.5: Downregulated genes shared between low and high dose treatment. The 12 differentially expressed genes that were downregulated in kochia leaf tissue three weeks after treatment with with 52.5 (Low) and 105 (High) g ae h-1 glyphosate but did not change in expression in the untreated plants. Plants were grown in a greenhouse setting. Gene ID, annotation in the kochia genome, and the Log2Fold Change and adjusted P values for both Low and High treatment are reported. 59 Figure 2.8: Biological Process terms from heatmap cluster one. Gene Ontology terms in the Biological Process category represented by 1,650 genes that are decreasing expression in kochia leaf tissue three weeks after treatment with increasing doses of glyphosate (0 (untreated), 52.5 (low) and 105 (high) g ae h-1) summarized using Revigo. Plants in this experiment were grown in a greenhouse setting. The cluster was defined using a hierarchical clustering heatmap, whose genes are highly upregulated in the untreated, intermediately upregulated with low sublethal glyphosate treatment, and high dose treated plants have no change in expression. The X-axis indicates the number of times a particular GO term was represented in the cluster (number of genes with this GO annotation), while the Y-axis shows the term description. 60 Figure 2.9: Cellular Component terms from heatmap cluster one. Gene Ontology terms in the Cellular Component category represented by 1,650 genes that are decreasing expression in kochia leaf tissue three weeks after treatment with increasing doses of glyphosate (0 (untreated), 52.5 (low) and 105 (high) g ae h-1) summarized using Revigo. Plants in this experiment were grown in a greenhouse setting. The cluster was defined using a hierarchical clustering heatmap, whose genes are highly upregulated in the untreated, intermediately upregulated with low sublethal glyphosate treatment, and high dose treated plants have no change in expression. The X-axis indicates the number of times a particular GO term was represented in the cluster (number of genes with this GO annotation), while the Y-axis shows the term description. 61 Figure 2.10: Molecular Function terms from heatmap cluster one. Gene Ontology terms in the Molecular Function category represented by 1,650 genes that are decreasing expression in kochia leaf tissue three weeks after treatment with increasing doses of glyphosate (0 (untreated), 52.5 (low) and 105 (high) g ae h-1) summarized using Revigo. Plants in this experiment were grown in a greenhouse setting. The cluster was defined using a hierarchical clustering heatmap, whose genes are highly upregulated in the untreated, intermediately upregulated with low sublethal glyphosate treatment, and high dose treated plants have no change in expression. The X-axis indicates the number of times a particular GO term was represented in the cluster (number of genes with this GO annotation), while the Y-axis shows the term description. 62 Figure 2.11: Biological Process terms from heatmap cluster four. Gene Ontology terms in the Biological Process category represented by 1,286 genes that are increasing expression in kochia leaf tissue three weeks after treatment with increasing doses of glyphosate (0 (untreated), 52.5 (low) and 105 (high) g ae h-1) summarized using Revigo. Plants in this experiment were grown in a greenhouse setting. The cluster was defined using a hierarchical clustering heatmap and is grouped by differentially expressed transcripts that are downregulated in untreated, intermediately downregulated in low dose treated plants, and show little change in expression in high dose treatment. The X-axis indicates the number of times a particular GO term was represented in the cluster (number of genes with this GO annotation), while the Y-axis shows the term description. 63 Figure 2.12: Cellular Component terms from heatmap cluster four. G Gene Ontology terms in the Cellular Component category represented by 1,286 genes that are increasing expression in kochia leaf tissue three weeks after treatment with increasing doses of glyphosate (0 (untreated), 52.5 (low) and 105 (high) g ae h-1) summarized using Revigo. Plants in this experiment were grown in a greenhouse setting. The cluster was defined using a hierarchical clustering heatmap and is grouped by differentially expressed transcripts that are downregulated in untreated, intermediately downregulated in low dose treated plants, and show little change in expression in high dose treatment. The X-axis indicates the number of times a particular GO term was represented in the cluster (number of genes with this GO annotation), while the Y-axis shows the term description. 64 Figure 2.13: Molecular Function terms from heatmap cluster four. Gene Ontology terms in the Molecular Function category represented by 1,286 genes that are increasing expression in kochia leaf tissue three weeks after treatment with increasing doses of glyphosate (0 (untreated), 52.5 (low) and 105 (high) g ae h-1) summarized using Revigo. Plants in this experiment were grown in a greenhouse setting. The cluster was defined using a hierarchical clustering heatmap and is grouped by differentially expressed transcripts that are downregulated in untreated, intermediately downregulated in low dose treated plants, and show little change in expression in high dose treatment. The X-axis indicates the number of times a particular GO term was represented in the cluster (number of genes with this GO annotation), while the Y-axis shows the term description. 65 Figure 2.14: Biological Process terms from heatmap cluster five. Gene Ontology terms in the Biological Process category represented by 672 genes that are increasing expression in kochia leaf tissue three weeks after treatment with increasing doses of glyphosate (0 (untreated), 52.5 (low) and 105 (high) g ae h-1) summarized using Revigo. Plants in this experiment were grown in a greenhouse setting. The cluster was defined using a hierarchical clustering heatmap and is grouped by differentially expressed transcripts that are upregulated under high dose treatment, and downregulated in low dose and untreated, to a greater degree. X-axis indicates the number of times a particular GO term was represented in the cluster (number of genes with this GO annotation), while the Y-axis shows the term description. 66 Figure 2.15: Cellular Component terms from heatmap cluster five. Gene Ontology terms in the Cellular Component category represented by 672 genes that are increasing expression in kochia leaf tissue three weeks after treatment with increasing doses of glyphosate (0 (untreated), 52.5 (low) and 105 (high) g ae h-1) summarized using Revigo. Plants in this experiment were grown in a greenhouse setting. The cluster was defined using a hierarchical clustering heatmap and is grouped by differentially expressed transcripts that are upregulated under high dose treatment, and downregulated in low dose and untreated, to a greater degree. X-axis indicates the number of times a particular GO term was represented in the cluster (number of genes with this GO annotation), while the Y-axis shows the term description. 67 Figure 2.16: Molecular Function terms from heatmap cluster five. Gene Ontology terms in the Molecular Function category represented by 672 genes that are increasing expression in kochia leaf tissue three weeks after treatment with increasing doses of glyphosate (0 (untreated), 52.5 (low) and 105 (high) g ae h-1) summarized using Revigo. Plants in this experiment were grown in a greenhouse setting. The cluster was defined using a hierarchical clustering heatmap and is grouped by differentially expressed transcripts that are upregulated under high dose treatment, and downregulated in low dose and untreated, to a greater degree. X-axis indicates the number of times a particular GO term was represented in the cluster (number of genes with this GO annotation), while the Y-axis shows the term description. 68 Step GeneID Annotation Notes Unt Sig. Low Sig. High Si LogFC LogFC LogFC g. 1 Bs.00g253620 casein kinase 1-like protein HD16 isoform X1 (DAHP 0.02 NS 0.46 NS 0.63 NS synthetase, class II) 1 Bs.00g142530 PREDICTED: phospho-2-dehydro-3-deoxyheptonate -0.81 NS -0.33 NS 0.88 NS aldolase 1, chloroplastic (DAHP synthetase, class II) 1 Bs.00g311790 phospho-2-dehydro-3-deoxyheptonate aldolase 2, 0.55 NS 0.12 NS -0.02 NS chloroplastic-like (DAHP synthetase, class II) 1 Bs.00g459970 uncharacterized protein LOC110706710 (DAHP -0.06 NS 0.06 NS -0.74 NS synthetase, class II) 2 Bs.00g461240 PREDICTED: uncharacterized protein -2.09 DOWN 0.09 NS -0.19 NS LOC104906292 isoform X3 (3-dehydroquinate synthase) 2 Bs.00g484840 3-dehydroquinate synthase, chloroplastic-like (3- -0.22 NS 0.05 NS -0.05 NS dehydroquinate synthase AroB) 3,4 Bs.00g071560 bifunctional 3-dehydroquinate dehydratase/shikimate -0.07 NS -0.34 NS -1.60 NS dehydrogenase, chloroplastic-like 3,4 Bs.00g269590 bifunctional 3-dehydroquinate dehydratase/shikimate -0.64 NS -0.82 NS 0.07 NS dehydrogenase, chloroplastic 3,4 Bs.00g436510 bifunctional 3-dehydroquinate dehydratase/shikimate 0.14 NS 0.01 NS 0.13 NS dehydrogenase, chloroplastic 5 Bs.00g261010 probable inactive shikimate kinase like 2, -1.18 DOWN -0.20 NS -0.25 NS chloroplastic isoform X1 5 Bs.00g484150 probable inactive shikimate kinase like 1, -1.02 DOWN 0.33 NS 0.03 NS chloroplastic isoform X1 5 Bs.00g485500 shikimate kinase, chloroplastic -0.70 NS 0.34 NS 0.25 NS 6 Bs.00g132470 EPSPS; 3-phosphoshikimate 1- -0.16 NS -0.32 NS 0.29 NS carboxyvinyltransferase 2-like 7 Bs.00g449990 chorismate synthase 1, chloroplastic 0.16 NS -0.43 NS -0.15 NS 7 Bs.00g129760 PREDICTED: uncharacterized protein 0.00 NS -0.02 NS -0.61 NS LOC104887015 (CHORISMATE SYNTHASE) Table 2.6: Shikimate pathway. Expression of genes in the shikimate pathway in kochia leaf tissue three weeks after treatment with with 0 (Unt), 52.5 (Low) and 105 (High) g ae h-1 glyphosate. Plants were grown in a greenhouse setting. Each gene’s step in the 69 Table 2.6 (cont’d): seven step enzymatic shikimate pathway, gene ID, annotation in the kochia genome, and the Log2Fold Change and adjusted P values for each treatment are reported. Significance (Sig.) is indicated as upregulated (UP), downregulated (DOWN), or not significant (NS) based on the cutoffs |Log2Fold Change > 1| and adjusted P values < 0.05. 70 REFERENCES 71 REFERENCES Aarati, P., Krishnaprasad, B. T., Ganeshkumar, Savitha, M., Gopalakrishna, R., Ramamohan, G., and Udayakumar, M. (2003). Expression of an ABA responsive 21 kDa protein in finger millet (Eleusine coracana Gaertn.) under stress and its relevance in stress tolerance. Plant Science, 164(1), 25–34. Adu-Yeboah, P., Malone, J. M., Fleet, B., Gill, G., and Preston, C. (2019). EPSPS gene amplification confers resistance to glyphosate resistant populations of Hordeum glaucum Stued (northern barley grass) in South Australia. Pest Management Science, 76(4), 1214– 1221. Ahsan, N., Lee, D.-G., Lee, K.-W., Alam, I., Lee, S.-H., Bahk, J. D., and Lee, B.-H. (2008). Glyphosate-induced oxidative stress in rice leaves revealed by proteomic approach. Plant Physiology and Biochemistry, 46(12), 1062–1070. Alberto, D., Serra, A.-A., Sulmon, C., Gouesbet, G., & Couée, I. (2016). Herbicide-related signaling in plants reveals novel insights for herbicide use strategies, environmental risk assessment and global change assessment challenges. Science of The Total Environment, 569–570, 1618–1628. Anunciato, V. M., Bianchi, L., Gomes, G. L. G. C., Velini, E. D., Duke, S. O., and Carbonari, C. A. (2022). Effect of low glyphosate doses on flowering and seed germination of glyphosate-resistant and -susceptible Digitaria insularis. Pest Management Science, 78(3), 1227–1239. Aseel, D. G., Elkobrosy, D. H., Abdelsalam, N. R., El-Saedy, M. A., Shama, S., and Hafez, E. E. (2020). Effect of cyst nematode (Globodera rostochiensis) isolate ddh1 on gene expression in systemic leaves of potato plant. Journal of Microbiology, Biotechnology and Food Sciences, 10(1), 93–97. Banerjee, P. (2020). Plant Abiotic Stress Responses and MicroRNAs. Advanced Agriculture. New Delhi Publishers. Bela, K., Horváth, E., Gallé, Á., Szabados, L., Tari, I., and Csiszár, J. (2015). Plant glutathione peroxidases: Emerging role of the antioxidant enzymes in plant development and stress responses. Journal of Plant Physiology, 176, 192–201. Bewick, A. J., Niederhuth, C. E., Ji, L., Rohr, N. A., Griffin, P. T., Leebens-Mack, J., and Schmitz, R. J. (2017). The evolution of CHROMOMETHYLASES and gene body DNA methylation in plants. Genome Biology, 18(1), 65. 72 Busi, R., Gaines, T. A., Walsh, M. J., & Powles, S. B. (2012). Understanding the potential for resistance evolution to the new herbicide pyroxasulfone: Field selection at high doses versus recurrent selection at low doses. Weed Research, 52(6), 489–499. Busi, R., & Powles, S. B. (2009). Evolution of glyphosate resistance in a Lolium rigidum population by glyphosate selection at sublethal doses. Heredity, 103(4), 318–325. Caverzan, A., Piasecki, C., Chavarria, G., Stewart, C. N., & Vargas, L. (2019). Defenses Against ROS in Crops and Weeds: The Effects of Interference and Herbicides. International Journal of Molecular Sciences 2019, Vol. 20, Page 1086, 20(5), 1086. Chen, J., Huang, H., Wei, S., Huang, Z., Wang, X., and Zhang, C. (2017). Investigating the mechanisms of glyphosate resistance in goosegrass (Eleusine indica (L.) Gaertn.) by RNA sequencing technology. The Plant Journal, 89(2), 407-415. Chen, J., Huang, H., Zhang, C., Wei, S., Huang, Z., Chen, J., and Wang, X. (2015). Mutations and amplification of EPSPS gene confer resistance to glyphosate in goosegrass (Eleusine indica). Planta, 242(4), 859–868. Chen, S., Zhou, Y., Chen, Y., and Gu, J. (2018). fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics (Oxford, England), 34(17), i884–i890. Chen, Y., Lun, A. T. L., and Smyth, G. K. (2016). From reads to genes to pathways: differential expression analysis of RNA-Seq experiments using Rsubread and the edgeR quasi- likelihood pipeline. F1000Research, 5, 1438. Craig, N. L. (2020). Mobile DNA III. Danecek, P., Bonfield, J. K., Liddle, J., Marshall, J., Ohan, V., Pollard, M. O., Whitwham, A., Keane, T., McCarthy, S. A., Davies, R. M., and Li, H. (2021). Twelve years of SAMtools and BCFtools. GigaScience, 10(2), 1–4. Dao, T. T. H., Linthorst, H. J. M., & Verpoorte, R. (2011). Chalcone synthase and its functions in plant resistance. Phytochemistry Reviews 2011 10:3, 10(3), 397–412. Deng, B., Wang, W., Ruan, C., Deng, L., Yao, S., and Zeng, K. (2020). Involvement of CsWRKY70 in salicylic acid-induced citrus fruit resistance against Penicillium digitatum. Horticulture Research 2020 7:1, 7(1), 1–12. Duhoux, A., Carrère, S., Gouzy, J., Bonin, L., and Délye, C. (2015). RNA-Seq analysis of rye- grass transcriptomic response to an herbicide inhibiting acetolactate-synthase identifies transcripts linked to non-target-site-based resistance. Plant Molecular Biology, 87(4–5), 473–487. Duke, S. O., and Powles, S. B. (2008). Glyphosate: a once-in-a-century herbicide. Pest Management Science, 64(4), 319–325. 73 Dupeyron, M., Singh, K. S., Bass, C., and Hayward, A. (2019). Evolution of Mutator transposable elements across eukaryotic diversity. Mobile DNA, 10(1), 12. Eceiza, M. V., Gil-Monreal, M., Barco-Antoñanzas, M., Zabalza, A., and Royuela, M. (2022). The moderate oxidative stress induced by glyphosate is not detected in Amaranthus palmeri plants overexpressing EPSPS. Journal of Plant Physiology, 153720. Ernst, A. M., Jekat, S. B., Zielonka, S., Müller, B., Neumann, U., Rüping, B., Twyman, R. M., Krzyzanek, V., Prüfer, D., and Noll, G. A. (2012). Sieve element occlusion (SEO) genes encode structural phloem proteins involved in wound sealing of the phloem. Proceedings of the National Academy of Sciences of the United States of America, 109(28). Friesen, L. F., Beckie, H. J., Warwick, S. I., and van Acker, R. C. (2009). The biology of Canadian weeds. 138. Kochia scoparia (L.) Schrad. Canadian Journal of Plant Science, 89(1), 141–167. Gaines, T. A., Lorentz, L., Figge, A., Herrmann, J., Maiwald, F., Ott, M. C., Han, H., Busi, R., Yu, Q., Powles, S. B., and Beffa, R. (2014). RNA-Seq transcriptome analysis to identify genes involved in metabolism-based diclofop resistance in Lolium rigidum. The Plant Journal, 78(5), 865–876. Gaines, T. A., Zhang, W., Wang, D., Bukun, B., Chisholm, S. T., Shaner, D. L., Nissen, S. J., Patzoldt, W. L., Tranel, P. J., Culpepper, A. S., Grey, T. L., Webster, T. M., Vencill, W. K., Sammons, R. D., Jiang, J., Preston, C., Leach, J. E., and Westra, P. (2010). Gene amplification confers glyphosate resistance in Amaranthus palmeri. Proceedings of the National Academy of Sciences of the United States of America, 107(3), 1029–1034. Geddes, C. M., and Sharpe, S. M. (2022). Crop yield losses due to kochia (Bassia scoparia) interference. Crop Protection, 157, 105981. Ghodke, P., Khandagale, K., Thangasamy, A., Kulkarni, A., Narwade, N., Shirsat, D., Randive, P., Roylawar, P., Singh, I., Gawande, S. J., Mahajan, V., Solanke, A., and Singh, M. (2020). Comparative transcriptome analyses in contrasting onion (Allium cepa L.) genotypes for drought stress. PLoS ONE, 15(8). Giacomini, D. A., Gaines, T., Beffa, R., and Tranel, P. J. (2018). Optimizing RNA-seq studies to investigate herbicide resistance. Pest Management Science, 74(10), 2260–2264. Gill, S. S., Anjum, N. A., Gill, R., Yadav, S., Hasanuzzaman, M., Fujita, M., Mishra, P., Sabat, S. C., and Tuteja, N. (2015). Superoxide dismutase—mentor of abiotic stress tolerance in crop plants. Environmental Science and Pollution Research, 22(14), 10375–10394. Gomes, M. P., and Juneau, P. (2016). Oxidative stress in duckweed (Lemna minor L.) induced by glyphosate: Is the mitochondrial electron transport chain a target of this herbicide? Environmental Pollution, 218, 402–409. 74 Green, J. M., and Owen, M. D. K. (2010). Herbicide-Resistant Crops: Utilities and Limitations for Herbicide-Resistant Weed Management. Journal of Agricultural and Food Chemistry, 59(11), 5819–5829. Hastings, P. J., Bull, H. J., Klump, J. R., and Rosenberg, S. M. (2000). Adaptive Amplification. Cell, 103(5), 723–731. Heap, I. (2022). The International Herbicide-Resistant Weed Database. www.weedscience.org Herrmann, K. M. (1995). The Shikimate Pathway: Early Steps in the Biosynthesis of Aromatic Compounds. The Plant Cell 7. American Society of Plant Physiologists. Jaswal, R., Rajarammohan, S., Dubey, H., Kiran, K., Rawal, H., Sonah, H., Deshmukh, R., & Sharma, T. (2021). A kiwellin protein-like fold containing rust effector protein localizes to chloroplast and suppress cell death in plants. BioRxiv. Jensen, M. K., Kjaersgaard, T., Nielsen, M. M., Galberg, P., Petersen, K., O’Shea, C., and Skriver, K. (2010). The Arabidopsis thaliana NAC transcription factor family: structure– function relationships and determinants of ANAC019 stress signaling. Biochemical Journal, 426(2), 183–196. Jha, U. C., Nayyar, H., Jha, R., Khurshid, M., Zhou, M., Mantri, N., and Siddique, K. H. M. (2020). Long non-coding RNAs: emerging players regulating plant abiotic stress response and adaptation. BMC Plant Biology, 20(1), 1–20. Jiang, L. X., Jin, L. G., Guo, Y., Tao, B., and Qiu, L. J. (2013). Glyphosate effects on the gene expression of the apical bud in soybean (Glycine max). Biochemical and Biophysical Research Communications, 437(4), 544–549. Josephs, E. B., van Etten, M. L., Harkess, A., Platts, A., & Baucom, R. S. (2021). Adaptive and maladaptive expression plasticity underlying herbicide resistance in an agricultural weed. Evolution Letters, 5(4), 432–440. Jugulam, M., and Shyam, C. (2019). Non-Target-Site Resistance to Herbicides: Recent Developments. Plants, 8(10), 417. Kim, D., Paggi, J. M., Park, C., Bennett, C., and Salzberg, S. L. (2019). Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nature Biotechnology 2019 37:8, 37(8), 907–915. Kolde, R. (2019). Package “pheatmap.” https://cran.r- project.org/web/packages/pheatmap/index.html Kumar, V., Jha, P., Jugulam, M., Yadav, R., and Stahlman, P. W. (2019). Herbicide-resistant kochia (Bassia scoparia) in North America: A review. Weed Science, 67(1), 4–15. 75 Law, J. A., and Jacobsen, S. E. (2010). Establishing, maintaining and modifying DNA methylation patterns in plants and animals. Nature Reviews Genetics, 11(3), 204–220. Leslie, T., and Baucom, R. S. (2014). De Novo Assembly and Annotation of the Transcriptome of the Agricultural Weed Ipomoea purpurea Uncovers Gene Expression Changes Associated with Herbicide Resistance. G3 Genes|Genomes|Genetics, 4(10), 2035-2047. Li, F., and Wang, A. (2018). RNA decay is an antiviral defense in plants that is counteracted by viral RNA silencing suppressors. PLOS Pathogens, 14(8), e1007228. Liu, N., Zhong, G., Zhou, J., Liu, Y., Pang, Y., Cai, H., and Wu, Z. (2019). Separate and combined effects of glyphosate and copper on growth and antioxidative enzymes in Salvinia natans (L.) All. Science of The Total Environment, 655, 1448–1456. Lorentz, L., Gaines, T. A., Nissen, S. J., Westra, P., Strek, H. J., Dehne, H. W., Ruiz-Santaella, J. P., and Beffa, R. (2014). Characterization of glyphosate resistance in Amaranthus tuberculatus Populations. Journal of Agricultural and Food Chemistry, 62(32), 8134– 8142. Malone, J. M., Morran, S., Shirley, N., Boutsalis, P., and Preston, C. (2016). EPSPS gene amplification in glyphosate resistant Bromus diandrus. Pest Management Science, 72(1), 81–88. Martin, S. L., Benedict, L., Wei, W., Sauder, C. A., Beckie, H. J., and Hall, L. M. (2020). High gene flow maintains genetic diversity following selection for high EPSPS copy number in the weed kochia (Amaranthaceae). Scientific Reports 2020 10:1, 10(1), 1–11. Matzke, M. A., and Mosher, R. A. (2014). RNA-directed DNA methylation: an epigenetic pathway of increasing complexity. Nature Reviews Genetics 2014 15:6, 15(6), 394–408. McCarthy, D. J., Chen, Y., and Smyth, G. K. (2012). Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Research, 40(10), 4288–4297. Meng, D., Yu, X., Ma, L., Hu, J., Liang, Y., Liu, X., Yin, H., Liu, H., He, X., and Li, D. (2017). Transcriptomic response of Chinese yew (Taxus chinensis) to cold stress. Frontiers in Plant Science, 8, 468. Moldes, C. A., Medici, L. O., Abrahão, O. S., Tsai, S. M., and Azevedo, R. A. (2008). Biochemical responses of glyphosate resistant and susceptible soybean plants exposed to glyphosate. Acta Physiologiae Plantarum, 30(4), 469–479. Nandula, V. K., Wright, A. A., Bond, J. A., Ray, J. D., Eubank, T. W., and Molin, W. T. (2014). EPSPS amplification in glyphosate-resistant spiny amaranth (Amaranthus spinosus): a case of gene transfer via interspecific hybridization from glyphosate-resistant Palmer amaranth (Amaranthus palmeri). Pest Management Science, 70(12), 1902–1909. 76 Neve, P., and Powles, S. (2005). Recurrent selection with reduced herbicide rates results in the rapid evolution of herbicide resistance in Lolium rigidum. Theoretical and Applied Genetics, 110(6), 1154–1166. Ngo, T. D., Malone, J. M., Boutsalis, P., Gill, G., and Preston, C. (2018). EPSPS gene amplification conferring resistance to glyphosate in windmill grass (Chloris truncata) in Australia. Pest Management Science, 74(5), 1101–1108. Pan, L., Yu, Q., Han, H., Mao, L., Nyporko, A., Fan, L., Bai, L., and Powles, S. (2019). Aldo- keto Reductase Metabolizes Glyphosate and Confers Glyphosate Resistance in Echinochloa colona. Plant Physiology, 181(4), 1519–1534. Patterson, E. L., Saski, C. A., Sloan, D. B., Tranel, P. J., Westra, P., and Gaines, T. A. (2019). The Draft Genome of Kochia scoparia and the Mechanism of Glyphosate Resistance via Transposon-Mediated EPSPS Tandem Gene Duplication. Genome Biology and Evolution, 11(10), 2927–2940. Perrone, A., and Martinelli, F. (2020). Plant stress biology in epigenomic era. Plant Science, 294, 110376. Pertea, G., and Pertea, M. (2020). GFF Utilities: GffRead and GffCompare. F1000Research, 9, 304. Pettinga, D. J., Ou, J., Patterson, E. L., Jugulam, M., Westra, P., and Gaines, T. A. (2018). Increased chalcone synthase (CHS) expression is associated with dicamba resistance in Kochia scoparia. Pest Management Science, 74(10), 2306–2315. Pettengill, E. A., Pettengill, J. B., and Coleman, G. D. (2013). Elucidating the evolutionary history and expression patterns of nucleoside phosphorylase paralogs (vegetative storage proteins) in Populus and the plant kingdom. BMC Plant Biology, 13(1), 118. Preston, C., Belles, D. S., Westra, P. H., Nissen, S. J., and Ward, S. M. (2009). Inheritance of Resistance to The Auxinic Herbicide Dicamba in Kochia (Kochia scoparia). Weed Science, 57(1), 43–47. Radwan, D. E. M., and Fayez, K. A. (2016). Photosynthesis, antioxidant status and gas-exchange are altered by glyphosate application in peanut leaves. Photosynthetica, 54(2), 307–316. Reams, A. B., and Roth, J. R. (2015). Mechanisms of Gene Duplication and Amplification. Cold Spring Harbor Perspectives in Biology, 7(2), a016592. Ritz, C., Baty, F., Streibig, J. C., and Gerhard, D. (2015). Dose-Response Analysis Using R. PloS One, 10(12). 77 Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics (Oxford, England), 26(1), 139–140. Salas, R. A., Dayan, F. E., Pan, Z., Watson, S. B., Dickson, J. W., Scott, R. C., and Burgos, N. R. (2012). EPSPS gene amplification in glyphosate-resistant Italian ryegrass (Lolium perenne ssp. multiflorum) from Arkansas. Pest Management Science, 68(9), 1223–1230. Seraj, Z. I., Elias, S. M., Biswas, S., and Tuteja, N. (2018). Helicases and Their Importance in Abiotic Stresses. Salinity Responses and Tolerance in Plants, (2), 119–141. Sergiev, I. G., Alexieva, V. S., Ivanov, S. v., Moskova, I. I., and Karanov, E. N. (2006). The phenylurea cytokinin 4PU-30 protects maize plants against glyphosate action. Pesticide Biochemistry and Physiology, 85(3), 139–146. Seymour, D. K., and Becker, C. (2017). The causes and consequences of DNA methylome variation in plants. Current Opinion in Plant Biology, 36, 56–63. Shaner, D. L. (2009). Role of Translocation as A Mechanism of Resistance to Glyphosate. Weed Science, 57(1), 118–123. Sharon, A., Amsellem, Z., and Gressel, J. (1992). Glyphosate Suppression of an Elicited Defense Response Increased Susceptibility of Cassia obtusifolia to a Mycoherbicide. Plant Physiology, 98(2), 654–659. Singh, S., Kumar, V., Datta, S., Wani, A. B., Dhanjal, D. S., Romero, R., and Singh, J. (2020). Glyphosate uptake, translocation, resistance emergence in crops, analytical monitoring, toxicity and degradation: a review. Environmental Chemistry Letters 2020 18:3, 18(3), 663–702. Soares, C., Pereira, R., Spormann, S., and Fidalgo, F. (2019). Is soil contamination by a glyphosate commercial formulation truly harmless to non-target plants? – Evaluation of oxidative damage and antioxidant responses in tomato. Environmental Pollution, 247, 256–265. Supek, F., Bošnjak, M., Škunca, N., and Šmuc, T. (2011). REVIGO Summarizes and Visualizes Long Lists of Gene Ontology Terms. PLoS ONE, 6(7), e21800. Tehranchian, P., Norsworthy, J. K., Powles, S., Bararpour, M. T., Bagavathiannan, M. v., Barber, T., & Scott, R. C. (2017). Recurrent Sublethal-Dose Selection for Reduced Susceptibility of Palmer Amaranth ( Amaranthus palmeri ) to Dicamba. Weed Science, 65(2), 206–212. Tian, T., Liu, Y., Yan, H., You, Q., Yi, X., Du, Z., Xu, W., and Su, Z. (2017). agriGO v2.0: a GO analysis toolkit for the agricultural community, 2017 update. Nucleic Acids Research, 45(W1), W122–W129. 78 Toledo-Ortiz, G., Huq, E., and Quail, P. H. (2003). The Arabidopsis Basic/Helix-Loop-Helix Transcription Factor Family[W]. The Plant Cell, 15(8), 1749–1770. Tzin, V., and Galili, G. (2010). The Biosynthetic Pathways for Shikimate and Aromatic Amino Acids in Arabidopsis thaliana. The Arabidopsis Book, 8, e0132. Underwood, C. J., Henderson, I. R., and Martienssen, R. A. (2017). Genetic and epigenetic variation of transposable elements in Arabidopsis. Current Opinion in Plant Biology, 36, 135–141. van Etten, M. L., Soble, A., & Baucom, R. S. (2021). Variable inbreeding depression may explain associations between the mating system and herbicide resistance in the common morning glory. Molecular Ecology, 30(21), 5422–5437. Wang, H., and Wang, H. (2015). Multifaceted roles of FHY3 and FAR1 in light signaling and beyond. Trends in Plant Science, 20(7), 453–461. Wickham, H. 2016. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. Wiersma, A. T., Gaines, T. A., Preston, C., Hamilton, J. P., Giacomini, D., Robin Buell, C., Leach, J. E., and Westra, P. (2015). Gene amplification of 5-enol-pyruvylshikimate-3- phosphate synthase in glyphosate-resistant Kochia scoparia. Planta, 241(2), 463–474. Yu, Y., Jia, T., and Chen, X. (2017). The ‘how’ and ‘where’ of plant microRNAs. New Phytologist, 216(4), 1002–1017. Zanatta, C. B., Benevenuto, R. F., Nodari, R. O., & Agapito-Tenfen, S. Z. (2020). Stacked genetically modified soybean harboring herbicide resistance and insecticide rCry1Ac shows strong defense and redox homeostasis disturbance after glyphosate-based herbicide application. Environmental Sciences Europe, 32(1), 1–17. Zhu, J., Kapoor, A., Sridhar, V. v., Agius, F., and Zhu, J.-K. (2007). The DNA Glycosylase/Lyase ROS1 Functions in Pruning DNA Methylation Patterns in Arabidopsis. Current Biology, 17(1), 54–59. 79 CHAPTER III EPIGENOMIC RESPONSE TO GLYPHOSATE IN KOCHIA SCOPARIA Summary The epigenome is a key regulator of gene expression in plants and DNA methylation is the most common mechanism. Methylation alters gene expression by preventing transcriptional machinery from accessing the DNA and silences transposons. The role of the epigenome in transient responses to abiotic stresses such as temperature, drought, and salinity has been characterized in depth in model species. However, herbicides are another stress that plants encounter and their effect on the epigenome is not well understood, particularly in weed species. In this study, we investigated the effects of two different sublethal doses of glyphosate on DNA methylation patterns of Kochia scoparia through whole genome bisulfite sequencing before and three weeks after treatment. We characterized methylation levels of cytosines in different sequence contexts (CG, CHG, and CHH) genome wide as well as exclusively gene bodies, transposons, and other features of the genome. Differentially methylated regions were also identified to determine treatment-specific directed changes in methylation and were correlated to gene expression data from the same plant population. Our results indicate that there are no global changes in cytosine methylation in kochia following glyphosate treatment. Overall responses were stochastic, but we observed more variability in the epigenome of high dose treated plants. Additionally, differentially methylated genes tended to be dose-dependent rather than treatment specific. There was also no direct correlation between changes in methylation and changes in gene expression, so other mechanisms may be at play. This work aids in understanding the role of the epigenome in weedy 80 species resilience to herbicide stress and suggests the need for further research connecting DNA methylation and herbicide resistance evolution. Introduction The epigenome - DNA methylation maintenance and regulation of gene expression The epigenome consists of histone variants, histone post-translational modifications, and non-coding RNAs that alter the chromatin state (Zhang et al., 2018). It plays a significant role in regulating plant development and stress response and the silencing of transposons (Perrone and Martinelli, 2020). DNA methylation is the most common and easily studied mechanism, which occurs when a methyl group is added to the 5th position of cytosines. DNA methylation influences the ability of transcriptional machinery to access the DNA and can either recruit or prevent transcription factors; therefore, modulating gene expression (Law and Jacobsen, 2010). Much of what we know about epigenomics in plants comes from Arabidopsis thaliana in which epigenetic recombinant lines (epiRILs) were inbred for multiple generations to study epigenetic stability and inheritance (Zhang et al., 2018). Inheritance of epigenetic alleles (epialleles), or DNA sequence with differing methylation can sometimes lead to drastic changes in the phenotypes of the plants, however not all methylation events equally contribute to phenotypic changes. DNA methylation is a highly conserved process and each cytosine context (CG, CHG, CHH where H = A, T, C) is established and maintained by different mechanisms (Bewick et al., 2017). Each context serves a different purpose as well, where CG methylation is commonly found on gene bodies, CHG and CHH methylation work to silence transposons and repetitive regions of the genome (Zhang et al., 2018). In plants, de novo methylation in all three contexts is established by the RNA-directed DNA Methylation (RdDM) pathway. RdDM involves small 81 interfering RNAs that target specific sequences and recruit DNA methyltransferases such as Domains Rearranged Methyltransferase 2 (DRM2) for transcriptional silencing, particularly the edges of transposons (Matzke and Mosher, 2014; Niederhuth and Schmitz, 2017). CG methylation is symmetrical, meaning it is maintained on both the forward and reverse strand of DNA, and is established by methyltransferase 1 (MET1) (Niederhuth and Schmitz, 2017). This is important to note, as methylation in this context is maintained across cell divisions. CG methylation in gene bodies often results in increased expression, though there is a notable decrease of methylation near the transcription start and termination sites (TSS and TTS, respectively) (Law and Jacobsen, 2010; Niederhuth and Schmitz, 2017). CHG methylation is also symmetrical but relies on a pathway involving Chromomethylase 3 (CMT3) to maintain methylation after DNA replication. CHH on the other hand is asymmetrical and must be continually re-established after each cell generation by Chromomethylase 2 (CMT2) (Niederhuth and Schmitz, 2017; Sasaki et al., 2019; Wang et al., 2015). Both non-CG contexts are related to plant development and environmental responses (Niederhuth and Schmitz, 2017). Interestingly, findings in Arabidopsis suggest that CG sites tend to be bimodal, either they are highly methylated (80-100%) or not (Cokus et al., 2008; Wang et al., 2015). In the CHH context, sites are either methylated at just 10% or remain unmethylated, while CHG can range from 20-100% (Cokus et al., 2008). This is important to consider when analyzing bisulfite sequencing data. Genome wide methylation levels by context also vary between species (Niederhuth et al., 2016), and specific cell types have their own distinct patterns of methylation as well (Kawakatsu et al., 2016). Active demethylation is performed both by DNA glycosylases in the Demeter (DME) family and by the Repressor of Silencing (ROS1) 82 family of proteins, sometimes to counteract methylation from members of the RdDM protein family (Law and Jacobsen, 2010). Epigenetics and plant stress Under stressed conditions, such as extreme temperature, water stress, salinity, or ultraviolet (UV) stress, the methylation state across the genome may be altered to modulate gene expression; more genes can be “turned on” to counteract the stress and allow the plant to acclimate in novel ways (Bilichak et al., 2012). Abiotic stress has been reported to induce stochastic methylation (Ashapkin et al., 2020). The change in methylation pattern is highly dependent upon the stress and the species being tested. For instance, Arabidopsis responded with hypermethylation to phosphate starvation (Bilichak et al., 2012). Epigenetic priming and stress memory can also be observed, where plants exposed to a particular stress are more prone to a faster or more rigorous response when encountering the same stimulus again, which can also be passed onto progeny (Lämke and Bäurle, 2017). Transposable elements The epigenome contributes heavily to the silencing of transposable elements (TEs), which are constitutively highly methylated in the genome as their activity can be harmful (Weil and Martienssen, 2008). However, demethylation can result in increased transposition of TEs, which may alter genome variability and the expression of other genes as well (Steward et al., 2002). In fact, Arabidopsis mutants with a loss of CG methylation maintenance displayed increased expression of transposons and increased transposition of Copia and Gypsy elements, as well as some Mutator TEs (Underwood et al., 2017). In response to heat stress, a particular copia-like TE ONSEN became transcriptionally active and produced multiple copies in extrachromosomal DNA (Ito et al., 2011). This may indicate that transposon activation is part of 83 the general plant stress response. Additionally, methylation levels of TEs tended to be higher in Gypsy and Copia elements above other TE families (Wang et al., 2015). Herbicide resistance and epigenomics Previous studies have shown that herbicides alter the plant methylome. In the case of Oryza sativa exposed to the photosystem II inhibitor atrazine, increased global methylation levels by approximately 2%, mostly attributed to the CHG and CHH contexts (Lu et al., 2016). In addition, differentially methylated regions were preferentially hypermethylated, especially in regions upstream and downstream of gene bodies. The results of this study indicated epigenetic regulation of genes associated with atrazine degradation and detoxification (Lu et al., 2016). Stress specific and dose-dependent methylation patterns were observed in Arabidopsis following sublethal glyphosate exposure, further supporting the role of the epigenome in plant response to herbicides (Kim et al., 2017). Only 7% of the differentially methylated regions detected from glyphosate treatment were shared with phosphate starvation and biotic stresses, meaning glyphosate treatment induced unique changes in methylation (Kim et al., 2017). More recently, a comparison of DNA methylation responses to glyphosate between tolerant and susceptible inbred Zea mays lines revealed a stark increase in methylation in the susceptible line 6 hours after treatment (>18%) (Tyczewska et al., 2021). Though, this increase in methylation level was not detected seven days after treatment where levels returned to control conditions. Among the differentially methylated regions they identified were transporter proteins, methyltransferases, and transposons (Tyczewska et al., 2021). There is a lack of research eliciting herbicide response patterns in methylation on their primary target: weeds. The mechanism of glyphosate resistance in multiple plants species is due to copy number variation whereby multiple copies of the target site allow the plants to survive 84 field rates (Patterson et al., 2019). In Kochia scoparia, this is caused by TE activity which commonly results in gene duplication events after transposition (Hastings et al., 2000; Patterson et al., 2019; Reams and Roth, 2015). Whether or not the epigenome is involved in regulating this process is particularly important for understanding the evolution of herbicide resistance, especially to glyphosate. Kochia (Bassia scoparia (L.) A.J. Scott syn. Kochia scoparia (L.) Schrad.) is an invasive weed that has ravaged grasslands, desert drylands and agricultural fields throughout the U.S. Great Plains and Canada (Friesen et al., 2009). Yield losses reported from kochia can be greater than 90% in sugar beets, corn, sorghum, and sunflower (Geddes and Sharpe, 2022). This C4 tumbleweed can withstand drought and salinity, is a prolific seed producer, and emerges earlier in the spring than other species, giving it a competitive advantage in the climates it invades (Friesen et al., 2009). It has also evolved resistance to multiple herbicides including photosystem II (PSII) inhibitors, acetolactate synthase (ALS) inhibitors, synthetic auxins, and glyphosate (Heap, 2022; Kumar et al., 2019). Glyphosate resistant kochia populations have been reported in 10 states in the U.S. and multiple locations in Canada (Kumar et al., 2019). Therefore, to understand the effect of herbicides on the methylome of the important weed species kochia, we used whole genome bisulfite sequencing to determine transient changes three weeks after sublethal glyphosate application. We hypothesized that glyphosate treatment would induce hypomethylation of the genome to allow for the activation of abiotic stress responsive genes and transposons, with the higher dose leading to more drastic changes in methylation (Busi and Powles, 2009). Our work aims to address how herbicides shape the epigenome to fully capture the physiological implications of modern weed control. Discoveries into the plasticity of 85 weed genomes in response to human selection will have ramifications for weed managers combating resistance and for researchers studying plant epigenomes alike. Methods Plant material and sublethal dose determination for experiment The previously characterized glyphosate susceptible kochia line ‘7710’ was used in this experiment to ensure homogenous genetic background (Patterson et al., 2019; Pettinga et al., 2018; Preston et al., 2009). Kochia seedlings were grown in a growth chamber at 25°C/23°C (day/night) with a photoperiod of 16 hours, and a light intensity of 600 µmol photon m-2 s-2. To determine appropriate dosage levels, nine glyphosate doses from 0-2x the standard field rate of 840 g ae ha-1 (0, 13, 26, 52.5, 105, 210, 420, 840 and 1,680 g ae ha-1) were assayed with three replicates of 18 seedlings for each dose (Wiersma et al., 2015). The formulation utilized in this study was Roundup PowerMAX (540 g ae L−1, Bayer Crop Science) applied with ammonium sulfate (2% w/v) was applied to plants averaging 10 cm in height using the Generation 4 DeVries Research Track Sprayer in the Michigan State University Plant and Pest Diagnostics Pesticide Application Laboratory. The sprayer was calibrated to deliver 187 L ha-1 spray volume at 207 kPa pressure. Plants were then moved to the greenhouse at a temperature of 29.4°C/23.8°C (day/night) with a 16 h photoperiod. All aboveground tissue was collected three weeks after treatment and dried in an Equatherm Oven (Curtin Matheson Scientific, Inc.) at 70°C for 48 hours. Dry weight was the response variable used to generate the dose response curve using a four-parameter log-logistic model in the drc package in R, which was selected using the mselect function (Fig. 2.1) (Ritz et al., 2015). We tested both three- and four-parameter models and chose the one with the smallest information criterion and residual standard error, as well as the largest 86 P value for lack-of-fit test. The median effective dose (ED50) was 92.5 g ae ha-1 and ED25 was 75.7 g ae ha-1. We saw a hormetic effect around the ED25 (Anunciato et al., 2022). The two sublethal doses for this study were chosen as 52.5 and 105 g ae ha-1 (1/16th and 1/8th of a field rate) to ensure we had two “extreme” doses, without impairing the plants’ ability to set seed for future studies (Fig. 2.1). Experimental design and leaf tissue collection Kochia plants from the same line ‘7710’ were grown in individual pots in the greenhouse at 29.4°C/23.8°C with a 16 h photoperiod until plants averaged 10 cm tall. Approximately 120 plants per treatment were sprayed as described above with either 0, 52.5 (“low”), or 105 (“high”) g ae ha-1 of glyphosate. Young leaf tissue was collected from four individuals per treatment pre- application (on the date of spraying) and was flash-frozen in liquid nitrogen and stored at -80°C until use. After application, plant height was recorded up to two weeks after treatment (Fig. 2.2A). Since visual glyphosate injury symptomology has a long effect, newly expanded leaf tissue three weeks post-application was collected from four individuals per treatment and handled in the same manner as pre-treatment samples. DNA sample preparation and sequencing Genomic DNA was extracted from both pre and post glyphosate tissue samples using the DNeasy Plant Mini kit (Qiagen). DNA concentration, purity and quality were assessed with a NanoDrop One spectrophotometer (Thermo Scientific), and quantity was assessed with a Qubit 1.0 fluorometer (Invitrogen). DNA samples (at least 1.0 µg per sample) were sent to Novogene Co., Ltd., where bisulfite treatment and library construction took place. Illumina NovaSeq6000 was used for paired end (150bp) whole genome bisulfite sequencing. 87 Data processing, mapping, and methylation extraction Sequence reads were filtered and trimmed using fastp with default parameters (Chen et al., 2018). Next, we evaluated two different pipelines for calling cytosine methylation sites across the genome including Bismark Bisulfite Mapper and the Python-based pipeline methylpy (Krueger and Andrews, 2011; Schultz et al., 2015). Bismark mapping efficiencies for each sample were only about 52%, and there was a remarkable difference in the methylation calls when compared to methylpy. For example, Bismark reported samples with as high as 81.0%, 62.0%, and 9.9% methylation of cytosines in the CG, CHG, and CHH contexts, respectively. The same sample when using the methylpy pipeline had 69.3%, 42.7%, and 8.3% methylation of cytosines in the same contexts. When looking at the number of methylated cytosines, Bismark called twice as many sites compared to methylpy. This is a stark difference which we attributed to more robust methylation calls, and we decided to use methylpy for the remainder of this study to process our paired-end bisulfite sequencing data. The existing kochia reference genome was used to build a converted reference for the aligner Bowtie within methylpy (Schultz et al., 2015). We then used the paired end pipeline to align the reads from each sample to the converted kochia genome and call the methylation level of all cytosines (Schultz et al., 2015). The unmethylated control was set to lambda DNA sequence that was spiked in during sequencing, which is used to estimate the bisulfite non- conversion rate (Schultz et al., 2015). Weighted methylation levels genome-wide and by different genomic features The output of the methylpy paired end pipeline is an ‘allc’ file containing the position, sequence context, number of methylated reads, and total read coverage of every cytosine in the genome. We created a custom python script to calculate the weighted methylation levels of all 88 cytosines in each sequence context by dividing the number of methylated reads by the total number of reads at each position. We also were interested in different genomic features. To analyze these, we separated individual annotation files into just genes, coding sequence (CDS), mRNA, exons, introns, and a file which included genes plus 2kb upstream and downstream defined using bedtools slop (Quinlan and Hall, 2010). We then used bedtools intersect to extract cytosines in the allc files for each sample that correspond with the various genome features. Weighted methylation levels of each feature were calculated using the same python script described above. Treatment effect differences in weighted methylation levels were determined in Microsoft Excel by a two-tailed student’s T-Test assuming heteroscedastic sample variance to calculate P values (Microsoft Corporation, 2018). Pattern of CG, CHG and CHH methylation in genes Allc files were used to visualize cytosine methylation level patterns upstream, downstream, and at the TSS and TTS of all genes in the entire genome by different sequence contexts. This was done using Custom python scripts modified from Dr. Chad Niederhuth at test- data/functions.py at master · niederhuth/test-data · GitHub. Analysis of differentially methylation regions Differentially methylated regions (DMRs) were called by methylpy DMRFind, which first identifies differentially methylated sites (DMS) and collapses those to define the DMR (Schultz et al., 2015). Average methylation of each DMR across the four replicates was used to calculate Log2Fold Change (Log2FC) between pre and post samples. Next a two-tailed student’s T-Test assuming heteroscedastic variance was used to determine P values (Microsoft Corporation, 2018). The significance threshold we set for DMRs was P value < 0.05 and 89 absolute Log2FC > 0.5, where Log2FC > 0.5 is hypermethylated, and Log2FC < -0.5 is hypomethylated. Volcano plots generated with ggplot2 in R were used to visualize the hypermethylated and hypomethylated DMRs in each sequence context (Wickham, 2016). We then defined differentially methylated genes (DMGs) as having one or more significant DMRs within 2kb upstream or downstream of the genes and used Venn diagrams to determine shared versus unique DMGs between treatments (VIB / UGent Bioinformatics & Evolutionary Genomics, 2022). The kochia genome annotation file was used to analyze the function of the top DMGs specific to glyphosate treatment and those that were dose dependent (Patterson et al., 2019). High dose specific karyotypes To visualize the locations of significant DMGs specific to high dose treatment, the package RIdeogram in R was used to generate three different karyotypes of the kochia genome, with differentially methylated sites in each context (Hao et al., 2020). Average methylation was converted into binary values for better visualization where zero is hypomethylated and one is hypermethylated. Unscaffolded chromosomes are not included in the karyotypes. Correlation with previously generated gene expression data The fold change in methylation was correlated with expression data generated from the same generation of kochia plants (see Chapter 2 for methods) to determine whether increased methylation resulted in decreased gene expression and vice versa. We compared differentially methylated genes to all differentially expressed genes, as well as only genes with significant differential expression (adjusted P value < 0.05 and |Log2FC > 1.0|). To test the null hypothesis that the slope of the regression line is 0, we used a T-Test of the slope in Excel to determine P values. 90 Transposon families Additionally, we used RepeatModeler with the LTRStruct flag to analyze the kochia genome for TEs (Smit and Hubley, 2008). We split them by family to analyze Copia, Gypsy, and MULE transposons specifically. We used bedtools intersect and the custom weighted methylation python script once again to determine changes in methylation over time by treatment. We then ran DMRFind to make comparisons between only post samples of untreated, low dose and high dose plants. We then intersected these DMRs with the transposon annotations to determine differential methylation of Copia, Gypsy or MULE TEs between treatments. Results Genome wide DNA methylation patterns in kochia As the first study to characterize the kochia methylome, we needed to quality check our bisulfite sequencing data and determine baseline methylation levels. For each sample, between 135 and 199 million reads were generated averaging approximately 22 Gb, which equates to 22x coverage of the kochia genome. The total number of cytosine sites in each sequence context ranged from 18.7 – 20.2 million CG’s, 17.9 – 19.3 million CHG’s and 120.3 – 129.6 million CHH sites per sample. The bisulfite conversion rate was 99.7% across all samples. All sequencing summary statistics are reported in Table 3.1. To determine the effect of glyphosate on global methylation, we compared the weighted methylation levels of all cytosines in the genome by their sequence context across treatments. In the pre samples, the highest methylation levels of cytosines were found in the CG context (68.5%), followed by CHG (41.2%) and CHH (8%) (Fig. 3.1). There were no significant differences in CG or CHG methylation between the pre and post samples in any treatment. However, we observed a significant increase in CHH methylation by 2.5% in the untreated post 91 samples (P < 0.01), which we are associating with natural developmental changes that occurred over the period of three weeks between sampling times. We see a similar significant increase in CHH methylation by 2% (P < 0.001) when plants are sprayed with the low dose of glyphosate; but there is no significant difference in CHH methylation in the high dose treatment. There seems to be higher variability in the high dose, post samples based on increased standard error in all sequence contexts. Since methylation in each sequence contexts plays key roles in different parts of the genome, we also analyzed weighted methylation values of gene bodies, genes plus 2kb upstream and downstream, CDS, mRNA, exons, and introns (Fig. 3.2). There were no differences in CG or CHG methylation between pre and post samples regardless of treatment in any of these regions. Additionally, we observed the same trend of slight but significant increases in CHH methylation in both untreated and low dose treated plants regardless of the genomic region being analyzed. Though, there were no significant changes in any sequence context in the high dose treatment. DNA methylation patterns of genes To visualize the methylation levels of all genes plus their upstream and downstream regions shown in Fig. 3.2A, we used metaplots to observe the distinct pattern of methylation in each sequence context across treatments (Fig. 3.3). Across all treatments, CG methylation is high across the genome and remains high in the gene bodies, but there is a sharp decrease at transcription start sites which is necessary for that gene to be expressed (Fig. 3.3A). There are few differences in this pattern across treatments or between time points, though the high post samples seem to continue to have greater variability. On the contrary, we see that CHG methylation is high upstream of the gene (about 45%), low across the entire gene body (10- 20%), and at intermediate levels downstream (35%) across all treatments (Fig. 3.3B). CHH 92 methylation is constitutively maintained at extremely low levels and decreases at gene bodies, which is seen regardless of treatment (less than 6%) (Fig. 3.3C). We also see an increase in CHH methylation across treatments post application, especially upstream and downstream from the genes, which was observed by the weighted methylation level (Fig. 3.2A), though we can see that one high post plant is behaving similar to the pre plants. This might explain why there is no statistical significance for the increase in CHH methylation under high treatment. Analysis of differentially methylated regions and their associated genes In order to determine if there was a directed response in DNA methylation to glyphosate treatment, we compared differentially methylated regions between pre and post timings for each dose and determined the genes these DMRs were associated with. In the CG context, there were 5,244 DMRs in untreated, 3,561 in the low dose treatment, and 4,619 in the high dose (Fig. 3.4). DMRs in this context were preferentially hypomethylated in post samples compared to pre across all three treatments (Fig. 3.4). With our differentially methylated gene (DMG) criteria and since there can be multiple DMRs within a gene, these represented 1,881 DMGs in untreated, of which 1,848 (98.2%) were hypomethylated and 33 (1.8%) were hypermethylated (Fig. 3.5A, B). There were 1,202 DMGs in low, of which only 24 (2.0%) were hypermethylated. Interestingly, a greater proportion of DMGs in high treatment were hypermethylated in the CG context, 95 out of 814 total (11.7%) (Fig. 3.5A, B). Next, we wanted to determine the genes that were most regulated through DNA methylation following treatment. There were 88 hypermethylated genes in the CG context unique to high treatment (Fig. 3.5A). The top hypermethylated genes with the greatest increase in Log2FC included transcription factors, a cytochrome P450, and hydroxyproline O- galactosyltransferase involved in growth and development (Table 3.2). The most hypomethylated 93 genes were related to auxin transport, potassium transport, very-long-chain-fatty acid elongation and a transcription factor involved in abiotic stress response (Table 3.2). Next, we analyzed differential methylation in the CHG context, which is responsible for silencing both genes and TEs. In the CHG context, there were less significant DMRs overall compared with the CG context across all treatments (Fig. 3.6). These represent a total of 2,911 DMRs in untreated, 2,193 in the low dose treatment of glyphosate, and 2,283 DMRs in the high dose. The total number of DMGs represented were 407, 265, and 204 in untreated, low, and high treatments (Fig. 3.5C, D). Of these, 16% were hypermethylated in untreated, 18.5% in low dose treatment, and 47.1% in high. There were 83 and 84 genes exclusively hypermethylated and hypomethylated, respectively, in the high dose treatment (Fig. 3.5C, D). Of these genes, those that were hypermethylated included multiple uncharacterized proteins, protein kinases, a heat shock protein and O2 thioredoxin (Table 3.3). Those that were most hypomethylated also consisted of uncharacterized genes, but also a psbP-like protein involved in photosystem II repair, calcium and potassium transport, and a glutathione-S-transferase. Finally, in the CHH context which would primarily silence TEs, cytosines are preferentially hypermethylated across all treatments as shown in the volcano plot (Fig. 3.7). Significant DMRs hypermethylated in the CHH context represented the greatest number of genes exclusive to high dose treatment (1,672, Fig. 3.5E). Those genes with the greatest increase in methylation unique to high treatment included an ABC transporter, a condensin protein, and a DNA repair enzyme (Table 3.4). Interestingly, there were no hypomethylated DMGs shared between untreated and high sublethal glyphosate treatment in this context (Fig. 3.5F). 94 Hypomethylated genes unique to high treatment included an F-Box protein, ethylene responsive transcription factor and zinc finger protein (Table 3.4). To determine the locations in the genome where these DMRs are located, we created karyotypes of the high dose specific, significant DMRs. We saw the preferential hypomethylation of CG sites, which are enriched in the euchromatic arms, particularly on chromosomes 3 and 4 (Fig. 3.8). There were far fewer DMRs in the CHG context, and they are dispersed more evenly throughout the genome. The preferential hypermethylation of CHH is obvious, and it is the most common in heterochromatic regions where it is responsible for silencing TEs (also the most common context overall based on number of sites). CHH methylation displays the most dynamic response to treatment where it can be continually re- established de-novo even after treatment. Changes in DNA methylation do not directly correlate to a change in expression We further identified the correlation between the change in methylation and differential gene expression using data from the same population of kochia plants sprayed with sublethal glyphosate. When comparing all significantly differentially methylated genes to their subsequent change in expression from our data set, regardless of whether or not they were signifigantly differentially expressed, there was no correlation with high dose treatment (Fig 3.9). Additionaly, the only regression line significantly different from zero was in the CHH context. However, when we excluded the comparison to only those DEGs that also met our significance thresholds for change in expression of absolute Log2FC > 1.0 and P value < 0.05, there was an increase in correlation where the strength of the relationship is greater according to the R2 value. The regression line in the CG context was significanlty different from zero in this case. 95 Methylation of transposons We were especially interested in methylation of Copia, Gypsy and MULE transposon families as they are agents of structural variation generation; however, there were no significant changes in methylation in any of the three contexts in the high dose treated plants between pre and post timings for these TE-families. It was surprising to find that, in all three TE-families, untreated plants exhibited significant changes (pre to post) in methylation in all three sequence contexts. In the case of Copia transposons, methylation increased by 1.7% in CG, 2.2% in CHG, and by 3.2% in CHH (Fig. 3.10A). Low dose plants only responded to treatments with an increase in CHH methylation by 3.1%, but changes in the other contexts were not significant. We also directly compared post-treatment plants and there was also no significant difference between post samples across all treatments. In the case of both Gypsy and MULE transposons, untreated and low dose plants responded with increased methylation in all three contexts over time, but high dose treated plants did not (Fig. 3.10B, C). There were no significant differences between post treatments when looking at Gypsy TEs, but the CHG comparison between untreated post and low post samples was significant when analyzing MULEs. The CHG methylation of MULEs increased by 2.1% in untreated, compared to 1.4% in low treatment (P <0.05). We then analyzed the number of DMRs between untreated and treated plants that intersect with TEs in each family. Unsurprisingly, DMRs in the CHH context were most common when analyzing TEs, particularly the Gypsy retrotransposons (Table 3.5). The greatest amount of treatment specific DMRs were in the Copia family, where more than 11% of DMRs were shared between low and high sublethal glyphosate treatments in the CG and CHG contexts. 96 Only 5% of those in the CHH context in each family were treatment specific, indicating a potential directed response in methylation patterns to glyphosate treatment. Discussion Here we have described the constitutive DNA methylation status of Kochia scoparia leaves and determined the effects of sublethal glyphosate application on the epigenome three weeks after treatment. The conversion rate of our bisulfite sequencing was exceptionally high (99.7%) (quality studies >99.5%), and when combined with 22x coverage of the genome, this gives us high confidence that our dataset truly captured changes in methylation due to glyphosate treatment for most of the genome (Table 3.1). Cytosine methylation sites in the CHH context were by far the most abundant, as it is every cytosine not in the CG or CHG contexts. Also, the TE density of the kochia genome is high, approximately 58% (Hall et al., unpublished), where CHH methylation is critical for TE silencing (Weil and Martienssen, 2008). Average methylation levels before treatment were 68.5% in the CG context, 41.2% in the CHG context and 8% in CHH (Fig. 3.1), which are similar to a close relative, Beta vulgaris, whose standing methylation levels were reported averaging 66.8% (CG), 41.7% (CHG), and 9.4% (CHH) (Gutschker et al., 2022). During untreated plant development, we observed a significant genome wide increase in CHH methylation in untreated (2.5%) and low dose treated plants (2%), but not in high dose plants (Fig. 3.1). Analyzing on specific genome features, such as genes or introns, did not change this observation significantly (Fig. 3.2). We did however see increased variance in the post samples in the high dose treatment compared to the low or untreated plants, mostly driven by one plant whose methylation did not change before to after treatment. Increased variability in plants after treatment might explain why there are no significant increases in methylation of cytosines in any sequence context as standard T-tests rely 97 on changes being consistent as well as large. In our data, we see stasis of the epigenome under high sublethal glyphosate treatment rather than normal developmental changes taking place, which is like the transcriptomic response we observed. Previous work attempting to understand glyphosate’s effect on the epigenome of treated plants have had mixed results. Findings in Arabidopsis show few glyphosate-specific long-term effects of methylation, however, sampling time in this study shifted based on the differing growth rates of the plants after treatment (Kim et al., 2017). This contrasts heavily with results from susceptible Z. mays treated with glyphosate where they observed a significant 18.65% increase in overall methylation when sampling just 6 hours after treatment (Tyczewska et al., 2021). Another herbicide, Atrazine, sprayed on O. sativa also displayed significantly increased global methylation levels by 2% six days after treatment (Lu et al., 2016). Contrary to other literature, we did not find significant changes in global methylation in response to even relatively high levels of sublethal glyphosate. This could mean that the high dose treatment stunted normal development in the plants and/or that our sampling time did not capture the full epigenomic response. Perhaps if we had taken samples during a time course (1h, 6h, 24h, 3 days, 7 days, 14 days) we could have seen a difference in the timing of these development shifts. Normal patterns of methylation were observed in genes and their promoters when averaging across all genes, where CG methylation decreases at the TSS but remains high in the gene body, and both nonCG contexts are depleted (Fig. 3.3) (Law and Jacobsen, 2010). We did not narrow down to promoter regions exclusively, but instead they were analyzed simultaneous to the genes. It is important to note that if there are a few important specific changes in methylation, they may be diluted by looking at genome wide averages and could only be uncovered by careful consideration of smaller genome subsets. The next steps would be to look 98 at specific genes and gene families. Additionally, note that methylation is often cell-type specific, and here we analyzed whole leaves with multiple cell types (Kawakatsu et al., 2016). When looking at differential methylation in specific sequence contexts, we found that CG sites were preferentially hypomethylated in all treatments, with untreated having the most differentially methylated genes (DMGs), followed by low, then high dose (Fig. 3.4, 3.5B). However, high dose treated plants had more hypermethylated DMGs than low or untreated, and when considered as a proportion of the total DMGs more still, 88 of which are high treatment specific (Fig. 3.5A)! Many TFs were in the high dose only DMR gene list and as we observed stunting and reduced developmental changes, the impacts of hypermethylation of these transcription factors are particularly interesting (Table 3.2). TFs are a good place to look for master regulators of transcriptomic response and epigenomic regulation of TF expression could have much larger effects on the transcriptome than expected (Hong, 2016). CHG sites had fewer overall changes and are the least interesting in terms of gene body methylation (Fig. 3.6, 3.5C, D). Many of the high dose-specific genes here were annotated as uncharacterized proteins (Table 3.3). 5-methylcytosines in the CHH context were preferentially hypermethylated across all treatments (Fig. 3.7). This is interesting because due to the mechanism of CHH methylation maintenance, it is prone to be passively lost over time if not continually re-established de novo after cell divisions (Sasaki et al., 2019; Zhang et al., 2018). Our results show that glyphosate treatment is not affecting this reestablishment (Fig. 3.8), which seems to go against our hypothesis, as we expected to see hypomethylation of the genome to allow for transposon activation. However, in both Arabidopsis and O. sativa, increases in CHH and CHG methylation are observed during the transition from vegetative to reproductive maturity, perhaps to prevent 99 TE activity while producing progeny (Lloyd & Lister, 2021). When taking this into account, seeing no change in CHH methylation may be attributed to the fact that high dose plants are stunted in their development. There were no CHH hypomethylated DMGs shared between untreated and high sublethal glyphosate treatment in this context (Fig. 3.5F), so these genes would be good candidates to select for continued study of glyphosate’s effect on activation of the genome. Changes in gene expression are not highly correlated with changes in methylation when looking at all significantly differentially methylated genes (Fig. 3.9). However, if we exclude only those DMGs which are also significantly differentially expressed, that correlation increases. Our data shows that an increase in CG methylation decreases expression however, the location of these CG sites with respect to the gene need to be verified. It has been shown that CG methylation in gene bodies increases gene expression but decreases expression when present in the promotor region (Zhang et al., 2018). On the contrary, CHH islands in Z. mays upstream of a gene may also indicate increased expression (Niederhuth and Schmitz, 2017). These rules are not hard and fast, and the impact of any change in methylation is case-by-case. For example, CG methylation at the promoter serves to enhance binding of transcription activators, which is the case in tomato fruit ripening (Zhang et al., 2018). Due to the complexity of epigenomics regulation there are potential limitations to correlating expression in this way. We need to look at specific genes and specific methylation sites for clear results as to what is causing significant differential expression. To our surprise, TEs were not significantly hyper- or hypo-methylated due to glyphosate treatment. Instead, the largest effect was developmental time, where untreated plants display an increase in methylation of TEs across the board as they develop (Fig. 3.10). Low dose treated 100 plants respond similarly except for Copia TEs, where they only increase CHH, not CG or CHG (Fig. 3.10A). The reasons we are not seeing changes in high dose treated plants could be due to increased variability diluting the results of the T-tests, or it could be that high plants are still at an early developmental stage and are not responding in the same manner yet. Transposons have been shown in kochia to effect glyphosate resistance evolution, especially a MULE-like FHY3/FAR1 homologue (Patterson et al., 2019). Methylation or activation of TEs is not strongly supported with our dataset, as the particular FHY3/FAR1 involved in the EPSPS gene duplication event in kochia did not display differential methylation following glyphosate treatment. Still, the differentially methylated TEs that are treatment specific would be interesting as follow ups especially if their expression changes (Table 3.5). We encountered challenges with these analyses since these TEs are not “annotated as genes” and would need special consideration to determine changes in expression. Conclusion This study elucidates the epigenomic response of kochia to sublethal glyphosate, where there are a number of dose-dependent and treatment specific differentially methylated regions that would be interesting to follow up on in future research. Overall, there are no global shifts in methylation levels in any sequence context as a result of glyphosate treatment. However, sublethal doses of glyphosate repress normal developmental shifts in CHH methylation specifically, which partially supports our hypothesis that glyphosate induces hypomethylation of the genome that results in the activation of transposons. In this case it is a lack of change in comparison to a 2.5% increase in untreated plants. We also did not find that changes in expression and methylation correlate, so future studies should work to incorporate epigenomics into their omics analyses, with this caveat in mind. It does not appear that methylation alone is 101 sufficient to capture the full picture of how plants modulate gene expression through the epigenome, and further research is needed to confirm this in other non-model species. 102 APPENDIX 103 Raw BS Total raw Effective Q20 Q30 GC Sample data conversion reads (%) (%) (%) (%) (G) rate (%) HIGH1 143,042,608 21.5 99.70 95.68 89.05 22.11 99.73 HIGH2 141,971,242 21.3 99.66 95.56 88.85 22.38 99.69 HIGH3 154,998,284 23.2 99.59 95.89 89.45 22.68 99.72 HIGH4 135,283,440 20.3 99.75 95.51 88.78 22.21 99.71 HIGHP1 142,224,710 21.3 99.66 95.06 87.90 22.77 99.74 HIGHP2 136,696,144 20.5 99.65 95.28 88.27 22.33 99.70 HIGHP3 152,556,078 22.9 99.60 95.37 88.60 23.37 99.73 HIGHP4 145,391,452 21.8 99.62 95.46 88.64 23.11 99.70 LOW1 143,177,490 21.5 99.52 95.80 89.36 22.63 99.75 LOW2 153,668,808 23.1 99.63 95.78 89.21 22.52 99.69 LOW3 199,080,600 29.9 99.64 95.76 89.22 22.28 99.74 LOW4 142,005,010 21.3 99.55 95.66 89.00 22.72 99.70 LOWP1 134,969,678 20.2 99.68 95.69 89.07 23.05 99.73 LOWP2 139,609,186 20.9 99.63 95.82 89.34 23.10 99.69 LOWP3 138,823,406 20.8 99.65 95.57 88.81 23.26 99.71 LOWP4 148,154,394 22.2 99.60 95.33 88.46 23.00 99.70 UNT1 137,413,418 20.6 99.57 95.73 89.15 23.08 99.74 UNT2 146,138,646 21.9 99.48 95.32 88.33 22.86 99.73 UNT3 161,238,672 24.2 99.62 95.69 89.03 22.20 99.73 UNT4 147,776,954 22.2 99.45 95.30 88.46 22.35 99.70 UNTP1 140,090,870 21.0 99.52 95.67 89.10 23.07 99.73 UNTP2 146,375,552 22.0 99.59 95.66 89.03 23.08 99.69 UNTP3 136,042,082 20.4 99.63 95.72 89.17 22.74 99.73 UNTP4 148,796,848 22.3 99.54 95.73 89.38 23.20 99.70 Table 3.1: Summary statistics from bisulfite sequencing data. Summary of Kochia scoparia whole genome bisulfite sequencing quality check report adapted from Novogene. DNA was extracted from young leaf tissue from four replicates each in the following samples: UNT = Untreated Pre, UNTP = Untreated Post, LOW = Low dose pre, LOWP = Low Dose Post, HIGH = High Dose Pre, and HIGHP = High Dose Post treatment. Pre samples were from 10 cm tall plants, which were then sprayed with 0, 52.5, or 105 g ae h-1 of glyphosate and post samples were collected three weeks after treatment. Total raw reads is the number of reads generated, raw data is the number of raw reads multiplied by their sequence length (PE150), effective is the percentage of clean reads out of the total raw reads, Q20 and Q30 are the base count of Phred value > 20 or 30 out of the total base count, GC content is listed as a percentage of total bases, and the BS conversion rate is the percentage of T changed from C by bisulfite treatment. 104 Figure 3.1: Genome wide methylation. Genome wide methylation level of all cytosines in the CG, CHG, and CHH contexts identified with bisulfite sequencing of kochia leaf tissue before (Pre) and three weeks after treatment (Post) with 0 (Untreated), 52.5 (Low) and 105 (High) g ae ha-1 of glyphosate. Plants were grown in a greenhouse setting and included four replicates from each of the treatments both pre and post application timings. A two-tailed student’s T-Test was used to determine the statistical significance between pre and post timings within treatment group, as denoted by *P value <0.05, **P value <0.01, ***P value <0.001. Error bars represent the mean ± standard error. 105 Figure 3.2: Methylation level of genes and their features. Methylation levels of all cytosines in various sequence contexts (CG, CHG, CHH) identified through bisulfite sequencing of kochia leaf tissue before (Pre) and three weeks after treatment (Post) with 0 106 Fig 3.2 (cont’d): (Untreated), 52.5 (Low) and 105 (High) g ae ha-1 of glyphosate for A) Genes plus 2kb upstream and downstream, B) only genes, C) mRNA, D) coding sequence, E) exons, and F) introns. Plants were grown in a greenhouse setting and included four replicates from each of the treatments both pre and post application timings. A two-tailed student’s T-Test was used to determine the statistical significance between pre and post timings within treatment group, as denoted by *P value <0.05, **P value <0.01, ***P value <0.001. Error bars represent the mean ± standard error. Note that axes are not scaled the same across plots. Figure 3.3: Methylation patterns of genes. Metaplots of A) CG B) CHG and C) CHH methylation levels averaged across all genes in the kochia genome, sequenced from leaf tissue before (Pre) and three weeks after treatment (Post) with 0 (Untreated), 52.5 (Low) and 105 (High) g ae ha-1 of glyphosate. Plants were grown in a greenhouse setting and included four replicates from each of the treatments both pre and post application timings. Note that y-axes are not scaled the same across plots. The x-axis indicates the transcription start site (TSS) and transcription termination site (TTS), as well as 2kb regions upstream and downstream from the gene body. Different colored lines are used for each treatment, with the solid line representing constitutive methylation of all cytosines prior to treatment, and dotted lines representing three weeks after treatment. 107 Figure 3.4: Volcano plot of DMRs in the CG context. Volcano plots representing all differentially methylated regions (DMRs) in the CG context between kochia methylomes generated from leaf tissue before and three weeks after treatment with A) 0 (Untreated), B) 52.5 (Low) and C) 105 (High) g ae ha-1 of glyphosate. Plants were grown in a greenhouse setting and included four replicates from each of the treatments both pre and post application timings. Each point represents a DMR, and significant (P value < 0.05 and |Log2FC > 0.5| based on a two-tailed student’s T-Test assuming heteroscedastic variance) differential methylation is represented by point color, where red indicates hypermethylation, blue indicates hypomethylation and gray indicates no significant changes. These plots represent 5,244 total DMRs in untreated, 3,561 in the low dose treatment of glyphosate, and 4,619 DMRs in the high dose. Note that axes are scaled the same across plots. 108 Figure 3.5: Venn diagrams of differentially methylated genes. Venn diagrams of differentially methylated genes (DMGs) that are significantly hypermethylated or hypomethylated between 109 Figure 3.5 (cont’d): kochia methylomes generated from leaf tissue before and three weeks after treatment with 0 (Untreated), 52.5 (Low) and 105 (High) g ae ha-1 of glyphosate. Plants were grown in a greenhouse setting. Overlaps show a shared treatment response. DMGs were defined as having one or more differentially methylated regions within 2kb upstream or downstream of the gene, with a P value < 0.05 and absolute Log2FC > 0.5 calculated with a student’s T-Test assuming heteroscedastic variance. A) hypermethylated CG, B) hypomethylated CG, C) hypermethylated CHG, D) hypomethylated CHG, E) hypermethylated CHH, and F) hypomethylated CHH. 110 Methylation Genes Annotation Log2FC Pvalue status Hyper Bs.00g110970 PREDICTED: uncharacterized protein LOC104890603 5.85 0.045 (Invertase/pectin methylesterase inhibitor domain superfamily) Bs.00g425320 PREDICTED: pentatricopeptide repeat-containing 5.80 0.031 protein At3g02330 Bs.00g252180 Unknown Protein (Region of a membrane-bound protein 5.37 0.013 predicted to be outside the membrane, in the extracellular region.) Bs.00g049890 Unknown Protein 4.81 0.034 Bs.00g049160 axial regulator YABBY 5 3.88 0.026 Bs.00g026340 hydroxyproline O-galactosyltransferase GALT2 3.85 0.047 Bs.00g359090 cytochrome P450 94A2-like 3.46 0.045 Bs.00g540350 PREDICTED: desiccation-related protein PCC13-62 3.39 0.007 Bs.00g194740 metacaspase-9-like 3.32 0.000 Bs.00g193680 uncharacterized protein LOC110792946 (Nucleic acid- 3.29 0.041 Bs.00g193690 binding, OB-fold), PREDICTED: dCTP pyrophosphatase Hypo Bs.00g360260 Unknown Protein, transcription initiation factor TFIID -5.96 0.0072 Bs.00g360270 subunit 4b-like Bs.00g517480 PREDICTED: lipid transfer-like protein VAS -5.62 0.0265 Bs.00g517490 potassium transporter 4 Bs.00g124580 PREDICTED: AP2-like ethylene-responsive -5.53 0.0462 transcription factor PLT2 Bs.00g039890 auxin efflux carrier component 2-like -5.13 0.0329 Bs.00g192850 protein ECERIFERUM 26-like -5.06 0.0101 Bs.00g226060 PREDICTED: putative B3 domain-containing protein -4.68 0.0350 Os03g0621600 Bs.00g118850 WAT1-related protein At1g09380 -4.62 0.0368 Table 3.2: High dose hypermethylated and hypomethylated CG DMRs. Top hyper and hypomethylated genes containing differentially methylated regions in the CG sequence context between kochia methylomes generated from leaf tissue before and three weeks after 111 Table 3.2 (cont’d): Methylation Genes Annotation Log2FC Pvalue status Hypo Bs.00g320630 protein trichome birefringence-like 42 -4.58 0.0441 Bs.00g320640 Unknown Protein (Ribonuclease H-like superfamily) Bs.00g284930 PREDICTED: ACT domain-containing -4.56 0.0002 protein ACR10 isoform X1 Bs.00g424610 polygalacturonase QRT2-like -4.48 0.0405 -1 treatment with 105 (High) g ae ha of glyphosate. Plants were grown in a greenhouse setting and included four replicates. DMRs are defined as having |Log2Fold Change > 0.5| and a P value < 0.05 calculated with a student’s T-Test assuming heteroscedastic variance. 112 Figure 3.6: Volcano plot of DMRs in the CHG context. Volcano plots representing all differentially methylated regions (DMRs) in the CHG context between kochia methylomes generated from leaf tissue before and three weeks after treatment with A) 0 (Untreated), B) 52.5 (Low) and C) 105 (High) g ae ha-1 of glyphosate. Plants were grown in a greenhouse setting and included four replicates from each of the treatments both pre and post application timings. Each point represents a DMR, and significant (P value < 0.05 and |Log2FC > 0.5| based on a two-tailed student’s T-Test assuming heteroscedastic variance) differential methylation is represented by point color, where red indicates hypermethylation, blue indicates hypomethylation and gray indicates no significant changes. These plots represent 2,911 DMRs in untreated, 2,193 in the low dose treatment of glyphosate, and 2,283 DMRs in the high dose. 113 Methylation GeneID Annotation Log2FC Pvalue status Hyper Bs.00g055630 PREDICTED: uncharacterized protein LOC104902397 4.86 0.046 Bs.00g139940 uncharacterized protein LOC110709754 isoform X1 (Checkpoint protein 4.55 0.004 Rad17/Rad24) Bs.00g441690 SAC3 family protein A isoform X7 4.31 0.045 Bs.00g352850 transcription factor TCP7-like 3.74 0.002 Bs.00g386890 PREDICTED: wall-associated receptor kinase-like 20 3.59 0.043 Bs.00g451600 Unknown Protein (Region of a membrane-bound protein predicted to be 3.00 0.039 embedded in the membrane) Bs.00g481010 heat shock protein 90-5, chloroplastic-like 2.93 0.026 Bs.00g481020 PREDICTED: glucomannan 4-beta-mannosyltransferase 9 Bs.00g171460 PREDICTED: uncharacterized protein LOC104892105 isoform X1 2.86 0.048 Bs.00g044380 thioredoxin O2, mitochondrial isoform X3 2.63 0.039 Bs.00g380040 Unknown Protein (MobiDBLite|mobidb-lite|consensus disorder prediction) 2.60 0.021 Hypo Bs.00g367780 psbP-like protein 1, chloroplastic -4.07 0.048 Bs.00g264960 E3 ubiquitin-protein ligase RFI2-like -3.81 0.040 Bs.00g122060 PREDICTED: uncharacterized protein LOC104907727 -3.74 0.012 Bs.00g196410 PREDICTED: uncharacterized protein LOC104889490 -3.51 0.005 Bs.00g172720 PREDICTED: calcium-transporting ATPase 3, endoplasmic reticulum-type -3.00 0.040 Bs.00g286270 Unknown Protein:(MobiDBLite|mobidb-lite|consensus disorder prediction) -2.94 0.025 Bs.00g286260 glutathione S-transferase T1-like -2.94 0.025 Bs.00g231980 PREDICTED: mitochondrial carrier protein MTM1 isoform X1 -2.83 0.033 Bs.00g498250 PREDICTED: probable potassium transporter 17 -2.58 0.030 Bs.00g356050 transcription factor-like protein DPB -2.51 0.006 Table 3.3: Top CHG DMGs in high treatment. Top hyper and hypomethylated genes containing differentially methylated regions in the CHG sequence context between kochia methylomes generated from leaf tissue before and three weeks after treatment with 105 (High) g ae ha-1 of glyphosate. Plants were grown in a greenhouse setting and included four replicates. DMRs are defined as having |Log2Fold Change > 0.5| and a P value < 0.05 calculated with a student’s T-Test assuming heteroscedastic variance. 114 Figure 3.7: Volcano plot of DMRs in the CHH context. Volcano plots representing all differentially methylated regions (DMRs) in the CHH context between kochia methylomes generated from leaf tissue before and three weeks after treatment with A) 0 (Untreated), B) 52.5 (Low) and C) 105 (High) g ae ha-1 of glyphosate. Plants were grown in a greenhouse setting and included four replicates from each of the treatments both pre and post application timings. Each point represents a DMR, and significant (P value < 0.05 and |Log2FC > 0.5| based on a two-tailed student’s T-Test assuming heteroscedastic variance) differential methylation is represented by point color, where red indicates hypermethylation, blue indicates hypomethylation and gray indicates no significant changes. These plots represent 12,870 total DMRs in untreated, 8,470 in the low dose treatment of glyphosate, and 8,979 DMRs in the high dose. 115 Methylation GeneID Annotation Log2FC P value status Hyper Bs.00g402840 PREDICTED: uncharacterized protein LOC104893067 5.69 0.032 Bs.00g338680 PREDICTED: ABC transporter A family member 7 5.62 0.015 Bs.00g370560 uncharacterized protein LOC110690478 (MobiDBLite|mobidb- 5.22 0.025 lite|consensus disorder prediction) Bs.00g196810 condensin complex subunit 1 5.02 0.040 Bs.00g301500 PREDICTED: tricyclene synthase TPS4, chloroplastic-like isoform X1 4.95 0.009 Bs.00g302140 REDICTED: (6-4)DNA photolyase isoform X1 4.94 0.020 Bs.00g054660 PREDICTED: uncharacterized protein LOC104900192 isoform X1 4.92 0.004 Bs.00g054670 Unknown Protein (Region of a membrane-bound protein predicted to be 4.92 0.004 outside the membrane, in the cytoplasm) Bs.00g452190 Unknown Protein 4.80 0.024 Bs.00g452200 isoleucine--tRNA ligase, chloroplastic/mitochondrial isoform X1 4.80 0.024 Bs.00g255710 Unknown Protein 4.80 0.011 Hypo Bs.00g226520 F-box protein At3g07870-like isoform X2 (Galactose oxidase/kelch, beta- -5.23 0.00002 propeller) Bs.00g258800 Unknown Protein (PANTHER|PTHR34538:SF4|EXPRESSED PROTEIN) -3.81 0.040 Bs.00g196410 PREDICTED: uncharacterized protein LOC104889490 (Sigma factor -3.12 0.008 binding protein 1/2) Bs.00g054940 zinc finger CCCH domain-containing protein 23-like -2.40 0.001 Bs.00g232460 uncharacterized protein LOC110696143 (Modifying wall lignin-1/2) -2.39 0.002 Bs.00g491000 PREDICTED: probable mannitol dehydrogenase -2.32 0.0003 Bs.00g189650 ethylene-responsive transcription factor RAP2-11-like -2.18 0.020 Bs.00g431710 protein NRT1/ PTR FAMILY 5.1-like (Proton-dependent oligopeptide -2.17 0.008 transporter family) Bs.00g225730 protein SHORT-ROOT-like -2.00 0.030 Bs.00g225740 probable transcriptional regulator SLK2 -2.00 0.030 Table 3.4: Top CHH DMGs in high treatment. Top hyper and hypomethylated genes containing differentially methylated regions in the CHH sequence context between kochia methylomes generated from leaf tissue before and three weeks after treatment with 105 116 Table 3.4 (cont’d): (High) g ae ha-1 of glyphosate. Plants were grown in a greenhouse setting and included four replicates. DMRs are defined as having |Log2Fold Change > 0.5| and a P value < 0.05 calculated with a student’s T-Test assuming heteroscedastic variance. Figure 3.8: Karyotypes. Karyotypes of significant differentially methylated regions (|Log2Fold Change > 0.5|, P value < 0.05 calculated with a student’s T-Test assuming heteroscedastic variance) in the CG, CHG and CHH contexts between kochia methylomes generated from leaf tissue before and three weeks after treatment 105 g ae ha-1 of glyphosate. Plants were grown in a greenhouse setting and included four replicates from each of the treatments both pre and post application timings. Methylation level of all four replicates was averaged and converted into binary values based on whether they were hypo- or hypermethylated after treatment. Hypomethylation is ‘Low’ on the scale while hypermethylation is ‘High’. 117 Figure 3.9: Correlation of methylation to gene expression. Scatter plots generated in Microsoft Excel showing the correlation of the Log2Fold Change in methylation of significant differentially methylated genes (DMGs) (genes containing one or more differentially 118 Fig 3.9 (cont’d): methylated regions with |Log2Fold Change > 0.5|, P value < 0.05 calculated with a student’s T-Test assuming heteroscedastic variance) in kochia leaf tissue three weeks after treatment with 105 (high) g ae ha-1 of glyphosate to the Log2Fold Change in expression of plants in the same population. Plants were grown in a greenhouse setting and included four replicates. All plots on the left show correlation of DMGs with all differentially expressed genes (DEGs) regardless of significance, and those on the right only show significant DEGs (P value < 0.05, |Log2FC in expression > 1.0|). Black points indicate a gene exclusively differentially methylated after high dose treatment, and purple are DMGs shared between high and low dose (52.5 g ae h-1 of glyphosate). Slope of the linear regression, R2 value, and P value from testing the null hypothesis that the line is zero using a t-test are all reported for each plot where *** is P value <0.01. 119 Figure 3.10: Weighted methylation of transposons. Methylation level of all cytosines in A) Copia, B) Gypsy, and C) MULE transposons by sequence contexts (CG, CHG, CHH) in kochia 120 Figure 3.10 (cont’d): methylomes generated from leaf tissue before (Pre) and three weeks after treatment (Post) with 0 (Untreated), 52.5 (Low) and 105 (High) g ae ha-1 of glyphosate. Plants were grown in a greenhouse setting. Significance between pre and post samples within treatments was evaluated with a two-tailed student’s T-Test assuming heteroscedastic variance is denoted as *P value <0.05, **P value <0.01, and ***P value <0.001. 121 TE Sequence # Of DMR # Of DMR Shared DMR % Of total family context UNTP vs LOWP UNTP vs HIGHP UNTP vs Glyphosate Glyphosate specific Copia CG 28 26 6 11.1 CHG 39 31 8 11.4 CHH 313 352 36* 5.4 Gypsy CG 46 25 6 8.5 CHG 75 40 7 6.1 CHH 556 595 60** 5.2 MULE CG 9 7 1 6.3 CHG 3 4 0 0 CHH 98 123 12 5.4 Table 3.5: DMRs with transposons. Number of differentially methylated regions (DMRs) compared between kochia methylomes generated from leaf tissue three weeks after treatment with 0 (Untreated), 52.5 (Low) and 105 (High) g ae ha-1 of glyphosate which intersect with Copia, Gypsy, or MULE transposons. Plants were grown in a greenhouse setting. Comparisons were made between post treatment samples, where UNTP = Untreated Post, LOWP = Low Dose Post, and HIGHP = High Dose Post. Shared DMRs between both glyphosate doses were also compared to UNT and reported as the proportion of “glyphosate specific” DMRs. In the Copia family, * denotes there are an additional ten transposons with DMRs shared between high and low dose treatment, but they included sites in multiple contexts. For Gypsy family, the ** denotes an additional eight transposons with multiple contexts as well. 122 REFERENCES 123 REFERENCES Anunciato, V. M., Bianchi, L., Gomes, G. L. G. C., Velini, E. D., Duke, S. O., and Carbonari, C. A. (2022). Effect of low glyphosate doses on flowering and seed germination of glyphosate-resistant and -susceptible Digitaria insularis. Pest Management Science, 78(3), 1227–1239. Ashapkin, V. v., Kutueva, L. I., Aleksandrushkina, N. I., and Vanyushin, B. F. (2020). Epigenetic Mechanisms of Plant Adaptation to Biotic and Abiotic Stresses. International Journal of Molecular Sciences, 21(20), 7457. Bewick, A. J., Niederhuth, C. E., Ji, L., Rohr, N. A., Griffin, P. T., Leebens-Mack, J., and Schmitz, R. J. (2017). The evolution of CHROMOMETHYLASES and gene body DNA methylation in plants. Genome Biology, 18(1), 65. Bilichak, A., Ilnystkyy, Y., Hollunder, J., and Kovalchuk, I. (2012). The Progeny of Arabidopsis thaliana Plants Exposed to Salt Exhibit Changes in DNA Methylation, Histone Modifications and Gene Expression. PLoS ONE, 7(1), e30515. Busi, R., and Powles, S. B. (2009). Evolution of glyphosate resistance in a Lolium rigidum population by glyphosate selection at sublethal doses. Heredity, 103(4), 318–325. Chen, S., Zhou, Y., Chen, Y., and Gu, J. (2018). fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics, 34(17), i884–i890.Cokus, S. J., Feng, S., Zhang, X., Chen, Z., Merriman, B., Haudenschild, C. D., Pradhan, S., Nelson, S. F., Pellegrini, M., and Jacobsen, S. E. (2008). Shotgun bisulphite sequencing of the Arabidopsis genome reveals DNA methylation patterning. Nature, 452(7184), 215–219. Gaines, T. A., Barker, A. L., Patterson, E. L., Westra, P., Westra, E. P., Wilson, R. G., Jha, P., Kumar, V., and Kniss, A. R. (2016). EPSPS Gene Copy Number and Whole-Plant Glyphosate Resistance Level in Kochia scoparia. PLOS ONE, 11(12), e0168295. Gutschker, S., Corral, J. M., Schmiedl, A., Ludewig, F., Koch, W., Fiedler-Wiechers, K., Czarnecki, O., Harms, K., Keller, I., Martins Rodrigues, C., Pommerrenig, B., Neuhaus, H. E., Zierer, W., Sonnewald, U., and Müdsam, C. (2022). Multi-omics data integration reveals link between epigenetic modifications and gene expression in sugar beet (Beta vulgaris subsp. vulgaris) in response to cold. BMC Genomics, 23(1), 144. Hao, Z., Lv, D., Ge, Y., Shi, J., Weijers, D., Yu, G., and Chen, J. (2020). RIdeogram: Drawing SVG graphics to visualize and map genome-wide data on the idiograms. PeerJ Computer Science, 6, 1–11. Hastings, P. J., Bull, H. J., Klump, J. R., and Rosenberg, S. M. (2000). Adaptive Amplification. Cell, 103(5), 723–731. 124 Hong, J. C. (2016). General Aspects of Plant Transcription Factor Families. Plant Transcription Factors: Evolutionary, Structural and Functional Aspects, 35–56. Ito, H., Gaubert, H., Bucher, E., Mirouze, M., Vaillant, I., and Paszkowski, J. (2011). An siRNA pathway prevents transgenerational retrotransposition in plants subjected to stress. Nature 472(7341), 115–119. Kawakatsu, T., Stuart, T., Valdes, M., Breakfield, N., Schmitz, R. J., Nery, J. R., Urich, M. A., Han, X., Lister, R., Benfey, P. N., & Ecker, J. R. (2016). Unique cell-type-specific patterns of DNA methylation in the root meristem. Nature Plants, 2(5), 1–8. Kim, G., Clarke, C. R., Larose, H., Tran, H. T., Haak, D. C., Zhang, L., Askew, S., Barney, J., and Westwood, J. H. (2017). Herbicide injury induces DNA methylome alterations in Arabidopsis. PeerJ, 5(7), e3560. Lämke, J., and Bäurle, I. (2017). Epigenetic and chromatin-based mechanisms in environmental stress adaptation and stress memory in plants. Genome Biology, 18(1), 124. Law, J. A., and Jacobsen, S. E. (2010). Establishing, maintaining and modifying DNA methylation patterns in plants and animals. Nature Reviews Genetics, 11(3), 204–220. Lloyd, J. P. B., & Lister, R. (2021). Epigenome plasticity in plants. Nature Reviews Genetics, 23(1), 55–68. Lu, Y. C., Feng, S. J., Zhang, J. J., Luo, F., Zhang, S., and Yang, H. (2016). Genome-wide identification of DNA methylation provides insights into the association of gene expression in rice exposed to pesticide atrazine. Scientific Reports, 6(1), 18985. Matzke, M. A., and Mosher, R. A. (2014). RNA-directed DNA methylation: an epigenetic pathway of increasing complexity. Nature Reviews Genetics, 15(6), 394–408. Microsoft Corporation. (2018). Microsoft Excel. Retrieved from https://office.microsoft.com/excel Niederhuth, C. E., Bewick, A. J., Ji, L., Alabady, M. S., Kim, K. do, Li, Q., Rohr, N. A., Rambani, A., Burke, J. M., Udall, J. A., Egesi, C., Schmutz, J., Grimwood, J., Jackson, S. A., Springer, N. M., and Schmitz, R. J. (2016). Widespread natural variation of DNA methylation within angiosperms. Genome Biology, 17(1), 194. Niederhuth, C. E., and Schmitz, R. J. (2017). Putting DNA methylation in context: from genomes to gene expression in plants. Biochimica et Biophysica Acta (BBA) - Gene Regulatory Mechanisms, 1860(1), 149–156. Patterson, E. L., Saski, C. A., Sloan, D. B., Tranel, P. J., Westra, P., and Gaines, T. A. (2019). The Draft Genome of Kochia scoparia and the Mechanism of Glyphosate Resistance via 125 Transposon-Mediated EPSPS Tandem Gene Duplication. Genome Biology and Evolution, 11(10), 2927–2940. Perrone, A., and Martinelli, F. (2020). Plant stress biology in epigenomic era. Plant Science, 294, 110376. Pettinga, D. J., Ou, J., Patterson, E. L., Jugulam, M., Westra, P., and Gaines, T. A. (2018). Increased chalcone synthase (CHS) expression is associated with dicamba resistance in Kochia scoparia. Pest Management Science, 74(10), 2306–2315. Preston, C., Belles, D. S., Westra, P. H., Nissen, S. J., and Ward, S. M. (2009). Inheritance of Resistance to The Auxinic Herbicide Dicamba in Kochia (Kochia scoparia). Weed Science, 57(1), 43–47. Quinlan, A. R., and Hall, I. M. (2010). BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics (Oxford, England), 26(6), 841–842. Reams, A. B., and Roth, J. R. (2015). Mechanisms of Gene Duplication and Amplification. Cold Spring Harbor Perspectives in Biology, 7(2), a016592. Ritz, C., Baty, F., Streibig, J. C., and Gerhard, D. (2015). Dose-Response Analysis Using R. PloS One, 10(12). Sasaki, E., Kawakatsu, T., Ecker, J. R., and Nordborg, M. (2019). Common alleles of CMT2 and NRPE1 are major determinants of CHH methylation variation in Arabidopsis thaliana. PLoS Genetics, 15(12). Schultz, M. D., He, Y., Whitaker, J. W., Hariharan, M., Mukamel, E. A., Leung, D., Rajagopal, N., Nery, J. R., Urich, M. A., Chen, H., Lin, S., Lin, Y., Jung, I., Schmitt, A. D., Selvaraj, S., Ren, B., Sejnowski, T. J., Wang, W., and Ecker, J. R. (2015). Human body epigenome maps reveal noncanonical DNA methylation variation. Nature, 523(7559), 212–216. Secco, D., Wang, C., Shou, H., Schultz, M. D., Chiarenza, S., Nussaume, L., Ecker, J. R., Whelan, J., and Lister, R. (2015). Stress induced gene expression drives transient DNA methylation changes at adjacent repetitive elements. ELife, 4. Smit, A. F. A., and Hubley, R. (2008). RepeatModeler Open-1.0. Steward, N., Ito, M., Yamaguchi, Y., Koizumi, N., and Sano, H. (2002). Periodic DNA Methylation in Maize Nucleosomes and Demethylation by Environmental Stress. Journal of Biological Chemistry, 277(40), 37741–37746. Tyczewska, A., Gracz-Bernaciak, J., Szymkowiak, J., and Twardowski, T. (2021). Herbicide stress-induced DNA methylation changes in two Zea mays inbred lines differing in Roundup® resistance. Journal of Applied Genetics, 62(2), 235–248. 126 Underwood, C. J., Henderson, I. R., and Martienssen, R. A. (2017). Genetic and epigenetic variation of transposable elements in Arabidopsis. Current Opinion in Plant Biology, 36, 135–141. VIB / UGent Bioinformatics & Evolutionary Genomics. (2022). Calculate and draw custom Venn diagrams. Wang, H., Beyene, G., Zhai, J., Feng, S., Fahlgren, N., Taylor, N. J., Bart, R., Carrington, J. C., Jacobsen, S. E., and Ausin, I. (2015). CG gene body DNA methylation changes and evolution of duplicated genes in cassava. Proceedings of the National Academy of Sciences, 112(44), 13729–13734. Weil, C., and Martienssen, R. (2008). Epigenetic interactions between transposons and genes: lessons from plants. Current Opinion in Genetics and Development, 18(2), 188–192. Wickham, H. 2016. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. Wiersma, A. T., Gaines, T. A., Preston, C., Hamilton, J. P., Giacomini, D., Robin Buell, C., Leach, J. E., and Westra, P. (2015). Gene amplification of 5-enol-pyruvylshikimate-3- phosphate synthase in glyphosate-resistant Kochia scoparia. Planta, 241(2), 463–474. Zhang, H., Lang, Z., and Zhu, J. K. (2018). Dynamics and function of DNA methylation in plants. Nature Reviews Molecular Cell Biology 2018 19:8, 19(8), 489–506. 127