Skip to main content

Transcriptome analysis reveals major transcriptional changes during regrowth after mowing of red clover (Trifolium pratense)

Abstract

Background

Red clover (Trifolium pratense) is globally used as a fodder plant due its high nutritional value and soil improving qualities. In response to mowing, red clover exhibits specific morphological traits to compensate the loss of biomass. The morphological reaction is well described, but the underlying molecular mechanisms and its role for plants grown in the field are unclear.

Results

Here, we characterize the global transcriptional response to mowing of red clover by comparing plants grown under greenhouse conditions with plants growing on agriculturally used fields. Unexpectedly, we found that biotic and abiotic stress related changes of plants grown in the field overlay their regrowth related transcriptional changes and characterized transcription related protein families involved in these processes. Further, we can show that gibberellins, among other phytohormones, also contribute to the developmental processes related to regrowth after biomass-loss.

Conclusions

Our findings show that massive biomass loss triggers less transcriptional changes in field grown plants than their struggle with biotic and abiotic stresses and that gibberellins also play a role in the developmental program related to regrowth after mowing in red clover. Our results provide first insights into the physiological and developmental processes of mowing on red clover and may serve as a base for red clover yield improvement.

Background

Trifolium pratense (red clover) is an important worldwide forage crop and thus of great economic interest. This perennial plant offers several advantages like a high protein content and soil improving characteristics, which reduce the use of artificial nitrogen fertilizer and enhance livestock intake. Well-known disadvantages of red clover include poor persistence under several land use scenarios, like grazing or cutting [1,2,3]. T. pratense is a member of the Fabaceae (or legumes), which are, due to their economic value, among the most examined families in the plant kingdom with genome sequences available for species like Medicago truncatula (barrel clover) [4], Lotus japonicus (birdsfoot trefoil) [5], Glycine max (soy) [6], Phaseolus vulgaris (common bean) [7], Cicer arietinum (chickpea) [8], Vigna unguiculata (cowpea) [9], Trifolium subterraneum (subterranean clover) [10] T. medium (zigzag clover) [11], and T. pratense (red clover) [12, 13].

Facing today’s challenges such as an increased demand on food production in an era of global climate change, together with the aim to solve these problems in an environmental friendly and sustainable way, requires improvement of forage crops like T. pratense [14, 15]. T. pratense breeding aims to offer genotypes with improved key agronomic traits (dry matter yield, high quality, resistance to diseases and abiotic/biotic stress, persistency), while improving its regrowth ability [2, 16]. Unfortunately, the morphological investigations of several T. pratense populations showed a correlation of persistency with non-favorable traits, like small plant size and prostrate growth habit [17]. Moreover, most T. pratense cultivars or accessions are locally adapted and require their specific local conditions to show the favored traits [18, 19], which decreases the stability for individual traits in breeding efforts [20]. T. pratense exhibits significant intraspecific variation due to high intrapopulation genetic diversity, thus, persistence and performance in response to mowing or cutting depends on the variety as well as on the developmental stage at the moment of damage [21,22,23,24].

Persistency can be defined as a sustained forage yield over several growing periods [25] and is a complex trait influenced by a variety of abiotic and biotic factors, and the regrowth ability of a plant [26]. Plants with high regrowth ability can survive more frequent and intense biomass loss and are therefore more persistent. Decapitation or biomass loss due to herbivory or mowing triggers a complex reaction affected by environmental conditions, plant morphology, architecture, developmental stage and genotype [21]. After decapitation, the first known stress response in other legumes like M. sativa and P. sativum involves the production of phytohormones: cytokinins, auxins, and strigolactones [27,28,29]. In addition, the mobilization of energy reserves is activated [30]. Phenotypic plasticity of plant architecture in combination with alterations of hormone concentrations can be observed in P. sativum and T. pratense after decapitation [24, 29, 31]. However, the molecular processes allowing plants to thrive even after an enormous loss of biomass remain still unclear, even in Arabidopsis thaliana [32, 33].

Here, we compare the transcriptomes of mown (cut) vs. unmown (uncut) T. pratense plants from two different, well investigated field locations on the Biodiversity Exploratory “Heinich-Dün” [34] and greenhouse grown plants. Our field samples were subjected to standard agricultural treatment and we can thus discriminate transcriptional changes caused by abiotic factors and biotic interactions in the field from those that regulate regrowth. We present the identification and in silico characterization of putative developmental regulators differentially expressed in the regrowth phase after mowing in the field and in the greenhouse that may contribute to the regrowth response of T. pratense and demonstrate that gibberellin is a major regulator of specific aspects of the regrowth morphology in red clover.

Results

RNA-Seq results, de novo assembly, and functional description of contigs

The RNA-Seq produced a total number of short reads between 44.7 and 58.1 million for each library with two exceptions (Table S4) totaling 608,041,012 raw reads. The de novo assembly of the reference transcriptome of T. pratense produced 44,643 contigs, of which 41,505 contigs were annotated and 29,781 contigs were identified as plant specific. The minimal length of the contigs was 124 bp, the maximal length 15,551 bp (Table S5). After the de novo assembly of the T. pratense transcriptome, each library was mapped back against the reference transcriptome to determine the overall alignment rate, which was between 77.85 and 90.32% (Table S6).

