RNA-Seq transcriptome analysis of Amaranthus palmeri with differential tolerance to glufosinate herbicide

Amaranthus palmeri (Amaranthaceae) is a noxious weed in several agroecosystems and in some cases seriously threatens the sustainability of crop production in North America. Glyphosate-resistant Amaranthus species are widespread, prompting the use of alternatives to glyphosate such as glufosinate, in conjunction with glufosinate-resistant crop cultivars, to help control glyphosate-resistant weeds. An experiment was conducted to analyze the transcriptome of A. palmeri plants that survived exposure to 0.55 kg ha-1 glufosinate. Since there was no record of glufosinate use at the collection site, survival of plants within the population are likely due to genetic expression that pre-dates selection; in the formal parlance of weed science this is described as natural tolerance. Leaf tissues from glufosinate-treated and non-treated seedlings were harvested 24 h after treatment (HAT) for RNA-Seq analysis. Global gene expression was measured using Illumina DNA sequence reads from non-treated and treated surviving (presumably tolerant, T) and susceptible (S) plants. The same plants were used to determine the mechanisms conferring differential tolerance to glufosinate. The S plants accumulated twice as much ammonia as did the T plants, 24 HAT. The relative copy number of the glufosinate target gene GS2 did not differ between T and S plants, with 1 to 3 GS2 copies in both biotypes. A reference cDNA transcriptome consisting of 72,780 contigs was assembled, with 65,282 sequences putatively annotated. Sequences of GS2 from the transcriptome assembly did not have polymorphisms unique to the tolerant plants. Five hundred sixty-seven genes were differentially expressed between treated T and S plants. Of the upregulated genes in treated T plants, 210 were more highly induced than were the upregulated genes in the treated S plants. Glufosinate-tolerant plants had greater induction of ABC transporter, glutathione S-transferase (GST), NAC transcription factor, nitronate monooxygenase (NMO), chitin elicitor receptor kinase (CERK1), heat shock protein 83, ethylene transcription factor, heat stress transcription factor, NADH-ubiquinone oxidoreductase, ABA 8’-hydroxylase, and cytochrome P450 genes (CYP72A, CYP94A1). Seven candidate genes were selected for validation using quantitative real time-PCR. While GST was upregulated in treated tolerant plants in at least one population, CYP72A219 was consistently highly expressed in all treated tolerant biotypes. These genes are candidates for contributing tolerance to glufosinate. Taken together, these results show that differential induction of stress-protection genes in a population can enable some individuals to survive herbicide application. Elevated expression of detoxification-related genes can get fixed in a population with sustained selection pressure, leading to evolution of resistance. Alternatively, sustained selection pressure could select for mutation(s) in the GS2 gene with the same consequence.

Introduction glufosinate-resistant crops is achieved using transgenic methods to insert the bar or pat (phosphinothricin acetyltransferase) gene from a bacterium to the plant's genome, allowing detoxification of glufosinate by acetylation [27].
Differential responses to glufosinate in weeds have been attributed to several factors including light, temperature, humidity, growth stage, application rate, application timing, species, and variation in the level of herbicide absorption, translocation, and metabolism [28,29]. Control of A. palmeri by glufosinate is variable [30][31][32]. A previous study reported higher uptake, mobility, and metabolism of glufosinate in A. palmeri compared to the more susceptible Ipomoea lacunosa [28]. As is commonly observed, differential tolerance to herbicides are often due to non-target-site (NTS) mechanisms, involving the detoxification of herbicide by biochemical modification and/or the compartmentation of the herbicide and its metabolites [33]. Cases of weed resistance to herbicides due to NTS mechanisms have been increasing (www. weedscience.org). The genetic bases of NTS mechanisms are not fully understood due to the complex interaction of biochemical processes and limited genomic information on weedy species [33][34][35]. In this study, we investigated A. palmeri accessions with differential tolerance to glufosinate.
The genome of A. palmeri is not yet deciphered although the genome and transcriptome of its cultivated relative grain amaranth (Amaranthus hypochondriacus) was completed in 2014 [36]. Also recently, the transcriptome of two weedy species Lolium rigidum [37] and Echinochloa crus-galli [38] were sequenced to identify genes involved in herbicide resistance. Understanding the molecular mechanisms underlying herbicide resistance could be used to mitigate and manage resistance evolution and reduce weed impact on crops. This study assembled the transcriptome sequence of A. palmeri exposed to glufosinate compared to controls to elucidate candidate genes involved in differential tolerance to glufosinate.

Plant materials
Amaranthus palmeri seed samples from 120 fields were collected in Arkansas, USA between 2008 and 2014. The collection of plant samples from the field was done with permission of farm owners, managers, consultants, or Extension Agents. In the process of collecting samples, no endangered species were affected. Inflorescences of at least 10 female plants per field were harvested, dried, threshed, and cleaned for herbicide bioassays in the greenhouse. One accession of interest (08-Lee-C) was collected from a field that was planted with glyphosate-tolerant (Roundup Ready 1 ) soybean in 2008 and glyphosate-tolerant cotton in 2006 and 2007. Although this field had no record of being sprayed with glufosinate, some plants survived exposure to glufosinate (0.55 kg ha -1 ) in the greenhouse. The survivors were grown and allowed to cross-pollinate to produce the first cycle of purified (intercrossed) progeny (C1).
To study the potential survival mechanisms, seeds of 08-Lee-C and the C1 progeny were planted in 4-cm-diameter pots using commercial potting soil mix (Sunshine Mix, Bellevue, WA, USA). Seedlings (100) were grown at one plant per pot in a growth chamber maintained at 32/26˚C day/night temperature with a 16-h photoperiod. Plants were watered daily and fertilized with a water-soluble all-purpose plant food containing 15-30-15% NPK (Miracle-Gro 1 , Marysville, OH, USA), every 2 wk. Fifty plants per accession (9-cm tall) were sprayed with glufosinate (0.55 kg ai ha -1 ) (Liberty 1 , Bayer) mixed with 20 g L -1 ammonium sulfate to identify S and T plants. Susceptible and T plants from the non-treated control were identified by ammonia accumulation assay. Six confirmed S plants from 08-Lee-C and T plants from C1 were used for ammonia accumulation assay, determination of chloroplast-encoded glutamine synthetase (GS2) copy number, and RNA-Seq experiment.

