Skip to main content
  • Research article
  • Open access
  • Published:

Transcriptomics analysis of the flowering regulatory genes involved in the herbicide resistance of Asia minor bluegrass (Polypogon fugax)

Abstract

Background

Asia minor bluegrass (Polypogon fugax, P. fugax), a weed that is both distributed across China and associated with winter crops, has evolved resistance to acetyl-CoA carboxylase (ACCase) herbicides, but the resistance mechanism remains unclear. The goal of this study was to analyze the transcriptome between resistant and sensitive populations of P. fugax at the flowering stage.

Results

Populations resistant and susceptible to clodinafop-propargyl showed distinct transcriptome profiles. A total of 206,041 unigenes were identified; 165,901 unique sequences were annotated using BLASTX alignment databases. Among them, 5904 unigenes were classified into 58 transcription factor families. Nine families were related to the regulation of plant growth and development and to stress responses. Twelve unigenes were differentially expressed between the clodinafop-propargyl-sensitive and clodinafop-propargyl-resistant populations at the early flowering stage; among those unigenes, three belonged to the ABI3VP1, BHLH, and GRAS families, while the remaining nine belonged to the MADS family. Compared with the clodinafop-propargyl-sensitive plants, the resistant plants exhibited different expression pattern of these 12 unigenes.

Conclusion

This study identified differentially expressed unigenes related to ACCase-resistant P. fugax and thus provides a genomic resource for understanding the molecular basis of early flowering.

Background

Common annual Asia minor bluegrass (Polypogon fugax) is a weed that is both distributed across China and associated with winter crops. This weed has evolved resistance to clodinafop-propargyl, an acetyl-CoA carboxylase (ACCase) herbicide [1]. The mechanisms of herbicide resistance have been intensively studied in the past twenty years and at least three have been identified: 1) target site change; 2) closure or translocation of herbicides; and 3) alteration in the rate of herbicide metabolism [2,3,4]. Nevertheless, these three mechanisms alone often fail to explain the development of herbicide resistance from an evolutionary and ecological perspective [5]. Herbicide resistance will increase the fitness of resistant individuals and hence their ability to produce the next generation. Identifying which biological characteristics play a major role on fitness and interactions with environmental factors is essential for predicting herbicide resistance.

It is generally believed that the initial occurrence of major resistance genes in weed populations is the main factor that influences the dynamic evolution of a resistance under herbicide selection [6]. In Lolium rigidum, carrying the common Leu-1781 in ACCase affects the fitness of resistant mutants, and the germination rate of the resistant biotype is lower than that of the sensitive biotype [7]. To palliate this lower germination rate, the mutant Setaria viridis produces more seeds than the sensitive population [8]. In addition, herbicide-resistant Setaria flowers earlier than the susceptible population, with more tillers and panicle, and with lighter seeds [9]. These mechanisms contribute to a more important spread of the resistant populations compared with susceptible ones.

Stress-induced flowering has recently received increased attention. Early flowering and seed setting of resistant plants allow the resistant seeds better access to resources [9]. Photoperiodic flowering and vernalization have been well characterized, and the regulatory mechanisms are well known [10,11,12]. In addition, flowering physiologists had reported that plants tend to flower when grown under unsuitable conditions, indicating that stress is a flower-inducing factor [13,14,15,16]. Stress-induced flowering is now considered as the third category of flowering responses alongside regulated autonomous flowering and environment-induced flowering [17].

Nevertheless, the mechanisms underlying early flowering remain poorly understood, particularly among herbicide-resistant weeds. Stress-induced flowering changes the life cycle and might alter fitness. In addition, stress adaptation extends to the evolution of the flowering characteristics [17]. Transcription factors are proteins regulating gene expression and specific transcription factors selectively regulate the transcriptional expression of specific genes. Therefore, in the present study, we aimed to investigate the transcription factors that regulate plant flowering in order to elucidate the relationship between early flowering and selection pressure (herbicide application) of ACCase-resistant P. fugax, and to identify candidate genes responsible for early flowering in resistant plants.

Methods

Plants

The seeds of a putative resistant population of P. fugax (RP population) were collected from Qingsheng County (29° 54′ 1″ N, 103° 48′ 57″ E), Sichuan Province, China, where clodinafop-propargyl has been used for more than five years and has failed to control P. fugax growth. A sensitive population of P. fugax (susceptible plant (SP) population) was sampled from a non-cultivated area in Xichang City, Sichuan Province (27° 50′ 56″ N, 102° 15′ 53″ E). The plants were collected without permissions being sought for the nature of scientific research according to the law of the People’s Republic of China.

Because gene expression can differ due to genetic background, genetically homogenized plant material was generated by controlled pairings to narrow the difference. In brief, F1 plants were transplanted to individual 1-L pots in a greenhouse that has an 18/15 °C day/night temperature and a 14-h photoperiod. At the four-tiller stage, the plants were subjected to vegetative propagation: all individual tillers of each plant were separated and transplanted to individual pots. Four clones were therefore obtained. At the three-leaf growth stage, ACCase herbicide was applied. The herbicide sensitivity of each F1 plant was assessed by spraying with clodinafop-propargyl (45 g. active ingredient (a.i.)/ha). Sensitive and resistant F1 plants were then crossed, yielding an F2 population. Visual phenotype rating of the F2 plants was carried out by clodinafop-propargyl selection. The F2 plants with contrasting phenotypes were selected as resistant and sensitive plants. And the F2 generation was used for transcriptome sequencing.

Sample collection

