Deciphering Trifolium pratense L. holobiont reveals a microbiome resilient to future climate changes

Abstract The plant microbiome supports plant growth, fitness, and resistance against climate change. Trifolium pratense (red clover), an important forage legume crop, positively contributes to ecosystem sustainability. However, T. pratense is known to have limited adaptive ability toward climate change. Here, the T. pratense microbiomes (including both bacteria and fungi) of the rhizosphere and the root, shoot, and flower endospheres were comparatively examined using metabarcoding in a field located in Central Germany that mimics the climate conditions projected for the next 50–70 years in comparison with the current climate conditions. Additionally, the ecological functions and metabolic genes of the microbial communities colonizing each plant compartment were predicted using FUNGuild, FAPROTAX, and Tax4Fun annotation tools. Our results showed that the individual plant compartments were colonized by specific microbes. The bacterial and fungal community compositions of the belowground plant compartments did not vary under future climate conditions. However, future climate conditions slightly altered the relative abundances of specific fungal classes of the aboveground compartments. We predicted several microbial functional genes of the T. pratense microbiome involved in plant growth processes, such as biofertilization (nitrogen fixation, phosphorus solubilization, and siderophore biosynthesis) and biostimulation (phytohormone and auxin production). Our findings indicated that T. pratense microbiomes show a degree of resilience to future climate changes. Additionally, microbes inhabiting T. pratense may not only contribute to plant growth promotion but also to ecosystem sustainability.


