Investigation of Thermomorphogenesis-Related Genes for a Multi-Silique Trait in Brassica napus by Comparative Transcriptome Analysis

In higher plants, the structure of a flower is precisely controlled by a series of genes. An aberrance flower results in abnormal fruit morphology. Previously, we reported multi-silique rapeseed (Brassica napus) line zws-ms. We identified two associated regions and investigated differentially expressed genes (DEGs); thus, some candidate genes underlying the multi-silique phenotype in warm area Xindu were selected. However, this phenotype was switched off by lower temperature, and the responsive genes, known as thermomorphogenesis-related genes, remained elusive. So, based on that, in this study, we further investigated the transcriptome data from buds of zws-ms and its near-isogenic line zws-217 grown in colder area Ma’erkang, where both lines showed normal siliques only, and the DEGs between them analyzed. We compared the 129 DEGs from Xindu to the 117 ones from Ma’erkang and found that 33 of them represented the same or similar expression trends, whereas the other 96 DEGs showed different expression trends, which were defined as environment-specific. Furthermore, we combined this with the gene annotations and ortholog information and then selected BnaA09g45320D (chaperonin gene CPN10-homologous) and BnaC08g41780D [Seryl-tRNA synthetase gene OVULE ABORTION 7 (OVA7)-homologous] the possible thermomorphogenesis-related genes, which probably switched off the multi-silique under lower temperature. This study paves a way to a new perspective into flower/fruit development in Brassica plants.