63% of the 44,643 contigs could be mapped to a known locus of the T. pratense genome annotation [12, 35], 32% could be mapped to an unknown locus of the T. pratense genome and 5% could not be mapped to the T. pratense genome (Fig. S2). All plant-specific contigs were annotated with several databases (Table S3). To further verify the quality of our replicates, we identified the transcripts shared by the two replicates. We calculated TPM values for each transcript and discarded transcripts with TPM values < 1. The percentage of transcripts shared between the two replicates was between 90 and 94% for all treatments/localities, suggesting that the RNA-Seq data are highly reproducible (Table S8). In addition, we validated the RNA-Seq data via q-RT-PCR of four randomly chosen contigs tdn_146439 (LTP), k65_9861 (P5CS), tdn_69411 (PME44), and tdn_85889 (ENGase85A). The expression pattern of all four genes was congruent in qRT-PCR and RNA-Seq samples (Fig. S3). When single tissues (axillary meristem (AM) and leaves (L) were tested individually, only LTP showed expression similar to the RNA-Seq data.

Differentially expressed gene analysis reveals diverse subsets of genes involved in regrowth influenced by location and environmental conditions

To identify gene expression responses underlying the regrowth response after mowing, a digital gene expression analysis was performed comparing field A mown vs. field A non-mown (FaM vs. FaNM); field B mown vs. field B non-mown (FbM vs. FbNM); greenhouse mown vs. greenhouse non-mown (GM vs. GNM) to identify DEGs (Table S9) from mown plants. Interestingly, using the |log2foldchange| ≥ 2, the number of differentially expressed genes (DEGs) is rather similar in all comparisons, ranging from 119 (GM vs. GNM) to 142 (FaM vs. FaNM) (Table 1).

Table 1 The numbers of differentially expressed transcripts (contigs) between libraries with changes equal or above |log2foldchange| 2. Upregulation for each comparison is shown

GO enrichment

We performed a GO enrichment analysis with the DEGs of each group to obtain a differential view on the transcriptional changes occurring in relation to regrowth (Tables 2, 3, 4). The identified GO terms for each sample (FaM, FbM, GM, FaNM, FbNM, and GNM) were then compared in a mown vs. non-mown manner to identify GO terms specific for the respective treatments. The results revealed that GO terms involved mainly in general metabolic processes and pathways, as well as general reactions are enriched in non-mown plants including i.e. the GO terms “protein metabolic process”, “metabolic process”, “cellular process”, “catabolic process”, “biosynthetic process” (Fig. S4 and Table S11). Within mown samples we found the following GO terms enriched: “nucleic acid binding “(GM); GO terms related to photosynthesis (“photosynthesis”, “thylakoid“), cell components and protein transport (“Golgi apparatus”, “cytoplasm“) and related to regrowth and stress response (“generation of precursor metabolites and energy”, “cell growth”, “cell communication“) (Fig. S4 and Table S11). Within the GO term “cell growth” the contigs GIBBERELLIN-REGULATED PROTEIN 1 which is involved in cell elongation and ROOT HAIR DEFECTIVE 3, a protein involved in root hair growth are present. For FbM we found the GO term related to metabolic processes (“metabolic process”, “lipid metabolic process”), cell related (“cytoplasm”, “extracellular space”), enzymatic and catabolic processes (“enzyme regulator activity”, “catalytic activity”) and the GO term “binding”, which included a contig encoding for “V”, a protein involved in the ethylene biosynthesis.

Table 2 Twenty most strongly differentially expressed genes of the GM vs. GNM analysis. Shown are the transcript name, |log2foldchange| ≥ 2 of the corresponding transcript, the library in which the transcript is upregulated (pattern), gene name based on T. pratense genome annotation, corresponding Phytozome description, gene name and species name of the next homologs and A. thaliana gene name, and locus name based on information available on Tair
Table 3 Twenty most strongly differentially expressed genes of the FaM vs. FaNM analysis. Shown are the transcript name, |log2foldchange| ≥ 2 of the corresponding transcript, the library in which the transcript is upregulated (pattern), gene name based on T. pratense genome annotation, corresponding Phytozome description, gene name and species name of the next homologs and A. thaliana gene name, and locus name based on information available on Tair
Table 4 Twenty most strongly differentially expressed genes of the FbM vs. FbNM analysis. Shown are the transcript name, |log2foldchange| ≥ 2 of the corresponding transcript, the library in which the transcript is upregulated (pattern), gene name based on T. pratense genome annotation, corresponding Phytozome description, gene name and species name of the next homologs and A. thaliana gene name, and locus name based on information available on Tair

Interestingly, most functional groups differ between the field and greenhouse location (Fig. 1a-c), for example, more genes related to growth are upregulated in the non-mown Fa location but in the Fb and greenhouse location, they more genes are upregulated in the mown plants. Only genes related to biotic stress processes were upregulated in all unmown locations and more transposon-related genes are upregulated in mown plants (Fig. 1 a-c).

Fig. 1
figure1

Transcriptome analysis showing biological processes relevant after mowing. a-c DEGs with a |log2foldchange| ≥ 2 are shown in percentage share of each class to the corresponding DEG list in bar charts, dark grey bars show datasets from mown, light grey bars indicate datasets from non-mown plants, showing the changing distribution of DEGs in the various biological processes after mowing. a Classes of DEGs from field A, b Classes of DEGs from field B, c Classes of DEGs from greenhouse grown plants. d-e Venn diagrams showing the number of shared upregulated genes within the mown samples (d) and the number of shared genes within the unmown samples (e) with a |log2foldchange| ≥ 2. Blue circles indicate genes upregulated in the greenhouse, green indicates genes upregulated in field A, and red indicates genes upregulated in field B. f Number of DEGs belonging to the class “phytohormones” within the DEG list of the two field and the greenhouse transcriptomes. The pie chart shows the number of DEGs related to the different phytohormones (abscisic acid, auxin, genes common between the auxin and cytokinin pathway, cytokinin, ethylene, gibberellins, jasmonic acid, and salicylic acid)

The photosynthesis- and phytohormone-related genes of field A show a similar pattern to the field B plants as do the phytohormone- and signaling related genes. Genes related to development, general cell functions and transcription are also similar between field A and the greenhouse grown plants, such that more transcription- and development-related genes are upregulated in mown plants. And unexpectedly, senescence-related genes are upregulated in mown plants of field B. However, as our analysis cannot discriminate between activating and repressing factors of senescence, we cannot conclude from our data whether the mown plants have activated or repressed their senescence program.

The largest group of differentially expressed genes is the one related to biotic stress with up to 38% differentially expressed genes in one location (field B, Fig. 1 b). This suggests that biotic stresses play a prominent role in non-mown plants. A similar phenomenon can be observed for growth related processes, where up to 24% genes were upregulated in the mown and unmown plants indicating that different growth programs are active in mown vs. unmown plants.

Taken together we can state that mown plants in all locations change their transcriptional programs upon mowing suggesting that they massively change their metabolism and signaling processes. However, the molecular answer to substantial biomass loss differed between all three locations.

To find similarly regulated genes between the treatments and/or locations, Venn diagrams were generated to compare the number of shared DEGs within the mown samples and the non- mown samples (Fig. 1 d-e, Table S12). Within the mown samples, we detected no overlap between the groups with the exception of four upregulated DEGs in the two field transcriptomes (FbM and FaM (Fig. 1 d). Among these genes were two that could not be annotated, on gene encoding tubulin beta chain 2, and on encoding the GA responsive protein GAST1. Within the non-mown samples, also four genes were shared between the field transcriptomes (FbNM and FaNM)). Among these genes encoding a Chitinase A (class III), expressed only under environmental stress conditions and involved in plant immunity, a Leucin-rich repeat (LRR) family protein, a protein arginine methyltransferase (ATPRMT6), and one gene that we could not annotate. The MADS-box transcription factor-encoding gene SEP1 was the only gene shared between the field B and the greenhouse (Fig. 1 e). No genes were shared between all three samples, neither in the mown treatment, nor in the non-mown treatment.

Genes involved in developmental processes

We were then interested to identify transcriptional changes in genes directing developmental processes required for the regrowth process. Thus, the results of the DEG analysis were restructured such that the DEGs were grouped in 16 descriptive classes defined by database and literature mining (Tables S3 and S10). Those classes (abiotic stress, abiotic/biotic stress, biotic stress, development, general cell function, growth, metabolism, photosynthesis, phytohormone, secondary metabolite biosynthesis, senescence, signaling, symbiosis, transcription, transposon, no available annotation) describe major functional groups and serve to broadly categorize the potential role of a gene (Table S10).

We compared the top 20 DEGs of mown vs. not mown plants and observed that the greenhouse plants displayed more DEGs (classes growth, transcription, and phytohormone) involved in regrowth processes (Fig. 1c, Table S3). Three DEGs involved in growth, two phytohormone genes, and two transcriptional regulators are among the top 20 DEGs, while ten DEGs are related to biotic and abiotic stress in the greenhouse (Table S3). The top 20 DEGs of field A grown plants include four growth related, three development related and five stress-related DEGs. The top 20 DEGs of field B grown plants included only two growth related and six stress-related transcripts. Taken together, the greenhouse grown plants showed most DEGs related to growth, transcription, and phytohormone actions indicative of a regrowth reaction, as they grew under less stressful conditions than the field grown plants, for which stress related DEGs were more dominant (Fig. 1 a-c /Tables 2, 3, 4).

Phytohormone-related genes

We were interested in the contribution of individual phytohormones to the regrowth reaction in T. pratense, as they are known to play a major role in the regulation of development and stress response. We identified DEGs related to phytohormone synthesis, homeostasis, transport, and signaling for all major classes of phytohormones in the datasets. The four phytohormones with the most associated DEGs were: abscisic acid (8 DEGs), gibberellin (8 DEGs), salicylic acid (6 DEGs), and auxin (5 DEGs) (Fig. 1 f). Abscisic acid and salicylic acid are well-known to be involved in response to drought and abiotic/biotic stress, respectively. Auxin is the major phytohormone required for growth and cell division regulation and thus, we expected DEGs related to these phytohormones to be abundant in our analysis. However, unexpectedly, eight DEGs with gibberellin association were found. As gibberellins regulate growth in response to stresses but have so far not been associated with regrowth after biomass loss, we suggest gibberellin as a novel candidate phytohormone to influence the regrowth response.

Specific transcriptional regulator families are differentially expressed during the regrowth process

As the regulation of stress response, growth and development depends on differential activity of transcription factors, we aimed to identify transcriptional relevant to the biological processes occurring 2 weeks after mowing by mapping the transcriptome to the PlnTFDB [36]. All members of a specific transcriptional regulator family (TRF) were identified in silico and their expression was compared between mown and unmown plants (Table S13). Figure 2 shows TRFs with significantly differential expression between mown and unmown conditions in at least 10% of their members, Fig. S5 and Table S13 includes also those TRFs with 5% of their members regulated differentially upon mowing.

Fig. 2
figure2

Expressed TRF members in mown and non-mown T. pratense plants. The y-axis shows the number of expressed genes (TPM value over 5) that are members of the respective TRF based on PlnTFDB, indicating significant differences in the expression of TRF members upon mowing. Names of the transcriptomes and TRFs are given on the x-axis. Expression of TRF members was compared in a pairwise manner (GM vs GNM, FaM vs FaNM, FbM vs FbNM). Shown are only those plant TRFs in which at least one of the comparisons resulted in a difference of more than 10% of the genes significantly upregulated in either the mown or the unmown condition (orange bars)

17 TRFs were identified of which at least 10% of the members showed differential expression in mown versus unmown comparisons (Fig. 2): ABI3VP1, AP2-EREBP, C2C2-Dof, C2C2-GATA, GRAS, HSF, LOB, MADS, mTERF, MYB, NAC, PHD, SBP, SNF2, TCP, TRAF, WRKY

On field A, the AP2-EREBP, LOB, MADS, MYB, NAC, PHD, SBP, TCP and WRKY TRFs were more prominent in unmown plants, and only the HSF TRFs were more prominent upon mowing. On field B, ABI3VP1, C2C2-Dof, GRAS, HSF, LOB, MADS, mTERF, SNF2, TRAF, and WRKY TRFs were reduced upon mowing. In the greenhouse-grown plants, members of ABI3VP1, C2C2-Dof, C2C2-GATA, and GRAS, show increased numbers in response to mowing. In addition, ARF, C2H2, homeobox, MYB, NAC, and TRAF TRFs show changes in expression in all locations, albeit with only between 5 and 10% of the members being differentially expressed (Table S13).

Two TRFs showed a repression of expression upon mowing: 10% of the WRKY transcripts were less abundant in mown plants regardless of the provenance. In addition, MADS-box transcripts were found upregulated as well, but only in the field-derived transcriptomes. Generally, only four of the 17 TRFs analyzed here showed significant changes in expression towards mowing in the greenhouse-derived plants, suggesting that they react less strongly towards mowing than the field-derived plants. Six TRFs (AP2-EREBP, MYB, NAC, PHD, SBP, and TCP) showed transcriptional changes in reaction to mowing only in field A while only three TRFs (mTERF, SNF2, TRAF) showed this only in field B, suggesting that the combination of biotic and abiotic factors with mowing differed between the two field locations, and, in a similar way, between the field locations and the greenhouse.

Notably, only the C2C2-GATA TRF showed transcriptional changes in at least 10% of its members towards mowing under greenhouse but not under field-conditions, indicating that transcriptional changes in reaction to other biotic and abiotic factors may overlay the regrowth reaction. Taken together, the TRF analysis showed that the reaction towards mowing induces transcriptional changes in only a subset of TRFs, suggesting that those play a major role in relieving the stress of biomass loss and regrowth.

Gibberellins are also important regulators after mowing in red clover

We have shown previously (Fig. 1g) that genes related to gibberellins are also differentially expressed, even though GA is not well-known to regulate biological processes related to loss of biomass. We then wanted to know if GA is relevant for the regulation of regrowth and treated red clover plants exogenously with GA3.

A weekly gibberellin application during the regrowth process led to significant and specific changes in morphology (Fig. 3). Previous work suggested that regrowing plants produce smaller and rounder leaflets with shorter petioles than non-mown plants [24]. Thus, number of leaves, shoots and inflorescences, leaf area, and the roundness of leaflets were measured in this experiment (Fig. 3, Fig. S6). The first visible effects of gibberellin treatment were recognized after 1.5 weeks, showing a significant higher leaflet area of gibberellin treated. Later it was observed that the petioles of treated plants were on average twice as long as petioles of untreated plants (16.7 ± 1.9 cm and 8 ± 1.2 cm, respectively). Leaflets of gibberellin treated plants had with 4.7 ± 0.9 cm2 almost double the size when compared with those of untreated plants (2.4 ± 0.6 cm2). However, gibberellin treated plants produced only 30% more total leaf area than control plants. Other morphological traits such as number of inflorescences, leaves, and shoots remained unaffected by the gibberellin treatment (Fig. S6). In summary, mown plants normally produce leaves with shorter petioles, restrict their leaflet area and their leaves become rounder. Gibberellin treatment partially alleviated these developmental changes such that the mown, gibberellin treated plants produced larger leaves with longer petioles while the leaf shape was unaffected by gibberellin treatment.

Fig. 3
figure3

Gibberellin treatment affects regrowth of mown T. pratense plants after biomass loss. a) The development of leaflet area in cm2, b) the development of petiole length in cm. Regrowing plants sprayed with GA3 showed significantly bigger leaflets and longer petioles. The graphs show average values of the respective plant growth parameters for each sampling date during 4 weeks of observation and the 95% confidence interval. Blue, GA3 treated plants; orange, mock-treated plants