RPs and SPs (n = 18/group) at the seedling stage (3-leaf stage) were selected, and treated with clodinafop-propargyl (R: 45 g.a.i./ha; S: 0.2 g.a.i/ha) (n = 9 plants/group-); the remaining plants were treated with an equal volume of water. After 72 h, the aerial parts were collected (n = 3 plants/group) and immediately cryopreserved in liquid nitrogen. Afterward, samples were taken at the tillering (four tillers, n = 3 plants/group) and flowering (early flowering, n = 3 plants/group) stages in the same manner and stored in liquid nitrogen.

Genomic library construction and sequencing

The total RNA was extracted using TRIzol (Invitrogen Inc., Carlsbad, CA, USA) and DNase I (Takara, Otsu, Japan) in accordance with the manufacturer’s instructions. Magnetic beads with oligo (dT) (Takara, Otsu, Japan) were used to isolate mRNA, and the mRNA was mixed together with fragmentation buffer by a ThermoMixer (Thermo Fisher Scientific, Waltham, MA, USA), breaking the mRNA into short fragments. cDNA was synthesized using the RevertAid First Strand cDNA Synthesis Kit (Fermentas, Thermo Fisher Scientific, Waltham, MA, USA) in accordance with the manufacturer’s instructions. Short fragments were purified and resolved with EB buffer (10 mM Tris·Cl, pH 8.5) to repair their ends by the addition of a single adenine nucleotide. The short fragments were then connected with adaptors (BGI, Beijing, China). Suitable fragments were selected as templates for PCR amplification. The constructed sample library was quantified and characterized using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and an ABI StepOnePlus Real-Time PCR system (Applied Biosystems, Foster City, CA, USA) and then sequenced using a HiSeq 4000 system (Illumina, Inc., San Diego, CA, USA).

Sequencing and de novo assembly

To obtain high-quality clean reads for de novo assembly, the raw reads were generated from transcriptome sequencing in accordance with the following steps. First, the adaptor sequences were removed. Then, reads with more than 5% of unknown nucleotides were removed and reads with more than 50% of low-quality bases (base quality ≤ 20) were discarded. The clean reads that remained were assembled into unigenes using the Trinity software with an optimized K-mer length of 25 for de novo assembly, as previously published [18]. The expression of unigenes was calculated via the reads per kb per million reads (RPKM, ≥ 0.5), which is a general method of quantifying gene expression from RNA sequencing data by normalizing for total read length and the number of sequencing reads [19].

Data analysis

Genes expressed at different levels in RPs and SPs (i.e., differentially expressed genes (DEGs)) were subjected to Gene Ontology (GO) functional analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. We used getorf software (http://emboss.bioinformatics.nl/cgi-bin/emboss/getorf) to identify the open reading frame (ORF) of each unigene. Then, the ORFs were aligned to the transcription factor (TF) domains using hmmsearch software (http://hmmer.org/) [20].

The false discovery rate (FDR) statistical method was used in multiple hypothesis testing to correct the P-value. A smaller FDR and larger ratio indicate a greater difference in the expression level between two samples. In this analysis, we chose samples with an FDR ≤ 0.001 and a ratio greater than 2.

Differential gene expression analysis

Gene expression levels were estimated using RSEM for each sample, as previously described [21]. Clean reads were mapped back onto the assembled transcriptome, and read count for each gene was obtained from the mapping results and normalized as FPKM (fragments per kilobase of transcript per million mapped reads). Differential expression analysis of the two groups was performed using the DESeq 2. The resulting P-values were adjusted using the Benjamini and Hochberg’s approach for controlling the false discovery rate. Genes with an adjusted P-value < 0.05 (fold change ≥ 2) found by DESeq were determined as being differentially expressed. Differentially expressed genes (DEGs) were analyzed by GO and KEGG enrichment analysis.

GO and KEGG enrichment analyses

All DEGs were mapped to terms in the GO database (http://www.geneontology.org//) and the number of genes corresponding to each GO term was calculated. We established a gene list and gene numbers for each GO term and then used a hypergeometric test to identify DEG GO terms whose enrichment was significantly greater than that in the genome background.

The KEGG pathway database contains networks of molecular interactions in cells and variants specific to particular organisms. We used pathway enrichment analysis to identify DEG metabolic pathways or signal transduction pathways whose enrichment was significantly greater than that in the whole genomic background. After multiple corrections, we selected pathways with an FDR value ≤ 0.001 to represent pathways significantly enriched in DEGs.

Real-time PCR

The total RNA was extracted as described above and cDNA was synthesized using an M-MLV Rtase Kit (Thermo Fisher Scientific, Waltham, MA, USA) in accordance with the manufacturer’s instructions. The qRT-PCR mix (25 μl) contained 12.5 μl of SYBR Green Mix (Thermo Fisher Scientific, Waltham, MA, USA), 0.5 μl of each primer (10 μM), 2 μl of cDNA, and 9.5 μl of RNase-free water. The reaction was performed on an ABI 7300 real-time PCR system (Applied Biosystems, Foster City, CA, USA). The qRT-PCR program consisted of 95 °C for 10 min, 40 cycles of 95 °C for 15 s and 60 °C for 45 s, and finally 60 °C for 15 s. EF1 was used as a reference gene for normalization. GraphPad Prism 5 software (GraphPad Software, Inc., La Jolla, CA, USA) was used for data analysis. Expression was calculated in accordance with the 2-ΔΔCt method. Each experiment was repeated at least three times and consisted of three replicates. Primer sequences are listed in Additional file 1: Table S1.

Data access

The raw reads have been deposited in the NCBI Sequence Read Archive (SRA) database (BioProject PRJNA385696, SRP106591). The data are available at https://trace.ncbi.nlm.nih.gov/Traces/sra_sub/sub.cgi?acc=SRP106591&focus=SRP106591&from=list&action=show:STUDY.

Results

Polypogon fugax Transcriptome sequencing and data assembly

Compared with P. fugax SPs, RPs flowered 10-15 days earlier, and their inflorescence was morphologically altered, RPs produced 1.9 times more seed than did SPs Additional file 2: Figure S1). Therefore, the transcriptome of SPs and RPs were compared to explore the potential genes involved in this process.

Strand-specific RNA-Seq was applied to assess RNAs from three pairs of RPs and SPs at the seedling, tillering, and flowering stages with or without clodinafop-propargyl treatment (seedling stage) to comprehensively identify the unigenes associated with herbicide resistance.

The sequencing reads containing low-quality, adaptor-contaminated, or high contents of unknown base (N) reads were removed before downstream analyses. Afterward, 150-base single-end sequence raw reads were subjected to quality control using the Phred scaled quality score. Overall, 1.39 billion raw reads and 1.07 billion clean reads (average clean read ratio of 77.05%) were obtained; 92.8% of the clean reads had a quality score ≥ 30, and 97.7% of the clean reads were quality filtered and matched the Illumina’s quality requirements. Read quality metrics after filtering are shown in Additional file 3: Table S2. De novo assembly of the 150-base reads yielded 206,041 unique sequences ranging from 300 to 3000 nt in length (Table 1) (including 14,166 unigenes with sequences of up to 3000 nt in length) Additional file 3: Table S3). The length distribution of the assembled contigs is shown in Additional file 4: Figure S2. Among the detected 206,041 unigenes, 165,901 unique sequences were annotated based on BLASTX alignment (E-value <0.00001) searches of seven databases: the NCBI non-redundant (NR), NCBI non-redundant nucleotide (NT), Swiss-Prot protein, KEGG, Cluster of Orthologous Groups of proteins (COG), InterPro protein, and GO databases (Additional file 3: Table S4). The 153,591 unique sequences were annotated by reference to the NR database, and then compared to those encoded in the genomes of all grass (Poaceae) species whose genome is fully sequenced, i.e., Brachypodium distachyon, Hordeum vulgare, Aegilops tauschii, and Triticum urartu. (Additional file 5: Figure S3).

