# Poplar stem transcriptome is massively remodelled in response to single or repeated mechanical stimuli

Lise Pomies, Mélanie Decourteix, Jérôme Franchel, Bruno Moulia, Nathalie Leblanc-Fournier

► To cite this version:

Lise Pomies, Mélanie Decourteix, Jérôme Franchel, Bruno Moulia, Nathalie Leblanc-Fournier. Poplar stem transcriptome is massively remodelled in response to single or repeated mechanical stimuli. BMC Genomics, 2017, 18, 16 p. 10.1186/s12864-017-3670-1. hal-01528849

HAL Id: hal-01528849

https://hal.science/hal-01528849

Submitted on 26 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
RESEARCH ARTICLE

Open Access

<figure><img src="image_2.png" /></figure>

CrossMark

# Poplar stem transcriptome is massively remodelled in response to single or repeated mechanical stimuli

Lise Pomiès†, Mélanie Decourteix†, Jérôme Franchel, Bruno Moulia and Nathalie Leblanc-Fournier*

**Abstract**

**Background:** Trees experience mechanical stimuli -like wind- that trigger thigmomorphogenetic syndrome, leading to modifications of plant growth and wood quality. This syndrome affects tree productivity but is also believed to improve tree acclimation to chronic wind. Wind is particularly challenging for trees, because of their stature and perenniality. Climate change forecasts are predicting that the occurrence of high wind will worsen, making it increasingly vital to understand the mechanisms regulating thigmomorphogenesis, especially in perennial plants. By extension, this also implies factoring in the recurring nature of wind episodes. However, data on the molecular processes underpinning mechanoperation and transduction of mechanical signals, and their dynamics, are still dramatically lacking in trees.

**Results:** Here we performed a genome-wide and time-series analysis of poplar transcriptional responsiveness to transitory and recurring controlled stem bending, mimicking wind. The study revealed that 6% of the poplar genome is differentially expressed after a single transient bending. The combination of clustering, Gene Ontology categorization and time-series expression approaches revealed the diversity of gene expression patterns and biological processes affected by stem bending. Short-term transcriptomic responses entailed a rapid stimulation of plant defence and abiotic stress signalling pathways, including ethylene and jasmonic acid signalling but also photosynthesis process regulation. Late transcriptomic responses affected genes involved in cell wall organization and/or wood development. An analysis of the molecular impact of recurring bending found that the vast majority (96%) of the genes differentially expressed after a first bending presented reduced or even net-zero amplitude regulation after the second exposure to bending.

**Conclusion:** This study constitutes the first dynamic characterization of the molecular processes affected by single or repeated stem bending in poplar. Moreover, the global attenuation of the transcriptional responses, observed from as early as after a second bending, indicates the existence of a mechanism governing a fine tuning of plant responsiveness. This points toward several mechanistic pathways that can now be targeted to elucidate the complex dynamics of wind acclimation.

**Keywords:** Mechanotransduction, Time series, Acclimation, Accommodation, microarray, Thigmomorphogenesis, Abiotic stress, Mechanical stimuli

<footnote>* Correspondence: nathalie.leblanc@uca.fr</footnote>

<footnote>†Equal contributors</footnote>

<footnote>Université Clermont Auvergne, INRA, PIAF, F-63000 Clermont-Ferrand, France</footnote>
## Background

In their fluctuating environment, plants are constantly exposed to abiotic stimuli to which they are sensitive and responsive. Some of these stimuli, such as wind exposure, count a strong mechanical component that has a major influence on plant growth and development. Exposure to such mechanical stimulations results in the so-called thigmomorphogenetic syndrome characterized by a reduction in stem elongation [1, 2], local stimulation of radial growth [3, 4], and modification of the stem's mechanical properties [5]. These alterations of plant architecture are thought to improve plant acclimation to chronic wind regimes [6]. Using continuous monitoring techniques [4, 7], the early kinetics of these plant growth responses were studied in tomato and poplar by applying quantified stem bending. In poplar, a single transitory stem bending led first to a short period of secondary growth inhibition (4 h) followed by a massive stimulation of the growth rate over 3 days, and finally a relaxation to normal values [4].

To better understand tree acclimation to wind, it is important to first unravel the processes that regulate thigmomorphogenetic syndrome. However, the way plant cells perceive and transduce mechanical signals is still poorly understood. Two major classes of potential mechanosensors are thought to be involved: plant MechanoSensitive (MS) ion channels and Receptor-Like Kinases (RLK) inserted into the cell wall-plasma membrane-cytoskeleton continuum (see [8] for review). Prior to the advent of transcriptomics, studies had identified a handful of mechanoresponsive genes, including *TOUCH* genes (*TCH*) that mainly encode calmodulins or calmodulin-like proteins and Xyloglucan endo-Transglycosylase/ Hydrolase (XTH) [9], genes encoding protein kinases [10, 11], Transcription Factors (TF) [12], genes involved in Jasmonic Acid (JA) and ethylene synthesis [13, 14], and genes involved in antioxidative responses [15]. In 2005, a transcriptome analysis of touch-stimulated *Arabidopsis* rosette leaves allowed a more global insight into the molecular functions altered after a touch-stimulus. Over 700 genes presented regulated expression 30 min after the stimulus [16], representing over 2.5% of the genome. Among these genes, the vast majority (589/760) were up-regulated. Analysis of the functional categorization of these up-regulated genes revealed enrichment in genes encoding calcium-binding proteins, cell-wall proteins, disease resistance proteins, kinases and TF, and a decline in genes involved in general metabolism and the ubiquitin/protein degradation pathway. For the down-regulated genes, the “transcription factor” and “cell-wall-associated protein” categories were over-represented. In a subsequent transcriptomic study of the effect of a long-term exposure (8 weeks) to low-speed wind in *Populus nigra* leaves, Fluch et al. found at least 98 up- and 94 down-regulated genes [17], including genes encoding cell-wall modification

proteins, proteins with regulatory roles (e.g. kinases, calmodulin, etc.), Reactive Oxygen Species (ROS) producing or scavenging proteins, or constituents of microtubules. Unfortunately, these two transcriptomic studies concern a unique time-point after stimulation, giving a very static view of the molecular response to mechanical stimulus. In 2004, Kimbrough et al. found 1,691 mechanoresponsive genes on root apices that were transiently (5 s time-lapse) stimulated by moving them back and forth and harvested at different times (0, 2, 5, 15, 30, and 60 min post-stimulation) [18]. To our knowledge, this is the only in-depth kinetic (time-series) study ever made in plants on responses to mechanical stimuli. However, these data cannot be viably extrapolated to the responses of aboveground organs, as the soil and aerial mechanical environments, and the selective pressures they apply to the corresponding organs, are highly distinct.

Despite the tentative general flowchart of physiological and molecular responses to mechanoperception established by Telewski [19], we still have no general understanding of molecular mechanoresponse. The fact is that most molecular thigmomorphogenesis studies were realized on different plant species, submitted to different types of mechanical stimulations, with a huge majority of them using uncontrolled mechanical loadings, all of which makes cross-comparison difficult.

There is an even more striking lack of molecular information for recurring mechanical loads such as wind episodes [6, 20]. In 2010, Martin et al. assessed poplar responses to multiple successive stem bendings mimicking wind, and provided the first evidence of a process of desensitization to mechanical loads at the growth and molecular levels [21]. They found that four early mechanoresponsive genes —encoding two calmodulins, a C2H2 TF, i.e. *PtaZFP2*, and a XTH, respectively— responded with lower intensity to a second stem bending, applied 24 h after a first one. This phenomenon is now referred to as ‘accommodation’ and is thought to be crucial to avoid over-response to wind [6, 22]. In trees, transcriptomic analysis of transgenic poplars overexpressing *PtaZFP2* constituted a first step toward the characterization of this phenomenon at molecular level. This study suggested that part of the accommodation process could be explained by a regulation at transcriptional level [6], thus requiring a more global analysis.

Deciphering the molecular mechanisms driving poplar accommodation to bending stimulation could be crucial to help better understand how trees acclimate to wind. In this study, our aim was to identify new molecular actors involved in mechanoresponse and the accommodation process, and unravel their kinetics of contribution. We thus conducted a time-series search of the transcriptome changes in poplar stem tissue subjected to one or two successive bendings.
# Results

## Kinetic analysis of the transcriptomic responses to a single bending in poplar stem

To unravel the complexity and dynamics of the transcriptomic response of poplar stem to a bending stimulus, we performed microarray gene expression profiling at 0.5 h, 2 h, 24 h or 72 h after a tightly controlled bending.

Statistical analysis revealed 2,663 genes differentially expressed along the time-course (Fig. 1; Additional file 1). These Differentially Expressed Genes (DEG) represented 6% of the whole poplar genome, and 75.6% of them were regulated at just 0.5 h or 2 h post-bending (0.5 h PB or 2 h PB). A third of the DEG observed at 0.5 h PB keep being differentially expressed at 2 h PB. At these earlier time-points, a higher number of DEG was up-regulated. At 24 h PB, the transcriptomic response concerned a smaller set of genes (312 DEG), most of which were down-regulated. At 72 h PB, the number of DEG was equally distributed between up- and down-regulated categories.

