Nitrate application or P deficiency induce a decline in Medicago truncatula N2-fixation by similar changes in the nodule transcriptome

Nitrogen fixation of Medicago truncatula is regulated by the nitrogen status of leaves through inducing a repeatedly occurring 24-h nodule activity rhythm that reduces per day nitrogen fixation. The hypotheses of the present study were that (1) long-term moderate whole-plant P deficiency in Medicago truncatula induces an according daily rhythm in nitrogenase activity comparable to that induced by nitrate application and (2), the changes in the nodule transcriptome that go along with a strong nitrogenase activity decline during the afternoon would be similar under P deficiency or after nitrate supply. The nodules of plants in a low P treatment developed a rhythmic pattern of activity that resembled the pattern following nitrate application. A comprehensive, RNAseq-based comparative transcriptome profiling of nodules during a repeated part of the rhythm revealed similarities between P deficiency versus nitrate supply. Under both treatments, the formation of nitrogenase was targeted by a reduction in the expression of genes for nodule-specific cysteine-rich peptides (NCR), and possibly also by a disturbance of the inner cell iron allocation. A strong reduction in the expression of leghemoglobin is likely to have restricted the supply of oxygen for respiration.

growing shoots have an increasing N need 6,7 . In contrast, there is less understanding of the molecular processes by which existing nodules are downregulated when the N need decreases. Numerous studies have shown that nodules lose activity quickly (in hours to days) when alternative N sources became available or growth slows (caused for example by emerging nutrient deficiency) [8][9][10] . The basic principle appears to be a systemic N feedback regulation 11,12 . This means that a loss of nodule activity occurs when the plant leaves encounter available N beyond their requirement for the current or possible growth rate. This situation can be induced by various and diverse factors (emerging external or internal [senescing leaves] alternative N sources, restriction of growth). Under such conditions the nodules rather than showing a linear decline develop a 24-h activity rhythm that reduces the per day nitrogen fixation 13 . The continuous measurement of nitrogen fixation revealed that the daily pattern is apparently the primary reaction of nodule activity to various factors that have in common that the leaf N need is satisfied 13 . One consistently occurring part of the 24-h activity rhythm was what might be called an 'afternoon decline' . A strong decline in activity began at around noon and abated and turned into an increase in activity at around 5 pm. One pronounced downregulation of nitrogen fixation activity of existing nodules is known to occur when nitrate becomes available 14,15 . The molecular basis for this phenomenon has been shown to consist of various molecular effects that rather than primarily restricting some external supply of carbon or oxygen, initially target nitrogenase formation and respiration on the transcriptional level 16 .
Based on the striking similarity of the pattern in the 24-h rhythm of nodule activity between nitrate application and other treatments that reduce nitrogenase activity 10,13 , the hypotheses of the present study were that (1) long-term moderate whole-plant P deficiency in Medicago truncatula induces an according daily rhythm in nitrogenase activity comparable to that induced by nitrate application and (2), the molecular reason for a strong nitrogenase activity decline during the afternoon -a major part of the daily rhythm -would be similar under P deficiency or after nitrate supply.
To validate these hypotheses, a procedure was established to grow M. truncatula plants in nutrient solution with limited P supply to produce lower growth but keep the plants alive. The growth system was embedded in a measurement system that allowed continuous, non-invasive monitoring of nitrogenase activity by measuring nodule H 2 evolution 13 . In a further step the molecular reasons for the 'afternoon decline' induced by nitrate application or P deficiency should be revealed through a comprehensive RNAseq-based transcriptome profiling of nodules during that decline.