Discussion

RNA-Seq and assembly

The de novo assembly in combination with a reference-based approach for the annotation led to the identification of 44,643 contigs of which 29,781 were annotated as plant-specific (Fig. S2). With the prior de novo assembly, 4051 additional contigs were obtained that have not been found in the T. pratense 1.0 (GCA_000583005.2) genome [12, 35]. The estimated genome size of T. pratense is ~ 440 Mbp [27]. The T. pratense transcriptome data in this study was ~ 55 Mbp in size, corresponding to ~ 12.5% transcribed regions in the T. pratense genome, which is within the range of previously published transcriptomes (~ 10% (42 Mbp) [37]). Interestingly, we found plant-specific, previously unreported contigs suggesting that the T. pratense genome might need improvement in terms of sequencing coverage and protein coding sequence annotation.

Biotic and abiotic stresses contribute to differential gene expression

Plant grown on fields face different stressors when compared to greenhouse grown plants and we were interested in how the field conditions contributed to differential gene expression. The transcriptome comparisons between locations revealed that the mown greenhouse plants showed the highest percentage of DEGs possibly involved in regrowth processes (Fig. 1a-c). In contrast, the field transcriptomes displayed patterns of abiotic and biotic stress reactions. Comparisons of the top 20 DEGs of the unmown field transcriptomes showed that plants grown on field A and B faced more biotic stress than abiotic stress. One of the upregulated genes in field A is a chitinase homolog suggesting that those plants are under attack of fungi and/or insects. Follow-up analyses to correlate environmental conditions, as well as biotic and abiotic stresses monitored within the Biodiversity Exploratories with differential gene expression at the two field locations would be an interesting project but is beyond the scope of this work. In contrast, the top 20 DE transcripts of the greenhouse plants include phytohormone- and transcription-related genes, but also a high proportion of biotic and abiotic stress-related genes. This suggests that also these plants have to cope with stresses, but to a lesser extent emphasizing their regrowth reaction more strongly within the top 20 DEGs. Generally, the non-mown plants show a much higher number of upregulated biotic stress-related genes during a phase in their life when senescence commences and they become more susceptible to pathogen attacks. The mown plants during their regrowth phase are not senescing and their younger organs seem to be less affected by biotic stress.

Cell walls are remodeled after mowing

After massive biomass loss, like mowing inflicts on T. pratense, plants firstly need to seal wounded tissues. Several transcriptional regulators known to play a role in the tissue-reunion processes were identified in Solanum lycopersicum, Cucumis sativus, and A. thaliana (reviewed in [38]). Homologs of these genes were also identified to be differentially regulated in the T. pratense transcriptome after mowing (Table S14), such as several members of the Auxin Response Factor (ARF) family or the No Apical Meristem (NAM) family member ANAC071, the homolog of the most highly upregulated transcription factor in greenhouse grown plants after mowing (Table S14) [39]. suggested that high levels of auxin induce the expression of ANAC071 via ARF6 and ARF8 (in the upper part of incised stems), at the same time, reduced auxin levels directly after the cutting activate the expression of RAP2.6 L. In addition auxin signaling via ARF6 and ARF8 influences jasmonic acid synthesis via the activation of DAD1, thus together with LOX2 increases RAP2.6 L expression during tissue reunion in A. thaliana (Table S14) [39]. Further it was demonstrated that ANAC071 can initiate the expression of members of the xyloglucan endotransglucosylase/hydrolases family (XTH20 and XTH19) which recombine hemicellulose chains to drive the cell expansion during tissue reunion [40]. Interestingly, we found that all members of the cell wall remodeling pathway show distinct expression patterns- Some are upregulated in mown plants including for example XTH32 (k69_7012, upregulated in FbM, tdn_94651, upregulated in GM, FaM and FbM), XTH6 (tdn_91763, upregulated in GM), XTH8 (k71_5058, upregulated in GM, FbM), XTH9 (tdn_113578, upregulated in GM), XTHA (tdn_87930, upregulated in GM), LOX2 (tdn_156279, upregulated FbM), and ARF8 (tdn_156886 upregulated in GM, tdn_156890 upregulated in GM) (Table S9, S13 and S15). This implies that the early steps in the regrowth reaction are conserved in core eudicots and that the cell wall remodeling processes continue at least 2 weeks after mowing.

Loss of auxin control on axillary buds can be detected two weeks after decapitation

Axillary buds of decapitated P. sativum plants export auxin upon growth activation, a process mediated by the P. sativum PIN1 homolog PsPIN1. Upon the loss of apical auxin transport, PsPIN1 polarization changes and this new polarization is causing auxin export from dormant axillary buds and is required for their activation [41]. Subcellular targeting and polarization of PsPIN1 starts about 6 h after decapitation and then PsPIN1 expression increases. However, it remains unclear if the elevated PsPIN1 expression is maintained for a prolonged period [42]. Our data show a higher expression of the three PIN1 homologs in greenhouse-grown mown plants when compared to the non-mown control plants (Table S2), indicating a sustained expression of the homolog of PIN1 even after 2 weeks of biomass loss, which might help to activate the remaining dormant buds in T. pratense.

Gibberellin-related genes influence regrowth of T. pratense in concert with other phytohormones