| INTRODUC TI ON
Forage legume crops with high protein and fiber contents are a major livestock feed source. The integration of forage legumes into the cropping systems can have beneficial effects on soil health and fertility, as well as on controlling weeds, insect pests, and pathogens (Sheaffer & Seguin, 2008). Trifolium pratense L. (red clover), a forage legume crop in the temperate regions, is a key component of sustainable livestock farming systems (De Vega et al., 2015). In the 16 th century, T. pratense was used as a protein-rich fodder in livestock agriculture. T. pratense was further used as a "nitrogen-assimilating crop" in the 19 th century when the soil nitrogen content depleted in Europe (Kjaergaard, 2003;McKenna et al., 2018). Red clover efficiently fixes atmospheric nitrogen (N) due to its symbiotic association with N-fixing bacteria (Fustec et al., 2010). Additionally, the use of red clover increases soil fertility through the rhizodeposition of plant exudates containing soluble N compounds (Paynel et al., 2008). The decomposition of red clover residues results in the release of 40-70% of the total plant N into the soil within 5-10 weeks of decomposition (Lupwayi et al., 2006). Therefore, red clover is considered a "fertility-building crop" (McKenna et al., 2018). The incorporation of red clover in agricultural crop rotations is an alternative and sustainable method of introducing N into low-input agricultural practices. In addition to its application in agriculture, red clover has pharmacological applications as it exhibits oestrogenic, antispasmodic, and expectorant properties (Coon et al., 2007;Leung & Foster, 1996;Lin et al., 2000).
The events associated with climate change, including increased global temperatures and altered precipitation patterns, adversely affect plant health and productivity across different agroecosystems (Franklin et al., 2016;Schmidhuber & Tubiello, 2007). Recent studies have suggested that climate change has led to shifts in plant phenology, species distribution, and population dynamics and has contributed to the emergence of new potential fungal plant pathogens (Delgado-Baquerizo et al., 2020;Franklin et al., 2016;Wahdan et al., 2020). T. pratense is adapted to a wide range of soil types and pH levels in temperate regions. However, it has a limited capacity to adapt to increased temperatures and extreme drought events (Hanna et al., 2018). Previous studies have reported that red clover is resistant to a maximum temperature of 25°C but that prolonged exposure to 28°C decreases the crop yield (Hanna et al., 2018). Additionally, red clover exhibited some resistance to moderate drought, however, extreme drought highly impaired the yield that did not recover after a postdrought period (Hofer et al., 2016). Various studies have examined the cascading effects of climate change on T. pratense performance.
However, the response of the T. pratense microbiome to climate changes has not been examined.
The plant holobiont, which comprises the host plant and its endocellular and extracellular microbiome (Rosenberg & Zilber-Rosenberg, 2018), is considered a biological entity associated with stability, adaptation, and evolution, and not as individual biotic components (Vandenkoornhuyse et al., 2015). The host plant traits, such as resistance against pathogens, immune system priming, and growth, are dependent on the host's microbiome composition (Hartmann et al., 2007;Mendes et al., 2011;Ritpitakphong et al., 2016;Trivedi et al., 2020). In contrast to the highly conserved plant genome, the microbiome genome is prone to rapid genetic changes (Rosenberg & Zilber-Rosenberg, 2018). Therefore, the plasticity of the microbiome to adapt to environmental changes enables rapid host adaptation (Voolstra & Ziegler, 2020). Microbiome plasticity is a broad phenomenon that includes a dynamic reconstruction of microbiome composition by increasing and/or decreasing the abundance of specific microbes and/or by the colonization of novel microbes that facilitate the host adaptation to external stress (Bulgarelli et al., 2013;Haney et al., 2015). However, enhanced microbiome plasticity increases the risk of pathogen invasion and undesirable microbes enrichment with a concomitant loss of beneficial ones (Voolstra & Ziegler, 2020). Beneficial microbiome plasticity depends on the dynamics within useful microbes that maintain high levels of functional redundancy in the original microbial communities. In another scenario, the microbiome may respond to environmental changes by exhibiting resistance or by maintaining a constant community structure with a high potential to adapt to external stress (Allison & Martiny, 2008). The plasticity or resistance of host-associated microbiomes may contribute to host adaptation. Nevertheless, the adaptive strategies employed by the T. pratense microbiome in response to future climate conditions are so far unclear.
Within the host plant, microbial communities vary between the belowground and aboveground plant compartments, which are distinct ecological niches with variations in nutrients and oxygen availability in different tissue types (Beckers et al., 2017;Cregger et al., 2018;Pangesti et al., 2020;Zarraonaindia et al., 2015).
Microorganisms reach their host to form the indigenous microbiome through the following two pathways: vertical transmission via seeds and horizontal transmission from the surrounding atmosphere, rhizosphere, and bulk soils (Trivedi et al., 2020). The rhizosphere is the soil zone around the roots in which microbes are impacted by the presence of plant roots (Vandenkoornhuyse et al., 2015). The density of microbial populations in the rhizosphere is higher than that in the bulk soils; therefore, it is considered a hot spot for plant-microbiome interaction (Berendsen et al., 2012). The microbiome composition of the root endosphere depends on the ability of microbes to invade root tissues from the surrounding rhizosphere and rhizoplane (Pangesti et al., 2020;Vandenkoornhuyse et al., 2015). Soil is also a microbial reservoir for the aboveground plant compartments as some endophytic microbes (microbes that colonize the internal plant tissue showing no infection or negative effect on their host; Schulz and Boyle (2006)) of the aboveground plant compartments/niches are recruited from soil (Zarraonaindia et al., 2015). Additionally, the aboveground endophytic microbiomes are derived from microbes that first colonize the leaf and flower surfaces as epiphytes (Vandenkoornhuyse et al., 2015) and can passively or actively invade the plant tissues (De Vrieze et al., 2018).
The ability of the host plants to utilize beneficial microbes, such as plant growth-promoting bacteria (PGPB) determines their response to the environmental and climate changes through direct and indirect mechanisms. These mechanisms include nutrient solubilization, biological nitrogen fixation, and the production of plant growth regulators, organic acids, and volatile organic compounds (Ahkami et al., 2017). Therefore, the identification of the T. pratense microbiome functional profile is critical for developing new strategies to enhance plant health, growth, and resistance against future climate changes.
This study aimed to investigate the responses of the bacteriome and mycobiome (i.e., bacterial and fungal microbiomes) associated with four ecological niches/compartments of T. pratense and to evaluate their potential ecological and metabolic functions in responding to future climate conditions. The rhizosphere and the endospheres of the root, shoot system (leaves and stems), and flower were examined. The study was performed at grassland plots of the Global Change Experimental Facility (GCEF) established in central Germany . GCEF is one of the largest experimental platforms designed to investigate the effect of a future climate scenario mimicking the prediction for the next 50-70 years on ecosystem processes in plots under different land-uses .
The sampling was performed 4 years after starting the climate manipulation in summer as it represents the critical season in which the future climate scenario is expected to have the highest impacts on soil functions (Yin et al., 2019). The period of 4 years after the onset of the experiment was sufficient for T. pratense generation and their microbiome to be affected by climate manipulation and adapt through the vertical and horizontal transmission of new microorganisms. MiSeq Illumina sequencing of the 16S rRNA gene (V5-V7 region) and the nuclear ribosomal internal transcribed spacer region 2 (ITS2) was performed to characterize the bacterial and fungal microbiomes, respectively. We hypothesized that T. pratense-associated microbiomes would be shaped by the influence of both biotic (plant compartments/ecological niches) and abiotic (climate change) factors that varied in their relative importance.

