Comparative transcriptome analysis provides insights into seed germination in herbicide resistance Polypogon fugax

Herbicide resistant mutations are predicted to exhibit tness cost under herbicide-free conditions. Asia minor bluegrass (Polypogon fugax) is a common weed species in the winter crops. Our previous study established a P. fugax accession (LR) resistant to acetyl-CoA carboxylase (ACCase) inhibiting herbicides. Besides, LR also exhibited cost, like lower germination, relative to the sensitive plants (LS). But little is known about the gene expression of in


Background
Herbicide resistant mutations were considered as pre-existed or arise spontaneously with very low rate within weeds populations [1][2][3]. With frequent herbicide application, herbicide resistant mutations are rapidly selected and enriched in weed populations [4]. Without herbicide stress, herbicide resistant mutations are predict to exhibit an adaptation cost ( tness cost) due to the impaired enzyme ability or altered feedback inhibition resulting in insu cient or excessive product biosynthesis [5][6][7]. Many herbicide resistant weeds exhibited altered phenotypes, including seed germination ability, growth cycle, biomass production, competitive ability and reproduction [8,9]. For instance, the acetolactate synthase (ALS) and ACCase herbicides multiple-resistance Lolium rigidum biotypes showed lower vigor and competitive ability against wheat than that of the sensitive biotype [10]. The Alopecurus myosuroides ACCase gene mutant showed a reduction of biomass, height and seed production [11]. Two Hordeum glaucum ACCase gene mutants exhibited lower growth rate, biomass and seed production [12].
Asia minor bluegrass (Polypogon fugax) is a common winter weed species distributed across China and other Asian countries. This annual grass has an extended period of emergence from early November to late December. Its life cycle is highly close to several winter crops, including wheat (Triticum aestivum L.), rapeseed (Brassica napus L.) and some vegetables [13,14]. P. fugax is a competitive weed, especially in moist soil and has become an increasing problem in wheat or rapeseed elds in rotation with rice [15]. In our previous study, a P. fugax accession (LR) collected from Sichuan province in China was highly resistant to ACCase inhibiting herbicides. The resistance mechanism was conferred by an Ile-2041-Asn substitution in ACCase gene [16]. In addition, this herbicide resistance plants exhibited tness cost, including lower germination and emergence potential, earlier owering and tiller and panicle emergence and seed shedding, relaively to the sensitive plants [17,18]. A transcriptome study of the P.fugax resistant and sensitive plants were performed at the owering stage, which provide a genomic resource for understanding the molecular basis of early owering. The study identi ed twelve genes exhibited different expression pattern in herbicide resistant plants and sensitive plants [19]. But the molecular basis of lower germination ability in the herbicide resistance plants was still unclear. Whether the altered phenotypes were a consequence of pleiotropic effects of resistance mutations was di cult to gure out.
The tness cost occurred in herbicide resistance mutations are not universal, which depends on the particular mutation [11,[20][21][22], genetic background [23], and environmental conditions [24,25]. Thus, the sensitive plants used in this study were collected from the same habitat of resistant plants to minimize the genetic background variance. The resistance and sensitive plants reproduced for three generations under the same environmental condition to minimize the in uence of environmental conditions. These limitations allowed us to focus on the effect of ACCase mutation on tness cost.
Seed germination and dormancy are important components of plant tness with complex regulation network, including energy production, protein metabolism, phytohormone control and transcription factors regulation. With the development of genomics, transcriptomic analysis was one of the most effective and economic approaches to investigate the difference of germination process between the tested samples. A transcriptome analysis of two maize inbred lines with different germination ability identi ed abundant genes associated with plant hormones, nutrient reservoir, amino acid metabolism, and ribsome in driving seed germination. Genes related to plant hormones and amino acid metabolism showed higher expression level in the maize with high germination ability [26]. A comparative transcriptome analysis between two sheepgrass (Leymus chinensis) germplasm with various germination percentages at three developing stages indicated that genes involved in starch and sucrose metabolism, phenylpropanoid biosynthesis, plant hormone signal transduction, amino sugar and nucleotide sugar metabolism, and photosynthesis were signi cantly changed [27]. The aim of this study is to investigate the differential expression pro le of seed germination process in the herbicide resistant and sensitive plants, and try to gure out the relevant biological pathways or function genes contribute to germination delay in herbicide resistant plants.

Results And Discussion
Characterization of seed germination in LS and LR The maximum germination (gMAX) of LR was about 52% while that of LS reached 94% (Fig 1a), consistent with the previous study results [17], suggesting that the low germination character was steadily inherited in LR. According to the germination curve, the rapid growth period was from 3d to 8d after incubated (DAI). Seed germination is a process with complex physiological, biochemical and molecular biological basis [28]. With the absorption of water, the storage macromolecules such as starch, proteins and lipids were decomposed into small and soluble molecules that can be easily utilized and transformed to support the germination [29]. The content of soluble sugar in the 0 DAI seeds of LS was higher than LR. After incubation, the soluble sugar content deceased and then dramatically increased at 6 DAI in LS and LR seeds. The soluble sugar content of LR seeds was always lower than that of LS seeds until 8 DAI (Fig. 1b). Similarly, the content of soluble protein showed down rst and then up, and the content in LR was lower than that in LS during the germination process (Fig. 1c). These results indicated that more soluble sugar and protein could directly participate carbon oxidation and enzyme catalysis to supply energy in LS in comparison with LR during germination. Combined with the dynamic change of germination rate and soluble sugar and protein content, we chose 3 DAI and 6 DAI seeds for transcriptome analysis to compare the difference of germination process between LS and LR.
Hormones regulate the seed germination of LS and LR Gibberellin (GA) and abscisic acid (ABA) antagonistically regulate the transition of germination and dormancy, that GA promote germination and ABA promote dormancy [30] To investigate whether the germination delay of LR was related to the hormonal regulation, the seeds of LS and LR were treated with exogenous GA and ABA and their synthesis inhibitors, uridone (FL) and paclobutrazol (PA), to observed the gMAX values (Fig. 1d). The application of GA signi cantly evaluated the gMAX of LR from 52% to 78%, suggesting that exogenous GA could promote the germination of LR. The application of ABA had no signi cant effect on the gMAX of LS but decrease that of LR to 17%, suggesting that the inhibitory by ABA was stronger in LR than LS. PA and FL inhibits the GA and ABA biosynthesis had been proved in many plant species [31]. The application of PA decreased the gMAX of LS to 73% and that of LR decreased to 25%, suggesting that the GA biosynthesis in LR was more strongly inhibited by PA. The application of FL dramatically increased the gMAX of LR to 94%, suggesting that the dormancy of LR completely removed when ABA biosynthesis was inhibited.
de novo assembly of P.fugx reference transcriptome Two P. fugax accessions LS and LR were incubated for 3d and 6d, respectively, and three biological replicates were designed in each accession at two time points. Thus twelve samples were used to conduct the RNA-seq analysis. The number of total clean reads of each samples were shown in Table 1. A total of 476,220 transcripts and 149,330 genes were de novo assembled, and the distribution of length interval was shown in Fig. 2a.

Identi cation of differentially expressed genes (DEGs)
To determine the molecular basis involved in the germination process of LS and LR, two comparison settings, LS_6d vs LS_3d and LR_6d vs LR_3d, were analyzed. The biological samples had good reproducibility due to the high correlation values (Fig. 3c). The distribution of the DEGs in LS and LR were directly exhibited in volcano plots ( Fig. 3a and 3b). A total of 11,856 DEGs were identi ed and 8,936 were up-regulated and 2,920 were down-regulated in the LS comparison setting. A total of 23,123 DEGs were identi ed and 20,055 were up-regulated and 3,068 were down-regulated in the LR comparison setting. A venn diagram showed that 7,165 up-regulated DEGs and 1,167 down-regulated DEGs were common in LS and LR comparison settings (Fig. 3d).

GO and KEGG enrichment of DEGs
To understand the difference in biological function and processes related to germination between LS and LR from 3 DAI to 6 DAI, all the DEGs were enriched to GO and KEGG database. The up-regulated DEGs in LS were enriched to 109 GO terms and the down-regulated DEGs were enriched to 9 GO terms (Additional le 5). The up-regulated DEGs in LR were enriched to 224 GO terms and the down-regulated DEGs were enriched to 40 GO terms (Additional le 6). The up-regulated DEGs in LS and LR were highly enriched in "metabolic process", "single-organism process", and "catalytic activity". With regard to the down-regulated DEGs, the enriched GO terms were much less in LS than LR. These results indicated that most of the upregulated DEGs in both LS and LR had similar biological function. All the up-regulated DEGs in LS and LR were enriched to KEGG database and the top 20 enriched pathways were shown in Fig. 4a and 4b. The majority of up-regulated DEGs were related to lipid metabolism, carbohydrate metabolism, amino acid metabolism and secondary metabolites biosynthesis in LS and LR seed during germination. Here, we selected some pathways to investigate the difference of regulation network in LS and LR (Table 3).

DEGs related to hormones biosynthesis and signal transduction
Seeds germination and dormancy processes are regulated by diverse endogenous hormones. According to the results of exogenous hormones regulation, GA and ABA affected the seed germination of LS and LR at different levels. The genes (Cluster-37472.24032 and Cluster-37472.23672) coding two ent-kaurene synthases (KS) and gene (Cluster-37472.70387) encoding gibberellin 3beta-hydroxylase (GA3ox) were upregulated in both LS and LR comparison settings, and the increasing was much higher in LS than that in LR. The rst step in GA biosynthetic pathway is transformed geranylgeranyl pyrophosphate to entkaurene catalyzed by ent-kaurene synthetases. The KS de cient mutants of Arabidopsis showed strong seed dormancy and recovered germinate with exogenous GA treatment [32,33]. GA3ox catalyzes the nal biosynthetic step to produce bioactive GAs. Two Arabidopsis genes, GA4 and GA4H, encoding GA3ox were expressed during seed germination [34]. These GA biosynthesis related genes with higher expression in LS suggested that GA biosynthesis was more active in LS than LR. Consistent with this result, the inhibition of GA biosynthesis by PA was more effective in LR than LS. The gene (Cluster-51396.0) encoding phytochrome-interacting factor 4 (PIF4) was down-regulated in LS but up-regulated in LR during germination. PIF4 regulated the gibberellin-signaling pathway via DELLA proteins, which blocked the GA signal pathway and resulting in germination delay [35]. Thus, the increased expression of PIF4 might repress seed germination via inhibiting the GA signal transduction in LR.
ABA is the major hormone that inducing seed dormancy. The gene (Cluster-37472.23110) coding 9-cisepoxycarotenoid dioxygenase (NCED) was down-regulated in LS but up-regulated in LR during germination. The expression of gene (Cluster-37472.23992) coding phytoene synthase (PSY) and gene (Cluster-37472.76697) coding lycopene epsilon-cyclase (LcyE) were both elevated in LS and LR during germination, and the increasing was higher in LR than LS. The cleavage of 9-cis-epoxycarotenoids catalyzed by NCED is the key regulatory step of ABA biosynthesis. The Arabidopsis mutants of NCED genes, Atnced6 and Atnced9, reduced ABA content level in seeds [36]. PSY is the limiting step of carotenoids synthesis, the upstream of ABA biosynthesis. The PSY gene overexpressed Arabidopsis mutant exhibited delayed germination, and the degree of delay was positively associated with the increased levels of carotenoids and ABA [37]. The gene (Cluster-37472.6593) encoding serine/threonineprotein kinase SRK2 (SNRK2) and gene (Cluster-48355.0) coding ABA responsive element binding factor (ABF) were up-regulated in LS and LR, and increasing were higher in LR than LS. It has established in rice and Arabidopsis that SNRK2 is activated by ABA and phosphorylate ABF, which is important for the activation of ABA signal transduction [38,39]. Taken together, these genes associated with ABA biosynthesis and signal transduction pathways exhibited higher expression level in LR than LS, which might explain why the LR seeds were more sensitive to ABA and FL could restore the germination ability of LR.
DEGs related to fatty acid metabolism ACCase catalyzes the carboxylation of acetyl-CoA to malonyl-CoA, acting as the initial step of fatty acid biosynthesis [40,41]. ACCase has three subunits: biotin carboxylase carrier protein, biotin carboxylase, and carboxyltransferase domains [42]. Molecular and biochemical studies have proved that the carboxyltransferase domain is the target site of ACCase-inhibiting herbicides [43]. LR plants had an amino acid substitution in carboxyltransferase domain that caused the resistance to ACCase-inhibiting herbicide [16]. Unfortunately, the transcriptome data did not identify any DEGs annotated as ACCase carboxyltransferase domain in LS or LR. But gene (Cluster-39490.0) coding ACCase biotin carboxylase (ACACA) was identi ed to up-regulated with the higher increasing level in LS than LR. ACCase is the ratelimiting enzyme in fatty acid metabolism, the different regulation level between LS and LR might in uence the downstream reactions as well. The gene (Cluster-31226.0) coding 3-oxoacyl-[acyl-carrier protein] reductase (FabG), gene (Cluster-37472.4287) coding long-chain acyl-CoA synthetase (ACSL) showed higher regulation level in LS compared with LR during germination. These enzymes were respectively involved in the elongation cycle and termination of fatty acid synthesis. Besides, the genes involved in fatty acid degradation were also up-regulated higher in LS than LR, such as the gene (Cluster-33548.0) coding acetyl-CoA acyltransferase 1 (ACAA1), gene (Cluster-37472.85335) coding acyl-CoA dehydrogenase (ACADM), gene (Cluster-30044.0) coding acetyl-CoA C-acetyltransferase (ACAT). The fatty acid forms triacylglycerol (TAG) that degrades to provide carbon and energy during germination and early seedling growth. Some ACCase herbicide resistant weeds also showed higher level of dormancy, absence of germination in dark conditions, and delayed seed germination [44][45][46]. We assumed that the mutated ACCase might alter normal fatty acid metabolism and impact the germination.

DEGs related to carbohydrate metabolism
Acetyl-CoA, acts as the precursor of fatty acid biosynthesis, is also the product of glycolysis and TCA. According to the dynamic change of soluble sugar and protein content, more soluble sugar and protein could directly participate carbon oxidation and enzyme catalysis in LS than LR during germination. Thus, the genes involved in carbohydrate metabolism were selected to compare the differential expression pattern. The gene (Cluster-37472.6545) coding glucose-6-phosphate isomerase (GPI) gene (Cluster-37472.71503) coding fructose-bisphosphate aldolase (FBA), gene (Cluster-35171.2) coding phosphoenolpyruvate carboxykinase (PCKA), gene (Cluster-40718.0 ) coding pyruvate dehydrogenase E1 component alpha subunit (PDHA), genes (Cluster-37472.39187, Cluster-37472.47005) coding alcohol dehydrogenase (ADH), were up-regulated in LS and LR comparison settings, and the increasing levels were higher in LS than LR. All of these genes were involved in glycolysis, which produces energy by utilizing endosaccharides [47]. PCKA catabolizes the storage lipid and protein to produce soluble sugar via glycolysis. The PCKA de cient mutants of Arabidopsis and tomato were observed growth suppression of germinated seedlings [48]. succinate dehydrogenase (ubiquinone) iron-sulfur subunit (SDHB), and gene (Cluster-41938.0) coding pyruvate carboxylase (PYC) were up-regulated in LS and LR at 6 DAI compared with 3 DAI, and the increasing levels were higher in LS than LR. All of these genes were involved in TCA cycle, which produced high-energy phosphate compounds as the main source of cellular energy [49]. MDH catalyze the interconversion of malate and oxaloacetate in TCA cycle. The Arabidopsis MDH gene mutants showed slowly germination rate, higher content of free amino acids, different sugar levels, and lower content of 2-oxoglutarate [50]. The seeds of the CS de cient Arabidopsis were dormant and did not metabolize triacyglycerol [51].
One of the major changes during germination is a rapid increase in respiration, which involves glycolysis, oxidative pentose phosphate pathway, TCA cycle and oxidative phosphorylation. All of the DEGs mentioned above were related to TCA cycle, glycolysis or pentose phosphate pathway that producing energy during seed germination. These genes were up-regulated in LS and LR during germination and had higher expression in LS, indicating that more energy were supplied to promote germination in LS in comparison with LR.

Validation of DEGs by qRT-PCR
A total of sixteen DEGs were randomly selected to verify the accuracy and reproducibility of the transcriptome results by qRT-PCR (Fig. 5). The correlation between transcriptome results (FC) and qRT-PCR results (2 -ΔΔCt ) were calculated using log 2 fold variation measurements to produce a scatter plot.
The results showed that the expression pro les of these DEGs were consistent with the transcriptome results, with relative R 2 = 0.8751 and 0.7376 in LS and LR, respectively.

Conclusion
In this study, we investigated the expression pro le during seed germination in two P.fugax accessions at two time points. GA and FL promote the germination of LR plants, while ABA and PA showed greater inhibiting effect on germination of LR plants. Consistent with this, the transcriptome results identi ed four GA signal related genes with higher expression in LS and ve ABA signal related genes with higher expression in LR during germination. LS showed higher content of soluble sugar and protein than LR during germination. The expression of the genes related to carbohydrate metabolism and fatty acid metabolism were highly increased in LS than LR, suggesting that more energy was supplied during seed germination in LS (Fig. 6). This study provided novel insight into the gene expression pro le of seed germination in P. fugax and identi ed abundant function genes involved in germination delay in herbicide resistant plants. It should stress that this study could not summarize the effect of herbicide resistance mutation on seed germination. Because the differential expression pro le of seed germination was described only in two accessions. The pathways or function genes identi ed in this study should be veri ed in other resistant plants with germination delay.

Plant materials and growth condition
The P.fugax accessions were originally collected in 2012 in or nearby a grower's eld in Qingshen County, Sichuan Province, China, and identi ed to be resistant and sensitive to ACCase inhibiting herbicides, respectively. No speci c permission is required for the P. fugax accessions collection, and the collected location are not privately owned or natural protected areas. The plant materials were identi ed by Prof. The peeled seeds of LS and LR were surface sterilized with 70% ethanol for 1 min and 10% sodium hypochlorite for 10 min and then washed with sterile water for three times. Sterile seeds placed on the 9 cm petri dishes with one layer of lter paper that moistened with sterile water. The seeds were incubated at 20/10 ℃ day/night cycles with a 12 h photoperiod in the growth chamber. The seed was considered as germinated when the radicle emergence to the seed length. The germinated seeds were observed by 40× microscope (Olympus) at 1d, 2d, 3d, 5d, 6d, 8d, 9d, 11d after incubated. The non-germinated seeds were tested the viability by bromothymol blue and the dead seeds were excluded from calculation. The germination curve was built based on the percentage of germination. Each sample had ve dishes (30 seeds per dish) as replications.
0.1 g sterile seeds from each sample were collected after incubated for 0d, 4d, 6d and 8d, immediately frozen with liquid nitrogen and stored at -80℃. All the samples were homogenized in 1 ml cold distilled water and centrifuged at 12,000 rpm for 5 min and then collected the supernatants for measurement. The soluble protein content was measured as Bradford methods [53]. The soluble sugar concentration was determined by anthrone method [54]. Three replicates were analyzed for each treatment.

Effect of exogenous hormones on seed germination of AS and AR
The peeled seeds of LS and LR were surface sterilized and placed on the petri dishes with one layer of lter paper that moistened with 5 ml waters or tested solutions, including 100μM gibberellin 4+7 (GA), 1μM abscisic acid (ABA), 100μM uridone (FL) and 1μM paclobutrazol (PA). The seeds were incubated under the same condition as mentioned above. The maximum germination (gMAX) was calculated the number of the germinated seeds divided with the total number of the viable seeds. Each treatment had ve dishes (30 seeds per dish) as replications.

RNA isolation and transcriptome sequencing
The peeled seeds of LS and LR were incubated for 3d and 6d and collected for RNA extraction. Each sample had three biological replications named as LS_3d_1, LS_3d_2, LS_3d_3, LR_3d_1, LR_3d_2, LR_3d_3, LS_6d_1, LS_6d_2, LS_6d_3, LR_6d_1, LR_6d_2, LR_6d_3. RNA extraction and RNA-seq was performed by Novogene Corporation (Beijing, China). RNA degradation and contamination was monitored on 1% agarose gels. RNA purity was checked using the NanoPhotometer® spectrophotometer (IMPLEN, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). A total amount of 1.5 µg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA) following manufacturer's recommendations.
Transcriptome assembly and functional classi cation Raw reads of fastq format were rstly processed through in-house perl scripts. In this step, clean reads were obtained by removing reads containing adapter, reads containing ploy-N and low quality reads from raw data. At the same time, Q20, Q30, GC-content and sequence duplication level of the clean data were calculated. All the downstream analysis were based on clean data with high quality. All clean reads of the twelve libraries were obtained using Trinity software to assemble the full-length transcript sequences [55].
Gene function was annotated based on the following databases: NR (NCBI non-redundant protein sequences) NT (NCBI non-redundant nucleotide sequences) Pfam (Protein family) KOG/COG (Clusters of Orthologous Groups of proteins) Swiss-Prot (A manually annotated and reviewed protein sequence database) KO (KEGG Ortholog database) GO (Gene Ontology). The transcription factors were identi ed by iTAK program [56].