Table 1 Summary of Polypogon fugax transcriptome sequencing, assembly, and annotation

Annotation of assembled unigenes

To further examine the integrity and effectiveness of the annotation process, the number of unigenes (that have NR matches) with a COG classification was calculated. A total of 74,434 unigenes were identified with a COG classification (Additional file 3: Table S5). Among the 25 COG categories, the cluster of “General function prediction only” had the highest number (21,107, 28.36%), followed by “Transcription” (15,595, 20.95%), and “Function unknown” (14,816, 19.90%). Categories of “Extracellular structures” (86, 0.001%) and “Nuclear structure” (11, 1.48 e−4) had the fewest matching genes (Additional file 6: Figure S4).

GO and KEGG enrichment analyses were used to classify the functions of the predicted P. fugax unigenes. Based on homologous genes, 80,312 sequences (Additional file 3: Table S6) from all unigenes of 36 P. fugax libraries were categorized into 56 GO terms comprising three domains: biological process, cellular component, and molecular function (Fig. 1). Most were categorized in “cellular process”, “metabolic process”, “cell”, and “cell part”. A high percentage of genes were also assigned to “binding”, “catalytic activity”, “organelle”, and “membrane” as well as “biological regulation”, “development process”, “transporter activity”, and “reproductive process” (Fig. 1).

Fig. 1
figure 1

Histogram of GO classification. Based on the homologous genes, 80,312 sequences from all unigenes of 36 P. fugax libraries were categorized into 56 GO terms comprising three domains: biological process, cellular component, and molecular function. The majority of the terms were categorized as “cellular process”, “metabolic process”, “cell”, and “cell part”. A high percentage of genes were also assigned to “binding”, “catalytic activity”, “organelle”, and “membrane”, and as well as “biological regulation”, “development process”, “transporter activity”, and “reproductive process”

There were 121,526 unigenes that mapped onto KEGG pathways (Additional file 3: Table S7). A total of 53,337 annotated unigenes between NR, COG, KEGG, SwissProt, and InterPro databases were identified with Venn diagrams (Fig. 2).

Fig. 2
figure 2

Venn diagram of the NR, COG, KEGG, SwissProt, and InterPro databases. A total of 53,337 annotated unigenes between the NR, COG, KEGG, SwissProt, and InterPro databases were identified by Venn diagram

TF prediction of assembled unigenes

Next, we studied unigenes that encoded TFs. The list of unigenes that encode TFs is shown in Additional file 3: Table S8. We also performed TF family classification, and found that 5904 unigenes were classified into 58 TF families (Fig. 3). Among those families, the MYB family had the highest number (742, 12.57%), followed by the MYB-related family (604, 10.23%) and the AP2-EREBP family (453, 7.67%). Genes of the MADS-box family were also found (131, 0.02% %); these genes are associated with plant development and adversity responses.

Fig. 3
figure 3

TF family classification of unigenes. In total, 5904 unigenes were classified into 58 TF families. Among these families, the MYB family represented the highest number (742, 12.57%), followed by the MYB-related family (604, 10.23%) and the AP2-EREBP family (453, 7.67%). Genes of the MADS-box family, which are associated with plant development and adversity responses, were also identified (131, 0.02%). The X axis represents the number of unigenes. The Y axis represents the TF family. “All-Unigene” indicates that the unigenes were those assembled from the 36 samples of Polypogon fugax

qRT-PCR validation of P. fugax expression data

To verify the gene expression patterns, qRT-PCR analyses were performed on Unigene12462, CL1441.Contig20, CL21112.Contig3, CL4600.Contig9 and CL4600. Contig2, CL12188.Contig2, and EF1 served as candidate reference genes for RT-PCR normalization.