| Study site and experimental design
The study was conducted in GCEF at the field research station of the Helmholtz Centre for Environmental Research in Bad Lauchstädt, Saxony-Anhalt, Germany (51_22060 N, 11_50060 E, 118 m a.s.l.). The area is characterized by a subcontinental climate (mean temperature, 8.9°C and mean annual rainfall, 498 mm for the period 1896-2013; mean temperature, 9.8°C and mean annual rainfall, 516 mm for the period 1995-2014). During the study period (2018), the mean temperature was 10.8°C with an annual rainfall of 254 mm. The study field comprised the Haplic Chernozem soil, which was characterized by a high content of organic carbon till a depth of more than 40 cm and a high water-holding capacity (Altermann et al., 2005).
The GCEF field infrastructure ( Figure A1) was designed to comparatively investigate the consequences of future climate and current climate conditions on ecosystem processes in different land types . Furthermore, the GCEF comprises 50 field plots (400 m 2 each), which were equally divided and subjected to the current and future climate conditions. Future climate condition is a consensus scenario across three models (COSMO-CLM (Rockel et al., 2008), REMO (Jacob & Podzun, 1997), and RCAO (Döscher et al., 2002)) of climate change in Central Germany for the years 2070-2100. Hence, future climate plots ( Figure A2) are equipped with mobile shelters and side panels, as well as an irrigation system. The roofs are controlled by a rain sensor. The continuous adjustment of irrigation or roof closing has decreased the precipitation by approximately 20% in the summer months and increased the precipitation by approximately 10% in spring and autumn. To simulate the increase in temperature, the standard method "passive night-time warming" was used. The shelters and panels were automatically closed from sundown to sunrise to increase the mean daily temperature by approximately 0.55°C accompanied by a stronger increase in minimum temperatures (up to 1.14°C on average). Current climate plots are equipped with the same steel constructions (but without shelters, panels, and irrigation systems) to mimic the possible microclimatic effects of the experimental setup. The resulting changes in climate conditions, due to climate manipulation, before and during the study period are shown in Figure A3. For more details on the field station design, see Schädler et al. (2019). The experiment was performed in the extensively used meadow plots subjected to future climate conditions (5 plots) in comparison with the plots of current climate conditions (5 plots). The vegetation comprises 56 plant species that were chosen from multiple regional natural source populations located in Central Germany. Each source population is genetically different. T. pratense species is represented by 2 gene pools (Madaj et al., 2020). The vegetation was mowed twice a year without the application of herbicides or fertilizers. The experiment was conducted in mid-July 2018 (summer), which corresponded with the highest effect of future climate conditions on soil ecosystem function (plant residue decomposition) at the GCEF in other years (Yin et al., 2019).