As a first examination of the mechanical regulatory network, we considered for each time-point the transcripts presenting the higher amplitude of regulation (Additional file 2). At 0.5 h PB, the 10 most highly up-regulated genes encompassed genes already known to present a short-term response to mechanical stimuli, such as an ortholog of a XTH (Potri.018G095200) [17], or *PtaZFP2* [23], a poplar C2H2 TF (Potri.001G235800). Strongly affected expression was also found for two other TF: one that was related to ethylene signalling (Potri.004G220400) and the other that was an ortholog of RRTF (Potri.009G101900), a redox-responsive member of

the Ethylene Response Factor (ERF) TF subfamily [24]. Interestingly, highly up-regulated genes also included homologs of the Late Embryogenesis Abundant (LEA) hydroxyproline-rich glycoprotein family (e.g. Potri.006G2 04300, Potri.009G158900), including *NHL10/NDR1*, a plasma membrane-embedded protein involved in mechanisms maintaining plasma membrane-cell-wall adhesions in *Arabidopsis* [25].

At 2 h and 72 h PB, the most highly up-regulated genes include many genes related to cell-wall modifications or wood formation. For example, three of the 10 most highly up-regulated genes at 2 h PB (Potri.004G2 10600, Potri.012G015000, Potri.009G012100) encoded homologs of *Arabidopsis* fasciclin-like arabinogalactan proteins (FLA) putatively involved in stem biomechanics in plants [26]. At 72 h PB, Potri.009G012200 and Potri.011G161500, encode, respectively; a FLA protein and a berberine bridge enzyme shown to play a role in lignin monolignol metabolism in *Arabidopsis* [27].

At 24 h PB, the set of genes most strongly up- or down-regulated are predicted to encode proteins involved in cellular metabolism, such as a Homocysteine S-methyltransferase (Potri.008G155900), antioxidant components (Potri.003G066400, Potri.006G141400, Potri.011G140400) or elements of ribosome biogenesis (Potri.014G135600 and Potri.015G079600).

A gene highly down-regulated (Potri.010G069000) at 0.5 h, 2 h, 72 h PB (and to a lesser extent 24 h PB) matched to a homolog of *Arabidopsis PHOSPHATE1* (*PHO1*) encoding a Pi exporter [28], pointing, for the first time, to a role of Pi homeostasis in mechanoresponse.

<figure><img src="image_3.png" /><figcaption>**Fig. 1** Global analysis of DEG after bending in poplar stem. Venn diagram shows the distribution and the overlaps of DEG at 0.5, 2, 24 and 72 h PB. 2,633 different DEG were identified. Numbers in the central part indicate the number of repressed and induced genes at each time-point. Proportion of DEG corresponding to TF oscillates between 13 and 6%</figcaption></figure>
Transcriptional regulation of mechanoresponsive TF-
encoding genes

Among the 2,663 mechanoresponsive genes, 271 corre-
sponded to known or putative TF. As shown in Fig. 1,
the proportion of TF varied over the time-course, repre-
senting 12.66% of the DEG at 0.5 h PB but only 5.90% at
24 h PB. These mechanoregulated TF were classified
into 34 of the 58 TF families described in the *PlantTF*
database [29]. The TF families were differentially repre-
sented depending on time post-bending (see Fig. 2). At
0.5 h PB, NAC-type TF represented a significant fraction
of the regulated TF (22 members), whereas only 5 NAC-
type TF were differentially regulated at 72 h PB. Other
TF families involved in stress responses were highly rep-
resented at early time-points, as illustrated by the 16
WRKY-type TF and the 38 members of the ERF family.
Conversely, TF members of the TALE family were in-
volved at later time-points (24 h and 72 h PB). Members
of the MYB family were also highly represented among
the mechanoresponsive TF, being present both at early
and late time-points, and either up- or down-regulated.
Three members of this MYB family show an interesting
pattern of expression along the time-course by being
down-regulated at 0.5 h and up-regulated at 72 h. These
poplar genes (Potri.003G13200, Potri.007G106100 and
Potri.009G053900) respectively encode orthologs of
MYB103, MYB69 and MYB46 *Arabidopsis* TF, known to
be involved in the regulation of secondary cell wall
biogenesis [30].

Expression clustering and Gene Ontology analysis of
mechanoresponsive genes in poplar

To gain insight into the biological processes underlying
poplar response to bending, two approaches were
combined: gene expression clustering and Gene Ontol-
ogy (GO) enrichment. Using the k-means partitioning
method, the expression profiles of the 2,663 DEG were
distributed into 35 clusters, each represented in Fig. 3 by
their average gene expression. This clustering approach
revealed a high diversity of expression profiles, many of
them supported by a high number of genes (apart from
clusters 1 and 10). The most striking expression pattern,
shown by 11 of the clusters, was an early and transitory
induction at 0.5 h PB. The most highly up-regulated
genes are included in these same clusters (Fig. 3). A
smaller number of genes showed an initial down-
regulation of expression followed by reversion to basal
levels (e.g. cluster 35). Clusters 5, 13, 23, 34, and 35
regrouped genes that were differentially expressed at
more than two time-points across the time-series. Some
interesting expression patterns are also shown by cluster
31 regrouping genes that were down-regulated at
0.5 h PB then up-regulated at 72 h PB, and clusters 29,
30 and 33 regrouping genes showing a regulation

throughout the time-course. To find out whether these
co-expressed genes could be indicative of putative
biological function, each cluster was analysed for enrich-
ment of GO terms in the Biological Process (BP),

<figure><img src="image_4.png" /></figure>
<figure><img src="image_5.png" /></figure>

Molecular Function (MF) and Cellular Component (CC) categories (Additional file 3). A clear picture delivered by this GO analysis is the early activation of responses linked to “defence” or “abiotic stress” shown by different GO enriched-terms found for clusters 1, 2, 3, 4, 6, 7, 8, 12 and 14. Three of the four clusters regrouping the earlier, transitory and most up-regulated genes showed an over-representation of terms related to transcription regulation or TF in both the BP and MF categories (clusters 2, 3, and 4). In contrast, cluster 10 is annotated with GO terms related to “negative regulation of

transcription”. This GO analysis also argues in favour of an early triggering of a “signal transduction” (cluster 5 and 12) involving protein phosphorylation (cluster 8 and 9) and positive regulation of the MAPK cascade (cluster 5). The early-up-regulated genes counted a high prevalence of GO terms related to plant hormones, especially ethylene (see GO:0009873 “ethylene-activated signalling pathway” for clusters 4 and 5) and JA (for example, see GO:0009753 “response to jasmonic acid” for clusters 11 and 13). Interestingly, clusters composed of genes down-regulated 2 h PB (clusters 18 and 19) were enriched with
genes related to photosynthesis or plastid function. For clusters of late-up-regulated genes (72 h PB), the GO BP analysis suggests the involvement of genes contributing to remodelling of the cell wall. For example, clusters 24, 25 and 31 are annotated with the “plant-type cell-wall organization or biogenesis” (clusters 24 and 25, GO:0071669) or “regulation of secondary cell-wall biogenesis” (cluster 31, GO: 2000652) terms. Given its enrichment with terms related to “microtubule-binding” (GO:0008017) and the importance of microtubules in the build-up of the cell wall [31, 32], cluster 29 could also contain genes that are central to this cell-wall remodelling process.

## Time-course expression analysis of candidate genes by qRT-PCR

A RT-qPCRs experiment was performed with the dual aim of validating the reliability of our microarray experiment and refining the dynamics of the biological processes triggered in poplar after stem bending. Expression of 46 selected genes was analysed along a new time-series made with samples collected at regular intervals between 0.5 and 24 h PB (Fig. 4). Each gene was selected to be representative of a cluster on the basis of its expression profile and its putative biological function as determined by GO categorization (Additional file 3).

Globally, except for one gene (PtaEB1B), the PCR data backs up the microarray analysis. Moreover, the refined time-series revealed an even greater complexity of the responses. Given the predominance of GO terms related to “defence response” in our data, we studied a variety of genes involved in this biological process, typically members of the WRKY, NAC, GRAS, C2H2 and ERF TF families. For some of these, the up-regulation observed at early time-points (0.5 h and 2 h PB) is maintained until 9 h PB (e.g. PtaZFP2, PtaNHL10, PtaXET6, PtaWRKY53, PtaRAV) whereas others show a more transient induction confined to the very first hours after bending (e.g. PtaWRKY18, PtaNACO62, PtaGRAS, PtaBZIP53, PtaERF4). Results considering several molecular components of calcium/ROS/Phosphorylation also showed contrasted patterns. Whereas the tested genes involved in phosphorylation cascade signalling (PtaMEK4, PtaCRK4, PtaPP2C) or calcium signalling (PtaCAMTA3) were transiently induced at early stages (before 5 h PB), the up-regulation of PtaACA13, encoding a putative calcium-transporting ATPase enzyme involved in the export of Ca²⁺, was maintained until 9 h PB. An unexpected result obtained with this PCR analysis was the down-regulation of a high number of early-induced genes at late stages (after 13 h) (e.g. PtaHRD, PtaACA13 or PtaGXR3). Considering hormonal signalling, our results suggested that PtaIAA8, homologous to an auxin-regulated TF, and PtaCESTA, homologous to bHLH TF involved in Brassinosteroid

