Transcriptome profiling of laser-captured crown root primordia reveals new pathways activated during early stages of crown root formation in rice

Crown roots constitute the main part of the rice root system. Several key genes involved in crown root initiation and development have been identified by functional genomics approaches. Nevertheless, these approaches are impaired by functional redundancy and mutant lethality. To overcome these limitations, organ targeted transcriptome analysis can help to identify genes involved in crown root formation and early development. In this study, we generated an atlas of genes expressed in developing crown root primordia in comparison with adjacent stem cortical tissue at three different developmental stages before emergence, using laser capture microdissection. We identified 3975 genes differentially expressed in crown root primordia. About 30% of them were expressed at the three developmental stages, whereas 10.5%, 19.5% and 12.8% were specifically expressed at the early, intermediate and late stages, respectively. Sorting them by functional ontology highlighted an active transcriptional switch during the process of crown root primordia formation. Cross-analysis with other rice root development-related datasets revealed genes encoding transcription factors, chromatin remodeling factors, peptide growth factors, and cell wall remodeling enzymes that are likely to play a key role during crown root primordia formation. This atlas constitutes an open primary data resource for further studies on the regulation of crown root initiation and development.


Introduction
Crown roots (CR) are a type of adventitious roots making up most of the root system in rice [1,2]. These roots develop post-embryonically from the stem and play an important role in the adaptation of plants to the soil hydro-mineral status [3,4]. Genetics and functional genomics approaches revealed several genes involved in CR development in rice [5,6]. However, these approaches are impaired by functional redundancies, and by lethality of some mutants. a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 Transcriptomic approaches, based on the quantification of gene expression at wholegenome level, overcome these limitations by identifying all the genes specifically expressed during a biological process. This type of approach led to the identification of new genes involved in root development and root physiology in Arabidopsis [7][8][9]. Until now, most of the rice root transcriptome datasets are based on the comparison at the whole root systems level between stressed/unstressed material, developmental stages, genetic backgrounds (knock out mutant or overexpressing lines versus wild type) or different root types or root tissues. For instance, transcriptomic analysis of root response to soil potassium deficiency revealed the importance of an ion transporter and several protein kinases encoding genes in the adaptation of the plant to this constraint [10]. Comparison between CR, large and fine lateral roots revealed a significant enrichment in transcripts associated with secondary cell wall metabolism, as well as increased levels of transcripts encoding auxin-signaling related genes in CR [11]. Recently, exogenous hormonal treatments and transcriptome analysis of rice stem base allowed the identification of a new gene involved in CR development [12].
We previously used the conditional ectopic induction of the CROWN ROOTLESS 1 transcription factor in its respective mutant background to control the formation of CR in stem bases and generate a transcriptomic dataset capturing the kinetics of gene expression during very early stages of CR formation [13,14]. These approaches are powerful but lack spatial resolution. High spatial resolution transcriptomic analyses can be obtained using laser capture microdissection (LCM). For instance, Jiao et coll. produced a cell type transcriptome atlas that includes 40 cell types from rice shoot, root and germinating seed at several developmental stages [15]. One of the most up-to-date root-based spatialized dataset focused on dividing the CR into eight developmental stages along the longitudinal axis, and three radial tissue types namely: (i) epidermis, exodermis and sclerenchyma; (ii) cortex; (iii) endodermis, pericycle and stele [16]. However, while this dataset covers different tissues of emerged CR, it does not cover the early stages of CR development.
In this paper, with the aim to complete our knowledge about genes expressed during CR formation, we used laser capture microdissection (LCM) to generate an atlas of genes differentially expressed during three early stages of CR primordia development before CR emergence. Using this data resource, transcriptome cross analysis identified key genes already known to be involved in CR initiation and development, and pinpointed new candidate genes involved in gene expression regulation, signaling pathways and cell differentiation and likely controlling CR development.

