Transcriptomic analysis of Medicago truncatula calli with MtWOX9-1 overexpression Транскриптомный анализ каллусов Medicago truncatula со сверхэкспрессией гена MtWOX9-1

Somatic embryogenesis (SE) is the development of embryo-like structures from somatic plant tissues. This process rarely can be observed in nature, but for many plant species, in vitro protocols are developed, which allow to obtain somatic embryos formation directly from tissues of plant explant or from the embryogenic callus. SE is widely used for plant propagation and transformation; therefore, the search for SE stimulators and revealing of the mechanisms of their functioning are very important for biotechnology. Among the SE regulators, proteins of the WOX family play significant roles. WOX (WUSCHEL-RELATED HOMEOBOX) is a homeodomain-containing transcription factor family. Different WOX genes function in different plant organs and tissues, maintaining meristem activity and regulating cell proliferation and differentiation. Recently, we have shown that transcription factor MtWOX9-1, belonging to the WOX family, can stimulate SE in the Medicago truncatula callus culture. In this research, transcriptomic analysis of highly embryogenic calli with MtWOX9-1 overexpression was performed in comparison to wildtype calli. It was shown that MtWOX9-1 overexpression led to the activation of several groups of genes, including genes related to cell division, tissue diff erentiation, and seed development. Enriched GO pathways included several groups related to histone methyltransferase activity as well as DNA methylation and chromatin binding, suggesting major epigenetic changes that occur in call overexpressing MtWOX9-1 . Using Medicago Truncatula Gene Expression Atlas, we also iden tified a group of genes coding for transcription factors that were both coexpressed with MtWOX9-1 in different plant organs and differentially expressed in samples. These genes are putative targets of MtWOX9-1, and they may act in the same pathway with this regulator during SE.