biosynthesis, were down-regulated in a biphasic manner. For JA signalling, PtaJAM2- and PtaJAZ2-like, two predicted JA-induced TF, showed a high but very transient induction within the first 2 h PB. In contrast, the peak accumulation of transcripts of PtaAOS4, a gene encoding an allene oxide synthase probably involved in the JA biosynthetic pathway, was shifted to later stages (5 h–9 h PB). Another set of genes (e.g. PtaLAF6, PtaRPOD1, PtaGLK, PtaSKOR) was rapidly down-regulated (1 h PB) and still repressed until 17 h PB. These genes encompassed functions related to cellular and metabolic processes, such as regulation of photosynthetic processes (PtaLAF6, PtaRPOD1 and PtaGLK) or cellular homeostasis, as illustrated by PtaSKOR, an ortholog of an Arabidopsis gene encoding a Shaker family K⁺ channel. Finally, the genes related to secondary cell-wall biogenesis (PtaMYB69; PtaMYB46), lignin synthesis (PtaRCOMT1; PtaLAC4) or microtubule dynamics (PtaMAP65) were repressed until around the 17th hour before recovering to basal levels.

Combined with GO analysis of the clusters, this time-course analysis provides a new picture of the molecular responses occurring in poplar stems in response to mechanical stimulation that was drawn in a simplified model of over-represented GO terms (BP and MF) in Fig. 5.

## Genome-wide evaluation of the effect of successive mechanical stimulations

To unravel the genome-wide impact of successive mechanical stimuli, we analysed the effect of a second bending performed 24 h after a first one. As shown in Fig. 6a, the genes presenting a modified expression due to the second bending were revealed through a statistical test by comparing their expression levels 0.5 h after the second bending with their expression at 24 h after the first bending. We then compared the level of regulation of these DEG against to their regulation observed at 0.5 h after a single bending. The different gene profiles were then classified through a decision-tree approach (Fig. 6a).

Several regulation patterns were observed. The results showed that among the 1,104 genes transiently regulated at 0.5 h after a single mechanical stimulation, 986 were unaffected by the second bending and 112 were less regulated (e.g. PtaWRKY18 (Potri.006G263600) Fig. 6b, “accommodated gene”). Among the 41 early-regulated genes that were still differentially expressed at 24 h after the first bending, 38 were not regulated by the second treatment and, for 3 genes, the analysis was unable to classify them. Overall, this means that 96% of the DEG regulated after the first mechanical stimulation are not any more or less affected by the second one.

This analysis also revealed 19 newly-regulated genes – i.e. genes that had shown no impact after a single mechanical treatment but were differentially regulated 0.5 h after the
Quantitative transcript relative accumulation

<figure><img src="image_6.png" /><figcaption>**Fig. 4** Expression profiles of selected genes along a refined time-course. Total RNAs were extracted from the stems of control plants and from bent stems at different times (0, 1, 2, 5, 9, 13, 17, 21 and 24 h) after bending. Changes in the abundance of the transcripts were analysed by qRT-PCR using specific primers (Additional file 4). Fold-changes are quantified in comparison to control plants collected at identical time-points. Plain symbols correspond to a significantly different transcript accumulation in bent plants according to a Mann-Whitney-Wilcoxon test</figcaption></figure>

second bending (Additional file 5). Finally, among the 1,145 genes regulated by a single treatment, only 44 genes were similarly regulated after the second bending (e.g. *PtaHRD* (Potri.006G021000) in Fig. 6b). As shown in Table 1, these “non-accommodated genes” were mainly distributed into three of the expression clusters previously defined in the time-series analysis: i.e. clusters 9, 10 and 14, grouping 9, 12, and 10 of these genes respectively. These co-expression profiles could suggest a common regulation.

Interestingly, among the 44 “non-accommodated” genes, several are involved in maintaining the functional integrity of the plant cell wall after biotic/abiotic stresses in other species. Indeed, several genes are predicted to encode members of the CrRLK1 family of RLK in *Arabidopsis*, such as Potri.006G110000, an ortholog of FERO-NIA recently described as a putative mechanosensor in *Arabidopsis* [33], but also ANX2 (Potri.005G055200) and RLK1 (Potri.005G056800). Furthermore, among the
<figure><img src="image_7.png" /><figcaption>**Fig. 5** Time-course of the biological processes involved in response to bending. These annotations correspond to the GO terms that are enriched in clusters of gene expression. Biological Processes are time-placed to match the expression of representative genes</figcaption></figure>

**A**
Differentially expressed 24 hours after 1 bending?

Differential expression between 0.5 h after the 2nd bending and 24 h after the 1st bending?

**ACCOMMODATED** 948 genes

**NON-ACCOMMODATED** 44 genes

**ACCOMMODATED** 112 genes

**ACCOMMODATED** 38 genes

**ACCOMMODATED** 0 genes

**UNKNOWN** 3 genes

<figure><img src="image_8.png" /><figcaption>PtaWRKY18 (Potri.006G263600)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>PtaHRD (Potri.006G021000)</figcaption></figure>

