Rhythmic Control of The Plant-Microbiome Holobiont

Background: Many physiological and metabolic processes in plants are modulated by a circadian clock. Rhizospheric microorganisms fed by root exudates provide benecial functions to their plant host. The intricate connection between the circadian clock and the rhizospheric microbial community remains poorly understood. Results: We investigated the role of the Arabidopsis circadian clock in shaping the rhizospheric microbial community using wild-type plants and clock mutants (cca1-1 and toc1-101). We performed transcriptomic and metabolomic analyses and sequenced amplicons of the 16S ribosomal RNA gene to characterize gene transcription, root exudation and the bacterial communities, respectively, throughout the day (24 h). Deciencies of the central circadian clock led to abnormal diurnal rhythms for thousands of expressed genes and dozens of root exudates. Bacterial community failed to follow obvious patterns in the 24-h period, and lack of coordination with plant growth in the clock mutants. Conclusions: Our ndings suggested that the biological clock was an important force that drove plants to adjust their rhizospheric microbiomes for adapting to different growth stages.


Introduction
The endogenous circadian clock in plants regulates most of the transcriptome and provides temporal coordination to ensure the optimal e ciency of many biological functions, including photosynthesis, carbohydrate metabolism and reproduction [1,2]. Plant roots host communities of rhizospheric microbes, which play a role in plant growth by regulating soil fertility, plant health and productivity [3][4][5]. The composition of the rhizospheric microbial community is predominantly determined by soil properties and environmental factors (e.g. climate change) but is also strongly affected by the genotype of the plant host [6,7].
The current concept of holobionts in intestinal and plant microbiotas considers the multicellular host and its associated microbiota as a functional entity [8]. The circadian clock of the host regulates the circadian rhythmicity and composition of the intestinal microbiota, and the diurnal rhythmicity of the microbiota reprograms the oscillations of the host transcriptome [9,10]. Both plants and microorganisms have an endogenous timer [11,12], so we can safely infer that the diurnal rhythmicity of the plants and microbes within the holobiont would in uence each other, regulating plant growth and the composition of the microbial community.
The structure of the bacterial community in the rhizosphere of Arabidopsis differs substantially between samples collected during the light and dark cycle, and this diurnal variation has been associated with the metabolism and exchange of carbon in the plants and microbes [13,14]. The molecular mechanisms involving a strong relationship between changes to the microbial community and the plant circadian clock nevertheless remain poorly understood. Understanding these mechanisms has thus become necessary.
The circadian clock in Arabidopsis involves a posttranscriptional component and a well-studied transcriptional-translational system of feedback regulation, in which circadian clock associated 1 (CCA1) and timing of CAB2 expression 1 (TOC1) are two core genes [12,15]. Mutants of these two genes disrupt the circadian rhythm [16,17]. We cultivated cca1-1 and toc1-101 mutants with the wild type (Col-0) under the same conditions to gain a better understanding of the role that the circadian rhythm of plants plays in the control of the rhizospheric microbial community. We performed transcriptomic and metabolomic analyses and sequenced amplicons of the 16S rRNA gene to comprehensively investigate whether the circadian rhythm affected the composition of the microbial community and, if so, whether this effect was intimately linked with root exudates.

