Sasanishiki (lowland japonica rice) with a functional allele of qSOR1 and a near-isogenic line with a non-functional allele of qSOR1 (qsor1-NIL; Kitomi et al., 2020) with a Sasanishiki background was used. Hulled seeds were washed three times with sterile water. The seeds were then soaked in 1.0% (v/v) PLANT PRESERVATIVE MIXTURE (PPM™; Plant Cell Technology, Inc., USA) and incubated at 30°C to germinate for 24 hours. The germinated seeds were sown in a 0.4% (w/v) agarose gel (Sigma-Aldrich, USA) in a microplate-type Petri dish (Stem, Japan) ( Figure 1A ; Seeding). Roots protruding from the medium were susceptible to desiccation, rather than gravitropic stimulation, so the Petri dish was tilted 60 degrees and placed in the dark at 28°C for 14 h to allow the roots to elongate into the agar medium ( Figure 1A ; Preculture). They were then rotated with a 3D clinostat PMS-VI (AES, Japan) in the dark for 6.0 h at 28°C ( Figures 1A ; S1 ; Movie S1 ; Clinorotation). SMC was mimicked by adjusting the rotational speed of the 3D clinostat (X-axis 11.0 RPM/Y-axis 13.0 RPM). However, the influence of other physical stimuli that occurred by this device, such as vibrations, cannot be excluded as a possibility in this clinostat treatment. After simulated microgravity treatment (SMT), the Petri dish was placed vertically and gravistimulation was applied at 90 degrees ( Figure 1A ; Gravistimulation). Eighteen seeds were sown on each plate, and 10-16 seminal root tips with 2-3 mm were collected and pooled for RNA-seq. Seminal roots that were too short were excluded from the sampling ( Figure S2 ). The root tips were collected at 0.0, 0.5, 1.5, and 3.0 h after fGS. For the control plot without SMT, the Petri dish was tilted 60 degrees for 14 h, then placed vertically, and the root tips were sampled after 6.0 h in the dark ( Figure 1B ). Therefore, the control samples are the same age as the clinorotated sample at 0h after fGS.

Root tips were immediately frozen in liquid nitrogen and ground to a fine powder with a ShakeMaster BMS-A20TP (BMS, Japan). RNA was extracted with the RNeasy Plant Mini Kit (QIAGEN, Germany) according to the manufacturer’s instructions.

RNA-seq libraries were prepared using the NEBNext Ultra II Directional RNA Library Prep Kit for Illumina (E7760, New England Biolabs, USA) according to the manufacturer’s instructions. Sequencing of the libraries was performed with 150 bp paired-end reads on the S4 flow cells of the Illumina NovaSeq6000 platform at Macrogen Japan. The reads were mapped to the IRGSP-1.0 genome assembly with MSU7 gene model annotation using the STAR aligner program (ver. 2.7.3a) with options “–outFilterMultimapNmax 1 –quantMode GeneCounts” (Dobin et al., 2013). Differentially expressed genes (DEGs; |log2[fold-change]| >1, the false discovery rate <0.05) were called using the glmLRT of R package edgeR (ver. 3.26.8; Robinson et al., 2010) and transcripts per kilobase million (TPM) values were computed. Genes with TPM values >2 were considered expressed. Pearson’s correlation coefficients between samples were calculated using expression levels (log2[FPKM + 1]) of all expressed genes. Relative expression in Figures 2A ; S3A is shown as log2[FC to average FPKM of each gene in all samples]. Principal component analysis (PCA) was performed using TPM values of all expressed genes. All heatmaps were plotted using ComplexHeatmap (ver. 2.0.0; Gu et al., 2016). Gene ontology (GO) enrichment analysis was performed using clusterProfiler (ver 3.12.0; Yu et al., 2012) and visualized using corrplot (ver 0.84; Wei and Simko, 2017). An initial quality check using PCA and heatmap detected the wired behavior in one replicate of Sasanishiki, Sasa3 ( Figures S3A–C ). Therefore, Sasa3 was excluded from subsequent analyses, although we could not determine this reason.

Clustal Omega (Sievers et al., 2011) was used to calculate the homology between HSF proteins.

