Louisa L H Tsang
email: louisa.tsang@geoazur.unice.fr

Mathilde Vergnolle

Cedric Twardzik

Anthony Sladen

Jean- Mathieu Nocquet

Frédérique Rolandone

Hans Agurto-Detzel

Olivier Cavalié

Paul Jarrin

Patricia Mothes

Imaging rapid early afterslip of the 2016 Pedernales earthquake, Ecuador

Keywords: afterslip, postseismic, subduction, Ecuador, aftershocks, GPS 43

High-Rate (HR) GPS time series following the 2016 M w 7.8 Pedernales earthquake suggest significant postseismic deformation occurring in the early postseismic period (i.e. first few hours after the earthquake) that is not captured in daily GPS time series. To understand the characteristics of early postseismic deformation, and its relationship with the mainshock rupture area, aftershocks and longer-term postseismic deformation, we estimate the spatiotemporal distribution of early afterslip with HR-GPS time series that span ~ 2.5 minutes to 72 hours after the earthquake, and compare with afterslip models estimated with daily GPS time series spanning a similar postseismic time period and up to 30 days after the earthquake. Our inversion technique enables us to image the nucleation of afterslip in the initial hours after the earthquake, bringing us closer to the transition between the coseismic and postseismic phases.

The spatial signature of early afterslip in the region updip of the mainshock rupture area is consistent with longer-term afterslip that occurs in the 30-day postseismic period, indicating that afterslip nucleated updip of and adjacent to peak coseismic slip asperities, in two localized areas, and subsequently continued to grow in amplitude with time in these specific areas. A striking difference, however, is that inversion of the 72-hour HR-GPS time series suggests early afterslip within the mainshock rupture area, but which may have been short-lived. More interestingly, we find that postseismic slip starts immediately after the earthquake at a rapid rate. Indeed, we find that early afterslip represents a significant contribution to the postseismic geodetic moment -afterslip in the first 72 hours is ~60 % greater than that estimated with daily GPS time series that span the first three post-earthquake daily GPS positions (i.e. covering the same time window). The results of our study demonstrate that imaging the spatio-temporal evolution of afterslip using subdaily GPS time series is important for evaluating postseismic slip budgets, and provides additional insights into the postseismic slip behaviour of faults. 42

Introduction