Plant growth and P supply. Plants of M. truncatula previously inoculated with Sinorhizobium meliloti
(102F51) and with visible active nodules were grown in a quasi flow-through nutrient solution system, where the P concentration was adjusted daily. Using this system, a daily adjustment of the P concentration to 5 μ M proved optimal for growth of M. truncatula when depending solely on nitrogen fixation for N supply. A daily measurement of the P concentration before the adjustment to 5 μ M showed that the concentration was never depleted below 3 μ M P during the experiment (eight weeks of growth), thus allowing the plants a continuous optimal P uptake and N 2 -fixation. As a treatment with limited P supply (P deficiency), a daily adjustment of the nutrient solution to 1 μ M (20% of full P supply) was chosen. The daily P concentration before the adjustment of the solution showed that the plants reduced the concentration to a minimum value below which the plants cannot take up P anymore 17 (C Lmin value [0.2-0.3 μ M]) throughout the experiment and thus suffered a limited P supply. The preliminary experiment had shown that M. truncatula plants under these conditions experienced reduced growth, but maintained a certain level of nitrogen fixation and did not die during the planned eight-week experimental period. Dry matter (DM) formation, nodule number and N and P concentrations in the various organs of the plants are shown in Fig. 1. The plants of the fully nourished treatment showed vigorous growth in these experiments and developed more than 10 g DM during the eight-week experimental period (Fig. 1A). The DM of shoots, nodules and the shoot/root ratio were significantly reduced by the low P treatment. Active nodules were significantly lower in P-deficient plants and the number of senescent nodules was not different for both treatments (Fig. 1B). N concentration in shoots and roots was significantly higher in the low P treatment (Fig. 1C) and the nodule P concentration in the low P treatment was lower when compared to the nodules of the control (Fig. 1D). The relative reduction in shoot P concentration was stronger than that in nodules. The overall picture with respect to these parameters was therefore in accordance with moderate P deficiency (slower growth, reduced shoot/root ratio and increased N concentration in leaves [18][19][20][21] ).
Nitrogen fixation under low P. Nitrogen fixation was measured continuously during one day in the eighth week of growth. While the plants in the control showed a stable nitrogen fixation throughout a 24-h period, the low-P plants had developed the typical daily rhythm. In particular, a strong decline occurred during the time between about noon and 5 pm (Fig. 2). Nitrogen fixation efficiency was reduced under low-P conditions. A lower relative efficiency (electron allocation efficiency [EAC]) compared to the control plants ( Fig. 2) contributed to that. The low EAC did not change when measured at noon compared to measurement after eight hours of darkness during the night. Furthermore, the nodule-specific activity, expressed as μ g N fixed per day and mg nodule DM (active nodules), was about 35% lower when compared to the control plants (Table 1). Active nodules were distinguished from inactive nodules by a pinkish versus a greenish to brown colour of the nodules. Although we are not aware of any direct measurement that greenish nodules have totally lost their nitrogenase activity, the fact that the green colour results from a decay product of leghemoglobins, which are of vital importance in nodule functioning, makes this highly probable. In any event, when considering the entirety of the collected nodules, the gap between the control and low-P treatment would even have widened since the share of inactive nodules was higher in low-P plants. Plants were grown with a daily adjustment of the nutrient solution to either 5 μ M or 1 μ M P in the sufficient (+ P) and limited (− P) P-treatment respectively. The P concentration in the + P treatment was never reduced below 3 μ M by the plants during the experiment, while the concentration in the limited P treatment reached the C Lmin value every day. Consequently while the plants in the + P treatment enjoyed a continuous P supply, the plants in the -P treatments suffered limited supply. The leaves of the -P plants partially developed a reddish colour, known to be the result of anthocyanin formation 48 , a symptom of P deficiency 49 . At the same time, a particular deep green colour of the remaining parts of the leaves was indicative of the higher leaf N concentration in the -P treatment. Data are means of six replicates. The two treatments were compared by the t-test. Statistical differences in comparison to the + P treatment are indicated by stars (*) above the bars (*different with a probability of 95%, **different with a probability of 99%, ***different with a probability of 99.9%; n = 6; n.s. = not significantly different). showed no shifts within a treatment over the time period, but a lower EAC in the P deficiency treatment when compared to the control and the nitrate application treatment. Data were taken every 5 min as means of six replicates. Data for control and nitrate were taken from 16  Gene expression after nitrate application (Mt3.5v3 versus Mt4.0v1). The study was supported by the fact that the harvest of P-limited nodules for RNAseq analysis was performed at the exact date and time as that of the nitrate-treated plants and the control 16 . The RNAseq experiment was performed before the study's authors became aware of the importance of the repeated 24-h rhythm in nodule activity under treatments that reduce nitrogen fixation. The original reason for combining both experiments was to be able to use one control for both treatments (nitrate application and long-term P deficiency). However, the fact that both treatments were grown in the same experiment and harvested at the same age and time of day also allowed a valid comparison of the transcriptome of nodules from nitrate-treated plants versus those from plants grown under limited P supply. The data of the comparison nitrate versus control are published in Cabeza et al. 16 (www.plantphysiol. org, Copyright American Society of Plant Biologists), but were recalculated on the basis of the newly released, significantly enhanced cds annotation for M. truncatula (Mt4.0v1 versus Mt3.5v3 that was previously used). The reason for that approach was to form a basis for comparison of the molecular reasons for the 'afternoon decline' induced by long-term P deficiency versus nitrate application. The harvest in both treatments was done at the exact same time, at which in both cases the 'afternoon decline' was in full effect (Fig. 2).
The main described effects of nitrate on the nodule transcriptome 16 were confirmed and reinforced (for example by the number of downregulated Transcription Units [TU] for nodule-specific cysteine-rich [NCR] peptides) by revisiting the data with the new cds annotation (Supplementary Dataset File 1: sheet 2). A strong downregulation of numerous genes for NCR-peptides was confirmed along with the downregulation of all expressed genes for leghemoglobin. Furthermore, a putative effect on nodule iron turnover was confirmed. In particular, a strikingly strong downregulation of a gene for a 'nicotianamine synthase-like protein' (Medtr1g084050) (log 2 FC = − 5.0) was again found with the new cds annotation. The importance of the gene for efficient nitrogen fixation and for nodule acclimation to altered external oxygen pressure has recently been shown 22 . Gene expression in nodules from plants with moderate P deficiency. Nodules for the RNAseq analysis were harvested during the P deficiency induced 'afternoon decline' parallel with the nodules harvested 4 h after nitrate application (see Fig. 2). The sequencing yielded 10 to 15 million mapping reads per sample and thus a sufficient depth for the objective of our analysis. About 80% of the reads hit at least one time in the cds annotation Mt 4.0v1. About 33,000 TUs of the cds annotation were hit at least one time. When setting a threshold of > 20 DESeq normalised counts in either control, treatment or both (see Supplementary Dataset File 2: sheet 1) 16,725 TUs were found to be expressed. This corresponded to 26.8% of all TUs annotated in the cds annotation Mt4.0v1. A cluster analysis (Fig. 3) shows that within the control and P-deficiency treatment, the replicates resembled one another, while there was a clear difference between the treatment and the control. A DESeq comparison of the control and the treatment revealed 1,738 differentially abundant TUs (False Discovery Rate [FDR] < 0.01; n = 3) which represent a 10.4% of the expressed TUs (Supplementary Dataset File 2: sheet 2). Of these TUs 900 were differentially expressed with a log 2 FC greater than [1], meaning that the abundance of the mRNA was at least doubled or halved. A quantitative reverse transcription-polymerase chain reaction (qRT-PCR) validation proved the reliability of the RNAseq data (Supplementary Dataset File 2: sheet 5). Twelve genes were tested using the same RNA pools previously used for the next-generation sequencing. The qRT-PCR results were significantly correlated to the RNAseq data (r 2 = 0.89). The slope of the regression line was close to one (0.94).