**Fig. 6** Decision tree used to classify the genes into the "accommodated" or "non-accommodated" gene categories. **a** For each gene differentially regulated at 0.5 h PB, "YES" or "NO" status was determined by the statistical test resulting from the microarray analysis described above. **b** Gene examples for the "accommodated" or "non-accommodated" gene categories. The *PtaWRKY18* and *PtaHRD* genes are differentially expressed at 0.5 h PB but not at 24 h PB, whereas the expression of *PtaHRD* is similarly regulated by the two successive bending while, *PtaWRKY18* expression is less affected by the second bending than the first one. According to the decision tree, *PtaWRKY18* is classified into the "accommodated gene" category, and *PtaHRD* into the "non-accommodated gene" category
<table><caption>**Table 1** List of the "non-accommodated" genes after two successive bendings of young poplar stem</caption><thead><tr><th rowspan="2">Cluster</th><th rowspan="2">Id_Gene Poplar</th><th colspan="2">Log2 expression fold-change</th><th rowspan="2">Id_Gene Ath</th><th rowspan="2">Ath Description (Alias)</th></tr><tr><th>1B0.5/C</th><th>2B0.5/1B24</th></tr></thead><tbody><tr><td>4</td><td>Potri.010G041300</td><td>2.92</td><td>1.22</td><td>NULL</td><td>NULL</td></tr><tr><td>6</td><td>Potri.T131200</td><td>2.21</td><td>1.64</td><td>AT3G19680</td><td>Protein of unknown function (DUF1005)</td></tr><tr><td>6</td><td>Potri.008G159300</td><td>2.03</td><td>1.85</td><td>AT3G22930</td><td>Calmodulin-like 11 (CML11)</td></tr><tr><td>7</td><td>Potri.004G086900</td><td>1.28</td><td>0.78</td><td>NULL</td><td>NULL</td></tr><tr><td>7</td><td>Potri.006G055500</td><td>1.39</td><td>1.24</td><td>NULL</td><td>NULL</td></tr><tr><td>7</td><td>Potri.009G084000</td><td>1.38</td><td>1.24</td><td>AT1G10020</td><td>Protein of unknown function (DUF1005)</td></tr><tr><td>7</td><td>Potri.009G092900</td><td>1.44</td><td>1.2</td><td>AT2G14960</td><td>Auxin-responsive GH3 family protein (GH3.1)</td></tr><tr><td>7</td><td>Potri.016G132500</td><td>1.42</td><td>0.93</td><td>AT3G51710</td><td>D-mannose binding lectin protein with Apple-like carbohydrate-binding domain</td></tr><tr><td>8</td><td>Potri.001G309600</td><td>0.89</td><td>1.07</td><td>AT1G49230</td><td>RING/U-box superfamily protein</td></tr><tr><td>8</td><td>Potri.006G110000</td><td>1.2</td><td>0.78</td><td>AT3G51550</td><td>Malectin/receptor-like protein kinase family protein (FER)</td></tr><tr><td>8</td><td>Potri.006G141500</td><td>0.95</td><td>0.91</td><td>AT2G05940</td><td>Protein kinase superfamily protein</td></tr><tr><td>8</td><td>Potri.008G126800</td><td>1.12</td><td>0.89</td><td>AT2G01100</td><td>unknown</td></tr><tr><td>8</td><td>Potri.009G030600</td><td>0.97</td><td>0.81</td><td>AT2G28890</td><td>Poltergeist like 4 (PPL4)</td></tr><tr><td>8</td><td>Potri.010G081100</td><td>1.31</td><td>1.07</td><td>AT3G22910</td><td>ATPase E1-E2 type family protein / haloacid dehalogenase-like hydrolase family protein</td></tr><tr><td>8</td><td>Potri.012G126100</td><td>1.17</td><td>0.88</td><td>AT5G51460</td><td>Haloacid dehalogenase-like hydrolase (HAD) superfamily protein (TPPA)</td></tr><tr><td>8</td><td>Potri.018G093600</td><td>1.22</td><td>0.94</td><td>AT5G57480</td><td>P-loop containing nucleoside triphosphate hydrolases superfamily protein</td></tr><tr><td>8</td><td>Potri.019G128100</td><td>0.93</td><td>1.01</td><td>AT2G30360</td><td>SOS3-interacting protein 4 (CIPK11)</td></tr><tr><td>9</td><td>Potri.001G344100</td><td>0.52</td><td>0.82</td><td>AT5G14500</td><td>Aldose 1-epimerase family protein</td></tr><tr><td>9</td><td>Potri.005G158400</td><td>0.59</td><td>0.79</td><td>AT1G02070</td><td>unknown</td></tr><tr><td>9</td><td>Potri.006G235800</td><td>0.73</td><td>0.98</td><td>AT5G25930</td><td>Protein kinase family protein with leucine-rich repeat domain</td></tr><tr><td>9</td><td>Potri.008G103300</td><td>0.64</td><td>0.83</td><td>AT5G49520</td><td>WRKY DNA-binding protein 48 (WRKY48)</td></tr><tr><td>9</td><td>Potri.010G132800</td><td>0.64</td><td>0.85</td><td>AT3G25600</td><td>Calcium-binding EF-hand family protein</td></tr><tr><td>9</td><td>Potri.013G115700</td><td>0.51</td><td>0.79</td><td>AT1G34300</td><td>Lectin protein kinase family protein</td></tr><tr><td>9</td><td>Potri.015G063500</td><td>0.76</td><td>0.78</td><td>AT3G47730</td><td>ATP-binding cassette A2 (ABCA2,ATH1)</td></tr><tr><td>9</td><td>Potri.017G000700</td><td>0.81</td><td>0.89</td><td>AT2G44130</td><td>Galactose oxidase/kelch repeat superfamily protein</td></tr><tr><td>10</td><td>Potri.006G021000</td><td>3.15</td><td>3.4</td><td>AT2G36450</td><td>Integrase-type DNA-binding superfamily protein (HRD)</td></tr><tr><td>11</td><td>Potri.012G048200</td><td>1.28</td><td>1.07</td><td>AT1G18210</td><td>Calcium-binding EF-hand family protein</td></tr><tr><td>12</td><td>Potri.002G126300</td><td>0.7</td><td>1.09</td><td>AT4G25000</td><td>Alpha-amylase-like (AMY1)</td></tr><tr><td>12</td><td>Potri.004G112700</td><td>0.83</td><td>1.48</td><td>AT5G15720</td><td>GDSL-motif lipase 7 (GLIP7)</td></tr><tr><td>13</td><td>Potri.005G056800</td><td>1.24</td><td>0.85</td><td>AT5G60900</td><td>Receptor-like protein kinase 1 (RLK1)</td></tr><tr><td>13</td><td>Potri.018G076500</td><td>1.31</td><td>0.84</td><td>AT3G24440</td><td>Fibronectin type III domain-containing protein (VIL1,VRN5)</td></tr><tr><td>14</td><td>Potri.003G159800</td><td>0.99</td><td>1.2</td><td>AT5G47910</td><td>Respiratory burst oxidase homologue D (RBOHD)</td></tr><tr><td>14</td><td>Potri.005G015300</td><td>0.9</td><td>0.86</td><td>NULL</td><td>NULL</td></tr><tr><td>14</td><td>Potri.005G055200</td><td>0.88</td><td>1.06</td><td>AT5G28680</td><td>Malectin/receptor-like protein kinase family protein (ANX2)</td></tr><tr><td>14</td><td>Potri.006G187900</td><td>0.89</td><td>0.82</td><td>AT2G24300</td><td>Calmodulin-binding protein</td></tr><tr><td>14</td><td>Potri.006G219100</td><td>1.03</td><td>1.04</td><td>AT1G72510</td><td>Protein of unknown function (DUF1677)</td></tr><tr><td>14</td><td>Potri.009G002200</td><td>0.82</td><td>0.95</td><td>AT5G22860</td><td>Serine carboxypeptidase S28 family protein</td></tr><tr><td>14</td><td>Potri.010G240800</td><td>0.5</td><td>0.89</td><td>AT5G04760</td><td>Duplicated homeodomain-like superfamily protein</td></tr><tr><td>14</td><td>Potri.013G125000</td><td>0.51</td><td>1.04</td><td>AT3G54420</td><td>Homolog of carrot EP3-3 chitinase (CHITIV)</td></tr><tr><td>14</td><td>Potri.013G130000</td><td>1.02</td><td>0.97</td><td>AT1G53050</td><td>Protein kinase superfamily protein</td></tr></tbody></table>
<table><caption>**Table 1** List of the "non-accommodated" genes after two successive bendings of young poplar stem (Continued)</caption><tr><td>14</td><td>Potri.014G111800</td><td>0.67</td><td>1.23</td><td>AT4G01700</td><td>Chitinase family protein</td></tr><tr><td>14</td><td>Potri.017G013300</td><td>0.87</td><td>1.09</td><td>AT4G28400</td><td>Protein phosphatase 2C family protein</td></tr><tr><td>15</td><td>Potri.006G101400</td><td>1.75</td><td>1.5</td><td>AT5G01380</td><td>Homeodomain-like superfamily protein</td></tr></table>

These genes are differentially expressed 0.5 h after a single bending and their regulation is not different whether one or two successive bendings are applied. *Id gene* corresponds to the gene model found in the Phytozome database [55] for poplar or *Arabidopsis* (Ath). Log2 expression fold-change is shown for the different analyses: expression fold-change 0.5 h after a single bending compared to a control plant (1B0.5 h/Control) and expression fold-change after a second bending compared to the 24 h post-bending time-point (2B0.5 h/1B24h)

other "non-accommodated" genes, some were described as components of signalling cascades coordinating cellular responses to help maintain plant cell-wall integrity [34]. This is the case for Potri.003G159800 encoding the respiratory burst oxidase homolog D, responsible in *Arabidopsis* for apoplast ROS production [35] or for genes encoding calcium-related proteins such as calcium-binding EF-hand family protein (Potri.012G048200, Potri.010G132800), Calmodulin-binding protein (Potri.006G187900 Potri.008G159300) or CBL-interacting protein kinase (Potri.019G128100).

## Discussion

### Short- and long-term impacts of transient stem bending on transcriptional responses in poplar

Lee et al. studying the transcriptomic effect of a touch stimulus at a single time-point (30 min) in the herbaceous *Arabidopsis* showed that the expression of over 2.5% of the *Arabidopsis* genome was up-regulated at least twofold [16]. Here, in poplar, a single transitory bending modulated the expression of about 6% of the genome over 4 time-points, and strikingly, the majority of changes in gene expression were observed within the first two hours post-bending. Such massive early response has been seen in other transcriptomic time-course studies of responses to stress (e.g. defence response to *Botrytis cinerea* in [36]). Moreover, at 0.5 h PB, genes are mainly up-regulated (782 up- vs. 362 down-regulated genes). A high proportion of these early-regulated genes showed no differential regulation at 2 h PB and thereafter. The transcripts they produce could be classified as unstable transcripts according to Gutierrez et al. [37]. Actually, about 20% of the unstable transcripts described by these authors were up-regulated at 0.5 h PB in our experiment (data not shown).

Our study illustrates the massive extent of transcriptomic changes in response to a non-wounding mechanical stimulus, but also sheds new light on the duration and dynamics of the phenomenon. The refinement of the time-course by qPCR for 46 candidate genes brought key insight by revealing unexpected short or long-lasting expression profiles. Notably, this analysis highlighted a high number of early-induced genes presenting inhibited expression at later time-points. This suggests the existence of an active negative feedback mechanism.

### A rapid activation of defence signalling

The association of a GO term enrichment analysis with a clustering approach enabled us to distinguish different major components of the response. Numerous clusters are enriched with genes related to plant responses to stimuli, either biotic or abiotic. A majority of these genes were differentially expressed at early time-points, showing the ability of trees to react rapidly to mechanical stimuli by modifying the expression of genes involved in the general response to abiotic stresses. These genes encompassed for example an ortholog of the PIA1 PP2C protein phosphatase known to regulate defence responses both positively and negatively in *Arabidopsis* [38], and an ortholog of the *Arabidopsis* CAMTA3 calmodulin-binding TF, which was recently shown to regulate amplitude of Rapid Stress Response Element (RSRE) peak activity in the general stress response, "the initial rapid cascade preceding the induction of consequent stress-specific signalling activities" [39].

