Marie Bouih
email: bouih@ipgp.fr

Isabelle Panet

Dominique Remy

Laurent Longuevergne

Sylvain Bonvalot

Deep mass redistribution prior to the Mw 8.8 Maule Earthquake (Chile) revealed by GRACE satellite gravity

Keywords: Gravity gradients, GRACE, Earthquake, Signal separation

published or not. The documents may come

Deep mass redistribution prior to the 2010 Mw 8.8 Maule (Chile) Earthquake revealed by GRACE satellite gravity

Introduction

The February 27 th , 2010, M w 8. [START_REF] Chao | Gravity Changes Due to Large Earthquakes Detected in GRACE Satellite Data via Empirical Orthogonal Function Analysis[END_REF] Maule earthquake is one of the largest instrumental earthquakes instrumentally recorded earthquakes. It nucleated in the central region of the historic 1835 Concepcion event (Mw 8.5), matching a zone of high coupling previously characterized as a mature seismic gap [START_REF] Ruiz | Historical and recent large megathrust earthquakes in Chile[END_REF]. This event ruptured a 500-km length segment of the interface between the downgoing Nazca and the over-riding South American plates, releasing stresses accumulated over more than 175 years since the last M w 9 earthquakes in 1730 and 1751 [START_REF] Udias | The Large Chilean Historical Earthquakes of 1647, 1657, 1730, and 1751 from Contemporary Documents[END_REF]. It produced up to 7-12 meters of thrust slip in the 24-35km depth range. The largest slip (∼ 16m) was found in the northern portion of the ruptured zone, where a M w 7.7 earthquake occurred in 1928 [START_REF] Lay | Teleseismic inversion for rupture process of the 27 February 2010 Chile (Mw 8.8) earthquake[END_REF] [START_REF] Moreno | Maule earthquake slip correlates with pre-seismic locking of Andean subduction zone[END_REF].

The slip distribution of the 2010 Maule earthquake has been derived from seismological records [START_REF] Lay | Teleseismic inversion for rupture process of the 27 February 2010 Chile (Mw 8.8) earthquake[END_REF], tsunami data and space geodetic observations [START_REF] Lorito | Limited overlap between the seismic gap and coseismic slip of the great 2010 Chile earthquake[END_REF] [START_REF] Delouis | Slip distribution of the february 27, 2010 mw= 8.8 maule earthquake, central chile, from static and high-rate gps, insar, and broadband teleseismic data[END_REF]. While the seismological data are sensitive to the propagation of the rupture during the event, space geodesy detects the surface motions osets after the rupture and their slow post-seismic variations, continuously over years or decades in the case of the GNSS Global Navigation Satellite Systems (GNSS). Based on these two types of observations, the earthquake slip distribution is however not fully constrained at depth (e.g. [START_REF] Lin | Coseismic and postseismic slip associated with the 2010 Maule Earthquake, Chile: Characterizing the Arauco Peninsula barrier eect: CHARACTERIZING ARAUCO BARRIER EFFECT[END_REF]) and the post-seismic deformation processes remain debated, from localized afterslip in the seismogenic zone to viscous ow in the mantle. Ambiguities result in particular from a limited spatial distribution of the ground stations, mostly on land, and an imperfect knowledge of material properties of the Earth.

At medium spatial scales, co-seismic and post-seismic mass redistributions have been detected by the Gravity Recovery and Climate Experiment (GRACE) satellites. This mission measured the space-time variations of the Earth's gravity eld with a decadal to monthly temporal resolution and 250-400km spatial resolution from 2002 to 2017 [START_REF] Tapley | GRACE Measurements of Mass Variability in the Earth System[END_REF]. As for other giant ruptures monitored by GRACE, a co-seismic dipole marked by a predomi-nant gravity decrease on the continental side of the subduction was observed for the Mw 8.8

Maule earthquake [START_REF] Heki | Coseismic gravity changes of the 2010 earthquake in central Chile from satellite gravimetry: COSEISMIC GRAVITY CHANGES IN CHILE[END_REF] [15] [START_REF] Dai | Improved source parameter constraints for ve undersea earthquakes from north component of GRACE gravity and gravity gradient change measurements[END_REF]. It was followed by shorter and long-term post-seismic longand short-term post-seismic signals featuring a slow gravity increase around the trench [START_REF] Tanaka | Long-and short-term postseismic gravity changes of megathrust earthquakes from satellite gravimetry[END_REF] [START_REF] Chao | Gravity Changes Due to Large Earthquakes Detected in GRACE Satellite Data via Empirical Orthogonal Function Analysis[END_REF]. The homogeneous spatial coverage of satellite gravity provided key additional information in order to constrain the geometric parameters of the ruptured fault and its average slip [START_REF] Wang | Coseismic slip of the[END_REF] [START_REF] Dai | Improved source parameter constraints for ve undersea earthquakes from north component of GRACE gravity and gravity gradient change measurements[END_REF] and to discuss the nature of the post-seismic processes [START_REF] Tanaka | Long-and short-term postseismic gravity changes of megathrust earthquakes from satellite gravimetry[END_REF] [START_REF] Han | Regional gravity decrease after the 2010 Maule (Chile) earthquake indicates large-scale mass redistribution[END_REF]. Recently, regional-scale gravity variations have been detected in the GRACE geoids in the months before the giant rupture of the March 2011 Mw 9.1 Tohoku-Oki earthquake.

They have been attributed to slab deformation at mid-upper mantle depth, eventually leading to the seismic slip as the deeper motion propagated towards the surface [START_REF] Panet | Migrating pattern of deformation prior to the Tohoku-Oki earthquake revealed by GRACE data[END_REF]. These results have been corroborated by independent GNSS data exhibiting regional crustal deformations of a few millimeters from October 2010 to March 2011, which have been related to slab extension prior to the earthquake, near 50-100km depth [START_REF] Bedford | Months-long thousand-kilometrescale wobbling before great subduction earthquakes[END_REF]. Thus, geodesy and gravity open new ways to analyze the subduction process from depth to surface, including the occurrence of giant ruptures. The unique sensitivity of satellite gravity to deeper mass redistributions oered an information complementary to the surface displacements in order to monitor aseismic motions at all depths in the subduction system. Deciphering these tenuous solid Earth signals in the gravity eld variations however requires to resolve a separation challenge: the GRACE data integrate the total gravity change induced by all the mass variations within the near-surface uid layer and the solid Earth. We need to decipher the signals from dierent sources such as hydrological, atmospheric and oceanic mass variability or viscoelastic Earth deformation from post-glacial rebound [START_REF] Boergens | Assessment of the capabilities of the temporal and spatiotemporal ICA method for geophysical signal separation in GRACE data[END_REF] [START_REF] Panet | An Analysis of Gravitational Gradients in Rotated Frames and Their Relation to Oriented Mass Sources[END_REF], predominant in the total gravity signal. This is a major challenge in the application of satellite gravity data to track deeper deformations, calling for dedicated analysis techniques. To solve this separation challenge, we will analyze horizontal gravity gradients rather than the geoid. Indeed, gravity gradients help identify a source from the spatial shape of its gravity signal, which is nely described thanks to the double dier-entiation of the gravity potential [START_REF] Panet | Migrating pattern of deformation prior to the Tohoku-Oki earthquake revealed by GRACE data[END_REF]. This way we can unravel smaller signals if their geometry diers from that of the predominant ones.

Here, we investigate whether anomalous gravity variations preceding the 2010 Maule event can be detected in the GRACE data. We consider a broad space-time window around the earthquake, from January 2003 to July 2014 in a 90 • × 120 • wide region around Central Chile. We rst analyze dierent sets of GRACE gravity eld models to search for abnormal signals before the rupture. For that, we enhance small gravity variations using gravitational gradients reconstructed from GRACE at dierent spatial scales. Then, we evaluate the obtained signals with respect to independent estimations of water storage changes by hydrological models and in-situ observations. This analysis allows us to propose and discuss a deeper origin inside the solid Earth, involving slab deformations near 150-km depth prior to the rupture.

Data and methods

GRACE geoid models

To search for gravity signals associated to with the Maule earthquake in the South American subduction system, we apply a space-time analysis of the time series of the GRACE geoid models over the January 2003-September 2014 period. To assess the sensitivity of the signals to the North-South oriented striping noise that degrades the GRACE geoids, we considered three sets of geoid models obtained by dierent groups, provided in the form of spherical harmonics expansions: the CNES/GRGS Release 3v1 (GRGS) up to degree/order 80 [START_REF] Lemoine | Temporal gravity eld models inferred from GRACE data[END_REF], the ITSG-2016 solution up to degree/order 90 60 [START_REF] Mayer-Gürr | ITSG-Grace2016 -Monthly and Daily Gravity Field Solutions from GRACE[END_REF] and the CSR Release-06 solution (CSR) up to degree/order 90 60 [START_REF] Save | GRACE RL06 reprocessing and results from CSR[END_REF] the CSR Release-06 solution (CSR) up to degree/order 60 [START_REF] Save | GRACE RL06 reprocessing and results from CSR[END_REF] and the ITSG-2016 solution up to degree/order 60 [START_REF] Mayer-Gürr | ITSG-Grace2016 -Monthly and Daily Gravity Field Solutions from GRACE[END_REF]. For the studied area at the beginning of 2010, we indeed found less striping in the North-South gravity gradients in the ITSG-2016 solution as compared to the more recent ITSG-2018 release.

Due to a dierent analysis of the GRACE observations, the ITSG-2016 and CSR gravity models show a higher level of striping artefacts than the GRGS solution. To minimize the striping these artefacts in these last two elds, we truncated their spherical harmonics expansion at the degree and order 40. We have veried that at this 500-km resolution, the signal-to-noise ratio remains favourable.

Hydrological models and in-situ data

Separation of solid Earth and hydrological signals is based on both a model-driven and a data-driven approaches. We considered an ensemble of four complementary hydrological models: 1. GLDAS NOAH 2.1 land surface model [START_REF] Rodell | The Global Land Data Assimilation System[END_REF], 2. WGHM global hydrological model [START_REF] Hannes | The global water resources and use model watergap v2. 2d: model description and evaluation[END_REF], 3. ERA5-Land land surface model [START_REF] Dee | The ERA-Interim reanalysis: conguration and performance of the data assimilation system[END_REF] and 4. the regional MGB model for South America [START_REF] Vinícius | Toward continental hydrologichydrodynamic modeling in South America[END_REF]. For comparison with the GRACE observations, we reconstructed each month the geoid and the gravity gradients predicted by these dierent models, considering the direct newtonian attraction of the water loads and applying a thin layer approximation.