The expression of CL12188.Contig2, Unigene12462, and CL1441.Contig20 was higher in the treated sensitive population when flowering (TFS) than in the untreated sensitive population when flowering (UFS). CL4600.Contig2 and CL21112.Contig3 were higher in the UFS than in the TFS. These results confirmed the reliability of the data (Fig. 4).

Fig. 4
figure 4

qRT-PCR validation of RNA-Seq results. qRT-PCR analysis of six randomly selected genes was conducted to confirm the expression patterns indicated by sequencing. The expression of CL12188.Contig2, Unigene12462, and CL1441.Contig20 was higher in the TFS than in the UFS. CL4600.Contig2 and CL21112.Contig3 was higher in the UFS than in the TFS. These results confirmed the reliability of the data. The error bars represent the standard error of the mean. U: untreated; T: treated; F: flowering stage; S: sensitive population; R: resistant population. *significant difference p>0.05, **significant difference 0.01<P<0.05.

Functional analysis of DEGs

To screen the flowering regulatory genes related to resistance, we analyzed DEGs among the seedling, tillering, and flowering stages under different treatments. The genes at the seedling and tillering stages served as a background to identify the specific DEGs at the flowering stage. Cluster analysis was used to compare DEGs, and the parameter was the log2 ratio of gene expression in the difference comparison scheme. The Euclidean distance (calculation of gene distance) referred to genes differentially expressed among all groups. Inter- and inner-group comparisons were performed by the same methods (Fig. 5a-b). Fifty-eight TF families were ultimately predicted, and nine families were related to the regulation of plant growth and development and stress response. We analyzed the different expression values of the TF families (log2-fold change) of the different populations (resistant vs. sensitive) under the same treatment at the same time, including the MYB (Additional file 3: Table S9), MYB-related (Additional file 3: Table S10), MADS (Additional file 3: Table S11), NAC (Additional file 3: Table S12), mTERF (Additional file 3: Table S13), ABI3VP1 (Additional file 3: Table S14), AP2-EREBP (Additional file 3: Table S15), bHLH (Additional file 3: Table S16), and GRAS (Additional file 3: Table S17) families.

Fig. 5
figure 5

Cluster analysis of DEGs. Cluster analyses of inter (a) and inner (b) group comparisons were analyzed by heat maps. Fifty-eight TF families were ultimately predicted, nine families were related to the regulation of plant growth and development and stress responses. We analyzed the different expression values of the TF families (log2-fold change) of the different populations (resistant vs. sensitive) under the same treatment at the same time, including the MYB family, MYB-related family, MADS family, NAC family, mTERF family, ABI3VP1 family, AP2-EREBP family, bHLH family, and GRAS family

For the DEGs in the nine abovementioned TF families, the screening criteria were as follows: 1) the expression levels of the RPs at the three stages were all higher than those of the SPs regardless of herbicide application; 2) the gene expression of the RPs was higher than that of the SPs after spraying; 3) the DEGs were expressed only at the flowering stage and not at the seedling or tillering stage; and 4) the expression in the sprayed resistant population was higher than that in the sensitive population only at the flowering stage (Additional file 3: Tables S9-S17, sheet 1). Afterward, a total of 30 candidate genes were selected for screening resistance related flowering-regulated genes. qRT-PCR was then carried out to verify the expression of these 30 genes in four samples: the UFS, untreated resistant population at the flowering stage(UFR), TFS, and treated resistant population at the flowering stage(TFR). Twelve DEGs were related to the regulation of plant development, flowering, and stress response (Fig. 6a-b, Table 2); the remaining 18 genes were false positives (data not shown). The ABI3VP1, BHLH, and GRAS families each had one gene (CL18402.Contig2, CL6193.Contig3, and CL20691.Contig17), and the other nine genes belonged to the MADSbox family. The expression of four genes (CL4600.contig2, CL278.contig6, CL10951.contig2, and CL18402.contig2) was higher in the RPs than in the SPs under herbicide and water treatments. The expression of eight genes (CL15323.contig1, CL6626.contig8, CL10710.contig2, CL19935.contig11, CL7805.contig1, CL19935.contig11, CL6193.contig3, and CL20691.contig17) was slightly lower in the RPs than in the SPs under water treatment, but their expression levels increased rapidly after herbicide application and, consequently, were significantly higher than those of the SPs. Interspecific comparisons showed that the expression of 12 unigenes in the RPs was higher than that in the SPs under herbicide selective pressure, suggesting that these genes in the RPs likely promote reproductive growth (flowering and fruiting) under stress conditions: this phenomenon constitutes an unknown resistance mechanism (Fig. 6a).

Fig. 6
figure 6

a Comparison of the expression levels of 12 unigenes in differently treated samples and groups at the flowering stage. Three replicates were performed for each of the three biological replicates. The error bars represent the standard error of the mean. U: untreated; T: treated; F: flowering stage; S: sensitive population; R: resistant population. *significant difference p>0.05, **significant difference 0.01<P<0.05, ***significant difference P<0.01. b Heat map analysis of the 12 unigenes

Table 2 Differentially expressed genes (DEGs) as candidate transcription factors (TFs) between the susceptible plants (SPs) and resistant plants (RPs) of Polypogon fugax at the early flowering stagea