Wounding induces a first stress response, activating the interplay of the phytohormones jasmonic acid, salicylic acid, and ethylene. This allows an individual response to various abiotic and biotic stresses, and the differentiation between wounding inflicted by physical forces, pathogens, or herbivory [43,44,45,46,47,48]. Abscisic acid is required for the fine tuning of the jasmonic acid/salicylic acid/ethylene induced stress response by i.e. suppression of other phytohormones [49]. After the first stress response, additional phytohormones are involved in the regrowth of the plant. Auxin, cytokinin, strigolactone, and gibberellin become involved in a later stage. Following the initiation of shoot outgrowth induced by auxin and cytokinin, an increased gibberellin concentration allows for shoot elongation [31, 50,51,52]. In addition, auxin, cytokinin, and salicylic acid are involved in the shoot branching, where high levels of auxin and salicylic acid have a suppressing effect on lateral bud outgrowth. High levels of cytokinin promote shoot outgrowth which was shown in A. thaliana, O. sativa, and P. sativum [53,54,55,56]. As we were mainly interested in processes that happen approximately 2 weeks after cutting and as the role of those phytohormones was already studied, we concentrated on the role of gibberellin during regrowth, which was also found as one the phytohormones with the most associated genes in our transcriptomes. Gibberellins are involved in multiple aspects of plant development like cell elongation, flowering time regulation, and seed germination. Consequently, genes encoding for proteins involved in the synthesis, perception, and catabolism of the various gibberellins can be assumed to influence plant form. Our RNA-Seq data showed a high abundance of differentially expressed gibberellin associated genes (Fig. 1 F, Table S13)) which may be connected to the morphological changes after mowing, such as rounder leaves, temporary dwarf-like appearance, and higher cumulative biomass production in mown plants [24].

When analyzing the morphological effects of gibberellin application to mown plants (Fig. 3)), external gibberellin application led to the disappearance of specific traits typical for the mowing response. Mown plants developed shorter petioles and their leaf size area was smaller [24], but when treated with gibberellin, leaves and petioles grow up to the size seen in unmown plants.

The cell-expansion and proliferation promoting abilities of gibberellins via stimulation of the degradation of growth-repressing DELLA proteins are well established [57]. The length increase of petioles in gibberellin treated mown plants is in line with reported data from non-mown Pisum sativum plants, but in those, leaf sizes remained unchanged after gibberellin treatment [58], suggesting a more specific role for gibberellin in the regrowth reaction after biomass loss in red clover. Moreover, it was shown in A. thaliana that elevated gibberellin concentrations enhance cell-division rates in the distal end of leaves (reviewed in [59]). If these results are transferred to T. pratense, gibberellin treatment should result in longer leaflets after gibberellin treatment of mown plants. Interestingly, leaf shape did not change, but only the size increased suggesting a regrowth-specific shift of growth pattern which is unaffected by gibberellin but similar to leaf shape of juvenile plants [24]. However, one can also assume a cell-division pattern in red clover leaves that is distinct from the one reported for A. thaliana and may react more uniformly to enhance gibberellin concentrations.

Interestingly, gibberellin treatment of mown T. pratense plants does not generally lead to stronger longitudinal growth as leaves retained the round shape characteristic for untreated mown plants. These regrowth-specific characteristics can also be found in other species, for example in A. thaliana, Fragaria ananassa, Duchesnea indica, and G. max, gibberellin treatment causes elongated petioles and increased leaf sizes and a more erect growth habit [60,61,62,63]. This may suggest a new way to increase the accumulation of biomass, suitable for animal fodder. Previous experiments with the grasses Leymus chinensis and Lolium perenne showed gibberellin action to be limited by N fertilization [64, 65]. Red clover, living in symbiosis with nitrogen fixing bacteria, is not dependent on additional N fertilization and can produce high-protein content biomass without fertilizer on poor soils.

Methods

Plant growth conditions, gibberellin treatment, tissue sampling, RNA extraction, cDNA library construction, and RNA-Seq

Plant material for RNA-Seq was collected from three locations (fields and greenhouse, Fig. S1 and Table S1). Field plant tissue for RNA-Seq was sampled on 06.11.2014 within the area of the long-term open research platform Biodiversity Exploratory “Hainich-Dün” [34], located in Thuringia, Germany. Collection permits from farmers and local authorities were obtained centrally by the Biodiversity Exploratory research platform. T. pratense is an agriculturally used species native to Germany and does not fall under the Nagoya protocol. ITS sequencing with subsequent database comparisons identified the species collected in the field. The Material was sampled on four neighboring sites; two mown pastures (FaM and FbM) and two meadows that were non-mown (FaNM and FbNM) (Fig. S1 and Table S1). In the year of the sampling, the non-mown pastures and meadows were not grazed upon or mown. The two treatments were comparable as these were the closest experimental plots neighboring one mown and one un-mown plot. The experimental plots were managed as normal agricultural fields and were populated with comparable red clover, as these are wild, established populations. For the greenhouse samples, seeds of regional T. pratense populations (from a region covering mainly Thuringia, Saxony, Saxony-Anhalt, Thuringian Forest and Uckermarck, Germany) were obtained from the Rieger Hofmann seed company (Blaufelden, Germany). Plants in the greenhouse were grown in 23 °C with 16 h of light in pots of 12 cm diameter, watered daily, and compound fertilizer (8′8’6′+) was given every 10 days. After 122 days after sowing, half of the plants were cut to 5 cm (GM and GNM). Material from mown plants was sampled approximately 14 days after mowing/cutting, to avoid sequencing of the transcripts related to the first stress response [38]. After collection, the samples were snap frozen in liquid nitrogen.. For each group (FaM, FbM, GM, FaNM, FbNM, and GNM) shoot and leaf material of eight plants was collected. The respective eight plants were separated into two biological replicates of four plants whose RNA was later pooled. As reviewed in [38], the expected time for tissue reunion and wound closure accounts approximately 7 days (cucumber and tomato) to 14 days (A. thaliana). Based on this information, we assumed that the first stress response and the initiation of regrowth in T. pratense will be approximately 2 weeks after cutting/mowing.

RNA was extracted using NucleoSpin® RNA Plant Kit (Macherey-Nagel GmbH & Co. KG, Düren, Germany) according to the manufacturer’s instructions. For each replicate, equal amounts of RNA of four plants were pooled. Preparation of the cDNA libraries and the strand-specific sequencing were conducted by Eurofins Genomics (Ebersberg, Germany). The RNA-Seq libraries were sequenced on an Illumina Hiseq2000 platform with chemistry v3.0, creating 2 × 100 bp paired end reads.

Assembly of reference transcriptome and annotation