Discussion
The low-P treatment resulted in P-deficient plants that could be kept alive during eight weeks of growth in the present study's experimental system. The plants showed typical growth reactions to limited P supply 21,23 . A higher N concentration in the shoots supports the view that plant growth under these conditions is not limited by nitrogen fixation but rather restricts nodule activity, possibly by a shoot-related N-feedback effect 20,24 . The low-P treatment in this experiment resulted in N 2 -fixation that clearly showed a 24-h rhythm of nitrogenase activity. The pattern of the rhythm was consistent with that under other treatments 13 . The relative, specific and per plant activity of the nodules under P deficiency was significantly lower when compared with the control plants. There are scattered reports that nodule-specific activity remains unchanged or is even increased under low P 20,21 . This was clearly not the case under the present study's conditions, even when only undoubtedly active nodules (pink colour) were considered. Overall, the specific activity might be maintained depending on the age of the plants and the extent of the P deficiency 19,25 , but is generally lower under whole-plant P deficiency. In addition, strong variations in measurements of specific activity might have been the result of point measurements at different times during the 24-h activity cycle. For example, under our conditions a point measurement in the P-deficient treatment at 5 pm would have yielded a value about one third lower when compared to a corresponding measurement at noon. Comprehensive nodule transcriptome profiling based on RNAseq analysis revealed that during the 'afternoon decline' of the 24-h activity pattern, similar molecular mechanisms appeared to be active in nodules of treatments as diverse as low-P supply versus (over)-sufficient nitrate application (Supplementary Dataset Files 3 and 4). In a sense, P deficiency and external N supply (nitrate) can be considered as the most distinct treatments that reduce nitrogenase activity. While one is promoting plant growth (nitrate) and is known to reduce nodule oxygen permeability, moderate P deficiency reduces plant growth and increases nodule oxygen permeability. Nevertheless, similar molecular mechanisms are active during the 'afternoon decline' of nodule activity. This was first indicated by the fact that while both treatments showed strong shifts in the transcriptome when compared to the common control (low-P or high nitrate, 1,738 or 2,257 TUs respectively; [FDR < 0.01, n = 3]), a direct comparison between the transcriptome of the nodules in both treatments revealed only 663 TUs to be of different abundance (FDR < 0.01, n = 3) (Supplementary Dataset File 3: sheet 2). That means that many of the changes in both treatments in comparison to the control went in the same direction. Among those still different 663 TUs in the direct comparison of P deficiency and nitrate application were 86 cases in which significant changes also went in the same direction compared to the control. However, a different extent of the change left significantly different TU abundances when the two treatments were compared. Of the remaining TUs, numerous are originated from genes with a known nitrate-specific regulation (nitrate transporters, nitrate and nitrite-reductase). When the lists of differentially regulated TUs in the two diverse treatments versus their common control were compared, an intersection of 1,074 TUs was found (Fig. 4). This corresponded to 48 or 62% of all differentially regulated TUs in the nitrate-treated or P-deficient nodules respectively. A closer look at the log 2 FC rates revealed a striking similarity between the reaction of the regulated TUs in the two diverse treatments, with a slightly stronger effect overall of the nitrate (log 2 FC average for all differentially expressed TUs either up or down regulated: nitrate + 1.38, − 1.48; P-def. + 1.35, − 1.25) (Supplementary Dataset File 4: sheet 1). The log 2 FC numbers showed a close correlation (r 2 = 0.89) and the predicted regression line was close to a 1:1 correlation (Fig. 5).
The commonly regulated TUs of both treatments reflected all the effects described for nitrate 16 . When the log 2 FC of significantly and commonly regulated TUs of both treatments were weighted together by the RPKM (Reads Per Kilobase per Million mapped reads) values of the control (RPKM control × log 2 FC after 4-h nitrate impact × log 2 FC after eight weeks of low P supply) the above-mentioned Medtr1g084050.1, annotated as 'nicotianamine synthase-like protein' , was the most strongly affected TU (Supplementary Dataset File 4: sheet 2). Iron is of central importance for nodule functioning because it is an essential component of nitrogenase, leghemoglobin and ferredoxin 26 . Nicotianamine synthase (NAS) is of pivotal importance for the inner cell and inner plant allocation and re-allocation of iron and other metal ions 27,28 . Among the TUs annotated as NAS, Medtr1g084050.1 was the by far most strongly abundant in nodules (Supplementary Dataset File 4: sheet 3). It has recently been shown, partially on the basis of a Tnt1 mutant, that this gene is important for nodule efficiency and the necessary neoformation of nitrogenase after the impact of elevated oxygen on nodules 22 . Taken together, it appears possible that the iron supply of bacteroids and infected cells might function as a regulatory mechanism for nitrogenase formation. There were a total of 203 TUs for NCR-peptides that were downregulated under both treatments from a high level of expression in the nodules. In fact, when the log 2 FC in both treatments were weighted by the RPKM values of the control, TUs for NCR peptides made up 56 of the 100 most differentially regulated TUs. NCR peptides are formed in the nodules of legumes of the inverted repeat-lacking clade of legumes and targeted to bacteroids 29 . They are closely related to defensins and attach to the bacteria to keep them in the bacteroid stage. This stage basically means that the bacteria cease to divide and strongly reduce metabolic activity except that they form nitrogenase vastly in excess of their own needs 30 . The delivery of these peptides appears to be the way in which the higher plant forces (sanctions) the microsymbiont to return sufficient ammonia in exchange for energy and nutrients 31,32 . During the 'afternoon decline' under both treatments, this pressure was strongly reduced. Thus a concerted downregulation of the expression of numerous genes for these peptides appears to be among the initial mechanisms that induce the 'afternoon decline' in nodule activity. Recent reports show that even the loss of function of a single gene for an NCR peptide results in reduced nitrogen fixation. The mutants dnf4 (defective in nitrogen fixation) and dnf7 are impaired in the functioning of the gene Medtr4g035705 and Medtr7g029760 respectively. Both genes encode an NCR-peptide 33,34 and were expressed strongly in the nodules in the present  and significantly downregulated by both treatments (log 2 FC Medtr4g035705.1 − 0.8 or − 0.59 under nitrate or P deficiency respectively; log 2 FC Medtr7g029760.1 − 1.42 or − 0.9 under nitrate or P deficiency respectively). Nevertheless, the expression and also the weighted log 2 FCs in transcript abundance of both TUs for the two NCR-peptides were only to a medium extent when considering the numerous other affected TUs for NCR-peptides. This fact might support the specific individual importance of the different NCR-peptides. On the other hand, in an earlier study, the individual knockdown or overexpression of five different genes for NCR-peptides resulted in no detectable phenotype, suggesting a certain redundancy among the NCR-peptides 35 . Our data suggest that the regulation of the NCR-peptides genes occurs in a concerted, one might say 'module-like'-fashion. The targeting of the peptides to the symbiosome is ensured through a signal peptidase complex. A mutation in the gene Medtr3g027890 (dnf1), which encodes for a core unit of this complex, results in fix − nodules 36 . Expression of the gene was not affected during the 'afternoon decline' , however it was strongly reduced when the nitrate impact had lasted for 28 h and the nodules had lost almost all activity and showed molecular indications for senescence 16 .
A further clear common effect of the treatments was the concerted downregulation of almost all the expressed genes for TUs of leghemoglobin. Eleven of the 16 genes annotated as leghemoglobins in the cds annotation Mt4.0 v1 were expressed and nine of them were strongly downregulated by both treatments in the nodules in the present experiments (Supplementary Dataset File 4: sheet 3). The expression had reached a high level in the control. This could be illustrated by the fact that five of these TUs were among the 100 most strongly abundant TUs in the nodules, and is supported by leghemoglobin protein accumulating to mM concentrations in the cytoplasm of nodule cells 37 . Functionally, leghemoglobin provides a substantial buffer capacity for free oxygen in the nodule interior, contributing to the maintenance of free oxygen concentrations at a low nM level. Furthermore, the leghemoglobin protein ensures a constant flow of oxygen for the high-affinity bacterial cytochrome oxidases to maintain sufficient respiration. The two imposed treatments are known to have a different effect on nodule O 2 conductance. While nitrate has been shown to lead to lower O 2 influx into the nodule and lower partial oxygenation of leghemoglobin 38 , nodule O 2 conductance increases under P deficiency 25,39 . The strong downregulation of leghemoglobin TUs under both treatments consequently suggests that O 2 buffering is not the primary function of leghemoglobin. Instead the lower expression of leghemoglobins might be a mechanism for lowering nodule and bacteroid respiration and thus ATP production for driving nitrogenase. RNA i interference with a gene for leghemoglobin in Lotus japonicus resulted in a fix − phenotype 40 . The loss of function of the leghemoglobin gene was accompanied by an almost total loss of nitrogenase protein in the nodules. This finding is supported by the comparative transcriptome analyses data so far obtained in our lab, in that the concerted expression of all expressed leghemoglobin genes always accompanies the nitrogenase activity of the nodules and the expression of genes for NCR-peptides 10,16,22 . These observations support the view that leghemoglobins might have functions beyond O 2 binding and buffering 41 , possibly in oxygen concentration signalling and regulating nodule activity.
In conclusion, the data from the present study supported the fact that the downregulation of nitrogenase under conditions of over-sufficient whole-plant N availability (induced by nitrate application or P deficiency) initially occurs through a 24-h rhythmic pattern of nodule activity. This pattern is induced by common molecular mechanisms in the nodule directly involved in nitrogenase formation and respiration.