In specic regions, we estimate water storage changes from in-situ observations: river discharge (Q) provided by the Global Runo Data Centre (GRDC), precipitation (P) from the Global Precipitation Climatology Center (GPCC) [START_REF] Schneider | Gpcc full data monthly product version 2020 at 0.25 • : monthly land-surface precipitation from rain-gauges built on gts-based and historical data[END_REF] and actual evapotranspiration (E) provided by the Max Planck Institute [START_REF] Jung | Recent decline in the global land evapotranspiration trend due to limited moisture supply[END_REF] (see Appendix E for a more detailed description of the used datasets).

To separate solid Earth and hydrological signals, we designed both a model-driven and a data-driven approach to dene the impact of water redistribution on gravity. We considered four complementary hydrological models: 1. The global GLDAS NOAH 2.1 model (including soil moisture, snow and water stored in the canopy) at 0.25 • resolution [START_REF] Rodell | The Global Land Data Assimilation System[END_REF], 2. the global WGHM model (including soil moisture, snow, groundwater and surface water) at 0.5 • resolution [START_REF] Hannes | The global water resources and use model watergap v2. 2d: model description and evaluation[END_REF], 3. the global ERA5-Land model at 9-km resolution (including soil moisture and snow) [START_REF] Dee | The ERA-Interim reanalysis: conguration and performance of the data assimilation system[END_REF] and 4. the regional MGB model for South America (including canopy, soil moisture, ground water and surface water) at 10-km resolution [START_REF] Vinícius | Toward continental hydrologichydrodynamic modeling in South America[END_REF]. We reconstructed each month the geoid and the gravity gradients predicted by these dierent models, considering the direct newtonian attraction of the water loads and applying a thin layer approximation. Here, the model ensemble is used to better quantify errors arising from forcing data, model structure, and model spatial resolution.

In specic regions, we complete the model analysis with water storage changes inferred from in-situ observations. We considered observations of river discharge (Q), precipitation (P) and actual evapotranspiration (E). The precipitation is based on the Global Precipitation Climatology Center (GPCC) Full Data Monthly Version 2020 dataset [START_REF] Schneider | Gpcc full data monthly product version 2020 at 0.25 • : monthly land-surface precipitation from rain-gauges built on gts-based and historical data[END_REF]. The GPCC provides gridded gauge-analysis products derived from quality controlled station data, at 0.25 • resolution [START_REF] Rudolf | The Global Precipitation Climatology Centre (GPCC)[END_REF]. The actual evapotranspiration is provided by the Max Planck Institute [START_REF] Jung | Recent decline in the global land evapotranspiration trend due to limited moisture supply[END_REF]. It is estimated from a data-driven approach, based on a global monitoring network, meteorological and remote-sensing observations, and a machine-learning algorithm.

Finally, we used river discharge data and basin outlines provided by the Global Runo Data Centre (GRDC). The hydrological analysis is performed over 2005-2012 when discharge data is available, in order to remove properly annual and semi-annual signals. Furthermore, in order to remove the potential impact of systematic bias in the uxes data (e.g. [START_REF] Long | Uncertainty in evapotranspiration from land surface modeling, remote sensing, and GRACE satellites[END_REF]), a linear trend is tted on water storage changes over the 2005-2012 period.

Gravity gradients at dierent spatial scales

In order To separate signals associated with from mass sources of dierent sizes, shapes or orientations in the GRACE geoids, we reconstruct each month from these geoid models the Earth's gravity gradients at dierent spatial scales, expressed in spherical frames:

1) the distinction between signals of dierent sizes is made by a wavelet analysis of the GRACE gravity potential [START_REF] Holschneider | From global to regional analysis of the magnetic eld on the sphere using wavelet frames[END_REF]; 2) the source geometry is emphasized by computing horizontal gravity gradients. In cartesian coordinates, these gravity gradients result from a double dierentiation of the wavelet-ltered gravity potential with respect to the three directions of space [START_REF] Panet | An Analysis of Gravitational Gradients in Rotated Frames and Their Relation to Oriented Mass Sources[END_REF]; then they are expressed at each point in the local spherical frame through appropriate coordinate transformations. As they highlight gravity signals elongated orthogonal to the dierentiation direction, the obtained horizontal gradients provide us with a detailed description of the geometry of the gravity eld variations at each spatial scale, reecting the structure and spatial extent of the sources (Appendix Fig. S1). From a general point of view, rotating the spherical frame along the radial axis is well-suited to separate gravity variations along the orientation of a subduction zone which could be potentially related to an earthquake, from water mass redistribution signals following other orientations. As both GRACE noise and the South-American subduction zone follow a North-South orientation, we average the gradients over a range of orientation (φφ gravity gradients from -10 to 10 • clockwise spherical frame rotations) around the North-South direction to increase signal-to-noise ratio, and we used directions close to East-West, orthogonal to the GRACE noise. The same methodology is applied to hydrological models.

Piece-wise linear t of anomalous signals

We then analyze the time series of gravity gradients at the dierent scales and in the different orientations, in order to search for anomalously large gravity variations before the Maule earthquake and compare them with the consecutive co-seismic signals. We rst estimate and remove from the time series annual and semi-annual sinusoidal terms accounting for the seasonal variability and a long-term trend, all estimated over the 2003/01-2008/12 period to limit the potential impact of a precursor, and apply the correction to the whole time series. Then, we analyze the residual time series g(t) as follows:

• We perform a piece-wise linear t [START_REF] Panet | Migrating pattern of deformation prior to the Tohoku-Oki earthquake revealed by GRACE data[END_REF] of the time series g(t) with a free jump in March 2010, which is the rst month actually recording the co-seismic step in the GRACE data (Appendix Fig. S2). For that, we decompose the time series of gravity gradients into four March 2011. The rst interval is the reference before the earthquake. The second interval represents the variations in the months preceding the earthquake, potentially including a fast pre-seismic signal. It is xed to 8 months as a compromise between a too short interval (for which the trend estimate would be very sensitive to noise) and a too long one (for which the meaning of a sudden pre-seismic gravity variation would be lost). Actually, our conclusions do not change when we vary the length of this interval between 1 and 8 months.

The free jump in March 2010 highlights the co-seismic signal. For the two last intervals, we take t 2 equal to March 2011. This way, we account for variations of the post-seismic gravity signals between the rst year and the rest of the time series.

• Anomalous variations before the earthquake are identied by the combination of 1. a large trend (>0.1 mEötvös) over the 8 months interval preceding the rupture in the piecewise linear t, together with 2. an abnormal gravity gradient signal in the month before the earthquake (February 2010), marked by a very low probability of occurrence in the observations (≥ 5σ, i.e. probability below 2.5 10 -6 %). To detect these abnormal variations in the monthly gravity gradients, we assume that the time series of residuals g(t) follow a Gaussian distribution and calculate its parameters, at each spatial grid point, over the reference period 2003/01-2008/12. This way we detect anomalous signals before the rupture, both in the February 2010 monthly snapshot and at a timescale of a few months, without making any hypothesis on the behaviour of the rest of the time series after February 2010 (such as the occurrence of a co-seismic variation). This approach allows us to detect a large and monotonous variation progressively realized over a few months, which culminates in highly abnormal values at the end of the considered period.

• We search for the co-seismic signal, identied by the combination of 1. a large jump in March 2010 in the piece-wise linear t, together with 2. an abnormal gravity gradient signal in March 2010 (probability below 2.5 10 -6 %), and 3. a shift in the distribution of the residual time series g(t) after the earthquake as compared to the years before, indicating some degree of persistence over time of the co-seismic jump (step-like temporal variation). This last criterium is simply implemented by a threshold in the amplitude of a Heaviside function centered at the time of the earthquake. This method enables us to make a clear distinction between fast variations right before the rupture and the co-seismic variation itself. The co-seismic amplitude is given by the jump in the piece-wise linear t, while the linear trends after the rupture can approximate faster and slower post-seismic signals. Over these relatively short post-seismic time intervals, exponential and logarithmic behaviors, as may be expected in the presence of afterslip or visco-elastic deformations, can be approximated by linear evolution [START_REF] Ivins | Transient creep of a composite lower crust: 1. Constitutive theory[END_REF]. Fitting the postseismic behavior is required to better estimate the pre-and co-seismic signals but remains beyond the scope of the paper.

Results

Slow to fast gravity jumps near the epicenter

We present here the results obtained for an analysis scale of 800-km, commensurate with the rupture length. At larger spatial scales, the earthquake signals progressively decay as the scales become too large as compared to the spatial extent of the signal. The smallest spatial scale that can be reached (500-km) given the resolution of the GRACE geoids is presented in the Section 4.2.1.

Around the epicenter, the two-lobe gravity signal conrms the impact of the co-seismic jump during the month of March 2010 (Fig. 1b, zoomed in Fig. 1d) that was observed in previous studies. It comprises a negative gravity gradient anomaly over the ocean and a positive one over the continent with amplitudes up to 0.19 mEötvös, persistent in the consecutive years. For comparison, the co-seismic gravity gradient signature of the 2011 Tohoku earthquake is twice larger at this spatial scale. Furthermore, we nd that this co-seismic variation is preceeded by an anomalously large gravity gradient increase in the months before the rupture, located North of the epicentral area (Fig. 1a, zoomed in Fig. 1c)(Fig. 1a,c). It stands out as the most abnormal signal over the entire South American continent, manifested by the most widespread February 2010 anomaly. This signal exceeds the 5σ level of the long-term distribution, after a large increase over 8 months. Thus, among all the medium-scale gravity variations in the months before the earthquake, the largest one is the closest to the epicenter. It is also detected in the wavelet-ltered geoid, although less well separated from neighbouring hydrological sources (Appendix B). These variations are reected in time series of the gravity gradients at dierent locations within these signals (Fig. 1e). For points near the maximum of the co-seismic anomaly below the latitude -34 • N, a large and sudden jump is observed between February and March 2010, starting from a high value in the time series in February. As we move towards the North (above latitude -34 • N), away from the epicenter, the co-seismic jump decreases while an anomalously large positive trend is observed in the preceding months, leading to a slow jump in the time-series. Because of this gradual increase, ending up with two highly anomalous values in January and February 2010, we infer a duration of the signal of at least 2 months. At the monthly resolution of the used GRACE data, and in the presence of hydrological contributions, it remains dicult to point out the exact starting time of this signal, which might be earlier in 2009, and whether it develops continuously or through pulses at submonthly timescales.

