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

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. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-021-02867-0.

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.

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 nonmown (FbM vs. FbNM); greenhouse mown vs. greenhouse non-mown (GM vs. GNM) to identify DEGs (Table S9) from mown plants. Interestingly, using the |log2fold-change| ≥ 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).

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.
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).
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 transcriptionand 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 Leucinrich 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 Growth stress  Biotic stress   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.
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 fieldderived 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 C 2 C 2 -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 GA 3 .
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 cm 2 almost double the size when compared with those of untreated plants (2.4 ± 0.6 cm 2 ). 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.

RNA-Seq and assembly
The de novo assembly in combination with a referencebased approach for the annotation led to the identification of 44,643 contigs of which 29,781 were annotated as plantspecific (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 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 regrowthspecific 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 unmown 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 H 2 O, 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, evalue 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 GA 3 (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.