Reduced coenzyme Q synthesis confers non-target site resistance to the herbicide thaxtomin A

Herbicide resistance in weeds is a growing threat to global crop production. Non-target site resistance is problematic because a single resistance allele can confer tolerance to many herbicides (cross resistance), and it is often a polygenic trait so it can be difficult to identify the molecular mechanisms involved. Most characterized molecular mechanisms of non-target site resistance are caused by gain-of-function mutations in genes from a few key gene families–the mechanisms of resistance caused by loss-of-function mutations remain unclear. In this study, we first show that the mechanism of non-target site resistance to the herbicide thaxtomin A conferred by loss-of-function of the gene PAM16 is conserved in Marchantia polymorpha, validating its use as a model species with which to study non-target site resistance. To identify mechanisms of non-target site resistance caused by loss-of-function mutations, we generated 107 UV-B mutagenized M. polymorpha spores and screened for resistance to the herbicide thaxtomin A. We isolated 13 thaxtomin A-resistant mutants and found that 3 mutants carried candidate resistance-conferring SNPs in the MpRTN4IP1L gene. Mprtn4ip1l mutants are defective in coenzyme Q biosynthesis and accumulate higher levels of reactive oxygen species (ROS) than wild-type plants. Mutants are weakly resistant to thaxtomin A and cross resistant to isoxaben, suggesting that loss of MpRTN4IP1L function confers non-target site resistance. Mutants are also defective in thaxtomin A metabolism. We conclude that loss of MpRTN4IP1L function is a novel mechanism of non-target site herbicide resistance and propose that other mutations that increase ROS levels or decrease thaxtomin A metabolism could contribute to thaxtomin A resistance in the field.


Introduction
The control of weeds by herbicides is threatened by the evolution of resistance in the field [1,2]. This is analogous to how human and animal health is threatened by the evolution of antibiotic resistance in bacteria. The intense selection pressure resulting from herbicide treatments can lead to the rapid evolution of herbicide resistance in the agricultural landscape [1]. Resistance can result from genetic changes in the gene encoding the herbicide target which lead to conformational changes in the target protein blocking herbicide binding, or overexpression of the target [3,4]. This form of resistance, which involves mutations only in the gene encoding the target protein, is known as target-site resistance. Plants with a target-site resistance mutation can be very strongly resistant to herbicides [1]. For example, target-site resistant biotypes of the weed Euphorbia heterophylla survive more than 50 times the recommended field dose of the acetohydroxyacid synthase (AHAS) inhibitor imazamox [5]. Resistance can also result from genetic changes in a variety of genes which prevent inhibitory levels of herbicide reaching the target, or that alleviate the toxic effects of the herbicide downstream from the target [1,3,6]. For example, a mutation which causes constitutive expression of a glutathione transferase gene confers resistance in the weed species Alopecurus myosuroides. In this example, resistance is the result of increased herbicide metabolism and increased accumulation of flavonol metabolites that act as antioxidants. Accumulation of these molecules can confer resistance to herbicides whose toxicity depends on the accumulation of toxic reactive oxygen species [7][8][9]. These types of resistance-which do not involve mutations in the gene encoding the herbicide target-are known as non-target site resistance. Individual alleles that contribute to non-target site resistance usually confer weaker resistance to herbicides than alleles responsible for target-site resistance [10].
Both target-site and non-target site resistance can be caused by gain-of-function mutations, which cause an increase in expression of a gene or enhanced activity of a protein. For example, the 5-enolpyruvyl shikiminate-3-phosphate synthase (EPSPS) gene encodes a protein essential for aromatic amino acid biosynthesis which is the target of the herbicide glyphosate [11]. Gain-of-function mutations in the gene encoding EPSPS resulting in gene amplification (multiple rounds of gene duplication) confer target-site resistance to the herbicide glyphosate [12,13]. EPSPS gene copy number is greater in many glyphosate-resistant populations than in glyphosate-sensitive populations, resulting in higher steady state levels of EPSPS mRNA, protein, and total enzyme activity. The resulting elevated levels of EPSPS activity increase glyphosate resistance [13]. Non-target site resistance can also be conferred by spontaneous gain-offunction mutations in genes encoding enzyme activities that chemically modify or compartmentalize herbicides and prevent the herbicide reaching the target [1,6,14]. For example, constitutive overexpression of an ABC transporter causing glyphosate export into the apoplast results in glyphosate resistance in E. colona [15].
Loss-of-function mutations-mutations that result in reduced gene and protein functioncan also cause target-site resistance. This is because mutations in the gene encoding a herbicide target which cause conformational changes that prevent herbicide binding can also reduce the capacity of the protein to carry out its function. For example, some mutations in EPSPS that confer glyphosate resistance by preventing the binding of glyphosate also reduce the affinity of EPSPS for its substrate phosphoenolpyruvate, leading to reduced catalytic activity [16]. However, very few loss-of-function mutations causing non-target site resistance have been reported in weeds. Recent genome wide association studies have identified mutations associated with glyphosate resistance in the weed species Amaranthus tuberculatus and Ipomoea purpurea of which some are likely to be loss-of-function, suggesting that the contribution of loss-of-function mutations to non-target site resistance has been underestimated, however the contribution of these mutations to resistance has not yet been experimentally confirmed [17][18][19]. Only two mechanisms of non-target site resistance involving loss-of-function mutations have been proven to confer resistance, using forward genetics in the model species Arabidopsis thaliana [20,21]. As such, the role of loss-of-function mutations in non-target site resistance in weeds is unclear.
It is likely that both gain-and loss-of-function mutations contribute to non-target site resistance in herbicide resistant weed populations, but it is difficult to identify resistance mechanisms of non-target site resistance, including those conferred by loss of gene function. This is partly because non-target site resistance is usually a polygenic trait, and each single allele may contribute little to resistance. For example, in a sample of glyphosate resistant A. tuberculatus individuals, it was estimated that 25% of the variance in resistance was attributed to non-target site resistance associated with SNPs in 274 genes with a range of allele frequencies and effect sizes [17,19]. Identifying the mutant genes that contribute such a small amount to total resistance phenotypes is challenging in weeds. This makes defining the molecular basis of non-target site resistance difficult. The few alleles which have been identified are those which have a large effect on resistance. So far, these have all been gain-of-function alleles involving genes from a few key gene families such as cytochrome P450s, glutathione-S-transferases, glycosyltransferases, and ABC transporters [6]. There are likely a variety of undiscovered non-target site resistance causing mutations which each have a small effect on resistance but additively cause strong resistance. Identifying these diverse mechanisms of non-target site resistancewhich likely include loss-of-function mutations-remains a challenge.
We set out to identify mechanisms of non-target site resistance to thaxtomin A-a natural product being developed as a herbicide-caused by loss-of-function mutations. We chose thaxtomin A because it has been approved by the US Environmental Protection Agency as a herbicide (EPA registration number 84059-12) [22][23][24][25][26][27][28]. Thaxtomin A is thought to act by inhibiting plant cell wall biosynthesis, a process necessary for plant growth and the inhibition of which leads to plant death [29,30]. This plant specific mode of action makes thaxtomin A a favourable herbicide choice due to its low toxicity to animals. To date, only one mechanism of resistance to thaxtomin A has been identified (via forward mutagenesis in Arabidopsis thaliana) caused by loss-of-function mutations in PAM16, which encodes a member of the TIM complex, involved in protein transport into mitochondria [20,31,32].
To identify novel loss-of-function mutations that confer resistance to thaxtomin A, we UV mutagenized a large population of Marchantia polymorpha spores and screened for mutants which survived the lethal dose of thaxtomin A. We chose M. polymorpha as a model species because a single female sporophyte produces millions of haploid spores which can be mutagenized and directly screened for resistance. This enables rapid mutant screens in large populations that would be difficult to carry out at the same scale in diploid flowering plant models [33,34]. Furthermore, there is evidence that gene networks underlying bryophyte gametophyte physiology and development are similar to those in angiosperms [35,36]. We identified 13 thaxtomin A-resistant mutants. Three independent mutants carried candidate resistance-conferring mutations in the MpRTN4IP1L gene, which functions in the biosynthesis of coenzyme Q, a component of the mitochondrial electron transport chain. Mprtn4ip1l mutants accumulated higher levels of reactive oxygen species (ROS) than wild type controls and were defective in thaxtomin A metabolism. Taken together these data suggest that mutations that lead to ROS accumulation or which hinder thaxtomin A metabolism may confer resistance to thaxtomin A. We predict that loss-of-function mutations that increase ROS levels or decrease thaxtomin A metabolism in weeds will confer resistance to thaxtomin A in the field.