We have assessed the sensitivity of the February 2010 GRACE monthly geoid to the coseismic mass redistributions. The Maule earthquake took place at 6 : 30 UTC on February 27, 2010. When cCounting the number of February 2010 orbits in a 20 • vicinity of the epicenter, we nd less than 5% of the monthly orbits in the time interval after the rupture. This conrms that the February 2010 signal North of the epicenter is not signicantly impacted by the co-seismic signal, and that the latter which is rst recorded in the March 2010 geoid.

Investigation of other GRACE gravity solutions

To assess the robustness of the signals before the rupture, we tested whether they could also be found in two other sets of GRACE geoid models, the CSR06 and the ITSG-2016 solutions , in addition to the GRGS solution. For that, we extracted their common spacetime patterns of variability using a Singular Value Decomposition between pairs of models (expressed in terms of gravity gradients): 1. GRGS03 versus CSR06 and 2. GRGS03 versus ITSG-2016. We found a highly coupled behaviour of each pair of solutions, featuring a slow jump initiated a few months before the earthquake in the region of the pre-and co-seismic signals, in both the North-South and the East-West directions (Appendix Fig.

D.

2) (Appendix D.2, Appendix Fig. S5). As each individual solution, the average of these three datasets shows the same behaviour (Fig. 2a-d). It is illustrated by the following analysis allowing us to identify a slow jump near March 2010 completed to a large extent in February 2010. As the CSR and ITSG-2016 solutions show more a higher level of striping noise than the GRGS one, we use a more constrained time evolution model, with less degrees of freedom than done for GRGS in the previous section. For that, in each set of GRACE solution as well as and in their average, we estimated a Heaviside step function in March 2010. This estimate is not very sensitive to the exact timing of the step, due to the length of the time series: it will not distinguish between steps completed in February or in March. We dene anomalous signals in February 2010 such that their amplitude exceeds 50 % of the estimated step, and their probability of occurrence is low (assuming a Gaussian distribution of the monthly gravity gradient values, as done before). and for their average for two ranges of orientations (North-South and East-West) (in the North-South and East-West orientations). A good agreement is found between the three solutions: the North-South oriented gravity gradient exhibits an anomalous increase before the earthquake in the same area for each individual solution and their average (Fig. 2a-d). The above discussed temporal pattern of a slow jump in the gravity gradients is detected in both North-South and East-West directions (Fig. 2f-g). It is associated with a well-resolved spatial pattern in the same area for both orientations, pointing to anomalously large signals during that month (Fig. 2d-e). This behaviour appears unique over all South America (Fig. 2a-c). Finally, we have also We nally investigated anomalously large signals in February 2010 without any constraint on a step-like evolution of the time series. A few other signals are detected in the CSR and ITSG-2016 solutions, in addition to the gravity gradient increase before the Maule earthquake (Appendix Fig. S11). However, the singularity specicity of the February 2010 signal in the region of Maule is further evidenced from comparisons with hydrological models and in-situ data.

4 Sources of the gravity signal before the earthquake 4.1 Inaccuracies in the GRACE data processing

We have investigated whether the gravity signal before the Maule earthquake could result from inaccuracies in the data analysis: (1) striping artefacts, (2) errors in the atmospheric dealiasing model, or (3) over or under-correction of the seasonal cycle. First, we estimated empirically the level of striping errors in the monthly horizontal gravity gradients from 2003 to 2014, for the 500-1000km spatial scales. For that, we computed each month the rms of the gravity gradients over a wide oceanic area centered at the latitude of the Maule earthquake epicenter. We conclude that the February 2010 anomalous gravity gradient variations exceed the noise by a factor 4 to 8 depending on the spatial scale. Looking then at the atmospheric dealiasing model of the GRGS geoids, based on the ECMWF ERA-Interim reanalysis [START_REF] Dee | The ERA-Interim reanalysis: conguration and performance of the data assimilation system[END_REF], we nd that, in the region of Maule, the amplitude of the February 2010 non-seasonal atmospheric signal is hundred times smaller than the GRACEobserved anomaly before the rupture. This reects the fact that the modelled atmospheric contribution is almost purely seasonal in the studied area. We nally investigated whether the gravity signal before the earthquake could result from an over-or under-correction of the seasonal cycle in the GRACE data or in the atmospheric model. We rst remark that an error in a periodic correction should appear as a periodic residual in the time series. Such behaviour is absent from the gravity gradient time series shown in Fig. 1e until 2011 at least. Comparing the amplitude of the GRACE pre-seismic signal with that of the seasonal cycles tted in the GRACE data or in the modelled atmospheric contribution at the same location (north of the epicenter) north of the epicenter, we notice that the February 2010 anomalous GRACE gravity gradient signal is at least two times larger than the amplitude of the GRACE seasonal cycle and ten times larger than that of the atmospheric model there. Thus, neither the striping artefacts nor the atmospheric or seasonal corrections can explain the observed variations. Most of the gravity signal recorded by GRACE comes from continental water mass redistribution, so we investigated a possible hydrological source to the gravity gradient increase before the earthquake. Note that the GIA gravity gradient signal from the Patagonian Ice Field does not aect our results due to dierent timescale, and even more so as we have removed a long-term trend from the data. The horizontal gravity gradients increase corresponds to mass decrease, hence a drying signal. We rst compare the GRACE anomaly with the predictions of the hydrological models presented in the Section 2.2.

In addition to the scale 800-km, we present here the results at the 500-km scale. As shown in Fig. 3, this 500-km scale gives access to ner details including large hydrological signals that might be not resolved at the 800-km scale (Fig. 3c,e). This higher resolution is interesting because it provides a better spatial separation between sources in the Andean Cordillera and those in the watersheds of Argentina, especially the major La Plata basin (hydrological context map Appendix Fig. S7). This smaller spatial scale is closer to the characteristic scale of some drainage basins in the region than the scale 800-km ; at the same time it brings other challenges as we notice a higher number of anomalous signals in the GRACE gravity gradients. They include a large negative anomaly in the La Plata basin (labelled 3 in Fig. 3c), and a positive anomaly around the point (294.5 3c. This positive signal is located at the North-East of the 800-km scale GRACE anomaly before the earthquake, which integrates the 500-km scale anomalies 1 and 2 (Fig. 3e).

• E ; -31 • N), labelled 2 in Fig.
We rst compare the time evolution of the GRACE-observed signal with that predicted from the hydrological models at locations spanning the corresponding area. Fig. 3a shows that GRACE and the models are coherent in the North-Eastern lobe (anomaly 2). In

Comparison of dierent hydrological basins

In a last comparison, we estimated the GRACE anomalous signals in February 2010 in all the investigated gravity eld solutions (GRGS, CSR and ITSG-2016)the GRGS, CSR and ITSG gravity eld solutions without any hypothesis on a step-like evolution of the time series in March 2010. With a lower level of abnormality of the February 2010 signal than in Fig. 1, we detect anomalous gravity gradient variations also over the Orenoco, Chaco and La Plata basins (Appendix Fig. S7,S11). These anomalous gravity gradient variations are probably related to the 2009-2010 El Niño event, which resulted in large mass redistributions associated with droughts in the Amazon basin and oods in the La Plata basin. We notice that the GRACE signals agree well with the predictions of the hydrological models in all the drainage basins, except for the anomaly located north of the epicenter of the Maule earthquake (Appendix Fig. S11). There, El Niño brings increased precipitations in winter (from June to December), and has no direct impact on the summer rainfalls [START_REF] Cai | Climate impacts of the El NiñoSouthern Oscillation on South America[END_REF]. Thus, it seems hardly consistent with the GRACE mass variations in this area, contrary to the other basins.

In situ observations

As a complement to the hydrological models, we investigated in-situ observations of horizontal and vertical water uxes (river discharges, precipitations and evapo-transpiration).

The GRACE signal is located in the region of Mendoza in Argentina, considered arid with its 150 to 300 mm of annual rainfall (Fig. 4a), and part of the so-called Arid Diagonal of South America. Most of the water in the watersheds comes from the annual melting of snow and ice and is transported from the mountain by the rivers ; it is collected through articial dams by the regional population, organized in three oases. Therefore, the rivers discharge upstream of the dams is representative of the water inux in the zone and appears signicantly correlated with the regional snow accumulation at inter-annual timescales [START_REF] Masiokas | Snowpack Variations in the Central Andes of Argentina and Chile, 19512005: Large-Scale Atmospheric Inuences and Implications for Water Resources in the Region[END_REF].

Fig. 4b-f shows the variations of discharge at ve stations upstream of the dams in the

Central Andes over the 2000-2017 period. First, we notice that the annual discharge shows a large decrease starting at the end of 2010. This decrease coincides with the beginning of a mega-drought in Central Chile, manifested as a sequence of dry years with reduced annual precipitations [START_REF] René | The 20102015 megadrought in central Chile: impacts on regional hydroclimate and vegetation[END_REF]. However, Fig. 4 also shows that the eect of this mega-drought on the river ows is still limited in February 2010; it starts later that year, together with a strong 2010/2011 La Nina episode initiated in June [START_REF] René | The 20102015 megadrought in central Chile: impacts on regional hydroclimate and vegetation[END_REF]. Thus, this mega-drought cannot explain the GRACE mass decrease signal observed ve months before, in February 2010.

Second, we compared the mass anomaly explaining the GRACE signal before the earthquake with the water storage variation corresponding to these discharge data. Here we neglect the impact of the vertical water uxes, which is consistent with an inux of water mostly from snow and ice melting in this region. The 0.18 mEötvös (resp. 0.19 mEötvös) amplitude of the GRACE signal at the 800-km scale (resp. 500-km scale) can be modelled by a water mass source of width 4500-km, length 500-km and thickness 300 mm EWH (Equivalent Water Height). It corresponds to ∼ 60 km 3 of water storage decrease between the beginning of January 2010 and the end of February 2010, or a discharge of 30 km 3 .month -1 during two months. Distributing this mass transport over the four major rivers of the region (Mendoza, Desaguadero, Tunuyan and San Juan) still leads to a 7.5 km 3 .month -1 ow. This is fteen times more than the largest monthly ow recorded in the region for the period 2000-2017 (0.5 km 3 .month -1 for the Mendoza river in 2006, Fig. 4c), and still much larger than the maximum historical discharge of 1 km 3 .month -1 , recorded in 1987 [START_REF] Moreiras | Climatic eect of ENSO associated with landslide occurrence in the Central Andes, Mendoza Province, Argentina[END_REF].