Materials and Methods
Plant growth. Seeds of Medicago truncatula (Gaertn.) cv. Jemalong A17 were submerged in H 2 SO 4 (96%) for 5 min for chemical scarification, sterilised with 5% (v/v) sodium hypochlorite for 3 min and rinsed several times with deionised water. The seeds were subsequently kept at 4 °C for 12 h in darkness, submerged in tap water. The next step was a two-to-four day slight shaking of the submerged seeds at 25 °C and continuous light. When the seeds had developed a ca. 20-mm-long primary root, they were transferred to small growth boxes (170 mm × 125 mm × 50 mm) filled with aerated nutrient solution. Each box contained 20 seedlings, which were fixed by small x-shaped cuts in the adhesive tape on the upper side of the boxes.
The plants were grown for two weeks in these boxes in a growth chamber with a 16/8-hour light/dark cycle at 25/18 °C respectively. The nutrient solution level in the boxes was maintained by the addition of an appropriate amount of nutrient solution every other day. The light intensity at plant height was approximately 400 μ mol m −2 s −1 . Immediately after being transferred to the growth boxes, the seedlings were inoculated with 1 mL box −1 of a stationary Sinorhizobium meliloti (Sm) (102F51) YEM-culture of an approximate cell density of 10 9 mL −1 . The Sm strain induced good nodulation with the first macroscopic nodules visible after about 7-10 days. This strain does not contain an uptake hydrogenase 42 .
After two weeks, the plants were transferred to glass tubes, allowing the separate measurement of root/nodule H 2 and CO 2 evolution. The system is described by Fischinger and Schulze 43 and Fischinger et al. 44 . The setup was enhanced by connecting a group of six plants through the lower side of the glass tubes to a 20-L nutrient solution container. Thus, using gravity, the nutrient solution level in the glass tubes depended on the height of the container, and losses via plant transpiration could be adjusted by the addition of nutrient solution to the container, thereby not interfering with the measurements in the root/nodule compartment. In addition, a pump in the container drove a nutrient solution flow of about 10 mL min −1 into the upper part of each individual glass tube 13  0.5 mM NH 4 + concentration by the addition of (NH 4 ) 2 SO 4 since low concentrations of ammonium support growth and nodule formation in young M. truncatula plants 45 . During their growth in the glass tubes, the plants depended solely on nitrogen fixation for their N nutrition. Phosphorus (P) was added daily as KH 2 PO 4 to the nutrient solution to adjust a concentration of 5 μ M P for plants in the control and the nitrate experiment. The nutrient solution in the low-P treatment was adjusted daily to 1 μ M P.
The nutrient solution was changed every week. During this procedure, the pump in the container was switched off and the backflow from the glass tubes to the container was blocked. Using this method, the ongoing measurements in the root/nodule compartment were not affected. The procedure of replacing the nutrient solution in the container took about ten minutes, after which the nutrient solution turnover system was set back to normal. Dry matter, P and N concentration in plant material. Plants were dried to constant weight at 65 °C before the dry matter was determined. Subsamples of plant dry matter were digested in concentrated HNO 3 at 180 °C and the P concentration in the digest was measured colorimetrically by the molybdenum-vanadate method 46 . N concentration was determined by means of an elementary analyser (Vario EL, Elementar Analysen GmbH, Hanau, Germany) 47 .
Root/nodule gas exchange measurement. The system for measuring nodule H 2 evolution, including the determination of apparent nitrogenase activity (ANA), total nitrogenase activity (TNA) and the calculation of the electron allocation coefficient (EAC), is described by Fischinger and Schulze 47 . ANA was measured continuously over a two-day period. The H 2 concentration data in the continuous gas stream were taken every 5 min. TNA was measured on parallel plants at various points in time during the experiment.