| Sample collection and compartmentalization of the belowground and aboveground plant compartments
Each climate scenario was represented by five plots. At each plot, three healthy T. pratense L. (red clover) plants were randomly selected and their two belowground compartments (rhizosphere soil and root) and two aboveground compartments (leaf/stem and flower) were examined. In total, 30 plants (3 plants × 10 plots) were sampled, the two halves of which are representing current and future climate scenarios. For each plant, the bulk soil was removed by vigorous shaking for 10 min. The adhering rhizosphere soil was collected by vortexing the roots for 10 min in a sterile polymerase chain reaction (PCR) water (Barillot et al., 2012). The root was separated from the aboveground compartments and surface-sterilized to collect the endophytes. Briefly, the root was washed under running distilled water, followed by three washes with 0.1% Tween 20, a 3 min wash with 70% ethanol, and five washes with sterilized distilled water. Similarly, the endophytes were obtained from the aboveground compartments after surface sterilization. The leaves and stems were considered as one compartment, while the flowers were considered a separate compartment. The two compartments were washed twice with 0.1% Tween 20, followed by five washes with sterilized distilled water. The samples from the three plants of each plot were pooled into a single composite sample. The entire sterilized compartments (root, leaves/stems, and flowers) were crushed using liquid nitrogen and the resulting powder was used for DNA extraction.
The ITS2 region of fungi was amplified using the following primers: fITS7F forward (5′-GTGARTCATCGAATCTTTG-3′) (White et al., 1990) and ITS4 reverse (5′-TCCTC CGCTTATTGATATGC-3′) (White et al., 1990). The amplification was performed in a two-step process. The forward primer of the first PCR was constructed using the Illumina i5 sequencing primer (5′-TCGTCGGCAGCGTCAGATGTG TATAAGA GACAG-3′) and a specific forward primer. The reverse primer was constructed using the Illumina i7 sequencing primer (5′-GTCTCGTGGG CTCGGAGATGTGTATAAGAGACAG-3′) and the specific reverse primer. The amplification was performed in a 25 μl reaction volume comprising 1 μl (5 μM) of each primer and 1 μl of the template using the Qiagen HotStar hi-fidelity polymerase kit (Qiagen Inc.). PCR was performed using an ABI Veriti thermocycler (Applied Biosystems). The PCR conditions were as follows: 95°C for 5 min, followed by 35 cycles of 94°C for 15 s, 54°C for 60 s, and 72°C for 1 min and one step of 72°C for 10 min and 4°C hold. The amplicons from the first PCR, whose concentrations were quantitatively determined, were used for the second PCR. In the second PCR, dual indices were attached using the Nextera XT index kit. The conditions for the second PCR were the same as those used for the first PCR, except for the amplification cycles (10 amplification cycles used in the second PCR). The amplicons were visualized using eGels (Life Technologies), following the manufacturer's instructions. Equimolar concentrations of the products were pooled, and the size of each pool was selected in two rounds using Agencourt AMPure XP (BeckmanCoulter) in a 0.75 ratio for both rounds. The size-selected pools were then quantified using a Quibit 2.0 fluorometer (Life Technologies). Sequencing was performed using MiSeq (Illumina, Inc) with a 2 × 300 bp paired-end strategy, following the manufacturer's instructions.

| Processing of amplicon data
The primer sequences were trimmed from the demultiplexed raw reads using cutadapt (Martin, 2011). The pair-end raw reads of bacterial and fungal datasets were merged using the simple Bayesian algorithm with a threshold of 0.6 and a minimum overlap of 20 nucleotides as implemented in PANDAseq (Masella et al., 2012). All the assembled reads were filtered for high-quality sequence reads (minimum sequence length, 350 and 120 nucleotides for bacteria and fungi, respectively; maximum sequence length, 500 and 580 nucleotides for bacteria and fungi, respectively; minimum average Phred score of 25; maximum length of 20 homopolymers in the sequence and without ambiguous nucleotides). Potential chimeras were removed using UCHIME (Edgar et al., 2011) as implemented in MOTHUR (Schloss et al., 2009). The high-quality reads were clustered into operational taxonomic units (OTUs) using cd-hit-est 4.6.2 (Fu et al., 2012) at a threshold of 97% pairwise similarity. The bacterial 16S rRNA OTU representative sequences were assigned against the SILVA v132 reference sequence database (Quast et al., 2013) to obtain the respective OTU tables. Fungal ITS representative sequences were assigned against the UNITE v7 sequence database (Kõljalg et al., 2013) using the Bayesian classifier as implemented in MOTHUR (Schloss et al., 2009). Singleton and doubleton OTUs originating from sequencing errors were removed from the datasets. The sequences that were classified as "Cyanobacteria," "Chloroplast," or "Mitochondria" and those that were not classified at the kingdom level were removed from the bacterial dataset. The ecological and metabolic functions of bacterial OTUs were predicted using FAPROTAX (Louca et al., 2016) and the functional annotation tool of prokaryotic taxa v.1.1, whereas those of fungal OTUs were predicted using FUNGuild (Nguyen et al., 2016). Additionally, the Tax4Fun (Aßhauer et al., 2015) R package, which employs 16S rRNA gene-based taxonomic information, and the Kyoto Encyclopaedia of Genes and Genomes (KEGG) database were used to predict the metabolic functional attributes of bacterial communities in the rhizosphere and endosphere of T. pratense.
Tax4Fun converted the SILVA-labeled OTUs into prokaryotic KEGG organisms and normalized these predictions using the 16S rRNA copy number (obtained from the National Center for Biotechnology Information genome annotations).