To complete our analysis of the hydrological sources in this region, we carried out a more precise estimation of the water storage variations of the four watersheds intersecting the GRACE anomaly, for comparison with the GRACE-observed mass transport in February 2010. For that, we estimate the water storage change (S-S0) at basin scale from in-situ observations of uxes, namely precipitation grids (P), actual evapotranspiration (E) grids and river specic discharge (Q). These variables are related via the mass balance equation dS dt = P -E -Q, which is integrated to derive storage variations per unit surface S -S 0 = ∆Sdt, in mm EWH in each basin, which is transformed in volume in km 3 by multiplying by the basin area. In the above equation, the river discharge (Q) denes the basin response to the eective rainfall (P-E). Here, the river discharge stations are located downstream of each basin on which we apply the mass balance equation. As shown in Fig. 5, the obtained water storage variations amount to ∼ 1.2 km 3 over two months between January and February 2010. This value is far smaller than the dozen of km 3 of water needed to explain the GRACE signal. Thus, the observed gravity gradient anomaly is not likely to be explained by a water source in this regional context.

Implications for deep Earth pre-seismic processes

The above analysis supports a solid Earth origin of this gravity gradient signal, involving mass decrease at depth. Here, we discuss its possible origin in the context of the South American subduction.

Multiple lines of evidence indicate that a regional change in stress state inside the crust along the South Chile subduction zone occurred prior to the megathrust earthquake, in relation with deeper slab motions. Bouchon et al [START_REF] Bouchon | Potential slab deformation and plunge prior to the Tohoku, Iquique and Maule earthquakes[END_REF] report a pre-earthquake seismic activity which began in early January 2010. This activity was characterized by an initial burst of seismic activity at depth followed by shallow foreshocks. In the USGS catalog, a large (Mw 5.8) intermediate-depth earthquake (depth ∼150-km) occurred on 12 February 2010 in the region where the pre-seismic gravimetric signal is reported, suggesting that the change in stress state aected a very wide area around the epicenter region. Large-scale anomalous GNSS displacements were detected four months before the main event over the whole South Chile subduction zone, corroborating such a hypothesis [START_REF] Bedford | Months-long thousand-kilometrescale wobbling before great subduction earthquakes[END_REF]. The spatial and temporal correlation between the shallow and the deep seismicity activities, the extensional mechanism of the deep shocks and the trenchward motion of 4-8 millimeters lead the authors of these previous studies to the same conclusion. Bouchon et al. [START_REF] Bouchon | Potential slab deformation and plunge prior to the Tohoku, Iquique and Maule earthquakes[END_REF] and Bedford et al. [START_REF] Bedford | Months-long thousand-kilometrescale wobbling before great subduction earthquakes[END_REF] propose that these observations are related to a sudden increase of slab pull at depth interacting with shallow slow slip further updip.

All these observations support an extensional deformation of the slab along its subduction direction. Also observed in the case of the 2011 Tohoku-oki earthquake [START_REF] Panet | Migrating pattern of deformation prior to the Tohoku-Oki earthquake revealed by GRACE data[END_REF], such deformation is consistent with our GRACE pre-seismic mass decrease signal. To evaluate its magnitude, we modeled a pre-seismic intra-slab extension using a model of quasi-static normal faulting in a vertically stratied elastic medium. We found that the observations can be explained by 1.5 m of slip along a 45 • dipping, North-South striking plane of width 100-km and length 500-km, located between the depths 115 and 185-km (Fig. 6). This corresponds to a Mw 8.2 event over a few months. As in the case of the Tohoku-oki earthquake, where a precursor signal equivalent to a Mw 8.4 rupture was detected [START_REF] Panet | Migrating pattern of deformation prior to the Tohoku-Oki earthquake revealed by GRACE data[END_REF], the magnitude of the pre-seismic event is smaller than that of the consecutive co-seismic rupture and its spatial extent is commensurate with the length of the co-seismically ruptured area. This corresponds to a smaller amount of deep deformation distributed over a wider zone around 250-km depth before the Tohoku-oki 2011 earthquake, and a larger but more localized deformation around 150-km depth before the 2010 Maule earthquake.

Interestingly, for both Maule and Tohoku, this pre-seismic gravity signal occurs in a region where the geometry of the subduction changes. The Tohoku pre-seismic signal coincides with a change in the strike of the subduction, in the vicinity of a triple junction. In the case of Maule, the gravity pattern before the rupture is located ∼ 400-km North-East of the epicenter, in a transition section where the dip of the subducted Nazca plate changes sharply. North of -33 • N, the Nazca plate remains sub-horizontal at about 100-km depth for several hundred kilometers, before dropping steeply into the mantle. South of -33 • N, the slab plunges continuously at an angle of 25 • [6] [38]. This change of dip is associated with a slab hole at ∼ 200 -300km depth, at the latitudes of the February 2010 GRACE signal [START_REF] Lynner | Mantle ow through a tear in the Nazca slab inferred from shear wave splitting[END_REF] [START_REF] Daniel | The nature of subslab slow velocity anomalies beneath South America[END_REF]; mantle ow through this opening at depth in the at slab is suggested from seismic tomography. The rapid change of slab geometry causes an increase of the slab stresses, due to localized lateral deformations [START_REF] Anderson | Anisotropy and mantle ow in the Chile-Argentina subduction zone from shear wave splitting analysis: CHILE-ARGENTINA SUBDUCTION ZONE[END_REF]. In addition, mantle ow beneath the at slab and through the slab hole [START_REF] Lynner | Mantle ow through a tear in the Nazca slab inferred from shear wave splitting[END_REF] [13] may contribute to entrain the motion of the subducted plate towards the depth.

The spatial and temporal synchronization observed between the pre-seismic gravity signal with the megathrust earthquake suggests that these two events are linked. We thus propose that, between the end of 2009 and March 2010, the deep pre-seismic slab deformation migrated both upwards to the surface and laterally along the slab. From the pre-seismic to the co-seismic phase, the gravity gradient signal indeed shifts from the North/East (around the point 292 • E; -32.5 • N, corresponding to a 150-km depth for the top of the slab) to the South/West (around the point 290 • E; -36 • N). This lateral migration towards the surface is reminiscent of the variations of coupling of the subduction interface, at the shallower depths. The pre-seismic extension is indeed located downdip a transition area between a low coupling segment in the North and a high coupling one in the South [START_REF] Metois | Interseismic coupling, megathrust earthquakes and seismic swarms along the chilean subduction zone (3818 s)[END_REF], and the motion propagates towards the more coupled zone in the South, where the rupture occurred [START_REF] Metois | Interseismic coupling, megathrust earthquakes and seismic swarms along the chilean subduction zone (3818 s)[END_REF] in the historical seismic gap of central Chile [START_REF] Udias | The Large Chilean Historical Earthquakes of 1647, 1657, 1730, and 1751 from Contemporary Documents[END_REF]. Comparing our co-seismic gravity gradient signal with that predicted from a geodesy-based co-seismic slip distribution model [START_REF] Lorito | Limited overlap between the seismic gap and coseismic slip of the great 2010 Chile earthquake[END_REF], we found a reasonable agreement between the GRACE-derived co-seismic anomaly and the modelled one (Appendix Fig. S12).

Finally, the existence of a gravity gradient signal without large surface displacements had been noticed for both the pre-and co-seismic signals of the Tohoku-Oki earthquake, and is again observed in the case of the Maule earthquake. In the Maule pre-seismic phase, the slab deformation estimated above should have generated centimetric surface displacements, which is one order of magnitude greater than the displacements highlighted by Bedford et al. [START_REF] Bedford | Months-long thousand-kilometrescale wobbling before great subduction earthquakes[END_REF]. Even if the spatial coverage of the GNSS network is very sparse in this region far of the trench, it is unlikely that such ground surface displacement remains undetected. Considering other processes that may have contributed to the observed preseismic gravity variation, we have investigated the eect of transient uid release related to oceanic plate dehydration. Upon subduction, hydrous minerals release water from the slab over a range of depths [START_REF] Dragovic | Pulsed dehydration and garnet growth during subduction revealed by zoned garnet geochronology and thermodynamic modeling, Sifnos, Greece[END_REF] [37] [45] [20]; the mechanisms, dimensions and timescales of uid ow within the subduction zone are not well understood. Dehydration uid is generally considered as a continuous process occurring on a time scale of 10 5 -10 6 years [START_REF] Dragovic | Pulsed dehydration and garnet growth during subduction revealed by zoned garnet geochronology and thermodynamic modeling, Sifnos, Greece[END_REF] [20]; however, it has been suggested that channelized uid ows are highly localized, accumulating and releasing high uid volumes within short time interval (1-4 months) [START_REF] John | Volcanic arcs fed by rapid pulsed uid ow through subducting slabs[END_REF] [37] [START_REF] Taetz | Fast intraslab uid-ow events linked to pulses of high pore uid pressure at the subducted plate interface[END_REF]. The increasing slab pull force accommodated by extensional cracks and motion along normal faults may have promoted the creation of uid pathways, improving the drainage of the subducting plate. Even if the mechanisms of large-scale uid release from the slab at depth remain unclear, a simple model based on a slight change of the porosity ( 10-5 -10 -6 of relative volume variation) of the slab segment in extension at depth is able to fully or partially explain the the observed gravity gradient signal. We have modelled the gravimetric signal caused by a uid inltration in the cracks generated by a deep extension, i.e. a density variation of ρ fluids -ρ rocks = -2400kg.m 3 over a volume of width 200-km, length 500-km and thickness 80-km, with a porosity δV V = 5 10 -6 , located around 150-km depth. According to this model, the uid-related mass transfer at depth induced a widespread gravity variation (Fig. 6e-f ), able to explain the observed gravity gradient signal (Fig. 6a-b). We thus hypothesize that the absence of large surface displacements could be related to a signicant contribution to the observed gravity signal of such deeper mass redistributions associated with uid migration, accompanying the extensional deformation of the slab.

Conclusion

From a dedicated analysis of time series of GRACE-derived gravity gradients, we have detected an anomalously large gravity gradient increase in the months before the Mw 8.8 2010 Maule earthquake, in the north of the epicentral area, most likely caused by mass redistributions at depth within the solid Earth on the continental side of the subduction.

This gravity signal prior to the rupture can be explained by a deep extensional deformation of the slab along the subduction direction, equivalent to a Mw 8.2 normal faulting event.