Under all treatment conditions, expression levels of the 12 unigenes were significantly higher in the UFS than in the TFS, which means that the expression of these genes in SPs was inhibited by herbicide selection. In the RPs, the expression of four unigenes (CL6193.contig3, CL20691.contig17, CL18402.contig2, and CL19935. contig9) in the UFR did not differ from that in the TFR (Fig. 6a). The expression of four unigenes (CL4600.contig2, CL19935.contig11, CL15323.contig1, and CL6626.contig8) was slightly lower in the UFR than in the TFR (Fig. 6a), and the expression of four unigenes (CL10710.contig2, CL7805.contig1, CL278.contig6, and CL10951.contig2) was slightly higher in the UFR than in the TFR (Fig. 6a). The results of the intraspecific comparison showed that herbicide selection pressure did not significantly influence the expression of these genes in the RPs (Fig. 6a). The results indicate that these genes have genetically adapted in RPs, while herbicides did not affect the growth and development of those plants.

Discussion

In this study, a P. fugax population resistant to clodinafop-propargyl and a susceptible population were selected. The RPs flowered 10-15 days earlier than did the SPs, but the mechanisms remain unclear. The goal of this study was to establish a resource to study transcriptomic patterns in P. fugax resistant or sensitive to clodinafop-propargyl in the absence or presence of herbicide and at different stages (seeding, tillering, and flowering), using experimental conditions as similar as possible to field conditions. The results revealed DEGs related to clodinafop-propargyl resistance in P. fugax. The assembled, annotated transcriptomes provide a genomic resource for understanding the molecular basis of P. fugax herbicide resistance.

Because there is no genomic resource for P. fugax, the Illumina technology was selected for sequencing, as it is the technology of choice for de novo transcriptome deep sequencing and assembly when a reference genome is absent [22]. Three replicates of 12 samples were sequenced by an Illumina HiSeq 4000 in this study, generating approximately 160.52 Gb bases in total. After discarding improper sequences, 206,041 unigenes were obtained; the total length, average length, N50, and GC content of these unigenes were 275,659,060 bp, 1337 bp, 1851 bp, and 51.51%, respectively. The N50 size of the contigs in this study was 1851 bp, which is higher than that recently obtained for plant de novo transcriptome assemblies based on Illumina sequence reads [23, 24]. The average contig size was 1337 bp, which matched the average length of gene coding sequences in grasses (1000 to 1300 bp) [25]. The use of functional annotation results revealed 152,966 coding DNA sequences (CDS), ESTScan (v3.0.2) was then used to predict the remaining unigenes, after which 11,716 additional CDS were obtained. Given that only 48 nucleotide sequences from P. fugax had been deposited in GenBank on May 22nd, 2017, our work tremendously increased the sequence data available for P. fugax.

The predicted peptide content of the P. fugax transcriptome was compared to that of five fully sequenced grass genomes. The five grass species belong to three major subfamilies of the Poaceae: Pooideae (Brachypodium distachyon and H. vulgare), Panicoideae (Zea mays and Sorghum bicolor), and Ehrhartoideae (Oryza sativa) [26]. The five grass genomes shared 58.64% of the identified protein families. These proportions were in agreement with a previous genome-wide study showing that genome peptide contents are largely shared among grass species, including peptide family representations [27]. In addition, good representation of the transcriptome of the aerial parts of P. fugax at different stages and confirmation of the relevance of the RNA-Seq-based expression data using qRT-PCR make the sequences a reliable resource for investigating the transcriptomic response to herbicide stress promoting early flowering in P. fugax.

Stress-induced flowering is a response to stress and is the ultimate stress adaptation, as plants can survive as a species if they flower and produce seeds under severe stress even when they cannot survive as individuals [17]. To validate the function of the candidate genes with respect to flowering in the resistant population, gene expression at the seedling and tillering stages served as a background. The present study used herbicide selection pressure when the RPs and SPs were about to flower. qRT-PCR was subsequently applied to analyze the candidate gene expression at the flowering stage to understand the relationship and mechanism of flowering and response to herbicide application in P. fugax.

The replication and diversity of MADS-box genes may be the factors affecting the morphological diversity of land plants and some angiosperms [28]. This gene family encodes conserved TFs and plays an important role in both vegetative and reproductive development. Among these genes, APETALA1 (AP1) is one of the earliest and most intensively studied genes. In Arabidopsis, AP1 deletion mutations delay flowering time and show a high frequency of floral meristem to inflorescence meristem transformation after flowering. On the other hand, constitutive expression of the AP1 gene changes the floral bud meristem to a floral meristem, and ectopic expression of the AP1 gene can significantly promote early flowering [29, 30]. The AP1 gene in herbaceous [31] and woody [32] plants also plays an important role in the initiation and development of flowers, and constitutive expression of AP1-like genes also contributes to early flowering. In contrast to the AP1 gene, overexpression of some MYB TFs (such as EPR1 and AtMYB44) in Arabidopsis leads to delayed flowering time [33, 34].

Table 2 summarizes the differentially expressed TFs identified in the present study. CL4600.Contig2 was annotated to the JOINTLESS-like protein (ADK55060.1). JOINTLESS (J) is a MADS-box gene that belongs in the same clade as the Arabidopsis flowering time genes SHORT VEGETATIVE PHASE (SVP) and AGAMOUS LIKE 24 (AGL24) [35]. Loss of J function causes premature termination of flower formation during inflorescence and reversion to a vegetative sympodial growth [36]. In addition, the formation of an inflorescence in tomato requires the interaction of J and a target of SINGLE FLOWER TRUSS (SFT) in the meristem [37].

CL10951.Contig2 was annotated to SUPPRESSOR OF OVEREXPRESSION OF CONSTANS1 (SOC1), which plays an important role in the regulation of flowering by integrating multiple flowering signals in Arabidopsis thaliana [38]. In the photoperiod pathway, SOC1 is regulated by CONSTANS (CO) through FLOWERING LOCUS T (FT), causing early flowering [39]. In addition, similar to SOC1, SOC1-like genes promote flowering when overexpressed, but some are also involved in floral development. Ectopic expression of UNSHAVEN (UNS), a SOC1-like gene of Petunia hybrids, leads to ectopic trichome formation on floral organs and the formation of petals into organs that exhibit leaf-like features and an early flowering phenotype [40, 41]. CL10951.Contig2 may be related to the promotion of flowering in P. fugax under herbicide stress [42].