In plants, abiotic stresses usually trigger the accumulation of ROS and increase the expression of scavenging enzyme-encoding genes [40]. Here, early responses to bending also included the regulation of genes related to ROS scavenging, such as several peroxidase-encoding genes, or *PtaGrx3* (Potri.017G017300), a member of the glutaredoxin family that regulates protein redox state. It is known that plants can respond to stresses by a coordinated modulation of chloroplast and nuclear gene expression. Such rapid acclimation of plant cell metabolism and photosynthetic machinery is considered a key factor for plant survival. Here, the expression of several genes involved in chloroplast function or development was down-regulated for several hours after the bending. To our knowledge, these results are the first evidence of coordinated modulation following a non-wounding mechanical stimulus.

Taken together, these data suggest that in trees, a transient non-wounding mechanical treatment also quickly activates defence-signalling pathways.

Previous studies have argued that $Ca^{2+}$ may be involved in plant responses to mechanical stimuli [16]. Accordingly, several genes encoding proteins related to calcium such as calmodulins, Calcineurin B-Like calcium sensor proteins (CBLs) and Autoinhibited $Ca^{2+}$-ATPases (ACAs) were also deregulated in poplar stem tissues after bending. The
long-term maintenance of the induced expression of *PtaACA13* suggests its role in maintaining $Ca^{2+}$ homeostasis thereby preventing cell injuries or maintaining $Ca^{2+}$ second messenger functions.

Plant hormones are known to play major roles in development and defence against biotic and abiotic stresses, so it was unsurprising to find concomitant modifications in expression of genes related to defence and to hormones at early time-points. In *Arabidopsis*, ethylene and JA are known to strongly contribute to the mechanical signalling pathway in the first-stage responses. Many of the genes studied here encode TF related to ethylene signalling and many others encode proteins involved in JA biosynthesis (Allene Oxide Cyclases, Allene Oxide Synthases, JA-deconjugating enzymes, JA sulfo-transferases) or JA signalling (JAZs and “JA-regulated” genes). These observations suggest JA and ethylene are central to mechanoresponse in the poplar model tree, just as they are in ruderal herbs such as *Arabidopsis thaliana*.

## Stem bending elicits transcriptional responses linked to cell wall modification and/or wood development

Poplar stem bending caused a complete cessation of cambial growth during the first hours [4]. The coincidence of early growth cessation in bent poplar and transitory overexpression of genes related to plant defence is consistent with the reported antagonism between plant immunity and growth [41]. In this tradeoff, hormones, especially JA, play a significant role. After the early phase of growth cessation, cambial growth increases during 3 to 5 days before resuming a normal rate of growth. In our study, 10 probesets revealed an increased expression of ribosomal protein-encoding genes at the 24 h PB time-point, which was further confirmed by the clustering approach with cluster 20 showing an enrichment of GO terms related to ribosomes. Sustained accelerated growth, as observed in bent poplar, is heavily reliant on an efficient supply of proteins by the cell machinery [42]. Thus, the effect of bending on the expression of genes encoding ribosomal protein could accompany the strong metabolic activity necessary for accelerated cambial growth.

Although early responses of poplar to a bending stimulus seem to apply to defence and growth, this transcriptionomic study indicated that later responses strongly affect cell wall and/or wood development. Indeed, a significant portion of cell wall-related genes displayed differential expression at different time-points. FLA genes were up-regulated at 2 h and 72 h PB while some NAC and MYB encoding genes, known to be part of the transcriptional network that regulates secondary cell formation and wood cell differentiation [43, 44], were first down-regulated at 0.5 h then up-regulated at 72 h PB. Moreover, the most highly induced gene corresponded

to a homolog of *Arabidopsis* XTR6, a XTH. In the cell wall, these proteins are thought to specifically catalyse the endolytic cleavage and re-ligation of xyloglucan chains, thus allowing adjustments of cell wall mechanical properties. These observations may indicate the existence of a two-phase process: (i) a relatively rapid response that might be interpreted as a way for plants to modulate growth through adjustment of the expansibility of the primary cell wall followed by (ii) a later response that may set up tissue differentiation and changes in properties to enable acclimation to the encountered mechanical stress. Indeed, multiple stem bendings in poplar were previously shown to modify the mechanical properties (cell wall fraction, elastic stiffness, flexural strength) of the stems [3, 5]. Until now, PtaZFP2 was the only protein putatively involved in such modifications [23]. Our study thus represents a first tentative molecular explanation for this two-phase effect of bending loads.

## A single mechanical stimulus is sufficient to attenuate plant responsiveness to recurrent stimulus

In natural conditions, wind triggers successive bendings of aerial axes. In a previous study, we demonstrated that plants reduce their responsiveness to mechanical signals as a function of their mechanical history—a phenomenon named accommodation [6, 21]. Here, we aimed to unravel the impact of repeated stem bending at genome-wide level by comparing the magnitude of transcriptional responses between two successive bendings and a single bending of the stem. This kinetics is intended to mimic the alternation between windy and quiet days observed in temperate climates [45] more than the frequencies (between 1 and 5 Hz) at which plants oscillate in wind [46]. Our transcriptionic results highlighted the existence of three classes of transcriptional response patterns.

First, the huge majority of these early mechanorespondive genes (96%) presented a smaller amplitude of regulation 0.5 h after the second bending. Furthermore, the data showed that most of these early-regulated genes returned to their basal expression levels at 24 h PB, at the time when the second bending was applied. The attenuation of transcriptional responses to the second bending thus probably results from a reduction in transcriptional responsiveness and not from saturation of the transcriptional machinery. These data suggest a global desensitization of the first stages of the mechanotransduction pathway.

The description of genes presenting altered responses to successive stimulations has been made elsewhere for other stresses. For example, Ding et al. used RNA-seq analysis to classify drought-responsive genes as “memory genes” or “non-memory” genes as a function of their transcriptional responses to two successive dehydration treatments [47]. Plants can record repeated stimulus and
respond in consequence, a phenomena named “attenu-
ation”, “stress imprint” or, as here, “accommodation”.
However, contrary to mechanical stem loading, re-
sponses are usually enhanced after a second exposure in
most of the studies [48]. Whereas gene priming and
transcriptional memory are classically thought to help
provide an increased, or at least improved, response to a
biotic or abiotic stress, accommodation mostly switches
off or attenuates the response. Such fine-tuning of the
mechanotransduction pathway could help prevent un-
necessary thigmomorphogenetic growth reductions
while still allowing the plant to become more resistant
(acclimated) to the stimulus.

Two unexpected categories of transcriptional responses were observed in our study: (i) the “newly-regulated” genes (19 genes), whose expression was only regulated after the second stimulation, and (ii) the “non-accommodated” genes (44 genes), i.e. mechanoresponsive genes that respond similarly to the first and the second stimulus. These data raise the question of the nature of a mechanism, or set of mechanisms, involved in this phenomenon without affecting all the mechanoresponsive genes similarly. As stated in [20], several hypotheses can be put forward: alteration of perception through direct modification of mechanosensors, or a negative feedback control on the mechanotransduction pathway (e.g. through a negative transcription control or epigenetic modifications). The existence of three types of expression profiles makes the hypothesis of a direct control on mechanosensors difficult to consider, except if several types of mechanosensors are involved. Interestingly, a recent study suggested a putative role for FERONIA in mechanosensing [33]. However, it also revealed that the pattern of cytoplasmic changes in [Ca²⁺] observed after a mechanical stimulation cannot be explained solely by FERONIA, which requires other mechanosensors to be involved. In our case, 10% of the “non-accommodated” genes encode proteins possessing a kinase domain, among which the poplar orthologous gene of FERONIA. The non-accommodation of these genes could be a way for the cells to help record recurrent stimuli or to fine-tune cell sensitivity to mechanical stimulation. However, the hypothesis of an active negative feedback control at the transcriptional level also needs to be explored, as suggested by the active inhibition of expression observed at later time-points for most of the early-induced genes.

Conclusion

This genome-wide study of poplar transcriptional re-
sponsiveness to bending brings a novel temporal dimen-
sion to our knowledge of plant responses to mechanical
stimuli. Moreover, it highlights the extent of the accom-
modation process and the existence of diverse gene
regulation types in response to recurring mechanically
stressing events. The identification of the underlying

molecular mechanisms opens news perspectives for
gaining deeper understanding of how the intensity and
recurrence of a stimulus affect plant growth and devel-
opment. We anticipate this study as a start point to the
generation of predictive models of the gene regulatory
networks mediating poplar transcriptional responses to
mechanical stimulation.

Methods

Plant material, culture conditions and bending treatments

Wild-type poplars (*Populus tremula* × *P. alba* cv 717-1B4) were obtained by *in vitro* micropropagation on MS medium [49]. After acclimation, plants were grown in liquid nutrient solution [50] in a growth chamber at 22 °C, with a relative air humidity of 60% and a 16 h/8 h light/ dark cycle under photosynthetic active radiation (PAR) of 50 µmol.m⁻².s⁻¹.