| Physicochemical analyses of the rhizosphere soil
The rhizosphere soil samples (100-200 g wet weight) from each plot were dried and sieved. The pH of the rhizosphere soil was measured using WTW Multi 3510 IDS. The rhizosphere soil was subjected to dry combustion at 1000°C to determine the total carbon (TC) and total nitrogen (TN) concentrations using a CHNS-Elemental Analyzer (Elementar Analysensysteme GmbH), following the manufacturer's instructions. Soil carbon/nitrogen (C/N) stoichiometry was calculated based on TC and TN. Available soil phosphorus was extracted and measured according to the Bray 1 method (Gutiérrez Boem et al., 2011). Cations (K + , Mg 2+ , Ca 2+ , and Na + ) in the rhizosphere soil were determined using an atomic absorption spectrophotometer (Hitachi Z 5300, Hitachi-Science & Technology), following the manufacturer's instructions. Physicochemical properties of soil did not differ significantly between current climate and future climate plots (Table A1).

| Statistical analysis
All statistical analyses were performed using the PAST program version 2.17c (Hammer et al., 2001) and R environment version 3.6.1 (R-Development-Core-Team, 2019). All the analyses were conducted based on five independent replicate plots of the field experiment (n = 5) for each treatment. The datasets were normalized to the minimum number of sequence reads per sample (5360 and 10,338 sequence reads for bacterial and fungal OTUs, respectively) using the function "rrarefy" from the vegan (Oksanen et al., 2019) package in the R environment version 3.6.1 (R-Development-Core-Team, 2019). To provide an overview of the bacterial and fungal operational taxonomic units (OTUs) distribution among different plant compartments, the shared and unique OTUs were represented using a Venn diagram with the software available at http://bioin forma tics.psb.ugent.be. The microbial diversity indices (Simpson's diversity, observed OTU richness, and estimated richness (Chao-1)) were calculated for both bacteria and fungi. Variance homogeneity was examined using Levene's test. The normal distribution of data was examined using the Shapiro-Wilk test. Since some samples' diversity was skewed, we used log 10 -transformed diversity indices data for further statistical analysis while the original values were used only for data visualization ( Figure 2). To test the influence of climate, plant compartment, and their interaction on microbial diversity, a splitplot analysis of variance (ANOVA) was performed using the function "sp.plot" from the agricolae R package (de Mendiburu, 2016). In detail, the impact of climate (two levels) was analyzed at the mainplot level, while that of the plant compartment (four levels) and both plant compartment and climate were analyzed at the sub-plot level.
Based on split-plot ANOVA results, the least significant difference (LSD) test was applied, using the function 'LSD.test', to show differences between treatments.
Microbial (bacteria and fungi) community composition was assessed by computing Jaccard and Bray-Curtis dissimilarity matrices and then visualized using non-metric dimensional scaling (NMDS) ordinations using the function "metaMDS" in the vegan R package (Oksanen et al., 2019) to visualize compositional differences.
To test whether ecological niche (plant compartment), climate, or their interaction had a significant effect on community composition, permutational multivariate analysis of variance (NPMANOVA) (Anderson, 2001), and analysis of similarities (ANOSIM) based on Bray-Curtis and Jaccard dissimilarities between microbial communities (OTU level) were performed for 999 permutations. Additionally, NPMANOVA pairwise post hoc comparisons were performed to evaluate the effect of the tested factors on bacterial and fungal communities separately using the function "pairwise.adonis" in the vegan R package (Oksanen et al., 2019) Only microbial endophytes were included, and no indicator species analysis was performed for the rhizosphere microbiome.