Materials And Methods
Experimental design and seedling culture We used wild-type (Col-0) and cca1-1 and toc1-101 mutants (both de cient in central transcriptionaltranslational feedback loops of the circadian clock) of Arabidopsisthaliana to determine the role of plant circadian clock in determining the nature of the rhizospheric microbiome. The surfaces of Arabidopsis seeds were sterilized to prevent bacterial contamination of solid media and were vernalized after that.
The seeds were then cultured in Petri dishes containing Murashige and Skoog (MS) medium under sterile conditions at 25 °C at a light intensity of 300 μmol photons/m 2 /s and a 12:12-h light:dark photoperiod.
The MS medium contained 3% sugar and 0.5% agar and was autoclaved at 115 °C for 30 min before use.
Two-week-old aseptic seedlings were transplanted into polycarbonate pots (400 mL) containing autoclaved potting soil (Sun Gro Horticulture, Avagam, Massachusetts, USA) and soil slurry. Approximately 30 g of grassland soil collected near the Zhejiang University of Technology, China (30°17′45.11″N, 120°09′50.07″E) was mixed with 200 mL of sterile water by vigorous shaking for 60 s.
Eighteen milliliters of soil slurry were added to each pot before transplantation. Plants were grown in an arti cial greenhouse at 25 ± 0.5 °C and 80% relative humidity under cool-white uorescent light (300 μmol photons/m 2 /s) with a 12:12-h light:dark cycle. Samples for transcriptomic and metabolomic analyses and 16S rRNA gene sequencing were collected 14 days after transplantation. The timepoints for the analysis of circadian rhythmicity were 07:00 (one hour before the start of the light exposure, Zeitgeber time (ZT) 0), 13:00 (ZT6), 19:00 (one hour before the dark cycle, ZT12), 01:00 (next day, ZT18) and 07:00 (next day, ZT24). The details of the sample collection are provided below.
Sampling of plant shoots and RNA isolation About 0.1 g of aboveground plant tissue was weighed at each timepoint and placed into a 1.5-mL centrifuge tube, frozen in liquid nitrogen and stored at -80 ℃. Three replicate samples were collected for each genotype for transcriptomic analysis.
Library preparation and RNAsequencing RNA puri cation, reverse transcription, library construction and sequencing were all performed at Shanghai Majorbio Bio-pharm Biotechnology Co., Ltd. (Shanghai, China) following the manufacturer's instructions (Illumina, San Diego, USA). One microgram of total RNA for each genotype was used to prepare the RNA-seq transcriptomic library using the TruSeqTM RNA kit (Illumina, San Diego, USA) for sample preparation. Messenger RNA was isolated by using the polyA-selection method and oligo(dT) beads and then fragmented in a fragmentation buffer. A SuperScript double-stranded cDNA synthesis kit (Invitrogen, Carlsbad, California, USA) with random hexamer primers (Illumina) was used to synthesize double-stranded cDNA. The synthesized cDNA was end-repaired and phosphorylated, and adenines were added to the 3′ end following Illumina's protocol for constructing libraries. Libraries were size-selected for cDNA target fragments of 200-300 bp on 2% Low Range Ultra Agarose followed by PCR ampli cation using Phusion DNA polymerase (New England Biolabs, Boston, Massachusetts, USA) for 15 cycles with the default parameters mentioned in the instruction of Phusion DNA polymerase. A paired-end RNA-seq library was sequenced using the Illumina Novaseq 6000 system (2×150 bp read length) after quanti cation using a TBS380 mini-uorometer.