Analyses were conducted on 3-month-old poplars.
Bendings were realized according to [23]. Briefly, each
basal stem part was bent transiently (10 s) against a plas-
tic tube so as to apply controlled and quantified flexural
strains as described in [4] (Additional file 6).

Two different types of bending treatment were applied. The first was a single bending treatment, for which the bent portion of the stem was collected from bent and control plants at 0.5, 2, 24 and 72 h after the mechanical stimulation, thus yielding a time-series expression dataset.

The second type of treatment consisted of two succes-
sive transient bendings separated by a 24 h time interval,
for which the bent portion of the stem was collected at
0.5 h after the second mechanical stimulation.

Microarray analysis

Three independent experiments (i.e. three biological repli-
cates) on two individuals per experiment and per condition
were analysed using Affymetrix GeneChip® Poplar Genome
Array oligonucleotide microarrays (Santa Clare, Canada).

RNA was isolated from 100 mg of bent or unbent stem
using the RNeasy Plant Mini kit followed by a DNAsel
treatment according to the manufacturer's instruction
(Qiagen,Germany). RNA sample integrity was checked on
an Agilent 2100 Bioanalyzer (Agilent Technologies,
Germany). The amplification, labelling, hybridization
and imaging procedures were performed according to
the manufacturer's instructions [51] at the URGV
Transcriptomics Platform (Evry, France).

For each independent experiment, arrays were hybrid-
ized with complementary RNA obtained from 2 µg of
total RNA made of a pool from two individuals. Arrays
were scanned with the GeneChip 3000-7G Scanner
driven by the GeneChip Operating Software (GCOS).
Raw data were imported into R software for analysis and
normalized by the GC-RMA algorithm [52] available in
the Bioconductor package.
A two-group t-test assuming equal variance between groups was performed to determine which genes were differentially expressed in the following comparisons: (i) unbent vs. bent stems at each time-point of the time-series, (ii) stems bent once and harvested 0.5 h after the mechanical stimulation vs. stems bent twice, (iii) stems bent once and harvested 24 h after the mechanical stimulation vs. stems bent twice. To fit the assumption of equal variance of gene expression between groups, genes displaying extreme variances (too small or too large) were excluded. The raw P-values were adjusted by the Benjamini–Hochberg procedure [53], which controls false discovery rate. A probeset was declared as representing differential expression of a gene if Benjamini–Hochberg P-value was < 0.05. All the raw and normalized data are available through the Gene Expression Omnibus repository at the NCBI [GEO submission GSE44321; [54]]. In order to limit the tendency of this Benjamini–Hochberg correction to lead to false negatives, the data was pre-treated prior to statistical analysis. When a probeset presented a log-transformed intensity under the threshold of 4 for each repeated sample in one comparison, it was considered as unexpressed and these samples were removed from the analysis.

## Microarray probeset selection and annotation up-dating

The Affymetrix 61 k poplar microarray (composed of 61,251 probesets) was designed on the first version of *Populus Trichocarpa* genome, but also on mRNA and EST from different species of *Populus* (e.g. *P. tremula*, *P. deltoids*, *P. alba*) [51]. Consequently, some genes are targeted by more than one probeset. In the Affymetrix poplar GeneChip, only 585 probesets corresponded to sequences of the *P. alba* x *P. tremula* hybrid used in our analysis. To reduce multiple gene probe elements without impacting the global analysis, probesets designed from *P. alba*, *P. tremula* and *P. trichocarpa* (the sequenced genome) and hybrids from these 3 species were also included in this analysis —representing 57,654 probesets of the 61,251 initial probesets of the Affymetrix 61 k poplar microarray.

To update the Affymetrix probeset annotations to the third version of the Phytozome *P. trichocarpa* genome, a BLASTn search was performed between the probeset sequences and the third version of *P. trichocarpa* sequencing (ncbi-blast-2.2.28 – 2013.04.01 version). 50,691 of the selected probesets were successfully associated to 31,121 different gene IDs among the 41,335 genes of the third version of *P. trichocarpa* genome. The *Arabidopsis thaliana* ortholog of each *P. trichocarpa* gene targeted by a probeset was retrieved in the Phytozome database of the third version of the *P. trichocarpa* genome (file: Ptrichocarpa_210 annotation info.txt) as "best *Arabidopsis* TAIR10 hit symbol" [55].

## Clustering of the microarray data

To identify groups of genes presenting similar patterns of expression in response to one bending and/or similar expression behaviour toward a second bending, we performed k-means clustering on the DEG (i.e. genes that are differentially expressed in at least one condition). We used the Euclidian metric to calculate the distances between the DEG, as follows. Let *x* and *y* two genes defined by *n* values (and *i* one of these *n* values). The distance between *x* and *y* is given by $d(x, y) = \sqrt{\sum_{i=1}^{n}(x_i - y_i)^2}$. Three types of *n* values were chosen. Two of them allowed genes to group according to their pattern of expression along the time-course: (i) the speed of transcript accumulation between the different time points, given by $v(x_{\Delta t}) = \frac{a_{t+1}+a_t}{\Delta t_{(t,t+1)}}$ where $a_{t+1}$ and $a_t$ are the differences of expression measured for gene *x* at *t* + 1 and *t*, and $\Delta t$ the lapse of time between the two time-points. The 0 h PB and 120 h PB time-points were also used and their level of differential expression was set to 0; (ii) a scale grade allowing discrimination between significant and non-significant differential expression values and expressing the type of regulation: 1 for a significant up-regulation, -1 for a significant down-regulation, and 0 for a non-significant regulation. The third type of *n* value, corresponding to the subtraction of expression ratios observed after one bending to expression ratios observed after two bendings, aimed to group genes according to their behaviour in response to the second bending.

## GO enrichment analysis

To categorize the DEG based on their biological functions, GO enrichment was searched for each expression cluster. In the third version of the poplar genome, genes are associated to a narrow numbers of annotations – e.g. *PtaZFP2* (Potri.001G235800) is annotated with 2 GO terms in the Phytozome database whereas its homologous *A. thaliana* gene has 10 GO annotations on Phytozome. Consequently, the GO enrichment of each cluster was performed with the term-enrichment tool from Amigo (1.8) using the *A. thaliana* orthologs and the TAIR database as a background. The maximum-thresholds P-value was set at 0.05 but most of the clusters presented P-values under 0.01.

## Selection of genes as representatives of each cluster

For each cluster, one gene representing the other genes of the cluster was selected based on two criteria: (i) its expression should be as close as possible to the median expression value of a cluster; (ii) it should be representative of either the GO functional category enrichment found for the cluster or a well-known component of the mechanical signalling pathway.
## Real-time quantitative RT-PCR experiments on candidate genes

The time-course analysis of candidate genes expression was performed via two independent experiments (i.e. two biological replicates) where young poplars were treated following the conditions described for microarray analysis. For each time-point, each independent experiment includes two mechanical-stimulated plants and one control. The bent portion of the poplar stem was collected 1, 2, 5, 9, 13, 17, and 21 h after mechanical stimulation.

Total RNAs was extracted from about 100 mg of bent stems using CTAB extraction buffer as described in Chang et al. [56] then treated with RNase-free RQ1 DNase (Promega, Charbonnières-les-Bains, France). RNA was quantified spectrophotometrically using a NanoDrop apparatus (Thermo Fisher Scientific Inc., Waltham, MA) and checked by agarose gel electrophoresis. First-strand cDNA was synthesized from 1 µg total RNAs using oligodT and SuperScript III (Invitrogen, Cergy-Pontoise, France). The real-time quantitative RT-PCR amplifications were performed using an iCycler IQ (Bio-Rad) with SYBR green as fluorescent dye. Each PCR reaction (15 µL) contained the following: cDNA (4 µL of a 1:40 dilution of the first cDNA strands), MESA GREEN qPCR MasterMix Plus for SYBR® Assay w/ fluorescein (Eurogentec, Angers, France) (1x) and primers (200nM of each). After a heat step at 95 °C for 3 min, PCR cycling conditions were 40 cycles of denaturation (95 °C, 10 s), annealing (temperature according to primers, Additional file 4, 10 s) and elongation (72 °C, 15 s), ending with a melt curve (increment of 0.5 °C each cycle from primer tm to 95 °C). Transcripts of each studied gene and reference genes were amplified using the primers described in Additional file 4. Specificity of amplification was confirmed by (i) determining the melt curves for the PCR products, (ii) gel electrophoresis and (iii) sequencing. An Index corresponding to the expression of 5 reference genes (EF-1α, UP1, UP2, TIP41 and UBC) was calculated. The Quantitative Relative Normalized (QRN) abundance of each gene transcript was calculated by comparison between bent and control plants and by normalizing with the expression of the index using the delta-delta method mathematical model [57].

Each run of the real-time PCR amplifications was carried out in triplicate. RNAs samples used for microarray, arising from control plants or plants collected 0.5, 2, or 24 h after bending, were also included in this Q-PCR analysis.

## Additional files

**Additional file 1: Table S1.** List of genes differentially expressed probset in response to one or two bendings, as identified by microarray analysis. Affymetrix 61 k poplar microarrays (composed of 61,251 probesets) were used to detect differences in expression. Four comparisons between bent plants (b) vs. control plants (c) at different