Quantitative reverse transcription-polymerase chain reaction (qRT-PCR).
For the quantitative reverse transcription-polymerase chain reaction (qRT-PCR), gene-specific primers were designed using Primer-Blast (http://www.ncbi.nlm.nih.gov/tools/primer-blast). To validate the RNAseq results, RNA from the RNA pools previously used for Illumina sequencing was used. The samples and reference genes were run in triplicate using the Fast SYBR Green Master Mix protocol (Applied Biosystems, Darmstadt, Germany) on a StepOne ™ Real-Time PCR System (Applied Biosystems, Darmstadt, Germany) following the manufacturer's recommendations.
qRT-PCR was performed in a 20-μ L reaction mix containing 10 μ L Fast SYBR Green Master Mix, 4 μ L ddH 2 O, 2 μ L forward primer (2 pmol/μ L), 2 μ L reverse primer (2 pmol/μ L) and 2 μ L template cDNA (100ng). The PCR conditions according to the Fast SYBR Green Master Mix protocol were 20 seconds of pre-denaturation at 95 °C, 45 cycles of 3 seconds at 95 °C and 30 sec at 60 °C, followed by steps for the dissociation curve generation (15 sec at 95 °C, 60 sec at 60 °C and a stepwise increase of 0.3 °C up to 95 °C). StepOnePlus software (Applied Biosystems, CA, USA) was used for data collection. Relative transcript levels were obtained using the comparative C T method (Δ Δ C T ).
Library preparation for RNAseq was performed using the TruSeq RNA Sample Preparation Kit (Illumina, Cat. N°RS-122-2002) starting from 500 ng of total RNA. Accurate quantitation of cDNA libraries was performed using the QuantiFluor ™ dsDNA System (Promega, Mannheim, Germany). The size range of the final cDNA libraries was determined by applying the DNA 1000 chip on the Bioanalyzer 2100 from Agilent (Böblingen, Germany) (280 bp). The cDNA libraries were amplified and sequenced using the cBot and HiSeq2000 from Illumina (SR; 1 × 50 bp; 6 GB ca. 30-35 million reads per sample).
The sequence images were transformed with the Illumina software BaseCaller to bcl-files, which were demultiplexed to fastq files with CASAVA v1.8.2. Quality check of raw data was performed with the FastQC Software (v.0.10.0, Babraham Bioinformatics, Cambridge, UK), in order to assess the overall and per-base quality of reads. Illumina adapters were removed with the Trimmomatic Software 48 . The sequences were aligned by Bowtie2 (v2.0.2, Johns Hopkins University, Baltimore, Maryland, USA) to the J. Craig Venter Institute (JCVI) cds annotation Mt4.0v1.
Gene expression and statistical analyses. The expression level of TUs in each library was calculated by quantifying the number of Illumina reads that mapped to the Mt4.0v1 cds annotation using the Bowtie program, counting only 'unique hits' . The term 'gene expression' or 'expression level of the TU' is used, although what is measured by RNAseq is a net result of expression and decay intensity of mRNA. The raw counts were normalised using the RPKM method 49 . RPKM values were used to compare the abundance of TUs within one treatment. Transcript units showing differential abundance were identified using the DESeq method for pair-wise differential expression analysis since this has proven to be the most reliable method for RNAseq based differential expression analysis 50 . Differentially abundant transcript units identified by DESeq were required to have an FDR value < 0.01. Additionally, either the control or treatment tissues had to have a DESeq value above 20 before a gene was considered as having been expressed. Figures 1, 2 and 5 were generated using GraphPad Prism 6. Statistical differences among treatments were calculated using a t-test. Cluster analysis for the dendrogram (Fig. 3) was calculated with the R package "hclust" through a complete-linkage clustering. The Venn diagram was generated for the set of TUs differentially expressed in both treatments at http://bioinformatics.psb.ugent. be/webtools/Venn/.