Analysis of the RNA-seq data
The sequences were trimmed and quality controlled using SeqPrep (https://github.com/jstjohn/SeqPrep) and Sickle (https://github.com/najoshi/sickle) using default parameters, and clean reads were separately aligned to the A. thaliana reference genome using TopHat (http://tophat.cbcb.umd.edu/, version2.1.1) in orientation mode. Only reads uniquely matched to the genome, allowing ≤2 mismatches without insertions or deletions, were mapped. The numbers of fragments per kilobase of exon per million mapped reads (FPKM) were calculated to represent the expression level. The ARSER algorithm [18] was used to detect circadian rhythmic genes using the parameters: minimal period, 20; maximal period, 28; default period, 24 and p = 0.05. A functional-enrichment analysis of the selected gene sets was performed to detect signi cantly enriched genes in the Kyoto Encyclopedia of Genes and Genomes (KEGG) relative to the whole-transcriptome background. The KEGG pathway analysis was carried out using KOBAS 2.1.1 (http://kobas.cbi.pku.edu.cn/download.php).

Preparation of exudate samples
Six replicates were taken from each genotype, and six seedlings were taken from each parallel. The plant was removed by cutting the polycarbonate pots which Arabidopsis planted in and its rhizosphere was gently nondestructively removed by shaking. The remaining soil attached to the roots was removed by washing with ddH 2 O, and the roots were then placed in a prepared petri dish (with an ori ce plate) containing 40 mL of sterile water. Each root was allowed to secrete for 6 h each replicate. The secretions were then collected in a 50-ml centrifugal tube using a 0.45-μm aqueous lter head, quickly frozen in liquid nitrogen and stored at -80 ℃. The weight of the plant and the volume of the secretions were recorded. The samples were thawed on ice, and the exudates were extracted in a 50% methanol buffer.
The exudates were centrifuged at 4000 g for 20 min, and the supernatants were transferred to new 96well plates. Samples for quality control (QC) were prepared by pooling equal volumes (10 μL) from all samples. The samples were stored at -80 °C prior to analysis by liquid chromatography and mass spectroscopy (LC-MS).
Measurement and identi cation of root exudates All chromatographic separations were performed using an ultraperformance liquid chromatograph (UPLC) (SCIEX, UK). An ACQUITY UPLC T3 column (100 mm×2.1 mm, 1.8 µm, Waters, UK) was used for the reversed-phase separation. The column oven was maintained at 35 °C, and the ow rate was 0.4 ml/min. The mobile phase consisted of solvents A (0.1% formic acid) and B (acetonitrile, 0.1% formic acid) with the following conditions for gradient elution: 0-0.5 min, 5% B; 0.5-7 min, 5 to 100% B; 7-8 min, 100% B; 8-8.1 min, 100 to 5% B and 8.1-10 min, 5% B. Exudates eluted from the column were detected using a TripleTOF5600plus high-resolution tandem mass spectrometer (SCIEX, Boston, USA) in both positive-and negative-ion modes. The Ionspray voltage oating was set at 5000 and -4500 V for the positive-and negative-ion modes, respectively. The MS data were acquired in IDA mode, and the TOF mass ranged from 60 to 1200 Da. The survey scans were acquired in 150 ms, and as many as 12 product ion scans were collected if exceeding a threshold of 100 counts per second and with a +1 charge. Total cycle time was xed to 0.56 s. Four-time bins were summed for each scan at a pulse frequency of 11 kHz by monitoring the 40-GHz multichannel TDC detector with four-anode/channel detection. Dynamic exclusion was set for 4 s. Mass accuracy was calibrated during the acquisition every 20 samples. A QC sample was collected after every 10 samples to evaluate the stability of the LC-MS throughout the acquisition. Raw LC-MS data les were converted to mzXML format and processed using the XCMS, CAMERA and metaX toolboxes in R. Each ion was identi ed by combining the retention times and m/z data. Exudates with mass differences between the observed and database values of <10 ppm were annotated using the online KEGG and HMDB databases. The molecular formulae of the exudates were determined and validated by measurements of isotopic distribution.

Analysis of metabolomic data
The intensity of the peaks of the exudates was calculated by the area of the peaks and preprocessed using metaX, and peaks were removed when detected in <50% of the QC samples or <80% of the biological samples. The remaining peaks with missing values were calculated with the k-nearest neighbor algorithm to further improve the quality of the data. The preprocessed data set was used for detecting outliers and evaluating batch effects using a principal component analysis (PCA). A robust QC-based LOESS signal correction was tted to the QC data in the order of injection to minimize the drift of signal intensity over time. Metabolic features with relative standard deviations >30% across all QC samples were removed. A principal coordinate analysis (PCoA) combined with a permutational analysis of variance (PERMANOVA) were performed to evaluate the dissimilarity between genotypes using Bray-Curtis distances. The ARSER algorithm was used to detect circadian rhythmic metabolites using the parameters: minimal period, 20; maximal period, 28; default period, 24 and p = 0.05. Spearman's rank correlations (r) between genus were calculated using the psych package in R. Only robust (r > 0.8 or r <−0.8) and statistically signi cant (P-value < 0.05) correlations were incorporated into co-occurrence network. Network visualization and modular analysis were made with Gephi (v 0.9.2). Node-level topological properties were also calculated using Gephi.

Collection of rhizospheres and bulk soils
Exposed topsoil (about 1 cm thick) was discarded, and the soil around the roots was gently removed to minimize mechanical damage to the roots. The roots were then placed in a 50-mL centrifuge tube lled with 20 ml of PBS buffer (137 mM NaCl, 10 mM phosphate buffer, 2.7 mM KCl, pH 7.3-7.5). The tube was then placed on a shaker at 180 rpm for 20 min, and the roots were discarded using sterile tweezers. The washing solution was centrifuged at 3900g for 20 min. The supernatant was discarded, retaining the rhizosphere [19]. The bulk soil (soil without roots) from the same pot as each sample of rhizosphere was also collected. Samples were stored at -80 °C prior to 16S rRNA gene sequencing. Four parallel samples were collected for each genotype for sequencing.
16S rRNA gene sequencing A FastDNA SPIN Kit for Soil (MP Biomedicals, California, USA) was used to extract total DNA from the microorganisms in the rhizosphere and bulk soil following the manufacturer's instructions. The V3-V4 region of the 16S rRNA gene was ampli ed using the primers 338F (5′-ACTCCTACGGGAGGCAGCAG-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) and a PCR thermocycler (GeneAmp 9700, ABI, USA).
Puri ed amplicons were sequenced on an Illumina MiSeq platform (Illumina, San Diego, California, USA). Sequences with ≥97% similarity were assigned to the same operational taxonomic unit (OTU). Taxa were assigned using the SILVA database (https://www.arb-silva.de/) based on the RDP classi er (Version 2.11, https://sourceforge.net/projects/rdp-classi er/) algorithm. The microbial compositions were also assessed using the UniFrac distance, which is a phylogenetically based measure of the degree of similarity between microbial communities [20]. Pairwise distances were determined for all pairs of samples based on either the taxonomic representation (unweighted) or the types and relative abundances (weighted). A PCoA of the UniFrac distances was used to visualize the differences among samples. A PERMANOVA was performed to evaluate the variance between genotypes. The ARSER algorithm was used to detect genera with circadian rhythms using the parameters: minimal period, 20; maximal period, 28; default period, 24 and p = 0.05.

Analysis of transcriptomes, root exudates and 16S rRNA gene sequences
The Procrustes test was used to analyze the correlations between rhythmic genes in the wild type and all exudates using the vegan package in R (3.6.2). Multiple linear regression was used to analyze the regulation of all exudates from the 10 most abundant phyla (Proteobacteria were divided into three classes) using SPSS (v 20.0.0).

Results
Physiological differences between the wild type Arabidopsis andclock mutants To explore differences in plant growth caused by an abnormal circadian clock, we measured the fresh weights and bolting rate of the wild type (Col-0) and the two mutants (cca1-1 and toc1-101). The fresh weights of 28-day-old seedlings of the wild type and the two mutants were similar and did not differ signi cantly ( Figure S1). The bolting rate after 15-28 d of culture was slightly lower (5-10%) for the toc1-101 mutant than the wild type, and the bolting rate was much higher (22-68%) for the cca1-1 mutant than the two other lines ( Figure S2). The phenotype of extremely early owering in the cca1 mutant was consistent with previous reports [21]. The bolting rate determined the timing of owering and was an indicator of the transition of Arabidopsis from the vegetative to the reproductive phase. These results indicated that the toc1-101 and cca1-1 mutants grew as well as the wild type, but their reproduction was intervened by the abnormal circadian rhythm.
The composition of the rhizospheric microbial community of plants is strongly regulated by plant development [6,22], implying that plants select for specialized microbial communities that change in response to plant growth, and this selection is believed to be conducive to plant growth or reproduction [4,6,23]. A previous study has reported an apparent discrepancy of microbial community structure in rhizosphere of Arabidopsis between the vegetative and bolting stage [6], indicating that the microbial community of the rhizosphere is continually changing, during this shift, so we explored whether the plant circadian clock in uences the rhizospheric microbial community during the transition from vegetative to reproductive phase.

Diurnal transcriptomic variation between wild type Arabidopsis andclock mutants
To further investigate the detailed cycling dynamics of transcription, we sampled leaves every 6 h over a 24-h period (from 07:00 to next 07:00), for a total of ve timepoints, hereafter designated as ZT 0-24, from the wild type and mutants for a transcriptomic analysis (Figure 1a). The circadian rhythmicity of the genes in both the wild type and mutants were identi ed using the ARSER algorithm [18,24] (Figure 1b). Heatmaps of the patterns of cycling genes (ARSER, p < 0.05) clearly differed between the two mutants and the wild type ( Figure S3a, c). The cca1-1 mutant lost 4902 cycling genes and gained 4940 cycling genes compared to the wild type, similar to the toc1-101 mutant. The detailed classi cation of the lost and gained cycling genes in the cca1-1 and toc1-101 mutants ( Figure S3b, d) indicated that the rhythmic operation of multiple pathways, including many types of metabolic and transport pathways, differed greatly from the wild type. The cycling genes shared between the wild type and the mutants were involved in photosynthesis, the TCA cycle, oxidative phosphorylation and nucleotide excision repair, indicating that the mutants retained the rhythmicity of the most basic activities for plant life. The growth (represented by fresh weight) and phenotypes of the mutants thus did not differ signi cantly from those of the wild type ( Figure S1).
The wild type and the two mutants shared 1800 cycling genes, but their oscillating loci were generally distinct (Figure 1c-e). The peak expression of most cycling genes in the wild type ranged between ZT15 and ZT18. The expression of cycling genes peaked equally between ZT12 and ZT24 in the cca1-1 mutant, but more genes peaked between ZT18 and ZT21 in the toc1-101 mutant, i.e. many genes shifted their phase in the mutants, although the oscillating amplitude (maximum diurnal variation in gene expression) of these genes was not signi cantly affected in the mutants. The detailed classi cation of genes with phase changes in the mutants indicated that their rhythmic patterns, including amplitude and phase, were strongly affected by the disruption of the central circadian clock (Figure 1f), although most basic activities for plant life retained their rhythmicity.
The circadian clock in Arabidopsis plays an important role in regulating the biosynthesis of central and secondary metabolites [25]. A Procrustes analysis indicated that the variation of the cycling in metabolic pathways, and changes to the phase shift and amplitude of oscillating genes, affected the synthesis and transport of material in the plant, which may further affect the root exudates ( Figure S4).

Distinct root exudates between wild type Arabidopsis and clock mutants
Plant-root exudates account for 10-40% of the photosynthetically xed carbon and ~15% of total plant nitrogen, including many kinds of primary and secondary metabolites [26], and act as key substrates or signaling molecules that affect microbial composition [6,[27][28][29]. We collected the root exudates of the three Arabidopsis lines (wild type and cca1-1 and toc1-101 mutants) at ve timepoints, with timepoint representing the continuous secretion of exudates for the previous 6 h, e.g. the exudates sampled at ZT12 consisted of the exudates secreted between ZT6 and ZT12 (Figure 2a). A total of 196 exudates were identi ed (Data set S1), 42% of which were lipids and lipid-like molecules. The daily production (total of all timepoints) of some exudates differed between the wild type and the mutants. Daily production differed from >2 to <0.5-fold for 9-11 exudates between the mutants and the wild type.
A global principal coordinate analysis (PCoA) combined with an Adonis analysis (Figure 2b) found that both genotype (R 2 =0.04, p < 0.01) and timepoint (R 2 =0.33, p < 0.001) could be linked to the pro le of exudates in Arabidopsis. The exudates in the wild type differed signi cantly between ZT0 and ZT6 ( Figure S5). A previous study reported that the expression of many photosynthetic genes peaked near ZT4 (the light cycle in our study started one hour after ZT0), when large amounts of sugar and other products of carbon xation are synthesized in Arabidopsis [30]. The exudates collected at ZT6 (continuously secreted by the roots from ZT0 to ZT6) would therefore bene t from the products of carbon xation. In sharp contrast to ZT6, exudates at ZT0 were secreted at the end of the dark cycle when the energy source in the plant was nearly depleted. These metabolic discrepancies may account for the dissimilarity between the exudates collected at ZT0 and ZT6. The differences between timepoints were not as obvious in the cca1-1 mutant as the differences in the wild type, although the Adonis analysis identi ed a signi cant difference between timepoints (R 2 =0.55, p < 0.001). Exudates for the toc1-101 mutant at ZT0, ZT6 and ZT12 differed from those secreted at other timepoints.
We detected 50 rhythmic circadian exudates in the wild type using the ARSER algorithm (Figures 2c, S6) but fewer in the cca1-1 (39) and toc1-101 (45) mutants. Some exudates retained rhythmicity in the mutants but had different rhythmic features (only four metabolites were concurrently rhythmic in all genotypes). The amplitude and phase for the cycling exudates shared between the mutants and the wild type varied with genotype (Figure 2d). These results were similar to those from the transcriptomic analysis, indicating that the plant could change the regularity of secretion for some exudates when mutations in clock genes induced large-scale changes in gene transcription and subsequently affected metabolism ( Figure S3b, d). The rhythm of exudation in both mutants and the wild type, however, was notably less important than the rhythm of gene expression. This difference may have been due to i) the joint control of metabolite synthesis by multiple genes, not all of which were rhythmically expressed, ii) a time lag between gene expression and metabolite synthesis and secretion or iii) the collection of exudates after six hours of accumulation, which may have reduced the importance of the rhythm.
The heatmap in Figure S7a represents the number of cycling exudates shared between the cca1-1 mutant and the wild type (10), with some exudates unique to the wild type (40) and the cca1-1 mutant (29). The heatmap in Figure S7b represents the number of cycling exudates shared between the toc1-101 mutant and the wild type (17), with some exudates unique to the wild type (33) and the toc1-101 mutant (28). Exudates that varied their rhythmicity in the mutants were mainly lipids and lipid-like molecules. The ve most rhythmic lipids and lipid-like molecules in the wild type were thus chosen based on their pvalues calculated using the ARSER algorithm for demonstrating the effects of the central circadian clock on metabolism. The levels of these ve exudates varied rhythmically in the wild type, but the exudates in the mutants lost their rhythmicity (dihydroterpineol, phosphatidylethanolamine (15:0/18:0) and piperochromenoic acid) or had abnormal amplitudes and phases (acylcarnitine 14:2 and citronellyl acetate). Both the abnormal amplitude and phase of exudation would lead to different trends in development of the rhizospheric microbiome and would lead to distinct community structures in the three Arabidopsis lines.