The raw-read-quality of the RNA-Seq data was analyzed with FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc). Illumina adapter and low quality regions were trimmed using Trimmomatic [42] with ILLUMINACLIP, SLIDINGWINDOW:5:20 and MINLEN:50 options. Quality trimmed reads were pooled and digitally normalized [66]. Multiple de novo assemblies were computed using Trinity [67] and Oases [68] with all odd k-mer parameters between 19 and 85. In addition, a genome guided assembly was performed with Trinity using the draft genome of T. pratense 1.0 (GCA_000583005.2) [12, 35]. The resulting contigs were screened for potential coding sequences (CDS) using TransDecoder (https://transdecoder.github.io/). The EvidentialGene pipeline (http://arthropods.eugenes.org/EvidentialGene/about/EvidentialGene _trassembly_pipe.html) was used to merge and filter the contigs based on the TransDecoder CDS prediction. Completeness of the final contig was confirmed by computing the mapping-rate of the non-normalized reads to the contigs. The raw sequence reads can be found at NCBI: PRJNA561285.

The contigs were uploaded to the “Sequence Analysis and Management System” (SAMS) [69] for functional annotation with the SwissProt [70], TrEMBL [71] and Phytozome [72] (e-value cutoff of 1e-5) databases. Additionally, attributes like gene names or functional descriptions were extracted from the blast hits. Contigs were mapped to the T. pratense reference genome using gmap [73]. The mapped contigs were compared to the reference annotation of the draft genome. Contigs with similar exon/intron chains (including contained contigs) to a reference gene were assigned the reference gene ID. All non-Viridiplantae contigs were discarded. Transcription factors were identified using a blastp search of the protein sequences against the plant transcription factor database Potsdam (PlnTFDB) ([36], version 3.0 protein database with an e-value cutoff of 1e-20. All functional annotations of transcripts can be found in the Table S2).

qRT-PCR confirmation of RNA-Seq

qRT-PCR samples included three biological replicates of leaves and axial shoot meristems from common garden experiment-grown plants (for growth conditions of the see Herbert et al. (2018)), as well as the 12 samples used for RNA-Seq. Sample treatment and RNA extraction were performed as described above. First strand cDNA was synthesized by using the RevertAid™ H-Minus First Strand cDNA Synthesis Kit (Thermo Scientific) following the manufacturer protocol.

PerlPrimer software (Marshall, 2004) was used to design qRT-PCR primer (Table S6), and primer efficiency tests were carried out with the same cycler settings, with a standard cDNA dilution series as template. The Roche LightCycler® 480 II system (Roche, Basel, Switzerland) was used for qRT-PCR. The total reaction volume of each sample was 20 μl consisting of 5 μl of 1:50 diluted cDNA template, 1 μl of each primer (10 μM), 3 μl sterile H2O, and 10 μl SYBR Green I Master (Roche). The qRT-PCR was carried out with the following cycler settings: 95 °C for 5 min followed by 45 cycles of 95 °C for 10 s, 60 °C for 10 s and 72 °C for 10 s. Contig k65_5754 was used as reference gene [37]. Three biological replicates with two technical replicates were used for each analyzed gene and primer combination. Only primers with an amplification efficiency between 1.8 and 2 were used. For the transcriptome library samples, only two biological replicates were used in the qRT-PCR analysis. Calculations were carried out as described in Pfaffl (2001).

Differential gene expression analysis, enrichment analysis, and classification of differentially expressed genes

Read counts for each contig of the final assembly in each sample were computed using RSEM [74] with bowtie mapping. To identify differentially expressed (DEGs) T. pratense genes, a pairwise comparison of all treatments was preformed using the DESeq2 [75] tool with FDR ≤ 0.01 and |log2foldchange| ≥ 2 between FaM and FaNM, FbM and FbNM; GM and GNM respectively. The top 20 DEGs were determined for each comparison based on the expression strength (log2foldchange). Homologs in the next closest species and A. thaliana for each T. pratense candidate gene were identified based on the T. pratense genome sequence deposited in Phytozome [72]. TPM (transcript per million) values were calculated to estimate contig expression level [76] (for RSEM read counts see Table S16).

We used the description and gene names obtained from TrEMBLE and SwissProt to search the UniProt [77], NCBI [78] and TAIR [79] databases to obtain further information (Table S3). Raw reads that were assembled to contigs, exhibiting a gene structure (ORF) and attained a putative annotation referred to below as genes.

GO enrichment analysis and Blast2Go analysis of T. pratense genomes

To further explore the digital gene expression results and to find more candidate genes/ to identify differentially expressed gene clusters, an enrichment analysis with Gene Ontology (GO) terms [80,81,82] was performed. For each pairwise comparison, the up-regulated genes were screened for enriched and depleted GO terms using the GOSeq package [83] separately for each treatment. Identified GO terms for each pairwise comparison were then also compared in a mown vs. non-mown manner to show treatment specific GO terms. The results of this analysis were visualized with the program GOplot [84] implemented in RStudio [85] with the program R [86].

Two local BLAST searches [87] with word-size of 3, e-value of 1.0e-3 and HSP length cutoff of 33 were performed against the PlnTFDB using Blast2GO [88]. Only the blast hits with the highest similarity were used for further comparisons (number of BLAST hits = 1), sequences with a similarity below 50% and an e-value higher than 1.0e-4 were omitted. The Blast2GO output was compared with an in-house python3-script utilizing NumPy (https://numpy.org/), Pandas (https://pandas.pydata.org/) and Seaborn (https://seaborn.pydata.org/) applying the list of transcription factors (TF) downloaded from PlnTFDB. We searched the Uniprot database hits for development and phytohormone related genes. Subsequently, we searched for gene IDs of gibberellin genes in our annotated T. pratense transcriptomes. Matches were filtered based on TPM values and classified based on biosynthesis and its regulation, catabolism, activation/repression or signaling/response, and the corresponding expression patterns within the transcriptome were identified.

Gibberellin treatment

To assess the effect of gibberellin during the regrowth reaction of T. pratense, 14 red clover plants were mown as described in [24]. Of these plants, seven were used as control plants and seven plants were sprayed with 100 μM GA3 (Duchefa Biochemie B. V, Haarlem, The Netherlands) once per week as described in [89]. Different morphological characters (leaf number, length/width of leaflets, petiole length, number of inflorescences, and number of main shoots) were measured every second day for 4 weeks.

Availability of data and materials

The datasets generated and/or analysed during the current study are available in the NCBI database repository, https://www.ncbi.nlm.nih.gov/search/all/?term=PRJNA561285

References

  1. 1.

    Isobe S, Klimenko I, Ivashuta S, Gau M, Kozlov NN. First RFLP linkage map of red clover (Trifolium pratense L.) based on cDNA probes and its transferability to other red clover germplasm. Theor Appl Genet. 2003;108:105–12. https://doi.org/10.1007/s00122-003-1412-z .

    CAS  Article  PubMed  Google Scholar 

  2. 2.

    Isobe S, Sawai A, Yamaguchi H, Gau M, Uchiyama K. Breeding potential of the backcross progenies of a hybrid between Trifolium medium × T. pratense to T. pratense. Can J Plant Sci. 2002;82:395–9. https://doi.org/10.4141/P01-034 .

    Article  Google Scholar 

  3. 3.

    Eriksen J, Askegaard M, Søegaard K. Complementary effects of red clover inclusion in ryegrass-white clover swards for grazing and cutting. Grass Forage Sci. 2014;69:241–50. https://doi.org/10.1111/gfs.12025 .

    CAS  Article  Google Scholar 

  4. 4.

    Young ND, Debellé F, Oldroyd GED, Geurts R, Cannon SB, Udvardi MK, et al. The Medicago genome provides insight into the evolution of rhizobial symbioses. Nature. 2011;480:520–4. https://doi.org/10.1038/nature10625 .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  5. 5.

    Sato S, Nakamura Y, Kaneko T, Asamizu E, Kato T, Nakao M, et al. Genome structure of the legume. Lotus japonicus Theor Appl Genet. 2008;15:227–39. https://doi.org/10.1093/dnares/dsn008 .

    CAS  Article  Google Scholar 

  6. 6.

    Schmutz J, Cannon SB, Schlueter J, Ma J, Mitros T, Nelson W, et al. Genome sequence of the palaeopolyploid soybean. Nature. 2010;463:178–83. https://doi.org/10.1038/nature08670 .

    CAS  Article  PubMed  Google Scholar 

  7. 7.

    Schmutz J, McClean PE, Mamidi S, Wu GA, Cannon SB, Grimwood J, et al. A reference genome for common bean and genome-wide analysis of dual domestications. Nat Genet. 2014;46:707–13. https://doi.org/10.1038/ng.3008 .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  8. 8.

    Varshney RK, Song C, Saxena RK, Azam S, Yu S, Sharpe AG, et al. Draft genome sequence of chickpea (Cicer arietinum) provides a resource for trait improvement. Nat Biotechnol. 2013;31:240–6. https://doi.org/10.1038/nbt.2491 .

    CAS  Article  Google Scholar 

  9. 9.

    Lonardi S, Muñoz-Amatriaín M, Liang Q, Shu S, Wanamaker SI, Lo S, et al. The genome of cowpea (Vigna unguiculata L. Walp.). Plant J. 2019;98:767–82. https://doi.org/10.1111/tpj.14349 .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  10. 10.

    Kaur P, Bayer PE, Milec Z, Vrána J, Yuan Y, Appels R, et al. An advanced reference genome of Trifolium subterraneum L. reveals genes related to agronomic performance. Plant Biotechnol J. 2017;15:1034–46. https://doi.org/10.1111/pbi.12697 .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  11. 11.

    Dluhošová J, Ištvánek J, Nedělník J, Řepková J. Red Clover (Trifolium pratense) and Zigzag Clover (T. medium) - A Picture of Genomic Similarities and Differences. Front Plant Sci. 2018;9:724. https://doi.org/10.3389/fpls.2018.00724 .

    Article  PubMed  PubMed Central  Google Scholar 

  12. 12.

    Ištvánek J, Jaros M, Krenek A, Řepková J. Genome assembly and annotation for red clover (Trifolium pratense; Fabaceae). Am J Bot. 2014;101:327–37. https://doi.org/10.3732/ajb.1300340 .

    Article  PubMed  Google Scholar 

  13. 13.

    JJ de V, Ayling S, Hegarty M, Kudrna D, Goicoechea JL, Ergon Å, et al. Red clover (Trifolium pratense L.) draft genome provides a platform for trait improvement. Sci Rep. 2015;5:17394. https://doi.org/10.1038/srep17394 .

    CAS  Article  Google Scholar 

  14. 14.

    Jahufer MZZ, Ford JL, Widdup KH, Harris C, Cousins G, Ayres JF, et al. Improving white clover for Australasia. Crop Pasture Sci. 2012;63:739. https://doi.org/10.1071/CP12142 .

    Article  Google Scholar 

  15. 15.

    Barrett BA, Faville MJ, Nichols SN, Simpson WR, Bryan GT, Conner AJ. Breaking through the feed barrier: options for improving forage genetics. Anim Prod Sci. 2015;55:883. https://doi.org/10.1071/AN14833 .

    Article  Google Scholar 

  16. 16.

    Řepková J, Nedělník J. Modern methods for genetic improvement of Trifolium pratense. Czech J Genet Plant Breed. 2014;50:92–9. https://doi.org/10.17221/139/2013-CJGPB .

    Article  Google Scholar 

  17. 17.

    Dias PMB, Julier B, Sampoux J-P, Barre P, Dall’Agnol M. Genetic diversity in red clover (Trifolium pratense L.) revealed by morphological and microsatellite (SSR) markers. Euphytica. 2008;160:189–205. https://doi.org/10.1007/s10681-007-9534-z .

    CAS  Article  Google Scholar 

  18. 18.

    Annicchiarico P, Proietti S. White clover selected for enhanced competitive ability widens the compatibility with grasses and favours the optimization of legume content and forage yield in mown clover-grass mixtures. Grass Forage Sci. 2010;140:no-no. https://doi.org/10.1111/j.1365-2494.2010.00749.x .

  19. 19.

    Ford JL, Barrett BA. Improving red clover persistence under grazing. Proceedings of the NZ Grassland Association. 2011;73:119–24.

    Article  Google Scholar 

  20. 20.

    Naydenova G, Hristova T, Aleksiev Y. Objectives and approaches in the breeding of perennial legumes for use in temporary pasturelands. Bio Anim Husb. 2013;29:233–50. https://doi.org/10.2298/BAH1302233N .

    Article  Google Scholar 

  21. 21.

    Tiffin P. Mechanisms of tolerance to herbivore damage:what do we know? Evol Ecol. 2000;14:523–36. https://doi.org/10.1023/A:1010881317261 .

    Article  Google Scholar 

  22. 22.

    Diaz S, Lavorel S, McIntyre SUE, Falczuk V, Casanoves F, Milchunas DG, et al. Plant trait responses to grazing ? A global synthesis. Glob Chang Biol. 2007;13:313–41. https://doi.org/10.1111/j.1365-2486.2006.01288.x .

    Article  Google Scholar 

  23. 23.

    van Minnebruggen A, Roldán-Ruiz I, van Bockstaele E, Haesaert G, Cnops G. The relationship between architectural characteristics and regrowth in Trifolium pratense (red clover). Grass Forage Sci. 2015;70:507–18. https://doi.org/10.1111/gfs.12138 .

    Article  Google Scholar 

  24. 24.

    Herbert DB, Ekschmitt K, Wissemann V, Becker A. Cutting reduces variation in biomass production of forage crops and allows low-performers to catch up: A case study of Trifolium pratense L. (red clover). Plant Biol (Stuttg). 2018;20:465–73. https://doi.org/10.1111/plb.12695 .

    CAS  Article  Google Scholar 

  25. 25.

    Conaghan P, Casler MD. A theoretical and practical analysis of the optimum breeding system for perennial ryegrass. Irish J Agricultural Food Res. 2011;50:47–63.

    Google Scholar 

  26. 26.

    Ortega F, Parra L, Quiroz A. Breeding red clover for improved persistence in Chile: a review. Crop Pasture Sci. 2014;65:1138. https://doi.org/10.1071/CP13323 .

    CAS  Article  Google Scholar 

  27. 27.

    Sato S, Isobe S, Asamizu E, Ohmido N, Kataoka R, Nakamura Y, et al. Comprehensive structural analysis of the genome of red clover (Trifolium pratense L.). DNA Res. 2005;12:301–64. https://doi.org/10.1093/dnares/dsi018 .

    CAS  Article  PubMed  Google Scholar 

  28. 28.

    Shimizu-Sato S, Tanaka M, Mori H. Auxin-cytokinin interactions in the control of shoot branching. Plant Mol Biol. 2009;69:429–35. https://doi.org/10.1007/s11103-008-9416-3 .

    CAS  Article  PubMed  Google Scholar 

  29. 29.

    Stafstrom J. Influence of bud position and plant ontogeny on the morphology of branch shoots in pea (Pisum sativum L. cv. Alaska). Ann Bot. 1995;76:343–8. https://doi.org/10.1006/anbo.1995.1106 .

    Article  Google Scholar 

  30. 30.

    Briske DD, Richards JH. Plant responses to defoliation: a physiological, morphological and demographicevaluation. In: Bedunah DJ, Sosebee RE, editors. Wildland plants: physiological ecology and developmental morphology. 1st ed. Society for Range Management: Denver, Colo; 1995. p. 635–710.

    Google Scholar 

  31. 31.

    Kotova LM, Kotov AA, Kara AN. Changes in Phytohormone status in stems and roots after decapitation of pea seedlings. Russ J Plant Physiol. 2004;51:107–11. https://doi.org/10.1023/B:RUPP.0000011309.47328.23 .

    CAS  Article  Google Scholar 

  32. 32.

    Li S, Strid Å. Anthocyanin accumulation and changes in CHS and PR-5 gene expression in Arabidopsis thaliana after removal of the inflorescence stem (decapitation). Plant Physiol Biochem. 2005;43:521–5. https://doi.org/10.1016/j.plaphy.2005.05.004 .

    CAS  Article  PubMed  Google Scholar 

  33. 33.

    Scholes DR, Wszalek AE, Paige KN. Regrowth patterns and rosette attributes contribute to the differential compensatory responses of Arabidopsis thaliana genotypes to apical damage. Plant Biol (Stuttg). 2016;18:239–48. https://doi.org/10.1111/plb.12404 .

    CAS  Article  Google Scholar 

  34. 34.

    Fischer M, Bossdorf O, Gockel S, Hänsel F, Hemp A, Hessenmöller D, et al. Implementing large-scale and long-term functional biodiversity research: the biodiversity Exploratories. Basic Applied Ecol. 2010;11:473–85. https://doi.org/10.1016/j.baae.2010.07.009 .

    Article  Google Scholar 

  35. 35.

    Ištvánek J, Dluhošová J, Dluhoš P, Pátková L, Nedělník J, Řepková J. Gene classification and Mining of Molecular Markers Useful in red clover (Trifolium pratense) breeding. Front Plant Sci. 2017. https://doi.org/10.3389/fpls.2017.00367 .

  36. 36.

    Pérez-Rodríguez P, Riaño-Pachón DM, Corrêa LGG, Rensing SA, Kersten B, Mueller-Roeber B. PlnTFDB: updated content and new features of the plant transcription factor database. Nucleic Acids Res. 2010;38:D822–7. https://doi.org/10.1093/nar/gkp805 .

    CAS  Article  PubMed  Google Scholar 

  37. 37.

    Yates SA, Swain MT, Hegarty MJ, Chernukin I, Lowe M, Allison GG, et al. De novo assembly of red clover transcriptome based on RNA-Seq data provides insight into drought response, gene discovery and marker identification. BMC Genomics. 2014;15:453. https://doi.org/10.1186/1471-2164-15-453 .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  38. 38.

    Asahina M, Satoh S. Molecular and physiological mechanisms regulating tissue Reunion in incised plant tissues. J Plant Res. 2015;128:381–8. https://doi.org/10.1007/s10265-015-0705-z .

    CAS  Article  PubMed  Google Scholar 

  39. 39.

    Pitaksaringkarn W, Ishiguro S, Asahina M, Satoh S. ARF6 and ARF8 contribute to tissue Reunion in incised Arabidopsis inflorescence stems. Plant Biotechnol. 2014;31:49–53. https://doi.org/10.5511/plantbiotechnology.13.1028b .

    CAS  Article  Google Scholar 

  40. 40.

    Pitaksaringkarn W, Matsuoka K, Asahina M, Miura K, Sage-Ono K, Ono M, et al. XTH20 and XTH19 regulated by ANAC071 under auxin flow are involved in cell proliferation in incised Arabidopsis inflorescence stems. Plant J. 2014;80:604–14. https://doi.org/10.1111/tpj.12654 .

    CAS  Article  PubMed  Google Scholar 

  41. 41.

    Balla J, Kalousek P, Reinöhl V, Friml J, Procházka S. Competitive canalization of PIN-dependent auxin flow from axillary buds controls pea bud outgrowth. Plant J. 2011;65:571–7. https://doi.org/10.1111/j.1365-313X.2010.04443.x .

    CAS  Article  PubMed  Google Scholar 

  42. 42.

    Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20. https://doi.org/10.1093/bioinformatics/btu170 .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  43. 43.

    Teichmann T, Muhr M. Shaping plant architecture. Front Plant Sci. 2015;6:233. https://doi.org/10.3389/fpls.2015.00233 .

    Article  PubMed  PubMed Central  Google Scholar 

  44. 44.

    Dar TA, Uddin M, Khan MMA, Hakeem KR, Jaleel H. Jasmonates counter plant stress: a review. Environ Exp Bot. 2015;115:49–57. https://doi.org/10.1016/j.envexpbot.2015.02.010 .

    CAS  Article  Google Scholar 

  45. 45.

    Schilmiller AL, Howe GA. Systemic signaling in the wound response. Curr Opin Plant Biol. 2005;8:369–77. https://doi.org/10.1016/j.pbi.2005.05.008 .

    CAS  Article  PubMed  Google Scholar 

  46. 46.

    Wasternack C. Action of jasmonates in plant stress responses and development--applied aspects. Biotechnol Adv. 2014;32:31–9. https://doi.org/10.1016/j.biotechadv.2013.09.009 .

    CAS  Article  PubMed  Google Scholar 

  47. 47.

    Huang H, Liu B, Liu L, Song S. Jasmonate action in plant growth and development. J Exp Bot. 2017;68:1349–59. https://doi.org/10.1093/jxb/erw495 .

    CAS  Article  PubMed  Google Scholar 

  48. 48.

    Mur LAJ, Prats E, Pierre S, Hall MA, Hebelstrup KH. Integrating nitric oxide into salicylic acid and jasmonic acid/ ethylene plant defense pathways. Front Plant Sci. 2013;4:215. https://doi.org/10.3389/fpls.2013.00215 .

    Article  PubMed  PubMed Central  Google Scholar 

  49. 49.

    Lee SC, Luan S. ABA signal transduction at the crossroad of biotic and abiotic stress responses. Plant Cell Environ. 2012;35:53–60. https://doi.org/10.1111/j.1365-3040.2011.02426.x .

    Article  Google Scholar 

  50. 50.

    Kebrom TH, Spielmeyer W, Finnegan EJ. Grasses provide new insights into regulation of shoot branching. Trends Plant Sci. 2013;18:41–8. https://doi.org/10.1016/j.tplants.2012.07.001 .

    CAS  Article  PubMed  Google Scholar 

  51. 51.

    Wang Y, Zhao J, Lu W, Deng D. Gibberellin in plant height control: old player, new story. Plant Cell Rep. 2017;36:391–8. https://doi.org/10.1007/s00299-017-2104-5 .

    CAS  Article  PubMed  Google Scholar 

  52. 52.

    Morris SE, Cox MCH, Ross JJ, Krisantini S, Beveridge CA. Auxin dynamics after decapitation are not correlated with the initial growth of axillary buds. Plant Physiol. 2005;138:1665–72. https://doi.org/10.1104/pp.104.058743 .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  53. 53.

    Domagalska MA, Leyser O. Signal integration in the control of shoot branching. Nat Rev Mol Cell Biol. 2011;12:211–21. https://doi.org/10.1038/nrm3088 .

    CAS  Article  PubMed  Google Scholar 

  54. 54.

    Dun EA, de Saint Germain A, Rameau C, Beveridge CA. Dynamics of strigolactone function and shoot branching responses in Pisum sativum. Mol Plant. 2013;6:128–40. https://doi.org/10.1093/mp/sss131 .

    CAS  Article  PubMed  Google Scholar 

  55. 55.

    Umehara M, Hanada A, Yoshida S, Akiyama K, Arite T, Takeda-Kamiya N, et al. Inhibition of shoot branching by new terpenoid plant hormones. Nature. 2008;455:195–200. https://doi.org/10.1038/nature07272 .

    CAS  Article  PubMed  Google Scholar 

  56. 56.

    Borghi L, Liu G-W, Emonet A, Kretzschmar T, Martinoia E. The importance of strigolactone transport regulation for symbiotic signaling and shoot branching. Planta. 2016;243:1351–60. https://doi.org/10.1007/s00425-016-2503-9 .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  57. 57.

    Achard P, Gusti A, Cheminant S, Alioua M, Dhondt S, Coppens F, et al. Gibberellin signaling controls cell proliferation rate in Arabidopsis. Curr Biol. 2009;19:1188–93. https://doi.org/10.1016/j.cub.2009.05.059 .

    CAS  Article  PubMed  Google Scholar 

  58. 58.

    DeMason DA, Chetty VJ. Interactions between GA, auxin, and UNI expression controlling shoot ontogeny, leaf morphogenesis, and auxin response in Pisum sativum (Fabaceae): or how the uni-tac mutant is rescued. Am J Bot. 2011;98:775–91. https://doi.org/10.3732/ajb.1000358 .

    CAS  Article  PubMed  Google Scholar 

  59. 59.

    Nelissen H, Gonzalez N, Inzé D. Leaf growth in dicots and monocots: so different yet so alike. Curr Opin Plant Biol. 2016;33:72–6. https://doi.org/10.1016/j.pbi.2016.06.009 .

    Article  PubMed  Google Scholar 

  60. 60.

    Guttridge CG, Thombson PA. The effect of gibberellins on growth and flowering of Fragaria and Duchesnea. J Exp Bot. 1964;15:631–46. https://doi.org/10.1093/jxb/15.3.631 .

    CAS  Article  Google Scholar 

  61. 61.

    Leite VM, Rosolem CA, Rodrigues JD. Gibberellin and cytokinin effects on soybean growth. Sci. agric. (Piracicaba, Braz.). 2003;60:537–41. https://doi.org/10.1590/S0103-90162003000300019 .

    CAS  Article  Google Scholar 

  62. 62.

    Tsukaya H, Kozuka T, Kim G-T. Genetic control of petiole length in Arabidopsis thaliana. Plant Cell Physiol. 2002;43:1221–8. https://doi.org/10.1093/pcp/pcf147 .

    CAS  Article  PubMed  Google Scholar 

  63. 63.

    Hisamatsu T, King RW, Helliwell CA, Koshioka M. The involvement of gibberellin 20-oxidase genes in phytochrome-regulated petiole elongation of Arabidopsis. Plant Physiol. 2005;138:1106–16. https://doi.org/10.1104/pp.104.059055 .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  64. 64.

    Cai Y, Shao L, Li X, Liu G, Chen S. Gibberellin stimulates regrowth after defoliation of sheepgrass (Leymus chinensis) by regulating expression of fructan-related genes. J Plant Res. 2016;129:935–44. https://doi.org/10.1007/s10265-016-0832-1 .

    CAS  Article  PubMed  Google Scholar 

  65. 65.

    Morvan-Bertrand A, Ernstsen A, Lindgard B, Koshioka M, Le Saos J, Boucaud J, et al. Endogenous gibberellins in Lolium perenne and influence of defoliation on their contents in elongating leaf bases and in leaf sheaths. Physiol Plant. 2001;111:225–31. https://doi.org/10.1034/j.1399-3054.2001.1110214.x .

    CAS  Article  Google Scholar 

  66. 66.

    Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, et al. De novo transcript sequence reconstruction from RNA-seq using the trinity platform for reference generation and analysis. Nat Protoc. 2013;8:1494–512. https://doi.org/10.1038/nprot.2013.084 .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  67. 67.

    Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29:644–52. https://doi.org/10.1038/nbt.1883 .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  68. 68.

    Schulz MH, Zerbino DR, Vingron M, Birney E. Oases: robust de novo RNA-seq assembly across the dynamic range of expression levels. Bioinformatics. 2012;28:1086–92. https://doi.org/10.1093/bioinformatics/bts094 .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  69. 69.

    Bekel T, Henckel K, Küster H, Meyer F, Mittard Runte V, Neuweger H, et al. The Sequence Analysis and Management System – SAMS-2.0: Data management and sequence analysis adapted to changing requirements from traditional sanger sequencing to ultrafast sequencing technologies. J Biotechnol. 2009;140:3–12. https://doi.org/10.1016/j.jbiotec.2009.01.006 .

    CAS  Article  PubMed  Google Scholar 

  70. 70.

    Boutet E, Lieberherr D, Tognolli M, Schneider M, Bairoch A. UniProtKB/Swiss-Prot. Methods Mol Biol. 2007;406:89–112.

    CAS  PubMed  Google Scholar 

  71. 71.

    Bairoch A, Apweiler R. The SWISS-PROT protein sequence data bank and its supplement TrEMBL. Nucleic Acids Res. 1997;25:31–6. https://doi.org/10.1093/nar/25.1.31 .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  72. 72.

    Goodstein DM, Shu S, Howson R, Neupane R, Hayes RD, Fazo J, et al. Phytozome: a comparative platform for green plant genomics. Nucleic Acids Res. 2012;40:D1178–86. https://doi.org/10.1093/nar/gkr944 .

    CAS  Article  PubMed  Google Scholar 

  73. 73.

    Wu TD, Watanabe CK. GMAP: a genomic mapping and alignment program for mRNA and EST sequences. Bioinformatics. 2005;21:1859–75. https://doi.org/10.1093/bioinformatics/bti310 .

    CAS  Article  PubMed  Google Scholar 

  74. 74.

    Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323. https://doi.org/10.1186/1471-2105-12-323 .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  75. 75.

    Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550. https://doi.org/10.1186/s13059-014-0550-8 .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  76. 76.

    Wagner GP, Kin K, Lynch VJ. Measurement of mRNA abundance using RNA-seq data: RPKM measure is inconsistent among samples. Theory Biosci. 2012;131:281–5. https://doi.org/10.1007/s12064-012-0162-3 .

    CAS  Article  Google Scholar 

  77. 77.

    The UniProt Consortium. UniProt: the universal protein knowledgebase. Nucleic Acids Res. 2016;45:D158–69. https://doi.org/10.1093/nar/gkw1099 .

    CAS  Article  Google Scholar 

  78. 78.

    Resource NCBI. Coordinators. Database resources of the National Center for biotechnology information. Nucleic Acids Res. 2016;44:D7–19. https://doi.org/10.1093/nar/gkv1290 .

    CAS  Article  Google Scholar 

  79. 79.

    Berardini TZ, Reiser L, Li D, Mezheritsky Y, Muller R, Strait E, Huala E. The Arabidopsis information Resource: making and mining the ‘gold standard’ annotated reference plant genome. Genesis. 2015;53:474–85. https://doi.org/10.1002/dvg.22877 .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  80. 80.

    Camon E, Barrell D, Brooksbank C, Magrane M, Apweiler R. The gene ontology annotation (GOA) project--application of GO in SWISS-PROT, TrEMBL and InterPro. Comp Funct Genomics. 2003;4:71–4. https://doi.org/10.1002/cfg.235 .

    Article  PubMed  PubMed Central  Google Scholar 

  81. 81.

    Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. Gene Ontol Consortium Nat Genet. 2000;25:25–9. https://doi.org/10.1038/75556 .

    CAS  Article  Google Scholar 

  82. 82.

    The Gene Ontology Consortium. Expansion of the gene ontology knowledgebase and resources. Nucleic Acids Res. 2016;45:D331–8. https://doi.org/10.1093/nar/gkw1108 .

    CAS  Article  PubMed Central  Google Scholar 

  83. 83.

    Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11:R14. https://doi.org/10.1186/gb-2010-11-2-r14 .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  84. 84.

    Walter W, Sánchez-Cabo F, Ricote M. GOplot: an R package for visually combining expression data with functional analysis. Bioinformatics. 2015;31:2912–4. https://doi.org/10.1093/bioinformatics/btv300 .

    CAS  Article  PubMed  Google Scholar 

  85. 85.

    Racine JS. RStudio: a platform-independent IDE for R and Sweave. J Appl Econ. 2012;27:167–72. https://doi.org/10.1002/jae.1278 .

    Article  Google Scholar 

  86. 86.

    R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. 2020.

    Google Scholar 

  87. 87.

    Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215:403–10. https://doi.org/10.1016/S0022-2836(05)80360-2 .

    CAS  Article  Google Scholar 

  88. 88.

    Götz S, García-Gómez JM, Terol J, Williams TD, Nagaraj SH, Nueda MJ, et al. High-throughput functional annotation and data mining with the Blast2GO suite. Nucleic Acids Res. 2008;36:3420–35. https://doi.org/10.1093/nar/gkn176 .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  89. 89.

    Blazquez MA, Green R, Nilsson O, Sussman MR, Weigel D. Gibberellins promote flowering of arabidopsis by activating the LEAFY promoter. Plant Cell. 1998;10:791–800. https://doi.org/10.1105/tpc.10.5.791 .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We thank Andrea Weisert for excellent technical help and Dietmar Haffer for skillfully raising the plants. We also thank Volker Wissemann and Birgit Gemeinholzer for continuous discussions on the project. Bioinformatics support by the BMBF-funded project “Bielefeld-Gießen Center for Microbial Bioinformatics—BiGi (Grant Number 031A533)” within the German Network for Bioinformatics Infrastructure (de. NBI) is gratefully acknowledged. This work was funded by a grant from the German Research Foundation (DFG, BE 2547/12–1) to A. B.

Third party submission

The manuscript is not in submission elsewhere.

Funding

This work was funded by a grant from the German Research Foundation (DFG, BE 2547/12–1, − 2) to A. B.. Open Access funding enabled and organized by Projekt DEAL.

Author information

Affiliations

Authors

Contributions

DH acquired, analyzed, and interpreted the data and drafted the manuscript, TG analyzed and interpreted the data a drafted the manuscript, OR analyzed the data, AB conceived and designed the work, analyzed the data, and revised the manuscript All authors read and approved the submitted version.

Corresponding author

Correspondence to Annette Becker.

Ethics declarations

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable.

Competing interests

Not applicable.

Additional information

Publisher’s Note

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

Supplementary Information

Additional file 1: Figure S1.

Map of sample locations within the Biodiversity Exploratories.

Additional file 2: Figure S2.

Annotation Overview: A: Distribution of transcripts that could be mapped to the T. pratense genome, to a known locus and were annotated with T. pratense genome identifier. B: Distribution of transcripts that could be mapped to an unknown T. pratense gene locus. C: Distribution of transcripts that could not be mapped o the T. pratense genome. D: Distribution of transcripts of whole transcriptome representing all 12 libraries.

Additional file 3: Figure S3.

qRT-PCR analysis of selected T. pratense genes: A) tdn_146439/ENGase85A (encodes a cytosolic beta-endo-N-acetyglucosaminidase (ENGase), B) tdn_69411/LTP (Lipid transfer protein, a pathogenesis-related protein), C) tdn85889/PME44 (pectin methylesterase 44), D) K65_9861 P5CS (DELTA1-PYRROLINE-5-CARBOXYLATE SYNTHASE 1). Two additional T. pratense tissues (axial meristem (AM) and leaves (L)) were included and also the six RNA extractions also used for constructing the transcriptome libraries. Transcriptome data are marked with TPM, all samples are from mown (M) or not mown (NM) plants. Gene expression levels were normalized to the housekeeping contig k65_5754 (not annotated, but with similar expression throughout the transcriptomes). Shown are mean values of the fold change of the respective contig in relation to the expression of contig K65_5754 with error bars that represent standard deviation. On the right side of each graph, the log2 values of the respective gene in each transcriptome are shown for better comparison of expression patterns.

