Frédéric Armetta
email: frederic.armetta@univ-lyon1.fr

Anthony Baccuet
email: anthony.baccuet@gmail.com

Mathieu Lefort
email: mathieu.lefort@liris.cnrs.fr

Algorithmic learning, a next step for AI. An application to arithmetic operations

Keywords: Algorithmic Machine Learning, Deep Learning, Active learning

Deep learning achieved state of the art performances in multiple domains (image recognition, natural language processing, etc.) One of the next steps is to be able to learn algorithms, as a way to provide some new forms of generalization for AI systems. This is currently a hard and challenging problem as it involves algorithmic recurrence, memory management and combination of subtasks, which leads to trainability problems. We show that even a simple algorithm of multiplication manifests trainability problems for neural networks. In this article, we present an original training method applied to multi-digit multiplication learning, called Unrolling Algorithmic Training (UAT). To learn the global algorithm, we use additional supporting tasks consisting in the successive subtasks composing the global algorithm (in our case the 1-digit multiplications and the final addition). The global end-to-end and the subtasks learning are then balanced with an active learning mechanism. This multi-task learning allows to overcome the problem of trainability encountered when learning directly the global algorithmic task. Our interpretation is that the network is able to somehow combine the subtasks in order to learn the global task. Moreover, we show that the global algorithm can be bootstrapped, fine-tuned and is even resilient without retraining it from scratch when we vary the size of recurrence provided to the network.

Introduction

Since the first success of convolutional neural networks on image classification [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF], deep learning methods have pushed forward the state of the art in multiple domains [START_REF] Goodfellow | Deep learning[END_REF]. This flows in tasks of increasing complexity such as language translation [START_REF] Vaswani | Attention is all you need[END_REF] or game playing [START_REF] Tian | Elf opengo: An analysis and open reimplementation of alphazero[END_REF]. One of the next steps is algorithmic learning, including mathematical expression calculation, as it will open the way to learn any task that can be expressed in a Turing machine and so to provide greater autonomy to AI systems.

Algorithmic resolution through neural networks is still an emerging area of research. It is a challenging problem as it requires to memorize values for a long period of time, to learn inferences, to combine procedures, to extrapolate to unknown domains, etc. More fundamentally this needs to fill the gap between symbolic understanding and statistical learning. Neural Turing Machine [START_REF] Graves | Neural turing machines[END_REF] was proposed to learn end-to-end algorithmic procedures, such as list sorting. It aims to reduce the trainability problem faced by RNN, which are Turing complete [START_REF] Siegelmann | On the computational power of neural nets[END_REF], by adding specific mechanisms such as a memory and differentiable ways to access it. However, it also suffers from the trainability problem itself, i.e. the learning procedure is very sensitive to the hyperparameters and/or to the initial values and may lack generality, so that the performances are hard to reproduce and inconstant [START_REF] Collier | Implementing neural turing machines[END_REF]. These problems are due to the depth of the proposed architecture [START_REF] Kaiser | Neural gpus learn algorithms[END_REF] but also to the intrinsic complexity of the operations to learn, especially the long-term dependencies between data or variables to manipulate.

To overcome this trainability problem for an algorithm learning, an alternative way is rather to decompose this algorithm into its successive steps so that the neural network will compute and learn iteratively by feeding back its previous output thanks to a handmade recurrence external to the model. This procedure was applied to a transformer-inspired architecture for learning some algorithmic procedures [START_REF] Yan | Neural execution engines: Learning to execute subroutines[END_REF] and to a MLP for arithmetic operations [START_REF] Nollet | Learning Arithmetic Operations With A Multistep Deep Learning[END_REF][START_REF] Yan | Neural execution engines: Learning to execute subroutines[END_REF]. In this paper, we aim to address the trainability problem end-to-end.

We choose to focus and illustrate our proposal on arithmetic operations, especially the multidigit multiplication algorithm, which requires multiple one-digit multiplications to add properly in order to get the final result. This computation involves many operations and variables in the calculation stream. The complexity for such an operation is also related to the propagation of the carry e.g. [START_REF] Cho | Problem difficulty in arithmetic cognition: Humans and connectionist models[END_REF]. Thus, learning directly the end-to-end multiplication leads to poor performances [START_REF] Hoshen | Visual learning of arithmetic operations[END_REF]. To tackle this trainability problem, we propose to mix during learning both intermediate sub-tasks and the global end-to-end multiplication, weighted by an active learning strategy. This is somehow similar to a multi-task learning, using intermediate computation as relevant complementary tasks to help the network to learn the global multiplication. However, with our approach, the network will perform all these tasks without any layer dedicated to each one. Thus, the network will have to accommodate intermediate operations in order to learn the global multiplication. Another difference is that this accommodation has to follow the course of the algorithm, by sequentially connecting tasks and propagating intermediate values through recurrences. We will show that this training procedure can also be used for pre-training the network, as a kind of bootstrap, for an efficient fine-tuning relying on the algorithmic target only (the multiplication).

We present the related work in section 2. We show in this section that language models such as chatGPT are also highly concerned by algorithmic learning, and suffer from trainability problems in the same way. This can be observed when the query involves lengthy algorithmic reasoning, as for arithmetic operations. We detail the use case of multi-digit multiplications in section 3. The chosen task is the decimal multi-digit multiplication, that is one of the hardest of the four arithmetic operations as it implies a lot of steps. Then, the model on which our learning procedure is applied is detailed in section 4 and next the protocol and results are presented in section 5. We conclude and discuss the perspectives of our work in section 6.

State of the art

Over the years, multiple deep learning models have been proposed to solve different kinds of problems [START_REF] Goodfellow | Deep learning[END_REF]. Convolutional neural networks [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF] are able to classify images with performances sometimes overcoming the ones of humans in specific scenarios. Recurrent neural networks with dedicated cells such as Long Short Term Memory (LSTM) [START_REF] Hochreiter | Long short-term memory[END_REF] or Gated Recurrent Units (GRUs) [START_REF] Cho | Learning phrase representations using rnn encoder-decoder for statistical machine translation[END_REF] can deal with temporal data.

Deep learning models are universal function approximators [START_REF] Cybenko | Approximation by superpositions of a sigmoidal function[END_REF] but they are often hard to train. For instance, in theory even a shallow architecture is sufficient, but in practice, the ability of deep networks to learn relevant representations from data is much better [START_REF] Ba | Do deep nets really need to be deep?[END_REF]. Moreover, RNNs are Turing complete [START_REF] Siegelmann | On the computational power of neural nets[END_REF], but also face a trainability problem [START_REF] Graves | Neural turing machines[END_REF]. This limitation tends to be more pronounced for complex tasks, especially when long-term inferences are needed. As an example, the error rate of an arithmetic operation is related to the number of carries for a MLP [START_REF] Cho | Problem difficulty in arithmetic cognition: Humans and connectionist models[END_REF].

Some global strategies were proposed to make network training easier and more effective. Among them we can mention active learning strategies that consist in finding the right next example to improve the training progress [START_REF] Settles | Active Learning Literature Survey[END_REF][START_REF] Burda | Large-scale study of curiosity-driven learning[END_REF]. A common way to do it is curriculum learning where tasks are ordered by increasing complexity [START_REF] Bengio | Curriculum learning[END_REF]. Another strategy is multi-task learning that consists in combining different tasks with similar objectives so that the network can better generalize and learn the underlying structure of data [START_REF] Zhang | A survey on multi-task learning[END_REF]. Even if our proposal appears similar to multi-task learning, in our model the tasks are taken on the same network and are of algorithmic nature.

Neural Turing Machine (NTM) [START_REF] Graves | Neural turing machines[END_REF] was specifically designed to learn algorithmic tasks. To limit the problem of long time dependency encountered by LSTM, it mimics some of the principles of a Turing Machine by introducing a memory and reading/writing mechanisms. Differentiable Neural Computer [START_REF] Graves | Hybrid computing using a neural network with dynamic external memory[END_REF] improves these accesses to learn more complex tasks as some NLP problems of inference and reasoning. However, these models are difficult to train as they seem to be very sensitive to initial weights [START_REF] Collier | Implementing neural turing machines[END_REF][START_REF] Saxton | Analysing mathematical reasoning abilities of neural models[END_REF]. The neural GPU architecture [START_REF] Kaiser | Neural gpus learn algorithms[END_REF] shares similar ideas with NTM but uses convolutional GRU in order to obtain parallel computing. It is able to learn some algorithmic tasks such as binary addition and multiplication, sequence reversing, etc. Its main achievement is its ability to generalize the learning to inputs with longer size in testing, e.g. from 20-digit binary multiplication in learning up to 2000-digit ones in testing, without any error. However, this model completely fails to learn decimal multiplications [START_REF] Freivalds | Improving the neural gpu architecture for algorithm learning[END_REF]. This limitation may be related to carry propagation that is hard to train and appears more frequently with decimal coding of numbers and can be partially overcome by using curriculum learning (from binary to quaternary then to decimal) [START_REF] Price | Extensions and limitations of the neural gpu[END_REF].

