Bastien Fernandez
email: bastien.fernandez@cpt.univ-mrs.fr

Global synchronisation in translation invariant coupled map lattices

 where CML with reflection invariant couplings are considered, the condition only involves the linearised dynamics in the diagonal, namely for all points in the diagonal, the derivative must be contractive in all transverse directions. In addition to this result, a (weaker) condition that ensures the CML attractor to be composed of either 2-periodic or constant configurations, is also obtained.

Synchronisation is probably the most commonly observed dynamical phenomenon in interacting (coupled) nonlinear systems. Generally speaking, it is said to occur in a multidimensional dynamical system when the attractor lies inside a one-dimensional subset of the phase space so that the knowledge of the trajectory of a single coordinate allows to determine all coordinate trajectories [Boccaletti et al., 2006;[START_REF] Pikovsky | Synchronization. A universal concept in nonlinear sciences[END_REF]. In other words, all components asymptotically evolve "in phase", possibly in a chaotic motion [START_REF] Pecora | Synchronisation of chaotic systems[END_REF]. In practice, the phenomenon takes various forms upon the system under consideration. For instance, phase synchronisation takes place in systems of coupled oscillators [START_REF] Fujisaka | Stability theory of synchronized motion in coupled oscillator systems[END_REF]] and master-slave, or generalized synchronisation are possible phenomena in unidirectionally coupled systems [START_REF] Hunt | Differentiable generalized synchronization of chaos[END_REF][START_REF] Rulkov | Multivalued mappings in generalized chaos synchronization[END_REF][START_REF] Tresser | Master-slave synchronisation from the point of view of global dynamics[END_REF].

Another specific form of synchronisation happens in dynamical systems with symmetry where it consists in the convergence to the symmetry fixed point set (which is invariant under dynamics, see e.g. [START_REF] Ashwin | From attractor to chaotic saddle: a tale of transverse instability[END_REF]). Challenging problems in this framework then concern conditions on parameters (or on components) for synchronisation. These conditions often involve the notion of transverse Lyapunov exponent in order to specify the basin of attraction of the symmetry fixed point set. In particular, a seminal result [START_REF] Alexander | Riddled basins[END_REF] states that this basin has positive Lebesgue measure (in phase space) when all transverse Lyapunov exponents are negative, for Lebesgue almost every point in the invariant set. Moreover, it is likely to have a fractal "riddled structure" (see the chapter by P. Ashwin in [Chazottes & Fernandez, 2005] for a mathematical definition) and has actually such a structure in various examples. However, [START_REF] Ashwin | From attractor to chaotic saddle: a tale of transverse instability[END_REF] showed that this basin is indeed a neighbourhood of the invariant set when the supremum of all transverse Lyapunov exponents with respect to all ergodic measures supported in the invariant set, is negative. Naturally, and as observed in [START_REF] Ashwin | From attractor to chaotic saddle: a tale of transverse instability[END_REF], without any further specification of the dynamics, these results are optimal and one cannot expect to control the global behaviour. In some cases however, it is possible to determine the fate of every orbit in phase space from properties of the linearised dynamics in the invariant set.

CML have been introduced at the beginning of the eighties as simple discrete-time models of reactiondiffusion systems [Chazottes & Fernandez, 2005;[START_REF] Kaneko | Synchronization analysis of linearly coupled networks of discrete time systems[END_REF]. Their specificity resides in the definition of their mapping which is the composition of an individual map and of a linear coupling. This presents the double advantage of being well-adapted to numerical simulations and to mathematical analysis. Recently, [START_REF] Kaneko | Synchronization analysis of linearly coupled networks of discrete time systems[END_REF]] has established and analyzed synchronisation conditions in CML with arbitrary linear coupling operators, mostly in absence of symmetry but yet with invariant diagonal.

For reflection invariant coupled map lattices (CML), [START_REF] Jost | Spectral properties and synchronization in coupled map lattices[END_REF] proved that if all transverse eigenvalues of the jacobian matrix are contractive for all points in the diagonal (the invariant set in this case), then all points in phase space asymptotically approach the diagonal, a property called global synchronisation. Although slightly stronger than the previous one, this condition is much simpler to check in practice.

In this Letter, we focus on translation invariant (but not necessarily reflection invariant) CML on periodic lattices. Translation invariance is usually assumed in most studies [Chazottes & Fernandez, 2005;[START_REF] Kaneko | Synchronization analysis of linearly coupled networks of discrete time systems[END_REF] as it reflects the simplifying assumptions that the individual systems are all identical and that the coupling is of diffusive type. As in [START_REF] Jost | Spectral properties and synchronization in coupled map lattices[END_REF], we show that global synchronisation holds provided that all transverse directions are contractive for all points in the diagonal. Our condition is actually a bit more general than the one in [START_REF] Jost | Spectral properties and synchronization in coupled map lattices[END_REF]] and the technique is different. The results in particular complete a previous result in [START_REF] Lin | Synchronisation in coupled map lattices with periodic boundary condition[END_REF] on global synchronisation for lattices of 2,3 and 4 coupled logistic maps.

A translation invariant CML on the periodic lattice Z L := {s ∈ Z modL} (L > 1) is the dynamical system generated by the following induction relation in R

x t+1 s = n∈Z L c n f (x t s-n ), s ∈ Z L . (1)
Here the coefficients c n are non-negative (c n 0) and normalized (

n∈Z L c n = 1). The most frequent
examples are respectively, the asymmetric nearest neighbour coupling, the symmetric one (L > 2) and the global coupling for which the coefficients are respectively given by ( ∈ [0, 1])

c n =    1 - if n = 0 if n = 1 0 otherwise , c n =    1 - if n = 0 2 if n = ±1 0 otherwise and c n = 1 - if n = 0 L otherwise
More generally, a coupling on the periodic lattice with L sites can be defined from any normalised sequence {γ n } n∈Z of non-negative coefficients by summing over the periods, namely

c n = k∈Z γ n+kL .
The individual map f in (1) is assumed to possess an invariant interval I (which may be the whole R) on which it satisfies the following inequality for some K > 0 (i.e. f is Lipschitz continuous)

|f (x) -f (y)| K|x -y|, ∀x, y ∈ I.
We assume that K is the smallest of such numbers. Note that if f is continuously differentiable, then K is the maximum of |f (x)| on I. However, f needs not be differentiable for our purpose.

The CML (1) leaves invariant the diagonal in I Z L . For any point x s = x (s ∈ Z L ) on this diagonal, the eigenvalues of the jacobian matrix associated with the CML derivative (which exists when f is differentiable) can be easily computed. Indeed, since this matrix commutes with translations on the lattice Z L , its eigenvectors in In particular, f (x)ĉ 0 = f (x) is the eigenvalue along the diagonal and all other eigenvalues correspond to transverse directions. The quantities ĉk and ĉL-k are complex conjugate for k = 0, L/2 (and equal in the case of reflection invariant couplings, i.e. c -n = c n ). This implies that the contraction rate in the subset of R Z L normal to the diagonal and defined by

C Z L are the elements e k (k = 0, • • • , L -1)
{ae k + āe L-k : a ∈ C} if k = 0, L/2 (and by ae L/2 : a ∈ R if k = L/2), is equal to |ĉ k |.
As a consequence, the normal Lyapunov exponents for any point on the diagonal write λ(x)

+ log |ĉ k | (k = 1, • • • L -1)
where λ(x) is the individual map Lyapunov exponent evaluated at x ∈ I. According to Theorem 2.12 in [START_REF] Ashwin | From attractor to chaotic saddle: a tale of transverse instability[END_REF], the condition log(K) + log max

1 k L-1 |ĉ k | < 0
ensures the existence of a neighbourhood of the diagonal in I Z L in which the orbit of every point asymptotically approaches the diagonal. Due to the specific structure of the CML (1), the very same condition actually implies that convergence to the diagonal holds for every point in phase space in such systems.

Theorem. The condition K max

1 k L-1 |ĉ k | < 1, (2)
implies global synchronisation of the CML. That is to say, the following limit holds

lim t→∞ max s∈Z L |x t s -x t s+1 | = 0 for all initial configurations {x 0 s } ∈ I Z L .
(The proof is given below.) Several comments on this result can be made. First, since the definition of ĉk implies max

1 k L-1
|ĉ k | 1, the condition (2) may hold in the case where K > 1 and in particular for chaotic individual maps as it happens in the case of chaotic synchronisation of CML. (On the other hand, global synchronisation always hold when K < 1 because the whole CML is contracting in this case.)

Moreover, the coupling eigenvalues have the following limit (whose proof is given after the proof of the Theorem) lim

L→∞ max 1 k L-1 |ĉ k | = 1. (3)
This implies that, for any coupling coefficient sequence {γ n } n∈Z and given individual map (with K > 1), the condition (2) becomes stronger and is hardly satisfied when the lattice size increases (and actually exceeds few sites in common examples). Alternatively, given a coupling coefficient sequence {γ n } n∈Z and a number of sites L, the condition (2) may not hold whenever K exceeds some threshold. For instance, for the symmetric nearest neighbour coupling defined above, L = 2 or 3 and K > 1, the condition (2) holds for all in some interval contained inside [0, 1]. However, if L = 4 (resp. L = 6), the same condition holds -for at least one value of in [0, 1] -only if K < 1/3 (resp. K < 3/5).