AGAMOUS-like proteins (CL10710.Contig2, CL19935.Contig11, CL19935. Contig9, CL278.Contig6, and CL7805.Contig1) belong to the AG subfamily whose members are involved in the specification of floral reproductive organs and are also required for the normal development of carpels and fruits in Arabidopsis and Gossypium hirsutum [43] as well as for both drought stress responses [44] and regulation of post-germination growth [45]. AGAMOUS-like proteins under herbicide selection pressure were identified in this study, suggesting that these genes may promote early flowering and increased seed yield in resistant P. fugax plants.

The remaining three unigenes (CL18402.Contig2, CL6193.Contig3, and CL20691.Contig17) belong to the ABI3VP1, BHLH, and GRAS families, respectively. The genes of these three families play important roles in plant growth and development and stress responses [46,47,48], but their specific roles in P. fugax under herbicide stress need to be studied further.

The present study is not without limitations. Only two cultivars of P. fugax were studied, it is possible that cultivars from different regions could yield different results. In addition, transcriptome analysis is limited by available comparative data. It must be stressed that the present study does not provide any mechanistic data, but rather makes available transcriptomic data that could be used for the determination of those mechanisms. Additional studies are necessary to improve our knowledge and understanding of transcriptomics in plants.

Conclusions

In conclusion, the present study compared ACCase-resistant to ACCase-sensitive P. fugax, as RPs flower earlier and yield more seeds. The results revealed nine related resistance genes in the MADS-box family (which regulate flowering), and three genes involved in the regulation of growth and development. These data lay the foundation for the further exploration of the specific functions of these genes and for the study of transcriptomics in grasses.

Abbreviations

ACCase:

acetyl-CoA carboxylase

COG:

Cluster of Orthologous Groups of proteins

DEGs:

differently expressed genes

FDR:

false discovery rate

GO:

Gene Ontology

KEGG:

Kyoto Encyclopedia of Genes and Genomes

NR:

non-redundant

NT:

nucleotide

ORF:

open reading frame

RPKM:

reads per kilobase per million reads

SRA:

Short Read Archive

TFs:

transcription factors