Pairwise Pearson’s correlation coefficients between the genes were calculated and gene pairs with correlation coefficients >0.9 were extracted for further analysis. Co-expression modules were analyzed using the “fastgreedy.community” function of igraph (ver. 1.3.4; Csardi and Nepusz, 2006). Modules with less than five genes were discarded. Cis-motifs enriched in the upstream regions of genes within each co-expression module were searched using MEME (ver. 5.4.1; Bailey et al., 2015) using the option “-mod zoops -nmotifs 3 -minw 6 -maxw 13 -revcomp -markov_order 0” and the detected motifs were compared with the Arabidopsis DAP-seq motif dataset (Bailey et al., 2015) using TOMTOM (ver. 5.4.1; Kitomi et al., 2020) and the options “-no-ssc -verbosity 1 -min-overlap 5 -dist pearson -thresh 1e-4.”

An overview of the SMT using a 3D clinostat and fGS treatment is presented in Figure 1A . Images of the roots at control, 0.0 h, and 3.0 h after fGS are shown in Figure 1B . The qsor1-NIL showed weaker gravitropism in the seminal roots compared with Sasanishiki (Kitomi et al., 2020). To identify genes downstream of qSOR1, we compared Sasanishiki and qsor1-NIL. Root elongation tended to be suppressed under SMC ( Figure S4A ). In the control without SMT, the root tips of both Sasanishiki and qsor1-NIL elongated toward the direction of gravity ( Figures 1B ; S4B ). In contrast, the roots of Sasanishiki and qsor1-NIL at 0.0 h after fGS showed various directions of elongation ( Figures 1B ; S4 ). After 3.0 h of fGS, the root tips of both Sasanishiki and qsor1-NIL elongated in the direction of gravity ( Figures 1B ; S4B ).

To identify genes associated with the gravitropic response in rice roots, we performed a time-series RNA-seq analysis along with fGS after SMC. We considered 0.0 h after fGS as under SMC. A heatmap of 16,648 expressed genes showed that transcriptional changes mainly occurred at 0.0 h and 0.5 h after fGS in qsor1-NIL and Sasanishiki ( Figure 2A ). The correlation matrix of expression levels revealed that replicates for each condition exhibited similar patterns (r > 0.99; Mann-Whitney U test P = 3.8e-9) and the time-series gene expression profile for qsor1-NIL broadly resembled that of Sasanishiki (r > 0.99; Mann-Whitney U test P = 2.0e-10; Figure 2B ). The controls without SMT were clustered with later-stage samples (1.5 h and 3.0 h) rather than early-stage samples (0.0 h and 0.5 h). PCA based on TPM values of the top 5% genes with the highest variability revealed that the difference in time points, rather than genotype, had more of an effect on gene expression ( Figure 2C ). This indicates that SMT-induced transcriptional changes were similar in both genotypes.

The expression of qSOR1 tended to be slightly lower in qsor1-NIL compared with that in Sasanishiki. Although the expression levels of qSOR1 were constant throughout the time-series in qsor1-NIL, qSOR1 expression in Sasanishiki was down-regulated at 1.5 h after fGS ( Figure S5 ). In contrast, no difference in DRO1 expression between Sasanishiki and qsor1-NIL was observed ( Figure S5 ). This suggests that the expression of qSOR1 and DRO1 in the root tips under SMC and subsequent fGS are similar between functional and non-functional qSOR1 genetic backgrounds. We identified only a few DEGs between Sasanishiki and qsor1-NIL after fGS ( Figure S1D ; Tables S1 , S2 ). The maximum number of DEGs that were up-regulated in Sasanishiki compared with qsor1-NIL was seven at 0.0 h, whereas the maximum number of DEGs that were down-regulated was 11 at 0.5 h and 1.5 h. Among these, no rice genes known to be involved in gravitropism was found. We also examined auxin-related genes associated with gravitropism; however, the expression patterns of these genes were almost identical between Sasanishiki and qsor1-NI ( Figure S6 ). These results indicate that qSOR1 is not involved in the regulation of gene expression patterns in the root tip during SMC and subsequent fGS.

We identified 152 and 117 DEGs in which expression was up-regulated during fGS in Sasanishiki and qsor1-NIL, respectively ( Table S3 ). In addition, we identified 286 and 333 DEGs in which expression was down-regulated during fGS in Sasanishiki and qsor1-NIL, respectively ( Table S4 ). GO enrichment analysis revealed that the down-regulated genes that were common to both genotypes were associated with response to heat, response to temperature stimulus, response to abiotic stimulus, cellular response to heat, response to hydrogen peroxide, protein folding, chaperone-mediated protein folding, and response to unfolded protein ( Figure 3 ). We found that the three GOs (response to heat, response to temperature stimulus, and response to abiotic stimulus) were detected in qsor1-NIL, but not in Sasanishiki, at 0.5 h after fGS. Because the expression levels of the six genes in these GOs were almost the same between the two genotypes, we concluded that this discrimination was dependent upon less statistical power in Sasanishiki at 0.5 h after fGS because of the number of replicates ( Figure S7 ). Based on these results, we focused our analysis on the DEGs that were common to both genotypes.