Additional file 4: Figure S4.

Enrichment analysis of the DEG GO terms. A: Enriched GO terms in mown and non-mown greenhouse samples. B: Enriched GO terms in mown and non-mown field B samples. C: Enriched GO terms in mown and non-mown field A samples. Information about the respective GO number can be found in Table S11.

Additional file 5: Figure S5.

Expressed transcription factor family members in mown and not mown T. pratense plants. The y-axis shows the number of upregulated transcription factors that are members of the respective transcription factor family. Names of the transcriptomes (GM, GNM, FaM, FaNM; FbM, FbNM) and transcription factor families are given on the x-axis. Expression of transcription factor members were compared in a pairwise manner (GM vs GNM, FaM vs FaNM, FbM vs FbNM). Comparisons that resulted in a difference of more than 10% of the contigs significantly upregulated in either the mown or the unmown condition were marked red, differences between 5 and 9 % were marked orange.

Additional file 6: Figure S6.

Plant architectural characteristics and growth habit of gibberellin treated plants. A-E: Measured, counted or calculated plant characteristics during phenotypic monitoring experiments. Gibberellin treated plants, blue; control plants, orange. Graphs show average values and 95% confidence intervals. Time is shown in weeks. Growth habit of control plants (left side) vs. gibberellin treated plants (right side), after approximately 2 weeks of gibberellin treatment and regrowth (F), and after 4 weeks (G).

