Luca Dal

Romain Jolivet

Ylona Van Dinther

Segmentation of the Main Himalayan Thrust illuminated by Bayesian inference of interseismic coupling

We propose a probabilistic interseismic coupling model of the Main Himalayan Thrust

• Our Bayesian analysis suggests heterogeneous coupling along the Main Himalayan Thrust

• Spatial coincidence between low coupling regions and subsurface ridges suggests a possible segmentation of the Main Himalayan Thrust

Introduction

The magnitude (M w ) 7.8 Gorkha-Nepal earthquake of 2015 was a dramatic reminder of the importance of understanding seismic hazard in the Himalayas (e.g., Bilham et al., 2001) (Fig. 1). This densely populated region has been struck by some of the largest continental earthquakes, including the 1905 M w 7.8 Kangra earthquake, and the 2005 M w 7.6 Kashmir earthquake. It is also known to have produced larger historical earthquakes of magnitude larger than 8, such as the 1950 M w 8.5 Assam earthquake and the 1934 M w 8.2 Nepal-Bihar earthquake [START_REF] Chen | Source parameters of earthquakes and intraplate deformation beneath the shillong plateau and the northern indoburman ranges[END_REF][START_REF] Sapkota | Primary surface ruptures of the great himalayan earthquakes in 1934 and 1255[END_REF].

