Léo Bouscarrat

Antoine Bonnefoy
email: antoine.bonnefoy@euranova.euleo.bouscarrat

Cécile Capponi

Carlos Ramisch
email: carlos.ramisch@lis-lab.fr

AMU-EURANOVA at CASE 2021 Task 1: Assessing the stability of multilingual BERT

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Event extraction is becoming more and more important as the number of online news increases. This task consists of extracting events from documents, especially news. An event is defined by a group of entities that give some information about the event. Therefore, the goal of this task is to extract, for each event, a group of entities that define the event, such as the place and time of the event.

This task is related but still different from named entity recognition (NER) as the issue is to group the entities that are related to the same event, and differentiate those related to different events. This difference makes the task harder and also complicates the annotation.

In the case of this shared task, the type of events to extract is protests (Hürriyetoglu et al., 2021a,b). This shared task is in the continuation of two previous shared tasks at CLEF 2019 [START_REF] Hürriyetoglu | Overview of clef 2019 lab protestnews: Extracting protests from news in a cross-context setting[END_REF] and AESPEN [START_REF] Hürriyetoglu | Automated extraction of socio-political events from news (AESPEN): Workshop and shared task report[END_REF]. The first one deals with English event extraction with three sub-tasks: document classification, sentence classification, and event information extraction. The second focuses on event sentence coreference identification, whose goal is to group sentences related to the same events. This year, task 1 is composed of the four aforementioned tasks and adds another difficulty: multilinguality. This year's data is available in English, Spanish, and Portuguese. Thus, it is important to note that there is much more data in English than in the other languages. For the document classification sub-task, to test multilingual capabilities, Hindi is available on the testing set only.

We have mainly focused on the last sub-task (event information extraction), but we have also submitted results for the first and second sub-tasks (document and sentence classification). We used multilingual BERT [START_REF] Devlin | Bert: Pre-training of deep bidirectional transformers for language understanding[END_REF], henceforth M-BERT, which is a model known to obtain near state-of-the-art results on many tasks. It is also supposed to work well for zero-or-few-shot learning on different languages [START_REF] Pires | How multilingual is multilingual bert?[END_REF]. We will see the results on these sub-tasks, especially for sub-task 4 where the training set available for Spanish and Portuguese is small. Thus, one of the issues with transformer-based models such as M-BERT is the instability on small datasets [START_REF] Dodge | Fine-tuning pretrained language models: Weight initializations, data orders, and early stopping[END_REF][START_REF] Ruder | Recent Advances in Language Model Fine-tuning[END_REF]. The instability issue is the fact that by changing some random seeds before the learning phase but using the same architecture, data and hyper-parameters the results can have a great variance. We will look at some solutions to mitigate this issue, and how this issue is impacting our results for sub-task 4. 1 2 Tasks and data Sub-tasks 1 and 2 can be seen as binary sequence classification, where the goal is to say if a given sequence is part of a specific class. In our case, a classifier must predict whether a document contains information about an event for sub-task 1 or if a sentence contains information about an event for sub-task 2.

Document and sentence classification tasks, subtasks 1 and 2, are not our main research interest. Moreover, the datasets provided for these tasks On the other hand, sub-task 4 not only has less training data available but also requires more finegrained token-based prediction. The goal of subtask 4 is to extract event information from snippets that contain sentences speaking about the same event. [START_REF] Hürriyetoglu | Overview of clef 2019 lab protestnews: Extracting protests from news in a cross-context setting[END_REF] have defined that an event has the following information classes (example in Figure 1):

• Time, which indicates when the protest took place,

• Facility name, which indicates in which facility the protest took place,

• Organizer, which indicates who organized the protest,

• Participant, which indicates who participated in the protest,

• Place, which indicates where the protest took place in a more general area than the facility (city, region, ...),

• Target, which indicates against whom or what the protest took place,

• Trigger, which is a specific word or group of words that indicate that a protest took place (examples: protested, attack, ...), Thus, not all the snippets contain all the classes, and they can contain several times the same classes. Each information can be composed of one or several adjacent words. Each snippet contains information related to one and only one event.

As the data is already separated into groups of sentences related to the same event, our approach consists of considering a task of named entity recognition with the aforementioned classes. Multilingual BERT has already been used for multilingual named entity recognition and showed great results compared to state-of-the-art models [START_REF] Hakala | Biomedical named entity recognition with multilingual BERT[END_REF].

The data is in BIO format [START_REF] Ramshaw | Text chunking using transformation-based learning[END_REF], where each word has a B tag or an I tag of a specific class or an O tag. The B tag means beginning and marks the beginning of a new entity. The tag I means inside, which has to be preceded by another I tag or a B tag, and marks that the word is inside an entity but not the first word of the entity. Finally, the O-tag means outside, which means the word is not part of an entity.

System overview

Our model is based on pre-trained multilingual BERT [START_REF] Devlin | Bert: Pre-training of deep bidirectional transformers for language understanding[END_REF]. This model has been pretrained on multilingual Wikipedia texts. To balance the fact that the data is not equally distributed between all the languages the authors used exponential smoothed weighting to under-sample the most present languages and over-sample the rarest ones. This does not perfectly balance all the languages but it reduces the impact of low-resourced languages.

The authors of the M-BERT paper shared the weights of a pretrained model that we use to do finetuning. Fine-tuning a model consists of taking an already trained model on a specific task and using this model as a starting point of the training for the task of interest. This approach has reached stateof-the-arts in numerous tasks. In the case of M-BERT, the pre-training tasks are Masked Language Modeling (MLM) and Next Sentence Prediction (NSP).

To be able to learn our task, we add a dense layer on top of the outputs of M-BERT and learn it during the fine-tuning. All our models are fine-tuning all the layers of M-BERT.

The implementation is the one from Hugging-Face's 'transformers' library [START_REF] Wolf | Transformers: State-of-the-art natural language processing[END_REF]. To train it on our data, the model is fine-tuned on each sub-task.

Sub task 1 and 2

For sub-tasks 1 and 2, we approach these tasks as binary sequence classification, as the goal is to predict whether or not a document (sub-task 1) or sentence (sub-task 2) contains relevant information about a protest event. Thus the size of the output of the dense layer is 2. We then perform an argmax on these values to predict a class. We use the base parameters in HuggingFace's 'transformers' library. The loss is a cross-entropy, the learning rate is handled by an AdamW optimizer [START_REF] Loshchilov | Decoupled weight decay regularization[END_REF] and the activation function is a gelu [START_REF] Hendrycks | Gaussian error linear units (gelus)[END_REF]. We use a dropout of 10% for the fully connected layers inside M-BERT and the attention probabilities.

One of the issues with M-BERT is the limited length of the input, as it can only take 512 tokens, which are tokenized words. M-BERT uses the wordpiece tokenizer [START_REF] Wu | Google's neural machine translation system: Bridging the gap between human and machine translation[END_REF]. A token is either a word if the tokenizer knows it, if it does not it will separate it into several sub-tokens which are known. For sub-task 1, as we are working with entire documents, it can be frequent that a document is longer than this limit and has to be broken down into several sub-documents. To retain contexts in each sub-documents we use an overlap of 150 tokens, which means between two sub-documents, they will have 150 tokens in common. Our method to output a class, in this case, is as follows:

• tokenize a document,

• if the tokenized document is longer than the 512-tokens limit, create different subdocuments with 150-tokens overlaps between each sub-document,

• generate a prediction for each sub-document,

• average all the predictions from subdocuments originated from the same document,

• take the argmax of the final prediction.

Sub-task 4

For sub-task 4, our approach is based on word classification where we predict a class for each word of the documents. One issue is that as words are tokenized and can be transformed into several sub-tokens we have to choose how to choose the prediction of a multitoken word. Our approach is to take the prediction of the first token composing a word as in [START_REF] Hakala | Biomedical named entity recognition with multilingual BERT[END_REF].

We also have to deal with the input size as some documents are longer than the limit. In this case, we separate them into sub-documents with an overlap of 150. Our approach is:

• tokenize a document,

• if the tokenized document is longer than the 512-tokens limit, create different subdocuments with 150-tokens overlaps between each sub-document,

• generate a prediction for each sub-document,

• reconstruct the entire document: take the first and second sub-documents, average the prediction for the same tokens (from the overlap), keep the prediction for the others, then use the same process with the obtained document and the next sub-document. As the size of each sequence is 512 and the overlap is only 150, no tokens can be in more than 2 different sequences,

• take the argmax of the final prediction for each word.

Soft macro-F1 loss

We used a soft macro-F1 loss [START_REF] Zachary C Lipton | Optimal thresholding of classifiers to maximize f1 measure[END_REF]. This loss is closer than categorical cross-entropy on BIO labels to the metric used to evaluate systems in the shared task. The main issue with F1 is its non-differentiability, so it cannot be used as is but must be modified to become differentiable. The F1 score is based on precision and recall, which in turn are functions of the number of true positives, false positives, and false negatives. These quantities are usually defined as follows:

tp = i∈tokens (pred(i) × true(i)) f p = i∈tokens (pred(i) × (1 -true(i))) f n = i∈tokens ((1 -pred(i)) × true(i))
With:

• tokens, the list of tokens in a document,

• true(i), 0 if the true label of the token i is of the negative class, 1 if the true label is of the positive class

• pred(i), 0 if the predicted label of the token i is of the negative class, 1 if the predicted label is of the positive class

As we use macro-F1 loss, we compute the F1 score for each class where the positive class is the current class and negative any other class, e.g. if the reference class is B-trigger, then true(i)=1 for B-trigger and true(i)=0 for all other classes when macro-averaging the F1.

We replace the binary function pred(i) by a function outputting the predicted probability of the token i to be of the positive class:

sof t tp = i∈tokens (proba(i) × true(i)) sof t f p = i∈tokens (proba(i) × (1 -true(i))) sof t f n = i∈tokens ((1 -proba(i)) × true(i))
With proba(i) outputting the probability of the token i to be of the positive class, this probability is the predicted probability resulting from the softmax activation of the fine-tuning network.

Then we compute, in a similar fashion as a normal F1, the precision and recall using the soft definitions of the true positive, false positive, and false negative. And finally we compute the F1 score with the given precision and recall. As a loss function is a criterion to be minimized whereas F1 is a score that we would like to maximize, the final loss is 1 -F 1.

Recommendation for improved stability

A known problem of Transformers-based models is the training instability, especially with small datasets [START_REF] Dodge | Fine-tuning pretrained language models: Weight initializations, data orders, and early stopping[END_REF][START_REF] Ruder | Recent Advances in Language Model Fine-tuning[END_REF]. [START_REF] Dodge | Fine-tuning pretrained language models: Weight initializations, data orders, and early stopping[END_REF] explain that two elements that have much influence on the stability are the data order and the initialization of the prediction layer, both controlled by pseudo-random numbers generated from a seed. To study the impact of these two elements on the models' stability, we freeze all the randomness on the other parts of the models and change only two different random seeds:

• the data order, i.e. the different batches and their order. Between two runs the model will see the same data during each epoch but the batches will be different, as the batches are built beforehand and do not change between epochs,

• the initialization of the linear layer used to predict the output of the model.

Another recommendation to work with Transformers-based models and small data made by [START_REF] Mosbach | On the stability of fine-tuning {bert}: Misconceptions, explanations, and strong baselines[END_REF] is to use smaller learning rates but compensating with more epochs. We have taken this into account during the hyper-parameter search.

Ruder (2021) recommend using behavioral finetuning to reduce fine-tuning instabilities. It is supposed to be especially helpful to have a better initialization of the final prediction layer. It has also already been used on named entity recognition tasks [START_REF] Broscheit | Investigating entity knowledge in BERT with simple neural end-to-end entity linking[END_REF] and has shown that it has improved results for a task with a very small training dataset. Thus, to do so, we need a task with the same number of classes, but much larger training datasets. As we did not find such a task, we decided to fine-tune our model on at least the different languages we are working with, English, Spanish and Portuguese. We used named entity recognition datasets and kept only three classes in common in all the datasets: person, organization, and location. These three types of entities can be found in the shared task.

To perform this test, the training has been done like that:

• the first fine-tuning is done on the concatenation of NER datasets in different languages, once the training is finished we save all the weights of the model,

• we load the weights of the previous model, except for the weights of the final prediction layer which are randomized with a given seed,

• we train the model on the dataset of the shared task.

4 Experimental setup

Data

The dataset of the shared task is based on articles from different newspapers in different languages. More information about this dataset can be found in (Hürriyetoglu et al., 2021a) For the final submissions of sub-tasks 1, 2, and 4 we divided the dataset given for training purposes into two parts with 80% for training and 20% for evaluation during the system training phase. We then predicted the data given for testing purposes during the shared task evaluation phase. The quantity of data for each sub-task and language can be found in Table 1. We can note that the majority of Sub-task English Spanish Portuguese Sub-task 1 9,324 1,000 the data is in English. Spanish and Portuguese are only a small part of the dataset. For all the experiments made on sub-task 4, we divided the dataset given for training purposes into three parts with 60% for training, 20% for evaluating and 20% for testing.

To be able to do our approach of behavioral finetuning, we needed some Named Entity Recognition datasets in English, Spanish and Portuguese. For English we used the CoNLL 2003 dataset [START_REF] Tjong | Introduction to the CoNLL-2003 shared task: Language-independent named entity recognition[END_REF], for Spanish the Spanish part of the CoNLL 2002 dataset [START_REF] Tjong | Introduction to the CoNLL-2002 shared task: Language-independent named entity recognition[END_REF] and for Portuguese the HAREM dataset [START_REF] Santos | Harem: An advanced ner evaluation contest for portuguese[END_REF]. Each of these datasets had already three different splits for training, development and test. Information about their size can be found in Table 2.

The dataset for Portuguese is pretty small compared to the two others, but the impact of the size can be interesting to study.

Hyper-parameter search

For sub-task 4, we did a hyper-parameter search to optimize the results. We used Ray Tune [START_REF] Liaw | Tune: A research platform for distributed model selection and training[END_REF] and the HyperOpt algorithm [START_REF] Bergstra | Making a science of model search: Hyperparameter optimization in hundreds of dimensions for vision architectures[END_REF]. We launched 30 different trainings, all the information about the search space and the hyper-parameters can be found in A.1. The goal is to optimize the macro-F1 on the evaluation set.

Our goal was to find a set of hyper-parameters that performs well to use always the same in the following experiments. We also wanted to evaluate the impacts of the hyper-parameters on the training.

Behavioral fine-tuning

For the first part of the behavioral fine-tuning, we trained an M-BERT model on the three NER datasets for one epoch. We only learn for one epoch for timing issues, as the learning on this datasets takes several hours. We then fine-tune the resulting models with the best set of hyper-parameters found with the hyper-parameter search.

Stability

To study the stability of the model and the impact of behavioral fine-tuning we made 6 sets of experiments with 20 experiments in each set:

• normal fine-tuning with random data order and frozen initialization of final layer,

• normal fine-tuning with frozen data order and random initialization of final layer,

• normal fine-tuning with random data order and random initialization of final layer,

• behavioral fine-tuning with random data order and frozen initialization of final layer,

• behavioral fine-tuning with frozen data order and random initialization of final layer,

• behavioral fine-tuning with random data order and random initialization of final layer, Once again it is important to note that what we called behavioral fine-tuning is different from behavioral fine-tuning as proposed by Ruder (2021), as we reset the final layer. Only the weights of all the layers of M-BERT are modified.

For each set of experiments we will look at the average of the macro-F1, as implemented in [START_REF] Nakayama | seqeval: A python framework for sequence labeling evaluation[END_REF], and the standard deviation of the macro-F1 on the training dataset, on the evaluation dataset, and on three different test datasets, one for each language. Thus we will be able to assess the importance of the instability, if our approach to behavioral fine-tuning helps to mitigate it and if it has similar results across the languages.

We can also note that in our implementation the batches are not randomized. They are built once before the learning phase and do not change, neither in content nor order of passage, between each epoch.

Results

Hyper-parameter search

The results of the hyper-parameter search can be seen in Figure 2. On the top pictures which represent the 30 experiments, we can see that a specific hyper-parameter seems to impact the worst results (in blue). This parameter is the learning rate, we can see it in the red box on the top image, all the blue lines are at the bottom, which means these experiments had a small learning rate. It seems that we obtain the best results with a learning rate around 5e-05 (0.00005), lower than 1e-06 seems to give bad results.

We can then focus on the bottom picture, with the same type of plot but with the worst results removed. Another hyper-parameter that seems to have an impact is the number of training epochs, 40 seems better than 20. We use a high number of epochs as recommended by [START_REF] Mosbach | On the stability of fine-tuning {bert}: Misconceptions, explanations, and strong baselines[END_REF] to limit the instability. Beyond the learning rate and number of epochs, it is then hard to find impactful hyper-parameters.

Finally, the set of hyper-parameters that has been selected is: For the stability experiments, the number of training epochs have been reduced to 20 for speed purposes. For the first part of the behavioral finetuning, the learning rate has been set to 1e-05 as more data were available.

Behavioral fine-tuning

The results on the test dataset of each model after one epoch of training can be found in Table 5.

We could not compare to state-of-the-art NER models on these three datasets as we do not take all the classes (classes such as MISC were removed before the learning phase). The metrics used on these datasets are not by classes, so the comparison cannot be made. However, the results are already much better than what a random classifier would output, thus the weights of the models should already be better than the weights of the base model.

Stability

The results of the different sets of experiments can be found in Table 3. First, we can see that the difference between behavioral fine-tuning and normal fine-tuning is not important enough to say one is better than the other. We can also note that the standard deviation is small for English, but not negligible for Spanish and Portuguese.

Final submission