Phenotypic evaluation of A. palmeri response to glufosinate
The response of A. palmeri, collected between 2008 and 2014, was evaluated in the greenhouse. A known herbicide-susceptible accession (SS) was included in each experiment as control [39]. Five-hundred mg of seeds from each field-collected plant were mixed to make a composite seed sample representing each accession. The experiment was conducted twice in a randomized complete block design with two replications. Each replication consisted of one cellular tray (28 X 54 cm) with 50 cells (Redway Feed Garden and Pet Supply, Reedway, CA, USA) filled with a commercial medium (Sunshine Mix, Bellevue, WA). Composite seeds from each accession were planted in each cell and seedlings were thinned to one per cell. Glufosinate was applied at 0.55 kg ha -1 when seedlings were 7.5-9 cm tall using a laboratory sprayer fitted with a flat fan nozzle tip (800067 TeeJet, Spraying Systems Co., Wheaton, IL, USA) delivering 187 L ha -1 at 269 kPa. The herbicide was applied with 20 g L -1 ammonium sulfate. The plants were assessed visually relative to the non-treated control 21 d after treatment (21 DAT) using a scale of 0 to 100, where 0 = no visible injury and 100% = complete desiccation. The number of survivors was recorded. Survivors from glufosinate treatment were grown to produce seed. Data were analyzed using hierarchal clustering in JMP Pro v. 12.

Herbicide resistance profiling of a selected A. palmeri accession
Data from the differential tolerance evaluation were used to select an accession for further study. Accession 08-Lee-C had the most number of survivors with minimum injury. Seeds from 08-Lee-C and SS accessions were planted as described in the previous section. Seedlings (7.5-9 cm tall) were treated with the recommended dose of fomesafen (264 g ha -1 ) (Flexstar 1 , Syngenta), glyphosate (870 g ha -1 ) (Roundup PowerMAX 1 , Monsanto), dicamba (280 g ha -1 ) (Clarity 1 , BASF), and ALS inhibitors pyrithiobac (73 g ha -1 ) (Staple LX 1 , DuPont) and trifloxysulfuron (8 g ha -1 ) (Envoke 1 , Syngenta). The ALS inhibitors were applied with 0.25% non-ionic surfactant (Induce 1 , Helena Chemical Co., Collierville, TN, USA), respectively. Herbicide treatments were applied as described in the previous section. Mortality was evaluated 21 d after treatment. The experiment was conducted in a randomized complete block design as in the previous section. Data were analyzed using ANOVA in JMP Pro v. 12.

Evaluation of tolerance level to glufosinate
A dose-response bioassay was conducted in the greenhouse to determine the tolerance level of 08-Lee-C and C1 to glufosinate. Seeds were sown in 15-cm diameter pots filled with commercial potting soil. Seedlings, 7.5-cm tall, were sprayed with 11 doses of glufosinate using a laboratory sprayer as described in the previous section. The 08-Lee-C and C1 accessions were sprayed with glufosinate at 0.0012 to 0.5940 kg ai ha -1 ; the SS accession was sprayed at 0.0006 to 0.5950 kg ha -1 . Non-treated checks were included for each accession. The herbicide was applied with 20 g L -1 ammonium sulfate. Shoot biomass was harvested 21 DAT, dried at 60˚C for 72 h, and weighed. The experiment was conducted in a randomized complete block design with four replications. Five plants were used per replication (20 plants total) per herbicide concentration.
Data were analyzed using SAS JMP Pro v. 13 in conjunction with SigmaPlot v.13 (Systat Software, Inc., San Jose, CA, USA) for nonlinear regression analysis. The percentage biomass reduction was fitted to a nonlinear, sigmoid, four-parameter logistic regression model defined by y ¼ c þ ½ðd À cÞ=ð1 þ efÀ aðx À bÞgÞ where y represents the biomass reduction expressed as percentage relative to the non-treated control, a is the growth rate, b is the inflection point, c is the lower asymptote, d is the upper asymptote, and x is the glufosinate dose. The herbicide doses that would cause 50% growth reduction (GR 50 ) were estimated using the fitted regression equation.

Ammonia accumulation assay
To identify S and T plants without glufosinate treatment, a leaf disc assay was conducted to measure ammonia accumulation caused by the inhibition of photorespiration by glufosinate [40]. The assay was conducted using 50 non-treated plants each from 08-Lee-C and the C1 progeny. In addition, leaf tissues from three 08-Lee-C plants that were controlled (S) and three C1 plants that survived (T) glufosinate application at the whole plant level were also tested. From each plant, two leaf discs (5-mm diameter) were excised from the youngest, fully expanded leaf of 6.4-cm tall seedlings. One leaf disc was placed per well in a microtiter plate containing 200 μM glufosinate. The plate was sealed with micropore tape and placed on a bench under light for 24 h. The plate was moved to a -80˚C freezer to stop the reaction. After two freeze-thaw cycles, ammonia content was measured in a spectrophotometer (Pharma Spec UV-1700, Shimadzu, Columbia, MD) at 630 nm using a modified method by Molin and Khan [41]. Leaf discs from S plants were expected to have higher ammonia content than those from T plants.