Afterslip plays an important role in redistributing stresses on and around megathrusts following subduction zone earthquakes, and can represent a significant amount of moment release [e.g. [START_REF] Hsu | Frictional Afterslip Following the 2005 Nias-Simeulue Earthquake, Sumatra[END_REF]. Thus, mapping out at various timescales where it occurs is vital for assessing seismic hazard on neighbouring sections of megathrusts. Most afterslip imaging studies are conducted using daily GPS time series spanning timescales ranging from weeks to years following the earthquake [e.g. [START_REF] Hsu | Frictional Afterslip Following the 2005 Nias-Simeulue Earthquake, Sumatra[END_REF][START_REF] Tsang | Afterslip following the 2007 M w 8.4 Bengkulu earthquake in Sumatra loaded the 2010 M w 7.8 Mentawai tsunami earthquake rupture zone[END_REF]. The first data point in these postseismic time series is at best taken on the first, or the second day, after the earthquake.

However, these postseismic daily GPS time series do not enable us to understand the nature of postseismic deformation occurring in the minutes to hours following the earthquake, hereafter referred to as the early postseismic period. Rather, early postseismic deformation is commonly captured within the static coseismic offset calculated with InSAR or daily GPS time series spanning a certain number of days before and after the earthquake. This manner of calculation can result in a contamination of the coseismic signal by early postseismic deformation, resulting in potential biases in the estimated coseismic source model. This likely explains why seismic moments for the same earthquake estimated from only daily GPS measurements can range from 1.5 to 2 times higher than those using only seismic data [e.g. [START_REF] Langbein | Coseismic and initial postseismic deformation from the 2004 Parkfield, California, earthquake, observed by Global Positioning System, Electronic Distance Meter, creepmeters, and borehole strainmeters[END_REF].

Furthermore, daily GPS time series do not allow us to resolve the transition between the coseismic and postseismic phases of deformation, and can hence potentially result in erroneous estimates of coseismic and postseismic slip budgets on faults. Thus, better understanding the early phase of postseismic deformation requires analysis of high-rate GNSS position time series that can capture the associated signal on the Earth's surface.

Thus far, only a handful of studies have focused on the early postseismic period. These studies focus on the 2004 M w 6.0 Parkfield earthquake in California [START_REF] Langbein | Coseismic and initial postseismic deformation from the 2004 Parkfield, California, earthquake, observed by Global Positioning System, Electronic Distance Meter, creepmeters, and borehole strainmeters[END_REF], 2003 M w 8.0 Tokachi-Oki [START_REF] Miyazaki | Coseismic and early postseismic slip for the 2003 Tokachioki earthquake sequence inferred from GPS data[END_REF][START_REF] Fukuda | Fault friction parameters inferred from the early stages of afterslip following the 2003 Tokachi-oki earthquake[END_REF] and 2011 M w 9.1 Tohoku-oki earthquakes [START_REF] Munekane | Coseismic and early postseismic slips associated with the 2011 off the Pacific coast of Tohoku Earthquake sequence: EOF analysis of GPS kinematic time series[END_REF] in Japan, 2012 M w 7.6 Nicoya earthquake in Costa Rica [Malservisi et al., 2012], and the 2010 M w 8.8 Maule, 2014 M w 8.3 Illapel and 2016 M w 7.8

Pedernales earthquakes in South America [START_REF] Twardzik | Unravelling the contribution of early postseismic deformation using sub-daily GNSS positioning[END_REF]. By analyzing HR-GPS position time series following these earthquakes, these studies demonstrate that the magnitude of early postseismic deformation is significant. For example, the geodetic moment in the early postseismic period following the 2012 Nicoya earthquake represents ~57 % of that released during the first 70 days [Malservisi et al., 2012]. Malservisi et al. [2012] and [START_REF] Munekane | Coseismic and early postseismic slips associated with the 2011 off the Pacific coast of Tohoku Earthquake sequence: EOF analysis of GPS kinematic time series[END_REF] estimated snapshots of the spatial distribution of early afterslip by inverting cumulative offsets calculated with a certain length time window of HR-GPS time series. However, this approach does not enable us to ensure coherency of afterslip from one time step to another. Also, it should be noted that cumulative offsets calculated with HR-GPS time series can vary depending on the level of noise in the time series and the method used to calculate these offsets. Of the studies mentioned above, only [START_REF] Miyazaki | Coseismic and early postseismic slip for the 2003 Tokachioki earthquake sequence inferred from GPS data[END_REF] inverted HR-GPS time series (spanning 4 hours after 2003 Tokachi-oki earthquake) to image the spatio-temporal evolution of afterslip. Their results show a complex evolution of early afterslip and suggest that it might have triggered the aftershock that occurred ~1.2 hours after the mainshock. They concluded that depth-dependent properties on the fault influence the propagation patterns of afterslip seen in their models. [START_REF] Fukuda | Fault friction parameters inferred from the early stages of afterslip following the 2003 Tokachi-oki earthquake[END_REF] investigated postseismic time series of this earthquake too, with a focus on the sudden acceleration of motions ~1.2 hours after the earthquake. They argue that the timing of this acceleration phase is driven by stress changes from the mainshock and the frictional parameters of the rate-and-state friction law, rather than due to the timing of the aftershock. Their results support theoretical studies demonstrating that faults experience an initial acceleration phase of afterslip, which is then followed by decelerating afterslip governed by steady-state velocity-strengthening friction [START_REF] Perfettini | Dynamics of a velocity strengthening fault region: Implications for slow earthquakes and postseismic slip[END_REF].

Clearly, investigating the spatio-temporal evolution of early postseismic deformation is relevant for understanding the physical mechanisms driving fault slip behaviour, as well as resolving its relationship with the mainshock rupture area, ensuing aftershocks, longer-term postseismic deformation, and quantifying its contribution to the postseismic slip budget.

Here, we model HR-GPS postseismic time series to estimate the spatio-temporal evolution of early afterslip following the 2016 M w 7.8 Pedernales earthquake (16 April 2016, 23:58:33) in Ecuador (Figure 1). The Ecuadorian megathrust hosts a diverse range of seismic and aseismic behaviour that has been studied in detail with regional GPS and seismic networks, as well as marine seismic studies [e.g., [START_REF] Collot | Subducted oceanic relief locks the shallow megathrust in central Ecuador[END_REF][START_REF] Font | Seismicity patterns along the Ecuadorian subduction zone: new constraints from earthquake location in a 3-D a priori velocity model[END_REF][START_REF] Vallée | Intense interface seismicity triggered by a shallow slow slip event in the Central Ecuador subduction zone[END_REF][START_REF] Mothes | Continuous GPS network operating throughout Ecuador[END_REF][START_REF] Chlieh | Distribution of discrete seismic asperities and aseismic slip along the Ecuadorian megathrust[END_REF][START_REF] Nocquet | Motion of continental slivers and creeping subduction in the northern Andes[END_REF]2016;[START_REF] Marcaillou | Seamount subduction at the North-Ecuadorian convergent margin: Effects on structures, interseismic coupling and seismogenesis[END_REF][START_REF] Rolandone | Areas prone to slow slip events impede earthquake rupture propagation and promote afterslip[END_REF][START_REF] Gombert | Strain budget of the Ecuador-Colombia subduction zone: A stochastic view[END_REF][START_REF] Vaca | Recurrent slow slip events as a barrier to the northward rupture propagation of the 2016 Pedernales earthquake (Central Ecuador)[END_REF][START_REF] Segovia | Seismicity distribution near a subducting seamount in the central Ecuadorian subduction zone, space-time relation to a slow-slip event[END_REF]Agurto-Detzel et al. (submitted); [START_REF] Twardzik | Unravelling the contribution of early postseismic deformation using sub-daily GNSS positioning[END_REF]. The optimal location of these networks have enabled detection of a number of slow slip events, repeating earthquakes and seismic swarms, as well as significant afterslip in the month following the Pedernales earthquake [START_REF] Vallée | Intense interface seismicity triggered by a shallow slow slip event in the Central Ecuador subduction zone[END_REF][START_REF] Rolandone | Areas prone to slow slip events impede earthquake rupture propagation and promote afterslip[END_REF][START_REF] Vaca | Recurrent slow slip events as a barrier to the northward rupture propagation of the 2016 Pedernales earthquake (Central Ecuador)[END_REF][START_REF] Segovia | Seismicity distribution near a subducting seamount in the central Ecuadorian subduction zone, space-time relation to a slow-slip event[END_REF]. These findings indicate that aseismic slip processes represent a significant mode of strain release along the megathrust, and it is therefore important to understand their contribution to megathrust slip budgets.

In this study, we capitalize on the availability of rich postseismic geodetic datasets from [START_REF] Twardzik | Unravelling the contribution of early postseismic deformation using sub-daily GNSS positioning[END_REF] and [START_REF] Rolandone | Areas prone to slow slip events impede earthquake rupture propagation and promote afterslip[END_REF] to image the spatio-temporal distribution of early afterslip in the first 72 hours after the earthquake, and compare this to afterslip estimated with daily GPS time series spanning the same time period (i.e. the first three post-earthquake daily GPS positions), as well as spanning the longer, 30-day post-earthquake time period (Figure 2). We also examine relationships between early afterslip, the earthquake rupture area, and the spatio-temporal distribution of relocated aftershocks.

Data

We use GPS data from 27 stations of the Ecuadorian continuous GPS network that were installed before the earthquake, and from which both HR-GPS and daily GPS time series were available (Figure 3). The horizontal components of the daily GPS time series from these stations were processed by [START_REF] Rolandone | Areas prone to slow slip events impede earthquake rupture propagation and promote afterslip[END_REF], with positions expressed with respect to the first daily GPS position (the day after the earthquake). In the days following the earthquake, 4 M w > 6 aftershocks occurred in the region updip of the mainshock rupture area. Coseismic offsets (and likely some ensuing postseismic deformation) associated with these aftershocks are visible in the time series of stations CABP and MOMP on the 22 and 20 April, respectively (Figure S1).

We did not attempt to correct the time series for these offsets, due to the absence of a technique to reliably estimate and remove the coseismic and postseismic signals associated with these aftershocks. No visible offsets associated with these aftershocks are observed at other stations.

The 30-second HR-GPS time series from these stations were processed by [START_REF] Twardzik | Unravelling the contribution of early postseismic deformation using sub-daily GNSS positioning[END_REF]. The time series span ~72 hours after the earthquake from 17 April 2016, 00:02:00 to 19 April 2016, 23:59:59. The first position in the HR-GPS time series was chosen to ensure that the time series are not contaminated by seismic waves associated with the mainshock (see [START_REF] Twardzik | Unravelling the contribution of early postseismic deformation using sub-daily GNSS positioning[END_REF] for details). We downsampled these time series to a position every 60 seconds, and used only the horizontal component of displacements, as the vertical component of displacements have relatively low signal-to-noise ratios and large uncertainties.

We examined the gCMT catalog of regional earthquakes, and found no significant signals of these earthquakes in the HR-GPS time series. The largest 72-hour cumulative displacements are observed at stations nearest the mainshock rupture area (Figure 3): MOMP, CABP and PDNS, with values of 5.5, 4.6 and 3.3 cm, respectively. The 30-day cumulative displacements for these stations, calculated from the daily GPS time series, are 12.6, 11 and 7.2 cm, respectively.

The 60-second HR-GPS time series are relatively noisy (Figures S2,S3). Hence, to verify that our inversion results are not significantly influenced by noise in the time series, we performed four additional inversions with time series low-pass filtered at various cut-off frequencies, to extract the reliable features of the models (Text S1, Figure S3) (see next section for details).

Inversion method to estimate the spatio-temporal distribution of afterslip

We assume that the signal in the early postseismic time series represents predominantly afterslip on the megathrust. This approximation seems valid considering that postseismic signals are recorded at stations covering a region greater than would be expected from poroelastic effects, the latter of which are likely to occur in more localized regions [e.g. [START_REF] Tung | Delayed Poroelastic Triggering of the 2016 October Visso Earthquake by the August Amatrice Earthquake, Italy[END_REF]. Also, the early postseismic time period considered is shorter than characteristic relaxation times associated with typical mantle viscosities previously reported at subduction zones (assuming a linear mantle rheology, mantle viscosities of 10 17 to 5*10 19 Pa s correspond to relaxation times ranging from ~38 days to 50 years). Even so, we do not exclude the possibility that our modelling approach might represent deformation associated with other postseismic mechanisms mapped onto the fault. We estimated the spatio-temporal distribution of afterslip that span the first: (1) 72 hours (with HR-GPS time series), (2) 2 days (with time series of the first three daily GPS positions), and (3) 30 days (with daily GPS time series) (Figure 2). While [START_REF] Rolandone | Areas prone to slow slip events impede earthquake rupture propagation and promote afterslip[END_REF] already published a spatio-temporal 30-day afterslip model, we estimated the spatio-temporal 2-day and 30-day models using the same inversion strategy as that used to estimate the 72-hour model, in order to enable a fair comparison between them.

Since the HR-GPS time series are relatively noisy compared to daily GPS time series, a method is required to filter for the main postseismic signal in the presence of considerable noise in the time series. We adopted the Principal Component Analysis-based Inversion Method (PCAIM) [START_REF] Kositsky | Inverting geodetic time series with a principal component analysis-based inversion method[END_REF] to invert for the spatio-temporal distribution of afterslip, as this method allows the main postseismic signal to be represented with a small number of principal components. The surface displacement patterns associated with each principal component are inverted to obtain a principal slip distribution; these are then linearly combined together with the corresponding time functions to obtain a model of the spatio-temporal evolution of afterslip (Text S2). While PCAIM has been validated on standard daily GPS time series, to the best of our knowledge it has never been implemented yet with noisy HR-GPS time series. However, it has been shown to be effective in estimating afterslip with noisy daily GPS time series [e.g. [START_REF] Gualandi | Space-time evolution of crustal deformation related to the M w 6.3, 2009 L'Aquila earthquake (central Italy) from principal component analysis inversion of GPS position time-series[END_REF]. To determine the optimum number of principal components to represent the postseismic signal, we analysed the time series fits as incrementally more principal components were employed, examined the time functions, surface displacement patterns associated with each principal component, and the corresponding principal slip patterns.

For the HR-GPS time series, we found that one principal component adequately represents the postseismic signal (Text S2, Figure S4), while for the daily GPS time series, two principal components were required (Figure S5). The additional component for the latter case is probably due to the presence of a slow slip event in the daily GPS time series (discussed in section 5.6). Offsets in the daily GPS time series seen at stations CABP and MOMP mentioned previously are not present in the time functions of the first two principal components, when the data are decomposed. Therefore, inclusion of these offsets does not significantly influence the afterslip model results.

We estimated afterslip on the same fault interface geometry employed by [START_REF] Rolandone | Areas prone to slow slip events impede earthquake rupture propagation and promote afterslip[END_REF]. This fault geometry spans latitudes ~2°S to ~1°N, extends to a depth of 80 km, and follows the SLAB1.0 model along this section of the Andean megathrust [START_REF] Hayes | Slab1.0: A three-dimensional model of global subduction zone geometries[END_REF].

The fault interface was discretized into quasi-equilateral triangles with 10-km-long edges. The horizontal displacements at each GPS station due to slip on each sub-fault patch were calculated with Okada's solutions for deformation due to point dislocations embedded in a homogeneous elastic half-space [START_REF] Okada | Internal deformation due to shear and tensile faults in a half-space[END_REF]. Similar to the approach of [START_REF] Rolandone | Areas prone to slow slip events impede earthquake rupture propagation and promote afterslip[END_REF], the rake angle on each triangular patch was fixed to a direction that is consistent with the Nazca/North Andean block relative velocity determined by [START_REF] Nocquet | Motion of continental slivers and creeping subduction in the northern Andes[END_REF] (although we note that allowing the rake angle to vary by +/-45° from 90° does not significantly affect the results; Figure S6). Positivity constraints were applied to ensure trench-directed slip.

The Laplacian constraints, which control the degree that slip is smoothed across neighbouring patches, were weighted spatially based on the sensitivity of each station to slip on each sub-fault patch (following the method of Ortega-Culaciati [2013]). The spatial variation of this weighting results in more smoothing of slip in poorly resolved areas (Text S3, Figure S7). In addition, λ, the value controlling the strength of the Laplacian constraints, was selected for each model based on the standard L-curve analysis [START_REF] Hansen | Analysis of discrete ill-posed problems by means of the L-curve[END_REF], corroborated by visual examination of the slip distributions associated with each value of λ. The optimum value of λ represents the best trade-off between the data-model misfits and the model norm, and which preserves the key spatial characteristics of the solution (Text S3, Figure S8). We note that the optimum value of λ differs by three orders of magnitude between the case of inverting the HR-GPS and the daily GPS time series, because in the case of the daily GPS time series, we took into account the data uncertainties (although we retrieve a similar spatial distribution of slip even if we do not take into account the data uncertainties). On the other hand, we did not take into account the uncertainties of the HR-GPS data in the inversion, as these could not be reliably estimated.

As mentioned in the previous section, for the 72-hour dataset, we performed five inversions, with each inversion using time series with various degrees of filtering applied, in order to extract the reliable features of afterslip. Our reference model is the afterslip distribution estimated with time series filtered with a cut-off frequency corresponding to a time period of ~2.7 hours, chosen as it yields visually smooth time series (Text S1).

Results

Early afterslip distribution from the 72-hour dataset

Six main early afterslip areas are evident in all five models (Areas A-F, Figure 3, inset, and Figure S9). The five models indicate geodetic moments ranging from 5.16E+19 to 6.35E+19 N m (M w 7.11-7.17), assuming a rigidity of 30 GPa (Table S1).

Our reference model, shown in Figure 3, yields a reduced chi-square of 0.65. Assuming a rigidity of 30 GPa, the estimated total geodetic moment is 5.22E+19 N m, which is equivalent to a moment magnitude of M w 7.11. This moment estimate is ~9.7% of the coseismic moment as estimated by [START_REF] Nocquet | Supercycle at the Ecuadorian subduction zone revealed after the 2016 Pedernales earthquake[END_REF]. The model shows peak cumulative afterslip of ~ 30 cm concentrated in two areas updip of and adjacent to asperities that experienced the largest coseismic slip (Areas A and B). The largest cumulative geodetic moment occurred in area A (9.25E+18 N m; M w 6.61), while less occurred in area B (8.28E+18 N m; M w 6.58). In addition, ~13-18 cm of cumulative afterslip is concentrated within the southern portion of the mainshock rupture area (area C) and immediately downdip of the mainshock rupture area (area D), with each area amounting to an estimated geodetic moment of 1.02E+19 N m (M w 6.63). ~1.67E+18 N m (M w 6.12) of geodetic moment is estimated north of the mainshock rupture area (area E).

Further south-east, ~190 km from the Pedernales earthquake hypocentre, ~ 28 cm of cumulative slip is estimated on a single patch located north of La Plata Island (area F, Figure 3). This patch of slip is located close to the remotely triggered slow slip event that was reported by [START_REF] Rolandone | Areas prone to slow slip events impede earthquake rupture propagation and promote afterslip[END_REF]; we discuss these results in section 5.6.

While we demonstrate that the main afterslip areas shown in Figure 3 (inset) are reliable based on inverting various groups of filtered time series, we additionally conducted synthetic tests (in the next section) to assess the resolution of slip in our reference model, which is related to the fault geometry that we employed and the spatial distribution of GPS stations.

4.2.

Model resolution

For each main afterslip area, 1 m of trench-perpendicular slip was assigned to each sub-fault patch. Then, horizontal displacements at each station were calculated with the Okada model [START_REF] Okada | Internal deformation due to shear and tensile faults in a half-space[END_REF]. Synthetic time series at each station were subsequently calculated by multiplying the horizontal displacements with the synthetic normalized time function (based on the first principal component obtained from decomposition of the 30-day GPS time series). For each station, noise was added to the synthetic time series that is on the order of the average standard deviation from the mean position estimated using the pre-earthquake HR-GPS time series (Figure S2). For the inversion, we employed the same value of λ as that in our reference model.

The magnitude of peak slip is well recovered in the two areas updip of the mainshock rupture area (Figure 4a), as well as within the mainshock rupture area (Figure 4b), as the inversion indicates that ~ 73-78 % of the input synthetic slip is recovered. On the other hand, the magnitude of peak slip is not well recovered downdip of the mainshock rupture area (Figure 4c), and near Esmeraldas (Figure 4d), with only ~ 30-38 % of the input synthetic slip recovered in the inversion, and significant smoothing of slip across adjacent sub-fault patches. These conclusions are supported by the recovered slip model based on an input model where all main afterslip areas were assigned trench-perpendicular synthetic slip of 1 m (Figure 4e); slip in the updip regions is well recovered, while slip downdip and near Esmeraldas is not. In the mainshock rupture area, slip is recovered, but is smeared across neighbouring sub-fault patches.

Importantly, Figure 4a demonstrates that slip in the southern updip area is not smeared into the mainshock rupture area, suggesting that slip in area C (Figure 3) is not an artifact of the spatial smoothing that we employ. In addition, slip in area C is evident in rougher afterslip models from all five inversions (Figure S8), suggesting that this is a reliable feature, since it is evident regardless of the temporal smoothing of the time series and spatial smoothing of slip.

We present in Text S4 and Figure S11 results of implementing an alternative method of assessing the resolution of our model. The results substantiate our conclusions drawn here.

Furthermore, [START_REF] Rolandone | Areas prone to slow slip events impede earthquake rupture propagation and promote afterslip[END_REF] reported that slip areas of > 40 km on this fault geometry are well resolved.

Discussion

Comparison of 72-hour and 2-day afterslip models highlights the significant contribution of early afterslip

Figure 5 and 5c), where the difference in geodetic moment is a factor of ~43 and 2.6, respectively. However, our model resolution tests suggest that slip is not well resolved in these two areas. In other areas (areas A-D), the geodetic moment between the two models differs by a factor of ~1.5 to 2.3. These findings suggest enhanced afterslip in these areas during the first 12 hours, the period that is not captured by the daily GPS time series.

South of the mainshock rupture area (area G, Figure 5c), more afterslip is estimated in the 2day model. According to our approach of extracting the robust features of the 72-hour model based on various filtered time series (Figure S9), slip here is not a reliable feature and may be due to high frequency noise in the time series. Slip here in the 72-hour model is thus not accurately resolved.

Overall, the estimated geodetic moment of the 72-hour model is a factor of ~1.6 greater than that of the 2-day model (Table 1) (this factor is ~1.68 in the case of excluding poorly resolved slip areas), suggesting that not accounting for the postseismic deformation recorded immediately after the coseismic phase could result in an underestimation of the postseismic geodetic moment by ~60 % -very early afterslip (i.e. before the first GPS daily position) therefore contributes significantly to the postseismic geodetic moment.

Our 72-hour afterslip model explains why coseismic slip distributions estimated with InSAR data place coseismic slip near the trench. For example, from InSAR data that span 6 days after the earthquake, [START_REF] He | Coseismic slip in the 2016 M w 7.8 Ecuador earthquake imaged from Sentinel-1A radar interferometry[END_REF] found coseismic slip that extends to regions near the trench, and a higher estimated seismic moment compared with studies that either exclude or put a low weighting on InSAR data [e.g. [START_REF] Nocquet | Supercycle at the Ecuadorian subduction zone revealed after the 2016 Pedernales earthquake[END_REF][START_REF] Ye | The 16 April 2016, M W 7.8 (M S 7.5) Ecuador earthquake: A quasi-repeat of the 1942 M S 7.5 earthquake and partial rerupture of the 1906 M S 8.6 Colombia-Ecuador earthquake[END_REF]. As the InSAR data likely contain deformation due to both the coseismic and early postseismic slip, their models are consistent with our results showing early afterslip located updip of the mainshock rupture area.

Spatio-temporal evolution of afterslip and aftershocks in the first 72 hours

Figure 6 shows how afterslip evolves in time within the 72-hour period. As pointed out by [START_REF] Rolandone | Areas prone to slow slip events impede earthquake rupture propagation and promote afterslip[END_REF], the spatial distribution of shallow afterslip shows little correspondence to that of stress changes induced by coseismic slip. Interestingly, our results suggest that this spatial distribution, where slip is localized in a few specific areas rather than broadly updip and downdip of the rupture area, is in place immediately after the earthquake. In the subsequent 72hour postseismic period, the spatial distribution of afterslip remains fixed, with the amplitude of afterslip increasing with time. This latter finding is expected given that only one principal component was used (and is sufficient) to represent the postseismic time series, such that the slip evolution on each patch is governed by the same time function. This specificity of using one principal component also limits our capability to discern possible patterns of accelerating afterslip in the different afterslip areas. Importantly, our spatio-temporal modelling approach enables us to analyze an enriched picture of the temporal evolution of afterslip that indicates that the highlighted updip areas are particularly prone to host afterslip, and enables us to compare our results with the spatio-temporal evolution of aftershocks.

Aftershocks occur due to the release of static and/or dynamic stress changes associated with the mainshock [START_REF] Dieterich | A constitutive law for rate of earthquake production and its application to earthquake clustering[END_REF][START_REF] Felzer | Decay of aftershock density with distance indicates triggering by dynamic stress[END_REF][START_REF] Stein | The role of stress transfer in earthquake occurrence[END_REF]. In addition, a number of studies have reported similar temporal evolution of afterslip and aftershocks following large earthquakes, which support a model whereby aftershocks are produced when rate-weakening asperities are loaded by afterslip and driven to coseismic failure [e.g. [START_REF] Hsu | Frictional Afterslip Following the 2005 Nias-Simeulue Earthquake, Sumatra[END_REF][START_REF] Perfettini | Postseismic relaxation driven by brittle creep: A possible mechanism to reconcile geodetic measurements and the decay rate of aftershocks, application to the Chi-Chi earthquake, Taiwan[END_REF]. To explore the spatio-temporal relationship between afterslip and aftershocks in the first 72 hours, we analyzed the aftershock catalogue of Agurto-Detzel et al.

(submitted). Errors in the aftershock locations at the 68 % confidence level are ~12 and 13 km in the horizontal and vertical directions, respectively. Based on their one-year-long aftershock catalogue, Agurto-Detzel et al. (submitted) reported a threshold magnitude of completeness of 2.5. However, the catalogue is most likely incomplete in the initial hours after the earthquake, as during this period the noise level is higher and events are harder to identify. To be more conservative, we therefore analyzed M L ≥ 3.5 aftershocks in this study (Figure S12).

Figure 6 shows that relocated M3.5+ aftershocks appear to concentrate in regions bordering the two updip peak afterslip areas (A and B), which is consistent with the results from other studies that reported limited occurrence of aftershocks within the regions of peak afterslip [e.g. [START_REF] Hobbs | Large and primarily updip afterslip following the 2012 M w 7.6 Nicoya, Costa Rica, earthquake[END_REF][START_REF] Hsu | Frictional Afterslip Following the 2005 Nias-Simeulue Earthquake, Sumatra[END_REF]. In area A, the aftershock-afterslip moment ratio, expressed as a percentage of the cumulative aftershock seismic moment to the cumulative afterslip geodetic moment, is 11%, while in area B, this ratio is ~ 6% (seismicity in each area is defined using the limits shown in Figure S12). These ratios indicate that the majority of the postseismic deformation is aseismic. A large number of aftershocks are also located within the mainshock rupture area; these occur in the region between the two peak coseismic slip asperities, and border the patches of peak afterslip in area C. However, the cumulative seismic moment released in area C is small -the aftershock-afterslip moment ratio is 0.93 %; around an order of magnitude smaller than in the two updip afterslip regions. In contrast, there is sparse seismicity in the region downdip of the mainshock.

In the first 72 hours, the similarity in the temporal evolution of afterslip and numbers of M L > 3.5 aftershocks in the two updip areas (A and B) (Figure 7) suggests that aftershocks in these regions are likely driven by afterslip, consistent with the findings of Agurto-Detzel et al.

(submitted) based on their one-year-long catalogue. In contrast, the evolution curves suggest that in area C, there is a much larger fraction of aftershock occurrence in the first ~48 hours compared to other areas (Figure 7). We hypothesize that a proportion of aftershocks represent the release of residual coseismically-induced stress on the megathrust interface as well as in the overlying crust (Figure S12), and are thus possibly controlled by two different processes. In addition, though our analysis may be limited by the completeness of the aftershock catalogue in the first 72 hours, the evolution curves (note that Figure 7 is in semilogarithmic scale) suggest that early afterslip and aftershocks do not appear to follow an Omoritype law, in contrast to findings from analyzing aftershocks over a longer time period of ~1 year [Agurto-Detzel et al. (submitted)].

Interestingly, the 1942 earthquake has been inferred to have roughly the same rupture area as the 2016 earthquake [START_REF] Ye | The 16 April 2016, M W 7.8 (M S 7.5) Ecuador earthquake: A quasi-repeat of the 1942 M S 7.5 earthquake and partial rerupture of the 1906 M S 8.6 Colombia-Ecuador earthquake[END_REF][START_REF] Nocquet | Supercycle at the Ecuadorian subduction zone revealed after the 2016 Pedernales earthquake[END_REF], and aftershocks following the 1942 earthquake were mainly located seaward of the hypocentre [START_REF] Mendoza | Seismicity associated with the great Colombia-Ecuador earthquakes of 1942, 1958, and 1979: Implications for barrier models of earthquake rupture[END_REF],

implying that these aftershocks may have been driven by afterslip updip of the mainshock rupture area too. If so, these observations may indicate that aseismic slip behaviour, and by inference, frictional properties, persist through at least two seismic cycles.

30-day afterslip distribution from daily GPS time series

Our 30-day afterslip results (Figure 8b) are consistent with that estimated by Rolandone et al.

[2018] (data-model fits at representative stations are shown in Figure 8d; the spatio-temporal evolution of afterslip in the first 30 days and data-model time series fits at all stations are shown in Figures S13 andS14, respectively). Estimated afterslip in the first 30 days amounts to a geodetic moment of 1.71E+20 N m, equivalent to an ~Mw 7.46. In the region updip of the mainshock rupture area, the spatial distribution of 30-day afterslip is consistent with that in the first 72 hours, indicating that afterslip nucleated in localized areas updip of and adjacent to patches that experienced the largest coseismic slip, and subsequently continued to grow in amplitude with time. Similar to the case within the first 72 hours, aftershocks in the subsequent 28 days mostly surround the updip afterslip areas (Figure 8b). Peak afterslip is concentrated in area B (~70 cm) and might have extended to the trench, though the model resolution is poor close to the trench. In both areas A and B, estimated afterslip in the first 72 hours is ~30 % of that accumulated in the first 30 days (Table 1).

On the other hand, there are differences between the 72-hour and 30-day afterslip models within, downdip and north of the mainshock rupture area (areas C, D and E, respectively). In area E, we have already pointed out that resolution of slip here in the 72-hour model is poor (see section 4.2). In area C, where slip is better resolved, a greater amount of afterslip likely occurred in the first 72 hours, but this afterslip may have been short-lived, as limited afterslip is imaged in these areas in the 30-day model. Also, notwithstanding the poorer model resolution in the downdip region, a comparison of the two models shows that afterslip here nucleated in the first 72 hours immediately downdip of the mainshock rupture at ~ 30-50 km depths, and 404 may have subsequently migrated southwards and extended to ~40-70 km depths. 405

Spatial overlap of early afterslip with the mainshock rupture area

Comparison of our 72-hour, 2-day and 30-day models suggests that afterslip in the mainshock rupture zone rapidly decayed with time, likely within the first 12 hours after the earthquake.

There is a spectrum of observations concerning afterslip occurring in mainshock rupture zones. Some studies have reported observations of coseismic deformation that is anticorrelated with afterslip and/or aftershock areas [e.g. [START_REF] Perfettini | Seismic and aseismic slip on the Central Peru megathrust[END_REF][START_REF] Sladen | Source model of the 2007 M w 8.0 Pisco, Peru earthquake: Implications for seismogenic behavior of subduction megathrusts[END_REF]. These studies support the 'rate-state asperity model' that stipulates that velocity-weakening patches on the fault that host seismic ruptures do not manifest significant afterslip. On the other hand, a number of recent studies report that afterslip in coseismic rupture zones is needed to explain geodetic observations, challenging this conceptual model [e.g. [START_REF] Barnhart | Coseismic slip and early afterslip of the 2015 Illapel, Chile, earthquake: Implications for frictional heterogeneity and coastal uplift[END_REF][START_REF] Bedford | A high-resolution, time-variable afterslip model for the 2010 Maule M w = 8.8, Chile megathrust earthquake[END_REF][START_REF] Miyazaki | Coseismic and early postseismic slip for the 2003 Tokachioki earthquake sequence inferred from GPS data[END_REF][START_REF] Johnson | Challenging the rate-state asperity model: Afterslip following the 2011 M9 Tohoku-oki, Japan, earthquake[END_REF]Salman et al., 2017]. It is possible that some of these observations can be explained by modelling biases. For example, coseismic slip distributions estimated with InSAR data often span the early and longer postseismic period, which potentially results in some amount of postseismic deformation being mapped as coseismic slip. Another bias may lie in the choice of model smoothing parameters for postseismic slip models [e.g. [START_REF] Fukuda | Fault friction parameters inferred from the early stages of afterslip following the 2003 Tokachi-oki earthquake[END_REF]. Despite these biases, several physical mechanisms have been proposed to explain the overlap of afterslip and coseismic slip areas. [START_REF] Helmstetter | Afterslip and aftershocks in the rate-and-state friction law[END_REF] stipulate that frictional properties on the megathrust are not steady state. Physical processes such as dynamic weakening and rupture directivity may push normally velocity-strengthening patches to participate in coseismic ruptures [e.g. [START_REF] Noda | Stable creeping fault segments can become destructive as a result of dynamic weakening[END_REF]Lapusta, 2013, Salman et al. 2017]. Along the Ecuadorian subduction zone, [START_REF] Kanamori | Variable rupture mode of the subduction zone along the Ecuador-Colombia coast[END_REF] proposed a model in which distinct 'asperities ' of the 1942, 1958 and 1979 ruptures are separated by weak zones that typically behave aseismically, but can slip abruptly when driven by failure of neighbouring asperities -this mechanism was proposed to explain the coordinated rupture of all three asperities in a M w 8.8 earthquake in 1906. Indeed, frictional properties may vary during the course of the earthquake cycle, and patches that are conditionally stable may participate in both coseismic and aseismic slip phases. [START_REF] Yabe | Why do aftershocks occur within the rupture area of a large earthquake?[END_REF] modeled a frictionally heterogeneous fault system and their quasi-dynamic numerical simulation results were able to reproduce afterslip and aftershocks occurring around and within the mainshock rupture zone. They proposed that if patches within the mainshock rupture zone do not fully release the accumulated slip deficit, afterslip and aftershocks on these patches are a plausible mechanism through which the residual slip deficit can be released. In addition, Agurto-Detzel et al. (submitted) proposed that mainshockreactivated processes in the damage zone above the mainshock rupture area might explain the large density of aftershocks located on and above the megathrust interface, in the first ~24 hours. From these perspectives, we hypothesize that the occurrence of rapidly decaying early afterslip in the mainshock rupture area could be explained by: (1) time-variable conditions that can alter the stress state, stability or frictional properties of patches, enabling patches in the mainshock rupture area to host afterslip, (2) release of the residual slip deficit during the short time frame after the mainshock, and/or (3) deformation within the volume surrounding the megathrust interface due to stress changes induced by the mainshock.

Spatial distribution of afterslip controlled by features on the incoming plate

The southern limit of updip afterslip (area A) overlaps with a highly coupled region (Figure 8c), and is marked by a high density of interseismic (see [START_REF] Font | Seismicity patterns along the Ecuadorian subduction zone: new constraints from earthquake location in a 3-D a priori velocity model[END_REF]) as well as aftershock seismicity (Figure 8b). The location of this seismicity cluster coincides with the northern flank of the Carnegie Ridge. Flanks of ridges have been proposed to be rough [START_REF] Bassett | Gravity anomalies, crustal structure, and seismicity at subduction zones: 1. Seafloor roughness and subducting relief[END_REF]; seismicity here may hence reflect the roughness of the subducting plate [Agurto-Detzel et al. (submitted)]. Also, south of this area, interseismic slip deficit maps indicate a 'creeping corridor', which may be related to the strike-slip Jama Fault Zone that extends down to the megathrust interface [START_REF] Chlieh | Distribution of discrete seismic asperities and aseismic slip along the Ecuadorian megathrust[END_REF], and that may have influenced the southern extent of afterslip.

In the north, a distinct cluster of aftershocks aligns with the northern extent of afterslip in area B (Figures 8b). The northern extent of afterslip in area B appears to be mostly confined to the region south of this cluster of aftershocks. North of this cluster lies the zone that ruptured during the 1958 earthquake, and where the interseismic slip deficit is higher in regions close to the trench (Figure 8c). Strikingly, less afterslip occurs in the updip region between areas A and B, where a double peaked seamount and associated pervasive fracturing was imaged by [START_REF] Marcaillou | Seamount subduction at the North-Ecuadorian convergent margin: Effects on structures, interseismic coupling and seismogenesis[END_REF] (Figure 8). We hypothesize that properties around the subducting seamount may inhibit the pervasive propagation of afterslip into this region. Collectively, these observations suggest that features on the incoming plate exert a strong control on the along-strike segmentation of

A slow slip event near La Plata Island

Daily time series at GPS stations located near La Plata Island (ISPT, MHLA, SLGO) are indicative of a slow slip event (SSE) in this region, with typical phases of increasing velocities then decreasing velocities of deformation (Figure 8d). This SSE was possibly remotely triggered by static stress changes caused by the mainshock rupture [START_REF] Rolandone | Areas prone to slow slip events impede earthquake rupture propagation and promote afterslip[END_REF]. Our inversion results suggest that this SSE occurred trenchward of the section of the megathrust under La Plata Island (area F, Figure 8b), consistent with the results of [START_REF] Rolandone | Areas prone to slow slip events impede earthquake rupture propagation and promote afterslip[END_REF].

The estimated geodetic moment of the SSE is 1.03E+19 (~M w 6.6). Intriguingly, our 72-hour model shows slip on a single patch north of La Plata Island amounting to a geodetic moment of 8.86E+17 N m (M w 5.93) (Figure 8a) that may suggest onset of the SSE in the early postseismic period, although this patch is poorly resolved. Coincidentally, seismicity is also concentrated in the region north of La Plata Island in the first three days, before then occurring in regions south and trenchwards of the island (Figure S13). If we suppose that this seismicity reflects an underlying SSE driving mechanism, we speculate that the SSE might have nucleated north of La Plata Island and migrated southwards and trenchwards with time. However, our speculation should be substantiated with studies of longer HR-GPS time series to better image its onset and migration.

Conclusions

We estimate the spatial and temporal distribution of afterslip spanning various timescales, ranging from ~ 2.5 minutes to 30 days after the earthquake, using HR-GPS and daily GPS time series. Although the HR-GPS time series are relatively noisy compared to the daily GPS time series, our employed inversion method enables us to invert the 72-hour HR-GPS time series to obtain a detailed description of the spatio-temporal evolution of early afterslip, and explore its spatial relation to longer-term afterslip, the mainshock rupture area, and ensuing aftershocks.

We find that the spatial signature of early afterslip in the first 72 hours within the region updip of the mainshock rupture area is consistent with that in the 30-day postseismic period, indicating that afterslip nucleated primarily updip of and adjacent to two peak coseismic slip asperities, and subsequently continued to grow in amplitude with time. The spatial pattern of afterslip appears to be localized in areas prone to host aseismic slip over several seismic cycles and controlled by features on the incoming plate. Interestingly, our results suggest that early afterslip may have occurred within part of the mainshock rupture area, but which may have decayed rapidly, as little afterslip is imaged here in the 30-day afterslip model. One of our most important findings is that early afterslip (in the first 72 hours) following the Pedernales earthquake represents a significant contribution, ~30 %, to the postseismic geodetic moment over the first 30 days. Furthermore, not accounting for afterslip before the first daily GPS position (in this case, in the 12 hours after the earthquake) strongly biases the postseismic geodetic moment, with ~60% of postseismic geodetic moment missing over the first 72 hours.

Figure 1 .[

 1 Figure 1. (a) Interseismic slip deficit, seismic and aseismic events along the

Figure 2 .

 2 Figure 2. Schematic of HR-GPS and daily GPS time series at one station

Figure 4 .

 4 Figure 4. Synthetic tests showing slip recovery in various areas. (a)-(e): Left: input

Figure 6 .Figure

 6 Figure 6. (a)-(f): Spatio-temporal evolution of cumulative afterslip in the first 72

  Table 1 show that in each of the six main afterslip areas, the geodetic moment of afterslip in the 72-hour model is greater than that in the 2-day model. The 72-hour model suggests afterslip in two areas that are not visually evident in the 2-day model: north of the mainshock rupture area near Esmeraldas and north of La Plata Island (areas E and F, respectively, Figure

Table 1 . Estimated geodetic moment in peak afterslip areas for each model. Calculations are based on a constant rigidity of 30 GPa.

 1

	A	Well	4.78E+18	1.59E+09	9.25E+18	3.07E+09	2.82E+19	9.33E+09	1.94	33	17
		resolved
	B	Well	5.62E+18	1.84E+09	8.28E+18	2.71E+09	3.14E+19	1.03E+10	1.47	26	18
		resolved
	C	Well	4.44E+18	1.37E+09	1.02E+19	3.15E+09	9.69E+18	2.99E+09	2.29	105	46
		resolved
	D	Less well	4.47E+18	8.77E+08	1.02E+19	2.00E+09	2.18E+19	4.27E+09	2.28	47	21
		resolved
	E	Less well	3.93E+16	6.85E+07	1.67E+18	2.91E+09	7.05E+17	1.23E+09	42.5	237	6
		resolved
	F	Poorly	3.45E+17	8.59E+07	8.86E+17	2.21E+08	1.03E+19	2.58E+09	2.57	9	3
		resolved
	Total	-	3.23E+19	-	5.22E+19	-	1.81E+20	-	1.62	29	18
			(Mw 6.97)		(Mw 7.11)		(Mw 7.47)

We advocate the importance of imaging afterslip in the minutes to hours following earthquakes, in order to better understand its contribution towards postseismic slip budgets on faults.

Acknowledgements

We thank Hugo Perfettini for insightful, fruitful discussions and advice with implementing PCAIM. We thank Francisco Ortega-Culaciati and Steven Tsang for insightful discussions.

Figures were produced with Generic Mapping Tools (GMT) [Wessel and Smith, 1998]. The acquisition of the data was carried out by IG-EPN and IRD in the frame of the Joint International Laboratory 'Seismes & Volcans dans les Andes du Nord'. We thank all the field operators who participated in the acquisition of this data. This research was supported by the Agence Nationale de la Recherche (ANR) project E-POST (M. Vergnolle, contract number ANR-14-CE03-0002-01JCJC 'E-POST'). L. Tsang's postdoctoral fellowship is funded by the ANR JCJC E-POST, and P. Jarrin's doctoral fellowship is funded by Senescyt -IFTH -Ecuador.

afterslip along this section of the Ecuadorian megathrust, in agreement with the hypothesis 470 proposed by Agurto-Detzel et al. (submitted). 471


Amin Golpour

Coralie Broquard

Sylvain Milla

Hadiseh Dadras

Abdul Rasheed

Taiju Saito

Martin Psenicka

Golpour

Dadras

Rasheed Baloch

Saito

Determination of annual reproductive cycle in male sterlet, Acipenser ruthenus

Keywords: Breeding season, Reproductive cycle, Spermatogenesis, Sterlet, Testicular development

Introduction

In studies on annual gonadal cycle of fish, it is crucial to examine gonadal cycle and spawning pattern [START_REF] Tan-Fermin | Ovarian development and serum steroid hormone profiles in hatchery-bred female catfish Clarias macrocephalus Gunther during an annual reproductive cycle[END_REF][START_REF] Manosroi | Annual sex hormonal profiles, gonad development and age determination of the Mekong giant catfish (Pangasianodon gigas, Chevey)[END_REF]Golpour et al., 2016a[START_REF] Golpour | Subcellular distribution of calcium during spermatogenesis of zebrafish, Danio rerio[END_REF][START_REF] Niksirat | Quantification of egg proteome changes during fertilization in sterlet Acipenser ruthenus[END_REF]. Several studies have already carried out regarding the artificial reproduction of wild sterlet fish [START_REF] Horvath | Successful sterlet, Acipenser ruthenus L., propagation with synthetic LH-RH hormone[END_REF][START_REF] Williot | Artificial spawning in cultured sterlet sturgeon, Acipenser ruthenus L., with special emphasis on hermaphrodites[END_REF]; however, knowledge on the reproductive cycle either in wild population or cultured condition is limited. Only, preliminary results dealing with gonadal development in cultured sterlet have been published previously [START_REF] Jackson | Anatomical, hormonal and histological descriptions of captive Russian sturgeon (Acipenser gueldenstaedtii) with intersex gonads[END_REF]. Development of staging designation in sturgeons is difficult due to the absence of morphological characteristics before sexual maturity, the gonad maturation cycle and the anatomical characteristics of the reproductive system (Lagler et al., 1977;[START_REF] Dettlaff | Sturgeon Fishes: Developmental Biology and Aquaculture[END_REF][START_REF] Doroshov | Observations on the reproductive cycle of cultured white sturgeon, Acipenser transmontanus[END_REF]. The gonadal cycles determination have been done in sturgeons including male and female of white sturgeon, Acipenser transmontanus [START_REF] Doroshov | Observations on the reproductive cycle of cultured white sturgeon, Acipenser transmontanus[END_REF], female siberian sturgeon, Acipenser baerii [START_REF] Williot | An analogue of GnRH is effective for induction of ovulation and spermiation in farmed Siberian sturgeon Acipenser baerii Brandt[END_REF], male and female of amur sturgeon, Acipenser schrenckii [START_REF] Qu | The relationships between gonad development and sex steroid levels at different ages in Acipenser schrenckii[END_REF], Sterlet, Acipenser ruthenus (Golpour et al., 2016b) and both sexes of shortnose sturgeon, Acipenser brevirostrum [START_REF] Colombo | Gonadal development and sex-specific demographics of the shovelnose sturgeon in the Middle Mississippi River[END_REF]. However, limited information is available on gonadal cycles and reproductive physiology of Acipenseriformes.

The determination of reproductive stages in sturgeons using diagnostic methods such as quantification of plasma sex hormones and histology have been widely used [START_REF] Feist | Early identification of sex in cultured white sturgeon, Acipenser transmontanus, using plasma steroid levels[END_REF][START_REF] Craig | Sex assignment of lake sturgeon (Acipenser fluvescens) based on plasma sex hormone and vitellogenin levels[END_REF][START_REF] Webb | Potential classification of sex and stage of gonadal maturity of wild white sturgeon using blood plasma indicators[END_REF][START_REF] Barannikova | Serum levels of testosterone, 11-ketotestosterone and oestradiol-17b in three species of sturgeon during gonadal development and final maturation induced by hormonal treatment[END_REF]. However, these techniques are time consuming and harmful for small sturgeons. Additionally, other methods such as determination of stage through an incision in the body wall has also been used [START_REF] Chen | Observations on gonadal development of 7-11 year old Chinese sturgeon Acipenser sinensis cultured in freshwater[END_REF]. The negative effects of these approaches on small captive broodstock cannot be ignored. Therefore, a reliable method for determination of maturity is critically needed in modern broodstock management to avoid the potential negative effects on gonads or muscles and consequently mortality. The non-invasive technique using ultrasound technology has been mostly tested for gender identification in sturgeons [START_REF] Wildhaber | Gender identification of shovelnose sturgeon using ultrasonic and endoscopic imagery and theapplication of the method to the pallid sturgeon[END_REF][START_REF] Bryan | Estimation of gonad volume, fecundity, and reproductive stage of shovelnose sturgeon using sonography and endoscopy with application to the endangered pallid sturgeon[END_REF][START_REF] Hurvitz | Use of endoscopy for sex and ovarian stage determinations in Russian sturgeon (Acipenser gueldenstaedtii) grown in aquaculture[END_REF][START_REF] Chebanov | Ultrasound diagnostics for sturgeon broodstock management[END_REF][START_REF] Divers | Field endoscopy for identifying gender, reproductive stage and gonadal anomalies in free-ranging sturgeon (Scaphirhynchus) from the lower Mississippi River[END_REF]. In addition, this non-invasive method is also useful to determine the maturity stage of small captive sturgeons.

The sterlet, Acipenser ruthenus, is a relatively small species of sturgeon which is widely distributed in central-eastern Europe and Asia [START_REF] Bemis | Sturgeon rivers: An introduction to acipenseriform biogeography and life history[END_REF]. Decline of natural spawning activity in its population from native habitat is mainly due to loss of spawning grounds, overfishing, poaching, and environmental pollution [START_REF] Jackson | Anatomical, hormonal and histological descriptions of captive Russian sturgeon (Acipenser gueldenstaedtii) with intersex gonads[END_REF]. Various studies have already been conducted regarding the artificial reproduction of wild sterlet fish [START_REF] Horvath | Successful sterlet, Acipenser ruthenus L., propagation with synthetic LH-RH hormone[END_REF][START_REF] Williot | Artificial spawning in cultured sterlet sturgeon, Acipenser ruthenus L., with special emphasis on hermaphrodites[END_REF]; however, knowledge on the reproductive cycle either in wild population or cultured condition is limited. The ultrasound method has been successfully used sex determination of various sturgeon species such as stellate sturgeon Acipenser stellatus [START_REF] Moghim | Determination of sex and maturity in Acipenser stellatus by using ultrasonography[END_REF], shovelnose sturgeon Scaphirhynchus platorynchus [START_REF] Colombo | Use of ultrasound imaging to determine sex of shovelnose sturgeon[END_REF]) and pallid sturgeon Scaphirhynchus albus [START_REF] Wildhaber | Gender identification of shovelnose sturgeon using ultrasonic and endoscopic imagery and theapplication of the method to the pallid sturgeon[END_REF]. Considering above statement, this tool solely is not exclusively used for sex determination and can be useful to ascertain the reproduction status throughout the year, particularly outside of the spawning season [START_REF] Blythe | Determination of sex and maturational status of striped bass (Morone saxatilis) using ultrasonic imaging[END_REF].

Thus, obtaining data on gonad maturity stages using the histology and non-invasive method provide reliable information for artificial culture and restocking programs. The objectives of the study was to evaluate the usefulness of ultrasonic method to assess stages of maturity of the male's gonads in sterlet throughout the year and validate/compare the results by histological analysis.

Materials and methods

Fish

Thirty six adult males of sterlet (~5 years old, weight 0.76-1.24 kg, length 52.5-64.6 cm) were provided by the Fischzucht Rhönforelle GmbH & Co. KG, Marjoss, Germany. Fish were kept in outdoor plastic tanks (4 m3) with a constant flow rate of 20 L/min pond water at hatchery of the Research Institute of Fish Culture and Hydrobiology (Vodnany, Czech Republic).

Sampling and analysis

Gonad samples were collected from three individuals every month from July 2014 to June 2015. Gonadal samples were kept in sampling bottles and fixed in Bouin's solution. After 24 hr, the samples were transferred into 70% ethanol and processed using the routine histological methods. Samples of gonad tissues were processed routinely by dehydration in a graded series of ethanol solutions of increasing concentration, followed by clearing in a series of xylenes solvent (three times), embedding in paraffin, and sectioning at 5 µm. Prepared slides were stained with hematoxylin and eosin (H&E), and then observed under light microscope (Leica DM 750, Heerbrugg, Switzerland). The percentage of spermatogenic cells at different phases of the reproductive cycle was quantified using light microscope (Leica DM 750). The developmental stages were identified based on descriptions given by [START_REF] Doroshov | Observations on the reproductive cycle of cultured white sturgeon, Acipenser transmontanus[END_REF]. In order to perform a quantitative assessment of spermatogenic cells, at least fifteen histology photos were randomly evaluated to count number of spermatogonia (SG), spermatocytes (SC), spermatids (ST) and spermatozoa (SZ) at each phase of the reproductive cycle were calculated in specific number of seminiferous tubules of gonad.

Ultrasound images were collected in the laboratory using a clinical portable ultrasound unit (Chison Eco 1 with a 7.5-MHz linear probe, China). Prior to scanning the body cavity of fish, the head of transducer was first covered with ultrasound gel and then, a plastic sheath was placed over the ultrasound transducer to be protected against rough contact with scutes. A right or left lateral side view of the fish between the pectoral and anal fins was obtained [START_REF] Chebanov | Ultrasound diagnostics for sturgeon broodstock management[END_REF]. During the scanning, constant focus, brightness and settings were used and the gained settings were fixed.

Data analyses

Computer statistical analysis was carried out using statistica software for Windows, v. 9.1 (StatSoft, Inc., Talsa, OK, USA). Data were expressed as mean ± SD. One-way analysis of variance (ANOVA) was followed by a Tukey's test. Statistical significance was tested at the p < .05 level.

Results

A total of 36 sterlets were examined from July 2014 to June 2015 for annual reproductive cycle assessment using histological descriptions and ultrasonography approach. Using histology, development and growth of testes during annual cycle were divided into four stages according to the different seasons. During testicular development, a variable number of spermatogenic cysts at different stages were found in each seminiferous lobule. The seasonal changes of the testes, weight, length, GSI and water temperature are presented in Table 1.

The ultrasonographic images of testes in fall season (Resting phase) showed crenate and narrow with a string-like appeared testis. In addition, the gonad was contained hypoechoic fat (Figure 1a). Ultrasonographic images of sterlet testes in winter (pre-spawning phase) appeared as a smooth grey region (Figure 1b), and the images in spring (Spawning phase) were homogeneous and smooth with distinct margin and well-defined actual shape of testis (Figure 1c). The testes in summer (post-spawning phase) were found to be with smooth dark grey structures and also without any actual shape in ultrasonography images (Figure 1d).

In resting phase, ultrasonographic images revealed that testes were small and surrounded by fat tissue. Gonadosomatic index was lowest throughout the year except in pre-spawning phase (Table 1). Also, this phase was dominated by type A spermatogonia (93.6 ± 0.76%) and few type B spermatogonia (5.2 ± 0.18%) (Figure 2a; Table 2). At the end of fall season, spermatogenesis started with proliferation of spermatogonia (Figure 2b).

In the next phase (Pre-spawning phase), the cyst was predominately occupied by spermatocytes and spermatids (Figure 3), which on average constituted 41 ± 0.63% and 34.8 ± 0.50% (Table 2), respectively. Furthermore, rapid increase of GSI was observed in pre-spawning phase of testis development compared to the other phases (Table 1). In spawning phase, GSI index followed a trend similar to previous months of pre-spawning phase (Table 1). In this phase, the bulged out testes in abdominal cavity appeared in white.

Due to rupture of spermatogenic cyst, spermatozoa uniformly filled the lumen of seminiferous lobules. In spawning phase, spermatids and spermatozoa were dominant (Figure 4), which on average constituted 5.6 ± 0.20% and 92.2 ± 62% (Table 2), respectively. In post-spawning phase, the testes mostly were empty. Further, the leftover meiotic cells and spermatozoa were reabsorbed (Figure 5 a and b). The spermatogonia populations were found to be abundant in cyst of seminiferous tubules (Figure 5b; Table 2). Gonadosomatic index was almost similar to resting and spawning phases (Table 1). The percentage of spermatogonia and spermatozoa were 87.2 ± 0.57% and 11.1 ± 0.43%, respectively (Table 2). In addition, 3 out of 36 males were found to have intersex gonads. The intersex gonads were typical testes with distribution of few pieces of ovary throughout each testis. All observed oocytes in intersex gonads were at the pre-vitellogenic stage (Figure 6 a,b).

Discussion

The present study revealed a seasonal pattern of gonadal development in male sterlet, with onset of testicular meiosis in fall season, and final gonadal maturation during spring. In the present study, the different seasonal water temperatures in holding tanks were associated with different patterns of gonadal development, emphasizing the importance of temperature as a major factor affecting testicular development. Using histology approach, a marked difference in the pattern of gonadal recrudescence in male sterlet was found. Spermatogenic activity was begun at the end of autumn and continued in winter, leading to the appearance of spermatozoa in lobule of lumen in early spring. In sampled testes from January onward, different spermatogenic cells were appeared and developed spermatozoa were observed frequently. We identified four stages of gonad development coinciding with those stages identified in different sturgeon species [START_REF] Doroshov | Observations on the reproductive cycle of cultured white sturgeon, Acipenser transmontanus[END_REF][START_REF] Flynn | Sex differentiation and aspects of gametogenesis in shortnose sturgeon Acipenser brevirostrum Lesueur[END_REF][START_REF] Matsche | Reproductive demographics, intersex, and altered hormone levels in shortnose sturgeon, Acipenser brevirostrum, from Delaware River, USA[END_REF][START_REF] Qu | The relationships between gonad development and sex steroid levels at different ages in Acipenser schrenckii[END_REF]. Further, a similar pattern of gonadal development and maturation of gonad has been observed in green sturgeon Acipenser medirostris (Reviewed by Alavi et al., 2010). It is also noteworthy to note that, all three sterlets with intersex gonads founded in the present study were contained small female gonad portions (testes-ova).

Occurrence of this phenomenon has been also observed in several sturgeons including shortnose sturgeon A. brevirostrum [START_REF] Williot | Artificial spawning in cultured sterlet sturgeon, Acipenser ruthenus L., with special emphasis on hermaphrodites[END_REF][START_REF] Henne | Synchronous hermaphroditism and self-fertilization in a captive shortnose sturgeon[END_REF], Russian sturgeon Acipenser gueldenstaedtii [START_REF] Jackson | Anatomical, hormonal and histological descriptions of captive Russian sturgeon (Acipenser gueldenstaedtii) with intersex gonads[END_REF], and sterlet Acipenser ruthenus [START_REF] Williot | Artificial spawning in cultured sterlet sturgeon, Acipenser ruthenus L., with special emphasis on hermaphrodites[END_REF]. This amount of intersexuality is consistent with the rate of 1% reported in Atlantic sturgeon Acipenser oxyrinchus [START_REF] Van Eenennaam | Reproductive conditions of the Atlantic Sturgeon (Acipenser oxyrinchus) in the Hudson River[END_REF], Russian sturgeon Acipenser gueldenstaedtiiand and Siberian sturgeon Acipenser baerii [START_REF] Rzepkowska | Intersex gonad differentiation in cultured Russian (Acipenser gueldenstaedtii) and Siberian (Acipenser baerii) sturgeon[END_REF].

In the present study, the different seasonal water temperatures in holding tanks were associated with four phases of gonadal development, emphasizing the importance of temperature as a major factor affecting testicular development. Ultrasound imaging has been successfully used to determine gender and maturity of different adult sturgeons including Chinese sturgeon Acipenser sinensis [START_REF] Wildhaber | Gender identification of shovelnose sturgeon using ultrasonic and endoscopic imagery and theapplication of the method to the pallid sturgeon[END_REF][START_REF] Du | Gender and gonadal maturity stage identification of captive Chinese sturgeon, Acipenser sinensis, using ultrasound imagery and sex steroids[END_REF], shovelnose sturgeon Scaphirhynchus platorynchus [START_REF] Colombo | Gonadal development and sex-specific demographics of the shovelnose sturgeon in the Middle Mississippi River[END_REF], beluga sturgeon Huso huso [START_REF] Masoudifard | High validity sex determination of three years old cultured beluga sturgeon (Huso huso) using ultrasonography[END_REF] and sturgeon hybrids Acipenser naccarii × Acipenser baerii [START_REF] Petochi | Sex and reproductive stage identification of sturgeon hybrids (Acipenser naccarii × Acipenser baerii) using different tools: ultrasounds, histology and sex steroids[END_REF] at different reproductive stages. In this context, some descriptions of ultrasonic images have been used for better characterization of fish gonad [START_REF] Jennings | Evidence for heritable preferences for spawning habitat between two walleye populations[END_REF][START_REF] Bryan | Estimation of gonad volume, fecundity, and reproductive stage of shovelnose sturgeon using sonography and endoscopy with application to the endangered pallid sturgeon[END_REF]. The ultrasound imaging also can be used to determine some parameters such as size of the ovary or testis and gonad shape [START_REF] Albers | Gonadosomatic index and fecundity of Lower Missouri and Middle Mississippi River endangered pallid sturgeon estimated using minimally invasive techniques[END_REF]. However, some other parameters including fecundity and GSI might be determined via further analyses [START_REF] Bryan | Estimation of gonad volume, fecundity, and reproductive stage of shovelnose sturgeon using sonography and endoscopy with application to the endangered pallid sturgeon[END_REF]. In some cases, application of ultrasound is not useful to characterize gonad stages, particularly in immature male sturgeons, due to the small gonads. In a study on Chinese Sturgeon, two sets of characterizations (i.e proportion of the gonad to the whole reproductive organs and thickness of the reproductive organs) using ultrasound were applied for sexing and determining sexual maturity stage [START_REF] Du | Gender and gonadal maturity stage identification of captive Chinese sturgeon, Acipenser sinensis, using ultrasound imagery and sex steroids[END_REF]. However, to date, no study has noted that ultrasound can be used determine the stage of maturity based on size of testicular tissue [START_REF] Webb | A practical guide for assigning sex and stage of maturity in sturgeons and paddlefish[END_REF]. Similar to our methodology and results, the density of testicular tissue observed on a sonogram as shades of white is being examined as an indicator of maturity stage [START_REF] Webb | A practical guide for assigning sex and stage of maturity in sturgeons and paddlefish[END_REF]. Considering other parameters of ultrasonography except morphological characteristics, however might be beneficial to determine errors associated with development of less invasive techniques.

The present information suggests a valuable guideline for seasonal gonad differentiation and sexual maturity in male sterlet. Our findings provide useful data for commercial purposes of propagation and comparison of reproductive development among sturgeon species.

Fig. 1 .

 1 Fig. 1. Ultrasound images taken from testes of sterlet during annual cycle: (a) undeveloped testis in autumn (resting phase), (b) testis in pre-spawning state in winter, (c) fully mature testis in spring and (d) testis in post-spawning state. F; fat tissue, M; muscle, T; testis.

Fig. 2 .

 2 Fig. 2. Transverse section of the sterlet testis during autumn (resting phase) shows the onset of testicular meiosis. a; arrows show spermatogonia type A enclosed in cysts (scale bar 50 µm) b; arrows show proliferation of spermatogonia type A (scale bar 100 µm). H&E staining was used for histology sections.

Fig. 3 .

 3 Fig. 3. Transverse section of the sterlet testis during winter (pre-spawning phase) exhibiting ongoing meiosis. a; early-spermatogenesis (meiosis); arrows show developing spermatocyte and spermatid in the clusters, b; late meiosis and spermiogenesis exhibit the cysts with developing spermatogenic cells and spermatozoa. H&E staining was used for histology sections. SG -spermatogonia, SCspermatocyte, ST -spermatid, SZ -spermatozoa. All scale bars are 100 µm.

Fig. 4 .

 4 Fig. 4. Transverse section of the sterlet testis during spring (spawning phase) shows cysts with differentiated spermatozoa. H&E staining was used for histology sections. SZ; spermatozoa. Scale bar is 100 µm.

Fig

  Fig. 5. Transver se section of the sterlet testis during summer (post-spawning phase) showing the absorption of unshed spermatozoa. a; spermatogonia type A located along the periphery of the lobule (scale bar 50 µm), b; Higher magnification shows the reabsorbing unshed spermatozoa and remained meiotic cells (scale bar 100 µm). H&E staining was used for histology sections. CS, Cyst; SGA, Spermatogonia type A; SGB, Spermatogonia type B; SZ, spermatozoa.

Fig. 6 .

 6 Fig. 6. Transverse section of sterlet intersex gonad (testis-ova). a; showing a pre-vitellogenic oocyte close to testicular lobules containing spermatozoa at the interface between ovary and testis (scale bar 50 µm), b; Higher magnification showing the close proximity between oocyte and spermatozoa (scale bar 200 µm). H&E staining was used for histology sections. OC: Oocyte, SZ: Spermatozoa.

Table 1

 1 Monthly growth parameters and gonadosomatic indices of sterlet males (n=3). Data are mean ± SD.

	Years	Developmental Stages	Months	Weight	Length	Gonadosomatic index (%)	Temperature
	2014	Post-spawning phase	July	1016 ± 27.6	64.3 ± 4.8	1.44 ± 0.1	24 ° C
		Post-spawning phase	August	1023 ± 28.4	64.6 ± 6.2	1.28 ± 0.2	22 ° C
		Post-spawning phase	September	1023 ± 25.8	64.3 ± 5.2	1.16 ± 0.3	20 ° C

Table 2

 2 Monthly percentage of spermatogonia (SG), spermatocytes (SC), spermatids (ST) and spermatozoa (SZ) during different phases of the reproductive cycle in sterlet.

	Resting phase	October		1040 ± 19.5	58.2± 7.3	1.3 ± 0.1	15 ° C
	Resting phase	November	1124.6 ± 16.3	60.3± 3.9	1.02 ± 0.02	11 ° C
	Resting phase Years Developmental stages December Month	760 ± 17.5 SG (%)	61 ± 4.6 SC (%)	1.01 ± 0.06 ST (%)	11 ° C SZ (%)
	2015 2014 Post-spawning phase Pre-spawning phase Post-spawning phase Pre-spawning phase Post-spawning phase	January February	1166.6 ± 20.8 86.5 ± 2.7 58.9 ± 5.7 1.2 ± 0.1 August July 87 ± 2 1.2 ± 0.2 1066.6 ± 27.5 60.7 ± 3.4 September 88.2 ± 2 1.5 ± 0.2	1.83 ± 0.1 0.3 ± 0.1 0.28 ± 0.01 2.04 ± 0.2 0. 18 ± 0.01	12.3 ± 0.05 10 ° C 11.2 ± 0.03 12 ° C 9.8 ± 0.07
	Pre-spawning phase Resting phase	March	850 ± 33.6 October 94.7 ± 1.5 53.3 ± 3.8 4.3 ± 0.1	1.95 ± 0.7 1.3 ± 0.1	13 ° C 0.59 ± 0.02
	Spawning phase Resting phase	April	900 ± 23.5 November 94.5± 1.4	52.5 ± 6.1 4.8 ± 0.3	1.78 ± 0.3 0.57 ± 0.01	15 ° C 0.67 ± 0.01
	Spawning Resting phase phase 2015 Pre-spawning phase Spawning phase Pre-spawning phase	May June	863.3 ± 25.9 December 92.1 ± 1 January 7.1 ± 0.09 54.8 ± 6.5 6.5 ± 0.5 42.3 ± 0.15 35.6 ± 0.2 0.8 ± 0.1 1.01 ± 0.06 800 ± 29.9 55.2 ± 5.3 0.79 ± 0.2 February 6.5 ± 0.07 41 ± 1 34.4 ± 0.6	0.72 ± 0.02 16 ° C 16.6 ± 0.4 19 ° C 18.2 ± 0.2

Acknowledgments

The study was financially supported by the Ministry of Education, Youth and Sports of the Czech Republic project CENAKVA (LM2018099) and Biodiversity (CZ.02.1.01/0.0/0.0/16_025/0007370) and the Czech Science Foundation (grant number 17-19714Y), by the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No. 642893 (IMPRESS).