| Richness and diversity of T. pratense microbiome under current and future climate conditions
The distribution of bacterial and fungal OTUs in the plant compartments (rhizosphere and the root, leaf/stem, and flower endospheres) under both current and future climate conditions was analyzed ( Figure 1). The rhizosphere soil harbored the highest number of unique OTUs (43.5% and 41.5% for bacteria and fungi, respectively).
A large proportion of the OTUs in the rhizosphere was shared with the root (30.6% and 39% for bacteria and fungi, respectively), followed by leaf/stem (23% and 32.4% for bacteria and fungi, respectively), and flower (3.3% and 4.6% for bacteria and fungi, respectively) endospheres (Figure 1b,c). Only 1.7% of bacterial and 2.6% of fungal OTUs were shared among all compartments. The future climate condition-specific OTUs were the highest for the root (40% and 36%  of all bacterial and fungal OTUs, respectively) endosphere, followed by leaf/stem (39% and 38% of all bacterial and fungal OTUs, respectively) endosphere, rhizosphere (29% and 22% of all bacterial and fungal OTUs, respectively), and flower (22% and 25% of all bacterial and fungal OTUs, respectively) endosphere ( Figure 1d).
The effects of plant compartments, climate, or both on alpha diversity indices (Shannon's diversity, observed richness, and estimated richness) of T. pratense microbiomes were examined. Plant compartments had the highest influence on shaping microbial diversity and richness (Table 1). Climate change influenced the diversity and estimated richness of fungi (Table 1). Fungal diversity and richness were significantly higher in the leaf/stem and root endospheres under the future climate conditions compared to current climate conditions ( Figure 2).

| Community composition and taxonomic structure of T. pratense microbiome
The composition of the bacterial and fungal microbiomes of T. pratense at the OTU level (97% identity) was examined. NPMANOVA  (Tables A3 and A4). The analysis of Bray-Curtis distance revealed similar findings (Table A2).
In total, 21 fungal classes were detected. Of these, the abundance of six classes, which accounted for more than 57% of the Fusarium (relative abundances of 15.6% and 5% among all fungal sequences, respectively). The compartment dissimilarity based on genera was calculated using SIMPER analysis (Table A5). Allorhizobium,

| Analysis of plant compartment/niche and climate indicator species
Indicator species analysis identified the bacterial and fungal taxa that significantly benchmark each plant compartment/niche and/or climate. We detected 35 bacterial indicator OTUs (   Figure A6; Table A7), while 47 fungal genera represented by 177 OTUs were assigned as plant pathogens ( Figure A7; Table A8). Interestingly, climate conditions did not affect these microbial functions (Table A9).

| Prediction of the metabolic functions of the bacterial community using Tax4Fun
The potential metabolic functional profiles of bacterial microbiomes were predicted based on the 16S rRNA genes of retrieved bacterial taxa using Tax4Fun according to the KEGG Ortholog groups (KOs) .
The highly abundant metabolic genes (>0.001% sequence relative abundance) belonged to the following four categories: metabolism, genetic information processing, environmental information processing, and signaling, and cellular processes ( Figure A8). Climate condi-     Cregger et al., 2018;Fonseca-Garcia et al., 2016;Gargouri et al., 2021;Zheng & Gong, 2019). The nicherelated differences in the microbiome composition can be attributed to variations in the microbial pools that invade different plant tissues through vertical transmission from seeds or horizontal transmission from soil and atmosphere (Cregger et al., 2018). The variations in the density of invading microbes and the unequal distributions of nutrients and oxygen among different plant tissues can also be a reason for microbial variations among different compartments (Vandenkoornhuyse et al., 2015).
Additionally, consistent with the results of other studies, the microbial diversity and richness varied between the plant compartments in this study. The analysis revealed that the microbial richness decreased from the rhizosphere to the endosphere tissues. This is due to the secretion of root exudates containing organic and amino acids, sugars, vitamins, hormones, and growth regulating substances in the rhizosphere, which promote microbial growth and colonization (Berg et al., 2016;Turner et al., 2013). In contrast, limited nutrients and available intercellular space in the plant endosphere limit microbial growth and colonization. The horizontal transfer of fungal communities from the rhizosphere to the endosphere was higher than that of bacterial communities. Among the rhizosphere fungal OTUs, 39% were transmitted to the root endosphere, 35% were shared with root and leaf/stem endospheres, and 6% were shared with all compartments. Similarly, among the rhizosphere bacterial OTUs, only 29% were transmitted to the root endosphere, 18% were shared with the root and leaf/stem endosphere, and only 4% were shared with all compartments. This suggested that the host genetic regulation of the bacterial composition is higher than that of the fungal composition and that the levels of host-specific selection factors are high in the aboveground compartments. In contrast, the high level of specificity in the flower endosphere may indicate specific microbiome recruitment through air or pollinators (Vannette, 2020).
The analysis of the taxonomic composition of the clover microbiome in different compartments revealed that Actinobacteria and Sordariomycetes were the predominant microbes in the rhizosphere.
In the root endosphere, Alphaproteobacteria (nitrogen-fixing Rhizobia) was the predominant microbe, which was consistent with the results of a previous study on red clover (Hartman et al., 2017).
Additionally, the root endosphere was less frequently colonized by other potential N-fixing bacteria, such as Bradyrhizobium, Devosia, Ensifer, Burkholderia, Mesorhizobium, Microvirga, and Phyllobacterium (Table A7). Previous studies have reported that the roots of Trifolium repens and Trifolium fragiferum comprised Rhizobium as the predominant microbe with decreased abundance of rhizobia species, such as Bradyrhizobium, Sinorhizobium, and Mesorhizobium (Liu et al., 2007;Marilley & Aragno, 1999). The T. pratense root endosphere was enriched in OTUs of various genera, such as Actinoplanes and Pseudomonas. To the best of our knowledge, the microbial composition of the leaf/stem endosphere has not been previously investigated. In this study, the leaf/stem compartment predominantly comprised Gammaproteobacteria and Dothideomycetes. The species or strains of the most dominant bacterial and fungal genera in the root and leaf endosphere, such as Actinoplanes (Lazzarini et al., 2000), Pseudonocardia (Mangamuri et al., 2016), Streptomyces (Gouda et al., 2016), and Cladosporium (Gouda et al., 2016) are reported to synthesize medicinally important natural products. Flowers provide a unique habitat for microorganisms because of their ephemerality and anatomy, which form distinct micro-niches (Aleklett et al., 2014).
This study investigated the red clover inflorescence microhabitats (calyx, corolla, pistil, and stamen) as one unit. The microbiome of the flower predominantly comprised Gammaproteobacteria (Pantoea) and Mollicutes (Candidatus Phytoplasma), while the most prevalent fungal community members remained unidentified. A recent study on the seed-borne endophytes of T. pratense revealed that the predominant bacterial taxa were Gammaproteobacteria (63% of relative sequences abundance, with a dominance of Pantoea) and unidentified fungi (70% of relative sequence abundance). This indicated that these taxa could be unique members of T. pratense flowers that are transmitted to the next generation via seeds. Candidatus Phytoplasma, which is the obligate bacterial pathogen of plant phloem, is transmitted through plant propagation materials and seeds, as well as by insect vectors (Kumari et al., 2019). In this study, microbial genera unique to the flower were detected, including the two insect symbionts, Arsenophonus and Rickettsia, which are transmitted by various arthropods (Caspi-Fluger et al., 2012;Novakova et al., 2009). Thus, the plant served as a reservoir for the horizontal transmission of both bacterial genera.