Glutamine synthetase (GS2) relative copy number
Leaf tissues were harvested from confirmed S and T plants (three each) without glufosinate treatment and stored at -80˚C until processing. Leaf tissues were harvested also from plants treated with glufosinate, 24 HAT. Upon evaluation of plant response 21 DAT, leaf tissues from three S and T plants were used also to determine the relative copy number of GS2 gene. Genomic DNA was extracted using the modified hexadecyltrimethylammonium bromide (CTAB) method [42]. Quantitative real-time PCR (qPCR) was used to determine the genomic copy number of GS2 relative to a housekeeping RNA dead box helicase gene GS2 in A. palmeri. The primer pair GS2-F (5'-ATACGGAGAAGGAAGGCAAAG -3') and GS2-R (5'-TGTGGGT TCCCAAAGTAGTG-3') were designed to amplify a region of the chloroplast GS. RNA dead box helicase gene primers A36-F (5'-TTGGAACTGTCAGAGCAACC-3') and A36-R (5'-GAACCCACTTCCACCAAAAC-3') were used as internal primers to normalize the samples for any differences in DNA quantities. Reactions were conducted in three technical replicates, and a negative control consisting of primer pairs with no template was included. An 8-fold serial dilution of genomic DNA samples, ranging from 0.00064 to 50 ng, was used to construct a standard curve. The slope of the standard curve was used to determine amplification efficiency (E). The qPCR reaction efficiency was 97% with an R 2 of 0.9907 and a slope of 3.271 indicating good assay validation. Genomic DNA templates (2 ng) were amplified in a 25-μL reaction containing 12.5-μL Bio-Rad iQ SYBR Green Supermix, 2-μL of primers (1:1 mix of forward and reverse primers), and nuclease-free water. Reaction conditions included 10 min incubation at 94˚C, then 40 cycles of 94˚C for 15 s and 60˚C for 1 min, followed by a melt-curve analysis to confirm single PCR product amplification. Data were analyzed using CFX Manager software (v.1.5). Relative GS2 copy number was calculated as ΔCt = (Ct, A36 − Ct, GS2) according to the method described by Gaines et al [43]. Increase in GS2 copy number was expressed as 2 ΔCt . Results were expressed as the fold increase in GS2 copy number relative to RNA dead box helicase.