This event is commensurate with the precursor signal detected before the Tohoku-Oki earthquake [START_REF] Panet | Migrating pattern of deformation prior to the Tohoku-Oki earthquake revealed by GRACE data[END_REF], and also located in a region of changes in the geometry of the subducted slab. In the case of 2010 Maule earthquake, the pre-seismic mass decrease signal highlights a larger amount of mass anomaly distributed over a more localized zone with a length commensurate to that of the co-seismically ruptured area. We notice that the modelled Maule pre-seismic normal faulting event should have generated centimetric surface displacements, which have not been observed by GNSS [START_REF] Bedford | Months-long thousand-kilometrescale wobbling before great subduction earthquakes[END_REF]. This leads us to propose that part of this gravity signal could reect deep mass redistributions from large-scale uid release promoted by extensional cracks and normal faults in the subducted slab -even if it remains dicult to decipher without ambiguity the physical processes at the origin of the gravity variations observed prior these two earthquakes. Nevertheless, the existence of these interactions be- Fig. S3 compares the February 2010 signals in the studied region, as obtained from the total geoid, from the wavelet-ltered geoid and from the horizontal φφ gravity gradients.

The contribution of the major hydrological sources from the La Plata basin predominates in the total geoid (Fig. S3a), partially masking the Maule pre-seismic signal. This smaller signal is emphasized in a high-resolution wavelet ltering of the geoid (Fig. S3b) and can be detected in the corresponding time series (Fig. S3e). However, it is not well separated from the nearby La Plata anomaly, which still perturbs the amplitude of the 500-km scale geoid signal in the region of the pre-seismic anomaly. The horizontal gravity gradients is reected by the fraction of the squared covariance explained by this mode (SCF k ). For the k th common mode, it is given by: SCF k = L(k,k) trace(L) .

D.2.2 Results

In the North-South direction, the rst common modes between GRGS and CSR, and between GRGS and ITSG-2016, represent more than 75 % of the squared covariance between each pair of solutions, in both cases (Fig. S6c,k). Thus, they point to a highly coupled behaviour of each pair of solutions. The associated spatial pattern covers the locations of the co-seismic and pre-seismic signals (Fig. S6a,b,i,j); the temporal pattern shows a progressive increase in the gravity gradients over three months up to March 2010, in each gravity gradient solution (Fig. S6d,l). In the East-West direction, the rst common modes explain a smaller amount of variance, at the level of 50 % of covariance (Fig. S6g,o

E Hydrological models and in-situ data

To separate solid Earth and hydrological signals, we designed both a model-driven and a data-driven approach to dene the impact of water redistribution on gravity. We considered four complementary hydrological models: 1. The global GLDAS NOAH 2.1 model (including soil moisture, snow and water stored in the canopy) at 0.25 • resolution [START_REF] Rodell | The Global Land Data Assimilation System[END_REF], 2.

the global WGHM model (including soil moisture, snow, groundwater and surface water) at 0.5 • resolution [START_REF] Hannes | The global water resources and use model watergap v2. 2d: model description and evaluation[END_REF], 3. the global ERA5-Land model at 9-km resolution (including soil moisture and snow) [START_REF] Dee | The ERA-Interim reanalysis: conguration and performance of the data assimilation system[END_REF] and 4. the regional MGB model for South America (including canopy, soil moisture, ground water and surface water) at 10-km resolution [START_REF] Vinícius | Toward continental hydrologichydrodynamic modeling in South America[END_REF]. We reconstructed each month the geoid and the gravity gradients predicted by these dierent models, considering the direct newtonian attraction of the water loads and applying a thin layer approximation. Here, the model ensemble is used to better quantify errors arising from forcing data, model structure, and model spatial resolution. Then, we use the gravity gradient signals predicted from these models for comparisons with the observed ones, in order to discuss the origin of the GRACE anomalies. Comparing these four dierent models also allow us to identify robust features and model-dependent errors arising from forcing data, model structure, and model spatial resolution.

In specic regions, we complete the model analysis with water storage changes inferred from in-situ observations. We considered observations of river discharge (Q), precipitation (P) and actual evapotranspiration (E). The precipitation is based on the Global Precipitation Climatology Center (GPCC) Full Data Monthly Version 2020 dataset [START_REF] Schneider | Gpcc full data monthly product version 2020 at 0.25 • : monthly land-surface precipitation from rain-gauges built on gts-based and historical data[END_REF]. The GPCC provides gridded gauge-analysis products derived from quality controlled station data, at 0.25 • resolution [START_REF] Rudolf | The Global Precipitation Climatology Centre (GPCC)[END_REF]. The actual evapotranspiration is provided by the Max Planck Institute [START_REF] Jung | Recent decline in the global land evapotranspiration trend due to limited moisture supply[END_REF]. It is estimated from a data-driven approach, based on a global monitoring network, meteorological and remote-sensing observations, and a machine-learning algorithm. Finally, we used river discharge data and basin outlines provided by the Global Runo Data Centre (GRDC). The hydrological analysis is performed over 2005-2012 when discharge data is available, in order to remove properly an annual and semi-annual signals.

Furthermore, in order to remove the potential impact of systematic bias in the uxes data (e.g. [START_REF] Long | Uncertainty in evapotranspiration from land surface modeling, remote sensing, and GRACE satellites[END_REF]), a linear trend is tted on water storage changes over the 2005-2012 period.

3. 2

 2 Singularity of the gravity gradient signal before the rupture We have investigated the unique character in space and time of the gravity gradient increase before the earthquake. The spatial unicity over the South American continent is illustrated from Fig.ure 1a, where the signal near Maule appears as the largest one. The unicity of this signal in time can be directly observed from the time series in Fig.ure 1e. These time series show that, in the considered region near the Maule 2010 epicentral area, a signal of a comparable amplitude had not been recorded before, nor in the consecutive years (note that the annual cycle correction may degrade in the end of the time series). This is conrmed when repeating the same time series analysis as described above, for hypothetic earthquake times t e spanning the [March 2004 -March 2010] interval with a monthly time step. Appendix Fig.ure S4 shows the obtained anomalous signals cumulated over the [t e -9 months to t e -1 month] intervals: there is no equivalent to the July 2009 -February 2010 gravity gradient signal over the whole period and the whole continent.

Fig.ure 2

 2 Fig.ure 2 presents the results obtained for each of the three gravity gradient solutions used in this study for one direction (North-South) (in the North-South direction)

4. 2

 2 Hydrological signals from global and regional models 4.2.1 Predicted signals in the vicinity of the epicentral area

  the southern lobe closer to the epicenter (anomaly 1, especially between latitudes -33 • N and -34.5 • N), GRACE and the models agree from 2004 to the end of 2009 but this consistency degrades starting from January 2010, after ve years of low variability. At the beginning of 2010, the GRACE time series indeed exhibit a fast and large increase whereas the hydrological signals remain in the continuity of the previous years. Thus, the North-Eastern component of the GRACE signal (anomaly 2) is probably impacted by a hydrological contribution near the western end of the Plata basin, a temperate climate area (Fig.3d, Appendix Fig.S8), while the South-Eastern component of the signal (anomaly 1) remains anomalous with respect to the hydrological models. Interestingly, it is located in an arid zone, as discussed later in this work.4.2.2 Spatial patterns of the modelled hydrological signalsThese dierences are conrmed when investigating the localization of the hydrological signals. Because the GRACE signal is probably not related to a variation in the seasonal cycle, we investigate the spatial patterns of the residual non-seasonal variability. During the 6 months period before March 2010 (September 2009 -February 2010) as well as in the 7 years before, the GLDAS, WGHM and MGB models do not predict any non-seasonal 500-km scale signal in the region of the GRACE anomaly 1, in both North-South and East-West directions (Appendix Fig. S8, S9). The hydrological signals are indeed controlled by the topographic reliefs of the Andes and the Chilean Coast Range, which localizes the rainfalls in a thin North-South elongated band on the Western ank of the mountains in Southern Chile, mostly south of the considered area. In contrast, the GRACE anomaly 2 is likely aected by a non-seasonal hydrological contribution from the western part of the La Plata basin. When comparing the hydrological models, we found a disagreement of ERA5-Land with the other models. This is mostly due to dierences in the seasonal cycle predicted by this model, leading to residual annual signals in the region of Maule (Appendix Fig. S10), unobserved by GRACE.

  tween slow mass variations at depth detected by GRACE and interplate seismicity, opens a new eld of research to better characterize and understand the dynamics of the seismic cycle at megathrusts. Observing again such interactions for large earthquakes in the future could lead to a paradigm shift in the study of the seismic cycle, which is today essentially based on the distribution of the recurrence times of large earthquakes for the estimation of the seismic hazard.

Figure 1 :

 1 Figure 1: Pre-seismic and co-seismic gravity gradient signals for the Maule earthquake. All panels show the 800-km scale, φφ GRGS gravity gradients, stacked for -10 to 10 • clockwise spherical frame rotations. Panels a and c: cumulated variation over the [July 2009 February 2010] interval with an absolute amplitude above 0.07 mEötvös, also shown in black contours in the map d. Panels b and d: March 2010 co-seismic variations. Red star: Maule earthquake epicenter; orange lines: plate boundaries [3]; violet lines: Pacic slab isodepth contours every 100-km [16]; brown lines: 3000 m topographic contours. Tectonic plates: NA = Nazca, SA = South-America, Red arrow = sudbuction direction. Panel e: time series of the gravity gradients after removing the annual and semi-annual cycles and a trend, at locations spanning the pre-and co-seismic anomalies, shown as black dots in the maps c and d. Blue dot in the time series: January 2010 ; green dot: February 2010 ; pink dot: March 2010.

Figure 2 :

 2 Figure 2: Anomalous gravity gradient signals before the Maule earthquake from the GRGS, CSR and ITSG-2016 gravity solutions and the average of the three solutions. Panels a-d and f : 800-km scale φφ gravity gradients in the local spherical frame, emphasizing North-South oriented signals. Panels e and g: 800-km scale θθ gravity gradients in the local spherical frame, emphasizing East-West oriented signals. Top panels: maps for each individual solution. Bottom line panels: maps and time series for the average of the three solutions. Panels a-e: Maps of anomalous gravity gradient signals in February 2010 (with a probability below 0.25% for ITSG-2016 , below 1% for CSR, below 2.5 10 -5 % for GRGS and below 0.01% for the average), persistent in time after March 2010 (see text). Panels f, g: time series of the gravity gradients at locations across the GRACE positive anomaly in February 2010, indicated by the black dots on the maps d, e respectively.

Figure 3 :

 3 Figure 3: Comparison of the GRACE GRGS gravity gradient signals with those predicted from four hydrological models in the region of Maule at two spatial scales. Panels a,c and d: scale 500-km, φφ gradients in the local spherical frame. Panels b and e: scale 800-km, same gradients. Panels a (resp. b): time series at points across the GRACE GRGS positive anomaly before the earthquake, marked in black dots in the map c (resp. e). GRGS time series in black ; time series from the hydrological models in colors as referred in the zoomed panel a. The gravity gradients from the ERA5-Land hydrology model have been scaled by a factor 0.5 for consistency with the other hydrological models (Appendix S10). Panels c: Map of the GRACE GRGS anomalous signals before the earthquake. This map is derived from the same analysis as in Fig. 1a, considering a lower level of abnormality of the February 2010 signals (outside the 0.02 -99.98 percentile range of the long-term distribution). Panel e: same as Fig. 1c. Panel d: Spatial patterns of the nonseasonal signals in the gravity gradients from the WGHM hydrological model, as expressed by the RMS of the 2009/09 2010/02 time series after subtraction of the annual, semi-annual and long-term trend components. For panels d and e, black lines are the contours of the positive GRACE GRGS anomaly shown in panel c.

Figure 4 :

 4 Figure 4: Panel a: map of the dierent climate zones of South America as reected by the annual amount of precipitations [23], from arid/semi-arid to humid, superimposed with the contours of the GRACE pre-seismic signal of Fig. 1a (black lines). Panels b-f: monthly river discharge at hydrometric stations in the vicinity of the GRACE pre-seismic signal and upstream of the dams (black dots in the zoomed map of panel a), reecting the annual supply of water in this arid region by snow and ice melting in the Andean Cordillera. Red dashed line: February 2010.

Figure 5 :

 5 Figure 5: Estimated water storage variations of the Desaguadero, San Juan, Mendoza and Tunuyan basins in the region of Maule. Panel a: mass changes over the basins from the beginning of January to the end of February 2010, superimposed with the contour of the GRACE pre-seismic signal of Fig. 1a (black line). Panels b-e: time series of the water storage variations for each basin, before and after removing a seasonal cycle.

Figure 6 :

 6 Figure 6: Comparison between the GRGS GRACE-observed (panels a-b) and the modeled (panels c-f) pre-seismic signals for the scale 800-km. Top panels: φφ gravity gradients, stacked for -10 to 10 • clockwise spherical frame rotations ; Bottom panels: θθ gravity gradients for a 10 • clockwise spherical frame rotation. Panel a: same as Fig.1c. The black lines depict the 0.10 mEötvös contours of the February 2010 gravity gradients. Panel b: same as Fig.1c for the considered θθ gravity gradients. Panels c-d: gravity gradient signals predicted from a model of quasi-static normal faulting in a vertically stratied elastic medium (see text). The violet lines are the slab isodepth contours every 100km [16]. Panels e-f: gravity gradient signals predicted from a model of density variations due to uid inltrations over a volume of 500 × 500 × 80-km at 150-km depth with porosity δV V = 5 10 -6 (see text).

Figure S2 :

 S2 Figure S2: Piece-wise linear t of the time series of the GRGS gravity gradients (800-km scale, φφ gravity gradients in the local south-east-up spherical frame, averaged for -10 • to 10 • clockwise frame rotations around the radial axis). Panel a-b: maps of the monthly gravity gradients in February and March 2010. Panel c: time series at the location marked by a black dot in the Panels a-b (thin grey line) and their piece-wise linear model (thick black line). Panel e d: cumulated variation over the [July 2009 -February 2010] interval, associated with abnormally large signals in February 2010 (probability below 2.5 10 -4 %) with an amplitude threshold of 0.1 mEötvös on the eight months trend. Panel d e: March 2010 co-seismic signal, associated with abnormally large signals in March 2010 (probability below 2.5 10 -4 %). Here, the level of abnormality of the February and March 2010 signals is lower than in Main Fig. 1.

  perform a better separation of these two signals, as shown by the comparison of the panel b with the panels c and d, where the amplitude of the La Plata anomaly has considerably decreased and that of the Maule pre-seismic signal starts to predominate. Indeed, the La Plata signal has a strong East-West component, which is ltered out in the North-South oriented gravity gradients.

Figure S4 :

 S4 Figure S4: Time series of anomalous gravity gradient signals over South America obtained by applying the same analysis as in Fig. 1a, for hypothetical earthquake times t e spanning the [January 2003 March 2010] time interval. Same scale and orientations of the spherical frame as in Fig. 1.

  )which is still relatively high. This is due to a larger contribution of the hydrological signals in the La Plata basin in this direction, in the 2 nd mode. The spatial and temporal patterns of this mode are consistent with those obtained in the North-South direction: the spatial pattern is localized in the epicentral region (Fig.S6e,f,m,n), and the temporal evolution shows a progressive step-like variation initiated months before the rupture, stabilized in March 2010 (Fig.S6h,p).