Proof of the Theorem. The proof is not as trivial as the presentation may suggest. The natural strategy which consists in controlling the dynamics of transverse modes by applying the triangle inequality, does not work for L > 2. Indeed, this method provides upper bounds for both the modulus of transverse modes at time t and at time t+1, but no direct relationship between the two quantities. The alternative strategy which consists in directly estimating the difference x t s -x t s+1 does not work either when L exceeds a threshold, because it implies the whole spectrum and leads to the condition K < 1. Therefore, one has to find a better combination of coordinates in order to find the appropriate expression.

We consider the euclidean scalar product and associated norm • in C Z L and the euclidean scalar product and associated norm

• R in R Z L . Let M 1 be the diagonal in I Z L and let M ⊥ 1 := {x ∈ I Z L : s∈Z L x s = 0}
be its orthogonal complement. Let also C denotes the coupling operator, i.e.

(C(x)) s = n∈Z L c n x s-n , s ∈ Z L
where x = {x s } s∈Z L . In a first step, we show that global synchronisation occurs under the condition x t -R(x t ) R = 0, for all x 0 ∈ I Z L . The CML dynamics commutes with R and since x -R(x) ∈ M ⊥ 1 for all x ∈ I Z L , we obtain the following inequality

K max x∈M ⊥ 1 : x R =1 C(x) R < 1. Let R
x t+1 -R(x t+1 ) R K max x∈M ⊥ 1 : x R =1 C(x) R x t -R(x t ) R
An induction then completes the first step.

In a second step, we show that max

x∈M ⊥ 1 : x R =1 C(x) R = max 1 k L-1 |ĉ k |.
We identify every vector {x s } ∈ R Z L with the vector {z s } ∈ C Z L where z s = x s for all s. It follows that every vector x ∈ M ⊥ 1 can be viewed as belonging to the subspace orthogonal to the diagonal in C Z L . Therefore, this vector writes x = L-1 k=1

x k e k with complex coordinates x k . Since the e k are the eigenvectors of C in C Z L with eigenvalues ĉk we have