RNA-Seq analysis
Sample preparation for RNA-Seq. This experiment used leaf tissues from non-treated and treated, confirmed S and T plants. These were the same plants used for ammonia assay and GS2 copy number determination. Tissues were collected 24 h after glufosinate application for RNA extraction. This collection time was selected to capture herbicide stress adaptation genes and because maximum absorption of glufosinate occurs 24 HAT [28]. Treatments were designated as non-treated S (susceptible without treatment, SWO), non-treated tolerant (tolerant without treatment, TWO), treated susceptible (SWT), and treated tolerant (TWT) plants with three biological replicates. Total RNA was extracted from young leaf tissues of S and T plants using PureLink RNA Mini Kit (Life Technologies, Carlsbad, CA, USA) following the manufacturer's protocol. The extracted RNA was treated with DNase (Invitrogen, Carlsbad, CA, USA) to remove potential genomic DNA contamination, according to the manufacturer's instructions. The samples were then sent to the Clemson University Genomics and Computational Biology Laboratory, South Carolina for sequencing the transcriptome.
Transcriptome sequencing and assembly. Total RNA was normalized, and converted to cDNA using the TruSeq RNA library kit v2.0 (Illumina). Final sequencing products were validated for size on an Agilent Bioanalyzer 2100 (Agilent Technologies, Waldbronn, Germany) and sequenced using a 2x125bp paired-end sequencing module on an Illumina HiSeq 2500 (Illumina). Raw sequence reads were assessed for quality using the FastQC software package (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and preprocessed to remove sequence adapters and low quality bases with the Trimmomatic software [44]. A reference unigene assembly comprehensive of developmental stage, tissue source, and experimental conditions was prepared by concatenating all preprocessed reads and assembling with the Trinity software package [45]. The resulting unigene assembly was filtered for genuine coding sequences (e.g., sequences without internal stop codons or chimeras) with the TransDecoder software, and clustered by identity with the CD-HIT software [46] in an attempt to collapse homologs, but not paralogs, at high identity thresholds.
Differential gene expression. Paired-end reads from each individual were aligned to the de novo transcriptome using the Subread package [47,48]. Samtools was used to convert alignments from sam to bam format, sort, and index [49, 50]. Subread's featureCounts counted the number of reads that aligned to each gene in the transcriptome [47,48]. The final gene counts were loaded into Bioconductor's edgeR package for statistical analysis [51][52][53][54][55]. Variance between samples was visualized by a multidimensional scaling (MDS) plot. Volcano plots were generated for each comparison of samples. The criteria for differential gene expression included a fold-change !2 between compared groups and statistical significance at P 0.05 [56]. Expression differences were compared between non-treated T and non-treated S (TWO vs SWO), treated T and treated S (TWT vs SWT), treated T and non-treated T (TWT vs TWO), and between treated S and non-treated S (SWT vs SWO).
Transcriptome annotation. The final reference assembly was annotated by blastx and blastp alignment to the non-redundant protein database at NCBI, the UniProt-swissprot database, and the uni-ref database to determine homology to known genes. Protein domains were determined by alignment to the HMMER database (http://hmmer.janelia.org). Signal peptides were determined with the SignalP software [57] and transmembrane regions predicted with tmHMM (http://www.cbs.dtu.dk/index.shtml). Gene ontology terms were derived from the best BLAST match [58] and clustered to determine enrichment using the Agbase tool [59].
The entire dataset was submitted as NCBI BioProject (PRJNA390774), which is a part of the U.S. National Library of Medicine of the National Institutes of Health. The 12 samples that were used to construct the transcriptome, and to run the differential gene expression comparisons were submitted as 12 separate BioSamples, SAMN07260017-SAMN07260028. The trimmed, paired-end .fastq files for each of the 12 samples were submitted to the Sequence Read Archive, SRR5759376-SRR5759387. Finally, the transcriptome was submitted to the Transcriptome Shotgun Assembly Sequence Database. The transcriptome, consisting of 72,780 transcripts, is under TSA submission: SUB2788796.
Sequence analysis of the glutamine synthetase 2 (GS2) gene. Glutamine synthetase 2 (GS2) gene sequences of the T and S biotypes were extracted and assembled from the transcriptome data. A 1296-bp GS2 gene (431 amino acids) from S and T plants was sequenced. The full length GS2 sequences of S and T plants were aligned using Sequencher 5.4.6 and BioEdit software packages to identify amino acid substitutions. Sequence alignment also included GS2 sequences of other Amaranthus species (A. viridis, A. albus, A. spinosus, A. hybridus, A. lividus and A. thunbergii) available in the database at http://www.weedscience.org.
Heat map analysis. Differentially expressed genes associated with abiotic stress response were subsampled and subjected to heat map analysis. Normalized read count averages were calculated to produce biological expression profiles followed by hierarchical clustering to recursively merge expression based on pair-wise distances between non-treated T (TWO), non-treated S (SWO), glufosinate-treated T (TWT), and glufosinate-treated S (SWT) samples. Digital expressions were visualized between rows for normalized read count numbers (minimum & maximum) expression patterns. The expression pattern was generated using GENE-E tool (http://www.broadinstitute.org/cancer/software/GENE-E/).
Selection of candidate non-target genes. Genes that were commonly expressed between treated S and T, and between treated T and non-treated T plants were selected for further evaluation based on their gene ontologies (GO). Genes that were assigned with GO molecular function and biological process related to metabolism and signaling pathways (oxidoreductase activity, nuclear acid binding transcription factor activity, hydrolase activity, transferase activity, transmembrane transporter activity, transferase activity, protein transporter activity, biosynthetic process, small molecule metabolic process, signal transduction, homeostatic process, immune system process, cell wall organization, secondary metabolic process, nitrogen cycle metabolic process) were evaluated based on UniProt and their fold change. Contig assemblies that were consistently upregulated in the treated T (relative to treated S and non-treated T) with a significant P-value in the DESeq analysis were selected, for a total of 49 contigs. A subset of this list was generated based on known gene function. Contigs with predicted annotations related to stress response, signaling, transcription factors, and herbicide metabolism were selected as potential candidate NTS genes involved in glufosinate tolerance.

Candidate gene validation by qRT-PCR
Two A. palmeri populations were treated with glufosinate at 0.37 kg ha -1 using the previously described procedure. Leaf tissues were collected three days before and 24 h after herbicide treatment. Tolerant and susceptible plants from each population were identified three weeks after herbicide treatment based on level of injury. Three biological replicates from the nontreated and treated samples from each biotype within each population were used for the validation experiment. Total RNA was extracted from leaf tissues using PureLink RNA Mini kit (Life Technologies, Carlsbad, CA). Genomic DNA was removed using DNAse I (Thermo Scientific, Waltham, MA). cDNA was generated from 5 μg total RNA using Reverse Transcription System first-strand cDNA synthesis kit (Promega). A 4-fold serial dilution of cDNA samples (1:1, 1:5, 1:25, 1:125) was used to construct a standard curve. Seven of the candidate NTS genes were subjected to real-time quantitative PCR with primers designed using Primer3 tool. Two genes (β-tubulin, RNA helicase) were used as internal controls for normalization of gene expression. Primers had an amplification efficiency of 96 to 110%.
The expression level of 7 candidate NTS genes was measured in 24 plants. Quantitative real-time PCR reactions were conducted in a 12-μL volume containing 6.25 μL of SyberGreen Master Mix, 1 μL of 1:25 diluted cDNA, and 0.5 μL of 10 μM primers (1:1 mix of forward and reverse primers). Amplification was performed in a Bio-Rad MiniOpticon System PCR machine (CFX96, Bio-Rad Laboratories, Inc., Hercules, CA, USA) under the following conditions: 10 min at 94˚C, 40 cycles of 94˚C for 15 s and 60˚C for 1 min, followed by a melt-curve analysis to confirm single PCR product amplification. Data were analyzed using CFX Manager software (v.1.5). Slopes for target and internal control genes were equivalent as observed in amplification plots. Comparative CT method was used to calculate relative expression levels as 2 −ΔΔCt where ΔCT = [CT target gene-geometric mean (CT internal control genes)] and ΔΔCT = [ΔCt tolerant − ΔCt susceptible]. Wilcoxon non-parametric test (α = 0.05) was used to determine statistical difference in gene expression between tolerant and susceptible biotypes.

Differential response of A. palmeri accessions to glufosinate
All field populations represented by the accessions tested were susceptible to glufosinate. The majority were killed 100% by 0.55 kg ha -1 glufosinate, except for some, which had few survivors. The 120 accessions had differential levels of susceptibility to glufosinate, separating into three groups based on mortality and levels of injury of survivors ( Table 1). The first group was composed of the 88 most sensitive accessions. The second group, composed of 28 accessions, were had 94 to 99% mortality with survivors incurring 60-99% injury. The third group consisted of four segregating accessions with 88 to 97% mortality. Survivors from these accessions incurred 30-80% injury. Of the possibly segregating accessions, only survivors from 08-Lee-C were able to produce sufficient seeds to conduct further experiments. Seven percent of 08-Lee-C survived glufosinate treatment, of which 4% of treated plants had <61% injury and produced seeds. Survivors from other recalcitrant accessions were not able to produce enough seeds due to either having high injury (>75%) or low number of survivors, which were all males. Considering that plants growing in the field tend to be more robust than those in the greenhouse, the likelihood of escapes in the field may be higher than that in the greenhouse. Plants in the field also do not receive uniform amounts of herbicide for various reasons such as partial coverage by other plants or differential plant size. In addition, plants maintained in the greenhouse that are well-watered and cultured under warm temperatures grow faster and reach the recommended spraying heights earlier than those growing in the field due to less variation in ambient conditions [60].

Response of 08-Lee-C to other foliar herbicides
The 08-Lee-C accession was susceptible to dicamba and fomesafen, but resistant to glyphosate and ALS inhibitors, which are commonly used herbicides. The mortality of 08-Lee-C was 98% Table 1. Cluster analysis of A. palmeri accessions treated with glufosinate at 0.55 kg ha -1 .

Resistance level to glufosinate
The response of A. palmeri to 11 doses of glufosinate fitted a sigmoidal, logistic function (Fig 1). The glufosinate doses required to reduce growth by 50% (GR 50 ) were 0.076, 0.110, and 0.214 kg ha -1 for SS, 08-Lee-C, and C1 accessions, respectively (Table 3). Based on these GR 50 values, the level of tolerance to glufosinate in 08-Lee-C and and C1 accessions was 1.4-and 2.8-fold, respectively, relative to the susceptible standard (SS). The GR 50 increased 2-fold, from 110 g ha -1 in 08-Lee-C to 214 g glufosinate ha -1 in C1.

Ammonia accumulation in response to glufosinate
Glutamine synthetase, the target site of glufosinate, is a nuclear-coded enzyme that catalyzes the conversion of L-glutamate to L-glutamine by assimilating ammonia in the cytoplasm and plastids, but predominantly in the chloroplast of green tissues [21]. Ammonia accumulation is a direct response to the inhibition of this pathway by glufosinate. The ammonia concentration in S plants was 830 μg g -1 fresh leaf tissue (±60) and was 394 μg g -1 fresh weight (±40) in T plants (Fig 2). The S plants accumulated 2X more ammonia than the T plants, indicating rapid depletion of functional glutamine synthetase as a consequence of glufosinate treatment.

Glutamine synthetase 2 (GS2) relative copy number
The relative GS2 copy number of S and T plants ranged from 1 to 3 (Fig 3). Similar GS2 copies were detected in both S and T plants, indicating that differential tolerance to glufosinate is not due to amplification of the GS2 gene. Transcriptome analysis also revealed that GS2 was not differentially expressed between T and S plants, which indicated that differential tolerance to glufosinate in A. palmeri was not due to changes in expression of the target enzyme.

Glutamine synthetase gene sequence analysis
Glutamine synthetase plays a primary role in plant nitrogen metabolism by catalyzing the conversion of glutamate to glutamine [21]. Glutamine synthetase in higher plants exists in two major isoforms: GS1 in the cytosol and GS2 in the chloroplast/plastids [61]. The cytosolic form (GS1) is the predominant isoform in roots and non-green tissues [62]. The chloroplast form of glutamine synthetase (GS2) is the major isoform in leaves, which is primarily responsible for recycling ammonia during photorespiration and synthesis of glutamine [63]. In our study, two different alleles of GS2 were observed in the T plants and one in the S plants.  37 Gln, Gly 39 Lys, Gln 54 Lys, Asp 56 Glu) in the upstream region were detected in one of the GS2 alleles of the T biotype (Fig 4). The second allele of the T biotype harbored only the Tyr 8-Asn substitution. These nonsynonymous substitutions identified in the A. palmeri T biotype also occur in herbicide-susceptible A. viridis, indicating that these substitutions do not contribute to tolerance to glufosinate. Other nucleotide polymorphisms between T and S plants were synonymous mutations.

Global transcriptional changes in A. palmeri 24 h after glufosinate application
A reference cDNA transcriptome consisting of 72,794 sequences was assembled (Table 4). Treatment samples were similar as indicated in the multidimensional scaling (MDS) plot ( Fig 5). Biological replicates from the same treatment clustered together indicating low bias and variation among treatment samples. One or more GO terms were assigned to 33,516 sequences with 76,455 GO assignments in total for biological process (31.9%), cellular component (10.8%) and molecular function (57.3%) categories.

Background differences between susceptible and tolerant plants.
Pairwise comparison between non-treated S and T plants showed 438 differentially expressed genes, 158 of which were downregulated and 280 were upregulated in the T plants relative to S plants ( Table 5). Genes that were notably more expressed in the T plants relative to the S plants without herbicide treatment included cytochrome P450s (Cyp72A219, Cyp86b, Cyt77A, Cyt71A, Cyt76A, Cyt86A), transporters (ABC transporter), transferases (glycosyltransferase, acylytransferases), antioxidants (glutathione-S-transferase, superoxide dismutase), and genes related to lipid metabolism (esterase lipase).  Amaranthus transcriptomes with and without glufosinate Genes induced by glufosinate application. Relative to the respective non-treated checks, 8154 genes were affected by glufosinate application in the T plants and 6034 genes in the S plants (Table 5, Figs 6 and 7). Comparison between treated T and S plants revealed 567 genes that were more repressed or more induced by glufosinate in the treated T plants. Overall, there  . The x-axis shows the log fold change or relative abundance. The P value (-log base 10) for differential gene expression is plotted on the y axis. Dots in black represent genes that did not achieve significant changes in expression; colored dots on the left indicate genes with significantly downregulated expression and colored dots on the right indicates genes with significantly upregulated expression. https://doi.org/10.1371/journal.pone.0195488.g006 were 210 upregulated and 357 downregulated genes in the treated T relative to the treated S plants (Figs 8 and 9). One hundred-five glufosinate-responsive genes were differentially expressed in both treated T (32 genes) and S plants (73 genes) (Fig 7). In addition, 239 genes that were differentially expressed between treated (52 genes) and non-treated (187 genes) T plants were more highly repressed or more induced in treated T plants than in treated S plants.
Of these 239 genes, the majority were related to biosynthetic process, cellular nitrogen compound and small molecule metabolic processes, response to stress, and oxidoreductase activity (Fig 9). Among the upregulated genes of this 239-gene subset, 91 were induced by glufosinate in T plants, including genes putatively annotated as NAC transcription factor, CYP94a1, and ABC transporter b. The majority of upregulated genes that were differentially expressed between treated T and S plants, and that were also differentially expressed relative to their respective non-treated counterparts, were related to nitrogen compound metabolic processes, oxidoreductase activity, nucleotide binding, and transferase activity. Two genes that were repressed exclusively in the treated T plants relative to the treated S plants were folypolyglutamate synthase and caffeic acid 3-o-methyltransferase (Fig 7). The functional classification of selected differentially expressed genes (DEGs) associated with stress response and herbicide metabolism was examined to investigate the pattern of transcriptome regulation that occurred during glufosinate treatment (Fig 10). Glufosinate treatment triggered the expression of genes related to stress response and xenobiotic detoxification as expected. Some genes associated with photosynthesis, structural stabilization, cell membrane binding, stress response, and detoxification were repressed. Increased expression of non-target site (NTS) genes related to stress response, stress signaling, detoxification, abiotic response, cell structure stabilization, and growth and senescence was observed in treated T plants. Genes that were exclusively induced in treated T plants were annotated as coding for transmembrane protein 45b, heat stress transcription factor B, hypersensitive-induced response protein, cytochrome P450 (Cyp72A219, 94A2, Cyt86b1-like), transcription factor, ethylene-responsive transcription factor, glutathione S-transferase (GST), zinc finger protein constans-like 10, and NAC transcription factor (Fig 8, Tables 5 and 6). These candidate genes likely play a role in the adaptation of A. palmeri to glufosinate and, possibly, also to other herbicides.
Of the 239 genes that were more highly induced or repressed in the treated T plants, 49 were consistently upregulated. These consistently upregulated genes are involved in biosynthetic process, cellular nitrogen compound metabolic process, nucleic acid binding transcription factor activity, oxidoreductase activity, stress response, signal transduction, transferase activity, transmembrane transporter activity, and transport (Table 6 and Fig 11). A subset of 13 glufosinate-inducible genes which are related to detoxification, stress signaling, and transport are candidate genes involved in conferring some tolerance to glufosinate. These include ABA 8'-hydroxylase, ABC transporter, chitin elicitor receptor kinase, cytochrome P450 72A, cytochrome P450 94A, GST, heat stress transcription factor, heat shock protein 83, ethylene response transcription factor, NAC transcription factor, NAC transcription factor 25, NADH ubiquinone oxidoreductase, and nitronate monooxygenase (NMO) ( Table 7). These genes were induced >2-fold in treated T plants relative to treated S plants and non-treated T plants.

Validation of selected genes using qRT-PCR
The expression of 7 candidate NTSR genes was measured in two A. palmeri accessions using quantitative real-time PCR. Gene expression was similar in the non-treated tolerant and susceptible plants. The 7 genes were induced in T and S plants upon glufosinate treatment. Five genes (HSP, NMO, ETF, ABC, NAC) were not differentially expressed between treated tolerant

Fig 8. Biological processes (A) and molecular functions (B) of upregulated genes in treated T relative to treated S plants.
https://doi.org/10.1371/journal.pone.0195488.g008 and treated susceptible plants. GST was differentially expressed in only one of the treated tolerant plants. On the other hand, CYP72A219 was expressed eight times higher in all the treated tolerant plants relative to the susceptible ones (Fig 12).

Differential tolerance to glufosinate among A. palmeri accessions
The differential response to glufosinate among 120 Amaranthus palmeri accessions demonstrates variation in herbicide efficacy on a weed species (Table 1). Several factors affect glufosinate activity even on a single species; these include temperature, light, relative humidity, time of day, plant size/age, and dose [30,[64][65][66][67][68][69]. Variation in environmental conditions was minimal in the greenhouse. Variation in plant factors was minimized by maintaining seedlings of the same size. The impact of time of day on glufosinate activity was eliminated by applying the herbicide at about the same time in each repetition of the experiment. Accessions with a few survivors reflect some heterogeneity within the population, as expected of this highly diverse species. Projecting such diversity to field conditions, where all the factors mentioned above can vary, we expect to see a higher frequency of individuals that would escape weed control activities in the field. Many such escapes could have been subjected to sub-lethal doses as a consequence of plant, environmental, and application variables. Controlling these escaped individuals, or adopting a weed management strategy that controls escaping genotypes, is a critical step in mitigating the accumulation of non-target-site genes that could eventually endow population-level resistance to herbicides [70].

Resistance profiling of the 08-Lee-C recalcitrant accession
Herbicides impose strong abiotic stresses to weeds in crop fields, managed turfgrass, gardens, roadsides, and other chemically managed areas. The evolution of resistance in populations of weedy species is an increasing problem worldwide. The recalcitrant A. palmeri accession, 08-Lee-C, was resistant to glyphosate and ALS-inhibiting herbicides, trifloxysulfuron and pyrithiobac ( Table 2). The occurrence of resistance to multiple herbicides in 08-Lee-C is not surprising because this field was sprayed with glyphosate and ALS inhibitors for several years. This field was planted with glyphosate-tolerant crops for more than three years and had been exposed to ALS inhibitors in the years prior when the grower was planting conventional soybean. Resistance to glyphosate and ALS inhibitors among A. palmeri in Arkansas is widespread [71]. Fomesafen, a PPO herbicide with soil and foliar activity, was first commercialized in the 1960s and had been used by farmers mainly for soybean. The usage of fomesafen, and almost all other herbicides in soybean, dropped when glyphosate-tolerant soybean was introduced in the mid-1990s. Upon the explosion of glyphosate-resistant Amaranthus species, soybean farmers reverted to using fomesafen and its use was expanded to cotton to control glyphosate-resistant Amaranthus species. The farmer of this field, like many others, had been using only glyphosate to control weeds. Although glufosinate had not been used in this field, some A. palmeri individuals had higher tolerance to glufosinate, allowing them to survive a single application. This population was susceptible to the field use rate of

Fig 11. Biological processes (A) and molecular functions (B) of differentially expressed genes that are common in treated T relative to treated S and non-treated T plants.
https://doi.org/10.1371/journal.pone.0195488.g011 glufosinate, with a few individuals that would tend to escape treatment. Being an obligate outcrossing species, A. palmeri exhibits high genetic diversity, which facilitates its tendency to evolve herbicide resistance. Intensive use of glufosinate in this field, in a manner that allows escapes to produce seed, will accelerate the evolution of resistance through accumulation of multiple low-impact tolerance genes as demonstrated already in some species, including Lolium rigidum [72] and A. palmeri [73]. Tolerance traits can accumulate and get fixed in the population as selection pressure continues.

Ammonia accumulation in response to glufosinate
Ammonia accumulation is directly related to glufosinate toxicity. Inhibition of glutamine synthetase and ammonia accumulation triggers a cascade of reactions, including inhibition of ribulose-1,5-bisphosphate carboxylase/oxygenase (RUBISCO) enzyme [74] and photosystem electron flow [75], affecting photosynthesis [30,76] leading to plant death. Ammonia reduces pH gradient across the membrane, which uncouples photophosphorylation [75]. Elevated levels of ammonia accumulated in glufosinate-treated rice and soybean cell cultures [77,78]. In our study, glufosinate-sensitive plants accumulated 2X more ammonia than the T plants ( Fig  2). Similarly, glufosinate-sensitive L. perenne ssp. multiflorum from Oregon accumulated 1.6X more ammonia than the resistant population [25]. Increased ammonia level in S plants after glusofinate treatment is the consequence of rapid depletion of functional glutamine synthetase. Reduced ammonia accumulation in T plants indicates the presence of mechanism(s) that reduce the impact of glufosinate on plant function. Such mechanisms could either be reduced binding affinity of glufosinate by target site modification or, mechanisms external to the herbicide-binding site (NTSM) including detoxification and others. The latter applies to the glufosinate-tolerant A. palmeri plants.

Glutamine synthetase (GS2) copy number
Gene amplification conferring herbicide resistance has been identified in GR weeds such as A. palmeri, A. tuberculatus, K. scoparia, and L. multiflorum [43,79,80]. The GR plants contain multiple copies of EPSPS, the target site of glyphosate, which results in increased production of EPSPS enzyme allowing the plant to function normally despite the presence of glyphosate. This mechanism has not been observed with other herbicide target genes either because it is exclusive to the EPSPS regulatory process, or simply because it has not been investigated in other herbicide target genes. Amplification of GS2, in glufosinate-resistant weeds, is not yet reported. However, a 4-to 11-fold amplification of GS2 in alfalfa cell culture lines resulted in increased GS enzyme production, endowing resistance to glufosinate [81]. In our study, S and T A. palmeri had similar copies (1-3) of GS2 indicating that tolerance to glufosinate was not due to GS2 amplification (Fig 3). This was supported by the fact that GS2 transcripts were not different between S and T plants.

Glutamine synthetase (GS2) gene sequence
A rare individual in a population may harbor a mutation in the herbicide-binding site that can alter the folding structure of the protein, resulting in reduced binding affinity of the herbicide. Resistance to ALS-inhibiting herbicides in weeds, in most cases, is due to mutation(s) in one or more of the binding domains in the ALS enzyme [82]. The higher frequency of SNPs in GS2 of T plants could predispose such individuals in the population to accumulate nonsynonymous nucleotide substitutions. However, genetic polymorphisms may not always get translated to protein polymorphisms, and only certain amino acid mutations will result in herbicide resistance [82]. In the current study, most of the polymorphisms observed in the nucleotide sequence of GS2 in T plants were synonymous. Although seven amino acid substitutions in the upstream region were detected in one of the alleles of the T biotype, these substitutions also occur in glufosinate-sensitive A. viridis and is in a region outside of the substrate-binding domain (Fig 4). A Tyr 8 Asn substitution was detected in the two GS2 alleles of the T biotype. Asparagine and tyrosine are both polar, uncharged amino acids, hence Tyr 8 Asn substitution may not alter the physiological and physicochemical stability of the plastid GS enzyme. The presence of this trivial amino acid substitution in glufosinate-tolerant plants as well as similar GS2 copies as that of the S plants suggests that non-target-site tolerance factors are involved. Transcriptome analysis could inform us on differential tolerance mechanisms. A glufosinateresistant L. perenne ssp. multiflorum from Oregon, USA, which showed similar level of ammonia accumulation to the tolerant A. palmeri in our study, harbors a Asn 171 Asp mutation (GAC TO AAC) in the GS2 gene which confers resistance to glufosinate [25].

Tolerance level to glufosinate in 08-Lee-C and C1
After one cycle of selection, the GR 50 for C1 increased from 1.44-fold to 2.80-fold relative to SS, reflecting increased frequency of tolerant plants in C1 (Table 3). Although a low frequency (<10%) of the plants survived exposure to glufosinate, the increase in GR 50 after one cycle of glufosinate selection is indicative that the population has become less sensitive to glufosinate after one cycle of purifying selection. Because glufosinate had not been used in the field, the frequency of survivors was low, and the tolerance level was only 1.44 to 2.80-fold, it is likely that these plant variants have low-level, non-target site resistance. We are possibly capturing an early phase of herbicide resistance evolution. Should the population continue to be under selection pressure from glufosinate, the probability of the population acquiring additional adaptive alleles and expressing resistance to field use rate of glufosinate would increase.

Transcriptome of A. palmeri and candidate NTS genes involved in glufosinate tolerance
A reference cDNA transcriptome consisting of 72,794 sequences was assembled for A. palmeri (BioProject PRJNA390774, TaxId 107608). The transcriptome of A. hypochondriacus had 57,658 assembled sequences [36]. Our data demonstrated broad effects of glufosinate on several metabolic pathways, as expected of a herbicidal compound, including nitrogen assimilation and metabolism similar to what is reported in Arabidopsis [83]. One of the apparent functional categories to which glufosinate-responsive genes belong is protein families known to participate in metabolism, stress response, and defense, with the majority of these genes potentially associated with abiotic stress response signaling and chemical detoxification pathways (Figs 10 and 11). Stress response genes are inducible by many other herbicides or stress factors. Generally, abiotic stress such as herbicide, salinity, and drought modulates the expression of genes that are involved in signaling cascades and in transcriptional control [84], genes that code for proteins involved in membrane protection [85], and genes that are involved in water and ion uptake and transport [86,87]. These stress-regulated genes are activated to counteract the stress effects, maintain homeostasis, and adapt. Cell membrane receptorkinases, stress signaling genes, detoxification-related genes, and antioxidants were activated upon glufosinate treatment in both S and T plants (Fig 10). Peroxidase and superoxide dismutase, for example, were upregulated to help counteract the oxidative stress caused by lipid peroxidation resulting from glufosinate treatment. This indicates that plants undergo extensive transcriptional adjustment in response to herbicide-induced stress. Activation of herbicidestress-response genes is hypothesized to be initiated by a herbicide sensor, which triggers the activation of regulator genes, which causes a cascade of reactions to either detoxify the herbicide or protect the plant from herbicide-mediated stress [35]. Two genes, folypolyglutamate synthase and caffeic acid 3-o-methyltransferase, were repressed in the treated T relative to S plants, but were not differentially expressed in other pairwise comparisons. Both are involved in one-carbon transfer and phenylpropanoid biosynthesis [88,89]. The phenylpropanoid pathway serves as a rich source of metabolites in plants, especially for lignin biosynthesis and the production of flavonoids, coumarins, hydrocinnamic acid conjugate, cutins and lignins [90,91]. Phenylpropanoids are involved in plant defense, structural support and survival [91,92]. Repression of genes involved in phenylpropanoid biosynthesis indicate less allocation of carbon resources to these plant products in the T plants after glufosinate treatment compared to the S plants, indicating a shift in carbon allocation to other intermediates, which are more critical for survival under herbicide stress.
Of the seven candidate genes subjected to qRT-PCR validation experiment, only two genes, (cytochrome P450 CYP72A219 and GST) were associated with tolerance to glufosinate ( Fig  12). The induction of GST and CYP72A219 suggests that T plants are able to deactivate glufosinate to some extent. Induction of cytochrome P450 and GST indicates possible conversion of the herbicide into a less toxic metabolite. The biochemical role of cytochrome P450-mediated herbicide metabolism has been well established in herbicide-resistant weed species. Plant cytochrome P450s facilitate the detoxification of toxic xenobiotics by catalyzing oxygen-and NADPH-dependent mono-oxygenation reactions which convert herbicide into a more hydrophilic metabolite [103]. Hundreds of P450 genes exist in higher plants. For example, Arabidopsis thaliana and Oryza sativa possess 272 and 458 putative P450 genes, respectively [104]. RNA-Seq transcriptome analysis of L. rigidum identified CYP72A genes to be involved in metabolic resistance to diclofop [37]. Non-target-site ACCase and ALS resistance in Alopecurus myosuroides [95,96], Stellaria media [105], Lolium [96,[106][107][108], Sinapis arvensis [109], Echinochloa phyllopogon [110], and Digitaria sanguinalis [111] were reported previously to be facilitated by cytochrome P450 enzymes. Upregulation of CYP72A and CYP94A was reported in a multiple-herbicide-resistant E. phyllopogon population [112]. Similarly, CYP94A1, a plant cytochrome P450-catalyzing fatty acid omega hydroxylase, was induced by chemical stress in Vicia sativa and by bentazon treatment in soybean [113,114]. Some cytochrome P450 genes in the CYP71A family were also demonstrated to be involved directly in herbicide metabolism in crops, such as O. sativa CYP71A31 and Zea mays CYP71A28 [115,116].
The involvement of glutathione S-transferases (GSTs) in herbicide resistance is reported in several weed species. Glutathione S-transferases are ubiquitous enzymes that catalyze the conjugation of harmful xenobiotics to reduced glutathione, facilitating their metabolism, sequestration or removal [117]. The primary factor for atrazine selectivity in corn is the activity of a soluble GST, which detoxifies atrazine by forming an atrazine-glutathione conjugate [118]. In a recent transcriptome study, increased expression of GST is associated with diclofop resistance in L. rigidum and nicosulfuron tolerance in Z. mays [37,116]. GST also functions as an antioxidant, protecting plants from herbicide-mediated oxidative stress by scavenging reactive oxygen species [119]. Increased expression of glutathione transferase gene (AmGSTF1) in multiple-resistant A. myosuroides led to accumulation of flavonoids which protects the plant from herbicide injury [120]. It has been reported that glutathione transferase orchestrate tolerance to abiotic stress through their ability to regulate redox signaling pathways that activate defense genes [121]. In glufosinate-tolerant A. palmeri, GST is possibly involved in converting glufosinate into a less toxic metabolite following possible minimal phase I detoxification by CYP72A219 as well as in protecting the plants against oxidative stress and lipid peroxidation from glufosinate phytoxicity.
Other candidate genes were identified by RNASeq, but were not differentially expressed in the validation experiment, including as the ABC transporter, NAC transcription factor, NMO, HSP, and ethylene transcription factor. Although these might not be involved in conferring some level of tolerance to glufosinate, their involvement in detoxification of toxic xenobiotic compounds and in stress response have been reported. For example, NMO in A. thaliana is associated with detoxification of the allelochemical benzoxazolin [122]. Increased expression of NADH ubiquinone oxidoreductase and induction of P450 genes were involved in resistance to pyriproxyfen insecticide in Bemisia tabaci [123]. Induction of heat shock proteins is associated with drought and oxidative stress in Brassica juncea [124] and A. thaliana [125], respectively. Plant ABC transporters have been associated with the movement of herbicide conjugates [126,127]. Modifications of ABC transporters have been suspected in some cases of weed resistance to glyphosate or paraquat [128,129]. Tolerance of Arabidopsis thaliana to paraquat is endowed by a mutation in the plasma membrane-localized ABC transporter, which resulted in reduced herbicide uptake in plant cells [130]. NAC transcription factors and ethylene-response transcription factors play an important role in the regulation of the transcriptional reprogramming associated with plant stress response such as cold, drought, and salinity [131][132][133][134].
Diversity in gene expression and regulation is an important factor driving herbicide resistance evolution [135]. As gene expression regulation also involves post-transcriptional and post-translational controls, protein expression of the identified genes may need to be investigated. Because of genetic diversity, plants have the potential to overcome herbicide stress through a concerted action of multiple genes. Plants with low-level tolerance showed greater induction of abiotic stress-protection-and detoxification-related genes than S plants. Thus, survival from glufosinate treatment is facilitated by stress-protection/stress-adaptation genes. Differential expression of stress-protection genes in a population can enable some individuals to survive herbicide application. Tolerance-related genes can get fixed in a population upon the exertion of sustained selection pressure. Selection pressure coupled with genetic diversity drives evolutionary processes leading to herbicide resistance.
A low frequency of plants with reduced sensitivity to glufosinate was observed across a twostage screening of survivors from a recalcitrant, segregating population. The surviving individuals were described as tolerant because there was no record of their exposure to glufosinate in this field and the plants can be controlled 100% with split application of glufosinate. However, a three-year record, which is what growers can provide generally, is insufficient to assert that the gene regulation demonstrated here is truly ancestral and that it would have been expressed in the total absence of glufosinate selection. Amaranthus palmeri is a prolific dioecious plant. Seed could have been transported to the field from elsewhere where glufosinate has been used. Similarly, pollen from other fields with a history of glufosinate use, could have blown through the field and fertilized female plants. Differential gene expression of CYP72A219 and GST and the presence of surviving progeny in intentionally intercrossed populations show that genetic mechanisms exist for the evolution of low-level, or potentially incipient, resistance to glufosinate in some A. palmeri populations. Vigilance will be required to detect elevated glufosinate tolerance, especially because it is likely multi-genic, impossible to eliminate from populations, and eventually might confer resistance to more than one class of herbicide chemistries and abiotic stresses.