Introduction
Somatic embryogenesis (SE) is a process of regeneration by which plants use somatic cells to grow embryo-like structures, which eventually can give rise to the new plant. This process isn't observed often in nature, but when plant explants are cultivated in vitro, several factors can induce direct SE or SE from callus tissue. Such factors include specific hormones, the concentration of nitrogen compounds (Reinert at al., 1967), the stress impact (Nic-Can et al., 2016), etc. Most of existing methods of SE induction in vitro include treatment with hormones and mechanical injury, and it is supposed that SE acts as the mechanism of defense against stressful in vitro conditions. Somatic embryo development occurs in general through the same stages as development of zygotic embryo, and therefore it is used as the model for studying embryogenesis. SE also has a lot of biotechnological applications in transformation of plants, artificial seeds production and micropropagation.
The way by which embryogenic cells are chosen among explant or callus cells is not fully investigated. In Medicago truncatula, one of SE model objects, somatic embryos are often derived from mesophyll cells near the damaged surface, which have to dedifferentiate, but some are derived from the stem-like vascular procambium cells (Wang et al., 2011;Rose, 2019). Auxin gradients are shown to play key role in determination of cells that will give rise to the embryo (Su et al., 2009).
Dedifferentiating cell, which will be capable to form somatic embryo later, undergoes a number of changes, including mitochondrial fusion and increase in peroxisomes (Tiew et al., 2015) and P-bodies (RNA processing bodies) numbers (Bhullar et al., 2017). The first one is used to provide a kind of quality control for mitochondrial populations in new generations (Rose, McCurdy, 2017). The second one is a kind of stress response, whereas the third one plays role in posttranscriptional gene regulation and cell reprogramming.
Genetic cascade that induces dedifferentiation of cell and development of somatic embryo is not fully uncovered. At the present time correlation with SE has been established for LEAFY COTYLEDON1 (LEC1), BABY BOOM (BBM), AGL15, SOMATIC EMBRYOGENESIS RECEPTOR KINASE (SERK) and other genes (Fehér, 2015). Among the regulators of SE, proteins of WOX family also play important roles. WOX (WUSCHEL-RELATED HOMEOBOX) is a homeodomain-containing transcription factors (TFs) family. Different WOX genes function in different plant organs and tissues, maintaining meristem activity and regulating cell proliferation and differentiation.
The most well-studied family members are WUS and WOX5 genes, which are expressed in the organizing and quiescent center cells of the shoot and root apical meristem (SAM and RAM), respectively, and regulate their development (Laux et al., 1996;Sarkar et al., 2007). The mechanism of WUS action in stem cell niche is related to stimulation of cytokinins activity, by which it represses the differentiation of SAM cells (Leibfried et al., 2005). WUS is also expressed in the floral meristem, where it stimulates AGAMOUS expression, providing termination of floral meristem activity (Lenhard et al., 2001).
WOX5 is functional analog of WUS gene in root. It is expressed in the quiescent center in the root apical meristem (Sarkar et al., 2007) and in different irregular meristems, such as nodule meristems (Osipova et al., 2012) and meristem-like structures of agrobacterial and spontaneous tumors Vinogradova et al., 2015).
The participation in zygotic and somatic embryogenesis was demonstrated for many WOX family genes. For example, WUS is an important stimulator of SE in different species, such as Arabidopsis thaliana (Su et al., 2009), Capsicum chinense (Solís-Ramos et al., 2009) and Gossypium hirsutum (Bouchabké-Coussa et al., 2013;Xiao et al., 2018). WOX5 is also involved in SE, participating in RAM development in somatic embryos (Su et al., 2015). Expression of WOX1 and WOX3 homologs was observed during SE process in different objects, such as C. chinense (Valle-Gough et al., 2015), Vitis vinifera (Gambino et al., 2011) and Picea abies (Alvarez et al., 2015). The WOX11 and WOX12 genes are expressed in the early stages of callus and adventitious roots development in Arabidopsis, stimulating the cambium cells proliferation (Liu et al., 2014). The expression of WOX11 homolog during SE was shown in V. vinifera (Gambino et al., 2011).
The genes WOX2, WOX8, and WOX9 play an important role in the zygotic embryogenesis, defining the differentiation of specific embryo domains: apical (WOX2), central (WOX9) and basal (WOX8) (Breuninger et al., 2008). In A. thaliana, WOX2 and WOX8 expression was detected in egg cell, though the Nicotiana tabacum homologs of these genes were shown to be de novo transcribed in zygote right after fertilization (Zhou et al., 2018). The WOX2 homolog in Larix decidua is expressed during early embryogenesis, both somatic and zygotic (Rupps et al., 2016). The WOX9 homologs are microspore embryogenesis markers in Brassica napus (Malik et al., 2007) and SE markers in V. vinifera (Gambino et al., 2011). Besides, expression of the WOX9 homolog MtWOX9-like was demonstrated during SE in M. truncatula (Kurdyukov et al., 2014).
In our previous studies, three new M. truncatula genes of the WOX family, that are expressed during SE, were found: MtWOX9-1, MtWOX11-like and STENOFOLIA (MtWOX1) . It was further shown that overexpression of STF or MtWOX9-1 stimulates the emergence of somatic embryos (Tvorogova et al., 2016(Tvorogova et al., , 2019. In the present work, we concentrated on studying the functions of the MtWOX9-1 gene, analyzing how its overexpression affects the expression profile of embryogenic callus.

Materials and methods
Plant growth, cultivation and sample collection. M. truncatula line R-108 (Hoffmann et al., 1997) and transgenic line with MtWOX9-1 overexpression were used in the analysis.
Line with MtWOX9-1 overexpression, obtained through argobacterial transformation of R-108 line plants with pMDC32 vector containing MtWOX9-1 coding sequence under the control of 35S promoter (Wolabu, 2015), was kindly provided by the laboratory of Dr Million Tadege.
M. truncatula seeds were sterilized in sulfuric acid for 10 minutes, rinsed 10 times with sterile water, and then put onto 1 % agar and left to germinate at 4 °C for 7 days. After germination, seedlings were transferred in the soil (Terra Vita, Russia) mixed with vermiculite (2:1). Plants were grown at 21 °C at 16 h photoperiod. Before in vitro cultivation, leaves of 30 day old plants were sterilized in 70 % ethanol for 1 minute, then in 50 ml solution of 0.5 % hypochlorite with two drops of Tween-20 for 10 minutes, and then rinsed 5-7 times with sterile water. In vitro cultivation and obtaining of embryogenic calli was performed as described previously (Tvorogova et al., 2016).
Two biological replicates for R-108 wildtype calli (wt from this point onward) and two biological replicates for MtWOX9-1 overexpressing calli (w9o from this point onward) were taken at 35th day of cultivation (5th day of cultivation of hormone-free medium). Plant material was divided into two parts, the first for transcriptomic analysis and the second for quantitative PCR (qPCR) analysis.
Transcriptome sequencing and bioinformatic processing. RNA extraction, library preparation and sequencing was performed by the Genoanalitica company (Moscow, Russia). Total RNA was extracted from calli with Trisol reagent and PureLink RNA Micro Kit (Invitrogen) according to manufacturer instruction. Quality was checked with BioAnalyser and RNA 6000 Nano Kit (Agilent). PolyA RNA was purified with Dynabeads ® mRNA Purification Kit (Ambion). Illumina library was made from polyA RNA with NEBNext ® Ultra™II RNA Library Prep Kit for Illumina ® (NEB) according to manual. Sequencing was performed on HiSeq1500 with 50 bp read length.

Results
To analyze the mechanisms of MtWOX9-1 functioning during SE, 35 day-old wt and w9o calli of M. truncatula were obtained. At this stage, somatic embryos, visible as green spots on callus surface, started to appear on w9o calli (Suppl. Fig. 1, a) 1 , but not on wt calli which usually start form embryos later during cultivation. The RNA was extracted, reverse transcribed and sequenced from embryogenic wt and w9oe calli (two biological replicates for each variant). According to qPCR analysis, expression level of MtWOX9-1 was several thousand times higher in w9o samples than in wt samples (see Suppl. Fig. 1, b). 4 complementary DNA libraries were sequenced with an average depth of approximately 15 millions of reads. After trimming and ribosomal RNA removing, about 12 millions of 35-bp reads were taken for analysis. After alignment with HISAT2, about 75 % of reads in each sample were uniquely mapped. Reads were counted by StringTie, and correlation analysis performed for DESeq normalized counts demonstrated high correlation between biological replicates (see Suppl. Fig. 1, c).
After the analysis of differential expression with DESeq package and imposition of 0.01 adjusted p-value and 1.0 log2 fold change cutoffs, 3133 genes out of 51628 analyzed were found to be differentially expressed (Suppl. Table S1), with 1608 and 1525 up-and downregulated, respectively, in w9o calli in comparison to wt calli. qPCR expression analysis of two differentially expressed genes (DEGs) supported transcriptome analysis data (Suppl. Fig. 2).
To find new potential stimulators and repressors of SE, we assessed expression levels of genes coding TFs among DEGs (selection of DEGs with GO annotation number 0006355, "regulation of transcription, DNA-templated"). We also added there five DEGs from WOX family, which happened not to be included in this GO group, but, according to numerous data, should have TF function. We found 173 DEGs coding TFs of which 94 gene was upregulated, including several TFs from NF-Y family, B3 domain TFs, MADS-box TFs etc (Suppl .  Table S2).
Transcriptomic analysis of Medicago truncatula calli with MtWOX9-1 overexpression As well as w9o calli appear to be more embryogenic than wt ones according to our data (Tvorogova et al., 2019), we supposed that MtWOX9-1 should specifically stimulate expression of genes related with embryo development and development of pods. To check this hypothesis, we selected the genes with specific expression in seedpod (expression level in seedpod is at least 4 times higher than average expression level measured in seedpod, leafblade, nodule, root, flower and bud) from the M. truncatula genome database (Krishnakumar et al., 2015). Among these genes, only 3.7 % (44 genes) were differentially expressed in our experiment, but 75 % of this group of DEGs were upregulated (Fig. 1), which is consistent with our hypothesis.
We also performed gene enrichment analysis with GSEA-Base package using GO terms (Ashburner et al., 2000). Upregulated and downregulated genes were found to be enriched with genes from 17 and 21 GO "Molecular Function" groups, respectively (Fig. 2, Suppl. Fig. 3, Suppl. Table S3).
Interestingly, upregulated GO pathways included several groups related with histone methyltransferase activity as well as DNA and chromatin binding, suggesting major epigenetic changes occuring in w9o calli. We analysed A. thaliana homologs of upregulated histone methyltransferases and found that most of them are responsible for histone repressive marks associated with DNA methylation (Table). 49 and 40 "Biological Process" GO terms were overrepresented among upregulated and downregulated genes, respectively (Fig. 3, Suppl. Fig. 4, Suppl. Table S4). In accordance with the data above, genes from DNA methylation and gene silencing GO groups were found among upregulated genes.
To find the genes which may work together with Mt-WOX9-1, we performed coexpression analysis using Benedito et al. (2008) data from Medicago Truncatula Gene Expression Atlas and WGCNA R package (Langfelder, Horvath, 2008). We found 55 coexpression modules, containing from 35 to 4914 different genes and isoforms (Suppl. Table S5). The MtWOX9-1 module contained 2337 genes, out of which 225 genes were differentially expressed in our calli samples according to transcriptome analysis (Suppl . Table S6). We searched for TF genes among this group using selection of genes with GO annotation number 0006355, "regulation of transcription, DNA-templated". We also added two DEGs from WOX family (including MtWOX9-1 itself), which happened not to be included in this GO group, but, according to numerous data, should have TF function. The resulting TF gene list contained 18 genes, up-or downregulated in w9o calli, including genes for three B3 domain proteins, several homeobox-containing factors, etc (Suppl .  Table S7).

Discussion
In this study, we analyzed transcriptome of embryogenic calli, affected by overexpression of SE stimulator MtWOX9-1. It is important to note that in this study we used a single line with MtWOX9-1 overexpression, therefore one cannot exclude the possibility that observed effects are related with specific location of transgenic insert in this line. However, according to our results (Tvorogova et al., 2019), the positive effect of MtWOX9-1 overexpression on embryogenic capacity was observed after several independent transformation events, therefore we can assume that the majority of the effects demonstrated in this study are not unique for specific transgenic line.
Analysis of gene groups activated in w9o calli allowed us to suggest that MtWOX9-1 overexpression and probably SE itself are associated with major epigenetic changes in embryogenic callus, such as DNA and histone methylation. Such changes are common for in vitro cultures (Kumar, Van Staden, 2017). Most of observed upregulated chromatin-related pathways were found to be repressive which is probably due to the start of differentiation of callus cells forming specific embryo tissues. Interestingly, such pathways as postembryonic development and seedling development were upregulated, suggesting important dif feren ces between zygotic embriogenesis and somatic embryogenesis, lacking dormancy stage.
As it was mentioned above, SE is the stress-induced process: temperature, wounding, starvation, heavy metal ions, and osmotic stress can lead to dedifferentiation of somatic cells and to other changes. In support of this, we also found some upregulated GO groups associated with stress response, including, for example, proline dehydrogenase and response to abscisic acid. Proline dehydrogenase enzyme is an important component of plants pathogen defense system which contributes to the hypersensitive response (HR) and disease resistance, and promotes the accumulation of reactive oxygen  (Cecchini et al., 2011;Monteoliva et al., 2014). Several downregulated GO groups associated with stress response were also found, including phenylalanine ammonialyase activity (Wada et al., 2014) and chitin binding groups. Several other SE and cell division associated pathways were also found to be upregulated in w9o calli, including phosphoenolpyruvate carboxylase and pectin methylesterase (PME) activity, cytokinetic process etc. Phosphoenolpyruvate carboxylase modulates cell division and elongation of cotton fibers, the ovule epidermal cells formed during flowering process (Li et al., 2010). Changes in the methylesterification status of pectins have been associated with cell wall remodeling that occurs during diverse plant developmental processes (Levesque-Tremblay et al., 2015). It was shown that PME expression level increased in heart-torpedo embryos and mature cotyledonary embryos (Pérez-Pérez et al., 2018). This tendency is kept during somatic embryogenesis process among both woody and herbaceous plants (Solis et al., 2016). Besides their roles in plant development, PMEs are also involved in stress response and pathogen defense (Ma et al., 2013).
WOX genes are master regulators of development, and their direct targets often encode for other TFs. The results of our coexpression analysis allowed to identify several TF genes which may be the direct targets of MtWOX9-1. Therefore, the next step of research will be to check the relations between MtWOX9-1 and these genes using molecular biology methods.

Conclusion
The observed differences between w9o and wt calli and control can be considered as the specific effect of MtWOX9-1 overexpression but, on the other hand, they may result from SE process itself, which tends to start earlier in w9o calli (Tvorogova et al., 2019). Thus, these data may be useful both for MtWOX9-1 target search and for the search for new SEassociated genes and SE stimulators.