References

  1. Tang W, Zhou F, Chen J, Zhou X. Resistance to ACCase-inhibiting herbicides in an Asia minor bluegrass (Polypogon Fugax) population in China. Pestic Biochem Physiol. 2014;108:16–20.

    Article  CAS  PubMed  Google Scholar 

  2. Xu HL, Zhu XD, Wang HC, Li J, Dong LY. Mechanism of resistance to fenoxaprop in Japanese foxtail (Alopecurus Japonicus) from China. Pestic Biochem Physiol. 2013;107:25–31.

    Article  CAS  PubMed  Google Scholar 

  3. Yu Q, Abdallah I, Han H, Owen M, Powles S. Distinct non-target site mechanisms endow resistance to glyphosate, ACCase and ALS-inhibiting herbicides in multiple herbicide-resistant Lolium Rigidum. Planta. 2009;230:713–23.

    Article  CAS  PubMed  Google Scholar 

  4. Shaner DL, Lindenmeyer RB, Ostlie MH. What have the mechanisms of resistance to glyphosate taught us? Pest Manag Sci. 2012;68:3–9.

    Article  CAS  PubMed  Google Scholar 

  5. Neve P. Challenges for herbicide resistance evolution and management: 50 years after Harper. Weed Res. 2007;47:365–9.

    Article  Google Scholar 

  6. Preston C, Powles SB. Evolution of herbicide resistance in weeds: initial frequency of target site-based resistance to acetolactate synthase-inhibiting herbicides in Lolium Rigidum. Heredity (Edinb). 2002;88:8–13.

    Article  CAS  Google Scholar 

  7. Vila-Aiub MM, Neve P, Steadmank J, Powles SB. Ecological fitness of a multiple-resistant Lolium Rigidum population: dynamics of seed germination and seedling emergence of resistant and susceptible phenotypes. J Appl Ecol. 2005;42:288–98.

    Article  Google Scholar 

  8. Wang T, Darmency H. Comparison of growth and yield of foxtail millet (Setaria italica) resistant and susceptible to acetyl-coenzyme A carboxylase inhibiting herbicides. In: 10th Colloque International Biologiedes Mauvaises Herbes de Dijon. Paris: Association Francaise de Protection des Plantes; 1996. p. 203-10.

  9. Wang T, Picard JC, Tian X, Darmency H. A herbicide-resistant ACCase 1781 Setaria mutant shows higher fitness than wild type. Heredity (Edinb). 2010;105:394–400.

    Article  CAS  Google Scholar 

  10. Bernier G, Perilleux C. A physiological overview of the genetics of flowering time control. Plant Biotechnol J. 2005;3:3–16.

    Article  CAS  PubMed  Google Scholar 

  11. Turnbull C. Long-distance regulation of flowering time. J Exp Bot. 2011;62:4399–413.

    Article  CAS  PubMed  Google Scholar 

  12. Song YH, Shim JS, Kinmonth-Schultz HA, Imaizumi T. Photoperiodic flowering: time measurement mechanisms in leaves. Annu Rev Plant Biol. 2015;66:441–64.

    Article  CAS  PubMed  Google Scholar 

  13. Yaish MW, Colasanti J, Rothstein SJ. The role of epigenetic processes in controlling flowering time in plants exposed to stress. J Exp Bot. 2011;62:3727–35.

    Article  CAS  PubMed  Google Scholar 

  14. Pieterse AH. Is flowering in Lemnaceae stress-induced? A review Aquatic Botany. 2012;104:1–4.

    Google Scholar 

  15. Riboni M, Robustelli Test A, Galbiati M, Tonelli C, Conti L. Environmental stress and flowering time: the photoperiodic connection. Plant Signal Behav. 2014;9:e29036.

    Article  PubMed Central  Google Scholar 

  16. Kazan K, Lyons R. The link between flowering time and stress tolerance. J Exp Bot. 2016;67:47–60.

    Article  CAS  PubMed  Google Scholar 

  17. Takeno K. Stress-induced flowering: the third category of flowering response. J Exp Bot. 2016;67:4925–34.

    Article  CAS  PubMed  Google Scholar 

  18. Garber M, Grabherr MG, Guttman M, Trapnell C. Computational methods for transcriptome annotation and quantification using RNA-seq. Nat Methods. 2011;8:469–77.

    Article  CAS  PubMed  Google Scholar 

  19. Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008;5:621–8.

    Article  CAS  PubMed  Google Scholar 

  20. Mistry J, Finn RD, Eddy SR, Bateman A, Punta M. Challenges in homology search: HMMER3 and convergent evolution of coiled-coil regions. Nucleic Acids Res. 2013;41:e121.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Ward JA, Ponnala L, Weber CA. Strategies for transcriptome analysis in nonmodel plants. Am J Bot. 2012;99:267–76.

    Article  CAS  PubMed  Google Scholar 

  23. Gardin JA, Gouzy J, Carrere S, Delye C. ALOMYbase, a resource to investigate non-target-site-based resistance to herbicides inhibiting acetolactate-synthase (ALS) in the major grass weed Alopecurus Myosuroides (black-grass). BMC Genomics. 2015;16:590.

    Article  PubMed  PubMed Central  Google Scholar 

  24. Chen J, Huang H, Wei S, Huang Z, Wang X, Zhang C. Investigating the mechanisms of glyphosate resistance in goosegrass (Eleusine Indica (L.) Gaertn.) by RNA sequencing technology. Plant J. 2017;89:407–15.

    Article  CAS  PubMed  Google Scholar 

  25. Carels N, Bernardi G. Two classes of genes in plants. Genetics. 2000;154:1819–25.

    CAS  PubMed  PubMed Central  Google Scholar 

  26. Bolot S, Abrouk M, Masood-Quraishi U, Stein N, Messing J, Feuillet C, et al. The 'inner circle' of the cereal genomes. Curr Opin Plant Biol. 2009;12:119–25.

    Article  CAS  PubMed  Google Scholar 

  27. Davidson RM, Gowda M, Moghe G, Lin H, Vaillancourt B, Shiu SH, et al. Comparative transcriptomics of three Poaceae species reveals patterns of gene expression evolution. Plant J. 2012;71:492–502.

    CAS  PubMed  Google Scholar 

  28. Kaufmann K, Melzer R, Theissen G. MIKC-type MADS-domain proteins: structural modularity, protein interactions and network evolution in land plants. Gene. 2005;347:183–98.

    Article  CAS  PubMed  Google Scholar 

  29. Mandel MA, Yanofsky MF. A gene triggering flower formation in Arabidopsis. Nature. 1995;377:522–4.

    Article  CAS  PubMed  Google Scholar 

  30. Weigel D, Nilsson O. A developmental switch sufficient for flower initiation in diverse plants. Nature. 1995;377:495–500.

    Article  CAS  PubMed  Google Scholar 

  31. Chi YJ, Huang F, Liu HC, Yang SP, Yu DY. An APETALA1-like gene of soybean regulates flowering time and specifies floral organs. J Plant Physiol. 2011;168:2251–9.

    Article  CAS  PubMed  Google Scholar 

  32. Wang J, Zhang X, Yan G, Zhou Y, Zhang K. Over-expression of the PaAP1 gene from sweet cherry (Prunus Avium L.) causes early flowering in Arabidopsis Thaliana. J Plant Physiol. 2013;170:315–20.

    Article  CAS  PubMed  Google Scholar 

  33. Kuno N, Moller SG, Shinomura T, Xu X, Chua NH, Furuya M. The novel MYB protein EARLY-PHYTOCHROME-RESPONSIVE1 is a component of a slave circadian oscillator in Arabidopsis. Plant Cell. 2003;15:2476–88.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Jung CK, Seo JS, Han SW, Koo YJ, Kim CH, Song SI, Nahm BH, Choi YD, Cheong JJ. Overexpression of ATMYB44 enhances stomatal closure to confer abiotic stress tolerance in transgenic arabidopsis. Plant Physiol 2008;146:623-35.

  35. Mao L, Begum D, Chuang HW, Budiman MA, Szymkowiak EJ, Irish EE, et al. JOINTLESS is a MADS-box gene controlling tomato flower abscission zone development. Nature. 2000;406:910–3.

    Article  CAS  PubMed  Google Scholar 

  36. Quinet M, Dubois C, Goffin MC, Chao J, Dielen V, Batoko H, et al. Characterization of tomato (Solanum Lycopersicum L.) mutants affected in their flowering time and in the morphogenesis of their reproductive structure. J Exp Bot. 2006;57:1381–90.

    Article  CAS  PubMed  Google Scholar 

  37. Thouet J, Quinet M, Lutts S, Kinet JM, Perilleux C. Repression of floral meristem fate is crucial in shaping tomato inflorescence. PLoS One. 2012;7:e31096.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Moon J, Suh SS, Lee H, Choi KR, Hong CB, Paek NC, et al. The SOC1 MADS-box gene integrates vernalization and gibberellin signals for flowering in Arabidopsis. Plant J. 2003;35:613–23.

    Article  CAS  PubMed  Google Scholar 

  39. Yoo SK, Chung KS, Kim J, Lee JH, Hong SM, Yoo SJ, et al. Constans activates SUPPRESSOR OF OVEREXPRESSION OF CONSTANS 1 through FLOWERING LOCUS T to promote flowering in Arabidopsis. Plant Physiol. 2005;139:770–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Ferrario S, Busscher J, Franken J, Gerats T, Vandenbussche M, Angenent GC, et al. Ectopic expression of the petunia MADS box gene UNSHAVEN accelerates flowering and confers leaf-like characteristics to floral organs in a dominant-negative manner. Plant Cell. 2004;16:1490–505.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Zhong X, Dai X, Xv J, Wu H, Liu B, Li H. Cloning and expression analysis of GmGAL1, SOC1 homolog gene in soybean. Mol Biol Rep. 2012;39:6967–74.

    Article  CAS  PubMed  Google Scholar 

  42. Liu S, Ma T, Ma L, Lin X. Ectopic expression of PVSOC1, a homolog of SOC1 form Phyllostachys Violascens, promotes flowering in Arabidopsis and rice. Acta Physiol Plant. 2016;38:166–74.

    Article  Google Scholar 

  43. Favaro R, Pinyopich A, Battaglia R, Kooiker M, Borghi L, Ditta G, et al. MADS-box protein complexes control carpel and ovule development in Arabidopsis. Plant Cell. 2003;15:2603–11.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Bechtold U, Penfold CA, Jenkins DJ, Legaie R, Moore JD, Lawson T, et al. Time-series Transcriptomics reveals that AGAMOUS-LIKE22 affects primary metabolism and developmental processes in drought-stressed Arabidopsis. Plant Cell. 2016;28:345–66.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Yu LH, Wu J, Zhang ZS, Miao ZQ, Zhao PX, Wang Z, et al. Arabidopsis MADS-box transcription factor AGL21 acts as environmental surveillance of seed germination by regulating ABI5 expression. Mol Plant. 2017;10:834–45.

    Article  CAS  PubMed  Google Scholar 

  46. Chen J, Tan RK, Guo XJ, Fu ZL, Wang Z, Zhang ZY, et al. Transcriptome analysis comparison of lipid biosynthesis in the leaves and developing seeds of Brassica Napus. PLoS One. 2015;10:e0126250.

    Article  PubMed  PubMed Central  Google Scholar 

  47. Ranjan R, Khurana R, Malik N, Badoni S, Parida SK, Kapoor S, et al. bHLH142 regulates various metabolic pathway-related genes to affect pollen development and anther dehiscence in rice. Sci Rep. 2017;7:43397.

    Article  PubMed  PubMed Central  Google Scholar 

  48. Song L, Tao L, Cui HP, Ling L, Guo CH. Genome-wide identification and expression analysis of the GRAS family proteins in Medicago truncatula. Acta Physiologiae Plantarum. 2017;39(4):93.