INTRODUCTION
Flower development, as well as the subsequent fruit formation, is vital to the crop life cycle. Siliques are important to rapeseed (Brassica napus, AACC = 38), the leading source of plant oil worldwide, which offers more than 13% of the global vegetable oil (Geng et al., 2016). In rapeseed, siliques supply nutrients from photosynthesis, transport carbohydrates from the vegetative organs to the seeds, and ensure their development (Shen et al., 2019). In addition, silique number is a crucial component determining seed yield per plant (Yang et al., 2012;Wang et al., 2016;Shen et al., 2019).
The environment can influence various aspects of plants. Temperature, a key environmental factor, affects the growth, development, and geographical distribution of plants, as well as the quality and productivity of crops (Ding et al., 2020). Therein, "thermomorphogenesis" is defined as the effect of temperature on the morphogenesis of plants (Casal and Balasubramanian, 2019). Take soybean (Glycine max) for example; a 4 • C rise in the temperature could increase its stem height more than 3-fold (Sionit et al., 1987); additionally, the temperature raised from 10 to 15 • C could increase the leaf number in wheat (T. aestivum) (Friend, 1965), whereas, in model plant Arabidopsis thaliana, higher temperatures reduced both silique number per plant and seed number per silique (Ibañez et al., 2017).
Similarly, as described in earlier studies (Chai et al., 2019(Chai et al., , 2020a, the multi-silique trait in zws-ms is significantly affected by the environment. Precisely, zws-ms showed stable multi-silique trait in Xindu (with an annual average temperature of 16.2 • C) for successive years; however, when grown in a colder area such as Ma'erkang (with an annual average temperature of 8.6 • C), zws-ms switched off the formation of multi-silique, and all siliques appeared normal. In previous studies, we identified two associated regions on chromosome A09 and C08 and screened out some potential candidate genes by the combination of bulked-segregant analysis and whole-genome re-sequencing; we also analyzed genes differentially expressed between NILs zws-ms (multisilique) and zws-217 (normal silique) from Xindu and selected potential underlying genes based on annotations. However, the transcriptome from Ma'erkang, where the multi-silique morphology is switched off, as well as the comparison of differentially expressed genes (DEGs) between Xindu and Ma'erkang, remains unclear. In other words, genes involved in the thermomorphogenesis of this multi-silique trait require further investigation.
Thus, in this study, we identified DEGs based on transcriptome sequencing (RNA-seq) from Ma'erkang and then compared them with those from Xindu. The variations in DEGs, combined with their annotations and information of orthologs in Arabidopsis, were then analyzed, and then some potential underlying genes related to thermomorphogenesis were identified.

Plant Materials, Growth Conditions, and Sample Collection
The rapeseed line zws-ms and its NIL, zws-217, were kept in the Crop Research Institute, Sichuan Academy of Agricultural Sciences, China. The multi-silique plant was originally discovered in progenies of B. napus × B. rapa and was successively selfpollinated for six successive generations until the homozygous multi-silique zws-ms line was obtained, whereas the singlesilique offspring were continuously backcrossed with zws-ms (current parent) for six generations, followed by six continuous generations of self-pollination until zws-217 with the singlesilique was obtained. The two lines were simultaneously grown in September in Xindu District of Chengdu in the Sichuan Basin under normal environmental conditions (with an annual average temperature of 16.2 • C). Moreover, both lines were also grown concurrently in early June in Ma'erkang, a mountainous area of western Sichuan, with an annual average temperature of 8.6 • C.
In each location, buds were sampled at their budding stage. In Xindu, three randomly selected individual plants of the zws-ms were assigned as T01, T02, and T03, and three plants in the zws-217 line were assigned as T04, T05, and T06, whereas in colder area Ma'erkang, plants in zws-ms were defined as T01 , T02 , and T03 , and plants in zws-217 were defined as T04 , T05 , and T06 . Eight to 10 buds were sampled from each plant and then quickfrozen and stored in liquid nitrogen.

RNA Isolation and the Library Preparation
The total RNA was isolated as Chai et al. (2019)

Transcriptome Sequencing
The transcriptome sequencing (RNA-seq) was performed on the Illumina HiSeq X-ten platform. In-house Perl scripts were used to remove adapter sequences and read containing poly-N, or low-quality reads, to process the initially generated raw reads into clean reads. Clean reads were then aligned to the B. napus "Darmor-bzh" reference genome 1 by using Tophat2 tools (Kim et al., 2013) to screen out the reads with a perfect match or one mismatch for the next investigation.

Differentially Expressed Gene Analysis
DEGs were detected by the DESeq R package (Love et al., 2014). The P-value was adjusted by controlling the false discovery rate (FDR) (Benjamini and Hochberg, 1995), and genes with an adjusted fold change (FC) > 4 (log 2 FC > 2) and an FDR < 0.001 were then identified as DEGs.

Annotation of Genes
Gene Ontology (GO) database 2 and GOseq R packages (Young et al., 2010) were used to provide gene annotations and calculate GO enrichment of the DEGs. The Kyoto Encyclopedia of Genes and Genomes (KEGG) database 3 and KOBAS software (Mao et al., 2005) were used to explore the high-level functions and utilities of the biological system (Kanehisa et al., 2007) and test the statistical enrichment of the DEGs in the various KEGG pathways. The TAIR database 4 was utilized to provide ortholog information from the model plant Arabidopsis. Sequences of genes from rapeseed were blasted on the website, and orthologs in Arabidopsis were then screened out. Orthologs in Arabidopsis were sufficiently annotated and able to give abundant information, as they have been well studied in the model plant and reported.

Multi-Silique Trait Under Normal Conditions and Its Absence in Colder Areas
Compared with its NIL zws-217 with normal siliques, the zwsms line displayed stable multi-pistil phenotype in successive years in Xindu, where a zws-ms plant showed approximately 32-53% of flowers with the multi-pistil trait, and that further developed into a multi-silique trait, precisely, two or three siliques on a carpopodium ( Figure 1A), like we previously described (Chai et al., 2019). When grown in Ma'erkang, a colder place in a mountain area of western Sichuan Province, zws-ms switched off this multi-silique phenotype; in other words, both lines produced normal siliques there ( Figure 1B).

Transcriptome Sequencing (RNA-Seq) From Two Environments
As described previously (Chai et al., 2019(Chai et al., , 2020a, the two lines were grown at Xindu and Ma'erkang. At each place, flower buds from three individual plants from each line were selected randomly for RNA isolation. Then, their total RNA was extracted, and sequencing libraries were generated, followed by being sequenced on a HiSeq X-Ten platform; the validation for the RNA-seq was also confirmed by real-time quantitative polymerase chain reaction (Chai et al., 2019(Chai et al., , 2020a.

Comparison of Differentially Expressed Genes Between Two Environments
Genes in an environment with expression level fold change > 4 (log 2 FC > 2) between zws-ms and zws-217 and FDR < 0.001 were identified as DEGs. In our earlier report (Chai et al., 2019), 129 DEGs were found between zws-ms and zws-217 in Xindu, among which 67 genes were upregulated, whereas 62 were downregulated. In this study, both lines were grown in colder area Ma'erkang, where they did not show any phenotypic differences to each other and were further subjected to RNAseq, and 117 genes were found expressed differentially in stamen and pistils between zws-ms and zws-217 (Table 1), including 63 upregulated and 54 downregulated in zws-ms (Supplementary Table 1). Samples from the two environments generated different DEGs between zws-ms and zws-217, but in either environment, chromosome A09 and C08, where two associated regions (Chai et al., 2019) underlying this trait were identified, provided most DEGs: in Xindu, 16 (12.4%) and 30 (23.26%) genes were located on chromosome A09 and C08, respectively, whereas 7 (5.98%) and 23 (19.66%) genes were found on chromosome A09 and C08, respectively, of samples from Ma'erkang.
Further analysis found that among these DEGs, there were some expressed only in zws-ms or zws-217, which were assigned "line-specific expressed genes" in this study. Herein, we discovered 21 genes that were line-specific expressed in zwsms and 25 in zws-217 from Xindu (Supplementary Table 2), whereas from Ma'erkang, 18 and one line-specific expressed genes were detected in zws-ms and zws-217, respectively (Supplementary Table 3).
Moreover, the comparison of DEGs between two environments identified 33 DEGs that represented the same or similar expression trends (Supplementary Table 4). In other words, these genes were up-(or down-) regulated in both environments, and they were assigned as "stable DEGs, " whereas the other genes showed different expression trends, which meant zws-ms and zws-217 displayed differential expression level of those genes when grown in Xindu, but no obvious difference while grown in Ma'erkang (Supplementary Table 5). These genes were defined as "environment-specific DEGs" in this study.

Annotations for the Differentially Expressed Genes
As mentioned earlier, we divided these DEGs into two classifications (stable and environment-specific DEGs) and then analyzed their GO annotations. GO terms were usually divided into three categories: biological processes, cellular components (CCs), and molecular functions (MFs). In classification (Figure 2A) for 33 stable DEGs, the biological process terms with the highest levels of enrichment included "metabolic process (GO:0008152)" and "cellular process (GO:0009987), " involving 17 and 16 DEGs; in the CC category, most enriched terms were "cell (GO: 0005623)" and "cell part (GO: 0044464), " whereas in the MF category, "catalytic activity (GO: 0003824)" and "binding (GO: 0005488)" accounted the top enriched terms, containing 13 and 10 DEGs, respectively. In the second classification ( Figure 2B) for 96 environment-specific DEGs, the terms "single-organism process (GO:0044699)" and "cellular process (GO:0009987)" got the highest numbers of a gene, at 38 and 37, respectively; the CC and MF categories showed similar most-enriched terms to those in the first classification.

DISCUSSION
Morphology of higher plants results from the interaction of genotype and environment. The multi-silique trait in rapeseed line zws-ms was found stable in Xindu for generations but absent in Ma'erkang where the climate is colder (Chai et al., 2019). In each place, there are some genes expressed differentially between zws-ms and its NIL zws-217. In Xindu, these DEGs may be the causal factors that distinguish zws-ms from zws-217 in the multi-silique trait, whereas in Ma'erkang, some of those genes turn to no significant difference in expression level between the two lines, which may be regulated by environment and switch off the multi-silique trait in zws-ms. To investigate the potential environment-regulated genes, which switch on/off the multi-silique trait, we compared DEGs generated in Xindu with those in Ma'erkang.
In our earlier report (Chai et al., 2019), we investigated the potential causal variation between zws-ms and zws-217 by the whole-genome re-sequencing, which then led to the identification of two associated regions on chromosome A09 and C08, respectively, as well as some candidate genes within them. It was based on the genomic level conferring inherent and stable genetic variations between the multi-silique and normalsilique lines, which would not be affected by environmental factors. The success of NIL construction, which conferred high genetic similarity between zws-ms and zws-217 and made them different from each other only in this multi-silique trait, ensured the accuracy of these studies. However, the expression level of transcripts can obviously change with external factors. In Chai et al. (2019), we also performed an RNA-seq and found some DEGs between zws-ms and zws-217 in Xindu, where they distinguished from each other in this multi-/single-silique trait; thus, we found some potential causal genes. However, which gene(s) switched off the multi-silique in the colder area remained unclear. Thus, in this study, we added another RNA-seq based on two lines grown in Ma'erkang, where the multi-silique disappears in zws-217, and both lines display normal siliques. By comparing the DEGs from the two environments, we can find out the changed factors.
Thus, among the 129 DEGs from Xindu, we found 33 were stable DEGs. It is worth noting that genes such as BnaA01g10540D, which was downregulated in both Xindu (log 2 FC = −4.53) and Ma'erkang (log 2 FC = −2.33), were defined as "having the same expression tendency, " whereas genes such as BnaA09g06740D, which was line-specifically expressed in zwsms (log 2 FC = +∞) from Xindu and strongly upregulated in it (log 2 FC = 7.57) from Ma'erkang, were defined as "having similar expression tendency." Two genes of them, BnaC08g39130D and BnaC08g42280D, got important annotations. BnaC08g39130D was line-specifically expressed in zws-ms and associated with "ovule development (GO:0048481)." Its ortholog in Arabidopsis, AT1G14980, encodes chaperonin 10, and this type of protein can be involved in physiological processes such as plant seed abortion (Hanania et al., 2007). BnaC08g42280D was strongly downregulated in zws-ms and annotated to "vegetative to reproductive phase transition of meristem (GO:0010228)." This implies that BnaC08g42280D probably serves as an inhibitor of the formation of multi-silique. Its ortholog, AT1G10930, is RECQ4A. Its mutant is hypersensitive to ultraviolet light and sensitive to methyl methanesulfonate. The two genes confer intrinsic and stable variations in zws-ms and are not influenced by the environment.

Gene ID GO annotation KEGG pathway
Stable DEG BnaC08g39130D copper ion binding (GO:0005507); calmodulin binding (GO:0005516); ATP binding (GO:0005524); mitochondrion (GO:0005739); cytosol (GO:0005829); gluconeogenesis (GO:0006094); glycolytic process (GO:0006096); protein folding (GO:0006457); tryptophan catabolic process (GO:0006569); response to heat (GO:0009408); response to cold (GO:0009409); chloroplast thylakoid membrane (GO:0009535); chloroplast stroma (GO:0009570); response to high light intensity (GO:0009644); response to salt stress (GO:0009651); chloroplast organization (GO:0009658); indoleacetic acid biosynthetic process (GO:0009684); chloroplast envelope (GO:0009941); isopentenyl diphosphate biosynthetic process, methylerythritol 4-phosphate pathway (GO:0019288); cysteine biosynthetic process (GO:0019344); response to endoplasmic reticulum stress (GO:0034976); response to hydrogen peroxide (GO:0042542); response to cadmium ion ( Frontiers in Genetics | www.frontiersin.org  Aiming to investigate the thermomorphogenesis-relative genes, we then turned to those environment-specific DEGs, which were defined as those genes such as BnaA05g21710D, which showed significant upregulation in Xindu (log 2 FC = 3.12) but no obvious change in Ma'erkang (log 2 FC = 0.93), that were considered as "environment-specific DEGs" herein. Among the 96 environment-specific DEGs, nine were noteworthy. BnaA09g45320D shared the same ortholog (AT1G14980) with BnaC08g39130D, but unlike the latter, it was only expressed in zws-217 in Xindu and showed no significant difference in Ma'erkang. Its annotation may explain this in response to heat (GO:0009408) and response to cold (GO:0009409). Moreover, it was also annotated to "ovule development (GO:0048481)." Take chaperonin 21 as an example: it was found differentially expressed in seedless and seeded grapes; its silencing resulted in seed abortion in transgenic tobacco (Nicotiana benthamiana) and seedless fruits in transgenic tomato (Lycopersicum esculentum) (Hanania et al., 2007). In Xindu, BnaA09g45320D was linespecific in the normal line, whereas in Ma'erkang, zws-ms and zws-217 showed no significant difference in it. This implied a potential relevance with the multi-silique. The ortholog of BnaC08g41780D, AT1G11870, encodes a seryl-tRNA synthetase, i.e., OVA7, of which disruption can result in ovule abortion in Arabidopsis (Berg et al., 2005). When zws-ms represents multisilique in Xindu, the expression of BnaC08g41780D was not detected in it, which indicates a possibility of it modifying the trait in a warm area. Thus, these two genes, due to the consistency of their expression verities with the changes of environment, as well as their annotations, are considered the most important thermomorphogenesis-related genes of all the candidates.
Although the other environmental-specific DEGs also implicated some indirect clues: BnaA09g44370D and BnaC08g40320D were annotated to MYB-like or MYB/SANTlike DNA-binding domain; and their orthologs, AT1G19000 and AT1G13450, were both homeodomain-like superfamily proteins. Some MYB proteins were found regulated by phytochromeinteracting factor 4, an important thermomorphogenesis factor in Arabidopsis (Wang et al., 2018). In addition to flower color, MYB transcription factors are also involved in the regulation of flower development: transgenic tobacco overexpressing MdMYB3 from apple (Malus × domestica) got longer peduncles of flowers and styles of pistils (Vimolmangkang et al., 2013). BnaA10g07970D and Cole_newGene_2073 were identified as heat shock proteins, which respond to temperature and are involved in signal transduction, protein trafficking, protein degradation, maintaining protein homeostasis, and so on (Raman and Suguna, 2015), but whether/how they are related to flower/fruit shape needs further research. AT2G21660 is homologous to BnaAnng35580D. It is ATGRP7 and encodes a small glycine-rich RNA binding protein. Loss-of-function mutations are late flowering in a non-photoperiodic manner. BnaC08g29060D was annotated to "photoperiodism, flowering (GO: 0048573)" and identified as homologous to AT1G12820, which is auxin signaling F-box 3 and involved in primary and lateral root growth inhibition (Dong et al., 2006). AT1G09970 is an ortholog of BnaC08g42450D and identified as RECEPTOR-LIKE KINASE 7. According to the existing knowledge, it is involved in the control of germination speed and oxidant stress tolerance.
There is another sort of gene, which was not discussed. Although these genes got annotations regarding environmentresponse, they showed no difference resulting from environmental changes. Take BnaC08g39990D for example; it was annotated to "response to cold (GO:0009409), " which seemed potentially related to what we were seeking. However, further analysis indicated that it was always upregulated in zws-ms in both Xindu and Ma'erkang, showing no variety as the environment changed. Thus, its expression level had a low correlation with environmental factors acting on the switching on/off of the multi-silique formation. So, this kind of gene was not emphasized herein.
Rapeseed is not the first crop conferring the multi-pistil trait reported. In fact, wheat contributes relatively sufficient studies about it. Researchers reported their multi-pistil wheat materials, in which the trait can be controlled by a recessive gene, two recessive non-complementary genes, or a single dominant gene, and the loci were located in different chromosomes (Zhu et al., 2019). Even so, no underlying gene in wheat has not been cloned so far (Zhu et al., 2019).
Notably, the multi-pistil trait was also reported in sweet cherry (P. avium L.), co-regulated by genes PaMADS3/4/5, members Frontiers in Genetics | www.frontiersin.org of the MADS-box family Wang et al., 2019). Interestingly, that was also induced by temperature: it occurred more frequently in the warm region of Shanghai than in the cool region of Dalian. This coincides with our discovery from rapeseed. Notably, there are many environmental factors, but the temperature is the most significant one, varying from Xindu to Ma'erkang. Although some other factors, such as photoperiod, are also vital for plant development, there has been no report showing it can change pistil number in flower; it is known to regulate flowering time, as well as some physiological characters such as hypocotyl length, cotyledon angle, chlorophyll content, etc., so far. Instead, a change in temperature, as mentioned earlier Wang et al., 2019), can promote the multipistil formation (in sweet cherry). Thus, we suppose it is the temperature that would switch on/off the multi-silique herein.
The multi-silique phenotype means different things to different plants. For some crops, the multi-silique phenotype was considered as an advantage: in wheat, for example, it is supposed to have the potential to increase yields (Wei, 2017;Zhu et al., 2019); on the other hand, for some other plants, such as sweet cherry, the multi-pistil is regarded as a disadvantage physiological disorder decreasing its commercial value . As to rapeseed, although how to use the multi-silique line to increase the seed yield, as well as the benefit of commercial use, still needs more investigations. We should explain that the mountain area (cold region such as Ma'erkang) is not the major rapeseed-producing area in Sichuan Province; we grow rapeseed there mainly for scientific purposes: the shuttle breeding makes it possible to produce two generations in 1 year, accelerating the population construction processes. Collectively, this study at least provides a new perspective into flower/fruit development in higher plants. In the next stage, the functional validation of these candidate genes will be carried out by investigating the phenotypic, biochemical, and molecular changes in the transgenic plants (both overexpression and knockdown).

CONCLUSION
By comparing DEGs between zws-ms and zws-217 in Xindu with those in Ma'erkang and referring to the gene annotations, we selected BnaA09g45320D (chaperonin 10-homologous) and BnaC08g41780D (OVA7-homologous) as the possible thermomorphogenesis-related genes, switching on/off the multi-silique under different environments.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI and accession number PRJNA736189.

AUTHOR CONTRIBUTIONS
LC and LJ conceived the experiment. LC and JZ performed the research. HL, CC, JJ, and BZ contributed to data analysis. LC wrote the manuscript. LJ and LW reviewed and revised the manuscript. All authors reviewed and approved this submission.