Figure S6 :

 S6 Figure S6: Common modes of variability in the region of Maule between the 800-km scale GRGS gravity gradients and the 800-km scale CSR (top and lower middle lines) or ITSG-2016 (upper middle and bottom lines) gravity gradients respectively, in the North-South direction (φφ gradients in the local spherical frame, top two lines) and in the East-West direction (θθ gradients in the local spherical frame, bottom two lines). Columns 1 and 2: non-dimensionalized spatial pattern of the rst common mode for each gravity gradient solution ; column 3: percentage of covariance explained ; column 4: dimensionalized time series of the rst common mode for each gravity gradient solution. Blue dot in the time series: January 2010 ; green dot: February 2010 ; pink dot: March 2010.

  consecutive segments separated by a free step in March 2010: [January 2003 -t 1 ], [t 1 -February 2010], [March 2010 -t 2 ] and [t 2 -September 2014] with t 1 = July 2009 and t 2 =

Acknowledgments

This work was supported by CNES. It is based on observations from the GRACE satellites. M.B. was funded by a CNES-IPGP fellowship. We thank Kerstin Schulze, Hannes Müller Schmied and Petra Döll for providing us access to the WGHM hydrology model. We would like to thank Ayan Fleischman, Nicolas Delbart, Emilie Lavie and Jean-Michel Lemoine for their kind support and constructive discussions. Most gures of this paper were generated using the free Generic Mapping Tools software (GMT) [50]. This is IPGP contribution number XXX.

D.2 Singular Value Decomposition of coupled elds D.2.1 Principle

To identify common spatio-temporal patterns between the three gravity models (GRGS, ITSG and CSR) taken two-by-two, we used a singular value decomposition (SVD). This method is well-suited to identify the coupled space-time patterns of variability between two elds. It is based on the decomposition of the cross-covariance matrix of the two space-time data matrices into a linear combination of orthogonal modes, each expressed by the multiplication of a spatial pattern with a time series. In more detail, the principle is as follows. We rst construct the temporal cross-covariance matrix (C) between two data elds e.g. GRGS (G) and ITSG (I). Each data eld is represented by a rectangular n × p matrix, where n is the number of epochs and p the number of locations (grid points in the case of a regular grid), such that G ij = g GRGS (t i , r j ) for the GRGS03 solution, I ij = g ITSG (t i , r j ) for the ITSG solution. Here, t i denotes the i-th epoch, r j the j-th position on the spatial grid, g GRGS (resp. g ITSG ) denotes the GRGS (resp. ITSG) gravity gradients. We have:

Then, we compute the SVD of the cross-covariance matrix C by solving the following equation, which can be seen as a generalization to rectangular matrices of the diagonalization of the square symmetric matrix:

We obtain two sets of spatially orthogonal singular vectors (the columns of U and V for matrices G and I respectively). The diagonal matrix L contains the singular values associated with each pair of singular vectors. Each common mode is represented by the product between a temporal mode a(t) and its associated spatial mode b(r). The i-th temporal mode is given by the i-th column of the matrix A = GU (resp. B = IV) for the GRGS data eld (resp. the ITSG data eld). The associated spatial patterns are given by the i-th columns of U and V respectively. Finally, the importance of each common mode When comparing the hydrological models, we noticed dierences between ERA5 and the other models. The ERA5 water storage amplitude appeared about twice larger than that of all the other models, hence the scaling by a factor of 0.5 applied for comparisons.

G.3 The seasonal cycle in ERA5-Land

In the Andean Cordillera close to the Maule region, contrary to the other hydrological models, the amplitude of its annual cycle varies irregularly by a factor up to 2.5 in the North-South gravity gradients, making the seasonal correction dicult (see Figure S10). This is why the RMS map of ERA5 in Fig. S8e comprises a small signal near Maule in the North-South direction (also present in the 2003-2009 period, see Appendix Fig. S9).

Nevertheless, this contribution cannot explain the GRACE anomaly before the earthquake due to dierent spatial and temporal patterns. The geometry of this ERA5 signal indeed follows the topographic reliefs of the Andes, leading to an absence of signal in the East-West direction and at the larger 800-km scale in the North-South direction. In addition, its time evolution is roughly periodic.


Sineewanlaya Wichit

Pauline Ferraris

Vale ´rie Choumet

Dorothe ´e Misse
email: dorothee.misse@ird.fr

The effects of mosquito saliva on dengue virus infectivity in humans

Arboviruses such as Dengue, Chikungunya, and Zika viruses represent a major public health problem due to globalization and propagation of susceptible vectors worldwide. Arthropod vector-derived salivary factors have the capacity to modulate human cells function by enhancing or suppressing viral replication and, therefore, modify the establishment of local and systemic viral infection. Here, we discuss how mosquito saliva may interfere with Dengue virus (DENV) infection in humans. Identification of saliva factors that enhance infectivity will allow the production of vector-based vaccines and therapeutics that would interfere with viral transmission by targeting arthropod saliva components. Understanding the role of salivary proteins in DENV transmission will provide tools to control not only Dengue but also other arboviral diseases transmitted by the same vectors.

Introduction

Mosquitoes are blood-feeding arthropods and, therefore, vectors of causative agents of important emerging infectious diseases. They are able to transmit several types of micro-organisms associated with vector borne disease, that cause more than 1 million deaths each year (WHO, 2016). Dengue virus (DENV) is the most important arbovirus in humans, infecting about 390 million individuals globally, with up to 96 million symptomatic infections. DENV infection results in a spectrum of illnesses, ranging from a flu-like disease (dengue fever) to dengue hemorrhagic fever (DHF) that can progress to dengue shock syndrome (DSS) and death. There are over 500 000 cases of DHF/ DSS reported worldwide each year, 25 000 of which lead to death, mostly among young children under the age of 15 [START_REF] Bhatt | The global distribution and burden of dengue[END_REF]. Epidemics with a high frequency of DHF/DSS continue to expand geographically in Asia and South America, highlighting the need to control the spread of this disease [START_REF] Guzman | Dengue: a continuing global threat[END_REF]. To date, there are no vaccines or effective treatments available for DENV infection.