| Red clover harbors various beneficial microbes for plant growth and system sustainability
The analysis of the predicted bacterial functional genes showed various genes involved directly or indirectly in plant growth initiation and adaptation to climate changes. For example, this study predicted the presence of bacterial genes involved in siderophore synthesis that indirectly induce plant systemic resistance by enabling bacteria to compete with pathogens through the removal of iron from the environment (Bakker et al., 2007). Moreover, genes involved in the production of phytohormones, such as auxin, indole-3-acetic acid (IAA), and 1-aminocyclopropane-1-carboxylic acid (ACC) that directly pro- In our study, Pseudomonas, Streptomyces, and Pantoea are three of the most abundantly detected genera in the rhizosphere and endosphere samples that are reported to promote plant growth and produce these bioactive compounds (Abbasi et al., 2019;Bakker et al., 2007;Jaemsaeng et al., 2018;Shariati et al., 2017). Furthermore, N-fixing and phosphate solubilization genes, which are involved in enhancing plant growth and nutrient release to the soil and reduce the need for N and P fertilization (Hayat et al., 2010), were predicted. Therefore, T. pratense is considered one of the most important soil biofertilizer forage crops that contribute to system sustainability.

| The impact of climate change on microbial community composition of T. pratense
Climate changes in terms of increasing temperature, summer drought, and altered precipitation patterns play a key role in shaping soil microbial communities (Mekala & Polepongu, 2019). However, few studies have investigated the effect of climatic conditions on plant-associated microbiomes. Drought conditions obstruct root development leading to the limitation of water and nutrients uptakes by plants and the diminishment of plant biomass (Al-Arjani et al., 2020;Hameed et al., 2014). In addition, severe drought may lead to over-accumulation of reactive oxygen species that result in extensive plant cell damage and death (Cruz de Carvalho, 2008).
Several plant-associated microbes were found to contribute to drought stress tolerance in plants by carrying out various strategies. For instance, arbuscular mycorrhizal fungi-plant associations lead to the induction of particular genes to elevated levels of expression such as P5CS involving in proline biosynthesis and genes coding for late embryogenesis abundant (LEA) proteins associated with ions and antioxidative stress system. Also, it regulates the abscisic acid (ABA) of plant content (Ahanger et al., 2014). Moreover, arbuscular mycorrhizal fungi could mitigate the negative effect of future climate conditions by altering the community composition and enhancing the richness of specific taxa (Wahdan et al., 2021).
Recent studies (Gargouri et al., 2021;Karray et al., 2020) performed on the genus Opuntia revealed that bacterial and fungal plant microbiomes changed in the rhizosphere and root endosphere along a climatic aridity gradient. Moreover, they identified specific biomarker taxa for each bioclimatic zone. Additionally, increasing the aridity resulted in highly cohesive soil microbial-root fungal networks. These microbial dynamics, biomarkers, and the highly correlative microbial networks could play a crucial role in the aridity stress and potentially promote the survival of Opuntia, one of the most xerophyte plants, across a wide range of arid zones (Gargouri et al., 2021). On the other hand, we have noticed that T. pratense resistance to cascading drought and rising soil temperature was limited. A marked reduction of T. pratense cover in the GCEF was detected after 4 years of growth under future climate conditions (unpublished results; Figure   A9). In contrast to previous studies, our results revealed that T. pratense harbored a highly conserved microbiome that did not provide plasticity to the host to acquire desirable microbes or reconstruct the community structure as observed in the bacterial community.
Fungal composition appeared to be more sensitive to environmental factors than bacterial composition, which was consistent with the results of previous studies Cregger et al., 2018;Fonseca-Garcia et al., 2016;Hacquard, 2016;Hamonts et al., 2018). We detected several dark septate endophytic fungal genera as indicators of future climate. Some of these indicators could be beneficial to plant growth and disease resistance, such as Cadophora, Myrothecium, and Stachybotrys (Banerjee et al., 2010;Busby et al., 2016;Yakti et al., 2019). However, these climate indicators are represented with few OTUs among the whole community.