Download references

Acknowledgements

Not applicable.

Funding

This research was supported by the National Natural Science Foundation of China Project (31501658), the National Key Research and Development Program of China (2016YFD0201305), and the Special Scientific Research Fund of Agricultural Public Welfare Profession of China (201403030).

Availability of data and materials

All data generated or analyzed during this study are included in this article and its supplementary information files. The raw reads have been deposited in the NCBI Sequence Read Archive (SRA) database (BioProject PRJNA385696, SRP106591). The data are available at https://trace.ncbi.nlm.nih.gov/Traces/sra_sub/sub.cgi?acc=SRP106591&focus=SRP106591&from=list&action=show:STUDY.

Author information

Authors and Affiliations

Authors

Contributions

FYZ and YZ conceived and coordinated the study: designed, performed and analyzed the experiments and wrote the paper. WT, MW and TCG carried out the data collection, analyzed the data, and revised the paper. All authors reviewed the results and approved the final version of the manuscript.

Corresponding author

Correspondence to Fengyan Zhou.

Ethics declarations

Ethics approval and consent to participate

The plants were collected without permissions being sought for the nature of scientific research according to the law of the People’s Republic of China.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1: Table S1.

Primers used for qRT-PCR. (DOC 63 kb)

Additional file 2: Figure S1.

Morphological characteristics in resistant (R) and susceptible (S) Asian minor bluegrass plants at different growth stages. (TIFF 1872 kb)

Additional file 3: Table S2.

Summary of sequencing reads after filtering. Table S3. Quality metrics of unigenes. Table S4. Summary of functional annotation result. Table S5. All-unigenes with a COG classification. Table S6. All-unigenes categorized by GO terns. Table S7. Unigenes mapped onto KEGG pathways. Table S8. TF prediction results. Table S9. MYB family. Table S10. MYB-related family. Table S11. MADS-box family. Table S12. NAC family. Table S13. mTERF family. Table S14. ABI3VP1 family. Table S15. AP2-EREBP family. Table S16. bHLH family. Table S17. GRAS family. (ZIP 3364 kb)

Additional file 4: Figure S2.

Length distribution of unigenes. (TIFF 5378 kb)

Additional file 5: Figure S3.

Distribution of annotated species. (TIFF 3199 kb)

Additional file 6: Figure S4.

Functional distribution of COG annotation. (TIFF 6509 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhou, F., Zhang, Y., Tang, W. et al. Transcriptomics analysis of the flowering regulatory genes involved in the herbicide resistance of Asia minor bluegrass (Polypogon fugax). BMC Genomics 18, 953 (2017). https://0-doi-org.brum.beds.ac.uk/10.1186/s12864-017-4324-z

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://0-doi-org.brum.beds.ac.uk/10.1186/s12864-017-4324-z

Keywords