A mechanism of non-target site resistance conferred by loss-of-function of PAM16 is conserved in M. polymorpha
To determine if M. polymorpha is a suitable model species for discovering novel mechanisms of non-target site resistance to thaxtomin A, we first tested if M. polymorpha is sensitive to thaxtomin A. Wild type gemmae from Takaragaike-1 (Tak-1) and Takaragaike-2 (Tak-2) accessions were grown on solid medium supplemented with different concentrations of thaxtomin A. The area of living tissue of these plants was measured after 14 days (Fig 1A and Tab A in S1 Table). Thaxtomin A treatment inhibits the growth of wild type M. polymorpha in a dose-dependent manner (Fig 1A and 1B). The IC 50 (dose at which average area of living tissue is reduced by 50%) was 56 ± 6 nM for Tak-1, and 135 ± 24 nM for Tak-2. The LD 100 (lethal dose at which no plant tissue is alive) for both Tak-1 and Tak-2 was 5 μM, although in rare cases Tak-1 can survive this dose (Fig 1A and 1B).
To validate the suitability of M. polymorpha to study non-target site resistance to thaxtomin A, we tested whether a known mechanism of non-target site thaxtomin A resistance discovered in Arabidopsis thaliana is conserved in M. polymorpha. Loss-of-function mutations in the AtPAM16 gene (AT3G59280) confer non-target site resistance to thaxtomin A in A. thaliana [20,32]. PAM16 encodes a subunit of the PAM/TIM complex that transports proteins across the inner mitochondrial membrane [31]. We tested the hypothesis that loss of PAM16 function confers thaxtomin A resistance in M. polymorpha.
The closest homologue of AtPAM16 in M. polymorpha was identified by BlastP as Mp3g09390.1 [37]. To determine the evolutionary history of the PAM16 gene and determine if Mp3g09390.1 is the only copy of PAM16 in M. polymorpha, a phylogeny was created of the homologues of PAM16 from 11 different species (Figs 1C and S1A and S1B and S2 Table and S1 Text). There were between one and three copies of PAM16 in each species and there was one copy in M. polymorpha (Mp3g09390.1). Mp3g09390.1 is therefore the only PAM16 gene in the M. polymorpha genome and will be referred to as MpPAM16.
To determine if loss of MpPAM16 function confers thaxtomin A resistance, Mppam16 lossof-function mutants were generated using CRISPR-Cas9 mutagenesis. Guide RNAs were designed to generate mutations in the highly conserved N-terminal predicted signal peptide which targets the PAM16 protein to the inner mitochondrial membrane (sgRNA 2), and a non-conserved region at the C-terminal end of the protein (sgRNA 1) (Fig 1D). Thirty-three Mppam16 mutants were generated. The highly conserved signal peptide was mutated in 28 Mppam16 mutants; the non-conserved C-terminal region of the protein was mutated in the 5 remaining Mppam16 mutants.
Gemmae from the Mppam16 mutant lines were grown on solid medium supplemented with 5 μM thaxtomin A (LD 100 ) (Fig 1E) or 0.1% DMSO (Fig 1F). The area of living plant tissue was quantified after 12 days of growth (Tabs B and C in S1 Table). Mppam16 mutant lines which were significantly larger than wild type plants on thaxtomin A were classed as resistant. Nine Mppam16 lines were significantly larger than both Tak-1 and Tak-2 on 5 μM thaxtomin A and therefore classed as thaxtomin A-resistant (p < 0.05) (Fig 1E).
Five of the nine thaxtomin A-resistant Mppam16 mutants were significantly smaller than Tak-1 and Tak-2 plants in control conditions, including the four Mppam16 mutants which grew largest on thaxtomin A-Mppam16 GE258 , Mppam16 GE255 , Mppam16 GE235 , and Mppam16 GE26 (Fig 1F). This suggests that growth defects are more likely in Mppam16 mutant lines with stronger resistance to thaxtomin A.
All 9 thaxtomin A-resistant Mppam16 mutants were mutated in the highly conserved region that comprises the predicted PAM16 signal peptide (S1C Fig). The mutations in the thaxtomin A-resistant mutants affect between 2 and 19 amino acids and involve a range of substitutions, insertions, and deletions (S1C Fig). A mutation in the signal peptide is predicted to disrupt PAM16 function by preventing its transport to the mitochondria. These mutations are therefore likely to be loss-of-function. Given that all 9 thaxtomin A-resistant Mppam16 mutant lines carry likely loss-of-function mutations, we conclude that loss-of-function of Mppam16 can confer thaxtomin A resistance. The resistance of Mppam16 mutants to thaxtomin A validates M. polymorpha as a model to identify novel mechanisms of non-target site resistance to this herbicide.