Common DEGs between Sasanishiki and qsor1-NIL in representative GOs were investigated ( Figure 4 ). One of the most prominent GO-enriched DEGs was protein folding. In Arabidopsis, HEAT SHOCK PROTEINs (HSPs) are involved in gravitropism through the regulation of auxin transport. HSP genes were enriched in the down-regulated DEGs ( Figure 4B ), which included five HSP20s, five HSP40s, eight HSP70s, and six HSP90s. The expression of most HSPs was higher at 0.0 h after fGS, but lower at 1.5 h and 3.0 h after fGS compared with those of the control without SMT ( Figure 4B ). These results indicate that the expression of HSP genes was induced by SMT, but repressed by prolonged fGS, likely reflecting an adaptive gravitropic response.

Transcription factors (TFs) regulate the expression of downstream genes and often play an important role in morphogenesis. Nine TFs were commonly up-regulated and 19 TFs were down-regulated during fGS in Sasanishiki and qsor1-NIL ( Figure 5 ). Of the 19 TFs, seven were HSFs, including HEAT SHOCK TRANSCRIPTION FACTOR A2D (HSFA2D), which acts upstream of LAZY1 and is involved in shoot gravitropism in rice (Zhang et al., 2018) ( Figure 5B ). The expression of these HSFs was higher at 0.0 h after fGS compared with those in the control without SMT, indicating that the expression of these HSFs was induced by SMT ( Figure 6 ). These HSFs belong to the HSFA2 and HSFB2 sub-groups. HSFs in other sub-groups were not induced by SMT, suggesting that HSFs in sub-groups HSFA2 and HSFB2 are specifically involved in the forced gravitropic response in rice ( Figure S8 ; Table S5 ).

The 3D clinostat is designed to minimize the effect of gravity by randomly and continuously changing the direction of gravity; however, the extent by which SMT reproduces real microgravity in space is not guaranteed. Transcriptome analysis of plants under microgravity in space have been conducted (Zupanska et al., 2013; Vandenbrink et al., 2019; Kruse et al., 2020; Barker et al., 2023). We compared genes detected as DEGs between our SMT experiment in rice and the previous space experiment in Arabidopsis (Zupanska et al., 2013). Many HSFs and HSPs were commonly detected DEGs in the two studies ( Table S6 ) as well as other studies (Correll et al., 2013; Choi et al., 2019; Barker et al., 2023; Paul et al., 2005; Paul et al., 2012; Kwon et al., 2015; Li et al., 2017). These results support that our SMT treatment properly mimics real microgravity in space.

To infer the gene regulatory networks under SMCs and following fGS, we performed a co-expression analysis using DEGs commonly expressed in Sasanishiki and qsor1-NIL. The DEGs were divided into five co-expression modules ( Figures 6A, B ; Table S7 ). The seven HSFs that were altered under the SMC were included in module 2 (M2), which was enriched with genes associated with “response to heat.” HSP40 genes (LOC_Os01g13760, LOC_Os03g57340, LOC_Os05g48810, LOC_Os06g02620), HSP70 genes (LOC_Os01g62290, LOC_Os03g16860, LOC_Os03g16920, LOC_Os05g38530), and HSP90 genes (LOC_Os04g01740, LOC_Os08g39140) were also included in M2, suggesting a possible interaction between HSFs and HSPs. We identified the enriched motif within a 200 bp upstream region of the genes belonging to M2 (55 out of 92 genes, 59.8%), which closely resembles the AtHSFA7 binding motif in the Arabidopsis cistrome dataset (q-value = 4.29e-9; Figure 6C ) (O’Malley et al., 2016), whereas TCP motifs (M1, q = 8.21e-6) and MYB-related motifs (M4, q = 2.45e-4) were identified in the other modules ( Figure S9 , q < 1e-4). Because Arabidopsis HSF-binding motifs are nearly identical to one another irrespective of class or sub-group, we anticipate that rice HSFs also recognize similar sequences ( Figure S10 ). These results suggest that HSFA2s and HSFB2 directly regulate the expression of genes in the M2 module and shape the M2 module.