Early cold stress responses in post-meiotic anthers from tolerant and sensitive rice cultivars

Background Rice grain production is susceptible to a changing environment that imposes both biotic and abiotic stress conditions. Cold episodes are becoming more frequent in the last years and directly affect rice yield in areas with a temperate climate. Rice is particularly susceptible to cold stress during the reproductive phase, especially in anthers during post-meiotic stages which, in turn, affect pollen production. However, a number of rice cultivars with a certain degree of tolerance to cold have been described, which may represent a good breeding resource for improvement of susceptible commercial varieties. Plants experiencing cold stress activate a molecular response in order to reprogram many metabolic pathways to face these hostile conditions. Results Here we performed RNA-seq analysis using cold-stressed post-meiotic anther samples from a cold-tolerant, Erythroceros Hokkaido (ERY), and a cold-susceptible commercial cultivar Sant’Andrea (S.AND). Both cultivars displayed an early common molecular response to cold, although the changes in expression levels are much more drastic in the tolerant one. Comparing our datasets, obtained after one-night cold stress, with other similar genome-wide studies showed very few common deregulated genes, suggesting that molecular responses in cold-stressed anthers strongly depend on conditions and the duration of the cold treatments. Cold-tolerant ERY exhibits specific molecular responses related to ethylene metabolism, which appears to be activated after cold stress. On the other hand, S.AND cold-treated plants showed a general downregulation of photosystem I and II genes, supporting a role of photosynthesis and chloroplasts in cold responses in anthers, which has remained elusive. Conclusions Our study revealed that a number of ethylene-related transcription factors, as putative master regulators of cold responses, were upregulated in ERY providing promising candidates to confer tolerance to susceptible cultivars. Our results also suggest that the photosynthesis machinery might be a good target to improve cold tolerance in anthers. In summary, our study provides valuable candidates for further analysis and molecular breeding for cold-tolerant rice cultivars.