thaxtomin A-resistant M. polymorpha mutants were isolated from a screen of 10 7 mutants
A genetic screen was carried out to identify mutants resistant to the herbicide thaxtomin A. To generate thaxtomin A-resistant mutants, 2 x 10 7 wild type M. polymorpha spores were plated on solid medium supplemented with 5 μM thaxtomin A and exposed to UV-B irradiation to induce mutations. A screening dose of 5 μM thaxtomin A (LD 100 ) was selected as screening at the lethal dose (as opposed to higher than the lethal dose) preferentially selects for non-target site resistance [39]. The dose of UV-B radiation used killed approximately 50% of spores (S2 Fig and Tab D in S1 Table). Ten million mutagenized spores were screened for resistance to thaxtomin A 14 days after mutagenesis (Fig 2A). Resistant mutants were transferred to fresh solid medium supplemented with 5 μM thaxtomin A to confirm their resistance. Mutants that survived the second round of selection were retained, and gemmae from each resistant mutant were plated onto fresh solid medium supplemented with 5 μM thaxtomin A to verify the retention of resistance through asexual reproduction (Fig 2B). Mutant gemmae that survived and grew significantly larger than wild type gemmae on 5 μM thaxtomin A were classed as thaxtomin A-resistant ( Fig 2C and Tab E in S1 Table). Using this protocol, 13 thaxtomin A-resistant A mechanism of non-target site resistance conferred by loss-of-function of PAM16 is conserved in M. polymorpha. A: Dose-response curves of the growth of wild type lines on thaxtomin A (TXTA). Wild type gemmae (Tak-1 and Tak-2) were grown for 14 days on solid medium supplemented with different concentrations of thaxtomin A (Tab A in S1 Table). The fitted curves and IC 50 values were calculated using the four-parameter log-logistic equation. Error bars represent ± standard deviation (n = 18). B: Wild type (Tak-1) gemmalings grown for 12 days on solid medium supplemented with DMSO, 100 nM thaxtomin A, or 5000 nM thaxtomin A. Images were taken with a Keyence VHX-7000. Scale bars represent 2 mm and were added in Inkscape v1.0.1. C: Phylogenetic analysis of PAM16 in plants. Proteins similar to AtPAM16 (At3G59280) were identified by protein BLAST search (E value < 1E-5) against the reference proteomes of various species (S2 Table and S1 Text) [37]. Orthologues were aligned and a maximum likelihood analysis was conducted. The tree was rooted with the PAM16 homologue from S. cerevisiae. A nonparametric approximate likelihood ratio test based on a Shimodaira-Hasegawa-like procedure was used to calculate branch support values [38]. D: Schematic representation of the MpPAM16 (Mp3g09390.1) gene constructed using Inkscape v1.0.1. CDS regions are represented in yellow. The regions of the gene targeted by guide RNAs are marked with green arrowheads. E and F: Wild type gemmae (Tak-1 and Tak-2) and gemmae from the nine thaxtomin A-resistant Mppam16 mutants were grown for 12 days on solid medium supplemented with 5 μM thaxtomin A (E) (Tab B in S1 Table) or 0.1% DMSO (F) (Tab C in S1 Table). The fitted curves and IC 50 values were calculated using the four-parameter log-logistic equation. Error bars represent ± standard deviation (n = 2-7). Stars represent the level of significance (as determined by Student's t-tests) of the difference between mutant and control lines (comparison to Tak-1 in red and Tak-2 in green): � = p < 0.05, �� = p < 0.01, ��� = p < 0.001. Mpthar (thaxtomin A resistant) mutant lines were identified out of the 10 7 mutagenized spores screened.

The MpRTN4IP1L gene is mutated in three Mpthar mutants
To identify mutations that confer thaxtomin A resistance, the genome of each Mpthar mutant was sequenced and the most likely resistance-conferring SNP in each mutant was identified using a modified version of a SNP calling bioinformatic pipeline adapted from [40]. Mismatches in the genomes of Mpthar mutants were compared to mismatches in the genomes of thaxtomin A-sensitive plants (Tak-1, Tak-2, and thaxtomin A-sensitive UV-B mutants [40]). Mismatches found only in Mpthar mutants and which passed further filtering steps were considered to be candidate resistance-conferring mutations. Between 6 and 27 candidate resistance-conferring mutations were identified for each Mpthar mutant. These SNPs were compared between Mpthar mutants. There were candidate resistance-conferring SNPs in the Mp3g19030.1 gene in 3 independent Mpthar mutants (Mpthar2, Mpthar4, and Mpthar6) ( Fig  3A and 3B). The presence of these SNPs was confirmed by Sanger sequencing (S3A, S3B and S3C Fig). Given that the Mp3g19030.1 gene carries a candidate resistance-conferring mutation in 3 independent Mpthar mutants, it is likely that a mutation in Mp3g19030.1 confers thaxtomin A resistance.
Mp3g19030.1 is annotated in the reference M. polymorpha genome (MpTak v6.1) as a homologue of the H. sapiens RETICULON 4-INTERACTING PROTEIN 1.1 (HsRTN4IP1.1) gene. HsRTN4IP1.1 encodes a mitochondrial NADH oxidoreductase involved in the biosynthesis of coenzyme Q, and its homologue RAD-8 was first discovered in C. elegans [41,42]. To determine the phylogenetic relationship between Mp3g19030.1 and its homologues, gene trees were generated from the homologs of HsRTN4IP1 from 10 species (Figs 3C and S3D and Table 1 and S1 Text). There were between 4 and 8 HsRTN4IP1.1 homologues in each of the species analyzed ( Table 1). The homologues group into a clade of 2-methylene-furan-3-one reductases (support value 1), a clade containing HsRTN4IP1.1 and its closest homologues (support value 0.9151), and a clade containing quinone oxidoreductases, alcohol dehydrogenases, and prostaglandin reductases (support value 0.8706) (Fig 3C). Based on the annotations, the identified homologues of HsRTN4IP1.1 therefore likely include enzymes with similar activities.
Mp3g19030.1 is the most closely related M. polymorpha homologue of HsRTN4IP1.1 gene based on BlastP. However, the topology of the phylogenetic tree suggests that Mp1g20380.1, a 13 thaxtomin A-resistant mutants were isolated from a screen of 10 7 mutants. A: Wild type spores from a cross between Tak-1 and Tak-2 grown for 14 days on solid medium supplemented with 0.1% DMSO or 5 μM thaxtomin A (TXTA), and UV-B mutagenesis selection plates supplemented with 5 μM thaxtomin A showing 14-day old thaxtomin A-resistant mutants Mpthar1 and Mpthar2 surrounded by dead mutagenized thaxtomin A-sensitive spores. Images were taken with a Leica DFC310 FX stereomicroscope. Scale bars represent 2 mm and were added in Inkscape v1.0.1. B: Wild type (Tak-1 and Tak-2) and Mpthar gemmalings grown for 12 days on solid medium supplemented with 0.1% DMSO or 5 μM thaxtomin A. Images were taken with a Keyence VHX-7000 and the maximum pixel value was adjusted to 188 in ImageJ. Scale bars represent 2 mm and were added in Inkscape v1.0.1. C: Gemmae from wild type lines (Tak-1 and Tak-2) and from Mpthar lines were plated on solid medium supplemented with DMSO (blue bars) or 5 μM thaxtomin A (orange bars) and grown for 21 days (n = 3-13) (Tab E in S1 Table). Error bars represent ± standard deviation. Stars represent the level of significance (as determined by Student's t-tests) of the difference between mutant and control lines subjected to the same treatment (comparison to Tak-1 in red and Tak-2 in green), or between individuals of the same genotype subjected to different treatments (black): � = p < 0.05, �� = p < 0.01, ��� = p < 0.001.

Loss-of-function of MpRTN4IP1L is a novel mechanism of herbicide resistance
To independently verify that loss of MpRTN4IP1L function confers thaxtomin A resistance, Mprtn4ip1l loss-of-function mutants were generated using CRISPR-Cas9 targeted mutagenesis. A guide RNA was designed to mutate the beginning of the gene to induce frameshifts (sgRNA 1), and 2 guide RNAs were designed to mutate the sites mutated in Mpthar2, Mpthar4, and Mpthar6 (sgRNA 2 and sgRNA 3) (Fig 3A). In total, 19 Mprtn4ip1l mutants were generated. Based on the nucleotide and predicted protein sequences, there were in frame indels in Mprtn4ip1l which affected 5 or fewer amino acids in the Mprtn4ip1l protein of Mprtn4ip1l Gemmae from the Mprtn4ip1l mutant lines and from lines with a wild type MpRTN4IP1L (Tak-1, Tak-2, and 3 lines which were transformed with the CRISPR construct but had no mutations in the MpRTN4IP1L gene; MpRTN4IP1L GE153 , MpRTN4IP1L GE314 , and MpRTN4IP1L GE356 ) were grown on solid medium supplemented with 5 μM thaxtomin A Proteins similar to the human RTN4IP1.1 protein were identified by protein BLAST search (E value < 1E-5) against the reference proteomes of various species (S1 Text) [37]. Orthologues were aligned and a maximum likelihood analysis was conducted. The tree was rooted with the RTN4IP1 homologue (Yim1p) from S. cerevisiae. A non-parametric approximate likelihood ratio test based on a Shimodaira-Hasegawa-like procedure was used to calculate branch support values [38]. D: Magnification of the red boxed branch from (C) showing that Mp3g19030.1 diverges into a monophyletic clade containing the human RTN4IP1.1 gene.

Picea abies Congenie 6
Brachypodium distachyon Ensembl 5 Azolla filiculoides Fernbase 7 Proteins similar to HsRTN4IP1.1 were identified via BlastP search of proteomes from these species [37]. Only homologues with an E value < 1 E-5 were included in the analysis.
https://doi.org/10.1371/journal.pgen.1010423.t001 (LD 100 ) ( Fig 4A and 4B and Tab F in S1 Table) or 0.1% DMSO (Fig 4A and 4C and Tab G in S1 Table). The average area of living plant tissue was quantified after 12 days of growth. Plant lines which grew significantly larger than Tak-1 and Tak-2 plants on thaxtomin A were classed as resistant. The 17 predicted loss-of-function Mprtn4ip1l lines were significantly larger than both Tak-1 and Tak-2 on 5 μM thaxtomin A so were classed as thaxtomin A-resistant (p<0.05) (Fig 4B). The 2 mutants with small in frame indels in Mprtn4ip1l and the lines with a wild type copy of MpRTN4IP1L however did not grow significantly larger than Tak-1 and Tak-2 on 5 μM thaxtomin A so were classed as thaxtomin A-sensitive ( Fig 4B). This suggests that loss-of-function of MpRTN4IP1L confers resistance to thaxtomin A.
When grown in control conditions, the areas of plant lines with a wild type copy of MpRTN4IP1L and the 2 Mprtn4ip1l mutants with small in frame indels are statistically indistinguishable from wild type ( Fig 4C). However, 13 of the 17 predicted loss-of-function Mprtn4ip1l mutants were significantly smaller than wild type lines in control conditions ( Fig  4C). Loss-of-function of MpRTN4IP1L therefore likely causes decreased growth in control conditions.
To quantify the thaxtomin A resistance of Mprtn4ip1l CRISPR mutants, gemmae from 2 independent Mprtn4ip1l CRISPR lines (Mprtn4ip1l GE118 and Mprtn4ip1l GE149 ) were grown on solid medium supplemented with different doses of thaxtomin A ( Fig 4D). The average area of living plant tissue was quantified after 12 days of growth (Tab H in S1 Table). Dose-response curves were fitted using the four-parameter log-logistic equation (Fig 4D). Three parameters of resistance were calculated from the curves: IC 50 (the concentration of thaxtomin A at which the area of living tissue is reduced by 50%), LD 100 (lethal dose at which there is no surviving tissue), and RI (resistance index; ratio between wild type and mutant IC 50 ). The IC 50 of Tak-2 was used to calculate the most conservative estimate of RI of each Mprtn4ip1l line. The IC 50 , RI, and LD 100 values of both Mprtn4ip1l mutant lines were significantly higher than either wild type line ( Fig 4D). Based on their RI, Mprtn4ip1l GE118 and Mprtn4ip1l GE149 are respectively at least 9.1 ± 2.95 times and 4.1 ± 2.29 times more resistant to thaxtomin A than wild type ( Fig 4D). Taken together, these data indicate that loss-of-function mutations in MpRTN4IP1L confer thaxtomin A resistance.

Mprtn4ip1l mutants produce lower levels of coenzyme Q and accumulate higher levels of reactive oxygen species than wild type plants
Having shown that loss-of-function mutations in the MpRTN4IP1L gene confer thaxtomin A resistance, we sought to determine how these mutations changed the physiology of plants to make them thaxtomin A resistant. HsRTN4IP1 interacts with proteins involved in coenzyme Q (CoQ) synthesis and knock-out mutants in the mouse RTN4IP1 gene accumulate less CoQ than wild type [42]. We therefore hypothesized that CoQ levels would be lower in Mprtn4ip1l mutants than in wild type M. polymorpha. To test this hypothesis, we measured the levels of coenzyme Q10 (CoQ) in the extracts of wild type (Tak-1) and Mprtn4ip1l GE149 mutants via LC-MS/MS. The levels of CoQ in Mprtn4ip1l GE149 are significantly lower than wild type levels in both control and TXTA-treated conditions, consistent with the hypothesis that MpRTN4IP1L is involved in CoQ synthesis (Fig 5A and Tab I in S1 Table).
Coenzyme Q is an electron acceptor that transports electrons from Complexes I and II to Complex III of the oxidative phosphorylation electron transport chain [46]. As electrons flow through the electron transport chain, electrons may leak from Complexes I, II and III and react with oxygen to form reactive oxygen species (ROS). The binding site of CoQ at Complex I is one of the sites of ROS production in the mitochondria [46]. Reduced levels of CoQ result in fewer electron carriers active in the oxidative phosphorylation transport chain leading to an excess of electrons leaking from the chain and reacting with oxygen to form ROS. Furthermore, non-mitochondrial CoQ acts as a lipid soluble antioxidant, protecting against oxidative stress by reacting with lipid peroxide radicals thereby preventing their ability to cause lipid peroxidation [47][48][49]. Therefore, we hypothesized that ROS levels would be higher in Mprtn4ip1l mutants than in wild type plants.
To test the hypothesis that Mprtn4ip1l mutants produce higher levels of ROS than wild type plants, 12-day old gemmalings of Tak-1, Tak-2, Mprtn4ip1l GE118 and Mprtn4ip1l GE149 were stained with 3,3'-diaminobenzidine (DAB). DAB reacts with hydrogen peroxide (H 2 O 2 )-a form of ROS-to form a brown precipitate. The amount of brown precipitate is proportional to the amount of H 2 O 2 in each sample [50]. The two Mprtn4ip1l mutants stained darker than the two wild type lines, indicating higher levels of H 2 O 2 in mutant thalli than in wild type ( Fig 5B). The concentration of various forms of ROS (including H 2 O 2 and lipid hydroperoxides) was also quantified in Mprtn4ip1l mutants using a modified ferric-xylenol orange (FOX) assay [51,52]. ROS concentrations were significantly higher in Mprtn4ip1l GE118 and Mprtn4ip1l GE149 than in wild type plants (Fig 5C and Tab J in S1 Table).
Since Mprtn4ip1l mutants produce more ROS than wild type, we hypothesized that Mprtn4ip1l mutants would be more sensitive (hypersensitive) to ROS treatment than wild type plants. To test this hypothesis the sensitivity of Mprtn4ip1l mutants to toxic levels of exogenous ROS was compared to that of wild type plants. Furthermore, a number of herbicides including paraquat, protoporphyrinogen oxidase inhibitors and glufosinate control weeds, at least in part, by causing ROS overproduction [53][54][55][56]. Mprtn4ip1l mutants produce more ROS than wild type, therefore we hypothesized that Mprtn4ip1l mutants would be hypersensitive to herbicides that cause ROS overproduction.
To test the hypothesis that Mprtn4ip1l mutants are more sensitive to ROS than wild type, the sensitivity of Mprtn4ip1l mutants and wild type was compared. Tak-1, Tak-2, Mprtn4ip1l GE118 and Mprtn4ip1l GE149 gemmae were grown on solid medium supplemented with a range of concentrations of H 2 O 2 or paraquat. The average area of living plant tissue was quantified after 10 days of growth (Tabs K and L in S1 Table). Dose-response curves were fitted using the four-parameter log-logistic equation, and resistance parameters IC 50 , LD 100 and RI were calculated from the curves (Fig 5D and 5E). Both Mprtn4ip1l mutants were more sensitive to H 2 O 2 than wild type lines. The RI values of Mprtn4ip1l GE118 and Mprtn4ip1l GE149 were 0.64 ± 0.26 and 0.88 ± 0.35 respectively, demonstrating that these lines are more sensitive to H 2 O 2 than either wild type line. The RI value of Mprtn4ip1l GE118 was statistically significant (Fig 5D). While the RI value of Mprtn4ip1l GE1549 was not statistically significant, the LD 100 of Mprtn4ip1l GE149 was lower than either wild type line ( Fig 5D). These data demonstrate that Mprtn4ip1l mutants are more sensitive to H 2 O 2 than wild type plants. Both Mprtn4ip1l mutants were also more sensitive to paraquat than wild type plants. The LD 100 of  (Tak-1 and Tak-2) and from Mprtn4ip1l mutants (Mprtn4ip1l GE118 and Mprtn4ip1l GE149 ) were grown for 12 days on solid medium supplemented with DMSO or 5 μM thaxtomin A (TXTA). Images were taken with a Keyence VHX-7000 and the maximum pixel value was adjusted to 200 in ImageJ. Scale bars represent 2 mm and were added in Inkscape v1.0.1. B and C: Wild type gemmae (Tak-1 and Tak-2) and gemmae from Mprtn4ip1l mutants were grown for 12 days on solid medium supplemented with 5 μM thaxtomin A (B) (Tab F in S1 Table) or 0.1% DMSO (C) (Tab G in S1 Table). Error bars represent ± standard deviation (n = 2-6). Stars represent the level of significance (as determined by Student's ttests) of the difference between mutant and control lines (comparison to Tak-1 in red and Tak-2 in green): � = p < 0.05, �� = p < 0.01, ��� = p < 0.001. Red arrowheads indicate lines with a wild type copy of MpRTN4IP1L; yellow arrowheads indicate lines with a small (< 5 amino acids) in frame indel in MpRTN4IP1L. Arrowheads were added in Inkscape v1.0.1. D: Dose-response curves of the growth of wild type and Mprtn4ip1l mutants on thaxtomin A, and table of corresponding IC 50 , LD 100 , and RI values. Gemmae from wild type lines (Tak-1 and Tak-2) and from two independent Mprtn4ip1l mutants generated via CRISPR-Cas9 mutagenesis (Mprtn4ip1l GE118 and Mprtn4ip1l GE149 ) were grown for 12 days on solid medium supplemented with different concentrations of thaxtomin A (Tab H in S1 Table). The fitted curves and IC 50 values were calculated using the four-parameter log-logistic equation. Error bars represent ± standard deviation (n = 3). https://doi.org/10.1371/journal.pgen.1010423.g004 Mprtn4ip1l GE118 and Mprtn4ip1l GE149 were 3.3 μM and 5 μM respectively, which are both lower than 7.5 μM, the LD 100 of Tak-1 and Tak-2 wild types ( Fig 5E). Neither of the RI values of Mprtn4ip1l GE118 and Mprtn4ip1l GE149 were statistically significant. Based on their LD 100 values, Mprtn4ip1l mutants are more sensitive to paraquat than wild type plants. Since paraquat kills plants by increasing ROS to toxic levels, these data demonstrate that Mprtn4ip1l mutants are more sensitive to ROS than wild type plants.
Together, these data indicate that ROS levels are higher in Mprtn4ip1l mutants than in wild type plants. Addition of H 2 O 2 to A. thaliana cells in culture confers partial thaxtomin A resistance [57]. The stiffening of the cell wall observed upon H 2 O 2 treatment is thought to prevent the toxic action of thaxtomin A on cell wall biosynthesis [57]. We hypothesize that the reduced levels of CoQ in Mprtn4ip1l mutants causes thaxtomin A resistance by increasing ROS levels which changes the physiology of the plant cell walls making thaxtomin A less effective. According to this hypothesis of ROS-induced cell wall alteration, loss-of-function of Mprtn4ip1l is a mechanism of non-target site resistance.

Loss-of-function of MpRTN4IP1L likely confers non-target site resistance
To test the hypothesis that loss-of-function of RTN4IP1L results in non-target site resistance, the strength of thaxtomin A resistance and the cross resistance to other herbicides were tested in Mprtn4ip1l mutants. Herbicide cross resistance and weak thaxtomin A resistance would be consistent with the hypothesis that RTN4IP1L loss-of-function mutations confer non-target site resistance.
The strength of thaxtomin A resistance in Mprtn4ip1l mutants was measured and compared to the strength of resistance to chlorsulfuron in chlorsulfuron target-site resistant mutants. Chlorsulfuron is a herbicide which targets the enzyme acetohydroxyacid synthase (AHAS), also known as acetolactate synthase (ALS). We generated 5 M. polymorpha chlorsulfuron-resistant target-site resistant mutants-Mpahas chlorsulfuron resistant (Mpahas chlr )-in a UV-B mutagenesis screen for resistance to chlorsulfuron. Each mutant carries a target-site resistance mutation in MpAHAS which confers chlorsulfuron resistance (S5 Fig). Mutations causing a Pro197Ser or Pro197Leu amino acid change, which have been observed to confer chlorsulfuron resistance in angiosperm weeds, were present in 4 out of the 5 Mpahas chlr mutants isolated in this screen [2,58,59].
To compare the strength of herbicide resistance in Mprtn4ip1l and Mpahas chlr lines, gemmae from 2 independent Mprtn4ip1l mutants (Mprtn4ip1l GE118 and Mprtn4ip1l GE149 ) were grown on solid medium supplemented with LD 100 or 5 x LD 100 of thaxtomin A (Fig 6A and 6B), and  Tab I in S1 Table). Tak-1 and Mprtn4ip1l GE149 gemmae were grown on solid medium supplemented with 0.1% DMSO for 14 days, then transferred to solid medium supplemented with 0.1% DMSO or 5 μM thaxtomin A for 2 days. Cellular fractions were extracted from the plants and analysed via LC-MS/MS. Stars represent the level of significance of the difference between Tak-1 and Mprtn4ip1l GE149 samples as determined by Student's t-tests: � = p < 0.05, �� = p < 0.01, ��� = p < 0.001. Error bars represent ± standard deviation (n = 3-6). B: DAB staining of Tak-1, Tak-2, Mprtn4ip1l GE118 , and Mprtn4ip1l GE149 . 12-day old gemmalings were stained with 3,3'-diaminobenzidine (DAB), which forms a brown precipitate upon reaction with H 2 O 2 . Stained plants were bleached to remove chlorophyll and imaged with a Keyence VHX-7000. Scale bars represent 500 μm and were added in Inkscape v1.0.1. C: Concentration of ROS in wild type and Mprtn4ip1l mutants. The concentration of ROS in 12-day old Tak-1, Tak-2, Mprtn4ip1l GE118 , and Mprtn4ip1l GE149 plants was determined by homogenising frozen samples in perchloric acid, followed by incubation with ferric xylenol-orange (FOX) and measurement of absorbance at 560 nm using an Ultrospec 3100 pro spectrophotometer. The absorbance measurements were converted to concentrations of ROS using a calibration curve (Tab J in S1 Table). Stars represent the level of significance of the difference between mutant and control lines (Tak-1 in red and Tak-2 in green) as determined by Student's t-tests: � = p < 0.05, �� = p < 0.01, ��� = p < 0.001. Error bars represent ± standard deviation (n = 3). D and E: Dose-response curves of the growth of wild type and Mprtn4ip1l mutants on H 2 O 2 (D) and paraquat (E). Gemmae from Tak-1 (orange), Tak-2 (green), Mprtn4ip1l GE118 (blue), and Mprtn4ip1l GE149 (purple) were grown for 10 days on solid medium supplemented with different concentrations of H 2 O 2 (D) (Tab K in S1 Table) or paraquat (E) (Tab L in S1 Table). The fitted curves and IC 50 values were calculated using the four-parameter log-logistic equation. Error bars represent ± standard deviation (n = 3). LD 100 is the lowest concentration at which area of all replicates was 0. RI is the ratio between each Mprtn4ip1l mutant's IC 50  gemmae from 5 independent Mpahas chlr mutant lines were grown on solid medium supplemented with LD 100 or 5 x LD 100 of CS (Fig 6C and 6D). The average area of living plant tissue was quantified after 12 days of growth (Tabs M and N in S1 Table).
Mprtn4ip1l mutants are significantly larger than wild type plants on the thaxtomin A LD 100 , confirming their thaxtomin A resistance (Fig 6A and 6B). Similarly, Mpahas chlr mutants are significantly larger than wild type on the chlorsulfuron LD 100 (Fig 6C and 6D). However, Mprtn4ip1l mutants grow significantly less on thaxtomin A LD 100 than in control conditions (Fig 6A and 6B). Conversely, Mpahas chlr mutants are either statistically indistinguishable or significantly larger when grown on chlorsulfuron LD 100 than in control conditions (Fig 6C  and 6D). This suggests that Mprtn4ip1l mutants are weakly resistant to thaxtomin A, whereas Mpahas chlr mutants are strongly resistant to chlorsulfuron. This is consistent with the hypothesis that Mprtn4ip1l mutants are non-target site resistant.
Mprtn4ip1l mutants die when grown on five times the LD 100 of thaxtomin A (Fig 6A and  6B), whereas Mpahas chlr mutants do not die on five times the LD 100 of chlorsulfuron (Fig 6C  and 6D). The inability of Mprtn4ip1l mutants to survive 5 x LD 100 is also consistent with the hypothesis that Mprtn4ip1l mutants are weakly resistant to thaxtomin A. Resistance from lossof-function mutations in Mprtn4ip1l is therefore more likely to be non-target site resistance than target-site resistance.
To further test the hypothesis that loss-of-function of Mprtn4ip1l function causes non-target site resistance against thaxtomin A, we tested if Mprtn4ip1l mutants were resistant to other herbicides (cross resistance). To quantify the cross resistance of Mprtn4ip1l mutants to different herbicides, the resistance of Mprtn4ip1l GE118 and Mprtn4ip1l GE149 to the herbicides isoxaben, dichlobenil, and 2,4-D was tested. Isoxaben targets cellulose biosynthesis, targeting the cellulose synthase subunits CESA1, CESA3, and CESA6 [60]. Dichlobenil is also classed as a cellulose biosynthesis inhibitor, but its target protein is still unknown [61]. 2,4-D is an auxin mimic, causing deregulation of the auxin signalling pathway [62]. Gemmae from Tak-1, Tak-2, Mprtn4ip1l GE118 , and Mprtn4ip1l GE149 were grown on solid medium supplemented with different doses of isoxaben (Fig 6E), dichlobenil (Fig 6F), or 2,4-D (Fig 6G). After 10 days of growth, the area of living tissue was quantified (Tabs O, P, and Q in S1 Table), and doseresponse curves were fitted using the four-parameter log-logistic equation. IC 50 , LD 100 , and RI for wild type and Mprtn4ip1l mutant lines were calculated from the resulting dose-response curves (Fig 6E, 6F and 6G).
Mprtn4ip1l mutant lines were resistant to isoxaben, Mprtn4ip1l GE118 with an RI of 4.22 ± 8.96 and Mprtn4ip1l GE149 with an RI of 3.94 ± 6.01. Although these RI values are not statistically significant, the LD 100 of both Mprtn4ip1l mutant lines was higher than wild type; both were alive when grown on 20 μM isoxaben which is lethal to Tak-1 and Tak-2 ( Fig 6E).

Fig 6. Loss-of-function of MpRTN4IP1L likely confers non-target site resistance. A and B:
Wild type gemmae (Tak-1 and Tak-2) and gemmae from Mprtn4ip1l mutants were grown on solid medium supplemented with LD 100 (5 μM) or 5 x LD 100 (25 μM) of thaxtomin A (Tab M in S1 Table). (A) Gemmalings were imaged using a Keyence VHX-7000. Scale bars represent 2 mm and were added in Inkscape v1.0.1. (B) Error bars represent ± standard deviation (n = 3). Stars represent the level of significance (as determined by Student's t-tests) of the difference between mutant and control lines (comparison to Tak-1 in red and Tak-2 in green), or between individuals of the same genotype subjected to different treatments (black): � = p < 0.05, �� = p < 0.01, ��� = p < 0.001. C and D: Wild type gemmae (Tak-1 and Tak-2) and gemmae from Mpahas chlr mutants were grown on solid medium supplemented with LD 100 (33 nM) or 5 x LD 100 (165 nM) of chlorsulfuron (Tab N in S1 Table). (C) Gemmalings were imaged using a Keyence VHX-7000. Scale bars represent 2 mm and were added in Inkscape v1.0.1. (D) Error bars represent ± standard deviation (n = 3). Stars represent the level of significance (as determined by Student's t-tests) of the difference between mutant and control lines (comparison to Tak-1 in red and Tak-2 in green), or between individuals of the same genotype subjected to different treatments (black): � = p < 0.05, �� = p < 0.01, ��� = p < 0.001. E,F, and G: Dose-response curves of the growth of Mprtn4ip1l mutants on herbicides with different targets and tables of corresponding IC 50 , LD 100 , and RI values. Gemmae from Tak-1 (orange), Tak-2 (green), Mprtn4ip1l GE118 (blue), and Mprtn4ip1l GE149 (purple) were grown for 12 days on solid medium supplemented with different doses of isoxaben (E) (Tab O in S1 Table), dichlobenil (F) (Tab P in S1 Table), or 2,4-D (G) (Tab Q in S1 Table). The fitted curves and IC 50 values were calculated using the four-parameter log-logistic equation. Error bars represent ± standard deviation (n = 3). LD 100 is the lowest concentration at which area of all replicates was 0. RI was calculated as the ratio between each Mprtn4ip1l mutant's IC 50  These data demonstrate that Mprtn4ip1l GE118 and Mprtn4ip1l GE149 are resistant to isoxaben. The cross resistance of Mprtn4ip1l mutants to isoxaben is consistent with the hypothesis that loss of MpRTN4IP1L function confers non-target site resistance.
Mprtn4ip1l mutant lines were not resistant to dichlobenil or 2,4-D. The RI values of Mprtn4ip1l GE118 and Mprtn4ip1l GE149 on dichlobenil are not statistically significant, and neither mutant can survive at a higher concentration of dichlobenil than wild type (Fig 6F). Mprtn4ip1l GE118 has an RI value of 0.33 ± 0.20 on 2,4-D so is significantly 2,4-D sensitive, while the RI of Mprtn4ip1l GE149 on 2,4-D is not statistically significant. Neither Mprtn4ip1l mutant has a higher LD 100 than wild type on 2,4-D (Fig 6G). Mprtn4ip1l lines are therefore not cross resistant to dichlobenil or 2,4-D.
Taken together, the weak thaxtomin A resistance and isoxaben cross resistance of Mprtn4ip1l mutants demonstrate that loss-of-function of MpRTN4IP1L is likely to confer non-target site resistance. It is possible that this non-target site resistance results from ROS-induced alterations to the cell wall that reduce the effect of thaxtomin A. Non-target site resistance is however often caused by chemical modification of a herbicide resulting in the formation of a non-toxic product, or by reduced herbicide uptake [1,6]. Therefore, an alternative hypothesis to the ROS-induced cell wall alteration hypothesis is that non-target site resistance to thaxtomin A in Mprtn4ip1l mutants is caused by a change in herbicide metabolism or uptake.

Mprtn4ip1l mutants are defective in thaxtomin A metabolism
We set out to test the hypothesis that non-target site resistance in Mprtn4ip1l mutants results from differences in the metabolism or uptake of thaxtomin A between wild type and mutants. To test this hypothesis, thaxtomin A and putative thaxtomin A metabolites in thaxtomin Atreated samples of Tak-2 wild type and Mprtn4ip1l GE118 mutants were detected via a precursor ion analysis for the two most abundant fragment ions of thaxtomin A using LC/MS-MS ( Fig  7A and 7B). There is a peak of pure thaxtomin A (retention time-RT-8.99 min) and a peak of an unknown compound (RT 9.25 min) in the chromatogram of wild type Tak-2 ( Fig 7A). Both peaks are absent in the analysis of the untreated samples (S6A Fig). The two prominent fragment ions of thaxtomin A are present in the MS/MS spectra of the unknown compound, suggesting that it is a thaxtomin A derivative (S6B and S6C Fig). The molecular mass of the unknown compound is 45 Da less than pure thaxtomin A. An explanation for these observations is that the unknown compound is a derivative of thaxtomin A that lacks the NO 2 group (delta mass of 45 Da) (Fig 7C and 7D). This suggests that wild type M. polymorpha metabolizes thaxtomin A by removing the NO 2 group of the 4-nitro-tryptophan moiety. These two peaks are present in the chromatogram from Mprtn4ip1l GE149 . However, the RT 9.25 peak is considerably smaller in Mprtn4ip1l GE149 than wild type (Fig 7B). This suggests that Mprtn4ip1l mutants accumulate less of the derivative and are therefore defective in thaxtomin A metabolism.
To quantify the potential defect in thaxtomin A metabolism in Mprtn4ip1l mutants, we measured the relative amounts of thaxtomin A (RT 8.99) and the putative thaxtomin A derivative (RT 9.25 min) in wild type and mutant plants using targeted LC-MS/MS employing selected reaction monitoring (Fig 7E and 7F, and Tabs R and S in S1 Table). The concentration of thaxtomin A is significantly higher in Mprtn4ip1l GE118 than in Tak-2 (p < 0.05) (Fig 7E). This is consistent with the hypothesis that the chemical modification of thaxtomin A is defective in the mutants (Fig 7E). Furthermore, this is inconsistent with the hypothesis that Mprtn4ip1l mutants uptake less thaxtomin A than wild type plants, suggesting that their resistance is not conferred by reduced herbicide uptake.
The ion counts of the putative thaxtomin A derivative (RT 9.25) are significantly lower in Mprtn4ip1l GE118 than in Tak-2 (p < 0.05), consistent with the hypothesis that the chemical modification of thaxtomin A is defective in the mutants (Fig 7F). Although the difference in the concentration of thaxtomin A between the samples cannot be directly compared to the difference in ion counts of the unknown compound, these data suggest that wild type plants metabolize thaxtomin A to form a derivative compound, but this metabolism is defective in Mprtn4ip1l mutants. Given that Mprtn4ip1l mutants are resistant to thaxtomin A, this is consistent with the hypothesis that thaxtomin A metabolism has a toxic effect on the plant.
To independently verify that the metabolism of thaxtomin A is defective in Mprtn4ip1l mutants, the presence of thaxtomin A and its putative metabolite were also detected via precursor ion analysis in Tak-1 wild type and Mprtn4ip1l GE149 mutants. The peak corresponding to the putative metabolite was stronger in Tak-1 than in Mprtn4ip1 GE149 , supporting the hypothesis that thaxtomin A metabolism is defective in Mprtn4ip1 mutants (S6D Fig). Together these data indicate that a thaxtomin A derivative that accumulates in wild type plants accumulates at much lower levels in the resistant Mprtn4ip1l mutants. It is hypothesized that the metabolism of thaxtomin A is toxic to plants. According to this hypothesis, Mprtn4ip1l mutants are resistant to thaxtomin A because they are defective in thaxtomin A metabolism. The high levels of ROS produced in Mprtn4ip1l mutants because of low levels of CoQ may inhibit the metabolism of thaxtomin A.

Discussion
We report the discovery of a novel mechanism of non-target site resistance caused by loss-offunction mutations in the MpRTN4IP1L gene using a forward genetic screen. We demonstrate that loss-of-function Mprtn4ip1l mutants are resistant to thaxtomin A. Resistance is relatively weak and mutants are also cross resistant to isoxaben which suggests that the relative insensitivity to thaxtomin A is not because of a form of target-site resistance. Instead, there is evidence that thaxtomin A resistance in the Mprtn4ip1l mutants is the result of non-target site resistance. Consistent with the hypothesis that the resistance is not due to target site mutation but instead an example of non-target site resistance is the observation that the mutants produce more ROS than wild type plants. Furthermore, the metabolism of thaxtomin A is different in the resistant mutants than in wild type, consistent with non-target site resistance.
While our data suggest that the resistance resulting from loss-of-function mutations in the MpRTN4IP1L gene is not caused by a thaxtomin A-insensitive target, there are a number of possible mechanisms which could explain non-target site resistance in these mutants. One hypothesis to explain thaxtomin A resistance is that metabolism of the herbicide is defective in Mprtn4ip1l mutants. If correct, this would be a rare instance of a decrease in herbicide metabolism conferring non-target site resistance, which is usually associated with increased metabolism [6].
If resistance were the result of defective thaxtomin A metabolism, then it is possible that thaxtomin A acts as a proherbicide. Proherbicides are non-toxic molecules that require

Fig 7. Mprtn4ip1l mutants are defective in thaxtomin A metabolism. A and B: LC-MS/MS analysis of thaxtomin A and its putative metabolite in Tak-2 (A)
and Mprtn4ip1l GE118 (B) extracts. Gemmae from each line were grown on solid medium supplemented with 0.1% DMSO for 14 days, then transferred to solid 1 medium supplemented with 5 μM thaxtomin A and grown for 2 days. Cellular fractions were extracted from thaxtomin A-treated samples and a precursor ion analysis was conducted via LC-MS/MS (n = 6). Chromatograms of the total ion currents of precursor ion scanning for m/z 247.1 are depicted. The peak eluting at a retention time of 8.99 min. is thaxtomin A (m/z 439.1), the peak eluting at a retention time of 9.25 (m/z 394.1) is its putative metabolite. C: Molecular structure of thaxtomin A D: Molecular structure of putative thaxtomin A metabolite E: Concentration of pure thaxtomin A (TXTA) in DMSO and herbicidetreated samples of Tak-2 and Mprtn4ip1l GE118 as quantified by targeted LC/MS-MS (Tab R in S1 Table). F: Ion counts (normalised by weight) of the putative thaxtomin A metabolite in DMSO and herbicide-treated samples of Tak-2 and Mprtn4ip1l GE118 as quantified by targeted LC-MS/MS (Tab S in S1 Table).
https://doi.org/10.1371/journal.pgen.1010423.g007 metabolism to form a herbicidally active, toxic molecule [63]. We demonstrated that a thaxtomin A metabolite is present at significantly lower levels in Mprtn4ip1l mutants than in wildtype plants. Therefore, if thaxtomin A is a proherbicide and non-toxic, the metabolite could be toxic and consequently herbicidal. If true, the decreased thaxtomin A metabolism in Mprtn4ip1l mutants would prevent the production of herbicidally active molecules that contribute to thaxtomin A non-target site resistance.
While our data are consistent with the hypothesis that thaxtomin A is a proherbicide, there is evidence that refute this hypothesis. The difference in molecular mass (45 Da) between thaxtomin A and the metabolite are consistent with the hypothesis that the metabolite lacks the NO 2 group, although this has not been verified. The NO 2 group of thaxtomin A has been shown to be required for phytotoxicity [22]. If the NO 2 group is required for phytotoxicity, and if the thaxtomin A metabolite lacks the NO 2 group, it would suggest that the observed thaxtomin A metabolite is not toxic. This would be inconsistent with thaxtomin A acting as a proherbicide. Therefore, while we propose that non-target site resistance in Mprtn4ip1l mutants results from defective thaxtomin A metabolism, we have not proved that thaxtomin A acts as a proherbicide.
It is formally possible that the defective thaxtomin A metabolism observed in Mprtn41ipl mutants may result from elevated levels of ROS in the mutant. Mprtn4ip1l mutants are defective in the biosynthesis of coenzyme Q, an electron transporter in oxidative phosphorylation. Consequently, mutants produce more ROS than wild type plants. The increased ROS levels in Mprtn4ip1l mutants may inhibit the chemical reactions involved in the metabolism of thaxtomin A observed in wild type. It is therefore possible that ROS-inhibited thaxtomin A metabolism could lead to thaxtomin A resistance.
It is also possible that high ROS levels in Mprtn4ip1l mutants confer herbicide resistance by altering the properties of the cell wall. Thaxtomin A is classed as a cellulose biosynthesis inhibitor [29,30]. It is therefore possible that physiological changes to the cell wall which hinder the action of thaxtomin A could confer non-target site thaxtomin A resistance. It is well established that ROS can cause physiological changes to the cell wall, such as cell wall stiffening upon abiotic stress by cross-linking or altering cell wall proteins [57,[64][65][66][67]. It is therefore possible that ROS-induced cell wall alterations could confer resistance to thaxtomin A in Mprtn4ip1l mutants.
The effects of ROS and thaxtomin A on A. thaliana cell walls are consistent with the hypothesis that ROS-induced cell wall alterations could confer thaxtomin A resistance. A. thaliana cell walls stiffen upon ROS treatment but loosen upon treatment with the cellulose biosynthesis inhibitors thaxtomin A and isoxaben [57,68]. Thaxtomin A and isoxaben may therefore have antagonistic effects to ROS on the cell wall. Furthermore, ROS pretreatment prevents thaxtomin A-induced cell wall loosening and confers resistance to thaxtomin A. [57]. This suggests that ROS-induced cell wall stiffening can protect against the toxic effect of thaxtomin A. Therefore, according to the ROS-induced cell wall alteration hypothesis, increased ROS levels in Mprtn4ip1l mutants cause cell wall stiffening which prevents the cell wall loosening action of thaxtomin A and isoxaben. The ROS-induced cell wall alteration hypothesis therefore also accounts for the cross-resistance of Mprtn4ip1l mutants to isoxaben.
According to the ROS-induced cell wall alteration hypothesis, cross resistance of Mprtn4ip1l to thaxtomin A and isoxaben would be due to the cell wall stiffening effect of ROS which counters the action of these two cellulose biosynthesis inhibitors. However, Mprtn4ip1l mutants are not resistant to the cellulose biosynthesis inhibitor dichlobenil. Thaxtomin A and isoxaben are classed as group 1 cellulose biosynthesis inhibitors and dichlobenil is classed as group 2 [30]. Group 1 cellulose biosynthesis inhibitors cause depletion of cellulose synthase complexes from the plasma membrane; group 2 cellulose biosynthesis inhibitors however cause accumulation of cellulose synthase complexes at the plasma membrane [30]. Given the opposite effects of Group 1 and Group 2 cellulose biosynthesis inhibitors on cellulose synthase enzymes, it is possible that some resistance alleles confer resistance to one group of cellulose biosynthesis inhibitors but not the other. Therefore, according to the ROS-induced cell wall alteration hypothesis, the cross-resistance pattern observed in Mprtn4ip1l mutants results from ROS-induced changes to the cell wall conferring resistance to group 1 but not group 2 cellulose biosynthesis inhibitors.
The mechanism of non-target site resistance reported here is also likely to operate in angiosperm weeds and not be restricted to M. polymorpha or other liverworts. First, the RTN4IP1L gene is found throughout the angiosperms and therefore agricultural weeds, where mutation would likely result in resistance. Second, our data indicate that loss-of-function mutations in the PAM16 gene confer weak thaxtomin A resistance in both M. polymorpha and A. thaliana [20]. Therefore, non-target site resistance caused by loss of PAM16 function in M. polymorpha also operates in angiosperms, demonstrating that some mechanisms of non-target site resistance are conserved between M. polymorpha and angiosperms. Furthermore, we report that four chlorsulfuron-resistant M. polymorpha UV-B mutants carry Pro197Ser or Pro197Leu mutations, which are among the most frequently reported mechanisms of chlorsulfuron resistance observed in chlorsulfuron-resistant angiosperm weeds [2,58,59]. Together, these data demonstrate that several herbicide resistance mechanisms are conserved between M. polymorpha and angiosperms, validating the use of M. polymorpha as a model with which to identify novel mechanisms of non-target site resistance which could evolve in weeds.
Most Mprtn4ip1l mutants are smaller than wild type plants in control conditions, suggesting that loss of MpRTN4IP1L function would incur a considerable fitness cost in the wild. Given this, the evolution of Mprtn4ip1l complete loss-of-function based resistance in the field is unlikely. However, some herbicide resistance alleles incur small fitness costs but are nevertheless maintained in weed populations due to the intense selection imposed by herbicide treatment [69][70][71][72]. We propose that weak loss-of-function alleles of RTN4IP1L (or mutations which affect the same molecular processes), which cause negligible or small fitness costs, could evolve in weed populations treated with thaxtomin A. If weak loss-of-function of RTN4IP1L is identified in resistant weeds in the field, our findings can inform weed management practices. Based on our characterization of phenotypes of Mprtn4ip1l loss-of-function mutants, weeds with RTN4IP1L loss-of-function mutations cannot be controlled by isoxaben. However, our data suggest that weeds that evolve thaxtomin A resistance through loss-of-function mutations in the RTN4IP1L gene, or through the overproduction of ROS, are hypersensitive to paraquat, and are therefore likely hypersensitive to other herbicides that kill plants by ROS overproduction such as protoporphyrinogen oxidase inhibitors or glufosinate [53][54][55][56].
Our findings demonstrate the potential to identify novel mutations conferring non-target site herbicide resistance which have been difficult to identify and characterize to date. Our methodology can uncover novel mechanisms of non-target site resistance caused by loss-offunction mutations which could inform efficient weed management.

Plant lines and growth conditions
Wild-type plants were M. polymorpha laboratory accessions Takaragaike-1 (Tak-1) and Takaragaike-2 (Tak-2) [73]. Mppam16 and Mprtn4ip1l lines were generated via CRISPR-Cas9 mutagenesis of wild-type M. polymorpha spores (from a cross between Tak-1 and Tak-2). All lines were maintained via asexual propagation of gemmae or thallus excision in the case of mutants which did not produce gemmae. Plants were grown on solid ½ Gamborg medium. All plant material used for experiments was grown in a growth chamber at 23˚C under 24-hour 10-30 μmol m -2 s -1 white light. Plant material to be stored long-term was kept in a growth chamber at 17˚C under 6 hours 10-30 μmol m -2 s -1 white light and 18 hours dark. To produce spores, plants were grown on a mixture of John Innes n. 2 compost and fine vermiculite at a ratio of 2:1. Plants were grown in growth chambers at 23˚C and under 16 hours white light and continuous far-red light to induce transition to the reproductive phase. Once mature, water was pipetted onto the antheridia produced by male plants to collect sperm, then transferred to archegonia produced by female plants. Sporangia were harvested before bursting and stored fresh at 4˚C in sterile dH 2 O.

Gemmaling dose-response assays
Gemmae were grown on solid ½ Gamborg medium supplemented with different concentrations of herbicide dissolved in DMSO (n = 2-18; see individual figure legends for experimentspecific biological replicate numbers). Gemmalings were imaged using a Berthold Nightowl II LB 983 In Vivo Imaging System (Berthold, Bad Wildbad, Germany). Images were taken after exposing gemmalings to 120 s white light. The imaging system detects chlorophyll autofluorescence (560 nm). The lateral area of autofluorescing (living) tissue was determined using the indiGo software package (Berthold, Bad Wildbad, Germany). Dose-response curves were generated using the ggplot2 and drc packages in R [74,75]; the four-parameter log-logistic equation with a fixed lower limit of 0 was used to fit a dose-response curve to the data.

Phylogenetic analysis of PAM16 and RTN4IP1 homologues
Homologues of the A. thaliana At3G59280 protein or of the H. sapiens RTN4IP1.1 protein were identified by protein BLAST search against the reference proteomes of various species [37]. Only homologues with an E value less than 1E-5 were used to construct the trees. Homologues were aligned via the L-INS-i strategy using MAFFT version 7 [76]. The alignments were manually trimmed using BioEdit7.2. A maximum likelihood tree was constructed with PhyML 3.0 using an estimated gamma distribution parameter and the LG model of amino acid substitution [38]. A non-parametric approximate likelihood ratio test based on a Shimodaira-Hasegawa-like procedure was used to calculate branch support values using PhyML 3.0 [38]. The trees were visualized in FigTree v1.4.4 and rooted with the respective Saccharomyces cerevisiae homologues (PAM16 or Yim1p).

CRISPR-Cas9 mutagenesis
CRISPR-Cas9 mutants were generated based on the protocols described in [77] and [78]. Guide RNAs (sgRNAs) annealing to different parts of the genes MpPAM16 (Mp3g09390.1) or MpRTN4IP1L (Mp3g19030.1) were designed to be 5' of an "NGG" site (PAM sequence) as required by the CRISPR-Cas9 system [77]. Guide RNAs were cloned into the pHB453 or pMpGE_En03 constructs [78] to generate the entry vectors using primers "sgRNA" (S3 Table). The expression vectors were generated by LR reaction of the destination vector pMpGE010 and the entry vectors. Vectors were transformed into Escherichia coli One Shot OmniMAX 2 T1R (Thermo Fisher Scientific: Cat. # C854003).
The expression vectors were transformed into Agrobacterium tumefaciens strain GV3101. Wild-type M. polymorpha spores from a cross between Tak-1 and Tak-2 accessions were transformed as per the protocol in [79]. Transformants were selected by growth on solid ½ Gamborg medium supplemented with 10 mg/l hygromycin B (Melford: CAS# 31282-04-9) to select for plants with the CRISPR-Cas9 construct insertion, and 100 mg/l cefotaxime (Melford: CAS# 64485-93-4) to kill any remaining A. tumefaciens.
To identify potential mutations in the MpPAM16 or MpRTN4IP1L genes, genomic regions including the loci targeted by the guide RNAs were amplified by PCR and Sanger sequenced using primers "PAM16_G" or "RTN4IP1L_G" (S3 Table). Sanger sequencing of the purified PCR products was conducted by Source BioScience or by the Molecular Biology service at Vienna BioCenter Core Facilities (VBCF), member of the Vienna BioCenter (VBC), Austria. Mutants were named according to the gene mutated (Mppam16 or Mprtn4ip1l), followed by "GE" (gene edited), the number of the guide RNA targeting the locus in which the mutation is found, and the number transformant that was genotyped (e.g. Mprtn4ip1l GE149 -a gene edited transformant with a mutation in the RTN4IP1L gene at the locus targeted by sgRNA 1 which was the 49 th transformant genotyped).

Generation of herbicide-resistant M. polymorpha lines via UV-B mutagenesis
Spores were mutagenized according to the protocol in [40]. Sterilised fresh wild type M. polymorpha spores were plated on solid modified Johnson's medium supplemented with thaxtomin A (5 μM-LD 100 ) or chlorsulfuron (140 nM-4 x LD 100 ). The spores were exposed to UV-B irradiation (302 nm) for an exposure time corresponding to or near the LD 50 (an exposure time of UV-B at which 50% of spores are killed-S2 Fig) using a UVP BioDoc-It. Plates of mutagenized spores were wrapped in aluminium foil and left in the dark overnight. Plates were then unwrapped and placed in a Sanyo growth chamber at 23˚C in 24-hour light. After 14 days of growth, plates were screened for surviving sporelings. Surviving sporelings were transferred to solid modified Johnson's medium supplemented with fresh herbicide. Plants which survived this second transfer onto a lethal dose of herbicide were maintained. To test for retention of resistance through asexual reproduction, gemmae from these plants were plated onto a lethal dose of herbicide. Lines whose resistance was inherited were classified as herbicide resistant lines and maintained. Lines with thaxtomin A resistance were named Mpthar (thaxtomin A resistant).
Chlorsulfuron-resistant lines were tested for target-site resistance by PCR amplification and Sanger sequencing of regions of the MpAHAS gene using primers GCS_Fw and GCS_Rv (S3 Table). Sanger sequencing of the purified PCR products was carried out by Source BioScience. Lines with target-site resistance to chlorsulfuron were named Mpahas chlr (Mpahas chlorsulfuron resistant).

Genomic DNA sequencing
Wild-type (Tak-1, Tak-2, OxTak1F, and OxTak2M) lines and herbicide-resistant lines (Mpthar and Mpahas chlr ) were grown on solid ½ Gamborg medium for three weeks on in a growth chamber at 23˚C under 24-hour 10-30 μmol m -2 s -1 white light. After three weeks of growth, plant material was harvested and flash frozen in liquid nitrogen. Samples were ground in liquid nitrogen and genomic DNA was extracted using the Qiagen DNeasy Plant Maxi kit (Cat. # 68163) according to the kit protocol. After elution, the gDNA was cleaned up and concentrated using the Zymo Genomic DNA Clean & Concentrator kit (Cat. # D4010). The gDNA was eluted in nuclease-free water.
The concentration of the gDNA was checked using a Nanodrop 1000 spectrophotometer and a Qubit 2.0 fluorometer according to the instruction manuals. The gDNA quality was checked by running 2 μl DNA (approximately 10 μg μl -1 ) on a 0.7% agarose gel at 70 V for 45 minutes and checking for degradation. Genomic DNA samples which had passed quality control checks were sent for sequencing to the Next Generation Sequencing Facility at Vienna BioCenter Core Facilities (VBCF), member of the Vienna BioCenter (VBC), Austria. DNA Libraries were prepared using the Westburg NGS Library Prep kit and fragment size was determined using a BioLabTech Fragment Analyzer. The DNA was sequenced on an Illumina NovaSeq 6000 SP flowcell using 150 bp paired-end reads.

Non-allelism based SNP discovery pipeline
The "non-allelism based SNP discovery" pipeline from [40] was used with slight adaptations to identify candidate resistance-conferring SNPs in herbicide resistant lines. Raw reads were trimmed to remove low-quality reads and NEB adaptors using Trimmomatic 0.38 [80]. The read coverage was then normalized using khmer 2.1.2 [81]. Reads were aligned to the Marchantia polymorpha reference genome (MpTak1 v5.1 plus the female chromosome from v3.1; now known as MpTak v6.1) using bowtie2. Reads were sorted by position and reads from different lanes were merged using samtools 1.10. The variant call analysis was carried out using samtools and bcftools (ploidy defined as haploid, multiallelic/rare variant call option selected) to generate bcf files listing mismatches between the sequencing reads and the reference genome for each line [82]. For each bcf file, only mismatches which were supported by 7-100 reads and where the number of reads supporting high quality alternative alleles � the number of reads supporting high quality reference alleles were retained (DP4[2]+DP4 [3])/sum(DP4)> 0.5). Non-canonical UV-B induced mismatches (not C ➔ T or G ➔ A) were filtered out using bcftools. Filtered mismatches were retained only if they were not present in other sequenced non-allelic lines without that phenotype. The resulting lists were filtered to retain only mismatches present in a CDS using bash scripting, then filtered manually in Interactive Genomics Viewer v 2.10 to remove mismatches which did not induce an amino acid change and which were misaligned or poor quality.
SNPs identified in the MpRTN4IP1L gene in Mpthar2, Mpthar4 and Mpthar6 by the pipeline were confirmed by PCR amplification and Sanger sequencing using primers "GRT" (S3 Table). Sanger sequencing was carried out by Source BioScience.

Metabolomic analysis of pure and modified thaxtomin A in M. polymorpha thallus
Gemmae were grown on autoclaved cellophane (AA Packaging Limited, Preston, UK) on solid ½ Gamborg medium supplemented with 0.1% DMSO for 14 days at 23˚C in 24 h light (n = 6). Gemmalings were then transferred (by transfer of the cellophane disc) onto solid ½ Gamborg medium supplemented with either 0.1% DMSO or 5 μM thaxtomin A and grown in these conditions for 2 days at 23˚C in 24 h light. Treated and untreated gemmalings were harvested by flash freezing in liquid nitrogen.
For thaxtomin A analysis, frozen tissue samples were ground in liquid nitrogen and homogenized in ice-cold extraction solvent (2:1:1 methanol:acetonitrile:H 2 O, v/v) for 2 min at 4˚C followed by incubation at -20˚C for one hour. Samples were centrifuged at full speed at 4˚C for 3 mins. Supernatants were collected and kept at -20˚C. Pellets were redissolved in ice cold 80% (v/v) methanol and homogenized for 1 min at 4˚C followed by incubation at -20˚C for one hour. Samples were centrifuged at full speed at 4˚C for 3 mins. Supernatants were collected and added to the supernatants from the first extraction step. Samples were incubated at -20˚C for 2 hours. Samples were then centrifuged at full speed at 4˚C for 10 mins. Supernatants were collected and shock frozen in liquid nitrogen. For coenzyme Q10 analysis, a chloroform/ methanol extraction from frozen plant tissue was carried out. Coenzyme Q10 was analyzed after chloroform/methanol extraction from powdered frozen plant tissue. 100 μl of the lower phase was diluted with 50 μl methanol and 1 μl was directly injected onto a Kinetex (Phenomenex) C8 column (100 Å, 100 x 2.1 mm) and employing a 5-minute-long linear gradient from 90% A (50% acetonitrile, 49.9% water, 0.1% formic acid, 10 mM ammonium formate) to 95% B (10% acetonitrile 88% isopropanol, 1.9% water, 0.1% formic acid, 10 mM ammonium formate) at a flow rate of 100 μl/min. Detection and quantification was done by LC-MS/MS, employing the selected reaction monitoring (SRM) mode of a TSQ Altis mass spectrometer (Thermo Fisher Scientific), using the transition 863.8 m/z to 197.1 m/z in the positive ion mode.

DAB staining
3,3'-diaminobenzidine (DAB) (Sigma Aldrich: D8001) was used to prepare a DAB staining solution according to [50]. Tak-1, Tak-2, Mprtn4ip1l GE118 and Mprtn4ip1l GE149 11-day old plants were incubated in 3 ml solution in 24 well plates and vacuum infiltrated for 5 minutes (n = 3). The plates were covered in aluminium foil and incubated for 1.5 hours shaking at 100 rpm. After incubation, the staining solution was replaced by bleaching solution prepared according to [50]. Plants were bleached for 2 hours. Bleached gemmae were imaged using a Keyence VHX-7000.

Ferric-xylenol orange (FOX) assay
A modified FOX assay was carried out according to the protocol in [52]. The FOX working solution was prepared as follows: 500 μM ammonium ferrous sulphate; 400 μM xylenol orange; 200 mM sorbitol; 25 mM H 2 SO 4 . Twelve-day old gemmalings of Tak-1, Tak-2, Mprtn4ip1l GE118 and Mprtn4ip1l GE149 were frozen and ground in liquid nitrogen and homogenized in 200 mM perchloric acid (n = 3). Samples were centrifuged at 1000g for 5 minutes at 4˚C. 500 μl of the supernatant were mixed with 500 μl FOX working solution. Samples were incubated at room temperature in the dark for 30 minutes followed by quantification of the absorbance at 560 nm using an Ultrospec 3100 pro spectrophotometer. The concentration of H 2 O 2 in each sample was calculated using a calibration curve of known concentrations of H 2 O 2 quantified using the same protocol.

Media
Modified Johnson's medium was prepared as described in [79].
½ Gamborg medium was prepared using the following components (Table 2) Wild type spores (from a cross between Tak-1 and Tak-2) were plated on solid medium and subjected to UV-B irradiation for different lengths of time using a UVP BioDoc-It. Spores were imaged after 14 days of growth using a Leica DFC310 FX stereomicroscope. The total area of spores on each plate was determined using ImageJ (Tab D in S1 Table). The fitted curve was calculated using the four-parameter log-logistic equation from the drc package in R. IC 50 (exposure time of UV-B at which the average sporeling area on the plate was reduced by 50%) was calculated from the dose-response curve as approximately 110 s. Error bars represent ± standard deviation (n = 3). (TIF)