Additional file 7: Table S1.

Overview of the sampling locations for the plant material. Names of the fields belonging to the Biodiversity Exploratory or greenhouse populations are shown. As well as the location, coordinates and conditions (mown/cut and not mown/uncut). Table S3. Sources and basis for description and classification of the top 20 DEG. Table S4. Number of reads for each sequenced library (transcriptome ID) before and after trimming. Table S5. General features of the transcriptome of T. pratense. Table S6. Overall alignment rate of the single transcriptomes to the references transcriptome, values above 80% are good. Table S7. Annotation of T. pratense plant-specific against different databases. Table S8. Quality of the replicates. For each library replicate the number of transcripts above TPM of 1 is shown. Further the number of transcripts shared between to related replica is shown, as well as the number of transcripts unique for a replica. The similarity of the replicas was evaluated by calculating he percentage of the shared transcripts compared to the total number of transcripts of a replica. Table S10. Main classes based on the DE contigs. 16 main classes were developed to group the DE contigs. Table S11. GO Terms specifically enriched in the individual transcriptomes. Table S12. Shared contigs with corresponding annotation. Table S13. Detailed information about transcripts described during the discussion part. Table S14. Primer sequences.

Additional file 8: Table S2.

Annotation of the T. pratense transcriptome.

Additional file 9: Table S9.

Analysis of the transcriptomes with Deseq2 of all transcriptome comparisons.

Additional file 10: Table S16.

TPM values for T. pratense transcripts.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. 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 in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Herbert, D.B., Gross, T., Rupp, O. et al. Transcriptome analysis reveals major transcriptional changes during regrowth after mowing of red clover (Trifolium pratense). BMC Plant Biol 21, 95 (2021). https://doi.org/10.1186/s12870-021-02867-0

Download citation

Keywords

  • Trifolium pretense
  • Red clover
  • RNAseq
  • Regrowth reaction
  • Biotic and abiotic stress
  • Field conditions
  • Gibberellins