The results of the final submissions can be found in Table 4. We can see that our results are lower than the best results, especially for sub-task 1 with a difference of between 30 to 50 macro-F1 score depending on the language, whereas for sub-tasks 2 and 4 the difference is close to 10 macro-F1 score for all the languages.

Conclusion

6.1 Sub-task 1 and 2

As we can see in Table 4, our final results for subtask 1 are much lower than the best results, but for sub-task 2 the difference is smaller. This is interesting as the tasks are pretty similar, thus expected the difference between our results and the best results to be of the same magnitude. One explanation could be our approach to handle documents longer than the input of M-BERT. We have chosen to take the average of the subdocuments, but if one part of a document contains an event the entire document does too. We may have better results looking if one sub-document at least is considered as having an event.

It is then hard to compare to other models as we have chosen to use one model for all the languages and we do not know the other approaches.

Sub-task 4

For sub-task 4 we have interesting results for all the languages, even for Spanish and Portuguese, as we were not sure that we could learn this task in a supervised fashion with the amount of data available. In a further study, we could compare our results with results obtained by fine-tuning monolingual models, where we fine-tune one model for each language with only the data of one language. This could show the impact of having data if using a multilingual model instead of several monolingual models improves or not the results. We do not expect good results for Spanish and Portuguese as the training dataset is pretty limited. The results seem to comfort the claim of [START_REF] Pires | How multilingual is multilingual bert?[END_REF] that M-BERT works well for few-shot learning on other languages.

The other question for sub-task 4 was about instability. In Table 3 we can see that the instability is way more pronounced for Spanish and Portuguese. It seems logical as we have fewer data available in Spanish and Portuguese than in English. The standard deviation for Spanish and Portuguese is large and can have a real impact on the final results. Finding good seeds could help to improve the results for Spanish and Portuguese.

Furthermore, our approach of behavioral finetuning did not help to reduce the instabilities. It was expected that one of the sources of the instability is the initialization of the prediction, and in our approach, the initialization of this layer is still random. In our approach, we only fine-tune the weights of M-BERT. This does not seem to work and reinforces the advice of [START_REF] Ruder | Recent Advances in Language Model Fine-tuning[END_REF] that using behavioral fine-tuning is more useful for having a good initialization of the final prediction layer.

On the two sources of randomness we studied, data order seems the most impactful for English, where we have more data. Nonetheless, for Spanish and Portuguese, the two sources have a large impact. In a further study, we could see how the quantity of data helps to decrease the impact of these sources of instabilities.

For the final submissions, the macro-F1 score for English and Portuguese is beneath the average macro-F1 score we found during our development phases. This could be due to bad seeds for randomness or because the splits are different. We did not try to find the best-performing seeds for the final submissions.

Figure 1 :

 1 Figure 1: Example of a snippet from sub-task 4.

Figure 2 :

 2 Figure 2: (Top) Parallel coordinates plot of the 30 experiments on sub-task 4 during the hyper-parameter search in function of the value of the hyper-parameters and the value of the F1 on the evaluation set. Each line represents an experiment, and each column a specific hyper-parameter, except the last which is the value of the metric. (Bottom) Same plot with the worst results removed to have a better view of the best results.

•

Table 2 :

 2 Number of elements for each dataset used in the behavioral fine-tuning in each split.

					1,487
	Sub-task 2 22,825	2,741	1,182
	Sub-task 4	808	33		30
	Table 1: Number of elements for each sub-task for each
	language in the data given for training purposes. Docu-
	ments for sub-task 1, sentences for sub-task 2, snippet
	(group of sentences about one event) for sub-task 4.
	Dataset	Train	Eval	Test
	CoNLL 2003 14,041 3,250 3,453
	CoNLL 2002 8,324 1,916 1,518
	HAREM	121		8	128

Table 3 :

 3 Average macro-F1 score, higher is better (standard deviation, lower is better) of the 20 experiments with the specified setup. N means normal fine-tuning and B behavioral fine-tuning. Data means data order and Init layer means initialization of the final layer. Rand means random, and fix refers to frozen.

	English	Spanish	Portuguese	Hindi
	Sub-task 1 53.46 (84.55) 46.47 (77.27) 46.47 (84.00) 29.66 (78.77)
	Sub-task 2 75.64 (85.32) 76.39 (88.61) 81.61 (88.47) /
	Sub-task 4 69.96 (78.11) 56.64 (66.20) 61.87 (73.24) /

Table 4 :

 4 Score of our final submissions for each sub-task, in parenthesis the score achieved by the best scoring team on each sub-task.

	Dataset	Test macro-F1
	CoNLL 2003	89.8
	CoNLL 2002	86.1
	HAREM	76.1

Table 5 :

 5 Macro-F1 score of the NER task on the test split of each dataset used in behavioral fine-tuning after training the base M-BERT for 1 epoch.

Acknowledgments

We thank Damien Fourrure, Arnaud Jacques, Guillaume Stempfel and our anonymous reviewers for their helpful comments.

A Appendix

A.1 Hyper-parameter search

The space search for our hyper-parameter search was:

• Number of training epochs: value in [20,25,30,40],

• Weight decay: uniform distribution between 0.001 and 1,

• Learning rate: value in [1e-5, 2e-5, 3e-5, 4e-5, 5e-5, 6e-5, 2e-7, 1e-7, 3e-7, 2e-8],

• Adafactor: value in "True", "False",

• Adam beta 1: uniform distribution between 0 and 1,

• Adam beta 2: uniform distribution between 0 and 1,

• Epsilon: value in [1e-8, 2e-8, 3e-8, 1e-9, 2e-9, 3e-10],

• Maximum gradient norm: uniform distribution between 0 and 1.

For the HyperOpt algorithm we used two set of hyper-parameters to help finding a good subspace. We maximized the macro-F1 on the evaluation dataset, and set the number of initial points before starting the algorithm to 5.


Loïc Francon
email: loic.francon@unige.ch

Christophe Corona

Irène Till-Bottraud

Philippe Choler

Erwan Roussel

Bradley Z Carlson

Samuel Morin

Brigitte Girard

Markus Stoffel

Shrub growth in the Alps diverges from air temperature since the 1990s

Keywords: dendroecology, alpine, Rhododendron ferrugineum, climate change, divergence, topography, shrubs

In the European Alps, air temperature has increased almost twice as much as the global average over the last century and, as a corollary, snow cover duration has decreased substantially. In the Arctic, dendroecological studies have evidenced that shrub growth is highly sensitive to temperature-this phenomenon has often been linked to shrub expansion and ecosystem greening. Yet, the impacts of climate change on mountain shrub radial growth have not been studied with a comparable level of detail so far. Moreover, dendroecological studies performed in mountain environments did not account for the potential modulation and/or buffering of global warming impacts by topography, despite its possible crucial role in complex alpine environments. To fill this gap, we analyzed a network of eight sites dominated by the dwarf shrub Rhododendron ferrugineum. The sites selected for analysis represent the diversity of continentality, elevation and slope aspect that can be found across the French Alps. We quantified annual radial increment growth for 119 individuals, assembled meteorological reanalyzes specifically accounting for topographic effects (elevation, slope and aspect) and assessed climate-growth relations using a mixed modeling approach. In agreement with a vast majority of dendroecological work conducted in alpine and arctic environments, we find that the number of growing degree days during the snow-free period snow-free growing degree days (SFGDDs) is a strong and consistent driver of R. ferrugineum growth across all sites since 1960 until the late 1980s. We also document a marked loss of sensitivity of radial growth to increasing SFGDD since the 1990s, with this decoupling being more pronounced at the driest sites. Our observations of the spatial and temporal variability of shrub sensitivity to limiting factors can be compared to the 'divergence' problem observed in tree-ring series from circumpolar and alpine regions and, accordingly, sheds light on possible future trajectories of alpine shrub growth in response to ongoing climate change.

Introduction