Differential gene expression and enrichment analysis
Gene expression levels were calculated using the fragments per kilobase per million fragments (FPKM) method [57]. Differential expression analysis was conducted at 6d compared to 3d in LS and LR using the DEGseq R package [58]. Q-value < 0.05 and | log2 (fold change) | > 1 as the threshold for signi cantly differential expression. All the DEGs were enriched to GO and KEGG database by GOseq and KOBAS software [59, 60]. The signi cant enriched GO terms and KEGG pathways were ltered by padj values < 0.05.

Quantitative real-time PCR validation
The quantitative real-time PCR (qRT-PCR) was performed to verify the DEGs identi ed from the transcriptome results. Total RNA was extracted from the LS and LR seeds under the same condition as the RNA-Seq samples using RNAprep Pure Plant Extraction Kit (Tiangen, China). Then, reverse transcription was performed with 1 μg of each RNA sample using cDNA Synthesis kit (TAKARA, China).
qRT-PCR was performed with iTaq Universal SYBR Green Supermix (Bio-Rad, USA) on QuantStudio TM 1 real-time PCR System and with the following thermal cycle conditions: denaturing at 95°C for 30 s followed by 40 cycles of 95°C for 5 s and 60°C for 1 min. The RT-PCR primers were designed based on the coding sequences (CDS) of the tested genes (Additional le 1.

Consent for publication
Not applicable.
Availability of data and materials.
The raw reads have been deposited in the NCBI Sequence Read Archive (SRA) database (BioProject: PRJNA625626).

Competing interests
The authors declare that they have no competing interests.    Values in a column followed by the different letters means signi cant differences (P ≤ 0.05).