In [START_REF] Saxton | Analysing mathematical reasoning abilities of neural models[END_REF], multiple algorithms (LSTM with or without attention, transformers) were tested on various mathematical inference tasks. The main objective was to propose a dataset and to compare the models, especially on the attentional aspect, so once again detailed performances are missing. One of the conclusions is that long-term dependencies are the harder to learn which can be compromising for the addition of multiple operands in specific cases. A transformer with operands given as variables in the text, achieves good performance except for subtraction and multiplication [START_REF] Wangperawong | Attending to mathematical language with transformers[END_REF]. [START_REF] Lample | Deep learning for symbolic mathematics[END_REF] solves mathematical equations, including differential ones, with Seq2Seq models. The originality lies in the equation being written as the prefix notation of its tree representation. Another approach proposes to solve arithmetic expressions by composing single-digit sub-tasks. This hierarchical combination is learned by reinforcement learning with curriculum [START_REF] Chen | Neural arithmetic expression calculator[END_REF]. The system manages to generalize to some extent to longer sequences, but the performance sometimes drops, especially for multiplication. Neural Arithmetic Logic Units [START_REF] Trask | Neural arithmetic logic units[END_REF] are cells that are able to perform arithmetic operations by introducing specific computation modules such as log, exp, etc. The aim here is not to learn arithmetic per se but to provide dedicated cells able to extrapolate learning with computations that extends to unknown domains. A direct extension dealing also with negative inputs was proposed in [START_REF] Schlör | inalu: Improved neural arithmetic logic unit[END_REF]. Other derived models can compute arithmetic operations on real numbers [START_REF] Madsen | Neural arithmetic units[END_REF]. As we saw, there is a huge variety of tasks that was explored in the literature, so that comparing different approaches is difficult.

Algorithmic reasoning is also part of the required material for Natural Language Processing. Thus, large Language Models which have been extensively developed in recent years, need to learn reasoning to answer appropriately. They can for instance learn to perform numerical reasoning when learning few examples in a few-shot setting (using a query such as "Q: What is 24 times 18?", taking GPT-J-6B as a base). It was shown that the performance is highly correlated with the frequency of the terms in the pretraining corpus [START_REF] Razeghi | Impact of pretraining term frequencies on few-shot numerical reasoning[END_REF]. When co-occurrence of terms is low, accuracy is also low. These observations underline the trainability problem and lack of generalisation to algorithms learning. As a consequence, performance decreases as the size of operands increases [START_REF] Brown | Language models are few-shot learners[END_REF]. Moreover, while addition seems to be less impacted by the size of operands, performance collapses for the multiplication which requires more reasoning or algorithmic processing. GPT-3 is among the most popular and large language model, using about 500 billion tokens for learning [START_REF] Ouyang | Training language models to follow instructions with human feedback[END_REF]. ChatGPT which is fine-tuned from GPT-3 also suffers from the limitations of GPT-based models. These observations suggest that model reasoning skills for such models are still limited. They highly rely on the size of the corpus available and statistics. The corpus will never contain all the combinations of terms or parameters of algorithms we have to learn. A better approach is then to enhance the ability to learn algorithms, which would allow to make better predictions for unknown configurations.

In this article, we choose to focus on the multi-digit multiplication of two decimal numbers. This is a simple enough task to not mix different problems but in the meantime it is challenging as multiple models are unable to learn it properly. This difficulty arises from long-term dependency due to the carry propagation, but also and more generally from the inherent complexity of algorithms which involves many operations and variables in the calculation flow.

Moreover, some articles precisely measured performance on this task. [START_REF] Hoshen | Visual learning of arithmetic operations[END_REF] proposes a MLP to learn addition and multiplication either from visual inputs or numerical encoding. In both cases, the accuracy achieved for multiplication is poor (see section 5).

Problem statement

In this article we will consider the decimal multi-digit multiplication. Formally, let 𝑛 be the (maximum) number of digits of any of the two operands. The multiplication of these two operands leads to 𝑛 + 1 sub-tasks: 𝑛 single-digit multiplications and 1 final addition of the partial multiplications (see figure 1 for an example). In our data representation, each operation will correspond to two lines of computation: one for the carries and one for the result. The maximum size of any intermediate operation is 𝑁 = 2𝑛 (this maximal length will be obtained for the final addition with a carry generated at the most significant digit position). All the operands will be padded with zero digits to match a 𝑁 size, yet there will be exactly 𝑛 + 1 intermediate operations even if the first operand has less than 𝑛 digits.

(1)

×0048

(2)

0012

(3) 2) are the two operands. Lines ( 3) and (4) (respectively ( 5) and ( 6)) represent the carries and result of the 1-digit multiplication of 8 (respectively 4) by [START_REF] Freivalds | Improving the neural gpu architecture for algorithm learning[END_REF]. Lines ( 7) and ( 8) correspond to the carries and result of the addition of lines ( 4) and ( 6), which is also the final result of the global multiplication of 48 by 23, i.e. 1104.