time points after the bending were realized: 0.5 h PB (1b0.5 h vs. C0.5 h), 2 h PB (1b2h vs. C2h), 24 h PB (1b24h vs. C24h) and 72 h PB (1b72h vs. C72h). To understand the effect of a second bending, we ran two comparisons (i) 0.5 h after the 2nd bending (with a time-lapse of 24 h between the two bendings) compared to 24 h after the 1st bending (2b0.5 h vs. 1b24h) and (ii) 0.5 h after the 2nd bending compared to 0.5 h after the 1st bending (2b0.5 h vs. 1b0.5 h). Expression was considered significantly different if the *P*-value from the Benjamini–Hochberg test (BH *P*-value) is below 0.05. (CSV 694 kb)

**Additional file 2: Table S2.** List of the 10 most up- or down-regulated genes at each time-point. *gene_ID corresponds to the gene model found in Phytozome database [55] for poplar or *Arabidopsis* (Ath). The log2 expression fold-change is shown for the different time-points: 0.5 h (exp 0.5 h), 2 h (exp 2 h), 24 h (exp 24 h) and 72 h (exp 72 h). (XLSX 15 kb)

**Additional file 3: Table S3.** Gene Ontology enrichment analysis applied to the clusters. The GO enrichment was performed using the Amigo (1.8) term enrichment tool [58]. The list of *A. thaliana* orthologs of each cluster was used, and the TAIR database was chosen as a background. The maximum *P*-value for the threshold was 0.01, except for clusters marked * and ** for which maximum *P*-value was set at 0.05 and 0.1, respectively. (XLSX 13 kb)

**Additional file 4: Table S4.** Description of the genes analysed in the refined time-course analysis. Primer sequences are described. *Id gene* corresponds to the gene model found in the Phytozome database [55] for poplar or *Arabidopsis* (Ath). **Poplar gene name was attributed arbitrarily following the gene function described for the *Arabidopsis* homologous gene in Phytozome. Gene membership to expression clusters is specified. (XLSX 14 kb)

**Additional file 5: Table S5.** List of "newly regulated" genes after two successive bendings. These genes are not regulated 0.5 h after a single stem bending, but they do become regulated 0.5 h after a second bending (occurring 24 h after the 1st one). *Id gene* corresponds to the gene model found in the Phytozome database [55] for poplar or *Arabidopsis* (Ath). Log2 expression fold-change is shown: expression fold-change 0.5 h after a single bending compared to a control plant (1B0.5 h/Control) and expression fold-change 0.5 h after a second bending compared to the 24 h PB time-point (2B0.5 h/1B24h). (XLSX 11 kb)

**Additional file 6: Figure S1.** Bending device. The stem is transitory – 10s back and forth– pushed against a plastic tube of known diameter; the stem is thus locally bent around the tube. Therefore, the stem is subjected to a quantified curvature in terms of strains. Locally, the applied strain is the product of the curvature of the central line and the stem radius ($r_{\text{stem}}$). In the case of small curvature, the curvature of the central line is given by the inverse of the sum of the stem radius and of the radius of the plastic tube ($r_{\text{tube}}$). (Modified from [4]). (PDF 747 kb)

## Abbreviations

BP: Biological process; DEG: Differentially expressed genes; ERF: Ethylene response factor; GO: Gene Ontology; JA: Jasmonic acid; MF: Molecular function; MS: MechanoSensitive; PB: Post-bending; RLK: Receptor-like kinase; ROS: Reactive oxygen species; TF: Transcription factor; XTH: Xyloglucan endo-Transglycosylase/Hydrolase

## Acknowledgements

The authors are grateful to Christelle Boisselet for her help on plant production and Sylvaine Labernia for her help on RNA extraction. We thank Stéphanie Huguet and Sandrine Balzergue (URGV, Evry) for assistance on the generation of Affymetrix microarray data and Marie-Laure Magniette (URGV, EVRY) for assistance with statistical analysis of the microarray data.

## Funding

This research was funded by a French Ministry of Research grant to Lise Pomiès and received support from the Agence Nationale de la Recherche [National Agency for Research] via grant ANR-09- BLAN-0245-01. Funding bodies were not implied in the design and analysis of the study.

## Availability of data and materials

The microarray datasets generated and analysed during the current study are available in the Gene Expression Omnibus repository at the
NCBI [GEO accession: GSE44321], [http://www.ncbi.nlm.nih.gov/geo/query/
acc.cgi?acc=GSE44321].

All qPCR data generated and analysed during this study are included in this
published article and its Additional files.

**Authors' contributions**

LP performed data mining, data analysis and wrote the manuscript. JF
performed the qPCR experiments and analysed the corresponding results.
MD and NLF conceived the study, coordinated the research and wrote the
manuscript. BM conceived the study and coordinated the research. All
authors have read and approved the final manuscript.

**Competing interests**

The authors declare that they have no competing interests.

**Consent for publication**

Not applicable.

**Ethics approval and consent to participate**

*P. tremula* x *P. alba*, cv 717-1B4 is an hybrid poplar, used for transgenic research, obtained from a cross between a female *Populus tremula* tree and a male *Populus alba* tree by the French research institute INRA. No special permission is required for its use.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in
published maps and institutional affiliations.

Received: 23 July 2016 Accepted: 29 March 2017

Published online: 17 April 2017

**References**

1. Telewski FW, Pruyn ML. Thigmomorphogenesis: a dose response to flexing in Ulmus americana seedlings. Tree Physiol. 1998;18:65-8.

2. Coutand C, Dupraz C, Jaouen G, Ploquin S, Adam B. Mechanical stimuli regulate the allocation of biomass in trees: demonstration with young Prunus avium trees. Ann Bot. 2008;101:1421-32.

3. Pruyn ML, Ewers BJ, Telewski FW. Thigmomorphogenesis: changes in the morphology and mechanical properties of two Populus hybrids in response to mechanical perturbation. Tree Physiol. 2000;20:535-40.

4. Coutand C, Martin L, Leblanc-Fournier N, Decourteix M, Julien J-L, Moulia B. Strain mechanosensing quantitatively controls diameter growth and PtaZFP2 gene expression in poplar. Plant Physiol. 2009;151:223-32.

5. Kern KA, Ewers FW, Telewski FW, Koehler L. Mechanical perturbation affects conductivity, mechanical properties and aboveground biomass of hybrid poplars. Tree Physiol. 2005;25:1243-51.

6. Moulia B, Der Loughian C, Bastien R, Martin O, Rodríguez M, Gourcilleau D, et al. Integrative Mechanobiology of Growth and Architectural Development in Changing Mechanical Environments. In: Wojtaszek P, editor. Mech. Integr. Plant Cells Plants. Berlin: Springer Berlin Heidelberg; 2011. p. 269-302.

7. Coutand C, Julien JL, Moulia B, Mauget JC, Guitard D. Biomechanical study of the effect of a controlled bending on tomato stem elongation: global mechanical analysis. J Exp Bot. 2000;51:1813-24.

8. Monshausen GB, Haswell ES. A force of nature: molecular mechanisms of mechanoperation in plants. J Exp Bot. 2013;64:4663-80.

9. Braam J, Davis R. Rain-, wind-, and touch-induced expression of calmodulin and calmodulin-related genes in Arabidopsis. Cell. 1990;60:357-64.

10. Botella JR, Arteca JM, Somodevilla M, Arteca RN. Calcium-dependent protein kinase gene expression in response to physical and chemical stimuli in mungbean (Vigna radiata). Plant Mol Biol. 1996;30:1129-37.

11. Mizoguchi T, Irie K, Hirayama T, Hayashida N, Yamaguchi-Shinozaki K, Matsumoto K, et al. A gene encoding a mitogen-activated protein kinase kinase kinase is induced simultaneously with genes for a mitogen-activated protein kinase and an S6 ribosomal protein kinase by touch, cold, and water stress in Arabidopsis thaliana. Proc Natl Acad Sci U S A. 1996;93:765-9.

12. Gilmour SJ, Zarka DG, Stockinger EJ, Salazar MP, Houghton JM, Thomashow MF. Low temperature regulation of the Arabidopsis CBF family of AP2 transcriptional activators as an early step in cold-induced COR gene expression. Plant J. 1998;16:433-42.

13. Mauch F, Kmec A, Schaffrath U, Volrath S, Görlach J, Ward E, et al. Mechanosensitive expression of a lipoxygenase gene in wheat. Plant Physiol. 1997;114:1561-6.

14. Arteca JM, Arteca RN. A multi-responsive gene encoding 1-aminocyclopropane-1-carboxylate synthase (ACS6) in mature Arabidopsis leaves. Plant Mol Biol. 1999;39:209-19.

15. Gadea J, Conejero V, Vera P. Developmental regulation of a cytosolic ascorbate peroxidase gene from tomato plants. Mol Gen Genet. 1999;262:212-9.