DENV is a single-stranded RNA virus belonging to the Flaviviridae family. The four dengue serotypes (DENV-1 to DENV-4) are transmitted to humans by the Aedes (Ae.) mosquitoes species (Ae. aegypti and Ae. albopictus), which are spreading throughout the world due to climate change and globalization [START_REF] Gubler | Dengue, urbanization and globalization: the unholy trinity of the 21(st) century[END_REF]. A dramatic increase in urbanization has also created ideal conditions for increased transmission of mosquito-borne dengue disease. Outbreaks of DENV have generally been restricted to resource-limited areas of tropical and subtropical latitudes, in which 40% of the global population live. However, due to the increasing range of the DENV vectors, the virus could spread outside of these regions [START_REF] Da Silva-Voorham | Dengue: a growing risk to travellers to tropical and sub-tropical regions[END_REF].

Understanding the transmission cycle of DENV is important for the control of the spread of this disease. The cycle is initiated when pathogen-containing fluids are ingested by the vector from an infected vertebrate during a blood meal. Once the virus crosses the midgut barrier and has replicated in the mosquito, it reaches the salivary glands which leads to high levels of viral titers in the saliva of infected arthropods [START_REF] Luplertlop | Induction of a peptide with activity against a broad spectrum of pathogens in the Aedes aegypti salivary gland, following Infection with Dengue Virus[END_REF][START_REF] Salazar | Dengue virus type 2: replication and tropisms in orally infected Aedes aegypti mosquitoes[END_REF][START_REF] Vazeille | Orally co-Infected Aedes albopictus from La Reunion Island, Indian Ocean, can deliver both dengue and chikungunya infectious viral particles in their saliva[END_REF]. An infected mosquito will then transmit the virus to a host by delivering most of the viral content through the saliva into both the epidermis and dermis where resident and migratory cells encounter the pathogen [START_REF] Briant | Role of skin immune cells on the host susceptibility to mosquito-borne viruses[END_REF].

Though it is known that the DENV contained in salivary glands is closely associated with mosquito saliva, little is known about the function of the Aedes mosquito's salivary proteins and their role in the transmission of the virus.

Mosquito saliva helps to successfully obtain a blood meal [9] and could also contribute to an optimal transmission of the pathogen by modulating the host immune response [START_REF] Valenzuela | High-throughput approaches to study salivary proteins and genes from vectors of disease[END_REF]. It is as yet unknown, however, how salivary proteins may affect the human immune system during the transmission of the virus. This review will summarize the effect of mosquito saliva and salivary proteins on the transmission of DENV in humans.

Human innate immune response against DENV

Human skin cells represent the first physical and immunological barrier after being infected by DENV through a mosquito bite [START_REF] Briant | Role of skin immune cells on the host susceptibility to mosquito-borne viruses[END_REF]. Skin cells may contribute to the proinflammatory and antiviral microenvironment that is established in the early stages of DENV infection. Pattern recognition receptors (PRRs) such as membranebound toll like receptors (TLRs), cytoplasmic helicases (retinoic-acid-inducible gene I (RIG-I) and melanoma differentiation-associated gene 5 (MDA5)) and protein kinase R (PKR) are the first lines of defense in the immune recognition of DENV [START_REF] Tsai | Human TLR3 recognizes dengue virus and modulates viral replication in vitro[END_REF][START_REF] Loo | Distinct RIG-I and MDA5 signaling by RNA viruses in innate immunity[END_REF][START_REF] Nasirudeen | RIG-I, MDA5 and TLR3 synergistically play an important role in restriction of dengue virus infection[END_REF]. These PRRs trigger signaling cascades that induce the activation of transcription factors including interferon regulatory factors (IRF3 and IRF7). This will lead to the production of type I interferons (IFN-a/b) resulting in the expression of hundreds of interferon-stimulated genes (ISGs) with various levels of antiviral activity [START_REF] Chen | The roles of IRF-3 and IRF-7 in innate antiviral immunity against dengue virus[END_REF]. TLR3, RIG-I, MDA-5 and PKR are induced during DENV infection as well [START_REF] Surasombatpattana | Dengue virus replication in infected human keratinocytes leads to activation of antiviral innate immune responses[END_REF]. DENV infection of human fibroblasts and keratinocytes induces the production of antimicrobial peptides such as LL37, b-defensin, RNase 7 and S100 family members, that are essential components of the innate immune response (Figure 1) [START_REF] Bustos-Arriaga | Activation of the innate immune response against DENV in normal non-transformed human fibroblasts[END_REF][START_REF] Surasombatpattana | Aedes aegypti saliva enhances dengue virus infection of human keratinocytes by suppressing innate immune responses[END_REF]. Dendritic cells (DCs) are sentinels of the immune system that 140 Virus-vector interactions detect pathogens at sites of entry such as the skin. Following this initial contact with the virus, DCs mature and migrate to the draining lymph nodes, becoming effective stimulators of T cell responses [START_REF] Banchereau | Dendritic cells and the control of immunity[END_REF]. In addition to the ability of DCs to control infections directly via their innate immune functions, these cells can also become infected [START_REF] Navarro-Sanchez | s P: Dendritic-cellspecific ICAM3-grabbing non-integrin is essential for the productive infection of human dendritic cells by mosquitocell-derived dengue viruses[END_REF]. Hence, DENV was shown to infect epidermal Langerhans cells (LCs) in healthy human skin explants in vitro [START_REF] Wu | Human skin Langerhans cells are targets of dengue virus infection[END_REF][START_REF] Limon-Flores | Dengue virus inoculation to human skin explants: an effective approach to assess in situ the early infection and the effects on cutaneous dendritic cells[END_REF]. Another cell type that plays an important role in host defense against viruses through its capacity to modulate both innate and adaptive immune responses is the mast cell [START_REF] Galli | Mast cells as sentinels of innate immunity[END_REF][START_REF] Abraham | Mast cells in infection and immunity[END_REF][START_REF] Henz | Mast cells as initiators of immunity and host defense[END_REF]. Following infection with DENV mast cells have been shown to produce a spectrum of chemokines that are important for the trafficking of natural killer cells, monocytes and T cells [START_REF] Brown | RNA sensors enable human mast cell anti-viral chemokine production and IFN-mediated protection in response to antibody-enhanced dengue virus infection[END_REF].

The RNA interference (RNAi) pathway has also been shown to exert antiviral activity against DENV. The knockdown of Dicer, Drosha, Argo1, and Argo2 proteins of the RNAi pathway during DENV infection of mammalian cells results in increased DENV replication [START_REF] Kakumani | Role of RNA interference (RNAi) in dengue virus replication and identification of NS4B as an RNAi suppressor[END_REF].

However, DENV has developed several ways to directly inhibit cellular signaling cascades by evading innate antiviral responses [START_REF] Green | Innate immunity to dengue virus infection and subversion of antiviral responses[END_REF].

Modulation of vertebrate immune response by mosquito saliva

It has been shown in the literature that mosquito saliva is able to modify the host immune response, although the underlying molecular mechanisms are unclear as yet.

Results from initial studies, carried out with C3H/HeJ mice on which Ae. aegypti fed for a blood meal, showed a significant down regulation of IFN-g production while that of (interleukin (IL) IL-4 and IL-10) was up-regulated [START_REF] Zeidner | Mosquito feeding modulates Th1 and Th2 cytokines in flavivirus susceptible mice: an effect mimicked by injection of sialokinins, but not demonstrated in flavivirus resistant mice[END_REF]. Further studies explored the effect of Aedes saliva or salivary proteins on the transmission of arboviruses. One study reported that viral replication of Rift Valley Fever Virus, as well as the pathogenesis associated with infection with this virus was strongly enhanced by factors contained in saliva and salivary gland extract (SGE) of Aedes mosquitos [START_REF] Coupanec | Aedes mosquito saliva modulates Rift Valley fever virus pathogenicity[END_REF]. Similarly, in the context of DENV infection, co-infection with A. aegypti saliva, saliva proteins or SGE was also reported to boost viral replication [17,30,31 ,32 ,33 ].

Pathological investigation in fatal human cases revealed spleens with significant destruction of the germinal centers and atrophy of lymphoid follicles, which may be associated with a decrease in the frequency of T cells [START_REF] Po ´voa | The pathology of severe dengue in multiple organs of human fatal cases: histopathology, ultrastructure and virus replication[END_REF]. Infection of AG129 mice with dengue serotypes triggered gross pathology and histopathological changes in spleens [START_REF] Sarathy | A lethal murine infection model for dengue virus 3 in AG129 mice deficient in type I and II interferon receptors leads to systemic disease[END_REF]. And in a study using murine splenocytes, a reduction of T-cell and splenocyte proliferation as well as cytokine production (Th1 (IFN-g) and Th2 (IL-10)) was observed in the presence of Aedes SGE [START_REF] Wanasen | Differential modulation of murine host immune response by salivary gland extracts from the mosquitoes Aedes aegypti and Culex quinquefasciatus[END_REF]. Production of IFN-g was also observed in inbred A/J mice infected intravenously with DENV-2 [START_REF] Shresta | Early activation of natural killer and B cells in response to primary dengue virus infection in A/J mice[END_REF]. The effects of Ae. aegypti saliva on DENV infection in human keratinocytes have also been investigated. Both the saliva and recombinant saliva proteins were associated with an increase of mRNA level in the cells and a decrease of type I IFNs, IRF-3 and IRF-7 (Figure 1) [START_REF] Surasombatpattana | Aedes aegypti saliva enhances dengue virus infection of human keratinocytes by suppressing innate immune responses[END_REF][START_REF] Surasombatpattana | Aedes aegypti saliva contains a prominent 34-kDa protein that strongly enhances dengue virus replication in human keratinocytes[END_REF].

It has furthermore been shown that SGE of Ae. aegypti inhibits the release of tumor necrosis factor alpha (TNFa) from rat mast cells, but has no effect on antigeninduced histamine secretion [START_REF] Bissonnette | Extracts of mosquito salivary gland inhibit tumour necrosis factor alpha release from mast cells[END_REF]. A previous study demonstrated that Ae. aegypti SGE also does not affect the viability of a fetal skin-derived DC line, nor the production of IL-12 by these cells [START_REF] Wasserman | Saliva of the Yellow Fever mosquito, Aedes aegypti, modulates murine lymphocyte function[END_REF]. Consequently, Ae. aegypti SGE has no effect on the basal expression of IFN-b by DCs, but instead decreases the production of this cytokine in West Nile virus-infected cells [START_REF] Schneider | Aedes aegypti saliva alters leukocyte recruitment and cytokine signaling by antigen-presenting cells during West Nile virus infection[END_REF]. There is reportedly no effect of Ae. aegypti SGE on the maturation and function of mice dendritic cells [START_REF] Bizzarro | Effects of Aedes aegypti salivary components on dendritic cell and lymphocyte biology[END_REF]. On the other hand, Ader et al. found that Ae. aegypti saliva had an effect on myeloid DCs by inhibiting their susceptibility to infection by DENV and by enhancing their production of both IL-12p70 and TNF-a (Figure 1). Moreover, presensitization of DCs with saliva prior to infection enhanced inhibition of DENV replication [START_REF] Ader | Modulation of dengue virus infection of dendritic cells by Aedes aegypti saliva[END_REF]. Variations in the results of these two observations may thus be due to the origin of the DCs, that is, whether they belong to the mouse or human lineage. Together, these observations suggest that Ae. aegypti saliva is involved in the regulation of both innate and adaptive immunity in vertebrate hosts.