Temporal changes of microbial community in rhizosphere and bulk soil
We sampled the rhizosphere soil of the wild type and the cca1-1 and toc1-101 mutants every 6 h for 24 h to assess the in uence of the circadian behavior of plants on the rhizospheric microbiome (Figure 3a). An equal number of samples of rhizosphere and bulk soils were collected from the same pots, so each genotype at each timepoint had a bulk soil as a control. A PCoA based on unweighted UniFrac distances clearly indicated that the rhizosphere and bulk soils had distinct microbial communities (Figure 3b). The initial soil microbial communities in the pots containing the cultivated cca1-1 and toc1-101 mutants and the wild type were from the same soil collected from a grassland. The composition of the microbiota in the bulk soil, however, was still very similar among the mutants and the wild type after 14 d of cultivation ( Figure S8a), but the microbial communities in the rhizosphere soils differed signi cantly between the mutants and the wild type ( Figure S8b), indicating that the circadian behavior of the plants adjusted the rhizospheric microbiome during cultivation [7].
The composition of the microbial community varied signi cantly over time (unweighted UniFrac, p = 0.002) in the rhizosphere of the wild type but not the mutants ( Figure S9). The ARSER algorithm identi ed more rhythmic OTUs in the rhizosphere of the cca1-1 (19.80%) and toc1-101 (22.67%) mutants than the wild type (14.43%) ( Figure S10). Most of these rhythmic OTUs belonged to rare taxa, as described by Xue et al [31] (Figure 3c). The distinct patterns of dynamic change of the microbiomes during the 24-h period in the rhizosphere soil of the mutants and wild type differed between ZT0 and ZT6, ZT12, ZT18 and ZT24 based on unweighted UniFrac distances (Figure 3d). The rhizospheric microbiome of the wild type was clearly coordinated with plant growth during the 24-h period as the dissimilarities to ZT0 continuously increased. The rhizospheric microbiomes for both the cca1-1 and toc1-101 mutants were more similar between ZT0 and the other timepoints, indicating that the rhizospheric microbiomes showed randomness and had no obvious temporal pattern during the 24-h period.
Proteobacteria and Actinobacteria were the two most abundant phyla of the rhizospheric microbiomes for all genotypes (Figure 3e), consistent with the core root microbiome of Arabidopsis de ned by Lundberg et al [7]. The relative abundances of both Proteobacteria and Actinobacteria decreased for the wild type during the 24-h period. The relative abundance of Proteobacteria for the mutants uctuated slightly during the 24-h period, but the relative abundance of Actinobacteria remained stable. The cooccurrence network analyses indicated that the rhizospheric microbiome was more complex for the wild type (average of 16.542 nodes) than the cca1-1 (average of 9.321 nodes) and toc1-101 (average of 7.96 nodes) mutants (Figure 3f, Table S1). The stronger species interactions in the wild type likely decreased the stability of the microbial communities, which could account for the highly dynamic development of the rhizospheric microbiome during the 24-h period [32].