16. Lee D, Polisensky DH, Braam J. Genome-wide identification of touch- and darkness-regulated Arabidopsis genes: a focus on calmodulin-like and XTH genes. New Phytol. 2005;165:429-44.

17. Fluch S, Olmo CC, Tauber S, Stierschneider M, Kopecky D, Reichenauer TG, et al. Transcriptomic changes in wind-exposed poplar leaves are dependent on developmental stage. Planta. 2008;228:757-64.

18. Kimbrough JM, Salinas-Mondragon R, Boss WF, Brown CS, Sederoff HW. The fast and transient transcriptional network of gravity and mechanical stimulation in the Arabidopsis root apex. Plant Physiol. 2004;136:2790-805.

19. Telewski FW. A unified hypothesis of mechanoperation in plants. Am J Bot. 2006;93:1466-76.

20. Leblanc-Fournier N, Martin L, Lenne C, Decourteix M. To respond or not to respond, the recurring question in plant mechanosensitivity. Front Plant Sci. 2014;5:401.

21. Martin L, Leblanc-Fournier N, Julien J-L, Moulia B, Coutand C. Acclimation kinetics of physiological and molecular responses of plants to multiple mechanical loadings. J Exp Bot. 2010;61:2403-12.

22. Bonnesoeur V, Constant T, Moulia B, Fournier M. Forest trees filter chronic wind-signals to acclimate to high winds. New Phytol. 2016;210:850-60.

23. Martin L, Decourteix M, Badel E, Huguet S, Moulia B, Julien J-LL, et al. The zinc finger protein PtaZFP2 negatively controls stem growth and gene expression responsiveness to external mechanical loads in poplar. New Phytol. 2014;203:168-81.

24. Khandelwal A, Elvitigala T, Ghosh B, Quatrano RS. Arabidopsis transcriptome reveals control circuits regulating redox homeostasis and the role of an AP2 transcription factor. Plant Physiol. 2008;148:2050-8.

25. Knepper C, Savory E, Day B. Arabidopsis NDR1 is an integrin-like protein with a role in fluid loss and plasma membrane-cell wall adhesion. Plant Physiol. 2011;156:286-300.

26. MacMillan CP, Mansfield SD, Stachurski ZH, Evans R, Southerton SG. Fasciclin-like arabinogalactan proteins: specialization for stem biomechanics and cell wall architecture in Arabidopsis and Eucalyptus. Plant J. 2010;62: 689-703.

27. Daniel B, Pavkov-Keller T, Steiner B, Dordic A, Gutmann A, Nidetzky B, et al. Oxidation of monolignols by members of the berberine bridge enzyme family suggests a role in plant cell wall metabolism. J Biol Chem. 2015;290: 18770-81.

28. Hamburger D, Rezzonico E, MacDonald-Comber Petétot J, Somerville C, Poirier Y. Identification and characterization of the Arabidopsis PHO1 gene involved in phosphate loading to the xylem. Plant Cell. 2002;14:889-902.

29. Jin J, Zhang H, Kong L, Gao G, Luo J. PlantTFB 3.0: a portal for the functional and evolutionary study of plant transcription factors. Nucleic Acids Res. 2014;42:D1182-7.

30. Ambawat S, Sharma P, Yadav NR, Yadav RC. MYB transcription factor genes as regulators for plant responses: An overview. Physiol Mol Biol Plants. 2013;19:307-21.

31. Wasteneys GO. Progress in understanding the role of microtubules in plant cells. Curr Opin Plant Biol. 2004;7:651-60.

32. Oda Y. Cortical microtubule rearrangements and cell wall patterning. Front Plant Sci. 2015;6:236.

33. Shih H-W, Miller ND, Dai C, Spalding EP, Monshausen GB. The Receptor-like Kinase FERONIA Is Required for Mechanical Signal Transduction in Arabidopsis Seedlings. Curr Biol. 2014;24:1887-92.

34. Engelsdorf T, Hamann T. An update on receptor-like kinase involvement in the maintenance of plant cell wall integrity. Ann Bot. 2014;114:1339-47.

35. Liu Y, He C. Regulation of plant reactive oxygen species (ROS) in stress responses: learning from AtRBOHD. Plant Cell Rep. 2016;35:995-1007.

36. Windram O, Madhou P, McHattie S, Hill C, Hickman R, Cooke E, et al. Arabidopsis defense against Botrytis cinerea: chronology and regulation deciphered by high-resolution temporal transcriptomic analysis. Plant Cell. 2012;24:3530-57.

37. Gutierrez RA, Ewing RM, Cherry JM, Green PJ. Identification of unstable transcripts in Arabidopsis by cDNA microarray analysis: rapid decay is
associated with a group of touch- and specific clock-controlled genes. Proc Natl Acad Sci U S A. 2002;99:11513-8.

38. Widjaja I, Lassowskat I, Bethke G, Eschen-Lippold L, Long H-H, Naumann K, et al. A protein phosphatase 2C, responsive to the bacterial effector AvrRpm1 but not to the AvrB effector, regulates defense responses in Arabidopsis. Plant J. 2010;61:249-58.

39. Bjornson M, Benn G, Song X, Comai L, Franz AK, Dandekar AM, et al. Distinct roles for mitogen-activated protein kinase signaling and CALMODULIN-BINDING TRANSCRIPTIONAL ACTIVATOR3 in regulating the peak time and amplitude of the plant general stress response. Plant Physiol. 2014;166:988-96.

40. Chehab EW, Eich E, Braam J. Thigmomorphogenesis: a complex plant response to mechano-stimulation. J Exp Bot. 2009;60:43-56.

41. Huot B, Yao J, Montgomery BL, He SY. Growth-defense tradeoffs in plants: a balancing act to optimize fitness. Mol Plant. 2014;7:1267-87.

42. Weis BL, Kovacevic J, Missbach S, Schleiff E. Plant-Specific Features of Ribosome Biogenesis. Trends Plant Sci. 2015;20:729-40.

43. Grant EH, Fujino T, Beers EP, Brunner AM. Characterization of NAC domain transcription factors implicated in control of vascular cell differentiation in Arabidopsis and Populus. Planta. 2010;232:337-52.

44. Nakano Y, Yamaguchi M, Endo H, Rejab NA, Ohtani M. NAC-MYB-based transcriptional regulation of secondary cell wall biosynthesis in land plants. Front Plant Sci. 2015;6:288.

45. Stull RB. An Introduction to Boundary Layer Meteorology. Netherlands: Springer; 1988.

46. Rodriguez M, De Langre E, Moulia B. A scaling law for the effects of architecture and allometry on tree vibration modes suggests a biological tuning to modal compartmentalization. Am J Bot. 2008;95:1523-37.

47. Ding Y, Liu N, Virlouvet L, Riethoven J-J, Fromm M, Avramova Z. Four distinct types of dehydration stress memory genes in Arabidopsis thaliana. BMC Plant Biol. 2013;13:229.

48. Bruce TJA, Matthes MC, Napier JA, Pickett JA. Stressful "memories" of plants: Evidence and possible mechanisms. Plant Sci. 2007;173:603-8.

49. Murashige T, Skoog F. A revised medium for rapid growth and bioassays with tobacco tissue cultures. Physiol Plant. 1962;15:473-97.

50. Morizet J, Mingeau M, Ann Agron. 1976;27:183-205.

51. Affymetrix. GeneChip® Poplar Genome Array. 2005 [cited 2014 Oct 2]; Available from: http://media.affymetrix.com/support/technical/datasheets/ poplar datasheet.pdf

52. Irizarry RA, Hobbs B, Collin F, Beazer-Barclay YD, Antonellis KJ, Scherf U, et al. Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics. 2003;4:249-64.

53. Ge Y, Dudoit S, Speed T. Resampling-based multiple testing for microarray data analysis. Test. 2003;12:1-77.

54. Barrett T, Wilhite SE, Ledoux P, Evangelista C, Kim IF, Tomashevsky M, et al. NCBI GEO: Archive for functional genomics data sets - Update. Nucleic Acids Res. 2013;41:991-5.

55. Goodstein DM, Shu S, Howson R, Neupane R, Hayes RD, Fazo J, et al. Phytozome: a comparative platform for green plant genomics. Nucleic Acids Res. 2012;40:D1178-86.

56. Chang S, Puryear J, Cairney J. A simple and efficient method for isolating RNA from pine trees. Plant Mol Biol Report. 1993;11:113-6.

57. Pfaffl MW. A new mathematical model for relative quantification in real-time RT-PCR. Nucleic Acids Res. 2001;29, e45.

58. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene Ontology: tool for the unification of biology. Nat Genet. 2000;25:25-9.

Submit your next manuscript to BioMed Central
and we will help you at every step:

* We accept pre-submission inquiries

* Our selector tool helps you to find the most relevant journal

* We provide round the clock customer support

* Convenient online submission

* Thorough peer review

* Inclusion in PubMed and all major indexing services

* Maximum visibility for your research

Submit your manuscript at
www.biomedcentral.com/submit

BioMed Central