(C(x) = L-1 k=1
x k ĉk e k , and then C(x) max

1 k L-1 |ĉ k | x . Since C(x) ∈ R Z L and x = x R , we have actually showed that max x∈M ⊥ 1 : x R =1 C(x) R max 1 k L-1 |ĉ k | .
It only remains to shows that the inequality is indeed an equality. Let k ∈ {1, • • • , L -1} be such that The sequence {γ n } is assumed to be summable. Thus for every δ > 0, there exists n δ such that

|ĉ k | = max 1 k L-1 |ĉ k | . The normalised vector v = 1 √ 2 (e k + e L-k ) (resp. v = e L/2 if k = L/
|n|>n δ γ n e -2iπn L |n|>n δ γ n < δ.
Now, let δ < 1 and define L δ to be sufficiently large so that for all L L δ , we have cos 2πn L 1 -δ for all |n| n δ . The two estimates imply that the following inequality hold for L L δ

|n| n δ γ n e -2iπn L |n| n δ γ n cos 2πn L (1 -δ) |n| n δ γ n (1 -δ) 2 It results that |ĉ 1 | (1 -δ)
2 -δ for L L δ and the desired limit (3) follows. 2

To conclude this Letter, we mention that global synchronisation can be seen as a special case of convergence to a spatially periodic configuration subset with prescribed period, namely period 1 in the present case. The arguments of the proof above can be adapted in order to obtain a condition such that the CML attractor is composed of periodic configurations with (not necessarily minimal) period, a given divisor of L. We only consider here the situation where this period is 2 and we provide the corresponding condition.

Proposition. Assume that L is even and that the following condition holds (-1) s x s = 0 is the hyperplane orthogonal both to the diagonal and to the linear subspace generated by e L/2 . As in the proof of the Theorem, we conclude that the condition K max

x∈M ⊥ 2 : x R =1
C(x) R < 1 implies the limit lim t→+∞ x t -R 2 (x t ) R = 0 which is equivalent to the desired result. Now, by definition of M ⊥ 2 , every configuration x ∈ M ⊥ 2 (viewed as an element of C Z L ) writes x =

1 k L-1, k =L/2
x k e k . As in the proof of the Theorem, this implies that max

x∈M ⊥ 2 : x R =1 C(x) R = max 1 k L-1 : k =L/2 |ĉ k |
The Proposition then easily follows. 2

  of the Fourier basis, where (e k ) s = 1 √ L e 2iπks L (s ∈ Z L ). A simple calculation shows that the eigenvalues are given by f (x)ĉ k where ĉk = n∈Z L c n e -2iπkn L .

  denotes the (right) translation operator. Since the norms • R and max s∈Z L |x s | are equivalent, global synchronisation is equivalent to the limit lim t→∞

  2) has real coordinates, belongs to M ⊥ 1 and satisfies C(v) = |ĉ k | since ĉk and ĉL-k are complex conjugate. In other words, there exists a normalised vector v ∈ M ⊥ 1 such that C(v) R = max1 k L-1 |ĉ k |.The Theorem is proved.Proof of relation (3). Firstly, the triangle inequality implies that for every n > 0 we have n∈Z γ n e -

  We start by observing that for all x ∈ I Z L , we have x -R 2 (x) ∈ M ⊥ 2 whereM ⊥ 2 = x ∈ I Z L : s∈Z Lx s = 0 and s∈Z L

Acknowledgments

I warmly thank V. Afraimovich for stimulating and fruitful discussions. I am also indebted to P. Guiraud for a careful reading of the manuscript which allowed to improve the presentation.


Diana Vaclavkova

Alex Delhomme

Clément Faugeras

Marek Potemski

Aleksander Bogucki

Jan Suffczynski

Piotr Kossacki

Andrew R Wildes

Benoît Grémaud

Andrés Saúl

Jan Suffczyński

Magnetoelastic interaction in the two-dimensional magnetic material MnPS3 studied by first principles calculations and Raman experiments

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

I. INTRODUCTION

The research dedicated to the so called twodimensional (2D) materials [1] has been lately expanded to studies of a large class of layered magnetic materials and their ultrathin crystals [2]. This, to great extent, is stimulated by the search for novel effects in the magnetic ordering in the 2D limit and possible design of new devices which take advantage of the spin degree of freedom. Relevant research efforts are, in particular, being focused on antiferromagnetic transition-metal thiophosphates [3], with their representative example, the MnPS 3 compound. Changes in the magnetic ordering are often traced with Raman scattering, a technique which appears to be especially well suited to study small size and ultrathin samples. The elementary excitations which are most often probed with Raman scattering are phonon excitations. Answering the question how and which phonon modes appear to be sensitive to the magnetic order is central for drawing conclusions on the applicability of Raman scattering methods to study magnetic phase transitions.

In this paper we present a first-principle study of magnetoelastic interaction in MnPS 3 antiferromagnets and compare our theoretical expectations with the results of temperature dependent Raman scattering measurements performed on these crystals. The aim of our work is to identify the phonon modes which are sensitive to magnetic order and to enlighten the mechanism of the phonon-spin coupling that is at the origin of this sensitivity. * saul@cinam.univ-mrs.fr Experimentally, we confirm [4] a robust sensitivity of the Raman active phonon mode at 155.1 cm -1 (at low temperatures) to the antiferromagnetic to paramagnetic transition in MnPS 3 and clarify a similar property for the other, but weaker, phonon mode observed at 120.2 cm -1 . These experimental observations are well accounted for by our theoretical calculations.

This paper is organized as follows. In Sections II we present the crystallographic structure of MnPS 3 and the details of the Raman experiments. In Section III we discuss the methods that we used for the calculation of the effective exchange interactions and the phonons modes. In Sections IV we present and discuss our Raman experiments versus temperature and in Section V the results of our calculations. Section VI is devoted to the discussion and conclusion.

II. EXPERIMENTAL DETAILS

A. MnPS3 structure MnPS 3 has a monoclinic structure, space group C2/m with lattice parameters a = 6.077 Å, b = 10.524 Å, c = 6.796 Å, and β = 107.35 • [6]. The atomic structure is shown in Figures 1(a) and (b). There are 4 MnPS 3 units in the monoclinic C2/m unit cell, giving a total of 20 atoms. The C-centring means that the unit cell is not primitive because there are two equivalent positions, one at (0, 0 ,0) and the other at (1/2, 1/2, 0). A primitive triclinic unit cell containing 10 atoms and having the minimal volume V = 1/2 abc sin β can be constructed with, for example, Bravais lattice vectors a 1 = (ax -bŷ)/2, a 2 = (ax + bŷ)/2, and a 3 = c (cos β x + sin β ẑ). In the paramagnetic phase all the Mn and P atoms are equivalent by symmetry, while there are two kinds of S atoms. The corresponding Wyckoff positions of the C2/m cell are 4g, 4i, 4i, and 8j for the Mn, P, S1, and S2 atoms, respectively.

As can be seen in Fig. 1(a), MnPS 3 is formed by the stacking along the c axis of PS-Mn-PS sandwiches. The Mn atoms form a distorted honeycomb lattice and the S atoms two distorted triangular lattices. The P atoms can be described as forming P 2 dimers penetrating the Mn hexagons (see Fig. 1(a)). Each Mn atom is octahedrally coordinated to six S atoms which can be associated to super-exchange paths between the Mn atoms (see Fig. 1(b) and the text below).

B. Raman scattering experiments

Raman scattering measurements were carried out on single crystals of manganese phosphorus trisulphide (MnPS 3 ) using, as the excitation source, either a 515 nm line of a continuous wave laser diode or a 514.5 nm line of a continuous wave Ar + laser. The crystals used for experiments were either home-made (see Ref. [7] for more details) or commercially available, both types of specimens displaying essentially identical spectra. In accordance to earlier studies [4], the Raman scattering in our samples has been found to be particularly efficient when using the laser excitation wavelengths around 515 nm. The collected signal was dispersed through a 0.5 m monochromator and detected by a CCD camera.

The temperature dependent Raman scattering response was measured at zero magnetic field with the aid of a continuous flow cryostat mounted on x-y motorized positioners. The sample was placed on a cold finger of the cryostat. The excitation light was focused by means of a 50x long-working distance objective with a 0.5 numerical aperture producing a spot of about 1 µm and the scattering signal was collected via the same objective.

Magneto-Raman scattering was measured in Faraday and Voigt configurations using an optical-fiber-based experimental setup. The sample was mounted on top of an x-y-z piezo-stage kept in gaseous helium at T= 4.2 K and inserted into a magnet. The excitation light was coupled to an optical fiber with a core of 5 µm diameter and focused on the sample by an aspheric lens (spot diameter around 1 µm). The signal was collected by the same lens, injected into a second optical fiber of 50 µm diam-eter. In all experiments, the excitation power was kept below 100 µW, to avoid the sample heating by the laser beam. Multiple test measurements of the power dependence of Raman scattering response have not shown any changes in the shape of spectra in the range of excitation powers up to 500 µW.

III. CALCULATION DETAILS

For the calculations we have used the Quantum Espresso [8] code based on density functional theory, Optimized Norm-Conserving Vanderbilt (ONCV) pseudopotentials and the PBE functional [9] with a planewave and charge-density cutoff of 90 Ry and 200 Ry, respectively. An electronic structure calculation with the experimental atomic parameters [6] and a PBE functional gives no gap for the ferromagnetic (FM) order and a gap of 1.30 eV for the experimentally observed AFM order to be compared with the experimentally reported gap of 2.96 eV [10].

Including an on-site Hubbard U eff on the Mn atoms using the simplified method developed by Cococcioni and de Gironcoli [11] increases the AFM gap up to a saturation value of 2.60 eV when U eff reaches 7 eV. The same dependence of the electronic gap with the value of U eff has has been reported by other authors, who use either U eff = 4 eV [12] or 5 eV [13,14]. A value of U eff = 5 eV (gap of 2.40 eV) has been used for the calculations of the effective exchange interactions and vibrational modes reported below.

The experimental lattice parameters and internal coordinates [10] has been used to calculate the effective exchange interactions. As mentioned above, the non primitive 1×1×1 base centered monoclinic cell contains 20 atoms. Double monoclinic super-cells (1×1×2, 1×2×1, and 2×1×1) containing 40 atoms and 8 Mn have been used to be able to separate the different exchange interactions (see below). We have used a 8×4×2 Monkhorst-Pack [15] grid for the first Brillouin zone sampling of the 40 atoms monoclinic 1×1×2 unit cell and adapted equivalent samplings of 8×2×4 and 4×4×4 for the 1×2×1 and 2×1×1 unit cells.

The phonon modes have been calculated in the primitive 10 atoms triclinic cell with the experimental lattice parameters but fully relaxed internal coordinates under the constraint of keeping the unit cell symmetry invariant. After relaxation, the calculation has been performed using the frozen phonon procedure with a 8×8×8 Monkhorst-Pack [15] grid for the first Brillouin zone sampling.

IV. EXPERIMENTAL RESULTS

The representative Raman scattering spectra of our MnPS 3 crystals, measured at low (5 K) temperature in the range of 100 to 650 cm -1 from the laser line, are presented in Fig. 2. A number of relatively sharp Raman scattering peaks, labeled X1 . . . X8 are observed, in fair agreement with previously reported observations of the specific phonon modes in MnPS 3 crystals [4,[16][17][18][19]. The temperature evolution of all but the X1 and X2 peaks is rather smooth. Indeed, a conventional red shift and broadening of X3 . . . X8 phonon modes is observed upon increase of temperature (see Supplementary Material Figures S1 andS2). In contrast, the X1 and X2 Raman scattering peaks appear to be sensitive to the magnetic ordering in our crystals. These peaks display a singular behavior around the expected, critical Néel temperature, T N = 78 K, for the antiferromagnetic to paramagnetic transition in MnPS 3 compound [20]. This is illustrated in Fig. 3. As can be seen, the center positions of the X1 and X2 Raman scattering peaks display an abrupt blue shift in close vicinity of T N . In agreement with previous reports [4,19], the observed blue shift amounts 5 cm -1 in the case of the X2 mode and we report here a similar behavior but with a smaller blue shift of about 2 cm -1 for the X1 Raman scattering peak (see Fig. 3(b) and (c)). More details on temperature evolution of X1 and X2 peaks as well as of all other X3 . . . X8 peaks are shown in the Supplementary Material's Figures S1 to S3.

Here we note that the X2 peak is relatively intense and well pronounced in the spectra measured in the whole temperature range, though it appears to be significantly broadened (by a factor of 4) around T N (the maximum in the X2 line-width is observed at temperature of 65 K, slightly below T N ). The X1 peak is relatively weak. It does not display any obvious broadening around T N but displays a relatively abrupt broadening (by a factor of 2) at T 110 K. The intensity of X1 peak initially weakens but significantly increases when temperature is raised above 100 K.

It is interesting to note that, as illustrated in Fig. 4, both X1 and X2 peaks are not sensitive to the applied magnetic field. As already perceived in a recent work [4] the low temperature Raman scattering spectra of MnPS 3 crystals show an additional weak but intriguing feature at energies slightly above the X2 phonon peak. This additional feature marked as X' in Fig. 2(b), might be seen as a continuous-like spectrum with a peak-like, highenergy onset around 183.3 cm -1 at low temperatures. A weak, Raman scattering feature centered at 185 cm -1 has been previously observed in MnPS 3 compounds studied at room temperature [16], though it is not observed in our experiments. Instead, as shown in Fig. 5, our Raman scattering studies performed with longer accumulation time and fine temperature intervals, show that the X' feature exhibits a characteristic narrowing effect (red shift of the high energy onset) upon the increase of temperature. It merges with the X2 phonon peak at temperatures around T N and it is not observed in the spectra measured above T N .

The X' feature is likely of a magnetic origin, and the first guess would be that it is a due to a coupled phonon (X2)-spin-wave (magnon) excitation. Such cou-pling would be favored for k ≈ 0 spin-wave modes whose energy as measured by the spin gap in the magnon dispersion relation [7,21] is only a few cm -1 (∼ 0.5 meV ∼ 4.03 cm -1 ). This is in contradiction to the rather large extension (∼ 20 cm -1 ) of the X' feature as measured at low temperatures (see Fig. 5). The above hypothesis can also be ruled out by our experiments performed as a function of the magnetic field (see Fig. 4) in which we do not observe any magnetic-field-induced change in the form of the X' feature. Very likely, and in accordance with the recent report [4], the X' feature is due to twomagnon excitations. Indeed, as deduced from neutron scattering studies [7] the (single) magnon excitations at the Brillouin zone boundary saturate at energies ∼ 11.5 meV ∼ 92 cm -1 , thus the two-magnon density of states is expected to peak at around 184 cm -1 , in accordance to the energy position of the onset of the X' feature measured at low temperatures. Notably, the magnon dispersion in the long-wavevector limit, i.e., close to the zone-boundary, may persist in the paramagnetic phase (above T N ) but the Raman scattering response due to two-magnon excitations can be modified at the antiferromagnetic or paraferromagnetic phase transistion [22,23]. Nevertheless, it remains surprising that our experiments do not show any trace of these two-magnon excitations at temperatures above T N [22,23].

FIG. 5. False-color map of Raman scattering response (logarithmic scale for the intensity), together with few selected spectra, measured with long acquisition times and fine temperature interval of 2 K in the spectral region 100 to 200 cm -1 and temperature range from 5 K to 100 K. The presented data emphasize a characteristic "narrowing" of the X' feature with temperature (red shift of its high energy onset), which eventually merges with the X2 phonon mode at TN . This behavior suggests that the X' feature is due to spin (possibly two-magnon) excitations in the antiferromagnetic phase of MnPS3.

V. CALCULATION RESULTS

A. The magnetic Hamiltonian

To understand the possible link between the vibrational modes and the magnetic order in MnPS 3 one needs to determine the model magnetic Hamiltonian of the interacting Mn 2+ ions.

If one considers that these ions are in a 3d 5 electronic configuration corresponding to a high spin 6 S (S = 5/2, L = 0) ground state [20], the effect of the spin-orbit coupling would be negligible [24]. The observation of a highly isotropic magnetic susceptibility above T N sustains this idea [7,[24][25][26][27] and suggests that the effective exchange interactions J ij between the magnetic centers are isotropic in nature.

Nevertheless, the AFM order below T N is anisotropic. It was first reported with the Mn spin moments pointing along the reciprocal lattice vector c * , i.e., perpendicular to the honeycomb lattice in the ab plane [7,20,25,26,28] but later on, a more careful analysis of the neutron diffraction data has shown a small tilt (∼ 8 • ) of the Mn 2+ magnetic moments with respect to c * [29]. Two posible energetic contributions can be proposed to explain the anisotropic order, the dipolar interactions between the magnetic moments or a single ion anisotropy arising from a non negligible spin-orbit coupling. Theoretical analysis [30][31][32][33] and recent magnon band measurements [21] suggest that the dipolar interaction would be the leading source of anisotropy responsible of the observed longrange order. While other experimental results (electron spin resonance [25,34], critical behavior [20,35], and the observation of a non negligible orbital contribution to the total magnetic moment [26,28]), support the idea that the single ion anisotropy is also playing an important role.

In conclusion the full magnetic Hamiltonian should contain a spin-independent contribution Ĥ0 , a rotational invariant Heisenberg part ĤHei , a dipolar term Ĥdip , and a single ion anisotropy contribution Ĥion :

Ĥ = Ĥ0 + ĤHei + Ĥdip + Ĥion (1) ĤHei = i>j J ij Ŝi • Ŝj (2) Ĥdip = - µ 0 4π i>j 1 r 3 ij [( mi • rij )( mj • rij ) -mi • mj ] (3) Ĥion = D i ( Ŝi • n) 2 (4)
where J ij are the magnetic couplings, Ŝi and Ŝj are the S = 5 2 spin operators localized on the Mn atoms at site i and j, respectively, mi = gµ B Ŝi and mj = gµ B Ŝj the corresponding spin only magnetic moments, and n the direction of the single ion easy axis. Notice that our Heisenberg interaction is defined as J S i • S j , similar to Le Flem et al [36], Sivadas et al [13], and Chittari et al [12], while other authors use -2J S i • S j [7,[24][25][26][27].

B. Effective exchange interactions

The full analysis of the magnetic Hamiltonian will be discussed in a forthcoming paper. To understand the link between the vibrational and magnetic properties we will only focus here on the rotational invariant Heisenberg contribution ĤHei .

We have classified the effective exchange interactions according to the Mn-Mn distances and we have calculated 9 different interactions up to interaction distances of 7.647 Å (see Table I and Figures 1(c) and(d)). It I. Effective exchange interactions. The nine interactions calculated in this work between the Mn atoms obtained using density functional theory are listed in the first column. The Mn-Mn distances appear in the second column.

dMn-Mn [ Å] Type Ji [K] Ji[K] J1 3
The third column gives the interaction type, either, intra or interlayer. In the fourth column, the effective interactions are given in units of K. A positive value is associated to an antiferromagnetic interaction (see text and Eq. ( 2) for the definition of the exchange interactions used in this work). In the last column we show some reported values of the corresponding interactions considering the system as having a 2D undistorted honeycomb lattice. The leading exchange interactions are in-plane and correspond to the first J1 (J1, J2) and third J3 (J7, J8) neighbors interactions of the undistorted honeycomb lattice. Interestingly, in this system, the in-plane second neighbor interactions J2 (J3, J4) is almost negligible.

is important to note that the distinction made here between the six (three pairs) of intralayer interactions is a consequence of the small distortion of the 2D honeycomb lattice. Indeed, J 1 and J 2 can be associated to a single intralayer first neighbors interaction in an undis-torted honeycomb lattice J 1 , J 3 and J 4 to a single second neighbors interaction J 2 , and J 7 and J 8 to a third neighbors interaction J 3 . Concerning the interlayer interactions, J 6 is the one strictly parallel to the c axis that can be associated to the single interlayer interaction J considered by other authors. The calculations of the exchange interactions were performed using a broken-symmetry formalism, i.e., by mapping total energies corresponding to various collinear spin arrangements within a supercell onto a Heisenberg Hamiltonian [37][38][39][40][41][42]. Different supercells were needed to distinguish the exchange interactions between two different Mn atoms which are connected by translation symmetry. For example, the 40 atoms 2×1×1 unit cell, whose total energy can be written as :

E 211 = E 0 + 5 2 2 [a 1 (J 1 + J 5 ) + a 2 (J 2 + J 7 ) + a 3 J 3 + a 4 (J 4 + J 6 ) + a 8 J 8 + 8 J 9 ]
neither allow the calculation of J 9 nor the separation of J 1 and J 5 ; J 2 and J 7 ; and J 4 and J 6 . Similar limitations arise with the 40 atoms 1×2×1 and 1×1×2 unit cells, whose total energies are :

E 121 = E 0 + 5 2 2 [b 1 J 1 + b 2 (J 2 + J 7 + 2J 8 ) + b 3 J 3 + 8 J 4 + b 5 J 5 + b 6 (J 6 + J 9 )
] and

E 112 = E 0 + 5 2 2 [c 1 (J 1 + J 5 ) + c 2 (J 2 + 2J 8 ) + c 3 J 3 + 8 J 4 + 8 J 6 + c 7 J 7 + 8 J 9 ]
The coefficients a j , b j , and c j depend on the spin arrangements of the Mn atoms. As every double supercell contains 8 Mn atoms there are a total of 256 distinct spin configurations. However, on taking crystal and spin reversal symmetries into account this number can be reduced to 22, 28, and 25 for the 2×1×1, 1×2×1, and 1×1×2 supercells respectively. The application of a least-squares minimization procedure to the three supercells allows us to calculate the nine exchange interactions. A summary of the calculated values is shown in Table I. There are four leading interactions (J 1 , J 2 , J 7 and J 8 ), which are all positive (AFM). They do not lead to frustration and they are compatible with the experimentally observed AFM order at low T [25]. As mentioned above, for an undistorted lattice they would correspond to the intralayer interactions J 1 12.7 K, J 2 0.3 K, and J 3 5.3 K which are in good agreement with reported values in the literature.

The distinction between the exchange interactions due to the small distortion will be important in the analysis of the magnetoelastic interactions because the distortion splits the energy of some two-fold degenerate Raman modes.

C. Phonon modes

To understand the reasons of the peculiar dependence with temperature of some Raman modes presented in Section IV, the 30 phonon modes of MnPS 3 have been calculated using the frozen phonon method in the 10 atoms primitive cell. The irreducible representations at the Γ-point are Γ tot = 8A g + 6A u + 7B g + 9B u , where only the 8 A g and 7 B g modes are Raman Active. The calculated 27 optical phonon frequencies that span from 46 to 557 cm -1 are shown in Fig. 6. In the lower panel we show schematically the energy of the modes, their symmetry and the relative participation weight of the three atom's types. For every mode, the corresponding bars are calculated with the normalized eigenvectors q of the standard eigenvalue problem:

M -1/2 H M -1/2 q = ω 2 q (5)
where q = M 1/2 x, M is the mass matrix, H the Hessian of the potential energy, and x the vector with the atomic displacements. As can be seen from the histogram, the high energy modes are mostly localized in the P and S atoms while the low energies ones concern vibrations of the Mn atoms (see Fig. 6).

The symmetry and energy of the calculated Raman active modes are also shown in Table II. They compare reasonably well with several reported values obtained from previous calculations [4,17]. As can be seen, the energies of some pairs of modes, separated by horizontal lines in Table II, are almost equal because they would be double degenerate E g modes in a perfect honeycomb lattice. Experimentally, it is not always easy to resolve the splitting, even though the splitting has been sucessfully resolved on feature X3 in Fig. 2(b) (see also Supplementary Material Fig. S1).

Calculations Experiments Symmetry

Energy (cm -1 ) Energy (cm 2 and Supplementary Material Fig. S1).

The correspondence between the experimental features X1 ... X8 and the calculated Raman modes is also shown in Table II. There is an overall good agreement between the corresponding energies, the discrepancy becoming slightly larger for the high energy modes. For the modes above 200 cm -1 , the energy difference increases with energy giving rise to an overall difference of 5% for the majority of the modes while the deviation of the mode around 110 cm -1 corresponds to a larger relative difference of about 9%. These relative deviations of about 5 to 10% are frequently found when comparing stateof-the-art DFT calculations with experimental measurements. It is important to note here that the approximations used in our calculations have been chosen to have a good overall description of the magnetic, electronic and structural properties of the system. A different choice of the DFT input parameters (pseudopotentials, Hubbard Hamiltonian, functional), may slightly change the calculated phonon energies and/or the distribution of deviations between the theoretical and experimental energies of different modes. The two calculated modes with irreducible representation Bg around 46 and 200 cm -1 are not visible in the experimental spectra.

VI. DISCUSSION AND SUMMARY

The two Raman modes whose frequencies increase when crossing the Néel temperature at about 120 (X1) and 155 cm -1 (X2) are associated to the two pairs of calculated Raman active modes at 109.8 (Bg) and 110.4 (Ag) for X1 and 152.5 (Bg) and 153.6 cm -1 (Ag) for X2 (see Table II). As mentioned above, these two pairs of almost degenerate modes are split in energy due to the small distortion of the honeycomb and triangular lattices. For every pair, the eigenvectors of these almost degenerate modes have very similar weights in the Mn, S, and P atoms (see lower panel in Fig. 6) and the relative displacements of the Mn and S atoms for the four modes is shown in Fig. 7. If we consider the six Mn atoms forming the hexagons in the honeycomb lattice as six consecutive Mn pairs. The effective exchange interactions of the pairs are either the first or the second neighbors intralayer interactions J 1 and J 2 respectively. These AFM interactions, which participate to stabilize the Néel AFM order at low temperature, can be associated to superexchange paths between the Mn atoms through the S atoms (see Fig. 1(b)).

One can imagine that the modification with temperature of these superexchange paths can be the mechanism at the origin of the magnetoelastic interactions responsible for the temperature dependence of the X1 and X2 modes. Above T N , due the thermal excitations, the spins localized in some Mn pairs would orient parallel to each other. As the interaction is AFM, the system would try to decrease the superexchange angle in order to decrease the interaction. The slight modification of the angles would explain the energy change of the Raman modes.

In order to visualize this effect, we have considered three virtual phonon modes : each one affecting the superexchange angle associated to two opposite exchange interactions (see Supplementary Information's Fig. S4). The upper panel in Fig. 6 shows the projections of the 27 optical modes on the three virtual phonons. The projections are siginificative for four pairs of Raman active modes and clearly larger for the two pairs of modes asscociated to peaks X1 and X2.

In conclusion, we have presented a comprehensive picture which largely accounts for the observed sensitivity of the particular phonon modes to the magnetic ordering in layered MnPS 3 antiferromagnet. More specifically, we have identified, both in the experiment and theory, that there are essentially two, phonon modes of MnPS 3 which appear at 120.2 cm -1 and 155.1 cm -1 in the low temperature Raman scattering spectra, and which show a pronounced change in the vicinity of the critical Néel temperature. The present approach can be extended to other magnetic materials including their thin layers, in reference to tracing the magnetic properties with Raman scattering response due two phonon excitations.

FIG. 1 .

 1 FIG. 1. (a) Atomic structure of MnPS3. Mn atoms are represented in violet, S in yellow, and P in gray. (b) Detail of the Mn -S bonds in the Mn-S-Mn sandwiches. The S1 sulfur atoms are associated to the super-exchange path that controls the effective magnetic interaction J2. Similarly the S2 sulfur atoms are associated to the effective magnetic interaction J1. (c) Detail of the six intralayer effective exchange interactions considered in this work. The difference between J1 and J2, J3 and J4, and J7 and J8 is a consequence of the distortion of the honeycomb lattice. (d) Detail of the three interlayer effective exchange interactions considered in this work. Atomic structures were drawn with VESTA [5].

FIG. 2 .

 2 FIG. 2. Characteristic Raman scattering spectrum of the MnPS3 crystal, measured under 514 nm laser line excitation, at low temperature (T = 5 K) in a wide spectral range 100 to 650 cm -1 from the laser line. The spectrum shown in the upper (a) panel is expanded in vertical direction in the bottom (b) panel to better visualize the less intense Ramanscattering modes. The observed narrow X1 . . . X8 Raman scattering peaks are due to phonons whereas a broad X' feature with the high energy onset at 183.3 cm -1 is of different origin.

FIG. 3

 3 FIG. 3. a) False-color map of Raman scattering response (logarithmic scale for the intensity) together with a few characteristic spectra measured as a function of temperature and focused on the spectral region encompassing the X1 and X2 phonon modes which are sensitive to the magnetic order in MnPS3 crystals. Temperature evolution of central positions of the (a) X1 and (b) X2 peaks, as extracted from fitting the shapes of these peaks with Lorentzian functions. TN denotes the critical Néel temperature for the antiferromagnetic to paramagnetic transition in MnPS3.

FIG. 4 .

 4 FIG.4. Evolution of the X2 and X' peaks with an applied magnetic field oriented perpendicular (left) and parallel (right) to the reciprocal lattice vector c * . The experiment was performed at 5 K.

FIG. 6 .

 6 FIG. 6. Upper panel: projection of the optical phonon modes in the three virtual modes shown in Fig. S4 in Supplementary Material. Lower panel: relative atomic participation in the optical phonon modes of MnPS3 of the three atom's types. The colors in the bars represent the contribution of the two Mn atoms (violet), the six S (yellow), and the two P atoms (gray). The symmetry of the modes is shown in the upper part of the figure.

FIG. 7 .

 7 FIG. 7. Schematic representation of the displacements of the Mn and S atoms in MnPS3 corresponding to the two pairs of Raman active modes associated to the experimental peaks X1 and X2. Phonon mode visualizations were drawn with VESTA [5].

FIG. S3. Temperature evolution of the amplitude of phonon-related, X1, ..X8, Raman scattering lines in MnPS3 crystals.

ACKNOWLEDGMENTS

The work has been supported by the EU Graphene Flagship project, the ATOMOPTO project (TEAM programme of the Foundation for Polish Science, co-financed by the EU within the ERDFund)

Supplementary Material: Magneto-elastic interaction in the two-dimensional magnetic material MnPS 3 studied by first principles calculations and Raman experiments Additional results concerning the analysis of the Raman scattering spectra of MnPS 3 , measured as a function of temperature are shown below. All Raman scattering lines, X1 ... X8, (for notation, see the main text) which are visible in our spectra and which have been attributed to phonon modes has been fit to Lorentzian functions and the central position, linewidth, and amplitude of these lines have been extracted. The resulting temperature evolution of these three parameters is illustrated in Figures S1,S2, and S3 for each X1 ... X8 resonance.

Figure S4 show the schematic atom displacement, i.e., phonon modes, that modifies the superexchange angles associated with the effective exchange interactions J 1 and J 2 .


Xiaoguo Xiang

Kunli Xiang

Rosa Del

C Ortiz

Florian Jabbour

Wei Wang

Integrating palaeontological and molecular data uncovers multiple ancient and recent dispersals in the pantropical Hamamelidaceae

Keywords: biogeography, divergence time estimates, fossil, Hamamelidaceae, phylogeny, tropical and subtropical forests

Aim: The integration of palaeontological and phylogenetic data can improve our understanding in the spatio-temporal evolutionary processes of living organisms. However, how best to use fossil data in divergence time estimation and ancestral range reconstruction remains challenging. Here, we integrated palaeontological and molecular data to investigate the historical biogeography of Hamamelidaceae, a pantropical angiosperm family with abundant fossils outside its present distribution.

Location: Global tropical/subtropical areas.

Methods: Using seven DNA regions (>7,500 bp) from plastid and nuclear genomes, we reconstructed a robust phylogenetic framework for Hamamelidaceae with the first complete genus-level sampling. We used the tip-dating method with the 22 fossils to estimate divergence times for the family, and inferred the ancestral range of lineages under the dispersal-extinction-cladogenesis model by incorporating the fossils.

Results: Our biogeographic analysis indicates that extant Hamamelidaceae most likely originated in tropical Asia during the mid-Cretaceous, and in the family 20 dispersals occurred during three major time intervals: the Upper Cretaceous (c. 93-69 Ma), Palaeocene-Eocene (c. 63-39 Ma) and late .

Main conclusions: Overland migrations through available land bridges and island chains may have been mainly responsible for hamamelidaceous range expansions during these three episodes. This study contributes to our knowledge on the assembly and evolution of angiosperm-dominated tropical and subtropical forests.

INTRODUCTION

A central goal of biogeography is to understand the geographic assembly of biodiversity through time [START_REF] Cox | Biogeography: An ecological and evolutionary approach[END_REF]. Two relatively independent research fields, palaeontology and phylogenetics, have been used to address this theme. It is now commonly acknowledged that a better understanding of the history of living organisms and their ecological and evolutionary processes can be accomplished through the integration of palaeontological and phylogenetic data [START_REF] Benton | Exploring macroevolution using modern and fossil data[END_REF][START_REF] Hunt | Integrating paleontological and phylogenetic approaches to macroevolution[END_REF]. The dominant role of fossil data in combination with molecular phylogenetics has been to serve as calibration points in dating analyses based on the assumed position of fossil taxa in the phylogenetic tree of extant taxa [START_REF] Parham | Best practices for justifying fossil calibrations[END_REF]. Currently, the Bayesian uncorrelated relaxed clock method [START_REF] Drummond | Relaxed phylogenetics and dating with confidence[END_REF], implemented in beast [START_REF] Drummond | BEAST: Bayesian evolutionary analysis by sampling trees[END_REF], is widely used. This method may rely on fossil data to calibrate one or several nodes of the tree and is accordingly regarded as node-dating [START_REF] Ronquist | MrBayes 3.2: Efficient Bayesian phylogenetic inference and model choice across a large model space[END_REF]. However, multiple fossil taxa associated with a particular node in a tree of extant taxa need to be reduced to a single calibration point [START_REF] Parham | Best practices for justifying fossil calibrations[END_REF][START_REF] Ronquist | MrBayes 3.2: Efficient Bayesian phylogenetic inference and model choice across a large model space[END_REF]. Another important aspect to consider when integrating phylogenetic and palaeontological data into biogeography is the reconstruction of ancestral ranges by combining information of fossil and living taxa. For a group with fossil taxa outside its modern distribution, the inclusion of fossil taxa can heavily influence its biogeographic reconstructions (Manchester & Tiffney, 2001;[START_REF] Mao | Distribution of living Cupressaceae reflects the breakup of Pangea[END_REF][START_REF] Meseguer | Integrating fossils, phylogenies, and niche models into biogeography to reveal ancient evolutionary history: The case of Hypericum (Hypericaceae)[END_REF]. Generally, fossil taxa can be directly placed in the phylogenetic trees of extant taxa based on the appearance time and can thereby be included in inferences of ancestral ranges (e.g. [START_REF] Clayton | Recent long-distance dispersal overshadows ancient biogeographical patterns in a pantropical angiosperm family (Simaroubaceae, Sapindales)[END_REF][START_REF] Mao | Distribution of living Cupressaceae reflects the breakup of Pangea[END_REF][START_REF] Zhang | Integrated fossil and molecular data reveal the biogeographic diversification of the eastern Asian-eastern North American disjunct hickory genus (Carya Nutt.)[END_REF]. However, this strategy can sometimes be problematic. For example, if two or more fossils belonging to a single clade are contemporaneous, based only on the appearance time, their placements will result in a polytomy for the lower clade (zero-length branch).

The recently developed tip-dating (TD) method estimates divergence times with all available fossil taxa as terminal tips based on morphological or combined morphological and molecular data [START_REF] Pyron | Divergence time estimation using fossils as terminal taxa and the origins of Lissamphibia[END_REF][START_REF] Ronquist | MrBayes 3.2: Efficient Bayesian phylogenetic inference and model choice across a large model space[END_REF]. This method can thereby generate a dated phylogeny containing extinct and extant taxa, or even entirely extinct taxa and overcome zero-length branches to a certain extent. Moreover, in comparison with node dating, TD method may establish direct calibrations using the known age of the fossils, and accommodates the uncertainty of their phylogenetic position through the co-estimation of time and topology [START_REF] O'reilly | Dating tips for divergence-time estimation[END_REF]. Therefore, TD can avoid some arbitrary settings of prior calibration points, and yield more accurate age estimates while providing a coherent measure of statistical uncertainty [START_REF] Heath | The fossilized birth-death process for coherent calibration of divergence-time estimates[END_REF].

Tropical and subtropical forests are well-known as the most species-rich and productive terrestrial biomes on Earth. How the past biogeographic events have shaped their present distribution pattern, diversification and endemism has long been the quest in biology and ecology [START_REF] Emerson | Species diversity can drive speciation[END_REF]. Several studies focused on the temporal patterns of the origin and evolution of angiosperm-dominated forests (e.g. [START_REF] Couvreur | Origin and global diversification patterns of tropical rain forests: Inferences from a complete genus-level phylogeny of palms[END_REF][START_REF] Davis | Explosive radiation of Malpighiales supports a mid-Cretaceous origin of modern tropical rain forests[END_REF][START_REF] Wang | Rosid radiation and the rapid rise of angiosperm-dominated forests[END_REF][START_REF] Yu | Insights into the historical assembly of East Asian subtropical evergreen broadleaved forests revealed by the temporal history of the tea family[END_REF]. A few studies have investigated spatial patterns of tropical angiosperm lineages in the early Tertiary, a warming period (e.g. [START_REF] Bartish | Vicariance or long-distance dispersal: Historical biogeography of the pantropical subfamily Chrysophylloideae (Sapotaceae)[END_REF][START_REF] Davis | Laurasian migration explains Gondwanan disjunctions: Evidence from Malpighiaceae[END_REF][START_REF] Lian | A multi-locus phylogeny for the Neotropical Anomospermeae (Menispermaceae): Implications for taxonomy and biogeography[END_REF]. Evidence from fossil and molecular data suggests that angiosperm-dominated tropical forests arose in the mid-Cretaceous (e.g. [START_REF] Couvreur | Origin and global diversification patterns of tropical rain forests: Inferences from a complete genus-level phylogeny of palms[END_REF][START_REF] Davis | Explosive radiation of Malpighiales supports a mid-Cretaceous origin of modern tropical rain forests[END_REF][START_REF] Morley | Cretaceous and Tertiary climate change and the past distribution of megathermal rainforests[END_REF][START_REF] Upchurch | Cretaceous vegetation of the Western Interior and adjacent regions of North America[END_REF][START_REF] Wang | Rosid radiation and the rapid rise of angiosperm-dominated forests[END_REF]. To better understand spatio-temporal patterns and processes underlying the distribution of biodiversity in tropical forests, one needs to investigate biogeographic events of an ancient group on a worldwide basis.

In this study, we focus on the historical biogeography of the angiosperm family Hamamelidaceae. During the past years, molecular phylogenetic studies have re-delimitated this family, where the subfamily Altingioideae was excluded and raised to the familial rank as Altingiaceae (e.g. [START_REF] Ickert-Bond | A taxonomic synopsis of Altingiaceae with nine new combinations[END_REF][START_REF] Jian | Resolving an ancient, rapid radiation in Saxifragales[END_REF][START_REF] Stevens | Angiosperm Phylogeny Website. Version 14[END_REF]. As currently circumscribed, Hamamelidaceae contains 28 genera and c. 120 species [START_REF] Judd | Plant systematics: a phylogenetic approach[END_REF], which are small shrubs or trees of great ecological importance inhabiting tropical and subtropical forests worldwide, with Hamamelis extending to temperate regions [START_REF] Ohsawa | Structural comparison of tropical montane rain forest along latitudinal and altitudinal gradients in south and east Asia[END_REF][START_REF] Wu | The families and genera of angiosperms in China[END_REF][START_REF] Ying | Plant geography of China[END_REF]; Figure S1). Hamamelidaceae is well-known for its broad and scattered geographic distribution and its endemics [START_REF] Endress | The families and genera of vascular plants[END_REF][START_REF] Zhang | Hamamelidaceae: Geographic distribution, fossil history and origin[END_REF]. The majority of genera are geographically restricted except for Hamamelis which shows a disjunct distribution between eastern North America and eastern Asia (Figure 1). The fruits of Hamamelidaceae are woody dehiscent capsules [START_REF] Endress | The families and genera of vascular plants[END_REF]. Seeds are usually dispersed by a ballistic mechanism, and hence are not expected to disperse across long distances or oceanic barriers even though seeds of Exbucklandia and Rhodoleia have small narrow wings [START_REF] Endress | The families and genera of vascular plants[END_REF]. The oldest fossil flowers of the family come from the late Santonian of North America and Europe (c. 84 Ma; [START_REF] Endress | Archamamelis, hamamelidalean flowers from the Upper Cretaceous of Sweden[END_REF]. Recent molecular clock estimates suggest a stem age of c. 84-100 Ma for the family (Jian et al., 2008) or 95.93 Ma (90.99-104.04;[START_REF] Magallón | A meta-calibrated timetree documents the early rise of flowering plant phylogenetic diversity[END_REF]. Fossil evidence indicates that Hamamelidaceae were widely distributed in the Northern Hemisphere in the Upper Cretaceous (Figure 1; [START_REF] Zhang | Hamamelidaceae: Geographic distribution, fossil history and origin[END_REF]. The family has an abundant fossil record in the early Tertiary, with representation of many extant and extinct genera (reviewed by [START_REF] Manchester | Eastern Asian endemic seed plant genera and their paleogeographical history throughout the Northern Hemisphere[END_REF][START_REF] Zhang | Hamamelidaceae: Geographic distribution, fossil history and origin[END_REF]. Importantly, most of the fossils are distributed outside the modern ranges of the taxa to which they have been assigned (Figure 1). Thus, Hamamelidaceae provides an opportunity to investigate historical biogeography of living organisms through the integration of palaeontological and molecular phylogenetic data.

Here, we first infer the phylogeny of Hamamelidaceae with complete genus-level sampling. With the resulting phylogenetic framework, we then estimate divergence times for the family using the TD method with 22 fossils and investigate the historical biogeography of the family by incorporating the fossils. Given the Cretaceous age and widespread distribution of Hamamelidaceae, a weak seed dispersal ability, the restriction of most genera to particular tropical regions, as well as an abundant fossil record, this study may provide new insights into the past tropical floristic exchanges among different continents.

MATERIALS AND METHODS

Taxon sampling and DNA sequencing

We sampled 55 taxa of Hamamelidaceae, representing all 28 of the currently recognized extant genera of the family; the genus Embolanthera, endemic to southeastern Asia, was included in a molecular phylogenetic analysis for the first time. Based on [START_REF] Jian | Resolving an ancient, rapid radiation in Saxifragales[END_REF] and [START_REF] Soltis | Angiosperm phylogeny: 17 genes, 640 taxa[END_REF], we selected 14 species from close relatives of Hamamelidaceae as outgroups: three species from its sister clade Daphniphyllaceae-Cercidiphyllaceae, 10 species from Altingiaceae and one species of Paeoniaceae. Six plastid (atpB, matK, rbcL, atpB-rbcL, trnH-psbA and trnLUAA-FGAA) and one nuclear (ITS) DNA region were used in this study. Voucher information and GenBank accession numbers are listed in Table S2. Total genomic DNA was extracted from fresh silica gel-dried leaves or herbarium specimens using DNeasy Mini Plant Kits (Tiangen Biotech). The selected DNA regions were amplified with standard polymerase chain reaction (PCR). The primers used in this study are listed in Table S3. Sequencing reactions were conducted using the ABI Prism BigDye Terminator Cycle Sequencing Kit (Applied Biosystems, ABI). Products were analysed on ABI 3730xl automated DNA sequencers.

Alignment and phylogenetic analysis

Sequences were aligned using the default parameters in Clustal X 1.83 [START_REF] Thompson | The Clustal X windows interface: Flexible strategies for multiple sequence alignment aided by quality analysis tools[END_REF]) and manually adjusted with BioEdit 5.0.9 [START_REF] Hall | BioEdit: A user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT[END_REF]. After excluding ambiguous regions, we obtained a concatenated alignment of 7,549 bp for 69 taxa (Table S4). Phylogenetic analyses were carried out using maximum parsimony (MP) and Bayesian inference (BI) methods in paup 4.0b10 [START_REF] Swofford | PAUP*: Phylogenetic analysis using parsimony (*and other methods)[END_REF] and MrBayes 3.2.6 [START_REF] Ronquist | MrBayes 3.2: Efficient Bayesian phylogenetic inference and model choice across a large model space[END_REF] respectively. For MP analyses, heuristic searches were performed with 1,000 random sequence addition replicates, tree-bisection-reconnection (TBR) branch swapping, MulTrees in effect and steepest descent off. Gaps were treated as missing data, characters were equally weighted and their states were unordered. Internal branch support was estimated by using 1,000 bootstrap replicates as described above. For BI analyses, each DNA region was assigned its own model of nucleotide substitution (Table S4), as determined by the Akaike information criterion (AIC) in Modeltest 3.06 [START_REF] Posada | MODELTEST: Testing the model of DNA substitution[END_REF]. Four chains of the Markov chain Monte Carlo (MCMC) were run, sampling one tree every 1,000 generations for 3,000,000, starting with a random tree. Majority-rule (> 50%) consensus trees were constructed after removing the burn-in period samples (the first 25% of sampled trees).

Divergence time estimation

Data about the fossil records of Hamamelidaceae were gathered from the literature and from the Paleobiology Database (http://paleo db.org/). We selected only the fossils that could be confidently attributed to specific taxonomic groups. A total of 19 hamamelidaceous fossils were used to estimate divergence times and reconstruct biogeographic scenarios. One fossil from Cercidiphyllaceae and two fossils from Altingiaceae were also included. Justifications of the ages and placements of the 22 fossils are given in Table S1. A likelihood ratio test was first carried out to estimate whether Hamamelidaceae evolved at a homogeneous rate along all branches of the phylogeny. The result rejected a constant rate for the family (δ = 1503.14, df = 67, p < 0.0001). Thus, we estimated divergence times of Hamamelidaceae using TD method with a Bayesian relaxedclock model. We first constructed a matrix containing 69 extant and 22 fossil taxa. The molecular data for fossils are coded as missing, and in total 68 morphological characters for 91 taxa were included (Table S5). The majority-rule consensus tree from the BI analysis was used and fossil taxa were placed in the relevant clades based on apomorphy-based or phylogenetic methods [START_REF] Sauquet | Testing the impact of calibration on molecular divergence times using a fossil-rich group: The case of Nothofagus (Fagales)[END_REF]; Table S1). Divergence time estimation was performed in MrBayes 3.2.6 under a fossilized birth-death (FBD) prior process [START_REF] Heath | The fossilized birth-death process for coherent calibration of divergence-time estimates[END_REF], with a lognormal prior and an uncorrelated independent gamma rates model. The calibrated ages of the 22 internal fossils were assigned a uniform prior distribution with ranges corresponding to the age ranges of their respective geological strata. To avoid overestimation of root age, we set a maximum age of 112.99 Ma (sd = 0.5) for the root, which is the crown group age of Saxifragales [START_REF] Magallón | A meta-calibrated timetree documents the early rise of flowering plant phylogenetic diversity[END_REF]. Morphological data were analysed under the Mk model [START_REF] Lewis | A likelihood approach to estimating phylogeny from discrete morphological character data[END_REF]. For the molecular partitions, each DNA region was assigned an unlinked substitution model (Table S4). Two independent MCMC runs were carried out, every run with 100 million generations and sampling trees every 5,000 generations. Convergence was evaluated as described in the beast analyses, and in each MCMC run, the first 10 million generations of the samples were discarded as burn-in.

Biogeographic analysis

Ancestral area reconstructions (AARs) were conducted using a recently developed Statistical BioGeoBEARS (S-BioGeoBEARS), implemented in rasp 4.0 [START_REF] Matzke | Probabilistic historical biogeography: New models for founder-event speciation, imperfect detection, and fossils allow improved accuracy and model-testing[END_REF][START_REF] Yu | RASP (Reconstruct Ancestral State in Phylogenies): A tool for historical biogeography[END_REF]. [START_REF] Ree | Conceptual and statistical problems with the DEC+J model of founder-event speciation and its comparison with DEC via model selection[END_REF] suggest that the biogeographic models with +J parameter (founder-event speciation) have conceptual and statistical problems. The dispersal-extinction-cladogenesis model (DEC; [START_REF] Ree | Maximum likelihood inference of geographic range evolution by dispersal, local extinction, and cladogenesis[END_REF] was used here. Following [START_REF] Buerki | An evaluation of new parsimony based versus parametric inference methods in biogeography: A case study using the globally distributed plant family Sapindaceae[END_REF], we specified dispersal probabilities between pairs of areas for four separate time slices (Table S6). According to current and fossil distributions of Hamamelidaceae, species sampled were assigned to one or more of the following seven biogeographic regions: tropical Asia, East Palaearctic, Oceania, Africa, Europe (including western Asia), North America and Central America (including southern Mexico and northern South America) (Figure 1). We randomly sampled 1,000 timetrees generated from MrBayes as a 'tree file' and used the majority-rule consensus tree derived from the MCMC stationary sample as a final representative tree.

RESULTS

Phylogeny

MP and BI analyses resulted in highly congruent trees (Figure S2). The phylogeny we inferred for Hamamelidaceae is the most taxonomically complete to date, and includes all genera recognized at present. Hamamelidaceae and each of its four subfamilies are all strongly supported as monophyletic. Exbucklandioideae is the earliest-diverging lineage (BS = 100%, PP = 1.00), followed by Mytilarioideae (BS = 100%, PP = 1.00). Within the family, relationships among the subfamilies are consistent with the results of [START_REF] Li | Molecular phylogenetics of Hamamelidaceae: Evidence from DNA sequences of nuclear and chloroplast genomes[END_REF], but are usually resolved with greater support for clades found therein. Within Hamamelidoideae, each of Fothergilleae, Hamamelideae, Dicorypheae and Eustigmateae is strongly supported as monophyletic, while Loropetaleae is paraphyletic, with Embolanthera recovered as sister to Corylopsideae (BS = 54%, PP = 1.00). The majority of intergeneric relationships in the family are well resolved.

Divergence times

The dated phylogenetic tree is presented in Figure S3.

Biogeographic analyses

The DEC model estimated the dispersal rate (d) as 0.0008, and the extinction rate (e) as 0.0002. The distribution of Hamamelidaceae is inferred to be the result of 20 dispersal events (Figure 2): six times from tropical Asia to North America (nodes 5, 29, 31, 33, 52, 60), four of which led to extinct branches; four from tropical Asia to Europe (nodes 7, 12, 56, 61), all of which led to extinct lineages; one from tropical Asia to Central America (node 20); one from tropical Asia to Africa (node 18); one from Africa to Oceania (node 26); one from tropical Asia to East Palaearctic (node 63); one from Europe to Central America (node 57); two from North America to East Palaearctic (nodes 37, 40), one of which led to an extinct branch (node 40); one from North America to Europe (node 30); one from East Palaearctic to tropical Asia (node 44); and one from East Palaearctic to Europe (node 45). The MRCA of Hamamelidaceae as well as of each of its four subfamilies were inferred to be distributed in tropical Asia.

DISCUSSION

Our inferred age estimation of the hamamelidaceous stem group ranges between 108.43 and 113.68 Ma, which is much older than previous estimates (84-100 Ma in [START_REF] Jian | Resolving an ancient, rapid radiation in Saxifragales[END_REF]95.93 Ma in Magallón et al., 2015). Previous studies had a limited taxon sampling for Hamamelidaceae and/or no hamamelidaceous fossil calibrations were used. Similar older ages were produced for Fagales [START_REF] Xiang | Large-scale phylogenetic analyses reveal fagalean diversification promoted by the interplay of diaspores and environments in the Paleogene[END_REF] and Ranunculaceae [START_REF] Wang | The rise of angiosperm-dominated herbaceous floras: Insights from Ranunculaceae[END_REF] when more taxa and fossils were included. Internal calibrations and dense taxon sampling have been regarded to be vital to obtain accurate time estimates [START_REF] Jian | Resolving an ancient, rapid radiation in Saxifragales[END_REF][START_REF] Yang | Bayesian estimation of species divergence times under a molecular clock using multiple fossil calibrations with soft bounds[END_REF]. Our analyses indicate that the split of Hamamelidaceae and its sister group occurred at 108.43 Ma (95% HPD: 104.08-111.98; node 1) and the MRCA of Hamamelidaceae was likely present in tropical Asia by the mid-Cretaceous (node 2; Figure 2). This result supports the hypothesis of [START_REF] Takhtajan | Flowering plants: Origin and dispersal[END_REF] that this family originated in the southeastern Asian continent, rather than the eastern Asia-western North America origin as hypothesized by [START_REF] Wolfe | Fossil forms of Amentiferae[END_REF].

Our biogeographic reconstruction shows that at least 20 dispersals need to be invoked to explain the current distribution of Hamamelidaceae and that of the fossils attributed to it. Despite the uncertainty involved in any estimate using a molecular clock, based on our estimated median ages, these 20 dispersals occurred during three major episodes: five in the Upper Cretaceous (c. 93-69 Ma), nine in the Palaeocene-Eocene (c. 63-39 Ma) and six in the late Oligocene (c. 27-23 Ma) (Figure 3). Among the five migration events inferred in Hamamelidaceae during the Upper Cretaceous, two occurred from tropical Asia to North America (node 29, 52), and one from North America to East Palaearctic (node 37). Palynological data suggest that during the Upper Cretaceous the Northern Hemisphere could be divided into two distinctive floras with different climates: the Normapolles flora (including eastern North America and Europe) and the Aquilapollenites flora (including major parts of Asia and western North America) [START_REF] Wolfe | Some aspects of the plant geography of the Northern Hemisphere during the Late Cretaceous and Tertiary[END_REF][START_REF] Zhou | The Cretaceous and Tertiary Aquilapolles in China[END_REF]. It has been documented that the Turgai Strait might have been a barrier to plant migration between Asia and Europe by the middle Upper Cretaceous [START_REF] Lomolino | Biogeography[END_REF]. On the other hand, [START_REF] Graham | The role of land bridges, ancient environments, and migrations in the assembly of the North American Flora[END_REF] has suggested that the exchange of plant populations between North America and Asia might have occurred via the Bering land bridge (BLB) in the Upper Cretaceous. Hamamelidaceous remains have been found in Siberia and the Amuro-Zeya Depression of North Asia during the Turonian [START_REF] Maslova | Infructescences of Kasicarpa gen. nov. (Hamamelidales) from the late Cretaceous (Turonian) of the Chulym-Yenisey depression, western Siberia, Russia[END_REF]. Thus, the BLB may have played an important role in biotic exchange between Asia and North America for Hamamelidaceae during the Upper Cretaceous. Our results indicate that one migration occurred from North America to Europe (node 30) during the same period. Both of these two areas belong to the Normapolles flora [START_REF] Wolfe | Some aspects of the plant geography of the Northern Hemisphere during the Late Cretaceous and Tertiary[END_REF], and the migration route across the North Atlantic was supported by the Cretaceous fossil records of Hamamelidaceae (Figure 1). In addition, one migration occurred from tropical Asia to Europe (node 7) in the late Upper Cretaceous, perhaps facilitated by the land bridge along the closing Tethys seaway that connected southern Europe and southwestern and southeastern Asia [START_REF] Scotese | A continental drift flipbook[END_REF]Tiffney & Manchester, 2001). The fossil record of Hamamelidaceae from southwestern Europe supports this route [START_REF] Endress | Archamamelis, hamamelidalean flowers from the Upper Cretaceous of Sweden[END_REF].

Nine migration events were inferred to have occurred during the Palaeocene-Eocene. Among them, three are from tropical Asia to Europe (nodes 12, 56, 61), two are from tropical Asia to North America (nodes 31, 60) and one is from East Palaearctic to Europe (node 45). The Asian precursors might have gradually migrated along the island chain of the Tethys seaway into southern Europe [START_REF] Heaney | Conservation biogeography in oceanic archipelagoes[END_REF]. This migration route was also suggested for Alangium sect. Marlea (Alangiaceae) during the similar period [START_REF] Feng | Phylogeny and biogeography of Alangiaceae (Cornales) inferred from DNA sequences, morphology, and fossils[END_REF]. Hamamelidaceae expanded further westwards to North and Central America through the North Atlantic Land Bridge (NALB), which extended south to roughly 45° and could have maintained sufficient warmth, light and humidity for thermophilic broadleaved evergreen plants to flourish (Tiffney & Manchester, 2001). In the early Palaeocene, one lineage migrated from Asia to Africa (node 18). Many 'stepping-stone' chains in the Tethys Sea [START_REF] Scotese | A continental drift flipbook[END_REF] could have facilitated floristic exchange between Laurasia and Africa. A similar migration route from southeastern Asia into Africa during the Palaeocene was suggested for Alangium [START_REF] Feng | Phylogeny and biogeography of Alangiaceae (Cornales) inferred from DNA sequences, morphology, and fossils[END_REF] and Cucurbitaceae [START_REF] Schaefer | Gourds afloat: A dated phylogeny reveals an Asian origin of the gourd family (Cucurbitaceae) and numerous oversea dispersal events[END_REF]. Within Asia, one migration is also inferred from East Palaearctic to tropical Asia (node 44). The Palaeocene-Eocene largely coincides with the time of dominance of the boreotropical flora [START_REF] Wolfe | Some aspects of the plant geography of the Northern Hemisphere during the Late Cretaceous and Tertiary[END_REF], which may have facilitated exchanges of thermophilic plants in the Northern Hemisphere.

The divergence time of the African-Oceanian disjunction in Hamamelidaceae was estimated to have occurred at c. 39.83 Ma (node 26), which is long after the fragmentation of eastern Gondwana (c. 95-84 Ma;[START_REF] Mcloughlin | The breakup history of Gondwana and its impact on pre-Cenozoic floristic provincialism[END_REF]. A similar African-Oceanian disjunction was reported in Monimiaceae [START_REF] Renner | Biogeography of the Monimiaceae (Laurales): A role for East Gondwana and long-distance dispersal, but not West Gondwana[END_REF] and Adansonia (Bombacaceae) [START_REF] Wickens | The baobab: Africa's upside-down tree[END_REF]. The seeds of Hamamelidaceous plants with a ballistic dispersal mechanism could be dispersed over only short distances. However, the very hard exocarps may be able to maintain the viability of the seed if they are dispersed over longer distances by water [START_REF] Endress | The families and genera of vascular plants[END_REF]. Therefore, we hypothesize that Hamamelidaceae could have migrated from Africa to Oceania by transoceanic dispersal in the late Eocene.

Six migration events occurred during the late Oligocene. Among these six events, two are from tropical Asia to North America (nodes 5, 33) and one is from North America to East Palaearctic (node 40). Two similar and contemporaneous biotic exchanges were also reported in Malpighiaceae [START_REF] Davis | High-latitude Tertiary migrations of an exclusively tropical clade: Evidence from Malpighiaceae[END_REF]. In addition, we inferred one migration from tropical Asia to Central America (node 20) and one from Europe to Central America (node 57) during this period. In Ampelopsis (Vitaceae), [START_REF] Nie | Evolution of the intercontinental disjunctions in six continents in the Ampelopsis clade of the grape family (Vitaceae)[END_REF] found that migration from North America into both Europe and Asia took place during the early Miocene, and suggested a role of the NALB in enabling those migrations. Evidence from geology and palaeontology also support the idea that the NALB might have acted as a dispersal corridor for tropical or warm temperate floras in the early Eocene to middle Miocene [START_REF] Tiffney | Perspectives on the origin of the floristic similarity between eastern Asia and eastern North America[END_REF]. Thermophilic broadleaved evergreen elements, such as Hamamelidaceae, Malpighiaceae and Vitaceae, could thus have migrated between Eurasia and tropical America through the NALB during the late Oligocene and early Miocene. Within Asia, one migration is also inferred from tropical Asia to East Palaearctic (node 63). In the late Oligocene, global climate became warming [START_REF] Zachos | An early Cenozoic perspective on greenhouse warming and carbon-cycle dynamics[END_REF], which may have driven range expansion of thermophilic lineages.

In conclusion, our integration of fossil and phylogenetic data indicate that extant Hamamelidaceae likely extend back to the mid-Cretaceous of tropical Asia and that shortly afterwards four major subfamilies diversified in tropical Asia about 105-97 Ma. A similar diversification pattern is found in the early evolutionary history of pantropical Arecaceae [START_REF] Couvreur | Origin and global diversification patterns of tropical rain forests: Inferences from a complete genus-level phylogeny of palms[END_REF] and Menispermaceae [START_REF] Wang | Menispermaceae and the diversification of tropical rainforests near the Cretaceous-Paleogene boundary[END_REF], both of which are important components of modern tropical forests. Thus, our data for Hamamelidaceae suggest that Laurasia, especially tropical Asia, is not only a key source region, but also an important area where some plant families currently inhabiting tropical and subtropical forests might have diversified. Our findings further indicate that the range expansions of Hamamelidaceae occurred during three major episodes: the Upper Cretaceous, Palaeocene-Eocene and late Oligocene. Several key angiosperm elements of tropical and subtropical forests, such as Annonaceae [START_REF] Thomas | The historical origins of palaeotropical intercontinental disjunctions in the pantropical flowering plant family Annonaceae[END_REF], Malpighiaceae [START_REF] Davis | High-latitude Tertiary migrations of an exclusively tropical clade: Evidence from Malpighiaceae[END_REF] and Menispermaceae [START_REF] Wang | Menispermaceae and the diversification of tropical rainforests near the Cretaceous-Paleogene boundary[END_REF], also appear to have expanded their ranges during one or two of these three episodes. Thus, we assume that these three episodes might also have been vital for the evolution and assembly of angiosperm-dominated tropical and subtropical forests. Yet, Hamamelidaceae is only one of many important components of angiosperm-dominated tropical and subtropical forests. This hypothesis needs to be further tested by studying other anciently thermophilic angiosperm groups in a broad phylogenetic context.

FIGURE 1

 1 FIGURE 1 Geographic distribution of both extant and extinct hamamelidaceous species. The extant members of widespread and endemic genera (the latter in brackets) are indicated for each region. The photograph inserts are Fothergilla major (A), Parrotia persica (B), Corylopsis spicata (C), Molinadendron sinaloense (D), Trichocladus crinitus (E) and Noahdendron nicholasii (F). The fossil distributions are referred in TableS1. Photographs by S.X.Y

FIGURE 3

 3 FIGURE 3 Age estimates for dispersal events in Hamamelidaceae. The inferred dispersal events between areas summarized from Figure2are indicated by colours (same colour coded as in Figure1). Black dots and horizontal bars represent the median ages and 95% highest posterior density intervals respectively. Node numbers are referred as those in Figure2. " †" indicates a fossil taxon. The three grey areas were hypothesized as important intervals for range expansion of Hamamelidaceae.

ACKNOWLEDGEMENTS

We sincerely thank Peter F. Stevens for kindly reading an earlier version of this manuscript, Shengxiang Yu for providing plant photos and Anmin Lu for his encouragement and valuable suggestions. We also thank the associate editor Richard Ree and four anonymous referees for their comments and invaluable suggestions that greatly improved our manuscript. This research was supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (XDB31030000) and National Natural Science Foundation of China (31470315, 31770231, 31300181, 31770233, and 31670212).

Funding information

Strategic Priority Research Program of the Chinese Academy of Sciences, Grant/Award Number: XDB31030000; National Natural Science Foundation of China, Grant/Award Number: 31470315, 31770231, 31300181, 31770233, and 31670212

CONFLICT OF INTEREST

The authors declared that they have no conflicts of interest in this work.

DATA AVAILABILITY STATEMENT

The sequences reported in this paper have been deposited in the GenBank database (accession nos. MG846141-MG846212).