Plant material and growth conditions
Oryza sativa L. ssp. Japonica cv. Taichung 65 seeds were disinfected (3 min in ethanol 70%, 90 min in sodium hypochlorite 3.8%:diluted in Milli-Q water plus 500 0L of Tween1 80 (Merck, Germany), quadruple rinsed in sterile Milli-Q water) and germinated under axenic conditions in round Petri dishes (Sarstedt, Germany) containing a filter pad (Whatman paper, GE Healthcare, UK) and 15 mL half strength Murashige and Skoog (MS/2) medium (Duchefa Biochemie, The Netherlands). Culture chambers were set at 27˚C, 70% relative humidity, and 14 hours daylight. After 5 days, the plantlets were transferred into 250 mL wide-collar Erlenmeyer flasks containing 30 mL MS/2 medium, and let to grow during 3 days before sampling.

Stem base sampling, two-photon imaging, tissue preparation, laser capture microdissection and RNA extraction
For two-photon imaging, stem bases were sampled and immersed in MS/2 medium at 4˚C overnight. They were included into 5% w:v agarose in sterile water (Duchefa Biochemie, The  [17]. For 3D reconstruction, aligned sections were used as an input in Imaris v9.1 (Bitplane AG, Germany). For LCM, stem bases were harvested in 2 mL ice-cold ethanol:acetic acid 3:1 (v:v) fixative solution, infiltrated thrice under a mild vacuum for 15 min and stored for 20 h in fresh fixative solution at 4˚C. Embedding for dissection was performed as described [18]. For post-embedding RNA quality assessment, stem bases were deparaffinized with three 10 min xylene baths at 4˚C, tissue disruption was done in liquid nitrogen using a TissueLyser II system (Qiagen, The Netherlands) with 3 mm steel beads for 15 sec at 30 Hz, and RNA isolation performed using the Plant RNEasy Kit (Qiagen, The Netherlands). RNA quality was assessed on a 2100 Bioanalyzer with the RNA 6000 Pico Kit (Agilent, Santa Clara, CA, USA), and RINs exceeded 9.3. Six embedded stem bases were cut into 15 μm sections on a Jung RM2055 (Leica Biosystems, Germany) or a Microm HM355S (Thermo Fisher Scientific, Waltham, MA, USA) microtome, and microdissected on a LMD 7000 laser dissector (Leica Biosystems). CR primordia belonging to the three successive rings present in stem bases at this stage were separately collected. Stem cortex tissues were sampled from the same sections where the CR primordia were collected from. Three pairs of two stem base were used as independant biological replicate. For each of the three rings observed in each pair of stem base, all crown root primordia (2~4 per ring) were collected, pooled, and used for RNA extraction in order to gather enough biological material. RNA extractions were performed using the Arcturus PicoPure RNA Isolation Kit (Thermo Fisher Scientific). One-round amplification and cDNA synthesis was done using the GeneChip WT Pico Kit (Applied Biosystems, Thermo Fisher Scientific) at the Transcriptomics platform of Montpellier University Hospital (CHRU Montpellier, France).

Array hybridization and analysis
All array hybridization steps were performed at the Transcriptomics platform of Montpellier University Hospital (CHRU Montpellier). Array hybridization on Gene Rice (Jp) 1.1 ST chips was done on the GeneAtlas system according to the manufacturer's instructions (Affymetrix, CA, USA). Probe intensities were normalized using Robust Multi-array Average [19] as implemented in Transcriptome Analysis Center software version 4.0.0.25 (Applied Biosystems, Thermo Fisher Scientific). Our dataset was deposited in NCBI's Gene Expression Omnibus [20] and is accessible through GEO Series accession number GSE133593.

Determination of differentially expressed genes, identification of coexpression clusters, and exploration of ontologies
Differentially expressed genes (DEG) were detected among the 29,664 genes present on the chip using the built-in Transcriptome Analysis Console (TAC) version 4.0.0.25 (Applied Biosystems, Thermo Fisher Scientific) function, using a fold-change cutoff of 1.5 and a p-value cutoff of 0.01. Co-expression clusters were identified by performing weighted gene coexpression network analysis (WGCNA) using the consensus method described by [21] on the normalized expression levels of the 3975 differentially expressed genes provided by TAC and for the four tissue types (cortex plus three stages for CR primordia, three replicates) considered as treatment factors. For our list of DEGs, we used the antilog normalized expression values, and chose a soft-thresholding power of P = 6. Coexpression cluster naming used the WGCNA convention of using color names. MapMan ontologies were assigned to detected DEG according to the ontology mapping BinTree RAPDB-IRGSP1.0 version 1.0 [22]. Granularity of the Map-Man bins, that are nested ontologies, was reduced to the second level (e.g. 27.2: RNA.transcription) for readability. The most represented ontologies were identified by considering the ontologies assigned to at least 2% of the total number of DEG, excluding the ontology 35.2 corresponding to "not assigned.unknown and unmapped genes". From the standard, most granular assignment of ontologies to DEGs, a parametric statistical enrichment in ontologies was computed using PageMan as implemented in MapMan v3.6.0-RC1 by performing a bin-wise Wilcoxon test, with Benjamini and Hochberg adjusted p-values [23]. Among the DEGs, the genes encoding transcription factors were identified using annotations from PlantTFDB v4.0 [24].

Three rings of CR primordia corresponding to 3 different developmental stages are present in 8 days old rice seedlings
The stem base of rice plants eight days after sowing was imaged to characterize the CR primordia developmental stages that were present. We observed three rings of CR primordia in the stem base (Fig 1A). Each ring corresponds to sequential events of CR initiation and therefore represents a different developmental stage [1]. CR primordia belonging to each ring are representative of a given CR developmental stage and were named as previously defined [13], from the top to the base of the stem, early, intermediate and late.
A progressive shift in transcriptional regulation occurring during crown root formation activates peptide mediated signaling cascades elements and repress specific cell wall metabolism genes Early, intermediate and late CR primordia, and cortex tissue neighboring the three stages of CR primordia were sampled using LCM from a total of six stem base (Fig 1B-1H). Transcriptome analyses were conducted on these different samples. A total of 3975 genes were detected as differentially expressed between cortex and CR primordia of all stages (S1 Table, Fig 2A).
2208 and 1767 genes were up-or downregulated, respectively. Nearly 30% of these differentially expressed genes (DEGs) were shared between the three developmental stages, and 10.5%, 19.5% and 12.8% of the DEGs were specific of the early, intermediate and late stages, respectively. 10% were shared between early and intermediate, and 15.3% between intermediate and late stages. Only 2% were shared between late and early stages, suggesting a progressive shift in genome expression profile during CR primordia formation. Furthermore, 10 co-expression clusters were identified, characterized by homogeneous changes in expression profiles across samples, and which can be used to explore genes implied in tightly regulated biological processes (S1 Fig).
The most strongly DEGs (detected in 1st and 99th centiles of the observed FC for at least one of the CR primordia developmental stage) are listed in S2 Table. Most of the genes are assigned to transcription or protein synthesis or of unknown function and their role in crown root formation should be further studied. It is interesting to note that XYLOGLUCAN ENDO-TRANSGLUCOSYLASE/HYDROLASES (XTH) OsXTH17 gene was among the most downregulated genes. OsXTH1, OsXTH6 and OsXTH10 were also significantly downregulated in CR primordia compared to cortex (S1 Table).
We next analyzed the ontologies associated with the different stages of development by assigning MapMan ontologies to the DEGs (Figs 2B and 3) and statistically tested their enrichment (S2 Fig). Altogether, this showed that the upregulated DEG were more frequently distributed in the protein synthesis, RNA/regulation of transcription, protein degradation, RNA processing, and DNA synthesis/chromatin structure categories, which can be expected from meristematic tissues. Notably, the protein synthesis category (with 89% of mapped genes corresponding to ribosomal proteins) was strongly represented and significantly enriched, indicating a dynamic protein synthesis and regulation in all three developmental stages. Also, we observed a high number of DEGs in the regulation of transcription category (377 up-and 257 downregulated).
Consistently, DEGs assigned to the "33.99 development.unspecified" MapMan ontology are also enriched in transcription factors (TF; S3 Table). Interestingly, the two genes that are the most upregulated at the early stage and that decrease in subsequent stages are PLETHORA 3/ OsPLT3 and PHYTOSULFOKINE 1 PRECURSOR/OsPSK1. OsPSK1 encodes for a sulfated peptide growth factor [25]. Its ortholog in Arabidopsis is involved in the regulation of the quiescent center cell number and of the differentiation and cell size of root cap and root meristem cells [25,26]. In Arabidopsis roots, maintenance of the quiescent center, meristem, and columella stem cell identity and cell division, is also dependent on TYROSYLPROTEIN The transcription factor genes AtPLT1 and AtPLT2 mediate these TPST/RGF responses [27,28]. TPST-dependent activation of RGFs was proposed to define the patterns and expression levels of PLTs that, in turn, regulate stem cell niche maintenance [27]. Because OsPLT1 and OsPLT3~6 were upregulated in CR primordia (S3 Table) and are in the same phylogenetic group as AtPLT1 and AtPLT2 [29], it is tempting to speculate that the OsPLTs DEG, and in particular OsPLT3 that has a similar expression profile with OsPSK1, could be involved downstream of a OsPSK1-mediated cascade that could contribute to CR primordia formation in rice.

Several transcription factors involved in root development are differentially expressed during crown root formation
Because of the large proportion of DEG encoding transcription factors in our dataset (S3 Fig), we investigated the TF differentially expressed during CR development associated to root development ( Table 1).
Interestingly and adding to the putative involvement of OsPLT3 in a OsPSK1-mediated cascade discussed in the previous paragraph, most of the genes of the PLT family (OsPLT1, OsPLT3~6) were upregulated in developing CR ( Table 1). They are all members of the clade A of the PLT family, and known for their specific expression in both primary root and CR meristem initials [29]. Expression pattern analysis of these genes by in situ hybridization showed that they were specifically expressed in different sectors of the CR meristem [29]. This suggests that these 5 TF might be important regulators in different phases of CR primordium development. PLT genes have been demonstrated to play important roles in root branching in Arabidopsis [48]. AtPLT5, along with AtPLT3 and AtPLT7, redundantly regulate the formative cell divisions in lateral root primordia and the proper establishment of gene expression programs leading to the establishment of a new growth axis [48]. AtPLT5 is the closest Arabidopsis phylogenetic relative of OsPLT3, OsPLT4 and OsPLT5; as is AtPLT4 for OsPLT1 and OsPLT6 [29]. In lateral root primordia, AtPLT3, AtPLT5 and AtPLT7 are required for the early and patterned expression of SHORTROOT (AtSHR) and WUSCHEL-RELATED HOMEOBOX 5 (AtWOX5) and the activation of SCARECROW (AtSCR), that specify the identity of the ground tissue and quiescent center cells [48]. By contrast, OsCRL5/OsPLT8, a member of the clade B of the PLT gene family which encodes an APETALA/ETHYLENE RESPONSE FACTOR (AP2/ERF) and possesses a mutant exhibiting a reduced number of CR due to impaired CR initiation [41], was detected in our study as downregulated during CR primordia development ( Table 1). However, this downregulation is coherent with previous RT-qPCR expression profiling experiments where it appeared to be more expressed in stem tissues than in roots [29]. Furthermore, pOsCRL5::GFP reporter expression pattern showed expression not only in CR primordia, but also in the stem cortex at the level where CR primordia are initiated [41]. This lack of specificity in expression may explain its detection as downregulated in the three analyzed considered stages of CR development versus cortex tissues ( Table 1). This may also indicate that this gene plays a key role during CR initiation but is not involved in later steps of CR primordia differentiation. In any case, further functional characterization of OsPLTs genes is needed to clarify their role during CR development.
In Arabidopsis, AtSHR and AtSCR are activated by AtPLTs [48]. Interestingly OsSCR1 and OsSHR1 were upregulated in early CR primordia (Table 1). OsSCR1, coding for a GRASdomain TF involved in the definition and stabilization of the quiescent center in primary roots [49] and in the establishment of cortex and endodermal cell lineages [50] is specifically expressed in the endodermis, whereas OsSHR1 is expressed in the stele. They are involved in quiescent center differentiation and maintenance [49], and they regulate the molecular events establishing the cortex and endodermal cell lineages [50]. Detection of their expression in the early CR developmental stage is coherent with previous observations [51,52] and their role in early stages of lateral root development in rice [16] and Arabidopsis [53]. These data suggest that OsSCR1 and OsSHR1 likely play a critical role during establishment of CR primordia in rice.
Other TF involved in hormonal signaling pathways were also differentially expressed ( Table 1) during CR primordia formation. This is the case of AUXIN RESPONSE FACTOR (ARF) OsARF5 and OsARF24 and ETHYLENE RESPONSE FACTOR (ERF) ERF99/OsEREBP2 and OsERF120. All are known to play a role in general modulation of root architecture and all are downregulated in CR primordia [54], whereas OsERF61, which cucumber homolog is involved in lateral root development is upregulated [31]. LSD1-LIKE ZINC FINGER PROTEIN OsLOL2/OsLSD5, that is upregulated in CR primordia, is known to be involved in the regulation of gibberellic acid synthesis. Interestingly, a crosstalk between gibberellic acid and auxin is involved in CR emergence and elongation in rice [55,56].

Genes encoding chromatin remodeling factors are differentially expressed during crown root formation
In addition to genes encoding TF, 14 genes encoding for chromatin remodeling factors were found among the upregulated DEG, suggesting that they also contribute to the progressive shift in transcriptional regulation we observed during CR formation (S4 Table). Four of those genes coded for SWItch/Sucrose Non-Fermentable (SWI/SNF) complexes including complex B (SWIB) domain containing proteins, and four other were part of the SBF2 family; both are elements of chromatin structure dynamics contributing to the control of root development [57][58][59][60].
CROWN ROOTLESS 6 (OsCRL6) encodes a member of the CHROMODOMAIN HELI-CASE DNA-BINDING (CHD) family protein and is involved in CR initiation during which it contributes to the regulation of meristem cells differentiation [61,62]. It was identified in a coexpression cluster prone to display high expression levels in the three primordia stages, and low in the cortex. Interestingly, most of the OsIAA genes are down-regulated in the crl6 mutant, suggesting that OsCRL6 regulate auxin signaling pathways [34,61,62]. In our data, a number of OsIAA genes were expressed differentially at the 3 different stages, indicating that different auxin signaling elements and pathways are sequentially activated during CR meristem formation and patterning.

Cross-analysis with other rice root related transcriptome data sets pinpoint key genes involved in de novo root formation in rice
In order to identify genes that may have a conserved function between the processes of CR and lateral root formation in rice, we crossed the list of DEG obtained in this study with a list of 204 genes (corresponding to the S5 Table from [16]) associated with lateral root initiation and development, identified from the RiceXpro RXP_4001 root gene expression dataset [16,63,64].
Seventy three DEG of our dataset (43 upregulated, 30 downregulated) crossed with the list of genes suggested to be involved in lateral root initiation and development [16] (S5 Table). Among these were genes encoding TF such as PLT family members or cell wall remodeling enzymes such as OsXTH10 or EXPANSINS (EXP). OsXTH10 is involved in the xyloglucan metabolism, that plays a central role in monocotyledon cell wall restructuring [65]. In rice, XTH gene expression is tissue-specific and growth stage-dependent [66]. In Arabidopsis, auxin accumulation in the tissues overlying the lateral root primordia induces the expression of a set of cell wall remodeling genes, such as that facilitate the growth of the lateral root through the cortex cells [54]. Similarly, the preferential expression of OsXTH1, OsXTH6, OsXTH10 and OsXTH17 genes in the cortex surrounding developing CR might play a role in the restructuring of cortical cell wall. Regarding the EXPANSINS, OsEXPA2, OsEXPA7 and OsEXPB16 were downregulated in CR primordia in comparison with cortex, and were identified as lateral root initiation-and development-specific [16]. OsEXPA2 is detected through in situ hybridization during the development of the seminal root [67].
Interestingly, most of the shared upregulated genes (33/43) are included in the co-expression cluster III of the RXP_4001 dataset which is suggested to play a role in regulating cell division during lateral root and CR initiation and development, particularly in the apical meristem [16]. These shared upregulated genes included in the co-expression cluster III were mainly assigned to the "midnightblue" (20/33, including OsSHR1, previously discussed) and "skyblue" (8/33) coexpression clusters in our study, characterized by a strongest expression level in the early and intermediate primordia stages, respectively. Downregulated genes are included into two other co-expression clusters identified by [16] (16/30 for cluster I including OsXTH10, 12/ 30 for cluster II). In particular, cluster I is proposed to play a specific role in the priming and initiation of the lateral root primordium. Also, 20 out of the 43 upregulated DEG were identified as involved in both lateral root initiation and development, at visible lateral root emergence stage ("overlapping gene" column), such as OsEXPA5 [16].
This cross-analysis showed that among genes likely involved in lateral root initiation and development, 36% (73/204) were also expressed in developing CR primordia. Over these 73 genes, most of the genes upregulated (33/43) are suggested to play a role in regulation of cell division, and a significant part of the genes downregulated (16/30) could play a role in priming and initiation of lateral roots. This suggests that several genes are involved both in CR and lateral root formation in rice, mostly during the post-initiation and cell proliferation stages.
Also, we compared our dataset with the list of 1753 DEG during the 45 hours following the ectopic expression of CRL1 in the crl1 mutant background, which is used to synchronously restore CR formation [14]. Three hundred and twenty-four DEG of our dataset (190 downregulated, 134 upregulated, representing about 8% of the total number of DEG) were induced or repressed for at least one time point during the 45 hours following CRL1 induction (S6 Table). 229/324 of these genes showed expression values that were coherent with the time series, i.e. genes upregulated or downregulated in the time series were also upregulated or downregulated in the developing primordia from our experiment, respectively. This suggested that the expression of these genes is likely under the control of CRL1 during CR formation. Interestingly genes encoding TF, chromatin remodeling factors and cell wall remodeling enzymes were retrieved.

Conclusions
In this study, we used LCM to generate tissue-level transcriptome gene expression atlas at 3 stages of rice CR primordia development between initiation and pre-emergence. Our results reveal a shift in transcriptional and post-transcriptional gene regulation occurring during CR development. Several genes already described for their involvement in CR development were retrieved in our dataset. New genes, such as members of the PLT gene family and OsPSK1, were also identified suggesting the involvement of peptide signal related pathway during CR formation similar to the pathways already described in Arabidopsis during lateral root formation. This suggest a conservation of some mechanisms during post-embryonic root formation in rice and Arabidopsis. Our analysis also pointed out a differential regulation of several cell wall remodeling enzyme belonging to XHT and EXP families during CR formation. Functional analyses will be needed to confirm their implication in this developmental process.
Our dataset constitutes an open valuable resource to identify new genes involved in CR development in rice, in particular via transcriptomic meta-analyses, similarly to what has been done in Arabidopsis for lateral root development using the VisuaLRTC tool [7] that revealed new key regulators of lateral root development that escaped mutant screening-based investigations [8]. Similarly to what has been successfully proposed on other species such as cotton for fiber quality [68], this resource will be also useful to help the identification of candidate genes in QTL related to rice root development.  Table. List of differentially expressed genes crossing with a list of genes related to lateral root initiation and development from crown root (S5 Table, Takehisa