Aedes mosquito saliva proteins involved in DENV infectivity

Since Ae. aegypti saliva can affect the transmission and the establishment of arboviruses in humans [START_REF] Briant | Role of skin immune cells on the host susceptibility to mosquito-borne viruses[END_REF], much effort has been devoted to characterize the sialome and sialotranscriptome of this vector [START_REF] Arca | An insight into the sialome of the adult female mosquito Aedes albopictus[END_REF]. The salivary gland protein profile from the Ae. aegypti mosquito was evaluated using a proteomic approach and approximately one hundred proteins were identified in this organ [START_REF] Almeras | Salivary gland protein repertoire from Aedes aegypti mosquitoes[END_REF]. Ae. aegypti saliva contains many pharmacologically active proteins that can facilitate blood feeding, pathogen transmission and modulate host immune responses [START_REF] Schneider | The enhancement of arbovirus transmission and disease by mosquito saliva is associated with modulation of the host immune response[END_REF][START_REF] Machain-Williams | Association of human immune response to Aedes aegypti salivary proteins with dengue disease severity[END_REF][START_REF] Ribeiro | An annotated catalogue of salivary gland transcripts in the adult female mosquito, Aedes aegypti[END_REF] (Table 1). A functional genomic and proteomic analysis of salivary glands of Ae. aegypti mosquitoes, identified four proteins that strongly enhance DENV replication in human keratinocytes: serpin-like protein (AT), adenosine deaminase (AD), a serine protease with a MW of 34 kDa and a putative secreted protein (VA) [START_REF] Surasombatpattana | Aedes aegypti saliva contains a prominent 34-kDa protein that strongly enhances dengue virus replication in human keratinocytes[END_REF]. All of these proteins have been shown to down regulate the expression of IRF-7, but only the 34 kDa protein was able to inhibit the transcriptional expression of IRF-3, resulting in a stronger abrogation of type I IFN production [START_REF] Surasombatpattana | Aedes aegypti saliva contains a prominent 34-kDa protein that strongly enhances dengue virus replication in human keratinocytes[END_REF]. These results are in line with a study showing that the combined action of IRF-3 and IRF-7 is required for efficient control of viral replication [START_REF] Chen | The roles of IRF-3 and IRF-7 in innate antiviral immunity against dengue virus[END_REF]. The mechanisms by which AT, AD and VA increase DENV infection in keratinocytes have not yet been studied. The 34 kDa protein also decreases the expression of antimicrobial peptides in the context of DENV infection, which may also explain the increase in viral replication in human keratinocytes. Another study also identified a serine protease, CLIPA3, responsible for DENV transmission enhancement [START_REF] Conway | Mosquito saliva serine protease enhances dissemination of dengue virus into the mammalian host[END_REF]. This protein is known for breaking down the extracellular matrix, which may lead to cell migration at the inoculation site.

Also, certain saliva factors may enhance, while others may inhibit, virus infectivity due to the complex nature of mosquito saliva. Aegyptin is a collagen-binding protein that inhibits platelet aggregation and a member of the GE-rich 30 kDa antigen family present in Ae. aegypti saliva [START_REF] Calvo | Aegyptin, a novel mosquito salivary gland protein, specifically binds to collagen and prevents its interaction with platelet glycoprotein VI, integrin alpha2beta1, and von Willebrand factor[END_REF]. Mice inoculated with Aegyptin, show a decrease in DENV titers, both at the inoculation sites and in the circulation, as well as increased serum concentrations of GM-CSF, IFN-g, IL-5, and IL-6, as compared to mice that were inoculated with DENV alone. These results suggest that Aegyptin impacts DENV infectivity through activation of the immune response [START_REF] Mccracken | Aedes aegypti salivary protein ''aegyptin'' co-inoculation modulates dengue virus infection in the vertebrate host[END_REF]. Reducing the abundance of Aegyptin through transgene-mediated dsRNA interference results in prolonged probing time and less successful uptake of a blood meal [START_REF] Chagas | Collagen-binding protein, Aegyptin, regulates probing time and blood feeding success in the dengue vector mosquito, Aedes aegypti[END_REF]. It has also been shown that a cecropin-like peptide in the salivary glands of Ae. aegypti decreased DENV infectivity in the C6/36 Ae. albopictus cell line [START_REF] Luplertlop | Induction of a peptide with activity against a broad spectrum of pathogens in the Aedes aegypti salivary gland, following Infection with Dengue Virus[END_REF].

Among the proteins up-regulated in DENV infected Ae. aegypti SGE, D7 proteins were shown to inhibit platelet aggregation, vaso-constriction and inflammation [START_REF] Calvo | Function and evolution of a mosquito salivary protein family[END_REF]. A 36 kDa member of the D7 protein family was found to be associated with dengue disease severity [START_REF] Machain-Williams | Association of human immune response to Aedes aegypti salivary proteins with dengue disease severity[END_REF]. The role of D7 proteins in DENV transmission, however, remains largely unknown. Another major protein found in Ae. aegypti saliva, Apyrase, which inhibits ADP-dependent platelet aggregation during blood-feeding and prevents neutrophils activation [START_REF] Sun | Expression of functional recombinant mosquito salivary apyrase: a potential therapeutic platelet aggregation inhibitor[END_REF], has also been largely uncharacterized in its role in DENV infection. SAAG-4 protein and sialokinin demonstrated the capacity to modulate the ratio of Th1/Th2 cytokine production by suppressing and simultaneously enhancing the expression of IFN-g and IL-4, respectively, outside the context of viral infection [START_REF] Zeidner | Mosquito feeding modulates Th1 and Th2 cytokines in flavivirus susceptible mice: an effect mimicked by injection of sialokinins, but not demonstrated in flavivirus resistant mice[END_REF][START_REF] Boppana | SAAG-4 is a novel mosquito salivary protein that programmes host CD4 T cells to express IL-4[END_REF]. The role of these proteins in DENV transmission remains to be studied.

Table 1 provides a list of saliva proteins that were identified in several studies. However, further investigation is needed to better understand their effect on the DENV infectivity in humans. It is of note that, even if all major epidemics of DHF have occurred only in areas where Ae. aegypti is found, Ae. albopictus is also a competent vector for DENV and has the ability to adapt to diverse environmental conditions in both tropical and temperate regions [START_REF] Lambrechts | Consequences of the expanding global distribution of Aedes albopictus for dengue virus transmission[END_REF]. To date, there have been no studies on the effect of salivary proteins of this vector in the transmission of DENV.

Conclusion

Arboviral diseases have become a great concern due to the expanding geographical range of many mosquito species. Unpredictable outbreaks of arboviral diseases, such as the recent expansion of ZIKA virus, highlight the need to control the spread of their vectors. DENV is arguably the most important arbovirus causing disease in the world today, and despite decades of research, a vaccine is not yet available. Increased understanding of the activity of saliva proteins on DENV infectivity could provide novel targets that may aid in the development of vector-protein based vaccines and therapeutics.

Current

  Figure 1

Table 1

 1 Identified proteins in Aedes aegypti and Aedes albopictus saliva

	142 Virus-vector interactions
	Name family	Gene Bank	Activity	DENV infectivity	Immunomodulation	Reference
		accession number
	Apyrase
	68 kDa	AAEL006347	Hydrolase	ND	ND	[46,47,56]
	AD	AAEL005672	Adenosine deaminase	Enhanced	Suppressed IRF7,	[31 ,47,55]
					IFN type-I and AMP
	D7
	D7L1	AAEL006417	Odorant binding	ND	ND	[56]
	D7s2	AAEL006423	Odorant binding	ND	ND	[44,51,56]
	D7L2/37kDA	AAEL006424	Odorant and biogenic	ND	ND	[44,46,51]
			amines binding
	Serine protease
	34 kDa	Q1HRW0_AEDAE	Serine protease	Enhanced	Suppressed IRF7,	[47,31 ]
					IRF3, IFN type-I and AMP
	CLIPA3	AAEL005718	Serine-type	Enhanced	Not changed Th1	[48 ]
			endonuclease		and Th2 expression
	Aegyptin
	30 kDa salivary gland	AAEL010235	Collagen binding	Inhibited	Increased GM-CSF,	[32 ,44,
	allergen				IFN-g, IL-5, and IL-6	46,47,49]
	Aed a 3
	Serpin-like protein
	serpin-4	AAEL002704	Serine protease inhibitor	ND	ND	[44]
	SRPN26	AAEL003182	ND	ND	ND	[56]
	Antithrombin (AT)	Q1HRTV7_AEDAE	Thrombin inhibitor	Enhanced	Suppressed IRF7,	[31 ]
					IFN type-I and AMP
	Vitellogenin-B	AAQ92367	Lipid transporter	ND	ND	[44]
	Vitellogenin-C	AAQ92366	Lipid transporter	ND	ND	[44]
	SGS1	AAV28546	ND	ND	ND	[44,47]
	34 kDa	AAL76018, Q8T9VI_AEDAE	ND	ND	ND	[44,47]
	Venom Allergen (VA)	Q8T9U5_AEDAE	ND	Enhanced	Suppressed IRF7,	[31 ]
					IFN type-I and AMP
	Cecropin
	3.8 kDa	AAEL000598	Antibacterial activities	Inhibited	ND	[5,55]
	Sialokinin	AAEL000229	Vasodilatory	ND	ND	[28,56]
	ND: not determined.

Mosquito saliva and dengue transmission Wichit et al. 141

Acknowledgements

We thank Dr Hans Yssel for critical discussions and Thanyakorn Rongsawat for technical assistance. This work was supported by grants from the Agence Nationale de la Recherche France (grants ANR-12-BSV3-0004-01 and ANR-14-CE14-0029). Sineewanlaya Wichit was supported by a fellowship from the Infectiopo ˆle Sud foundation.