Background
Rice is one of the major staple cereals in the world, providing essential caloric requirement for nearly half of the population (Khush 2005). Production of rice will need to be increased in coming decades to satisfy a steadily increasing demand from a fast-growing world population (Anderson et al. 2004). A sustainable increase in productivity requires intensified efforts to develop varieties with improved yield and greater stress tolerance. Global climate change projections pinpoint increasingly episodes of extreme temperatures. Cold temperatures affect rice growth during its whole life cycle, and especially during the seedling and reproductive stages (Andaya and Mackill 2003;Yoshida et al. 1996). In temperate regions, where rice cultivation experiences cold episodes, it is cultivated in a time-frame window of three-to-four months depending on the cultivar. In order to avoid cold episodes during seedling stages sometimes a delay in seed sowing is needed. Hence, rice plants during booting stage may encounter cold episodes at the end of the summer, highlighting the importance of planting time for the proper management of this crop (Biswas 2019).
In the last two decades an increasingly number of studies have helped to elucidate how plants respond to cold stress and the genes/proteins involved in this response. Cold is primarily sensed by cell membranes where a cold sensor triggers an influx of calcium, in order to activate downstream responses (Ma et al. 2015). The following cold signal cascade involves a number of calcium-binding proteins, mitogen-activated protein kinases and, finally, transcription factors that reprogram the transcriptome in the nucleus (Guo et al. 2018).
Cold stress produces dramatic outcomes during postmeiosis, affecting microspore generation and pollen development . Molecular analyses helped to understand how genomes reprogram the transcriptome and the proteins synthetized or silenced, in order to give a full molecular response to the stress. Six genome-wide transcriptional studies from cold-stressed anthers or spikelets during booting stage were reported (Bai et al. 2015;Ishiguro et al. 2014;Oda et al. 2010;Shimono et al. 2016;Suzuki et al. 2015;Yamaguchi et al. 2004). The obtained results from those studies were diverse, probably due to the use of different conditions, type and duration of treatments carried out. Nevertheless, samples analysed in those studies were collected after several days of cold stress conditions, suggesting that results obtained may not reflect the early molecular events that plants may accomplish in order to face the stress.
Here we performed RNA-seq experiments to study the genome-wide expression profiles in one-night cold-stressed post-meiotic anthers from tolerant (Erythroceros Hokkaido, ERY) and commercial susceptible (Sant'Andrea, S.AND) cultivars. Interesting genes and pathways putatively involved in cold-stress responses were identified. A detailed study of datasets obtained from this work allowed us to propose a candidate list of molecular processes and genes that might be responsible for the early differential physiological response to cold stress and, hence, be valuable for molecular breeding.

Results
Spikelet fertility in contrasting cultivars exposed to coldstress In order to study how different rice cultivars respond to low temperatures at the post-meiotic stage, the Italian accessions Arborio, Carnaroli, S.AND and Volano were chosen for analyses. Two cultivars putatively coldtolerant, ERY and Lijiangheigu, were analysed in parallel together with the reference cultivar Nipponbare. ERY is a cultivar from Poland (Zhao et al. 2011), one of the highest latitudes in the world where rice is cultivated, while Lijiangheigu is a well-known cold-tolerant cultivar at different growth conditions (Ye et al. 2009). All cultivars tested were grown in phytotrons until prebooting stage and subjected to cold stress (8°C), or control conditions (22°C) during the night in five consecutive days. Spikelet fertility was measured for all seven cultivars as the percentage of filled grains compared to total grains counted at maturity. ERY subjected to cold-stress shows spikelet fertility similar to control conditions (around 80%), behaving as a coldtolerant cultivar in our conditions, whereas fertility in Lijiangheigu was significantly affected by cold-stress ( Fig. 1). Among the Italian cultivars that we tested we found S.AND as most cold-susceptible showing a decrease of about half in spikelet fertility (Fig. 1). These results led us to choose ERY and S.AND as reference contrasting cultivars to compare physiological and molecular responses to cold-stress in detail. Fig. 1 Spikelet fertility in cold-stress tolerant and sensitive cultivars. Spikelet fertility as the percentage of filled/empty spikelets after grain maturation is shown. NB, Nipponbare. Bars represent mean ± SEM of at least two biological replicas. Asterisks indicate statistically different mean values compared with the control (Student's t-test, p < 0.05)

Physiological responses to cold stress in S.AND and ERY plants
As was reported previously, one of the most cold-sensitive stages during reproductive development in rice is the early microspore stage. This particular stage often correlates with the pre-booting stage but differences in reproductive development between cultivars might exist. Thus, we wanted to address the precise developmental stage that coincides with meiosis or post-meiosis to study cold responses in ERY and S.AND. We collected spikelets from both cultivars at different auricle distances (AD) between the flag leaf and the penultimate leaf (between − 5/+ 4 cm) during pre-booting stage. Transversal sections of anthers show that samples range between early post-meiotic and late vacuolated stages, corresponding to anther developmental stage 9 (Additional file 1: Figure S1a), between stage 9 and 10 (Additional file 1: Figure S1b and S1d), and stage 10 (Additional file 1: Figure S1c, S1e, and S1f) according to Zhang et al. (2011), before tapetum layer degraded. Next, we exposed ERY and S.AND panicles at AD between − 2 and + 2 cm to cold and control conditions and took samples 1, 3 and 5 days after the beginning of the treatments. Gross morphology of cold-treated panicles for both cultivars did not show any difference compared to their respective controls (Additional file 2: Figure S2A). Similar results were obtained for spikelets (Additional file 2: Figure S2b) and transverse sections of anthers (Additional file 2: Figure S2c). However, a pale green colour can be seen in cold-treated S.AND (Additional file 2: Figure S2b). In fact, many spikelets (around 10-15%) from these plants showed whitish palea and lemma (Additional file 3: Figure S3a), a phenotype that cannot be seen in control samples or in cold-treated ERY plants (Additional file 2: Figure S2). Furthermore, transversal sections of anthers from these spikelets clearly showed different levels of damaged morphology (Additional file 3: Figure S3b and S3c). These results might probably explain part of the decreased fertility observed in S.AND after cold exposure. For further experiments we only used anthers that showed normal anatomy in order to avoid biased results due to degradation of macromolecules in collected samples.

Genes differentially expressed in post-meiotic coldstressed anthers of S.AND and ERY cultivars
In order to compare the early molecular responses to coldstress from tolerant ERY and susceptible S.AND cultivars, anthers were isolated from one day cold-treated panicles specifically coinciding with post-meiotic stages as explained above (samples ERY COLD and S.AND COLD). The same type of tissues was collected from plants subjected to a normal growing temperature (24°C) and used as controls (ERY CONTROL and S.AND CONTROL). Three biological replicates from each treatment were obtained. Total RNA was extracted from these tissues and used for Illumina RNA-seq experiments. More than 20 million reads were obtained from each of the 12 samples and mapped to the Oryza MSU7.0 database using the commercially available CLC Genomics Workbench (CLC bio, Aarhus, Denmark). Gene expression was quantified as reads per kilobase of coding sequence per million reads (RPKM). After normalization by total reads, statistical analyses were carried out (Baggerley's test), and differentially expressed genes (DEGs) from control and cold samples were obtained, setting False Discovery Rate (FDR) < 0.05 and fold change (FC) > 2 as cut-offs (SSTF, statistically significant by two-fold). Nearly the same number of SSTF DEGs was obtained from each cultivar (1017 and 1015 DEGs for ERY and S.AND, respectively), although up-regulated genes were more abundant in ERY cultivar than in the cold sensitive one ( Fig. 2a and Additional file 5: Table S1). Both cultivars share common cold-responsive genes showing that part of the molecular response to cold-stress is similar, as expected. However, many genes were only deregulated in one specific cultivar and might be responsible for the differential cold tolerance during reproductive stages. In order to validate RNA-seq results, six genes with altered expression levels either between controls or between control and cold-stressed samples were chosen as representatives to quantify their expression by qRT-PCR. Relative expression analysis shown in Additional file 4: Figure S4A follow the same pattern obtained by RNA-seq analysis (Additional file 4: Figure S4b), validating the genome-wide transcriptional profiling.
Next, we wanted to understand which categories were overrepresented in the DEGs lists in comparison with the whole rice genome. DEGs were further analysed for gene ontology functional annotations with agriGO analysis tools (Du et al. 2010). Up-regulated DEGs from both cultivars are enriched in genes involved in stimulus responses, nitrogen compound and cellular metabolic processes, among others, belonging to the "Biological Process" category (Additional file 6: Table S2 and Fig. 2b). Also, a significant amount of transcription factors (TFs) in those lists stand for an enrichment in the categories "Molecular Function" and "Cellular component" (Fig.  2b, Additional file 6: Table S2 and Additional file 7: Table S3). We performed a scatterplot of log 2 FC from those TFs belonging to ERY and S.AND's DEGs lists (black circles from SSTF DEGs in both cultivars, grey squares from SSTF DEGs in only one cultivar; Fig. 2c). Strikingly, most of these transcription factor genes were more upregulated in the cold-tolerant ERY than in the susceptible S.AND, showing a positive correlation with tolerance to cold stress. On the other hand, genes down-regulated by cold in ERY do not stand for any enriched ontological group. However, a number of thylakoid-related genes putatively involved in photosynthesis are statistically enriched in the list of cold down-regulated S.AND genes ( Fig. 2d and Additional file 6: Table S2). Details of enriched gene ontology accessions and associated FDRs of DEGs are provided in Additional file 6: Table S2.

Expression of genes involved in anther development and starch accumulation
After meiosis, there is a number of crucial biological processes taking place in anthers to ensure pollen fertility. The transport of nutrients from tapetum to microspores, degeneration of tapetal cells via programmed cell  Table S2. c, scatterplot of log 2 FC from deregulated transcription factors in cold-stress ERY and S.AND anthers. Black circles correspond to SSTF (statistically significant by two-fold) genes in both cultivars, while grey squares correspond to genes in only one cultivar death, and pollen wall development should be well coordinated in order to produce viable pollen. Dozens of genes belonging to different pathways in these anther developmental steps have been characterized (Guo and Liu 2012;Shi et al. 2015). We reasoned that the expression of some of these genes could be deregulated between our cold tolerant and sensitive cultivars. Surprisingly, there were only three SSTF DEGs in ERY dataset, among 45 genes analyzed in both datasets (Additional file 8: Table S4 and Additional file 12: Data S1). Two of those genes showed similar FC in both cultivars while the expression of the third one, the recently characterized TDR INTERACTING PROTEIN 3 (TIP3) (Yang et al. 2019), was more induced by cold stress in ERY than in S.AND cultivar (Additional file 8: Table S4 and Additional file 12: Data S1).
Starch accumulation is another crucial step in the maturation and viability of pollen grains. A cell-wall invertase and two monosaccharide transporter genes were shown to be important for starch accumulation in pollen cells and their expression were affected by cold stress (Oliver et al. 2005). We tested whether the expression of these genes, together with those reported to be involved in starch synthesis (Ohdan et al. 2005) and carbon partitioning (see references in Data S1), was affected in our datasets. We found very few SSTF DEGs in each dataset with moderated changes in their expression (Additional file 8: Table S4). The only gene consistently induced by cold stress in both cultivars was Rice Starch Regulator 1 (RSR1), an AP2/EREBP transcription factor proposed to be a starch biosynthesis regulator (Fu and Xue 2010). On the other hand, only one gene, described as starch synthase IVa (Ohdan et al. 2005) was induced by cold stress in ERY but remained unaffected in S.AND (Additional file 8: Table S4).

Common molecular responses in cold-stressed rice anthers
In order to elucidate common molecular responses in anthers subjected to cold-stress we compared lists of deregulated genes from our work and other published transcriptomic analysis. Six other studies have shown a variable number of differentially expressed genes in cold-stressed anthers or spikelets during young microspore stage, although cultivars, duration of treatment, temperature, etc., were very variable (Additional file 9: Table S5). From these previous studies there were only three datasets available of common cold-induced DEGs (Bai et al. 2015;Ishiguro et al. 2014;Suzuki et al. 2015) to compare with our work. Multiple comparison shows very few common genes between two datasets, where only two genes were common in three datasets: LOC_ Os03g18870, encoding the heat shock protein DnaJ; and LOC_Os12g39630, encoding a calcium/calmodulin dependent protein kinase CAMK-like.49. No common deregulated genes were found in all four transcriptomic analysis (Fig. 3). These results suggest that molecular responses in cold-stress anthers are very variable and probably highly dependent on the experimental setup.

Specific cold-induced molecular responses in ERY and S.AND
To gain further insights into metabolic pathways or other biological processes differentially affected in cold-stress, ERY and S.AND samples we analysed SSTF genes specifically deregulated in each cultivar with Mapman Software (Thimm et al. 2004). Thus, a list of 444 and 442 genes specifically deregulated in ERY and S.AND, respectively, were mapped to the Oryza MSU7.0 database. An overview of deregulated genes falling in different hierarchical categories (Fig. 4) shows: most of the photosynthesis (bin 1), cell wall (bin 10), secondary metabolism (bin 16), and miscellaneous enzyme families (bin 26) genes were specifically downregulated in S.AND; while most of the hormones (bin 17), RNA (bin 27), protein (bin 29), and signalling (bin 30) genes were specifically upregulated in ERY. To test which bins exhibit statistically significant differences we applied Benjamini-Hochberg correction as statistical analyses, provided in Mapman software with p < 0.05. Two main categories were overrepresented in S.AND: bin 1.1 corresponding to "photosynthesis-light reaction", where all genes but one is downregulated (Additional file 10: Table  S6, Fig. 5) and; bin 29 corresponding to "protein", where there is no clear trend in the regulation of these genes. Additional file 10: Table S6 also shows the two main categories overrepresented in ERY: bin 17 and 27.3.3 corresponding to "hormone metabolism" (most of them related to ethylene) and "RNA-regulation of transcription-AP2/EREBP", respectively. Although some reports have shown, on one hand, the involvement of ethylene in male gametophyte development (Kovaleva 2011) or anther dehiscence (Rieu et al. 2003;Wang and Kumar 2007) and, on the other hand, the role of ethylene in cold stress responses, reviewed in (Kazan 2015), no reports have linked ethylene metabolism in anther tissues as part of the molecular responses to cold stress. Several genes putatively involved in ethylene biosynthesis, sensing, signal transduction and related transcription factors have been described (Rzewuski 2008). Table 1 highlights the transcriptionally active set of genes in ERY anthers during post-meiotic stages and how they respond to cold stress. Remarkably, the anther most expressed ACC oxidase genes, that catalyse a committed step in ethylene biosynthesis, OsACO2 and OsACO4 are upregulated by cold stress in ERY. Other proteins involved in the regulation of ethylene receptors (OsRTE2), and in signal transduction (OsEBF1, OsEBF2, and OsEIL2) are also upregulated. Finally, three well-known ethylene response factors are significantly upregulated as well (Sub1C, OsERF1 and osERF3). Table 1 also shows a number of ethylene-related

Discussion
Cold episodes during reproductive stages in plants are becoming more frequent due to many factors in temperate regions, including climatic change, leading to important yield losses in crops. Commercial rice cultivars were selected for high yield but poor resistance to many abiotic or biotic stresses. Many molecular studies were performed in reference rice cultivars providing important insights in the response to stress, although variability lead to the notion that not always conclusions are translatable to other cultivars with agronomic relevance. This work was focused on the cold-susceptible Sant'Andrea, one of the most cultivated commercial varieties from Italy, the main producer of rice in Europe. The Erythroceros Hokkaido cultivar served as the cold-tolerant variety for comparative purposes. This cultivar grows in Poland, one of the most northern latitudes in the world in which rice is cultivated. Its highly tolerant behavior to cold stress points to a good potentiality to serve as breeding resource to introduce tolerance in cold-sensitive commercial cultivars. Available transcriptomic studies from rice anthers subjected to cold stress led to very heterogenous results in terms of observed molecular responses, probably due to diverse conditions in which datasets were produced. In those studies, cold stress was applied at different developmental stages, at different temperatures and conditions (air or cold water). However, one of the most heterogenous factor was the duration of stress until samples were collected. This may explain the unexpectable low amount of common deregulated genes between the different studies. Many of the genome-wide transcriptomic studies were based on samples collected several days after the beginning of the cold treatment, suggesting that observed molecular changes might correspond to secondary or tertiary responses. In this study, RNA-seq data were obtained from samples subjected to cold stress during one night, where probably primary changes in the transcriptome are observed as a response to the stress applied. Changes in expression of a large number of TF encoding genes in our datasets are in agreement with that and provide a good candidate list of putative master regulators of cold responses, and may in the future be further studied to identify those genes that facilitate an increase in cold tolerance in sensitive cultivars. Remarkably, tolerant cultivar ERY showed a stronger molecular response in the activation of several TFs compared to S.AND. It might be due to a transcriptional higher rate of these master regulators, or to a faster response in the regulation of these by for instance reprogramming of the chromatin and thereby triggering the activation of the molecular machineries to successfully face stress before cell death is induced. S.AND plants were susceptible to cold episodes during reproductive stages as it was shown for the premature degeneration of anthers and a visible decay in pigments, giving as a result a decrease in spikelet fertility. Chloroplasts might be impaired in this cultivar since results show a global decay in the synthesis of photosystem components necessary for a successful photosynthesis in anthers. This may reflect the relevance of chloroplasts for a proper male gametophyte development within anthers. The role of chloroplasts in anther development is not clear, although they are present especially in the anther epidermis, endothecium and the middle layers during microsporogenesis (Clément 2001). They may have a role in the synthesis of carbohydrates to feed tapetum, the male gametophyte and pollen grains. ERY did not show a massive downregulation of photosynthesis genes, either because chloroplasts are protected by specific metabolites and the signal triggered to degrade plastids are inhibited, or upstream master genes that regulate transcription of photosynthesis machinery are not impaired.
Finally, ethylene metabolism may be activated in ERY cold-stress anthers as a result of the transcriptional increase of several genes involved in its biosynthesis, signaling or response in the nucleus. The role of ethylene in cold-stress anthers is not known although some reports lead to the conclusion that a proper balance of its levels must be achieved for programed cell death during pollen development (Takada et al. 2006;Wei et al. 2013). The expression of a number of ethylene-related TFs is upregulated after cold treatment in ERY anthers. Only three of them (Sub1C, OsERF1 and OsERF3) are well-known TFs involved in ethylene metabolism or responsive to this phytohormone (Fukao et al. 2006;Hu et al. 2008;Wan et al. 2011). However, no evidence showing a participation of these TFs in cold stress responses has been shown yet. This work also shows more than ten putative TFs with increased expression after cold treatment in ERY, some of them with high relative expression in anthers.

Conclusions
Our results suggest that active chloroplasts may be important during anther development to face cold episodes, making the photosynthetic machinery in anthers a target to improve cold tolerance in susceptible varieties. Our study also suggests that anther-specific ethylene metabolism, through the many transcription factors responding to cold stress in ERY, might be a good candidate for cold tolerance improvement in agronomical relevant rice cultivars as part of breeding programs or novel tools emerging for a rationale design of rice cultivars. Our studies provide an important step in the identification of possible pathways controlling cold tolerance but will of course need further validation by molecular and genetic analyses.

Plant material and growing conditions
Oryza sativa subspecies with contrasting tolerance to cold stress (Erythroceros Hokkaido, Lijiangheigu, Sant'Andrea, Carnaroli, Arborio, Volano and Nipponbare) were used in the study. Plants were grown in a temperaturecontrolled growth chambers maintained at 26/22°C day/night temperature and day/night relative humidity (RH) of 75-85% under LD (16/8 h light/dark) or heading-inductive SD conditions (10/14 h light/dark). Light intensity was set to ∼500 μmol/m 2 /s. Seeds were pre-germinated in seeding trays for 3-4 weeks and then seedlings were transplanted individually into pots. After 7-10 weeks, depending on the cultivar, plants were transferred to SD inductive conditions for 3 more weeks. Cold-stress treatments (8°C) were applied only during the night period for 1, 3 or 5 days when rice plants showed particular auricle distances between flag leaf and penultimate leaf (AD) as described in results.

Spikelet fertility
Between 2 and 13 tillers from at least three plants of each cultivar were used to estimate spikelet fertility for both control and stress treatments. At physiological maturity, the total number of spikelets and the number of grains formed were recorded and used to determine spikelet fertility. Each experiment was replicated between 2 and 4 times and mean ± SEM was obtained.

Light microscopy
Anther developmental stages were determined by histological analysis and light microscopy. Anther samples were fixed in 3.7% formaldehyde, 5% acetic acid and 50% ethanol for 30 min with vacuum, and therefore dehydrated in ethanol series. Infiltration and embedding were carried out with Technovit acrylic resin according to the manufacturer's instructions. 5-μm-thick sections were obtained with microtome, mounted in slides and stained with 0.1% toluidine blue. Images were captured with a Leica DM LB microscope coupled to a DFC280 camera.

Sample collection
Between 5 and 8 plants for each biological replica were used for sample collection. For cold stress or control treatments, plants were transferred to the growth chambers when the AD matched post-meiotic developmental stage for each cultivar, corresponding to stage 10 according to Zhang et al. (2011). Spikelets from the middle part of the panicles were collected immediately after one night of corresponding treatment (when the lights of the growth chamber switched on in the morning), dissected to separate the anthers and collected in tubes suspended in liquid nitrogen. Frozen anther samples were stored at − 80°C until further use.

Total RNA extraction and expression analysis
Total RNA was extracted from pooled anthers by using TRIzol® (Invitrogen). Three biological replicates were collected for each treatment (Cold-stress and Control). The quality and quantity of RNA samples were assessed by gel electrophoresis and Nanodrop quantification.
cDNA synthesis and quantitative RT-qPCR were performed as described previously (Gonzalez-Schain et al. 2016). Briefly, standard RT-qPCR was carried out to test efficiency of the primers (listed in Additional file 11: Table S7) and therefore large-scale RT-qPCR was used to study the expression of deregulated genes from anther samples by using Microfluidic Dynamic Array developed by Fluidigm Corporation (Spurgeon et al., 2008). Three biological replicates with three technical replicates were performed for each sample. Data were normalized using OsUBQ (LOC_Os02g06640) and NABP (LOC_Os06g11170) genes as reference.

Illumina sequencing
Upon treatment with TURBO DNAse I (AMBION), 4 μg of RNA from each sample were used to produce sequencing libraries with the TruSeq mRNA sample preparation kit (Illumina). Sequencing of poly(A) RNA samples was carried out in multiplex (6 samples per lane, single 50 bp reads) with Illumina Hi-seq 2000 platform at IGA Technology Services (Udine, Italy). Quality control of the raw sequence data was done using FastQC (Babraham Bioinformatics).
Mapping of short reads, quality analysis and assessment of gene expression analysis for RNA-seq Evaluation and processing of raw data was performed on the commercially available CLC Genomics Workbench v.4.7.1 as described previously (Gonzalez-Schain et al. 2016). Gene expression values were based on reads per kilobase of exon model per million mapped read (RPKM) values. Fold change (FC) and log 2 FC was calculated in terms of RPKM of the corresponding transcripts. To obtain statistical confirmation of the differences in gene expression, P and FDR values were computed using Baggerley's test on expression proportions. We applied a threshold value of P < 0.05 and FDR < 0.05 to ensure that differential gene expression was maintained at a significant level (5%) for the individual statistical tests. Absolute FC ≥ 2 was set as threshold limit to obtain the differentially expressed genes. Raw and processed data have been deposited in NCBI's Gene Expression Omnibus (Edgar et al. 2002) and are accessible through GEO Series accession number GSE137002 (https://www.ncbi.nlm.nih. gov/geo/query/acc.cgi?acc=GSE137002).
To gain insight into the biological processes associated with the regulated genes, we determined which GO annotation terms were over-represented, in both upand down-regulated gene lists. Gene set enrichment analysis was performed with the agriGO database (Du et al. 2010) using the Singular Enrichment Analysis (SEA). Mapman Software (Thimm et al., 2004) was also used to get an overview of deregulated genes falling in different hierarchical categories (or bins), applying Benjamini-Hochberg correction as statistical analyses provided in the software with p < 0.05.