Discussion
The compositions of the rhizospheric microbial communities of the two mutants and the wild type uctuated with distinct patterns, indicating that the mutant genes of the core circadian clock in Arabidopsis led to different programming effects on the rhizospheric microbial communities. The abundance of the core microbes for the wild type was regulated by rhythmic root exudation and developed continuously due to the strong species interactions. This dynamic pattern may account for the effect of plant development on the composition of rhizospheric microbiome [6]. The composition of the rhizospheric microbiome for the mutants was also regulated by rhythmic root exudation, but it recovered to the status of ZT0 during the 24-h period due to the weak species interactions. The more rhythmic the dynamic pattern of the rhizospheric microbiome, the slower the coordination of the rhizospheric microbiome with plant development.
The rhythmicity of gene transcription differed signi cantly between the wild-type and mutant lines ( Figure  1). The mutation of CCA1 and TOC1 did not reduce the number of cycling genes, but most cycling genes differed between the two mutants and the wild type, i.e. the mutation of the genes controlling the core circadian clock created another gene oscillating model in the mutants. The loss of genetic rhythmicity or the abnormality of its amplitude or phase in the cca1-1 and toc1-101 mutants were closely associated with the metabolism of lipids, carbohydrates, energy and amino acids. These metabolic pathways determine the synthesis of root exudates. A previous report found that the secretion of many chemicals from roots, such as phytochemicals, was not predominantly in uenced by the circadian clock [33], but some exudates in our study had diurnal rhythms, especially lipids (Figures 2b, S7).
We collected our samples during the critical transition of Arabidopsis from the vegetative to the reproductive phase. The root exudates at this stage contributed mainly to shaping the rhizospheric microorganisms (Figures 4, S11). A multiple linear regression found that 75 root exudates signi cantly regulated the relative abundance of some phyla or classes of Proteobacteria ( Figure S11, Data set S2). Half of these exudates were lipids and lipid-like molecules, which are critical exudates transferred from the plant to the microbes [34,35] and play important roles in the communication between them. The lipid exudates differed between the mutants and the wild type ( Figure S7), potentially due to the abnormal rhythmicity of the genes involved in lipid metabolism ( Figure S3). For example, the transcriptomic analysis indicated that the rhythmicity of genes involved in sphingolipid metabolism (cca1-1) and glycerolipid metabolism (toc1-101) varied more than the rhythmicity of these genes in the wild type ( Figure S3b, d). These two lipids are involved in establishing the membrane interface between plant and microbial cells [35].
Staley et al [13] found that the composition of organic matter differed signi cantly between rhizospheric samples during a light and dark cycle and that 13% of the microbial community was cycling, suggesting that the rhizospheric bacterial community was regulated by the circadian clock. We also found that some bacteria in the rhizospheres of all genotypes had a diurnal rhythm ( Figure S12). More interestingly, however, the entire microbial community in the rhizosphere of the wild type varied continuously from ZT0 to ZT24 (Figure 3d, e) rather than rhythmically. In contrast, the cca1-1 and toc1-101 mutants, which had abnormal circadian clocks, were weak drivers for adjusting the microbiome toward the "designated" target. Similarly, the abundance of some rhizospheric microorganisms in barley is not circadian but uctuates from light to dark timepoints [36]. The rhythmicity of genes may have led to the rhythmicity of some of the root exudates, but this rhythmicity could not be extrapolated to the microbial community. This nding was similar to the key ndings of a force analysis: a rhythmic force acting on an object does not lead to the oscillation of the object, but accelerates the object. The rhythmic secretion of organic matter by plants would thus be bene cial for transforming the rhizospheric microbiome into a special community. This process is strictly programmed and controlled mainly by a biological rhythm. The disturbance of this rhythm or a change to a new oscillating model in the mutants would slow the dynamic uctuation and development of the rhizospheric microbiome. We therefore hypothesize that the circadian clock in plants is an important driving force for the positive development of rhizospheric communities and thus plant development (Figure 4), which has been neglected in most previous studies.
Hubbard et al [14] found that the mutation of two clock genes (TOC1 and ZTL) had a strong impact on the A. thaliana rhizospheric microbiome, particularly the rare taxa. They also found that wild-type plants germinated earlier and were larger when inoculated with soils where the wild type had grown compared to genotypes with mutant clock genes. This nding implied that the microbial community in the rhizosphere of plants with a normal circadian clock would be bene cial to the germination and growth of their offspring. The role of the circadian clock is particularly important during owering, and the mutation of core genes directly affects the timing of owering [37]. Rhizospheric microorganisms can in uence the timing of owering [4], so an abnormal clock may affect the rhizospheric microbes and ultimately indirectly affect owering, in addition to the direct in uence of the clock oscillator.

Conclusions
Previous studies recently reported oscillating variations of the relative abundances of taxa in plant rhizospheric microbiomes [13,14,37]. Our study has made a step forward by connecting the circadian clock of the plant to the rhizospheric microbial composition and nding a strong in uence. We do not know, though, the speci c molecules dominating the regulation of microbial community by biological rhythms nor how does an abnormal microbial community generated by disturbed biological clock in turn affect plants. Future research will closely revolve around these two questions. Our results show that the