Transcriptome analysis reveals major transcriptional changes during regrowth after mowing of red clover (Trifolium pratense)
BMC Plant Biology volume 21, Article number: 95 (2021)
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.
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.
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.
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) , Lotus japonicus (birdsfoot trefoil) , Glycine max (soy) , Phaseolus vulgaris (common bean) , Cicer arietinum (chickpea) , Vigna unguiculata (cowpea) , Trifolium subterraneum (subterranean clover)  T. medium (zigzag clover) , 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 . 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 . 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  and is a complex trait influenced by a variety of abiotic and biotic factors, and the regrowth ability of a plant . 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 . 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 . 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”  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 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).
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 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).
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 . 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.
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 . 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.
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 . 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) ). 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 ). 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) . 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) . 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 . 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 . 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 . 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 . 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 .
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 , 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 . 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 , 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 ). 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 . 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.
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” , 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 . 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 , 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  with ILLUMINACLIP, SLIDINGWINDOW:5:20 and MINLEN:50 options. Quality trimmed reads were pooled and digitally normalized . Multiple de novo assemblies were computed using Trinity  and Oases  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)  for functional annotation with the SwissProt , TrEMBL  and Phytozome  (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 . 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) (, 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 . 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  with bowtie mapping. To identify differentially expressed (DEGs) T. pratense genes, a pairwise comparison of all treatments was preformed using the DESeq2  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 . TPM (transcript per million) values were calculated to estimate contig expression level  (for RSEM read counts see Table S16).
We used the description and gene names obtained from TrEMBLE and SwissProt to search the UniProt , NCBI  and TAIR  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  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  implemented in RStudio  with the program R .
Two local BLAST searches  with word-size of 3, e-value of 1.0e-3 and HSP length cutoff of 33 were performed against the PlnTFDB using Blast2GO . 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.
To assess the effect of gibberellin during the regrowth reaction of T. pratense, 14 red clover plants were mown as described in . 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 . 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
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
Ř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 .
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 .
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 .
Ford JL, Barrett BA. Improving red clover persistence under grazing. Proceedings of the NZ Grassland Association. 2011;73:119–24.
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 .
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 .
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 .
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 .
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 .
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.
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 .
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 .
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 .
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 .
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.
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
Teichmann T, Muhr M. Shaping plant architecture. Front Plant Sci. 2015;6:233. https://doi.org/10.3389/fpls.2015.00233 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
Boutet E, Lieberherr D, Tognolli M, Schneider M, Bairoch A. UniProtKB/Swiss-Prot. Methods Mol Biol. 2007;406:89–112.
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 .
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 .
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 .
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 .
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 .
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 .
The UniProt Consortium. UniProt: the universal protein knowledgebase. Nucleic Acids Res. 2016;45:D158–69. https://doi.org/10.1093/nar/gkw1099 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. 2020.
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 .
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 .
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 .
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.
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.
Ethics approval and consent to participate
Consent for publication
Map of sample locations within the Biodiversity Exploratories.
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.
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.
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.
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.
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).
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.
Annotation of the T. pratense transcriptome.
Analysis of the transcriptomes with Deseq2 of all transcriptome comparisons.
TPM values for T. pratense transcripts.
About this article
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