0184 (4) 0010 (5) +0920 (6

Model

Our model is very similar to networks used for natural language processing. We choose a recurrent network in order to provide a potential for algorithm unrolling. We use an agnostic encoding for digits (no binary encoding) in order to prevent any facilities associated with the selected problem and focus on the algorithmic concerns.

Data representation

Each digit 𝑑 is represented by a 10-dimension vector with a one-hot encoding, i.e. (𝛿 𝑑𝑖 ) 𝑖∈{0,...,9} . This vector can also have two other values. First, a null vector represents an empty value in the input to form empty lines (see table 1). Note that it will also be used as the starting character for the decoder presented below. Second, a one vector corresponds to the end of the line either when reading or writing the data.

Model architecture

We use a Seq2Seq model as proposed in [START_REF] Sutskever | Sequence to sequence learning with neural networks[END_REF] (see figure 2). It is composed of an encoder that will sequentially read the data and recurrently embed it in some hidden states. From it, a decoder will iteratively produce a sequence from previously outputted characters, beginning with some predetermined code, until an ending character is written. Both encoder and decoder are implemented with an LSTM neural network model. During training, we used teacher forcing, i.e. the next character is not produced from the previous output but from the ground truth. The current output is then compared to the expected one and the error is backpropagated through the decoder and the encoder.

A headband allows to read and write successive lines two at a time, digit by digit. The digits, encoded in one-hot vectors, are then read from right to left, with the <eol> vector at the end of each line. Similarly, encoded one-hot vectors are written at the outputs. Both the inputs of the encoder and output of the decoder are 2 × 10 in size. Each task can then be encoded individually in order to infer the ground truth value (see below). In practice each digit is encoded as a one-hot vector (see section 4.1).

Tasks definition

The multiplication algorithm involves sub-tasks that can be submitted to the model (see table 1). The encoder records the formulated task, the decoder outputs the target value. The multiplication algorithm can also be inferred end-to-end in the same way. As formalized in section 3, the multiplication operation of two n-digit numbers can be decomposed into n single-digit multiplications and 1 final addition. These n+1 operations and the global end-to-end multiplication (ie computing the final result and carries from the two operands) compose the n+2 tasks that the network has to learn. For each task, the encoder reads the corresponding inputs and the decoder outputs the target value of the corresponding operation as illustrated in table 1. So that the network can differentiate the task of the first single digit multiplication and the end-to-end multiplication, which both require only the two operands, in the latter case we add 2n empty lines (corresponding to unfilled intermediate results lines ). Sub-tasks and end-to-end inferences associated with a 4-digit multiplication (see numbered lines in figure 1, lines are read and written two at a time ("//"))

Active learning mechanism

We have so far described several tasks. One can choose to select a bundle of tasks to train concurrently on the model. In this article, different settings are compared (sub-tasks only, all the tasks, end-to-end only) for training.

In order to balance the learning effort between the chosen bundle tasks, in all cases we use the same active learning mechanism as proposed in [START_REF] Nollet | Learning Arithmetic Operations With A Multistep Deep Learning[END_REF]

+ (1 -𝜆) 1 𝑐𝑎𝑟𝑑(taskList)
where 𝜆 is an hyperparameter and taskList is the list of all tasks involved. The general idea is that the more difficult is a task to learn (first term), the more it is present in the next epoch with a lower bound depending on the 𝜆 value.

Experiments

Protocol

For training, 100000 unique couples of operands are generated. At each epoch, in order to update the dataset, we generate and attribute an operation to each of the couple so that the global distribution between operations matches the one decided by the active learning mechanism presented in section 4.4.

For the validation and test datasets, we use respectively 1000 and 10000 additional unique couples of operands. The size of the latent space of the encoder is set to 500 and the active learning parameter 𝜆 to 0.5. The model is trained using the ADAM optimizer with a learning rate of 10 -4 and a batch size of 10. All the results presented in the next sections are averaged over 4 runs that learned during 500 epochs.

Results

test = end-to-end

3 × 3 (6 output digits) restricted 4 × 4 (7 output digits) 4 × 4 (8 output digits)
Hoshen et al. [START_REF] Hoshen | Visual learning of arithmetic operations[END_REF] n.c 37.6 % n.c. UAT (train = end-to-end only) 4.05% (± 1.72%) 35.87% (± 28.10%) 92.42% (± 3.64%) UAT (train = sub-tasks and end-to-end) 3.33% (± 1.32%) 4.51% (± 1.21%) 23.34% (± 14.69%)

Table 2

Error rate on the end-to-end multiplication task

In this section, we want to quantify the effect of our proposition to mix the learning of the sub-tasks with the global end-to-end multiplication. For the purpose of comparison, we reproduce 4 × 4 multiplications, as presented in [START_REF] Hoshen | Visual learning of arithmetic operations[END_REF] and [START_REF] Nollet | Learning Arithmetic Operations With A Multistep Deep Learning[END_REF], with 4-digit inputs chosen so that the final outputs are restricted to 7-digit numbers.

Sub-tasks and end-to-end

The main purpose of this paper is to tackle the trainability problem encountered for the end-toend multiplication task we choose as a first step towards general algorithmic learning. For the so considered bundle of tasks, sub-tasks and the end-to-end task are trained simultaneously on the same network following the active learning mechanism. We report in table 2 the performances for the end-to-end task, and compare with [START_REF] Hoshen | Visual learning of arithmetic operations[END_REF] which uses an MLP while we use a Seq2Seq. We also vary the scale and complexity of problems in order to show the improvements of our proposal.

Results clearly demonstrate the contribution of sub-tasks for the end-to-end task, lowering the error rate (from 35.87% to 4.51 % for the 4 × 4 restricted multiplications). This suggests that the model is able to self-organize and take advantage of sub-tasks, which facilitate the training, validating that our proposed training procedure is the key element of our model. While results are similar without the support of sub-tasks, i.e. both architectures face trainability problems, our training approach clearly outperforms [START_REF] Hoshen | Visual learning of arithmetic operations[END_REF]. The standard deviation reduction shows that stability of learning is also improved (with the exception of UAT end-to-end only 8 output digits that can be excluded from the comparison because of homogeneous but poor efficiency).

We report below other significant improvements for the most difficult problem presented (the 4 × 4 (8 output digits) problem) thanks to fine tuning and additional recurrences. 3.84% (± 1.61%) 7

3.83% (± 1.29%)

Table 4

Influence of the recurrence provided to the model when fine tuning (pre-training = sub-tasks + end-toend, 4 free double lines)

End-to-end only (fine tuning)

To improve further the performance of our model, we propose to fine tune the end-to-end task only, during 500 additional epochs, for the hardest 4×4 problem. As we can see on table 3, when keeping the same recurrence to the network, the fine tuning leads to a significant improvement of the performance as the error drops from 23.34% to 6.68% for the 8 output digits problem. This result also shows that once the end-to-end task is bootstrapped (by simultaneous learning with subroutines), its performance can be further improved. On the contrary, table 3 shows that not including the end-to-end task in the pre-training bundle leads to trainability problems. This confirms that the model has to learn simultaneously the intermediate operations and the end-to-end one to overcome the trainability problem.

Adaptability

In our model, the input of the end-to-end multiplication is composed by some empty lines to match with the input size of the final addition subroutine (see table 1). These empty lines provide recurrence steps for the network (encoder). To test the influence of this factor, we ran multiple fine tuning of our model, setting each with a different number of empty lines (initial training done with 4 double lines). We can observe on table 4 that the error rates tend to decrease monotonically with the number of free recurrences. This may sound logical as increasing the number of recurrences also increases the computational power of the model.

Active learning dynamic

In section 4.4, we introduced the active learning mechanism that we used to balance the learning effort between the different type of operations involved. In figure 3 is represented the evolution of their distribution in the training dataset. We can observe that it is relatively constant after some time and that the hardest task is the end-to-end one. However, one can notice that the second hardest task is the addition procedure and not any of the intermediate multiplication, whose mean error rate is even close to 0. This is surprising as in the literature the addition appears to be a simpler task to learn. This may be due to the fact that in this case the addition involves many operands.

Sub-tasks only

We will train the network only on the sub-tasks (ie excluding the end-to-end operation). To obtain the global result, we will provide the network with the successive intermediate tasks in the right order (with intermediate lines filled iteratively with previous outputs of the network) and consider only the output of the final addition as the global result. This method was also used in [START_REF] Yan | Neural execution engines: Learning to execute subroutines[END_REF]. This task is easier to process (as the global operation is decomposed in its successive steps for its execution). Especially, we aim to estimate the effect of the recurrence available in our network knowing that [START_REF] Nollet | Learning Arithmetic Operations With A Multistep Deep Learning[END_REF] uses a network without recurrence.

Table 5 shows that our model reduces the error rate, in this context of execution. This shows that the recurrent model and data representation we use succeeds in capturing all the Error rate comparison between [START_REF] Nollet | Learning Arithmetic Operations With A Multistep Deep Learning[END_REF] and UAT, when combining sub-tasks only thanks to an handmade recurrence

Conclusion and perspectives

Despite its multiple successes, deep learning architectures still struggle to learn algorithms, such as arithmetic operations. Algorithmic learning is highly expected as it would allow new classes of problems to be addressed, such as algorithmic generalisation. This is nevertheless a hard task, due to the complexity of the learning, especially the long-term dependencies, but also a trainability problem faced by multiple models. This limitation can as well be observed for language models whose performances come from the corpus available, and do not apply well to algorithmic inferences when variability in parameter values is large and usual generalisation do not apply as observed for arithmetic calculations, which underlines the need for new inference mechanisms.

In this article, we propose an original way to learn the end-to-end multi-digit multiplication of two (decimal) operands by guiding the model via the introduction of sub-tasks (or subroutines) concurrently with the targeted end-to-end task while applying an active learning strategy. While in the literature the addition appears to be a simple task to learn, multiplication is especially challenging. We show through analysis of our experiments that the multiplication can nevertheless be learned end-to-end, we measure improvements directly caused by UAT, the learning procedure we introduce. Learning the target end-to-end task after the sub-tasks can improve the results but the best way, by a large margin, is to mix all the tasks together during training in order to bootstrap the target end-to-end training task. Once bootstrapped, the targeted task can efficiently be fine-tuned alone. This process can be described as a new algorithmic transfer, i.e. the tasks complement and support each other to achieve the overall arithmetic task.

By fine tuning the model on the final end-to-end task we show an interesting additional property of our learning method. Once the end-to-end multiplication learned, the network can adapt to a number of recurrence different from the one it was first been trained for. This seems to indicate that the model is able not only to combine the sub-tasks to resolve the global endto-end task but also to autonomously extract and to adapt some kind of high level algorithmic knowledge from it. Restricting the provided recurrences for the computing and maintaining accuracy looks like a constrained parallelization of the algorithmic task.

The main motivation for this work is not only the learning of the multiplication, but to provide a new method for alleviating the hard trainability problem that is observed for algorithms. This work raises multiple research questions. We want to investigate more precisely how the transfer from intermediate steps to the global task is achieved by the network. To overcome the trainability problem encountered by the classical training procedure, we provide all the intermediate steps to the network during learning. A specific case that we want to study is to provide only some of the supporting tasks and see if the network can complement the unknown tasks by itself. For that we will investigate how the transfer from intermediate steps to the global task is achieved in the network. This understanding may give us a way to control the flow of interaction between supporting steps and the general dynamic of the algorithmic transfer, and thus get a more detailed understanding of how an AI could exploit and adapt its algorithmic knowledge in order to expand its capabilities.

Figure 1 :

 1 Figure 1: Example of the representation of a multiplication of two 2-digit operands. Note that the signs (+ and ×) and lines are only shown for clarity. Lines (1) and (2) are the two operands. Lines (3) and (4) (respectively (5) and (6)) represent the carries and result of the 1-digit multiplication of 8 (respectively 4) by[START_REF] Freivalds | Improving the neural gpu architecture for algorithm learning[END_REF]. Lines (7) and (8) correspond to the carries and result of the addition of lines (4) and (6), which is also the final result of the global multiplication of 48 by 23, i.e. 1104.

Figure 2 :

 2 Figure 2: Seq2Seq architecture. The encoder receives successively the two digits (here boxed) of two lines (here the operands 23 and 48) from right to left. From the encoder embedding, the decoder recurrently produces the digits, from right to left, of the two next lines (here the carry and result of the first intermediate multiplication 8 × 23). In practice each digit is encoded as a one-hot vector (see section 4.1).

4x4 ( 8

 8 output digits) error rates initial UAT train = sub-tasks + end-to-end 23.34% (± 14.69%) UAT fine tuning (pre-training = sub-tasks + end-to-end) 6.68% (± 4.31%) UAT fine tuning (pre-training = sub-tasks) 73.31% (± 11.10%)

Figure 3 :

 3 Figure 3: (Top) Evolution of the error rate for each kind of operation (estimated from the training dataset) used by the active learning. (Bottom) Resulting evolution of the distribution of the various tasks in the training dataset due to active learning

Table 1

 1

	Task	Encoder inputs Decoder output
	st1	(1//2)	(3//4)
	st2	(1//2) (3//4)	(5//6)
	st3	(1//2) (3//4) (5//6)	(7//8)
	end-to-end (1//2) (empty line * 4)	(7//8)

  . It consists in measuring the error rate, denoted err task , of each involved task in the training dataset at the end of each training epoch. For the next epoch, the training dataset is constructed by randomly picking up examples from a global fixed set of examples so that each task is present with a fraction

	𝐹 task = 𝜆	∑	err t ta∈taskList err ta

Table 3

 3

	Influence of fine tuning
	free double lines	error rates
	0	88.67% (± 5.03%)
	1	50.19% (± 18.16%)
	2	14.39% (± 2.65%)
	3	4.44% (± 1.65%)
	4	6.68% (± 4.31%)
	5	5.53% (± 2.86%)
	6

Acknowledgment

This work was performed using HPC resources from GENCI-IDRIS and a GPU donated by @NVIDIA Corporation. We gratefully acknowledge this support.


Monique N Simon

Priscila S Rothier
email: priscilasrd@gmail.com

Colin M Donihue
email: colindonihue@gmail.com

Anthony Herrel
email: anthony.herrel@mnhn.fr

Jason J Kolbe
email: jjkolbe@uri.edu

Can extreme climatic events induce shifts in adaptive potential? A conceptual framework and empirical test with Anolis lizards

Multivariate adaptation to climatic shifts may be limited by trait integration that causes genetic variation to be low in the direction of selection. However, strong episodes of selection induced by extreme climatic pressures may facilitate future population-wide responses if selection reduces trait integration and increases adaptive potential (i.e., evolvability). We explain this counter-intuitive framework for extreme climatic events in which directional selection leads to increased evolvability and exemplify its use in a case study. We tested this hypothesis in two populations of the lizard Anolis scriptus that experienced hurricaneinduced selection on limb traits. We surveyed populations immediately before and after the hurricane as well as the offspring of post-hurricane survivors, allowing us to estimate both selection and response to selection on key functional traits: forelimb length, hindlimb length, and toepad area. Direct selection was parallel in both islands and strong in several limb traits.

Even though overall limb integration did not change after the hurricane, both populations showed a non-significant tendency toward increased evolvability after the hurricane despite the direction of selection not being aligned with the axis of most variance (i.e., body size).

The population with comparably lower between-limb integration showed a less constrained response to selection. Hurricane-induced selection, not aligned with the pattern of high trait correlations, likely conflicts with selection occurring during normal ecological conditions that favors functional coordination between limb traits, and would likely need to be very strong and more persistent to elicit a greater change in trait integration and evolvability. Future tests of this hypothesis should use G-matrices in a variety of wild organisms experiencing selection due to extreme climatic events.

Introduction

Despite the vast literature focusing on adaptation to climate, predictions for how species will respond to climate change generally ignore the degree of adaptive potential of populations in the face of these changes, especially when multiple traits are under selection (Nadeau and Urban 2019; but see [START_REF] Etterson | Constraint to Adaptive Evolution in Response to Global Warming[END_REF][START_REF] Etterson | Evolutionary potential of Chamaecrista fasciculata in relation to climate change. II. Genetic architecture of three populations reciprocally planted along an environmental gradient in the Great Palins[END_REF][START_REF] Peschel | Estimating the capacity of Chamaecrista fasciculata for adaptation to change in precipitation[END_REF]. Adaptive responses to climatic pressures can alter ecological responses of populations, lowering their vulnerability to future changes, or even rescue populations from extinction [START_REF] Quintero | What determines the climatic niche width of species? The role of spatial and temporal climatic variation in three vertebrate clades[END_REF]Wiens 2013, Diniz-Filho et al. 2019). Of course, these adaptive responses depend on standing genetic variation in the traits under selection, the frequency and severity of climatic perturbations that induce selection and even on interactions with other species [START_REF] Grant | Evolution caused by extreme events[END_REF]. For instance, we may expect adaptation to be especially challenging in response to extreme climatic events that induce very strong episodes of selection, leading to high mortality and small population sizes that may never rebound. Given that extreme climatic events are likely to become more frequent and intense with future climate change, such as accelerating warming rate [START_REF] Fischer | Increasing probability of record-shattering climate extremes[END_REF], the potential of populations facing extreme events to adapt and persist may be limited. However, the scenario may not be so daunting if populations can increase their adaptive potential in response to predicted extreme climatic events.

In this paper, we lay out a conceptual framework for the intriguing possibility that directional selection may increase the adaptive potential of populations by reorganizing genetic variance such that it aligns more with the direction of selection [START_REF] Pavlicev | Evolution of adaptive phenotypic variation patterns by direct selection for evolvability[END_REF], such as when selection is caused by extreme climatic pressures. Then we show how to apply this framework with a case study of hurricane-induced selection on functional traits in Anolis lizards. Although our data set has some important limitations which we discuss, such as only phenotypic (and no genetic) data and limited sample sizes, our goal is to stimulate future research on evolvability and extreme climatic events, a theme that is currently underappreciated in the literature.

Genetic integration, constraint and evolvability

Given that selection induced by climate is likely on several traits at once, an important determinant of population-level responses is phenotypic and genetic integration among traits [START_REF] Etterson | Constraint to Adaptive Evolution in Response to Global Warming[END_REF]. Theoretical models posit that if selection induced by changes in the climate acts on several traits (i.e., multivariate selection; [START_REF] Lande | Quantitative Genetic Analysis of Multivariate Evolution, Applied to Brain: Body Size Allometry[END_REF][START_REF] Lande | The measurement of selection on correlated characters[END_REF] that are tightly integrated-showing high genetic correlations-the response to selection may deviate from the direction of selection [START_REF] Hellmann | Constraints and reinforcement on adaptation under climate change: Selection of genetically correlated traits[END_REF]Pineda-Krch 2007, Chevin 2012). Genetic trait integration is measured by the additive genetic variancecovariance matrix, the G-matrix (G), which has trait variances along its diagonal and trait covariances in the off-diagonal and characterizes the pattern of trait co-inheritance [START_REF] Lande | Quantitative Genetic Analysis of Multivariate Evolution, Applied to Brain: Body Size Allometry[END_REF][START_REF] Arnold | Understanding the evolution and stability of the G-matrix[END_REF]. The response to selection is then determined by an interaction between directional selection and the G-matrix: ∆𝑧 = 𝐺𝛽; in which ∆𝑧 is the evolutionary response vector and 𝛽 is the selection gradient, which measures only the direct selection on traits by accounting for phenotypic trait correlations [START_REF] Lande | Quantitative Genetic Analysis of Multivariate Evolution, Applied to Brain: Body Size Allometry[END_REF]Arnold 1983, Cheverud 1984).

High genetic integration among traits results in a very stretched 𝐺 that accumulates additive variance in only one or a few trait combinations (a 'cigar-shaped' G-matrix, in contrast with a 'ball-shaped' G-matrix; see Figure 1), while most trait combinations will have little additive genetic variance (Arnold et al. 2001, Walsh and[START_REF] Walsh | Abundant Genetic Variation + Strong Selection = Multivariate Genetic Constraints: A Geometric View of Adaptation[END_REF]. Therefore, the shape of 𝐺 (i.e., how genetic variance is distributed along trait combinations) is crucial for determining how much additive variation will be available in the direction of 𝛽, that is, how high is the evolvability or the ability of a population to adapt (Houle 1992, Hansen and[START_REF] Hansen | Measuring and comparing evolvability and constraint in multivariate characters[END_REF]. In highly integrated systems, if the direction of 𝛽 is not aligned with the combinations of traits accumulating genetic variance, G will have low evolvability in the direction of 𝛽 and ∆𝑧 will deviate from the direction of selection [START_REF] Arnold | The adaptive landscape as a conceptual bridge between micro-and macroevolution[END_REF]; scenario B in Figure 1). Thus, whether populations will respond adaptively to climate-induced selection depends on evolvability: the higher the evolvability, the more likely the response to selection will be effectively adaptive. However, populations may be able to escape low evolvability if selection induced by climate events can change the shape of 𝐺 to become more aligned with the direction of selection (i.e., more additive genetic variance in the direction of selection; e.g., [START_REF] Assis | Directional selection effects on patterns of phenotypic (co)variation in wild populations[END_REF].

Rapid shifts in evolvability and strong selection induced by climatic extremes [START_REF] Pavlicev | Evolution of adaptive phenotypic variation patterns by direct selection for evolvability[END_REF] presented a theoretical model in which directional selection changes the pattern of trait integration (i.e., changing the shape of G), by acting on genes that regulate the degree of association between traits (the so-called 'rQTLs', relationship Quantitative Trait Loci; [START_REF] Cheverud | Pleiotropic effects on mandibular morphology II: Differential epistasis and genetic variation in morphological integration[END_REF][START_REF] Pavlicev | Genetic variation in pleiotropy: Differential epistasis as a source of variation in the allometric relationship between long bone lengths and body weight[END_REF]. The fundamental idea is that rQTLs change the effects of loci that control more than one trait (pleiotropic loci), which are one of the genetic bases for trait correlations together with linkage disequilibrium. Therefore, this model depends on genetic variation in pleiotropy, which was shown to be produced by interactions of pleiotropic loci with different rQTLs (i.e., differential epistasis; [START_REF] Pavlicev | Genetic variation in pleiotropy: Differential epistasis as a source of variation in the allometric relationship between long bone lengths and body weight[END_REF].

Importantly, selection on rQTLs can change genetic variation without an appreciable effect on trait means [START_REF] Cheverud | Pleiotropic effects on mandibular morphology II: Differential epistasis and genetic variation in morphological integration[END_REF][START_REF] Pavlicev | Genetic variation in pleiotropy: Differential epistasis as a source of variation in the allometric relationship between long bone lengths and body weight[END_REF]). Thus, different genotypes at rQTLs result in different G-matrices, and if they also show differences in mean fitness, selection will favor a specific G and may culminate in changes in the population covariance structure after several generations.

The most interesting result of this model is that G stretches in the direction of selection, that is, the combination of traits with the highest additive variance gets aligned with the direction of selection, increasing evolvability [START_REF] Pavlicev | Evolution of adaptive phenotypic variation patterns by direct selection for evolvability[END_REF]. This happens because the rQTL genotype with higher variance along the direction of selection has higher fitness.

The change in G is gradual across generations until the most advantageous genotype is fixed by selection (Figure 2A). Hence, this model posits that directional selection can theoretically change genetic trait integration to increase heritable variation in the direction of selection. This is different from what is predicted with a more traditional view of selection acting strictly on additive variation without epistasis, which would result in no change in genetic variance, and therefore no change in evolvability after selection when considering polygenetic traits and mutation-selection equilibrium (i.e., the genetic variation that is removed by selection is replenished by mutation; [START_REF] Hill | Rates of change in quantitative traits from fixation of new mutations[END_REF][START_REF] Burger | Prediction of the dynamics of a polygenic character under directional selection[END_REF]. However, both the model of [START_REF] Pavlicev | Evolution of adaptive phenotypic variation patterns by direct selection for evolvability[END_REF] and evolutionary simulations of the mutational matrix (Mmatrix; Jones et al. 2014, Melo andMarroig 2015) indicate that gene interactions (epistasis) are a crucial factor to allow the alignment between the M-matrix, the G-matrix and directional selection. More importantly, an increase in evolvability can happen quite fast, only after 100 years of directional selection on chipmunks in natural conditions [START_REF] Assis | Directional selection effects on patterns of phenotypic (co)variation in wild populations[END_REF]) and after 50 generations of artificial selection on mice [START_REF] Penna | The evolution of phenotypic integration: How directional selection reshapes covariation in mice[END_REF].

If climate-induced selection elicits a change in G such that more additive variation accumulates in the direction of selection, populations may be more prone to evolve adaptively to the same climatic factor in the future, even if traits are strongly integrated before the action of selection. This scenario may be even more relevant for population survival in face of rapid climate change if selection induced by climatic events is strong.

Extreme climatic events, which are historically rare but increasing in frequency, may induce strong episodes of directional selection [START_REF] Grant | Evolution caused by extreme events[END_REF]. Such selection has the potential to change not only trait means but also change the pattern of trait integration, culminating in a shift in evolvability. It is important to note that the [START_REF] Pavlicev | Evolution of adaptive phenotypic variation patterns by direct selection for evolvability[END_REF] simulations were done with very weak selection (selection coefficients of 0.001 or 0.05), and still the change in G could happen in less than 50 generations depending on the level of genetic linkage (Figure 2A,B). Yet, with very strong selection induced by extreme climatic pressures, a small change in G could happen very fast, perhaps in just a single generation, but being enough to increase evolvability in a biologically meaningful amount (Figure 2C).

We test this new hypothesis using lizard limb traits subjected to selection induced by hurricanes to demonstrate how this framework can be applied in future studies. Different from previous studies showing increased evolvability after selection [START_REF] Pavlicev | Genetic variation in pleiotropy: Differential epistasis as a source of variation in the allometric relationship between long bone lengths and body weight[END_REF][START_REF] Assis | Directional selection effects on patterns of phenotypic (co)variation in wild populations[END_REF][START_REF] Penna | The evolution of phenotypic integration: How directional selection reshapes covariation in mice[END_REF], we test this hypothesis on a system in which selection is not on the main axis of variation, which is body size here and in many other studies. It is unclear whether evolvability can increase with directional selection that is not aligned with the axis of greatest variation. In this case study, we use P-matrices as surrogates of Gmatrices, and our results suggest that extreme climatic events could lead to increased evolvability in some populations. However, we caution that the ultimate test needs to be performed with G-matrices and larger sample sizes. An analytical framework developed by [START_REF] Grabowski | How many more? Sample size determination in studies of morphological integration and evolvability[END_REF] can be used to estimate the minimum sample sizes needed to find significant changes in evolvability. Nevertheless, we emphasize that having ideal samples to study unpredictable extreme events in nature is quite difficult, and thus, we must work with the available data to reach an understanding of potential impacts of these increasingly common phenomena on biodiversity.

Case study: Hurricane-induced selection on lizard limb traits

We focused on two populations of the lizard Anolis scriptus that were exposed to Hurricanes Irma and Maria in 2017, showing strong selection on limb traits induced by this extreme climatic event [START_REF] Donihue | Hurricane-induced selection on the morphology of an island lizard[END_REF]. In both populations, toepad area became larger, relative humerus length became longer, whereas relative femur length was shorter in survivors of the hurricane. These morphological changes were shown to be associated with higher clinging capacity and reduced aerodynamic drag, performances that are likely relevant to increase survival in the face of hurricanes [START_REF] Donihue | Hurricane-induced selection on the morphology of an island lizard[END_REF][START_REF] Debaere | An aerodynamic perspective on hurricane-induced selection on Anolis lizards[END_REF]).

These populations were re-surveyed in 2019 and data on offspring of the survivors were collected, allowing for the assessment of the degree to which the evolutionary responses to selection in the offspring (∆𝑧) followed the direction of selection (𝛽). Toepad area is likely heritable because the offspring of survivors from the hurricane also had large toepads [START_REF] Donihue | Hurricane effects on Neotropical lizards span geographic and phylogenetic scales[END_REF]. Although this change could be due to phenotypic plasticity, limb traits are moderately to highly heritable. Narrow-sense heritability of hindlimb has been measured in Anolis sagrei, exhibiting quite high values (0.78 + 0.13, see [START_REF] Calsbeek | Postcopulatory fertilization bias as a form of cryptic sexual selection[END_REF]. Significant genetic variance on limb length was also measured in five species of Anolis, where heritabilities varied from 0.45 -0.54 for femur and tibia, and around 0.15 for humerus and ulna [START_REF] Mcglothlin | Adaptive radiation along a deeply conserved genetic line of least resistance in Anolis lizards[END_REF]. Therefore, the limb structures of Anolis that we are dealing with comprise heritable traits that show additive variance that can be molded by selection.

Because we know that hurricane-induced selection is multivariate and limb traits are expected to be highly integrated in lizards, we need to analyze more than just a single trait for a robust understanding of the response to selection and whether evolvability changes after this episode of strong selection. Morphological integration of limb traits has been characterized in different Anolis species, showing that phenotypic correlations are generally high among limb elements [START_REF] Kolbe | Convergent evolution of phenotypic integration and its alignment with morphological diversification in Caribbean Anolis ecomoprhs[END_REF]. Recent studies in Anolis species revealed that genetic correlations are also high between limb traits [START_REF] Mcglothlin | Adaptive radiation along a deeply conserved genetic line of least resistance in Anolis lizards[END_REF][START_REF] Mcglothlin | The Adaptive Radiation of Anolis lizards[END_REF]. High positive between-limb genetic correlations mean that it is easier to change all limb traits in the same direction (increasing or reducing all limb traits; directions that bear high genetic variation) than changing just fore-or hindlimb traits or changing fore-and hindlimb in different directions. We expect a similar pattern of high limb integration for A. scriptus populations (i.e., high positive correlations between fore-and hindlimbs), and that toepad area will be mostly correlated with the length of the corresponding longest finger or toe on which it was measured. Moreover, given that hurricane-induced selection resulted in an increase in the lengths of forelimb traits, but a reduction in the lengths of hindlimb traits [START_REF] Donihue | Hurricane-induced selection on the morphology of an island lizard[END_REF], we expect the evolutionary response to selection in the offspring to have not followed the direction of selection because fore-and hindlimbs are genetically correlated in anoles [START_REF] Mcglothlin | Adaptive radiation along a deeply conserved genetic line of least resistance in Anolis lizards[END_REF][START_REF] Mcglothlin | The Adaptive Radiation of Anolis lizards[END_REF]). Yet, we also expect to see an increase in evolvability in the offspring if hurricane-induced selection reduced correlations between the limbs (i.e., if selection reduced the overall integration of limb traits due to the opposing direction of selection on fore-and hindlimbs).

Methods

Sampling and morphological measurements

Anolis scriptus populations were sampled on two small islands in the Turks and Caicos: Pine Cay and Water Cay (350 and 250 ha, respectively). The islands are presently connected by a sandbar that was deposited by hurricane storm surge in the mid-1990s. However, the two populations are genetically distinct (unpublished data), and therefore can be treated as replicates of the response to extreme climatic events. Both cays are similarly covered by low vegetation (1 to 3 m in height). We sampled the anole populations along a 2-km transect on each island. The initial survey began on 28 August 2017 and lasted until 4 September 2017.

We did not permanently mark individuals on the islands, and consequently, direct measures of individual survival are impossible. Hurricane Irma hit the Turks and Caicos on 8 September 2017 and Hurricane Maria passed through the region on 22 September 2017.

Populations were revisited after the hurricanes between 16 October and 20 October 2017 and again approximately 18 months after the storms between 1 April 2019 and 8 April 2019. We caught 71 adult lizards in the initial survey, 93 during the 2017 revisit, and 117 lizards in 2019. This difference in sample size before and after the hurricane reflects an increased sampling effort after the hurricane. However, we do not expect these differences in sampling to strongly affect selection differentials and response to selection because these metrics are calculated using sample means of surveyed populations, and means are generally well estimated with just a few individuals (although standard errors are generally higher with smaller sample sizes).

To identify the most likely offspring of hurricane survivors in 2019, we calculated the maximum body size of a lizard that had hatched on 1 April 2018, one year before our revisit (following [START_REF] Donihue | Hurricane effects on Neotropical lizards span geographic and phylogenetic scales[END_REF]). To do this, we used a logistic-by-length model and estimated growth rates in ecologically similar species to characterize the growth curve of lizards (see Donihue et al. 2020 for details). By excluding lizards larger than this size threshold, we conservatively estimated the subset of the sampled population less than a year old and thus, offspring of hurricane survivors. We retained for further analyses only the sexually mature adults of each sampling period (i.e., pre-hurricane / post-hurricane / offspring) by using an adult body size threshold for each sex. See Supplementary Material and Methods for a full description of the body size thresholds used. Our final sample sizes for each population were 36, 53 and 57 for Pine Cay, and 39, 52 and 57 for Water Cay, for pre-hurricane, posthurricane, and offspring surveys, respectively.

All morphological measurements were taken by the same researcher on all three visits. Measurements included the snout-vent length (SVL), eight limb length segments, and toepad area on the fore-and hindlimbs (as described in [START_REF] Herrel | Rapid large-scale evolutionary divergence in morphology and performance associated with exploitation of a different dietary resource[END_REF]; see also [START_REF] Lowie | Do the relationships between hind limb anatomy and sprint speed variation differ between sexes in Anolis lizards[END_REF]. For further details on the morphological measurements, see Supplementary Material.

Selection analysis and response to selection

To measure total hurricane-induced selection (indirect plus direct selection), we estimated a proxy of selection differentials (S), the mean trait difference between post-and pre-hurricane sampling periods (S = [mean traits of post] -[mean traits of pre]), which represent the survivors after selection and the population before selection, respectively. The selection differential should be a difference between the mean of parents of the next generation, after the hurricane, and the mean of pre-hurricane populations, thus we are assuming that most survivors reproduced. We calculated separate selection differentials for each island, given that the populations are genetically distinct, and they differed in mean trait values in the pre-hurricane survey (Table S1). To measure only the direct selection on traits, we estimated the corresponding selection gradients (β) for each island. We used two different equations: 𝛽 = 𝑃 -1 𝑆 and 𝛽 = 𝑃 -1 𝑧 [START_REF] Lande | The measurement of selection on correlated characters[END_REF], where P -1 is the inverse of the phenotypic covariance matrix of limb traits in pre-hurricane populations and z is the response to selection. The second equation estimates 𝛽 as the realized selection because it is based on the realized response to selection. We calculated the response to selection (z) for each island as the mean difference in traits between offspring and pre-hurricane populations.

Inverting the P-matrix discounts the effects of trait correlations from the selection vector, retaining just direct selection on traits [START_REF] Lande | The measurement of selection on correlated characters[END_REF]). However, given that Pmatrices are estimated with sampling error, their inversion makes more error accumulate in the smallest eigenvalues (which dominate the inverted P-matrix), biasing estimates of selection gradients. To correct for this bias, we performed a sampling noise control approach following [START_REF] Marroig | Modularity, noise, and natural selection[END_REF], in which we replace the smallest unreliable eigenvalues of Pmatrices by the last reliable eigenvalue (which is the sixth one for the A. scriptus P-matrices, determined using the second derivative of the distribution of eigenvalues). This procedure has been shown to improve estimation of selection gradients (for more details, see [START_REF] Marroig | Modularity, noise, and natural selection[END_REF]. We explain below how we estimated P-matrices. Finally, we measured the vector correlation between both estimates of 𝛽 to check if they are similar, and thus ensure that our assumption of post-hurricane population comprising mostly parents of the next generation holds true. To determine the significance of vector correlations, we randomly sampled 1,000 vectors of 11 elements (number of traits) from a normal distribution of mean zero and unit SD and calculated their vector correlations. The critical values considering  = 0.05 for this random distribution is -0.57 and 0.57. Therefore, any vector correlation below or above these values were considered significant.

To estimate uncertainty due to sampling error in our selection analysis, we performed bootstrapping of individuals caught on each survey in each island. We resampled with replacement 1,000 times for each population and then we calculated the selection metrics on resampled populations, resulting in 1,000 values of S and β. We extracted 95% confidence intervals (CI) for the selection metrics and if they did not encompass zero, the values were considered statistically significant. We also calculated 95% CI for mean-standardized selection to determine if hurricane-induced selection is strong compared to the literature [START_REF] Hereford | Comparing strengths of directional selection: how strong is strong?[END_REF], and for z values to determine if the offspring showed significant changes.

P-matrices, overall integration and modularity of limb traits

We estimated P-matrices for each sampling period on each island to use as surrogates for Gmatrices. Environmental effects leading to phenotypic plasticity may cause the P-matrix to differ from the G-matrix [START_REF] Willis | Can One Predict the Evolution of Quantitative Characters Without Genetics?[END_REF]), but the potential for plasticity to impact traits in this study is likely low because growth after sexual maturity is limited. Moreover, Cheverud's conjecture has shown that well estimated G-matrices tend to resemble P-matrices [START_REF] Cheverud | A comparison of genetic and phenotypic correlations[END_REF], especially for morphological traits with generally high heritabilities (i.e., high contribution of genetic factors to overall variance), which is the case for lizard limb traits. Moreover, empirical examples of similarity between G and P are mounting in the literature (e.g., [START_REF] Roff | The evolution of genetic correlations: An analysis of patterns[END_REF][START_REF] Phillips | Hierarchical Comparison of Genetic Variance-Covariance Matrices. I. Using the Flury Hierarchy[END_REF][START_REF] Marroig | A comparison of phenotypic variation and covariation patterns and the role of phylogeny, ecology, and ontogeny during cranial evolution of New World monkeys[END_REF][START_REF] Sodini | Comparison of Genotypic and Phenotypic Correlations: Cheverud's Conjecture in Humans[END_REF][START_REF] Styga | Development of G: a test in an amphibious fish[END_REF]. Finally, we compared A. scriptus P-matrices with a G-matrix of the closest relative, A. cristatellus, extracted from [START_REF] Mcglothlin | Adaptive radiation along a deeply conserved genetic line of least resistance in Anolis lizards[END_REF], and matrices are substantially similar (see Supplementary Methods and Tables S2,S3). We also confirmed that P-matrices are reasonably well-estimated using a rarefaction analysis and we compared all P-matrices across populations of both islands (see Supplementary Methods and Figure S1). Given that P = G + E [START_REF] Falconer | Introduction to Quantitative Genetics[END_REF], if P-matrices are similar, we may assume that P and G-matrices are also similar, otherwise the environmental matrix (E) in all populations would have to compensate for differences in G-matrices, which is a very unlikely scenario [START_REF] Marroig | A comparison of phenotypic variation and covariation patterns and the role of phylogeny, ecology, and ontogeny during cranial evolution of New World monkeys[END_REF]. Even though it is preferable to use estimates of Gmatrices to test the hypothesis of increased evolvability with strong selection induced by climatic extremes, because the potential to adapt ultimately depends on heritable additive variation, estimating G-matrices in wild populations with appropriate sample sizes may be difficult. Despite this, the growing use of genomic data allowing the estimation of more precise G-matrices in a wide variety of organisms [START_REF] Gienapp | Genomic quantitative genetics to study evolution in the wild[END_REF]) will enable researchers to test this hypothesis more rigorously in the future.

To estimate mean overall integration of traits in both Pine and Water Cays, we used correlation P-matrices and calculated an integration index using the variance of eigenvalues (egv), which is less prone to bias due to differences in sample size compared to other integration metrics if negligible eigenvalues are excluded (O' Keefe et al. 2022). This index is a ratio between the variance of eigenvalues and its theoretical maximum, which simply corresponds to the number of traits minus one [START_REF] Pavlicev | Measuring morphological integration using eigenvalue variance[END_REF]). To test if egv values differed between surveys, we used the 1,000 bootstrapped populations to calculate 1,000 Pmatrices for each survey and extract a 95% CI of overall integration for each survey. If the 95% CI did not overlap between-surveys, we considered the differences in egv significant.

We performed this analysis for all traits (including SVL), and separately for forelimbs and hindlimbs (excluding SVL) to check if one of the limbs showed a stronger response to hurricane-induced selection than the other.

We also performed a modularity analysis to explore the extent of integration of foreand hindlimb components in relation to between-limb integration, and whether this pattern changed after the hurricanes. If the selection induced by the hurricane reduced the magnitude of between-limb correlations in the P-matrix, we expect to see an increase in limb modularity.

We constructed three hypothetical modules: forelimb, hindlimb and toepad, the latter being composed of toepad area of fore-and hindlimbs as well as the longest finger and longest toe.

Even though toepads are epidermal modifications that likely develop from an embryonic origin semi-independent from the skeletal limb elements [START_REF] Weatherbee | Chapter 7. Mechanisms of Chondrogenesis and Osteogenesis in Limbs[END_REF][START_REF] Liu | Gekko japonicus genome reveals evolution of adhesive toe pads and tail regeneration[END_REF][START_REF] Gamble | Duplications in Corneous Beta Protein Genes and the Evolution of Gecko Adhesion[END_REF], toepads are functionally linked to toes and fingers in locomotion [START_REF] Ruibal | The structure of the digital setae of lizards[END_REF][START_REF] Elstrott | Evolutionary correlations among morphology, habitat use and clinging performance in Caribbean Anolis lizards[END_REF], and therefore may constitute a functional module together with toes and fingers. We then tested whether average correlations within hypothetical modules (AVG+) were higher than between-modules (AVG-) using the empirical correlation P-matrices and a Mantel test to control for non-independence between correlations [START_REF] Cheverud | Morphological Integration in the Saddle-Back Tamarin (Saguinus fuscicollis) Cranium[END_REF][START_REF] Porto | Size variation, growth strategies, and the evolution of modularity in the Mammalian skull[END_REF]; using function 'TestModularity' in evolqg R package, Melo et al. 2015). Although the Mantel test is effective in detecting biologically meaningful modules when the modularity hypotheses are grounded in wellestablished developmental and functional studies as is our case here, we also tested modularity using a more recent test, the Covariance Ratio coefficient (CR), and permutations to calculate significance of effect sizes using the R package geomorph (Adams and[START_REF] Adams | geomorph: an R package for the collection and analysis of geometric morphometric shape data[END_REF][START_REF] Adams | Comparing the strength of modular signal, and evaluating alternative modular hypotheses, using covariance ratio effect sizes with morphometric data[END_REF]. The CR is based on a ratio of within-and betweenmodule covariances and gives standardized effect sizes (using Z-score transformations that take the scale-dependence of traits into account) that characterize the strength of the modular signal [START_REF] Adams | Comparing the strength of modular signal, and evaluating alternative modular hypotheses, using covariance ratio effect sizes with morphometric data[END_REF].

Evolvability analysis

Evolvability measures the amount of genetic variation in G in the direction of selection (𝑒(𝛽) = (𝛽 𝑇 𝐺𝛽)/|𝛽| 𝟐 , where T means transpose and |𝛽| 𝟐 is the squared strength of selection; [START_REF] Hansen | Measuring and comparing evolvability and constraint in multivariate characters[END_REF]). As we mentioned before, we substituted P-matrices for G-matrices to estimate evolvabilities. Thus, our interest was to compare evolvability in the directions of 𝛽 before and after the hurricane to see if the P-matrices increased the amount of variation in the direction of selection, and to compare evolvability in the directions of 𝛽 and 𝑧. If any of these vectors were aligned with trait correlations (i.e., aligned with pmax, the phenotypic equivalent of gmax), we expected their evolvability to be comparably high. To estimate the evolvability in the direction of selection, we then used the equation 𝑒(𝛽) = 𝛽 𝑇 𝑃𝛽 (e.g., [START_REF] Assis | Directional selection effects on patterns of phenotypic (co)variation in wild populations[END_REF]; using function 'Evolvability' from evolqg R package, Melo et al., 2015) with β normalized to unit length and covariance P-matrices. Final 𝑒(𝛽) were standardized using the geometric mean of all traits to account for scale differences between traits and populations (e.g., [START_REF] Assis | Directional selection effects on patterns of phenotypic (co)variation in wild populations[END_REF]. We chose to standardize using a single value, the average geometric mean of traits, instead of using the mean of each trait (e.g., [START_REF] Hansen | Measuring and comparing evolvability and constraint in multivariate characters[END_REF] to preserve the structure of covariation.

We calculated two values of evolvability in the direction of selection for each island:

(1) 𝑒(𝛽) before the hurricane: using pre-hurricane P-matrices and [postpre] selection gradients, and ( 2) 𝑒(𝛽) after the hurricane: using offspring P-matrices and [postpre] selection gradients. To estimate uncertainty in evolvability, we once again used the 1,000 bootstrapped values of 𝛽 and the 1,000 bootstrapped P-matrices to construct 95% CI. This procedure accounts for error in estimation of both selection and P-matrices. We interpreted values as significantly different when 95% of the distributions did not overlap.

All analyses were conducted in the R programming environment (R Core Team 2021, version 4.1). We provide R scripts showing the selection and evolvability analyses, and showing the calculations to estimate inaccuracy in evolvability and number of individuals needed to achieve a certain level of accuracy in the Supplementary Materials.

Results

Hurricane-induced direct selection on limb traits and response to selection

Selection (S) was mostly negative on traits, however, we only found significant negative selection on the hindlimb longest toe in Pine Cay (Table 1). In contrast, we found that several selection gradients (β) were significant on both island populations-mostly negative in SVL and hindlimb traits and positive for toepads. We also found positive direct selection on the humerus, but it was only significant on Water Cay (Table 1). Results were similar when using β vectors estimated with either S or 𝑧, given that both ways of calculating β resulted in significantly similar vectors (i.e., vector correlations of 0.86 for Pine Cay and 0.93 for Water Cay, both above the critical value of 0.57). Moreover, selection was parallel in both islands, showing a vector correlation of 0.85. Interestingly, the mean-standardized β vectors are similar between-islands in magnitude, ranging from weak selection (around 0.05) to strong selection (around 0.7) when compared to the average mean-standardized β (0.89) found in [START_REF] Hereford | Comparing strengths of directional selection: how strong is strong?[END_REF] for morphological traits (see their Table 1). Response to selection (𝑧) showed mostly a decrease in lengths of all traits, however, significant changes were only found in Pine Cay (Table 1). 𝑧 vectors were not similar between populations, showing a non-significant vector correlation of 0.49. While the response to selection was aligned with β in Water Cay (vector correlation of 0.84), this was not the case for Pine Cay (vector correlation of 0.48), even though both β vectors are not correlated with the axis of highest phenotypic variance (vector correlations of 0.0 and 0.05, respectively). These results indicate that direct selection was Table 1. Selection differentials (S), selection gradients (), and response to selection (z) separated by island. We show empirical values and the bootstrapped 95% CI. CI in bold are significant.

Table 2. Matrix similarity between all P-matrices. We used two methods to compare Pmatrices. Random Skewers (RS) indicates the similarity in the response to random selection, whereas Krzanowski correlation (Krz) indicates similarity based on the best matched pairs of eigenvectors. All values are significant at  = 0.05. Table 3. Empirical evolvabilities and associated 95% confidence intervals.

Discussion

Incorporating multivariate adaptive responses of populations to climate change is a necessary step to better predict how and which species may overcome these challenges (Peschel et al. 2020). However, if climate-induced selection acts on multiple traits simultaneously and if traits are highly integrated, the response to selection may be constrained, not closely following the direction of selection [START_REF] Etterson | Constraint to Adaptive Evolution in Response to Global Warming[END_REF][START_REF] Hellmann | Constraints and reinforcement on adaptation under climate change: Selection of genetically correlated traits[END_REF][START_REF] Chevin | Genetic constraints on adaptation to a changing environment[END_REF]. In this scenario, populations may not be able to cope with climate change even when responding to selection because the response would lead to a lower increase in mean fitness given that some traits will respond non-adaptively. This is especially likely if genetic correlations between positively selected traits are negative (e.g., [START_REF] Etterson | Constraint to Adaptive Evolution in Response to Global Warming[END_REF], but also if selection is not aligned with axes accumulating the most variance in the G-matrix (Arnold et al. 2001, Blows and[START_REF] Walsh | Abundant Genetic Variation + Strong Selection = Multivariate Genetic Constraints: A Geometric View of Adaptation[END_REF]). Yet, if selection also changes the pattern of trait integration and increases the amount of genetic variance aligned with selection, future response to this same selection may be facilitated [START_REF] Pavlicev | Evolution of adaptive phenotypic variation patterns by direct selection for evolvability[END_REF]).

These alternative scenarios ultimately depend on the extent to which trait integration results in only a few directions accumulating most of the genetic variance, and how much selection is aligned with these directions (showing high evolvability; scenario B in Figure 1) or whether selection changes trait integration to increase evolvability.

We investigated these possibilities in a case study of an extreme climatic event, hurricane-induced selection on lizard limb traits [START_REF] Donihue | Hurricane-induced selection on the morphology of an island lizard[END_REF][START_REF] Donihue | Hurricane effects on Neotropical lizards span geographic and phylogenetic scales[END_REF], in which we had samples of survivors and their offspring for two different populations. Given that selection was inferred to be strong on some of the limb traits [START_REF] Donihue | Hurricane-induced selection on the morphology of an island lizard[END_REF], we hypothesized that selection would change the pattern of limb integration and re-distribute genetic variation in a way that would increase the chance of adaptation to hurricanes in the future. Because we could not estimate G-matrices for the populations, our results are only suggestive, but indicate that both populations had higher evolvability in the direction of selection after the hurricane, although this result was not significant. This lack of significance may be more related to a low sample size than to a genuinely non-significant effect, as we calculated that a sample size of around 70 individuals with the same evolvability values would lead to enough accuracy to find significant difference in evolvability before and after the hurricane.

Interestingly, our results raise the possibility that evolvability may increase after a single episode of selection and even if the direction of selection is not aligned with size (or PC1, the direction accumulating most genetic variance). Still, using P-matrices as surrogates of Gmatrices likely may have overestimated evolvability values if environmental variation contributes to total variation in the direction of selection. Alternatively, environmental variation may be much higher along pmax than gmax given that this direction is related to size, and several environmental factors can affect growth. In this case, evolvabilities can be underestimated by using P-matrices instead of G-matrices.

The model of directional selection influencing evolvability predicts changes in the pattern of correlation even with fluctuating selection, but the change depends on the level of epistasis between rQTLs and other genes [START_REF] Pavlicev | Evolution of adaptive phenotypic variation patterns by direct selection for evolvability[END_REF]). Hence, it may be that only sustained or frequent directional selection, such as maintaining the same direction for several generations, would be able to increase evolvability in the direction of selection. Although we do not expect hurricane-induced selection to be this persistent because of its infrequency, lizard populations experiencing periodic hurricanes (up to four hurricanes over a 70-year period) show larger toepads than populations exposed to fewer hurricanes [START_REF] Donihue | Hurricane effects on Neotropical lizards span geographic and phylogenetic scales[END_REF], indicating that hurricanes can leave persistent marks on phenotypes. Moreover, hurricane frequency and magnitude are predicted to increase in the future due to climate change [START_REF] Bender | Modeled Impact of Anthropogenic Warming on the Frequency of Intense Atlantic Hurricanes[END_REF], Knutson et al. 2010).

Direction of selection, locomotor performance and variation in responses to selection

Direct selection on limb traits in Pine Cay and Water Cay populations was parallel, giving some confidence in their estimation with the sampling noise control approach when inverting P-matrices [START_REF] Marroig | Modularity, noise, and natural selection[END_REF]. Parallel selection was also found in other lizard systems inferred to have suffered hurricane-induced selection (e.g., [START_REF] Rabe | Post-hurricane shifts in the morphology of island lizards[END_REF]. Direct selection was strong in several limb traits, but not so strong on toepad area, although selection gradients were towards increased surface area, and this is likely related to increased clinging capacity in A. scriptus [START_REF] Crandell | Stick or grip? Coevolution of adhesive toepads and claws in Anolis lizards[END_REF][START_REF] Donihue | Hurricane-induced selection on the morphology of an island lizard[END_REF]). However, shifts in toepad area do not always lead to higher clinging capacity in lizards (e.g., [START_REF] Dufour | Parallel increases in grip strength in two species of Anolis lizards after a major hurricane on Dominica[END_REF]) and other traits beyond toepad area also change clinging capacity (such as setae length and limb length; [START_REF] Hagey | Variation in setal micromechanics and performance of two gecko species[END_REF][START_REF] Kolbe | Effects of Hind-Limb Length and Perch Diameter on Clinging Performance in Anolis Lizards from the British Virgin Islands[END_REF]. Direct selection to reduce femur length was stronger, supporting the evidence in favor of reduced aerodynamic drag in lizards that have smaller hindlimbs and consequently exposing less of the pelvic region to the strong winds induced by hurricanes [START_REF] Debaere | An aerodynamic perspective on hurricane-induced selection on Anolis lizards[END_REF]. Thus, selection induced by hurricanes, and possibly selection induced by other climatic events, may vary in direction and strength on different traits, showing the relevance of using a multivariate approach to estimate selection.

Another source of variation induced by climatic events is how much the responses are constrained by trait integration. Fore-and hindlimbs compose a single integrated unit in both populations. This high between-limb integration was expected from previous studies [START_REF] Kolbe | Convergent evolution of phenotypic integration and its alignment with morphological diversification in Caribbean Anolis ecomoprhs[END_REF][START_REF] Mcglothlin | Adaptive radiation along a deeply conserved genetic line of least resistance in Anolis lizards[END_REF]) and likely results from the need to coordinate both limbs in locomotion [START_REF] Aerts | Lizard locomotion: how morphology meets ecology[END_REF]. However, populations may still differ in the strength of between-limb integration. For instance, the modularity analysis showed that between-limb correlations were lower in Water Cay than in Pine Cay, and the former had a morphological response that followed much more closely the direction of selection than Pine Cay (vector correlation of 0.84 compared with 0.48, respectively). Because direct selection was not aligned with limb integration in either of the populations (i.e., selection did not increase or decrease all limb traits, as evidenced by the lower correlation between  and pmax), the lower correlations between limbs may have resulted in less correlated responses of limb elements in Water Cay, and therefore, a lower constraint in the direction of the response to selection. Finally, we must also consider that populations may differ due to chance and evolutionary forces unrelated to selection and having more than two populations in future studies will aid in evaluating the consistency of patterns of selection and responses to selection induced by hurricanes.

Therefore, although extreme climatic events like hurricanes are predicted to become more frequent and intense with future climate change [START_REF] Bender | Modeled Impact of Anthropogenic Warming on the Frequency of Intense Atlantic Hurricanes[END_REF], Knutson et al. 2010), the responses of lizard populations may vary because of differences in strength of selection, and differences in trait integration resulting in different degrees of constrained responses to selection. A broader implication of all these sources of variation for responses to climate-driven selective pressures is the resultant variation in how long these changes last in terms of evolutionary time. For traits experiencing stronger selection and having lower genetic correlations, a longer-lasting evolutionary response to extreme climatic events may be more common, as shown for populations and species more often hit by hurricanes showing larger toepads [START_REF] Donihue | Hurricane effects on Neotropical lizards span geographic and phylogenetic scales[END_REF]. On the other hand, for traits that show a more constrained response because of high trait integration, the evolutionary responses may be more ephemeral.

Selection on trait integration patterns and evolvability

Even though we found evidence for strong selection on some limb traits in both populations, integration patterns were not changed by this selection, suggesting that only very strong selection (or lesser magnitude of selection sustained over multiple generations) can redistribute variation in the G-matrix. For instance, in two different species of chipmunks that experienced selection induced by 100 years of climate change, only the species under stronger selection showed increased evolvability after selection [START_REF] Assis | Directional selection effects on patterns of phenotypic (co)variation in wild populations[END_REF]. Whether strong selection changes trait integration, and perhaps also evolvability, may be related to a conflict with multivariate stabilizing selection (Cheverud 1984, Roff andFairbain 2012) associated with maintenance of a pattern of higher correlations among traits that perform the same function. In this scenario, the G-matrix is kept evolutionary conserved by a stable pattern of correlational selection, that in the case of lizard limbs would be associated with locomotor tasks normally performed in the environment (e.g., sprinting, jumping, climbing; [START_REF] Irschick | A comparative analysis of the ecological significance of maximal locomotor performance in Caribbean Anolis lizards[END_REF][START_REF] Perry | Locomotor performance and social dominance in male Anolis cristatellus[END_REF][START_REF] Foster | How forelimb and hindlimb function changes with incline and perch diameter in the green anole, Anolis carolinensis[END_REF]. These performances demand a coordination between fore-and hindlimbs, and thus selection to enhance these performances would favor high integration between limbs [START_REF] Mcglothlin | Adaptive radiation along a deeply conserved genetic line of least resistance in Anolis lizards[END_REF]. This conflict between episodes of directional selection causing short-term, fast changes in G, and multivariate stabilizing selection promoting long-term evolutionary stability in G, must be common in nature given that different types of selection may act on populations [START_REF] Arnold | Understanding the evolution and stability of the G-matrix[END_REF]. Selective forces and other factors, such as mutation, all influence simultaneously the G-matrix and the balance between these factors ultimately shape the evolutionary potential of populations. Therefore, hurricane-induced selection on A. scriptus limbs in the Turks and Caicos appears to conflict with selection associated with normal environmental conditions with the latter being much more persistent. Hence, the P-and G-matrices of A. scriptus are likely shaped by multivariate selection related to normal limb functioning to a much greater extent than by hurricane-induced selection, even if these extreme selection events are much stronger than selection in normal ecological conditions. Even for populations subjected to hurricanes in a higher frequency, selection in normal ecological conditions may erase any potential long-lasting signal of lower limb integration due to hurricane-induced selection.

Conclusion

We showed how the hypothesis of directional selection leading to increased evolvability can be studied in the context of population-wide responses to extreme climatic events. Ideally, this hypothesis should be tested with G-matrices with appropriate sample sizes before and after the extreme event, which are likely above 70 individuals if studying around 10 traits.

Although we expected strong hurricane-induced selection to reduce limb trait integration, this did not happen in A. scriptus populations, potentially because the episodic selection caused by this hurricane was not persistent enough to contrast normal ecological conditions, despite hurricanes being a relatively strong selective force. Even so, our results indicate that evolutionary responses follow more closely the direction of selection if between-limb integration is lower. Therefore, persistence of populations facing strong environmental changes may be enhanced by lower trait integration, even if individuals perform normal functions at a lower level than more highly integrated individuals in normal ecological conditions. most variance accumulating in the direction of gmax. In this scenario, the evolutionary response may not follow the direction of selection, representing a case of genetic constraint, as shown for population 1 (Pop 1). However, if selection is aligned with gmax (Pop 2) or with PC2 (Pop 3), the evolutionary response will follow the direction of selection even with high genetic integration between traits. A set of traits is considered a module when within-set average correlation is higher than between-set average correlation. See table S2 for the results of the modularity tests.

  Figure adapted from Arnold et al. 2001.

Figure 2 .

 2 Figure 2. Conceptual framework of directional selection shifting evolvability. A) A

Figure 3 .

 3 Figure 3. Mean overall integration of limb traits and toepad area. The values are

Figure 4 .

 4 Figure 4. Average correlations within (w) and between (b) the three hypothetical modules.

parallel across islands, not aligned with the direction of most phenotypic variance, and the response in the offspring differed between the two islands.

Pattern of limb trait correlations and overall integration

Correlation P-matrices are shown in Table S4. All correlation P-matrices show a PC1 that is predominantly isometric size, with similar magnitudes of loadings for all traits and all loadings with the same sign (Table S5). Also, PC1 explains high amounts of proportional variation in all surveys, showing that variance accumulates in the direction of size. However, the prehurricane population in Water Cay shows less proportional variation accumulated in PC1 compared with Pine Cay. All P-matrices show a very high similarity when using Random Skewers (RS), with all values above 0.91, indicating that matrices respond similarly to random selection. Similarity of the best matched PC axes (Krzanowski correlation, Krz) is in general lower, ranging from 0.70 to 0.79, which indicates that some PCs diverge across populations (Table 2). Still, P-matrices can be considered quite similar across populations, giving confidence in our assumption that P is a good surrogate of G.

We found a tendency for overall integration to be higher in post-hurricane populations, but the integration distributions overlap considerably (Figure 3A). The same pattern occurs if separating forelimb (95% CI for Pine Cay: pre = 0.21, 0.50; post = 0.36, 0.61; offspring = 0.19, 0.42; for Water Cay: pre = 0.12, 0.45; post = 0.25, 0.52; offspring = 0.18, 0.45; Figure 3B) and hindlimb (95% CI for Pine Cay: pre = 0.38, 0.73; post = 0.47, 0.73; offspring = 0.30, 0.60; for Water Cay: pre = 0.16, 0.5; post = 0.35, 0.65; offspring = 0.28, 0.60; Figure 3C). Results are the same if keeping all eigenvalues or just the first six eigenvalues. Therefore, mean limb trait integration did not change with the hurricane-induced selection.

Modularity of limb traits

The hindlimb is the only hypothetical module that shows higher average correlations withinmodule than between-modules, especially for Pine Cay populations (Figure 4). However, the hindlimb showed a signal of modularity that is significant only at P = 0.10 (Table S6).

Therefore, forelimb, hindlimb and toepad traits can be seen as a single integrated unit. Still, we noticed that between-limb integration in pre-hurricane population was lower in Water Cay than in Pine Cay for both limbs (Figures 4B and4C).

When using the CR coefficient, we found that both pre-hurricane populations have a significant CR coefficient. However, CR values are above 1.0 indicating that covariances within-modules are lower than between-modules. Moreover, we found no differences in pairwise comparisons of CR effect sizes across populations (Table S7, Figure S2). Therefore, we confirmed our results from the Mantel test that the limbs are a single integrated unit in all populations and that the pattern did not change after the hurricane.

Evolvability

We found a similar trend of increased evolvability in the direction of β from pre-hurricane parents to offspring in both islands (Table 3). However, the simulated confidence intervals are quite large, making these differences non-significant. This lack of significance may be due to low sample size to detect an evolvability difference. To determine the level of accuracy in our evolvability values, we used the analytical framework developed by [START_REF] Grabowski | How many more? Sample size determination in studies of morphological integration and evolvability[END_REF], in which we could calculate the inaccuracy of evolvability given our sample size (lowest in pre-hurricane populations, around 35 individuals) and the overall integration of P-matrices (from 0.5 to 0.6). By looking at Figure S3A, our inaccuracy in evolvability ranged from 0.05 to 0.06. We also determined how many individuals we would need to find a significant difference (accuracy of 0.03 or lower) given the overall integration of P-matrices: we would need at least 65 to 75 individuals (Figure S3B).