Mountainous and Arctic regions are experiencing a temperature increase that is about two times larger than the global average; as a corollary, observations also show a general decline in snow cover [START_REF] Holmes | Tree growth responses to changing temperatures across space and time: a fine-scale analysis at[END_REF]). In the Arctic, higher temperatures during summer lead to intense greening and increased shrub productivity and cover, which ultimately drive climate feedbacks (Myers-Smith et al 2015a, [START_REF] Berner | Summer warming explains widespread but not uniform greening in the Arctic tundra biome[END_REF]. Likewise, the ongoing shrub expansion has been shown to modify a range of ecosystem processes (Myers-Smith et al 2011) by altering surface albedo, snowpack depth, energy, water balance, or permafrost conditions (Sturm et al 2001, Liston et al 2002[START_REF] Chapin | Role of land-surface changes in Arctic summer warming[END_REF]. By analyzing a network of circumpolar shrub rings, Myers-Smith et al (2015a) suggested that increasing summer temperature have enhanced shrub radial growth. Several arctic studies also showed that topography may have a positive buffering effect on changing summer temperatures by modifying soil moisture conditions (Myers-Smith et al 2015a). Ropars et al (2015) and [START_REF] Ackerman | Arctic shrub growth trajectories differ across soil moisture levels Glob[END_REF] thereby demonstrated that the response of shrub individuals to climatic parameters can be strongly modulated by local topography. However, [START_REF] Bär | Growth-ring variations of dwarf shrubs reflect regional climate signals in alpine environments rather than topoclimatic differences[END_REF] and [START_REF] Ackerman | Uniform shrub growth response to June temperature across the North Slope of Alaska[END_REF] did not detect any inconsistent responses of shrubs sampled in contrasting topographic situations to ongoing climate warming.

In the European Alps, increasing vegetation productivity and denser shrub cover above treeline have been documented through remote sensing or resurveys of vegetation plots [START_REF] L20703 Dullinger | Patterns of Shrub invasion into high mountain grasslands of the Northern Calcareous Alps, Austria Arctic[END_REF][START_REF] Cannone | Unexpected impacts of climate change on alpine vegetation Front[END_REF][START_REF] Carlson | Observed long-term greening of alpine vegetation-a case study in the French Alps[END_REF][START_REF] Filippa | Climatic drivers of greening trends in the Alps Remote Sens[END_REF], Malfasi and Cannone 2020). Dendroecological studies have highlighted the negative role of winter precipitation and subsequent, delayed melt-out dates on radial growth of shrubs, but have also demonstrated the positive effect that increasing summer temperatures can have on shrub growth (Pellizzari et al 2014[START_REF] Francon | Warm summers and moderate winter precipitation boost Rhododendron ferrugineum L. growth in the Taillefer massif (French Alps)[END_REF]. In the Italian Alps, a study suggested a positive, spatially consistent and stable effect of increasing growing season length on shrub radial growth (Juniperus communis nana; [START_REF] Carrer | Winter precipitation-not summer temperature-is still the main driver for Alpine shrub growth[END_REF]. By contrast, two recent, albeit spatially-limited studies emphasized that elevation and slope aspect can significantly modulate Rhododendron ferrugineum growth response to increasing temperature significantly (Francon et al 2020a(Francon et al , 2020b)). Slope aspect, as a local topographic factor, can mimic temperature differences of large elevational or latitudinal gradients over distances of just a few meters (Scherrer and Körner 2009) by modulating solar incidence angles, frost exposure and soil moisture conditions [START_REF] Aalto | Revealing topoclimatic heterogeneity using meteorological station data: revealing topoclimatic heterogeneity[END_REF]. Additionally, within mountain ranges, rain-shadow effects-controlling precipitation sums and ariditywere demonstrated to be of paramount importance for plant interactions and growth (Michalet et al 2003). However, no dendroecological study has so far specifically investigated the impacts of topoclimatic variability (defined as the local scale of a few hectares, where climate is driven by variations in topography; Aalto et al 2017) on shrub growth in mountain regions. This research gap is critical given that (a) topoclimatic heterogeneity in alpine environments [START_REF] Dobrowski | A climatic basis for microrefugia: the influence of terrain on climate: a climatic basis for microrefugia Glob[END_REF]) is likely to buffer the response of shrubs to global warming; and that (b) strong moisture gradients that exist over short distances due to the presence of rugged terrain (Scherrer and Körner 2011) could potentially modulate the productivity and performance of alpine vegetation more markedly than in the Arctic [START_REF] Ernakovich | Predicted responses of arctic and alpine ecosystems to altered seasonality under climate change Glob[END_REF]. Recent shrub growth and expansion in the European Alps is a critical and understudied topic, considering the various ecosystem 'disservices' expected to be conferred by the replacement of alpine meadows by dense shrub canopies, including: modifications of the microclimate (Mekonnen et al 2021) localized reductions in vascular plant diversity [START_REF] Anthelme | Shrub encroachment in the Alps gives rise to the convergence of sub-alpine communities on a regional scale[END_REF][START_REF] Boscutti | Shrub growth and plant diversity along an elevation gradient: evidence of indirect effects of climate on alpine ecosystems[END_REF] diminished pasture resource and quality (Klein et al 2007), and altered soil chemical composition (Pornon et al 1996).

This paper assesses the impacts of global warming on shrub growth for a dominant dwarf shrub species in alpine heathlands, R. ferrugineum (L.), with a focus on the possible mediating role of topoclimatic conditions could have on shrub growth responses to climate change. To this end, we (a) assembled annual growth records from a network of eight sites designed to be representative of the contrasts in continentality, elevation and slope aspect across the French Alps. Analyses were performed (b) with SAFRAN-Crocus (Durand et al 2009a, 2009b, Vionnet et al 2012, Vernay et al 2019), a meteorological and snow cover reanalysis providing an estimation of meteorological and snow cover variables for different slope, aspects and per 300 m altitude steps within mountainous regions called 'massifs' . In a last analytical step, we (c) combined bootstrap correlation functions with mixed effects models applied to single R. ferrugineum shrubs to quantify threshold effects and shifts in shrub-climate relations resulting from ongoing global warming in the French Alps.

Material and methods

Shrub-ring chronology development

R. ferrugineum is an evergreen dwarf shrub species with highly branched trailing that can reach up to 0.8 m in height. It has a wide geographic range of occurrence across the European Alps (Ozenda 1985) and is in expansion in the treeline ecotone (Pornon et al 1996). Complementary sexual and vegetative reproductive strategies [START_REF] Doche | Changes of humus features along with a successional gradient of Rhododendron ferrugineum (L.) Populations (Subalpine Level, Northwestern Alps, France) Arctic, Antarct[END_REF]. This strategy also allows R. ferrugineum to outcompete other plants in subalpine heaths (Theurillat and Schlüssel 2000), to accumulate large amounts of biomass per unit area and to form large heathlands that reach up to 90%-100% cover after 150-250 years (Pornon et al 1996, Pornon andDoche 1996)explain its abundance in the subalpine belt. Given its (a) wide distribution from c. 1600-2500 m asl across the European Alps, (b) its longevity, as well as (c) the formation of clearly identifiable annual growth rings, R. ferrugineum is a reliable model to study shrub response to long-term climate variability [START_REF] Francon | Warm summers and moderate winter precipitation boost Rhododendron ferrugineum L. growth in the Taillefer massif (French Alps)[END_REF]. We thus sampled R. ferrugineum 1. Umbro-thermal diagrams (Chartreuse, Oisans, Mont-Blanc and Queyras) were obtained from the Safran-CROCUS (S2M) reanalysis (Météo France, 1981[START_REF] Bonet | Modelling the production and species richness of wild mushrooms in pine forests of the Central Pyrenees in northeastern Spain Can[END_REF].

Table 1. Characteristics of the eight study sites with the name of the mountain massif (Massif), elevation of plots above sea level (Elev.), index of continentality (α), and aspect (Asp.) as well as growth ring chronology characteristics and period covered (Length), number of individuals used (Nb.), signal strength (expressed population signal or EPS), mean intercorrelation between individual series inter-series correlation (rbar) and autocorrelation of order 1 (AC) included in the dendroecological network. The EPS, rbar and AC values have been computed for detrended chronologies. Specimens of R. ferrugineum were sampled in heathlands or large patches at all sites except at CN2400 where individuals were growing isolated. When sampling, a minimum distance of 4 m was maintained between each sample selected for analysis so as to avoid replication of samples from the same individual [START_REF] Escaravage | Clonal diversity in a Rhododendron ferrugineum L. (Ericaceae) population inferred from AFLP markers[END_REF]. For each individual, one to three sections were taken from one or two stems to allow for a serial sectioning and repeat measurements of growth-ring series on the same individual (figure S2 (available online at stacks.iop.org/ERL/16/074026/mmedia)), but also to detect missing rings in growth-ring series as shrubs often produce very narrow, wedging growth rings (Kolishchuk 1990). From each section, a microsection was obtained with a sliding microtome; each section was then stained, colored with safranin and astra blue dye to enhance ring boundaries and finally mounted on slides (Schweingruber and Poschold 2005).

Cross-dating (i.e. the assignment of exact calendar years to the formation of each ring) was realized on the digitalized images of thin sections with CooRecorder and CDendro (CYBIS Elektonik & Data AB); this included a three-step procedure relying on the visual synchronicity between (a) elementary ringwidth (RW) series measured at three radii of the same cross section, (b) mean RW series obtained for each section of the same individual, and finally (c) between individual chronologies at a given site. Visual cross-dating was validated statistically with COFECHA [START_REF] Holmes | Tree growth responses to changing temperatures across space and time: a fine-scale analysis at[END_REF]. In a subsequent step, we detrended each section chronology using a negative exponential function to eliminate non-climatic biological trends using ARSTAN [START_REF] Cook | The decomposition of tree-ring series for environmental studies Tree Ring[END_REF]. After detrending, ring-width indices (RWIs) were averaged for each cross-dated section and individual with a bi-weighted robust mean; this step aimed at reducing the influence of outliers and at developing mean detrended site chronologies [START_REF] Cook | The smoothing spline: a new approach to standardizing forest interior tree-ring width series for dendroclimatic studies Tree Ring[END_REF]. Finally, the eight detrended site chronologies were compared with each other using a correlation matrix to check for radial growth consistency between the sites.

We computed descriptive statistics including standard deviation, first-order AC, mean rbar and the running 35 years EPS for each detrended site chronology with ARSTAN (table 1), whereby EPS quantifies the strength of the common climate signal in growth-ring proxies (Wigley et al 1984). We used a quality threshold EPS ⩾ 0.8 to evaluate reliability of our chronologies. Sample depth is shown on figure S2.

Meteorological time series and topoclimatic characteristics of network sites

We extracted series of hourly temperature, hourly precipitation totals and daily snowpack depths at each of the 8 network sites from SAFRAN-Crocus (S2M) reanalyzes for the period 1959-2017 (Durand et al 2009a, 2009b, Vionnet et al 2012, Vernay et al 2019). They provide continuous time series of meteorological variables at hourly resolution, for different elevation bands, slope aspects and angles within massifs, i.e. horizontally climatologically homogeneous regions. Furthermore, they also provide snow height and snow water equivalent values. Note however, that vegetation is represented in a simplified way in the simulations, assuming continuous grass cover to simulate typical snow conditions and underlying ground temperature for open areas (Vionnet et al 2012). Based on these reanalyzes, we derived snow melt-out dates at our study sites according to the elevation, slope aspect and angle that had been measured in the field for each plot. Melt-out dates correspond to the date when snow cover depth reaches 10 cm for the last time within a snow cover period of at least 7 d. From the date of melt-out to August 31 (considered here as the last day of growth, based on multiple dendrometric series, unpublished data), we summed all above-zero daily mean temperatures to obtain snow-free growing degree days (SFGDDs). SFGDD summarize thermal conditions experienced during the period during which growth occurs. We also generated monthly and seasonal mean temperature series and precipitation totals. To account for the potential impact of summer drought on radial growth, we computed series of standardized precipitationevapotranspiration indices (or SPEI, Vicente-Serrano et al 2010) for two month periods corresponding to the length of the growing season. Trends in summer precipitation and temperature, SFGDD and meltout dates are shown in the supplementary materials attached to this manuscript (figure S3). Topographic characteristics (elevation and slope aspect) and the Gams continentality index (α, Gams 1932) were employed as additional variables to characterize each site, whereas the Northness index-ranging from-1 for a southern to 1 for a northern orientation-was given at each site as the cosine of slope aspect (α) and is used here as an indicator of rain-shadow effects while controlling for the effect of elevation. α was calculated following the equation: cotg (α) = P/A where P and A represent the mean annual precipitation sum (mm) and the elevation (m asl), respectively, at each site (table 1). Elevation and slope aspect are known to affect temperature, light availability, solar incidence angles and soil conditions in mountainous terrain (Scherrer and Körner 2009[START_REF] Bonet | Modelling the production and species richness of wild mushrooms in pine forests of the Central Pyrenees in northeastern Spain Can[END_REF][START_REF] Aalto | Revealing topoclimatic heterogeneity using meteorological station data: revealing topoclimatic heterogeneity[END_REF].

Relationships between shrub radial growth and climate

Climate-growth analyses were performed in two steps: we first calculated bivariate correlation functions (BCFs) between the detrended site chronologies (RWI) and time series of monthly (and summer) temperature, precipitation, SPEI, melt-out dates and SFGDD. This preliminary analysis was intended to detect the main drivers of radial growth at each site. BCFs were computed separately for the periods 1960-1988 and 1989-2017 so as to detect potential shifts in climate-growth relationships related to global warming. The two periods of equal length cover almost 30 years each, in accordance with the definition of a climate normal. The second period also reflects the marked warming observed since the late 1980s (Lo and Hsu 2010, Reid et al 2016). This regime shift to warmer temperatures had a particularly significant impact on biophysical environments including the physiology, development, behavior and phenology of living organisms (Oberhuber et al 2008, Reid et al 2016). Subsequently, relationships between SFGDD and log-transformed scaled raw RW of R. ferrugineum were investigated further using linear mixed models. The advantage of these models relies in (a) their ability to manage intra-site variance (Myers-Smith et al 2015b); (b) the possibility to compare shrub response between different topoclimatic conditions, thereby allowing to account for potential interactions between variables, for example topography and climate (Myers-Smith et al 2015b); (c) the possibility to use raw ring measurements, thus avoiding detrending that could remove meaningful ecological information from growth RW series [START_REF] Ackerman | Uniform shrub growth response to June temperature across the North Slope of Alaska[END_REF]; (d) the detection of possible non-linear relationships and threshold effects in climate-growth relations by introducing, for example, quadratic effects; and in (e) the inference of predictions of general trends for sites that have not been sampled by using data of multiple sites.

Linear mixed models computed using the nlme package (Pinheiro and Bates 2000) in R (R Core Team 2016) are an extension of simple linear models to allow both fixed and random effects. The inclusion of random effects enables control of pseudo-replication, i.e. of the degree of dependence between observations within the same group. In our case, we chose multiple nested grouping levels for the random intercept structure of models: section nested within individual as well as individual within site as random terms [START_REF] Harrison | A brief introduction to mixed effects modelling and multi-model inference in ecology Peer J. 6 e4794 Hoch G[END_REF]. Fixed effects, equivalent to explanatory variables in standard linear regressions, are variables that we expect to affect RWs. These include cambial age, elevation, Gams continentality index, northness index and mean annual SFGDD (normalized and scaled). Cambial age (i.e. number of years) was introduced in the model to account for age-related trends in growth-ring series. Particular attention was also paid to interactions between SFGDD and other geographical variables (elevation, Gams continentality index and northness index) as they could modulate shrub response to increasing air temperature.

To detect potential threshold effects, we also explored second-order (non-linear) relationships between SFGDD and radial growth in the models. We therefore log-transformed and scaled mean raw RW measurements at the level of individual sections (hereafter referred to as RW) to obtain the required normality of residuals. An AC structure (AR1, firstorder autoregressive process) was included as well so as to eliminate temporal biological memory from growth-ring series (Weijers et al 2018b). The predictive capacity of a single independent variable was then determined by the absolute value of Student's t statistics for its linear or quadratic effect on RW. We calculated a pseudo-R 2 of the selected models comprising marginal (R 2 m) and conditional (R 2 c) values (Nakagawa and Schielzeth 2013). R 2 m and R 2 c account for the proportions of the variance explained by the fixed factors and by the whole model (i.e. fixed plus random factors), respectively. R 2 m and R 2 c values were calculated with the MuMIn package [START_REF] Barton | Package 'MuMIn' Model selection and model averaging based on information criteria[END_REF]. Elimination of non-statistically significant terms (at p ⩽ 0.05) was then used to select the most parsimonious model. The Akaike Information Criterion (AIC) and the Bayesian Information Criterionimposing a stronger penalty for increasing model complexity-were then used to compare the resulting models. We considered the overall climate sensitivity to be the comparison in AIC value between the best model and a null model (Myers-Smith et al 2015a). Finally, multicollinearity issues where checked based on the variance inflation factor with a threshold set at three (Zuur et al 2010) in car package available in R. To detect potential impacts of global warming on climate-growth relationships, three models were implemented for the periods 1960-88, 1989-2017 and 1960-2017.

Results

Figure 2 provides a synthesis of changes in SFGDDs and RWI values at the eight sites sampled across the French Alps. The SFGDD series do not exhibit any significant trend from 1960 until the late 1980s. By contrast, and starting in the early 1990s, a clear shift of SFGDD is observed at each site along with an increasing trend in all SFGDD series over the period 1989-2017. No comparable trend is observed in the RWI chronologies over the same period (figure 2).

Using the growth RW series obtained from the 225 cross-sections sampled from 119 individuals of R. ferrugineum, we constructed eight detrended site chronologies (see table 1 and figures S1 and S2 for their characteristics) that were subsequently used for the analysis of climate-growth relationships. All chronologies exceed the 0.8 EPS threshold over the 1960-2017 period with the exception of SAL1800 (1963-2017) (table 1). Over the period 1960-89, BCFs (see figure S4) computed between meteorological variables and detrended chronologies reveal strong (r values ranging between 0.39 and 0.75, figure 2) and statistically significant correlations (p < 0.05) between RWIs and SFGDD at all sites except at CHAR (r = 0.22, p = 0.24, figures 2 and S4). SFGDD was the main climatic driver of radial growth at SAL1800, CN2400, LORIAZ, QGL and QN. At SAL2000, comparable correlations (r = 0.46, p < 0.05) were computed between the shrub-ring chronology, SFGDD and summer temperature while at QW, correlation with SFGDD ranks second after the correlation with February precipitation. More generally, RWI correlates negatively with winter precipitation at a vast majority of sites (except for SAL1800 and SAL2000). Finally, mean late summer and early fall temperature in the year preceding growth-ring formation (September-October) correlated positively with RWI at SAL2000, CN2400, LORIAZ and QGL (figure S4).

Over the period 1989-2017, site-to-site climategrowth correlation profiles (figure S4) differ much more than over the period 1960-88. Interestingly, SFGDD remains the single most important driver of radial growth only at CN2400 (r = 0.42, p < 0.05, figure 2), whereas correlations between SFGDD and radial growth dropped close to zero at CHAR, LORIAZ, QGL, and SAL1800-and even become significantly negative at QW and QN (figure 2). Moreover, we find significant negative correlations between shrub chronologies and August air temperatures at CHAR, QN, QGL, and QW (figure S4). At QW, a negative impact of drought on R. ferrugineum radial growth is also evidenced by significant correlations between SFGDD (r = -0.29, p < 0.05) and summer SPEI (r = 0.55, p < 0.001). By contrast, detrended RW chronologies from SAL1800, SAL2000 and QGL are positively correlated with summer and early fall (September-October) temperatures. Finally, winter precipitation remains an important driver of radial growth at CN2400, LORIAZ, QGL and QN (figure S4).

The correlation matrices computed between site chronologies confirm the loss of consistency between both time periods (figure 3) and indicates that median intercorrelation decreased from 0.36 to 0.13 from 1960-88 to 1989-2017, respectively. Whereas all chronologies-with exception of LORIAZ-were significantly and positively correlated between sites until 1988, we observe a shift to low and mostly statistically insignificant values for a majority of correlation coefficients over the period 1989-2017.

To further explore potential drivers of shrubs growth, we used a linear mixed modeling approach to associate RW with SFGDD and a suite of geographical variables (i.e. elevation, Gams continentality index and northness index). Prior to computing mixed models, we investigated relationships between RW and SFGDD (figure S5). Conditional R 2 values reveals that, regardless of the period, all models can explain ∼50% of the interannual growth variability (44%, 45%, and 47% for 1960-2017, 1960-88, and 1989-2017, respectively; table 2). Marginal R 2 values (fixed factors) explain approximately half of the conditional R 2 showing that the proportions of variance explained by fixed and random variables remain comparable, regardless of the time period considered. According to the ∆AIC null vs model values, the models explain a greater part of the interannual growth variability than the null model. ∆AIC null vs model values yet are much higher for 1960-88 than for 1989-2017 (with values of 701.8 and 151.8, respectively).

For the period 1960-88, the model yields a positive, highly significant linear relation between and meteorological variables already encapsulate a great dependence on elevation.

R 2 c = 0.44: R 2 m = 0.23 R 2 c = 0.45; R 2 m =0.19 R 2 c = 0.47; R 2 m = 0.

Discussion

SFGDDs as the main driver of R. ferrugineum radial growth before 1989

Our results reveal strong, highly significant and consistent relationships between SFGDD and shrub radial growth until 1988 for all eight sites considered in the network. Our findings are consistent with the widely accepted concept that temperature-and more particularly growing degree days (Jochner et al Furthermore, previous studies considered sensitivity to temperature was considered to be consistent over contrasted topographic contexts in most studies [START_REF] Bär | Growth-ring variations of dwarf shrubs reflect regional climate signals in alpine environments rather than topoclimatic differences[END_REF][START_REF] Ackerman | Uniform shrub growth response to June temperature across the North Slope of Alaska[END_REF], Carrer et al 2019).

Temperatures and RWs diverge since 1989

Interestingly, the increasing number of growing degree days observed in our meteorological series after 1989 did not result in a release of radial growth rates from thermal limitations. On the contrary, our results suggest a divergence between increasing SFGDD and shrub radial growth and thus, point to non-linear relationships between the two variables.

In other words, and in agreement with smaller scale studies carried out in the Taillefer and Queyras massifs (Francon et al 2020a(Francon et al , 2020b)), our data suggest a loss of sensitivity of shrub growth to SFGDD beyond a certain threshold, and that this threshold is indeed modulated by continentality and slope aspect. We hypothesize that the observed non-linear response of shrub growth is affected by a combination of two effects. On the one hand, the strong increase in growing season length makes shrubs less constrained by the availability of energy to complete their life cycle (Francon et al 2020a). On the other hand, however, this beneficial effect of rising temperatures on growth also seems to be compromised by the negative effects of drought, with the latter depending strongly on the position of the shrub along moisture gradients (Francon et al 2020b). Such non-linearity was not observed by [START_REF] Carrer | Winter precipitation-not summer temperature-is still the main driver for Alpine shrub growth[END_REF] who reported a consistent and negative role of winter precipitation on J. communis nana growth in the Italian Alps, despite the fact that a persistent and significant warming trend could be observed at their site as well. and high-elevation sites [START_REF] Büntgen | Testing for tree-ring divergence in the European Alps Glob[END_REF]. This body of literature suggests-although with a focus on trees-complex and non-linear growth responses of trees to climate changes, thereby leading to a decrease of year-to-year sensitivity of tree growth in previously temperature-limited environments since the mid to late 20th century [START_REF] D'arrigo | Thresholds for warming-induced growth decline at elevational tree line in the Yukon Territory, Canada: thresholds for warminginduced growth decline Glob. Biogeochem[END_REF], 2008, Wilmking 2005[START_REF] Büntgen | Testing for tree-ring divergence in the European Alps Glob[END_REF], Oberhuber et al 2008, Leonelli et al 2009). Potential causes for this divergence include warming-induced thresholds of tree growth (D'Arrigo et al 2004) or increasing drought stress [START_REF] Büntgen | Summer temperature variations in the European Alps[END_REF]. Assessment of divergence phenomena for shrubs has been only mentioned once at the Arctic tundra biome in response to warming and drying climate [START_REF] Buchwal | Divergence of Arctic shrub growth associated with sea ice decline[END_REF]. In the same manner as for R. ferrugineum, Arctic shrubs display a pronounced growth response heterogeneity starting in the mid-1990s. However, a careful assessment of a divergence problem for alpine shrubs would require shrub ring chronologies covering longer periods, but also broader regions and accounting for other species.

Topographically driven moisture limitation as a possible driver of divergence

The fact that lower correlations are observed between the chronologies at all eight sites after 1989 and that greater heterogeneity emerges in the BCF profiles points to strong site effects that have become the key modulators of shrub growth responses to climate over the last few decades. Thus, we hypothesize that the recent loss of sensitivity and the non-linearity in the response of shrub growth to SFGDD, also referred to the 'divergence problem' in forestry and tree-ring research, could be explained by moisture limitations driven by topoclimate and microsite (Wilmking et al 2004, Wilmking 2005). The influence of topoclimate is confirmed by our mixed modeling approach which includes continentality and slope aspect variables in interaction with SFGDD as significant terms. Indeed, both variables directly control soil moisture both at the mesoscale (slope aspect) and at regional (continentality) scales. Therefore, shrubs growing on drier sites (characterized by higher Gams angles and lower northness indices) are more affected by the ongoing increase of air temperatures. Based on dendroecological studies realized in the Arctic, Myers-Smith et al (2015a) likewise postulated that topography, soil properties and snow-related moisture availability would modulate the positive effect of rising temperature on shrub growth. In other words, the sensitivity of shrubs to moisture stress has increased clearly over the last few decades, but this evolution still greatly depends on site characteristics. We are not aware of other studies reporting on topographical effects of moisture stress on shrub growth. In the Italian Alps, the absence of topographic effects on shrub growth reported by [START_REF] Carrer | Winter precipitation-not summer temperature-is still the main driver for Alpine shrub growth[END_REF] for J. communis nana can be explained by the its lesser sensitivity to drought as compared to R. ferrugineum [START_REF] Gracia | Overstory structure and topographic gradients determining diversity and abundance of understory shrub species in temperate forests in central Pyrenees (NE Spain) For[END_REF] or to less pronounced topographic variability within the Italian site network. Besides moisture limitation, an increasing frequency of late frost events has been shown to cause both short-term (through damage to aerial parts of the shrub, [START_REF] Bokhorst | Winter warming events damage sub-Arctic vegetation: consistent evidence from an experimental manipulation and a natural event[END_REF], Treharne et al 2019) and longer-term (through winter droughts, Phoenix and Bjerke 2016, Charrier et al 2017) stresses to shrub vegetation and could, thus, further exacerbate the observed divergence in growth responses. However, when comparing the RWI with the freezing degree days index (calculated as the sum of daily minimum temperature below-2 • C estimated from the moment of melt-out to the end of August), we do not detect any significant effects By contrast, at drier sites (bottom left panel), a non-linear response is anticipated and will likely limit radial growth above a certain SFGDD threshold.

of late frost on R. ferrugineum growth in our dataset. This absence of frost effects should yet be interpreted with caution as frost risk depends strongly on plant exposure and vulnerability, which in turn will be determined by microtopography, plant phenology and hardening [START_REF] Bréda | Vulnerability to forest decline in a context of climate changes: new prospects about an old question in forest ecology[END_REF]. On this specific point, studies coupling high-resolution microclimatic and ecophysiological monitoring are yet to be realized to allow for more robust conclusions about effects of frost on shrub growth [START_REF] Charrier | Monitoring of freezing dynamics in trees: a simple phase shift causes complexity[END_REF]. By contrast, the linear mixed model approach allowed determination, for the first time, of thresholds for the occurrence of divergence. Accordingly, values of 900 • C and 1600 • C represent SFGDD thresholds for R. ferrugineum radial growth in the dry (or eastern) and moist (or western) mountain massifs of the French Alps. The threshold also varies with changing aspect and could be set to 900

Conclusion

The approach developed in this study relies on a network of R. ferrugineum RW series from eight contrasted topoclimatic locations across the French Alps. The unprecedented size and spatial coverage of the dataset, along with the temporal resolution of the meteorological and snow cover reanalyzes data used allowed facilitated the disentanglement of variables driving shrub growth at the individual level across sites. We demonstrated that a clear divergence emerged between shrub growth across the French Alps and increasing air temperatures after the early 1990s, a phenomenon that has previously not been demonstrated in the alpine environment [START_REF] Carrer | Winter precipitation-not summer temperature-is still the main driver for Alpine shrub growth[END_REF], but has been detected in shrub rings in the arctic environment [START_REF] Buchwal | Divergence of Arctic shrub growth associated with sea ice decline[END_REF]. Our study also points to a likely modulation of the observed divergence by temperature thresholds, which in turn are controlled by differences in moisture availability resulting from contrasted topographical and geographical parameters. Given the increasing air temperatures and the earlier snow melt-out timing projected by climate model simulations [START_REF] Gobiet | 21st century climate change in the European Alps-a review[END_REF][START_REF] Beniston | The European mountain cryosphere: a review of its current state, trends, and future challenges[END_REF][START_REF] Holmes | Tree growth responses to changing temperatures across space and time: a fine-scale analysis at[END_REF]), these thresholds will likely be exceeded even more frequently during the 21st century. By synthesizing the findings of this study, we offer a conceptual model (figure 5) summarizing the potential response of R. ferrugineum radial growth to increasing SFGDD in a warmer climate for wetter and moisture-limited sites. Individuals should benefit from increasing temperatures on north-facing slopes in the north-westernmost regions of the French Alps. By contrast, we expect their productivity to decrease as a result of increasingly severe moisture limitations at the drier sites in the Southern Alps. Collectively, our results underline the high sensitivity of R. ferrugineum to climate change. We therefore consider this shrub to be a valuable sentinel species as it integrates the complex and multi-directional effects that climate change will continue to have on alpine ecosystems (Francon et al 2020a). Scaling up, we encourage future studies to investigate potential relationships between divergent growth responses demonstrated by dendroecological analyses and the greening or browning trends quantified with remote sensing methods [START_REF] Cunliffe | Aboveground biomass corresponds strongly with drone-derived canopy height but weakly with greenness (NDVI) in a shrub tundra landscape[END_REF]. At a more micro-scale, further work is also need to disentangle the underlying eco-physiological processes for the non-linear relationships observed between ring growth and temperature increase (de Swaef et al 2015, de Micco et al 2019). The latter should necessarily be complemented with high-resolution spatial imagery (Myers-Smith et al 2020) to account for the heterogeneity of alpine environments. We encourage future alpine studies to focus on the implications of the bidirectional evolution of shrub growth evidenced here on shrub expansion, carbon transfers, radiative feedbacks, and biodiversity.

Figure 1 .

 1 Figure 1. Location of the eight study sites of the R. ferrugineum network used in this study: CHAR refers to the site sampled in the Chartreuse massif; CN2400, SAL2000, SAL1800 to three sites with different slope aspects in the Taillefer massif at 2400, 2000 and 1800 m asl; LORIAZ to the site sampled in the Mont-Blanc massif; QGL, QN, QW refer to three sites with different slope aspects sampled at elevations ranging from 2350 to 2550 m asl in the Queyras massif. Details on site characteristics are given in table1. Umbro-thermal diagrams (Chartreuse, Oisans, Mont-Blanc and Queyras) were obtained from the Safran-CROCUS (S2M) reanalysis(Météo France, 1981[START_REF] Bonet | Modelling the production and species richness of wild mushrooms in pine forests of the Central Pyrenees in northeastern Spain Can[END_REF].

Figure 2 .

 2 Figure 2. Annual SFGDD series with positive and negative anomalies with respect to the 1960-1988 average indicated by red and blue colors, respectively as well as mean (RWIs, black line) (1960-2017) for each of the eight sites within the French Alps. SFGDD values were transformed into Z-scores with respect to their mean over the 1960-1988 reference period. Dashed and solid red lines represent linear trends (non-significant and significant at p ⩽ 0.05) in SFGDD series for the time periods 1960-1988 and 1989-2017. Correlation coefficients between SFGDD and RWI are given for each period.

Figure 3 .

 3 Figure 3. Correlation matrices between detrended R. ferrugineum site chronologies for 1960-1988 (left panel) and 1989-2017 (right panel). Color intensity and circle size are proportional to the correlation coefficients (see color scale). * (p < 0.05) indicates the level of significance.

Figure 4 .

 4 Figure 4. Linear mixed models of R. ferrugineum RWs (log-transformed and scaled RW values) as a function of SFGDDs and their interaction with continentality (A) and slope aspect (B) for the period 1960-2017. Ranges for northness and Gams continentality index are indicated to the right side of each panel.

  2018)-indeed are the single-most important drivers of plant growth in cold environments by driving cell division and phenological transitions(Körner 2003, Kudo and Suzuki 2003, Wipf 2010, Hoch 2015). Our findings are also in line with dendroecological studies realized in the Alps (Pellizzari et al 2014, Francon et al 2017, Carrer et al 2019). These studies demonstrated the positive effect of summer temperatures and the negative effect of snowpack duration on shrub radial growth at their upper elevational limit. In the Arctic, the vast majority of dendroecological studies found summer temperatures to be the main driver of shrub growth(Blok et al 2011, Jørgensen et al 2015, Myers-Smith et al 2015a, Andreu-Hayles et al 2020).

Figure 5 .

 5 Figure 5. Conceptual model synthesizing potential responses of R. ferrugineum radial growth along the moisture gradient to increasing SFGDD in a warmer climate. In wetter sites (top left panel), increasing SFGDD are expected to favor radial growth.By contrast, at drier sites (bottom left panel), a non-linear response is anticipated and will likely limit radial growth above a certain SFGDD threshold.

  • C and 1700 • C for south-facing and north-facing slopes, respectively. However, one should keep in mind that the linear mixed model only explains 45% of the interannual RW variability with 24% attributed to the fixed effects. In that respect, taking both microsite (microtopography, soil properties) and (inter-)individual characteristics (competition, morphology, genetic/phenotypic variability) into account will facilitate the definition of potential drivers of shrubs growth further (Wipf et al 2006, Myers-Smith et al 2015a, Boulanger-Lapointe et al 2016, Bjorkman et al 2018, Berner et al 2020) and thereby also enhance the reliability of our approach.

  By contrast, warming climate negatively impacted alpine juniper shrub growth and shrub recruitment in the central Himalayas and Tibetan Plateau due to increasing drought stress (Lu et al 2019, Pandey et al 2020).

	attrib-
	uted to global warming and its effects on increasing
	soil moisture limitation or more frequent late frost
	events. Similar negative effects of increasingly higher
	summer temperatures have been highlighted in terms
	of productivity of alpine meadows in the European
	Alps (de Boeck et al 2016, Cremonese et al 2017,
	Corona-Lozada et al 2019), underlining that, if com-
	bined with water deficit, heat waves will lead to veget-
	ation browning. The loss of sensitivity to temperature
	that we highlight in this study echoes the 'divergence
	problem' that has been reported for trees growing
	in certain regions of circumpolar northern latitudes
	(Briffa et al 1998, Driscoll 2005, D'Arrigo et al 2008)

Moreover, a loss of sensitivity to temperature has been observed since the 1990s in Betula nana and Salix glauca shrubs from Western Greenland (Gamm et al 2018), Betula nana and Empetrum hermaphroditum shrubs in the central Norwegian Scandes (Weijers et al 2018a) and in Rumex alpinus herbs in the Low Tatras of Slovakia (Dolezal et al 2020) which was

  120 years for Alnus viridis and Salix glauca in West Greenland. Michalet R (ed) J. Veg. Sci. 26 155-65 Klein J A, Harte J and Zhao X-Q 2007 Experimental warming, not grazing, decreases rangeland quality of the Tibetan Plateau Bligny R, Gout E and Doche B 1996 Growth rates and nutrition status of an open and a closed population of Rhododendron ferrugineum L . in the northwestern Alps (France). Trees 11 91-98 Pornon A and Doche B 1996 Age structure and dynamics of Rhododendron ferrugineum L. populations in the northwestern French Alps J. Veg. Sci. 7 265-72 R Core Team 2016 R: a language and environment for statistical computing. R foundation for statistical computing. R foundation for statistical computing Vienna, Austria (available at: www.R-project.org/) Reid P C et al 2016 Global impacts of the 1980s regime shift Glob. Change Biol. 22 682-703 Riihimäki H, Heiskanen J and Luoto M 2017 The effect of topography on arctic-alpine aboveground biomass and NDVI patterns Int. J. Appl. Earth Obs. Geoinf. 56 44-53 Ropars P, Lévesque E and Boudreau S 2015 How do climate and topography influence the greening of the forest-tundra ecotone in northern Québec? A dendrochronological analysis of Betula glandulosa J. Ecol. 103 679-90 Scherrer D and Körner C 2009 Infra-red thermometry of alpine landscapes challenges climatic warming projections: thermometry of alpine landscapes Glob. Change Biol. 16 2602-13 Scherrer D and Körner C 2011 Topographically controlled thermal-habitat differentiation buffers alpine plant diversity against climate warming: topographical control of thermal-habitat differentiation buffers alpine plant diversity J. Biogeogr. 38 406-16 Schweingruber F H and Poschold P 2005 Growth rings in herbs and shrubs: life span, age determination and stem anatomy For. Snow Landsc. Res. 79 195-415 Sturm M, Holmgren J, McFadden J P, Liston G E, Chapin F S and Racine C H 2001 Snow-shrub interactions in Arctic Tundra: a hypothesis with climatic implications J. Clim. 14 336-44 Theurillat J P and Schlüssel A 2000 Phenology and distribution strategy of key plant species within the subalpine-alpine ecocline in the Valaisan Alps (Switzerland) Phytocoenologia 30 439-56 Treharne R, Bjerke J W, Tømmervik H, Stendardi L and Phoenix G K 2019 Arctic browning: impacts of extreme climatic events on heathland ecosystem CO2 fluxes Glob. Change Biol. 25 489-503 Vernay M, Lafaysse M, Hagenmuller P, Nheili R, Verfaillie D and Morin S 2019 The S2M meteorological and snow cover reanalysis in the French mountainous areas (1958-present) (available at: https://en.aeris-data.fr/metadata/?865730e8edeb-4c6b-ae58-80f95166509b) (Accessed 8 December 2020) Vicente-Serrano S M, Beguería S and López-Moreno J I 2010 A multiscalar drought index sensitive to global warming: the standardized precipitation evapotranspiration index J. Clim. 23 1696-718 Vionnet V, Brun E, Morin S, Boone A, Faroux S, Le Moigne P, Martin E and Willemet J-M 2012 The detailed snowpack scheme Crocus and its implementation in SURFEX v7.2 Geosci. Model Dev. 5 773-91 Weijers S, Beckers N and Löffler J 2018a Recent spring warming limits near-treeline deciduous and evergreen alpine dwarf shrub growth Ecosphere 9 e02328 Weijers S, Pape R, Löffler J and Myers-Smith I H 2018b Contrasting shrub species respond to early summer temperatures leading to correspondence of shrub growth patterns Environ. Res. Lett. 13 034005 Wigley T M L, Briffa K R and Jones P D 1984 On the average value of correlated time series, with applications in dendroclimatology and hydrometeorology J. Clim. Appl. Meteorol. 23 201-13 Wilmking M 2005 Increased temperature sensitivity and divergent growth trends in circumpolar boreal forests Geophys. Res. Lett. 32 Wilmking M, Juday G P, Barber V A and Zald H S J 2004 Recent climate warming forces contrasting growth responses of white spruce at treeline in Alaska through temperature thresholds Glob. Change Biol. 10 1724-36 Wipf S 2010 Phenology, growth, and fecundity of eight subarctic tundra species in response to snowmelt manipulations Plant Ecol. 207 53-66 Wipf S, Rixen C and Mulder C P H 2006 Advanced snowmelt causes shift towards positive neighbour interactions in a subarctic tundra community Glob. Change Biol. 12 1496-506 Zuur A F, Ieno E N and Elphick C S 2010 A protocol for data exploration to avoid common statistical problems: data exploration Methods Ecol. Evol. 1 3-14

	Ecol. Appl. 17 541-57
	Kolishchuk V 1990 Dendroclimatological study of prostrate
	woody plant Methods of Dendrochronology, Applications in
	the Environmental Sciences (Dordrecht: Springer) pp 51-5
	Körner C 2003 Alpine Plant Life: Functional Plant Ecology of High
	Mountain Ecosystems 2nd edn (Berlin: Springer)
	Kudo G and Suzuki S 2003 Warming effects on growth,
	production, and vegetation structure of alpine shrubs: a
	five-year experiment in northern Japan Oecologia 135 280-7
	Leonelli G, Pelfini M, Battipaglia G and Cherubini P 2009
	Site-aspect influence on climate sensitivity over time of a
	high-altitude Pinus cembra tree-ring network Clim. Change
	96 185-201
	Liston G E, Mcfadden J P, Sturm M and Pielke R A 2002 Modelled
	changes in arctic tundra snow, energy and moisture fluxes
	due to increased shrubs Glob. Change Biol. 8 17-32
	Lo T-T and Hsu H-H 2010 Change in the dominant decadal
	patterns and the late 1980s abrupt warming in the
	extratropical Northern Hemisphere Atmos. Sci. Lett.
	11 210-5
	Lu X, Liang E, Wang Y, Babst F, Leavitt S W and Julio Camarero J
	2019 Past the climate optimum: recruitment is declining at
	the world's highest juniper shrublines on the Tibetan
	Plateau Ecology 100 e02557
	Malfasi F and Cannone N 2020 Climate warming persistence
	triggered tree ingression after shrub encroachment in a high
	Alpine Tundra Ecosystems 23 1657-75
	Mekonnen Z A et al 2021 Arctic tundra shrubification: a review of
	mechanisms and impacts on ecosystem carbon balance
	Environ. Res. Lett. 16 053001
	Michalet R, Rolland C, Joud D, Gafta D and Callaway R M 2003
	Associations between canopy and understory species
	increase along a rainshadow gradient in the Alps: habitat
	heterogeneity or facilitation? Plant Ecol. 165 145-60
	Myers-Smith I H et al 2011 Shrub expansion in tundra
	ecosystems: dynamics, impacts and research priorities
	Environ. Res. Lett. 6 045509
	Myers-Smith I H et al 2015a Climate sensitivity of shrub growth
	across the tundra biome Nat. Clim. Change 5 887-91
	Myers-Smith I H et al 2015b Methods for measuring arctic and
	alpine shrub growth: a review Earth Sci. Rev. 140 1-13
	Myers-Smith I H et al 2020 Complexity revealed in the greening of
	the Arctic Nat. Clim. Change 10 106-17
	Nakagawa S and Schielzeth H 2013 A general and simple method
	for obtaining R 2 from generalized linear mixed-effects
	models Methods Ecol. Evol. 4 133-42
	Oberhuber W, Kofler W, Pfeifer K, Seeber A, Gruber A and
	Wieser G 2008 Long-term changes in tree-ring-climate
	relationships at Mt. Patscherkofel (Tyrol, Austria) since the
	mid-1980s Trees 22 31-40
	Ozenda P G 1985 La végétation de la chaîne alpine: dans l'espace
	montagnard européen. Masson: Paris u.a
	Pandey J, Sigdel S R, Lu X, Salerno F, Dawadi B, Liang E and
	Camarero J J 2020 Early growing-season precipitation drives
	radial growth of alpine juniper shrubs in the central
	Himalayas Geogr. Ann. A 102 317-30
	Pellizzari E, Pividori M and Carrer M 2014 Winter precipitation
	effect in a mid-latitude temperature-limited environment:
	the case of common juniper at high elevation in the Alps
	Environ. Res. Lett. 9 104021
	Phoenix G K and Bjerke J W 2016 Arctic browning: extreme
	events and trends reversing arctic greening Glob. Change
	Biol. 22 2960-2
	Pinheiro J C and Bates D M 2000 Mixed-effects Models in S and
	S-PLUS (Berlin: Springer)
	Pornon A,

© 2021 The Author(s). Published by IOP Publishing Ltd

Acknowledgments

This research was supported by the program ProXyC-liM (2016ProXyC-liM ( -2017) ) funded by the CNRS (Mission interdisciplinarité-RNMSH) and by the ANR/FWF ACOUFOLLOW project of the French Agence Nationale de la Recherche (ANR-20-CE91-0008) and the Austrian Science Found (FWF I 4918). We wish to thank all who contributed to the meteorological series data used in the analyses, especially Carlo Maria Carmagnola (Météo-France-CNRS, CNRM, CEN). The S2M data are provided by Météo-France-CNRS, CNRM Centre d'Etudes de la Neige, through AERIS. We acknowledge Olivier Voldoire (GEOLAB, CNRS) for technical support. Rhododendron shrubs were sampled with the help of Séverine Finet, Sébastien Guillet, and Robin Mainieri. CNRM/CEN and LECA are part of LabEX OSUG@2020.

Data availability statement

The data that support the findings of this study are available upon reasonable request from the authors.