| CON CLUS IONS
To the best of our knowledge, this is the first study to examine the composition and functions of bacteria and fungi in the four compartments of the forage legume crop T. pratense under both current and future climate conditions. Although the T. pratense microbiomes did not differ at the community level, it is possible that the microbial communities changed at the genomic level that was not detected by our approach (16S and ITS sequencing data). Therefore, further studies on microbial functions involving the integration of the highresolution metagenome, metatranscriptome, and metaproteome ap-

CO N FLI C T O F I NTE R E S T
None declared.

E TH I C S S TATEM ENT
None required.              Tax4Fun  TABLE A11 List of the selected KEGG enzyme-encoding gene for plant growth-promoting traits involved in biofertilization (nitrogen fixation, phosphate solubilization, and siderophore synthesis) and biostimulation (indole acetic acid (IAA) production, 1-aminocyclopropane-1-carboxylate (ACC) deaminase activity, and general plant growth-promoting traits). All data were extracted from the Kyoto Encyclopaedia for Genes and Genomes (KEGG) database www.genome.jp/kegg/   Figure A3 Effects of climate manipulation on (a) total precipitation (sum of seasons), (b) air temperature in a height of 10 cm (daily mean temperature + standard error), and (c) soil temperature (daily mean temperature + standard error) in a depth of 1 cm in experimental plots managed as extensively used grassland in the GCEF. Precipitation is not manipulated during the winter months. Note that the effects of soil temperature are strongly modulated by indirect effects via the change of vegetation cover (see also Schädler et al., 2019). Thus, the values presented here, are the net result of the direct increase of temperature manipulation and the indirect modulation by the vegetation cover.