Tectonic shortening imposed by plate motion on interseismically locked, seismogenic faults leads to elastic strain accumulation, which is ultimately released by earthquakes (e.g., [START_REF] Savage | A dislocation model of strain accumulation and release at a subduction zone[END_REF]. Geodetic observations allow derivation of spatial variations of strain rates induced by interseismic locking of the main boundary faults.

Interseismic locking varies spatially along plate boundary faults, highlighting segments that are locked versus those that slip aseismically (e.g., [START_REF] Avouac | From geodetic imaging of seismic and aseismic fault slip to dynamic modeling of the seismic cycle[END_REF]. In the Himalayas, most large earthquakes have ruptured the Main Himalayan Thrust (MHT), the main plate boundary fault between the Indian and Eurasian plates (e.g., Bilham et al., 2001). Interseismic strain results from relative plate motion across the MHT, which seems to increase eastwards from ∼13 ±1.7 mm yr -1 to ∼21 ±2.0 mm yr -1 across the Himalayas [START_REF] Stevens | Interseismic coupling on the main himalayan thrust[END_REF].

There is now a relatively good coverage of Global Positioning System (GPS) stations [START_REF] Kreemer | A geodetic plate motion and global strain rate model[END_REF], which allows quantification of the degree of fault locking and its spatial distribution on the MHT. For the purposes of this study, fault locking is referred to as coupling, defined as the ratio of the deficit of slip rate in the interseismic period divided by the long-term slip rate imposed by tectonic shortening. This ratio varies from 0 for a portion of the MHT that slips aseismically at a rate equal to the plate convergence rate to 1 for a fully locked fault. Many models derived from geodetic data suggest that interseismic coupling on subduction megathrusts worldwide is spatially heterogeneous, defining locked patches where stress increases more rapidly surrounded by segments that slip aseismically (e.g., [START_REF] Chlieh | Heterogeneous coupling of the sumatran megathrust constrained by geodetic and paleogeodetic measurements[END_REF][START_REF] Moreno | Heterogeneous plate locking in the south-central chile subduction zone: Building up the next great earthquake[END_REF]. Megathrust earthquakes occur to first order on these interseismically locked patches (e.g., [START_REF] Loveless | Spatial correlation of interseismic coupling and coseismic rupture extent of the 2011 mw= 9.0 tohoku-oki earthquake[END_REF]. Conversely, previously published models in the Himalayan region suggest that the MHT is nearly uniformly locked from the surface to beneath the front of the high range over a width of 100 ±20 km [START_REF] Ader | Convergence rate across the nepal himalaya and interseismic coupling on the main himalayan thrust: Implications for seismic hazard[END_REF][START_REF] Li | Geodetic imaging mega-thrust coupling beneath the himalaya[END_REF][START_REF] Ponraj | Estimation of strain distribution using gps measurements in the kumaun region of lesser himalaya[END_REF][START_REF] Sreejith | Audit of stored strain energy and extent of future earthquake rupture in central himalaya[END_REF][START_REF] Stevens | Interseismic coupling on the main himalayan thrust[END_REF][START_REF] Yadav | Strong seismic coupling underneath garhwal-kumaun region, nw himalaya, india[END_REF]. Therefore, the MHT appears to be an anomaly with respect to subduction megathrusts worldwide. However, a recent and detailed GPS velocity solution over Bhutan highlights local variations in coupling, with a less wide fully coupled segment in eastern Bhutan and a partially unlocked deeper crustal ramp [START_REF] Marechal | Evidence of interseismic coupling variations along the bhutan himalayan arc from new gps data[END_REF].

Mapping the spatial distribution of interseismic coupling is key to understanding the segmentation of the MHT and its seismic potential (e.g., [START_REF] Stevens | Millenary mw>9.0 earthquakes required by geodetic strain in the himalaya[END_REF]. In addition, quantitative assessments of uncertainties on the slip deficit inferred from geodetic data are still missing. Here, we use a Bayesian approach to derive a probabilistic estimate of fault coupling along the MHT. Such probabilistic estimate does not rely on any spatial smoothing, an approach often used to regularize slip inversions. Contrarily to more classic, regularized slip and coupling inversion techniques, our approach explores all the possible distributions of coupling along the fault that explains the data. In doing so, we avoid explicit formulation of arbitrary constraints that restrict the range of possible models. Furthermore, we provide full posterior probability distributions of coupling made of the ensemble of plausible models that fit the observations and are consistent with prior constraints.

-3-manuscript accepted for publication at Geophysical Research Letters

Methods

Geodetic Observations

To assess the extent of interseismic coupling, we use an updated set of GPSderived velocities spanning the 1990s to 2015 [START_REF] Kreemer | A geodetic plate motion and global strain rate model[END_REF]. This dataset was recently used to estimate the width of coupling across the MHT through a series of two-dimensional transects [START_REF] Lindsey | Structural control on down-dip locking extent of the himalayan megathrust[END_REF] (Fig. S1). To better resolve the pattern of coupling we also include levelling measurements from the Survey of Nepal collected between 1977[START_REF] Chen | Source parameters of earthquakes and intraplate deformation beneath the shillong plateau and the northern indoburman ranges[END_REF][START_REF] Jackson | Constraints on himalayan deformation inferred from vertical velocity fields in nepal and tibet[END_REF].

Our fault model considers the total length of the arc over roughly 2000 km [START_REF] Stevens | Interseismic coupling on the main himalayan thrust[END_REF], from 73 • E to 96 • E (Fig. S1). The dip angle is set to 10 • as in previous studies [START_REF] Ader | Convergence rate across the nepal himalaya and interseismic coupling on the main himalayan thrust: Implications for seismic hazard[END_REF][START_REF] Stevens | Interseismic coupling on the main himalayan thrust[END_REF]. To account for the arcuate and irregular shape of the MHT, we divide the fault geometry into four planar sub-fault blocks. These sub-faults are discretized in 353 triangles with variable size, smaller near the surface (∼45-50 km) and larger down-dip (∼250 km).

As estimates of the extent of plate coupling depend on the direction and rate of the plate convergence, we set a specific convergence rate and azimuth for each sub-fault based on plate motion models and regional geodetic measurements (see Table S1).

We also assume that the otherwise rigidly moving Indian plate is fragmenting around the eastern end of the Himalayan arc where there is evidence for internal deformation [START_REF] Vernant | Clockwise rotation of the brahmaputra valley relative to india: Tectonic convergence in the eastern himalaya, naga hills, and shillong plateau[END_REF]. To account for this deformation, we consider two distinct blocks and related subnetworks: Indian and the Shillong block south of the eastern Himalaya (Fig. S1). The elastostatic Green's functions of each sub-fault and block, relating unit slip on each patch to surface displacements at the data locations, are thus computed assuming a different long-term convergence rate and azimuth. The inversion solves for the poles of rotation with respect to India, as well as the coupling pattern. We predict the response to subsurface fault slip at each data (GPS and levelling) location assuming a stratified semi-infinite elastic medium [START_REF] Simons | Coseismic deformation from the 1999 m w 7.1 hector mine, california, earthquake as inferred from insar and gps observations[END_REF][START_REF] Zhu | A note on the dynamic and static displacements from a point source in multilayered media[END_REF]. Our elastic half-space consists of 3 layers with elastic properties following the average properties of the crust in the region [START_REF] Pandey | Interseismic strain accumulation on the himalayan crustal ramp (nepal)[END_REF] (Fig. S2).

Bayesian analysis

We consider the forward problem d = Gm, where d is the vector containing the horizontal component of the GPS-derived velocities and of the vertical displacement of the levelling-derived measurements, m is the vector of model parameters, and G is the matrix of Green's functions. Green's functions consist of surface displacement of each data location (GPS and levelling) for unit slip on each point of the fault.

G includes elements of reference frame estimate for the GPS displacement field as well. We solve the inverse problem to infer the distribution of model parameters (m) consistent with our data (d). The solution to this inverse problem is non-unique and large uncertainties on the parameters are expected. Therefore, instead of deriving a single model, we explore all potential solutions using a Bayesian approach to derive -4-manuscript accepted for publication at Geophysical Research Letters the posterior probability density function (PDF) of our model given our set of data, p(m|d) [START_REF] Minson | Bayesian inversion for finite fault earthquake source models i-theory and algorithm[END_REF]. We thus derive a probabilistic estimate of the degree of fault locking along the MHT. Following Bayes's theorem, we write the posterior PDF as:

p(m|d) ∝ p(m) exp [- 1 2 (d -Gm) T C -1 χ (d -Gm)]; (1)
where p(m) is the prior PDF of the model that describes our state of knowledge before acquiring data, d is the data vector, G the Green's functions matrix, and m is the vector of model parameters. C χ is the covariance matrix in the data space, defined as the sum of the data covariance matrix C d , (i.e., measurement error matrix) and the prediction error matrix, C p , which describes uncertainties in our elastic model for the crust:

C χ = C d + C p . (2)
We account for prediction uncertainties (i.e., in the C p covariance matrix) due to inaccuracies in this layered elastic model using the approach of Duputel et al.

(2014). The uncertainty on the elastic structure, presented as grey standard deviation in Fig. S2, is estimated by comparing previously published models across the entire Himalayan region. We need a careful description of errors in order to not overfit the data and produce reasonable estimates of coupling uncertainties on each portion of the megathrust. Considering the variability and uncertainties on the 1D velocity models (Fig. S2), we test both 5% and 10% uncertainty on the elastic properties of each layer of the elastic half-space. We prescribed 10% uncertainty in the reference model.

Previous studies have shown that not including the prediction errors when sampling for posterior PDF of coupling may lead to biased results [START_REF] Jolivet | Aseismic slip and seismogenic coupling along the central san andreas fault[END_REF].

The solution of the inverse problem is the posterior ensemble of all plausible interseismic coupling models (m) that fit the GPS data (d) and that are consistent with our prior hypotheses. We test two different prior hypothesis: (1) a uniform PDF assuming no prior knowledge on the model parameters and (2) a binary prior (with 10% uncertainties) along-dip -that is, imposing a prior of 1 from the surface to 100 km downdip (18 km depth), followed by a linear decrease between 100 and 110 km (18-20 km depth), and a constant prior of 0 at depth >20 km. We generate 360'000 models corresponding to the posterior information on geodetic coupling given measured interseismic velocities. Our final solution consists of an ensemble of models that are statistically distributed according to the posterior PDF.

Interseismic coupling distribution

The posterior mean coupling model shows remarkable spatial variations in interseismic coupling, both down-dip and along-strike the MHT (Fig. 2a). We infer discrete highly-coupled segments along the fault extending throughout the arc. Coupling is high at the front, where the Main Frontal Thrust (MFT) -the youngest and southernmost thrust of the Himalayan wedge -reaches the surface and the resolu- regions with lower coupling (<0.3). These lateral variations of coupling raise two fundamental questions. First, are these variations a robust feature driven by the data?

Second, what difference between our approach and previously published methods leads to the emergence of these low coupling segments along the MHT?

We interrogate the posterior PDF of coupling to assess the robustness of our estimate by considering six representative local PDFs of coupling (nodes 1-6 in Fig. 2a). Marginal PDFs for nodes 2, 4, and 6 show a >75% probability that fault coupling exceeds 0.8. On the other hand, for node 1 and 5, we infer a 68% probability that fault coupling is lower than 0.4, confirming a significantly higher probability of low coupling.

In contrast, the shape of the cumulative distribution function of coupling on node 3 suggests a wider range of probable values, not far from the initial prior distributions (i.e., uniform probability of coupling) but still tending towards lower values of coupling.

This similarity of prior and posterior distributions suggests the data are less informative to derive fault coupling along this limited section of the model. Although our inversions include larger covariances compared to previous studies [START_REF] Ader | Convergence rate across the nepal himalaya and interseismic coupling on the main himalayan thrust: Implications for seismic hazard[END_REF][START_REF] Stevens | Interseismic coupling on the main himalayan thrust[END_REF], observed vs. predicted horizontal GPS displacements (Fig. 2b) and residuals (Fig. S3) from our mean model show a good agreement with the geodetic measurements. However, it is important to note that the regions with low coupling are spatially coincident with the block boundaries, and where the GPS data are less dense (Fig. S1). To check the reliability of our results, we perform two tests in which we assume (i) a fault model with a constant convergence rate (17±3 mm yr -1 )

and azimuth (8 • E), and (ii) 5% uncertainty on the elastic properties of each layer of the elastic half-space. Despite small variations, these fault models show very similar results as the regions with low coupling are still easily detectable (Fig. S4 andS5).

These results thus provide strong evidence that the inferred variations of interseismic coupling are consistent with the data and the given posterior uncertainties.

For comparison, instead of assuming a uniform non-informative prior, we consider a constrained prior PDF distribution. In particular, we assume that the MHT is highly-coupled from the surface to 18-km-depth and with negligible coupling at depths greater than 20 km. To do so, we consider a tight Gaussian prior PDF with a standard deviation of 10%, centered on 1 from the surface to 18-km-depth (i.e., fully locked) and on 0 for depths greater than 20 km (i.e., slip at plate convergence rate).

All other assumptions and hypothesis are kept identical to our reference model. This

strong assumption is based on previous geodetic studies using GPS campaign and levelling measurements [START_REF] Bollinger | Stress buildup in the himalaya[END_REF], a few continuous GPS stations [START_REF] Bettinelli | Plate motion of india and interseismic strain in the nepal himalaya from gps and doris measurements[END_REF], and a recent physics-based forward model of the Himalayan seismic cycle [START_REF] Dal Zilio | Bimodal seismicity in the himalaya controlled by fault friction and geometry[END_REF]. Also, the location of the down-dip end of the locked fault zone is consistent with the seismicity pattern observed in Nepal Himalaya (Cattin and -6-manuscript accepted for publication at Geophysical Research Letters [START_REF] Cattin | Modeling mountain building and the seismic cycle in the himalaya of nepal[END_REF], and the postseismic deformation following the 2015 M w 7.8 Gorkha earthquake [START_REF] Gualandi | Pre-and post-seismic deformation related to the 2015, mw 7.8 gorkha earthquake, nepal[END_REF].

The posterior mean model suggests the MHT is highly-coupled from its surface expression along the Himalayan foothills to beneath the front of the high range about 100 km to the north (Fig. S6). In general, coupling is nearly binary (1-to-0), with a sharp transition between the coupled and creeping zones, similarly to what has been observed in previous studies [START_REF] Ader | Convergence rate across the nepal himalaya and interseismic coupling on the main himalayan thrust: Implications for seismic hazard[END_REF][START_REF] Stevens | Interseismic coupling on the main himalayan thrust[END_REF]. This distribution of coupling is actually very close to the prior PDF we have considered.

Strikingly, residuals are of the same order of magnitude as in our reference model with uniform prior. This leads to the conclusion that both models, the one assuming a constrained prior and the one with a uniform prior, fit the data almost equally well and are both possible. However, examination of marginal posterior PDFs suggest that the model derived from highly-constrained prior is only a subset of potential models.

As shown by marginal PDFs of coupling on node 1, 3, and 5 (Fig. S6), models with high coupling are not excluded, but they are less probable than low coupled models given the available data. Homogeneous coupling along the MHT thus results mainly from a strong prior assumption. On the other hand, relaxing this hypothesis reveals a complex distribution of coupling and this complex distribution is more probable given the available data as shown by the posterior uncertainties directly translated from the data.

Segmentation of the Himalayan megathrust and seismicity

The spatial variability of interseismic coupling is an important factor to consider in interpreting relationships with the slip distribution of historical ruptures (Fig. 3). Along-strike, low coupling regions around nodes 1, 3, and 5 (Fig. 2a) coincide with three subsurface Indian basement ridges in the Indo-Gangetic plains: The Delhi-Haridwar ridge (DHR), the Faizabad ridge (FR), and the Munger-Saharsa ridge (MSR). These subsurface basement ridges extend northward across the Himalayan foothills and they might have affected the development of the arcuate shape and the morphology of the Himalayan arc. They have been detected first from basin thickness changes [START_REF] Raiverman | Basin geometry, cenozoic sedimentation and hydrocarbon prospects in northwestern himalaya and indo-gangetic plains[END_REF], and they have been confirmed by the analysis of arc-parallel gravity anomalies [START_REF] Hetényi | Segmentation of the himalayas as revealed by arc-parallel gravity anomalies[END_REF]. They typically include oblique to transverse synclines, folds, fracture zones and/or strike-slip faults oriented transverse to the strike of the Himalaya, which are associated with normal faulting and strike-slip seismicity [START_REF] Cannon | Segmented strain accumulation in the high himalaya expressed in river channel steepness[END_REF][START_REF] Murphy | Limit of strain partitioning in the himalaya marked by large earthquakes in western nepal[END_REF]. Details of the origin and nature of the lineaments, however, remain mostly obscure at present.

The subduction of these major tectonic structures may relate to a segmentation in the seismic behavior of the MHT. Despite significant uncertainties, the rupture extent of major and great Himalayan earthquakes for past two hundred years appears to correlate with the regions of ridge subduction and the pattern of coupling (Fig. 3b). For example, the eastern edge of the 1905 Kangra earthquake rupture appears to approximately coincide with the Delhi-Hardwar ridge [START_REF] Gahalaut | Possible influence of subducting ridges on the himalayan arc and on ruptures of great and major himalayan earthquakes[END_REF][START_REF] Wallace | Surface deformation in the region of the 1905 kangra mw= 7.8 earthquake in the period 1846-2001[END_REF]. The western edge of the 1803 rupture, though not well constrained, abuts the Delhi-Hardwar ridge [START_REF] Rajendran | The status of central seismic gap: a perspective based on the spatial and temporal aspects of the large himalayan earthquakes[END_REF]. The 1934 Nepal Bihar earthquake rupture propagated eastward until or into the Munger-Saharsa ridge [START_REF] Hough | Site response of the ganges basin inferred from re-evaluated macroseismic observations from the 1897 shillong, 1905 kangra, and 1934 nepal earthquakes[END_REF]. However, it should be noted that an error of about 50 km could exist in the estimation of the rupture extent and rupture dimensions, and in the northward extrapolation of the ridges under the Himalayan region [START_REF] Gahalaut | Possible influence of subducting ridges on the himalayan arc and on ruptures of great and major himalayan earthquakes[END_REF].

Recorded M≥4.5 earthquakes (NEIC catalog from 1964 to 2018) show that the background crustal seismicity varies along-strike (Fig. 3c). Cumulative seismic moment released reveals along-strike variations with less seismic moment released in the three regions corresponding to the subsurface ridges. In particular, the cumulative effect of all smaller events does not reach that of a single M6 event. Furthermore, low seismicity across the Faizabad ridge has been evidenced by a dense seismic network in Nepal [START_REF] Ader | Convergence rate across the nepal himalaya and interseismic coupling on the main himalayan thrust: Implications for seismic hazard[END_REF]. Relocated microseismic activity [START_REF] Ader | Convergence rate across the nepal himalaya and interseismic coupling on the main himalayan thrust: Implications for seismic hazard[END_REF][START_REF] Mahesh | One-dimensional reference velocity model and precise locations of earthquake hypocenters in the kumaon-garhwal himalaya[END_REF], which extends over two regions from 77 • E and 81 • E and from 81 • E to 87 • E, respectively, follow the down-dip edge of the highly-coupled patches quite well (Fig. S7). East of Bhutan, in Arunachal Pradesh, there is a seismically active band that extends into the orogen. Also, the cluster of seismicity extending into the Himalayan orogen at 77.5-78 • E turns northward as a possible continuation of the Delhi-Haridwar ridge.

Coupling vs. topography

If interseismic coupling is stable over geological time scales (Kyr), the resulting pattern of strain should influence the construction of topography locally. To first order, topography contour lines show how the Himalaya rise abruptly at a distance of about 120-150 km from the Main Frontal Thrust (Fig. 4a). Notably, the down-dip end of the highly-coupled section of the MHT is spatially coincident with the abrupt rise in topography. Such correlation could reflect anelastic strain, which contributes to topographic uplift in the long run [START_REF] Bilham | Gps measurements of present-day convergence across the nepal himalaya[END_REF][START_REF] Stevens | Interseismic coupling on the main himalayan thrust[END_REF].

Such anelastic deformation suggests that the down-dip limit of the brittle-ductile transition is located where the shear stress accumulation is the highest, and drops under topography greater than 3500 m of elevation [START_REF] Bollinger | Stress buildup in the himalaya[END_REF].

It is also interesting to notice that the largest across-strike graben of the Himalayan belt, the Kaurik-Chango Rifts, the Thakola Rifts and the Yadong Rifts, lie near the prolongation of the three ridges (Fig. 4a). In particular, normal faulting earthquakes continuously occur along N-S fault planes in the Kaurik-Chango Rift [START_REF] Gahalaut | Segmentation of seismicity along the himalayan arc due to structural heterogeneities in the under-thrusting indian plate and overriding himalayan wedge[END_REF], where the eastern limit of the 1905 Kangra earthquake is located. Possibly, the Kaurik-Chango Rift keeps this part of the Himalayan seismic belt in a dilated (or low stress) state so that earthquake ruptures along the MHT do not extend through it. It should be noted, however, that the assumption of a linear extension in the trend of these ridges below the Himalaya is uncertain and that the ridges may in fact deviate from their extrapolated linear trend.

A link between topography and coupling suggests a long-term influence of the seismic behavior of the megathrust and topography building. One possibility to reconcile the segmentation of topography and that of coupling is to consider the topographic effects of ridge subduction. In subduction zones, seafloor observations [START_REF] Singh | Aseismic zone and -13-manuscript accepted for publication at Geophysical Research Letters earthquake segmentation associated with a deep subducted seamount in sumatra[END_REF] and numerical experiments [START_REF] Ruh | Crustal deformation dynamics and stress evolution during seamount subduction: High-resolution 3-d numerical modeling[END_REF] suggest that subduction of large topographic structures (e.g., seamounts) result in a bulge above the underthrusting seamount and normal faults within the overriding plate. If this holds for the Himalayan orogen, the long-term passage of these three basement ridges could explain the topography cusps and the arc-perpendicular normal faults. More generally, the correlation between the pattern of interseismic coupling and the topography suggests that segmentation and effective frictional properties vary little in time, while irrecoverable strain of the Himalayan wedge is small during the interseismic period. As previously inferred from the colocation of the bulge of interseismic uplift with the front of the Higher Himalaya [START_REF] Meade | The signature of an unbalanced earthquake cycle in himalayan topography?[END_REF], some fraction of the interseismic geodetic strain, albeit small as discussed above, could be actually anelastic and contributes to topographic uplift over geological time scales.

To investigate lateral variations in collisional structure we also compare (i) the Arc-Parallel Topography Anomalies (APTA) [START_REF] Hetényi | Segmentation of the himalayas as revealed by arc-parallel gravity anomalies[END_REF], which define the topography differences from the average arc-perpendicular profile, (ii) the 0. This latter is certainly a major boundary in (or of) the Indian plate [START_REF] Hetényi | Segmentation of the himalayas as revealed by arc-parallel gravity anomalies[END_REF]. Despite some exceptions (e.g., between 89 and 92 • E), the correlation between the topography anomalies and coupling suggests that fault segments beneath positive APTA are likely more aseismic, as already observed in subduction zones [START_REF] Song | Large trench-parallel gravity variations predict seismogenic behavior in subduction zones[END_REF].

Conclusions

In summary, we propose that inherited structures along the Himalayan arc induce a zonation of interseismic coupling along the plate interface. Regardless of uncertainties, our probabilistic estimate of interseismic coupling along the MHT suggests that the Himalayan megathrust is segmented by regions of low coupling. These low-coupled patches overlap with inherited structures. Despite competing hypotheses, the segmentation of the orogen into these blocks is also visible in the along-strike rupture extension of historical earthquakes. Large earthquakes of the past millennium do not propagate across segment boundaries defined by lower plate inherited structure [START_REF] Gahalaut | Segmentation of seismicity along the himalayan arc due to structural heterogeneities in the under-thrusting indian plate and overriding himalayan wedge[END_REF][START_REF] Hetényi | Segmentation of the himalayas as revealed by arc-parallel gravity anomalies[END_REF]. Therefore, the here identified low-coupled patches potentially set three boundaries for the lateral extent of future megathrust earthquakes. From a general agreement between previous results [START_REF] Ader | Convergence rate across the nepal himalaya and interseismic coupling on the main himalayan thrust: Implications for seismic hazard[END_REF][START_REF] Stevens | Interseismic coupling on the main himalayan thrust[END_REF], geodetic data, and geological constraints, our results confirm that the fraction of anelastic deformation in the highly-coupled regions might be within uncertainties of the geodetic data (≤10%). These new results highlight important targets for future studies, which will be essential to question how these segment boundaries control interseismic strain accumulation. Future GPS campaigns should assist in answering this challenging but crucial question.

  5 contour line of coupling, and (iii) the location of the three subsurface ridges (Fig 4b). The variation of the APTA exceed 500 m with a number of alternating negative and positive patches throughout the arc. Neglecting important effects of climate conditions and surface processes (e.g., Whipple, 2009), the higher values at the front of the range are located where the coupling is low and the three subsurface ridges come across the Main Frontal Thrust. The western boundary near 76.5 • E, aligned with the Delhi-Haridwar ridge (DHR), is visible in the topography anomalies. The middle boundary near 82.5 • E, aligned with the Faizabad ridge (FR), as well as the eastern end of the Munger-Saharsa ridge (MSR) near 87.5 • E, are also visible in the topography anomalies.

Figure 1 .

 1 Figure 1. Great earthquakes in the Himalayas since 1500 A.D. Yellow rupture areas of magnitude >7.5 and larger events are schematic and represent the published along-arc extent estimates Berthet et al. (2014); Bilham (2009); Bilham and England (2001); Bollinger et al. (2014); Kumar et al. (2010). Background map shows the population density distribution Ciesin (2005). Inset locates main map.

  Geophysical Research Letterstion resulting from observations is highest. Down-dip, coupling tapers to negligible values around 100-120 km north from the MFT. From west to east, we observe a large lateral variability in geodetic coupling. Four large (ca. 300-450 km) and discrete highly-coupled segments (>0.8) are separated by three relatively smaller (ca. 100 km)

	manuscript accepted for publication at
	-5-

Acknowledgments

Work for this article has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 Research and Innovation Programme (Grant Agreement Number: 758210 -Geo4D). L.D.Z. was supported by the SNSF 2-77090-14 project Swiss-AlpArray SINERGIA and the Cecil and Sally Drinkward Postdoctoral Scholar at Caltech. We thank V.L. Stevens, J.-P. Avouac, G. Hetény, M.

Simons, L. Bollinger, C. Pranger, A. Gualandi, K. Wang, and S. Michel for constructive comments and discussions. We are grateful to G. Hetény for providing us with GMT scripts to plot Fig. 1 and the Arc-Parallel Topography Anomalies (Hetényi et al., 2016). We thank the Editor and two anonymous reviewers for providing insightful comments that helped to improve the quality of this paper. Data related to this paper can be downloaded from the following link: https://doi.org/10.22002/d1.1317

Author contributions

L.D.Z. and R.J. conceived the research, designed the study, performed and analyzed the fault modeling, and wrote the manuscript. Y.v.D. commented on the analysis. All authors took part in finalizing the manuscript.


PharmD Thuy Nguyen
email: thithu-thuy.nguyen@aphp.fr

Constance Delaugerre

Marc-Antoine Valantin

Corinne Amiel

Emmanuelle Netzer

Thomas L'yavanc

Michel Ohayon

Nadia Valin

Nesrine Day

Georges Kreplak

Gilles Pialoux

Vincent Calvez

Jean-Michel Molina

Anne-Geneviève Marcelin

Eve Todesco

Shared HCV Transmission Networks among HIV-1 Positive and Negative Men Having Sex with Men by Ultra-Deep Sequencing Running head: Shared HCV Transmission Networks

Keywords: shared HCV transmission chains, recent HCV infection, men having sex with men, ultra-deep sequencing

Objective: Several studies reported HCV transmission networks among men having sex with men (MSM) in Europe and the spread of HCV strains from HIV-HCV co-infected toward HCV mono-infected MSM. We aimed to investigate HCV transmission dynamics among HIV-positive and HIV-negative MSM by ultra-deep sequencing (UDS).

Design and Methods: NS5B fragment (388 bp) was sequenced from virus of 50 HIVpositive and 18 HIV-negative patients diagnosed with recent HCV infection. UDS data were analysed by Geneious (version 10.3.2). Phylogenetic trees were constructed by FastTree (version 2.1) and submitted to ClusterPicker (version 1.2.3) for transmission chain detection at different thresholds of maximum genetic distance (MGD) (3% for Sanger, 3% and 4.5% for UDS).

Results: Ten, seven-teen, and eight-teen HCV transmission chains were identified by Sanger at 3%, UDS at 3% and at 4.5% of MGD, respectively. Of 68 subjects enrolled, 38 (55.9%), 38 (55.9), and 43 (65.3%) individuals were involved in transmission networks found by Sanger at 3%, UDS at 3% and at 4.5% of MGD, respectively. Mixed transmission chains including HIV-positive and HIV-negative subjects were detected for 8/10 chains by Sanger at 3%, for 9/17 by UDS at 3%, and for 10/18 by UDS at 4.5% of MGD. Overall, the number of HIVnegative individuals clustering with HIV-positive ones was 9/18 by Sanger, 9/18 by UDS at 3%, and 10/18 by UDS at 4.5% of MGD.

Conclusions: HIV-positive and HIV-negative MSM shared HCV transmission networks, which emphasizes the need for HCV surveillance and prevention measures in these communities regardless of the HIV status.

INTRODUCTION

HCV infection among HIV-infected men who have sex with men (MSM) has become an outbreak since several years particularly in urban centres in Europe, Australia, and the United States [START_REF] Götz | A cluster of acute hepatitis C virus infection among men who have sex with men--results from contact tracing and public health implications[END_REF][START_REF] Van De Laar | Evidence of a large, international network of HCV transmission in HIV-positive men who have sex with men[END_REF][START_REF] Giraudon | Increase in diagnosed newly acquired hepatitis C in HIV-positive men who have sex with men across London and Brighton, 2002-2006: is this an outbreak?[END_REF][START_REF] Van De Laar | Increase in HCV incidence among men who have sex with men in Amsterdam most likely caused by sexual transmission[END_REF][START_REF] Gamage | Incidence of hepatitis-C among HIV infected men who have sex with men (MSM) attending a sexual health service: a cohort study[END_REF] . Indeed, a remarkable increase in HCV incidence among HIV-positive MSM was reported by studies of longitudinal cohorts. In France, data from the large Dat'AIDS cohort showed a regular increase of HCV incidence from 4.3 to 11.1 per 1000 person-years (PY) from 2012 to 2016 in French HIV-positive MSM despite a high HCV treatment coverage and cure rate [START_REF] Pradat | Incidence of new hepatitis C virus infection is still increasing in French MSM living with HIV[END_REF] . Data from a meta-analysis of 28 studies revealed a pooled incidence of HCV at 7.8 per 1000 person-year (PY) in HIV-positive MSM while it was only 0.4 per 1000 PY in HIVnegative MSM in resource-rich countries such as Europe, Australia, the United States, and Canada [START_REF] Ghisla | Incidence of hepatitis C in HIV positive and negative men who have sex with men 2000-2016: a systematic review and meta-analysis[END_REF] . Of note, a high incidence of HCV infection in HIV-negative MSM (14 per 1000 PY) was seen in individuals eligible for Pre-Exposure Prophylaxis (PrEP), probably because of their high-risk behaviours [START_REF] Gras | HCV RNA and Antigen Detection for Diagnosis of Acute Hepatitis C Among MSM on PrEP[END_REF] .

Sequencing and phylogenetic analyses are powerful tools to understand transmission dynamics at molecular level. Individuals are considered to share the same transmission chain if their viral populations are more genetically similar to each other than expected by chance and demonstrated by a tight cluster on phylogenetic trees satisfying requirements of branch support value and genetic distance threshold [START_REF] Hué | HIV-1 pol gene variation is sufficient for reconstruction of transmissions in the era of antiretroviral therapy[END_REF] . For example, a collaborative study by phylogenetic approach enrolling HIV-positive MSM recently diagnosed with HCV infection from England, Netherlands, France, Germany, and Australia (n= 226) highlighted a large European MSM specific HCV transmission network 2 . Moreover, several studies showed the spread of HCV strains from HIV-positive toward HIV-negative MSM [START_REF] Hoornenborg | MSM starting preexposure prophylaxis are at risk of hepatitis C virus infection[END_REF][START_REF] Charre | Hepatitis C virus spread from HIV-positive to HIVnegative men who have sex with men[END_REF] . HCV antibody testing is therefore recommended for MSM at high risk of HIV infection and included in PrEP programs.

To date, phylogenetic studies are often based on Sanger sequencing method to identify transmission chains [START_REF] González-Candelas | Molecular evolution in court: analysis of a large hepatitis C virus outbreak from an evolving source[END_REF][START_REF] Féray | A novel method to identify routes of hepatitis C virus transmission[END_REF] . However, Sanger sequencing, with only one bulk or consensus sequence generated, is unable to fully characterize intra-host genetic diversity especially for RNA viruses such as HCV [START_REF] Domingo | Viral quasispecies evolution[END_REF] . Moreover, HCV infections are considered to be frequently established through transmission of minority variants [START_REF] Wang | Hepatitis C virus transmission bottlenecks analyzed by deep sequencing[END_REF][START_REF] Bull | Sequential bottlenecks drive viral evolution in early acute hepatitis C virus infection[END_REF][START_REF] D'arienzo | Sequence and functional analysis of the envelope glycoproteins of hepatitis C virus variants selectively transmitted to a new host[END_REF] ; a consensus sequence cannot therefore reliably capture such transmissions. Ultra-deep sequencing (UDS) with a high throughput of sequencing data allows detecting minority viral populations down to 1% and is able to characterize in-depth viral population. Therefore, in this study, we aimed firstly to identify and characterize HCV transmission chains and secondly to detect closely related HCV transmission events by UDS among HIV-positive and HIV-negative MSM with recent HCV infection.

II. MATERIALS AND METHODS

2.1.Study design and patients

Fifty-five patients with recent HCV infection (50 HIV-positive and 5 HIV-negative), followed at the Pitié-Salpêtrière, Saint-Antoine and Tenon hospitals, Paris, France and 13 HIV-negative patients from the ANRS IPERGAY study (Intervention for prevention of HIV acquisition by antiretroviral therapy for PrEP among gay men at high risk of HIV-1 infection) were also enrolled and considered as possible recent HCV infections.

Sanger and UDS were performed on frozen plasma samples and HCV transmission network was constructed on the 2 datasets to compare quantity and characteristics of transmission chains identified by both techniques.

2.2.Extraction, amplification, and deep-sequencing

Eighty microliters of HCV RNAs were extracted from 1 ml of plasma using® easyMAG® (bioMérieux Clinical Diagnostics). Extracted RNAs were reverse transcribed in complementary DNAs, and NS5B fragment (position 8256 to 8644 compared to H77, fragment of 388 bp) was amplified by PCR in a one-step process (Superscript III One-step RT-PCR with platinum Taq kit; Invitrogen, Carlsbad, CA, USA) according to the manufacturers' protocol, by Forward primer: ATATGAYACCCGCTGYTTTGACTC-3' and Reverse primer : 5'-GCNGARTAYCTVGTCATAGCCTC-3'. Samples were then multiplexed with and subjected to standard Illumina Miseq paired-end sequencing at 2x250 bp.

2.3.UDS data analysis

UDS data were analyzed by Geneious software (version 10.3.2, http://www.geneious.com) [START_REF] Kearse | Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data[END_REF] .

Paired reads were firstly merged, primer-removed and quality-trimmed using quality threshold of 30. Sequences with good quality (quality scores of 30 on at least 95% of bases)

were error-corrected by BBNorm from the BBtools package included in Geneious. Corrected reads of each sample were de novo assembled by Geneious assembler with custom sensitivity following different thresholds of similarity to reduce the number of reads and time for further analysis while maintaining viral population diversity. Firstly, reads were assembled at 100% of similarity. Reads unable to assemble at this threshold were then assembled at 99% of similarity. The process continued and finished at threshold of 97% of similarity where almost all reads were assembled. All contigs and unassembled sequences produced in this step were grouped in one file used for phylogenetic analyses.

2.4.Phylogenetic analysis to study transmission chains

Contig sequences retrieved from clustering process were aligned to a reference sequence corresponding to the subtype of the sample by "Map to Reference" function in Geneious.

Phylogenetic trees were constructed on Sanger dataset (all genotypes in one tree) and on UDS dataset (separate tree for each genotype for better visualization) by FastTree software (version 2.1) [START_REF] Price | FastTree 2--approximately maximum-likelihood trees for large alignments[END_REF] ClusterPicker software was also used to detect individuals belonging to closely related transmission evens, defined as a less than 0.5% of MGD among numerous sequences of different individuals. Adjacent samples in the same PCR or sequencing plate were not considered and presented as closely related transmission events for risk of contamination.

Samples were considered only when they came from different experiments or when they were non-adjacent in the same experiment without any sign of contamination in other samples between them.

Trees were visualized in MEGA7 [START_REF] Kumar | MEGA7: Molecular Evolutionary Genetics Analysis Version 7.0 for Bigger Datasets[END_REF] .

III. RESULTS

3.1.Patients' characteristics and sequencing results

The patients' characteristics are presented in table 1. Briefly, the median age was 38.5 years (IQR, 30.5-46.0) and the majority were MSM (85.3%). Among the 68 patients enrolled, 15

were cases of HCV reinfection, three were considered as possible recent HCV infections, 50 (73.5%) were HIV-coinfected. Significant difference among HIV-positive and HIV-negative patients was observed only for age. HCV sequencing showed genotype 1a (47.1%), 4d

(41.2%), 3a (8.8%), and 2k (2.9%).

A median of 2389 sequences (interquartile range [IQR], 1851-2960) per sample was obtained after quality trimming step.

3.2.Comparison of HCV transmission chains identified by UDS and Sanger

At 3% of MGD, Sanger detected 10 transmission chains in which a median of 3 subjects (min-max = 2-6) was identified while UDS at 3% of MGD detected 17 chains (median = 2 subjects, min-max = 2-5) and UDS at 4.5% of MGD detected 18 chains (median = 2 subjects; min-max=2-6). The number of subjects identified within each transmission chain was not statistically different among Sanger and UDS at 3% and 4.5% of MGD (p value > 0.31 using the Wilcoxon signed rank test in JASP software) 24 .

In particular, UDS allowed detection of hidden transmission chains through minority variants.

UDS at 3% and 4.5% of MGD allowed detection of three and four additional transmission chains, respectively which were not detected by Sanger (table S2 of supplementary data). One transmission chain among these was formed from one individual living in Paris and another residing outside of Paris (from IPERGAY study). Moreover, four subjects were additionally detected by UDS to be included in transmission chains (table S1 of supplementary data).

However, Sanger sequencing also allowed detection of one transmission chain (chain 5 in table S1 of supplementary data) which was not noticed at all by UDS.

3.3.Individuals inside and outside HCV transmission chains

Out of 68 individuals enrolled, 38 (55.9%), 38 (55.9), and 43 (65.3%) were detected to be part of transmission chains by Sanger and UDS at 3% and at 4.5% of MGD, respectively.

Regarding characteristics of individuals inside and outside transmission chains, statistical

analyses showed no significant difference between two groups for age, HIV co-infection and HCV reinfection rates, proportions of MSM, and proportions of individuals infected with HCV-GT1a or HCV-GT4d whatever the technique used (table 2).

3.4.HCV transmission chains including HIV-positive and HIV-negative individuals

HCV transmission chains including HIV-positive and HIV-negative individuals were observed in 8/10 (80%) chains by Sanger, in 9/17 (52.9%) by UDS at 3%, and in 10/18

(55.6%) by UDS at 4.5% of MGD. Overall, among 18 HIV-negative MSM included in this study, the number of HIV-negative individuals clustering with HIV-positive ones was 9 by Sanger, 9 by UDS at 3%, and 10 by UDS at 4.5% of MGD. By UDS at 4.5% of MGD, 8 out of 13 HIV-negative individuals (61.5%) from IPERGAY trial enrolled in this study were detected to belong to transmission chains.

3.5.Closely related HCV transmission events

In a second analysis based on UDS data, we described individuals belonging to closely related transmission events because numerous sequences of different samples were identical or almost identical (MGD <0.5%). Five events considered as closely related transmission were detected in transmission chains number

IV. DISCUSSION

In this work, we identified HCV transmission chains in MSM either co-infected by HIV or at high risk of HIV acquisition in Paris by UDS and Sanger sequencing. Our study revealed a high HCV clustering rate (from 56% to 65%) whatever the techniques used signifying a dynamic transmission among them. Moreover, one patient under PrEP living outside Paris was enrolled and this patient was found to be part of a transmission chain by UDS. Therefore, in case of HCV infection, early initiation of treatment should be carried out in this population to rapidly prevent further spread of the virus.

Transmission chains were identified at cut-off of 3% of MGD by Sanger and at two different cut-offs of 3% and 4.5% of MGD by UDS. Indeed, few studies have conclusively established the cut-off of MGD to identify a transmission chain among HCV-infected people by UDS, varying from 2 to 4.5% [START_REF] Charre | Hepatitis C virus spread from HIV-positive to HIVnegative men who have sex with men[END_REF][START_REF] Olmstead | A molecular phylogenetics-based approach for identifying recent hepatitis C virus transmission events[END_REF][START_REF] Montoya | Deep sequencing increases hepatitis C virus phylogenetic cluster detection compared to Sanger sequencing[END_REF][START_REF] Campo | Accurate Genetic Detection of Hepatitis C Virus Transmissions in Outbreak Settings[END_REF][START_REF] Rose | Complex patterns of Hepatitis-C virus longitudinal clustering in a high-risk population[END_REF] . It may be difficult to compare Sanger and UDS techniques at the same cut-off of MGD because UDS allows a much deeper characterization of viral diversity. In this work, NS5B deep sequencing improved the discrimination of transmission chains versus Sanger sequencing which was in line with results from a study of Montoya et al. [START_REF] Montoya | Deep sequencing increases hepatitis C virus phylogenetic cluster detection compared to Sanger sequencing[END_REF] . UDS at both thresholds of MGD identified a median of two subjects within a transmission chain compared to a median of three subjects by Sanger sequencing but this difference was not statistically significant. Importantly, UDS allowed establishing more solid transmission events and the transmission dynamics among individuals within each chain could be further evaluated through detection of numerous clustered viral strains. For example, we detected five transmission events considered very closely related i.e. individuals harboured viruses with numerous overlapped sequences (MGD < 0.5%). Importantly, among them, some harboured viruses with multiple identical sequences (MGD = 0%) suggesting direct transmission events [START_REF] Romero-Severson | Phylogenetically resolving epidemiologic linkage[END_REF] . However, it is not possible to confirm direct transmission from one person to another using molecular data alone. Indeed, both could be infected from a third source, or they could be connected indirectly through a transmission chain including one or more intermediaries.

Although transmission directionality was not inferred due to lack of specific epidemiological data, patients included in these events should be followed more closely including the communities around them to assure a rapid intervention.

In this study, UDS also detected hidden transmission chains by identifying transmission linkages through minority viral strains. However, the deeper characterization of viral variability is also the reason why UDS did not detect one transmission chain found by Sanger technique. The MGD among sequences of the three individuals involved in this transmission chain was 2.64% with Sanger sequencing while the genetic distance among viral sequences of the three individuals is higher than the MGD threshold of 3% and 4.5% with UDS. Therefore, UDS did not capture the transmission linkage among these individuals as Sanger did. Further studies would be necessary to determine the most suitable MGD cut-off for transmission chain identification by UDS. Thereby, UDS would be interesting to deeply characterize transmission patterns such as directness or directionality among individuals; however, it is not more useful than Sanger sequencing in term of large-scale prevention and rapid intervention.

Importantly, depending on the techniques and MGD cut-off used, 53% to 80% of transmission chains identified included both HIV-positive and HIV-negative subjects and more than 50% of HIV-negative subjects enrolled in this study clustering with HIV-positive ones. The shared HCV transmission networks among HIV-positive and HIV-negative MSM were also observed in two studies conducted in Amsterdam, the Netherlands and in Lyon, France [START_REF] Hoornenborg | MSM starting preexposure prophylaxis are at risk of hepatitis C virus infection[END_REF][START_REF] Charre | Hepatitis C virus spread from HIV-positive to HIVnegative men who have sex with men[END_REF] . Our results raise an alert for better screening, monitoring, and surveillance of HCV infection in this high-risk community regardless of the HIV status.

Even though the HCV reinfection rate in subjects inside transmission chains was not statistically higher than in subjects outside them, the need of follow-up for possible HCV One limitation of the study is that 14.7% of patients had unknown sexual orientation.

However, they all engaged in risky behaviours for multiple HCV exposures as other MSM enrolled in the study, demonstrated by their HIV co-infection status, HCV reinfection rate, and other sexual transmission diseases discovered with their HCV infection. Furthermore, individuals with unknown sexual orientation were also involved in transmission chains and no significant difference was observed in proportions of MSM or unknown sexual orientation subjects among those inside and outside transmission chains. Another limitation in our study was the length of fragment NS5B sequenced. Indeed, a longer sequence could permit more accurate differentiation of linked or unlinked virus and thus more exactly identify transmission chains [START_REF] Lamoury | The Influence of Hepatitis C Virus Genetic Region on Phylogenetic Clustering Analysis[END_REF] . In our study, a quite short fragment of NS5B was amplified but it was counterbalanced by a high depth of coverage by UDS. Furthermore, this strategy had been applied in different settings [START_REF] Wang | Hepatitis C virus transmission bottlenecks analyzed by deep sequencing[END_REF][START_REF] Rossi | Multiregion deep sequencing of hepatitis C virus: An improved approach for genetic relatedness studies[END_REF] . In this study, we did not have enough epidemiological data to confirm the true transmission events among individuals. However, the fact that almost all individuals (67/68) were from Paris justified in some way their epidemiological connection as well as the fact that transmission chains identified by UDS were established through multiple clustered viral strains increased the likeliness of true transmission event identification.

In conclusion, in this study, a high clustering rate of HCV was observed in HIV-positive and HIV-negative MSM communities in Paris, particularly those engaged in PrEP program.

Furthermore, HIV-positive MSM shared HCV transmission networks with HIV-negative MSM. The more frequently screening and surveillance of HCV infection regardless of the HIV status is essential to prevent the spread of HCV in these high-risk communities.

  3, 4, 7, and 8 in table S1 of supplementary data. In detail, we detected this event between individuals 9 and 10 (2 months of difference in date of HCV infection); 13 and 14 in chain 3 (14 months), among individuals 15, 16, 18, and 19 in chain 4 (14 months), among individuals 27, 28, and 29 in chain 7 (5 months), and among

	individuals 30, 33, 34, and 35 (39 months) in chain 8.
	Examples of two phylogenetic trees constructed from UDS sequences from two and three
	individuals considered as closely related transmission events are shown in figures S1 and S2
	of supplementary data.

  reinfection and of patient support and education to prevent HCV reinfection and transmission arise in this high-risk population. Last but not least, 61.5% of HIV-negative individuals under PrEP enrolled in this study were detected to belong to transmission chains by UDS at 4.5% of MGD. Therefore, surveillance of HCV infection by HCV viral load instead of anti-HCV antibody test would be more advantageous to rapidly intervene and control transmission for those in PrEP programs.

An interesting study by Caro-Pérez et al. has showed an HCV outbreak in HIV-positive MSM in Barcelona related to a previously described European MSM transmission network

[START_REF] Caro-Pérez | Phylogenetic analysis of an epidemic outbreak of acute hepatitis C in HIV-infected patients by ultra-deep pyrosequencing[END_REF]

. For that reason, our high throughput of HCV sequencing data from HIV-positive and HIVnegative MSM in Paris will be further investigated to study the transmission network of these populations with HCV sequences of other MSM at European level 2 .

Acknowledgements

We would like to thank all the patients who agreed to participate in the study, all the participant doctors who followed the patients, in particular Drs. ROUDIERE, LIOTIER, GOSSET, CARDON, GRIVOIS, ISRAEL, KIRSTETTER, LAYLAVOIX, BOTTERO, WORMSER and Pr. KATLAMA and all the participant virologists, Drs. ELAERTS and SCHNEIDER.

We thank the INSERM SC10 and the Trial Scientific Committee for IPERGAY trial.

We thank ANRS AC43 Next Generation Sequencing and STIs working groups for their support.

Funding

This work received financial support from the Agence Nationale de Recherches sur le SIDA et les hépatites virales (ANRS) (decision number 2018-139).

Presented in part

This study was presented in part at the International AIDS Conference, Amsterdam, the Netherlands, 2018 (Abstract THAB0203, oral presentation).

Conflict of interest

All authors do not have any conflicts of interest to declare.