Decoupling Functional and Morphological Convergence, the Study Case of Fossorial Mammalia

Morphological similarity between biological structures in phylogenetically distant species is usually regarded as evidence of convergent evolution. Yet, phenotypic similarity is not always a sign of natural selection acting on a particular trait, therefore adaptation to similar conditions may fail to generate convergent lineages. Herein we tested whether convergent evolution occurred in the humerus of fossorial mammals, one of the most derived biological structures among mammals. Clades adapting to digging kinematics possess unusual, by mammalian standards, humeral shapes. The application of a new, computationally fast morphological test revealed a single significant instance of convergence pertaining to the Japanese fossorial moles (Mogera) and the North-American fossorial moles (Scalopini). Yet, the pattern only manifests when trade-off performance data (derived from finite element analysis) are added to shape data. This result indicates that fossorial mammals have found multiple solutions to the same adaptive challenge, independently moving around multiple adaptive peaks. This study suggests the importance of accounting for functional trade-off measures when studying morpho-functional convergence. We revealed that fossorial mammals, a classic example of convergent evolution, evolved multiple strategies to exploit the subterranean ecotope, characterized by different functional trade-offs rather than converging toward a single adaptive optimum.


INTRODUCTION
Independent evolution of similar morphological adaptations to similar environmental challenges is commonly regarded as evidence for convergence (Harmon et al., 2005;Losos, 2011) and convergent evolution underpins many hypotheses and processes that explain the repeated evolution of similar forms through time (Tseng and Flynn, 2018). However, species can adapt to similar environmental pressures in different ways Losos, 2011). One-to-many mapping, and functional trade-offs will drive the phenotypic evolution of a structure according to which traits affect fitness the most (Polly et al., 2016;Dickson and Pierce, 2019;Moen, 2019), and the evolution of a particular phenotype may result in the reduced fitness of other attributes (Alfaro et al., 2004;Moen, 2019). In addition, many different phenotypes can produce highly similar functional outputs, known as many-to-one mapping, and divergent morphologies may perform convergently in function Wainwright et al., 2005;Losos, 2011;Renaud et al., 2018), meaning that functional convergence can be recognized even in the absence of morphological convergence. In this case, species distribute around different adaptive peaks, precluding a direct relationship between the traits and the environment. Overall, both trade-offs and many-to-one mapping can provide insights into how organisms evolve under similar environmental pressures and can necessitate reassessment of the way we measure or recognize convergent evolution (Moen, 2019).
Fossorial mammals have often been reported as a textbook example of convergent evolution (Nevo, 1979;Hildebrand et al., 1985). The subterranean environment is relatively predictable in time, and requires a high degree of specialization, which in turn can lead to comparable morphological forms (Nevo, 1979;Marcy et al., 2016;Sansalone et al., 2019). Moreover, locomotion underground requires extreme morphological specialization and digging kinematic is likely to have a strong impact on other functional traits (Piras et al., 2012(Piras et al., , 2015Sansalone et al., 2018). For instance, the overall anterior limb mobility of the forelimb is sacrificed in favor of maximizing the abduction moment (Gambaryan et al., 2002). Here we will measure the extent of functional and morphological convergence among ecologically convergent fossorial mammals ranging in degree of fossorial specialization (cursorial, semi-fossorial, fossorial) and from across Mammalia (monotremes, marsupials and placentals) including a number of extinct forms. When comparing distantly related groups, the presence of similar morphological features could be interpreted as the result of phylogenetic relatedness (e.g., parallelisms) rather than convergence. For this reason, it is highly recommended to consider the phylogenetic history of the species examined. Here, we have combined 3D geometric morphometrics (GMM), finite elements analysis (FEA) and novel phylogenetic comparative methods  to quantify the extent of functional and morphological convergence in the humeri of fossorial mammals.
We managed to include four living species of Chrysochloridae, the well-preserved North-American proscalopid Mesoscalops montanensis, the extinct marsupial notoryctid Naraboryctes philcreaseri and the monotreme Tachyglossus aculeatus. Supplementary Table S1 summarizes details about all the specimens included. All specimens, except for N. philcreaseri, Chrysochloris sthulmanni, T. aculeatus and were scanned at the "Studio dentistico Moscato", Rome, Italy using a Kodak 9000 3D-Cone Beam Computed Tomography (CBCT) scanner. C. sthulmanni CT are available from the digimorph database 1 . M. montanensis has been scanned using a Skyscan 1076 microCT at the Small ANimal Tomographic Analysis Facility (SANTA) at the Seattle Children's Research Institute, Seattle, Washington. The humeri of N. philcreaseri and T. aculeatus were scanned using a Siemens Inveon MicroPET-CT scanner at the University of New South Wales. Each specimen scan comprised 1000 slices and the images were processed and stacked using GE Phoenix proprietary software.
The volume models possessed between 1.2 million and 1.4 million tetrahedral elements and between 317,549, and 401,341 nodes (i.e., tetrahedral vertices). The reconstructed volume meshes were then imported into Strand7 (v. 3.0) where finite element models were generated for each specimen. We used FEA in a comparative fashion, therefore we applied the same approximations to all the meshes. Material properties were assigned for cortical bone with a Young modulus of 10 GPa and Poisson ratio of 0.41 and homogeneously applied through the structure (Rayfield et al., 2001). According to published protocols (Stayton, 2009), simulations were conducted after scaling all models to a unit volume: 0.82 cm 3 (Piras et al., 2015).
Following Piras et al. (2015), we identified the different functional types occurring in our sample. For talpids, chrysochlorids and proscalopids we replicated the modeling approach presented in Piras et al. (2015) simulating the action of the main muscles involved in digging for each functional group. For notoryctids and tachyglossids, we utilized detailed anatomical descriptions of both myology and osteology in order to identify the main muscle groups involved in the digging kinematics (Jenkins, 1970;Warburton, 2003Warburton, , 2006Archer et al., 2010;Beck et al., 2016;Regnault and Pierce, 2018). Due to the comparative nature of this study, we applied the same resultant force to all the models. For this purpose, we used Talpa europaea as reference model since in vivo muscle force estimation is available for this species only (Gambaryan et al., 2002). We applied the same resultant force to the trochlear area for all the models following the methodology proposed by Piras et al. (2015).

Restraints
We applied on the models surface, over the selected restraint areas, an elastic solid mesh composed of beam elements (Structural Steel AS 4100-1998, Young's Modulus of 200,000 MPa). This strategy mitigates the impact of artifacts that can occur where single nodes are loaded (Attard et al., 2016;Tsang et al., 2019).
We applied the restraints on muscle insertion areas for each taxon. The restraints can perform along the axes of the reference system and were modeled as spring foundations. Talpids anatomical restraints were defined on the head of humerus and the clavicular articular facet. Chrysochlorids, proscalopids, notoryctids and tachyglossids restraints were identified only on the humeral head because these latter taxa either do not have a humerus clavicle articulation or it is extremely reduced (Jenkins, 1970;Puttick and Jarvis, 1977;Barnosky, 1981;Gasc et al., 1986;Warburton, 2006;Beck et al., 2016).

Loadings
As with restraints, the loadings have been applied on different selected areas according to the modeled muscles for the different taxa. For talpids, the loadings were defined on the areas of insertion of the m. teres major and m. pectoralis pars sternalis, two of the most important muscles used for burrowing (Dobson, 1882;Freeman, 1886;Edwards, 1937;Gambaryan et al., 2002). M. montanensis was loaded in the areas where m. teres major and either m. infraspinatus or m. spinodeltoideus insert (following Barnosky, 1981, Figure 28, p. 326;Piras et al., 2015; see also Supplementary Figure S1). For chrysochlorids, we marked the insertion areas of m. latissimus dorsi and m. triceps (Puttick and Jarvis, 1977;Gasc et al., 1986). For notoryctids we selected the areas corresponding to insertion of m. triceps, m. pectoralis and areas of the medial epicondyle corresponding to the attachment of muscles aimed at wrist and digit flexion (Warburton, 2006;Beck et al., 2016). For tachyglossids we selected the areas where m. latissimus dorsi and m. teres major insert (Regnault and Pierce, 2018). Orientation of muscles has been based on anatomical descriptions available (Jenkins, 1970;Puttick andJarvis, 1977, Gasc et al., 1986;Gambaryan et al., 2002;Warburton, 2006; see Figures 1C,D).

Stress Evaluation
We integrated the results from FEA with GMM outputs to provide a more direct comparison of the results. At this purpose we followed the protocol published by Parr et al. (2012). The mean von Mises (VM) stress values were extracted from the four nearest tetrahedral elements (each tetrahedral element has four nodes) at the point were the 38 landmarks have been placed during the digitization process (see GMM section).

Forearm Mobility
In fossorial mammals, abduction of the humerus is instrumental to put the manus in the correct position to generate a laterally directed force on the substrate (Gambaryan et al., 2002). The movement of the scapula and the flexion of the shoulder and elbow joint are coordinated with humeral abduction. Among mammals, the range of movement available at the shoulder joint varies greatly (Luo, 2015), however, in fossorial species these movements are commonly limited by factors such as reduction of movement of the shoulder joint. Abduction can take place at the sterno-clavicular joint, and the humeral head can be entirely displaced in the transverse dimension within the scapula glenoid cavity (Gambaryan et al., 2002;Warburton, 2006;Archer et al., 2010;Regnault and Pierce, 2018). The result is that all abduction at the shoulder can be extremely constrained in its range of movement. In other mammals the longitudinal axes of the glenoid fossa is parallel to the plane of the costal surface of the scapula, whereas in fossorial mammals it may be directed at different angles. In this context a measure of this angle provides an appropriate proxy for the mobility of the shoulder joint, where higher angles suggest a reduction in the possible movements at the shoulder. Consequently, we measured the angle between the longitudinal axes of the glenoid fossa and the costal surface of the scapula for all the species included in our sample (see Figures 1A,B) to predict the limitations imposed by the digging kinematics.
We then imported the landmark configurations into R version 3.5.3 (R Core Team, 2019) to further process the landmark coordinates. We rotated, translated, and scaled the landmark configurations to unit centroid size (CS, square root of the sum of squared distances of a set of landmarks from their centroid, Bookstein, 1991) performing a Generalized Procrustes Analysis (GPA, Rohlf and Slice, 1990;Goodall, 1991) on all landmark coordinates. The GPA is implemented in the ProcSym() function from the R package "Morpho" (Schlager, 2014). The multivariate ordination of the aligned Procrustes coordinates has been visualized through a principal component analysis (PCA).

Phylogeny
We built synthetic phylogenies, presented in Supplementary Figure S2, in Mesquite 3.02 (Maddison and Maddison, 2015). The stratigraphic range and phylogenetic position have been derived from the most updated literature (Piras et al., 2015;Sansalone et al., 2019). Supplementary Table S1 includes the complete literature corpus employed to build the phylogenetic trees. We built two different phylogenies because of the lack of agreement between molecular and morphological data for talpids (see Sansalone et al., 2019, for an extensive discussion about the topic). For Mesoscalops, we used the data presented by Gunnel et al. (2008) and Barnosky (1981). The position of chrysochlorids has been based on the data provided by Piras et al. (2015) and O'Leary et al. (2013). The divergence between the echidna, N. philcreaseri and placental mammals has been set following Luo et al. (2007Luo et al. ( , 2011, Rowe et al. (2008); Phillips et al. (2009) and Bi et al. (2018).

RRphylo and Convergence Analysis
We employed a novel Phylogenetic Comparative Method specifically designed to deal with phylogenies including fossil taxa. The RRphylo method (Castiglione et al., 2018) performs a phylogenetic ridge regression on tree and data and returns branch-wise evolutionary rates and ancestral character estimates (ACEs) at nodes. In the case of multivariate data, such a procedure is performed independently for each component of the phenotype, by applying a normalization factor essential to avoid extreme rate values and multicollinearity, under the hypothesis that rate variation is minimized within clades (Castiglione et al., 2018).
To identify instances of phenotypic convergence between either clades or between groups of species within specific categories (i.e., states) we used search.conv, a new method which has proven to be both fast and powerful in identifying morphological convergence . Under search.conv the phenotypic distance between species is quantified as the angle between their phenotypic vectors (i.e., multivariate phenotypes for each species). Under a Brownian Motion (BM) model of evolution, this angle should increase proportionally to the patristic distance between species. When morphological convergence is present, this assumption no longer applies as the angle (per unit time) become smaller than expected. The function search.conv is specifically designed to detect instances of convergence, as occurring either between entire clades ("clade case") or between selected species evolving according to specific states ("state case"). When comparing clades, the function computes the mean angle over all possible combinations of species pairs taking one species per clade, and divides this value by the patristic distance between the nodes subtending to the clades (i.e., the phylogenetic distance between the most recent common ancestors, MRCAs, of the clades). To assess significance, this value is contrasted to a random distribution of 1,000 angleby-time values. The function randomly samples two species within the tree, computes the angle between them and divides it by the patristic distance between their immediate ancestors.
Additionally, search.conv determines whether convergence occurred early in the clade (at divergence) or only applies at the level of terminal branches (i.e., species). Given two clades presumed to evolve under convergence, the function derives the ACEs for the MRCAs for both clades from the RRphylo results and computes the angle between them. Such angle is summed to the mean angle between species and divided by the patristic distance between the MRCAs. If the "summed" angle is smaller than expected by chance, it means the clades converged since their origin and subsequently followed parallel phenotypic evolutionary trajectories. The significance level is assessed as described above by randomizing phenotypes across the tree tips. Finally, if no specific hypothesis about converging clades is available, the function automatically scans the phylogeny to identify instances of convergence.
Under the "state case", search.conv computes the mean angle over all possible combinations of species pairs using one species per state. Each individual angle is divided by the patristic distance between the species. Significance is assessed by contrasting this value with a family of 1,000 random angles obtained by shuffling the state across the species.
Here, we applied search.conv under both "clade" and "state" cases on shape values (PC scores), by classifying species under fossorial/non-fossorial lifestyle. We repeated the analyses by collating the values of functional variables (i.e., angle and von Mises stress values) to PC values for each species and taking such "composite" vectors as the search.conv entry. All the functions used to compute rates, phenotypes and to test convergence are available as part of the R package RRphylo (Raia et al., 2019 available at 2 ).

Shape Analysis
3D shape analysis evidences a clear phylogenetic signal with the different clades clustering in distinct regions of the morphospace (see Figure 2). Specifically, chrysochlorids grouped on the negative end of PC1 and PC2, whereas the echidna and N. philcreaseri separate along positive values of PC2. Highly fossorial talpids (i.e., Talpini and Scalopini) and M. montanensis cluster along positive values of PC1, but fossorial talpids separates between consensus and slightly negative values of PC2, whereas M. montanensis is placed at highly negative values. The remaining taxa (cursorial, semi-fossorial and semi-aquatic moles) are distributed from consensus to negative values of the PC1 and positive values of the PC2 (see Figure 2).
At positive values of PC1 the humerus shows a robust and rounded configuration characterized by a medially displaced and enlarged teres tubercle, the proximal region is relatively developed and the greater tuberosity is highly expanded. At negative values of the PC1 the humeri possess a relatively slender morphology characterized by a less developed teres 2 https://github.com/pasraia/RRphylo tubercle, the proximal region is reduced and the articulation with clavicle is absent, features typical of chrysochlorids and N. philcreaseri. Slender, non-fossorial, talpids cluster at the positive values of PC2, whereas M. montanensis is positioned in a unique region of the morphospace at the negative extreme of the PC2. At positive values of this same axis the humerus presents a relatively slender configuration with an expanded proximal region, whereas the teres tubercle is poorly developed. At negative values of PC2, the teres tubercle is placed on the distal region of the humerus, the distal region broaden and bear an extremely elongated lateral epicondyle, peculiar of the proscalopid M. montanensis.

Finite Elements Analysis
The application of FEA to the humerus models, we revealed that the different functional groups were characterized by different VM stress magnitudes and distributions. The cursorial and semi-aquatic talpids showed the highest stress values, with a highly stressed humeral shaft. The semi-fossorial talpids display intermediate values with most of the stress again concentrated on the humeral shaft (see Figure 3). Fossorial talpids showed overall low stress values, with stress peaks mainly concentrated around the head of the humerus and the clavicular articular facet, whereas the areas of the muscles insertion showed lower VM stress values. The notoryctid showed averaged stress values lower than the cursorial and semi-aquatic talpids but higher than the semi-fossorial talpids and the chrysochlorids. N. philcreaseri displayed a highly stressed humeral shaft, in particular along the proximal half. The deltopectoral crest and the medial epicondyle showed remarkably lower stress values, reflecting the insertion of the main muscles involved during the power stroke (m. triceps longus, mm. pectorales, m. flexor carpi radialis, and m. flexor carpi ulnaris and m. flexor digitorum profundus). The echidna showed stress values lower than N. philcreaseri but comparable with the semifossorial talpids. The highest stress values were located on the proximal half of the humeral shaft, in particular on the region below the humeral head. The area of insertion of the digging muscle and, overall, the greatly enlarged distal area of the humerus showed remarkably lower stress values. In chrysochlorids, high stress values were found on the distal region of the humerus, whereas the lateral and medial epicondyles, having the direct m. latissimus dorsi and m. triceps insertions were particularly stressed. In contrast, the proximal humeral region in chrysochlorids is less stressed compared to talpids because the muscles involved in digging do not insert in that region (Puttick and Jarvis, 1977;Barnosky, 1982;Gasc et al., 1986;Piras et al., 2015). The humerus of M. montanensis showed moderate stress values on the proximal region compared to talpids, though the teres tubercle was more stressed because of its distally shifted position (Piras et al., 2015). PC1 and PC2 scores, VM stress values, and angles for each species are reported in Supplementary Table S3.

RRphylo and Convergence Analysis
By applying search.conv on shape data (PC scores) alone we found no significant instance of convergence, neither under "clade" (angle per unit time = 1.051 degrees/myr, p = 0.117) or under "state" (angle per unit time = 1.240 degrees/myr, p = 0.746) case. When integrating functional variables, we found Japanese moles (i.e., Mogera species exclusive of Mogera insularis) to significantly converge on the clade including the north American moles within Scapanus plus Scalopus genera (i.e., Scalopini, Figures 4A,B), starting from their respective common ancestors (angle per unit time, tips plus ancestors = 0.199 degrees/myr, p = 0.046).
FIGURE 4 | Convergence test on fossorial mammals humeral shape. The convergence test is implemented for entire clades (A,B) and for fossorial species as they appear sparsely on the talpid tree (C,D). PC plots appear in the upper row (A,C). (B) phenogram describing the path of the convergent clades. (D) angle between the fossorial species as compared to the random distribution of angles obtained shuffling the states across species in the tree. The asterisks represent the ancestral states of the clades found to converge. The pink and purple dots represent the mean phenotype of the species within the respective clades found to converge.
As for the "state" case, the mean angle between species within the "fossorial" category becomes 0.404 degrees/myr by adding performance variables, which is significant (p = 0.042; Figures 4C,D).

DISCUSSION
Shape analysis showed the presence of a mixture of taxonomic and functional signals (see Figure 2). All the functional groups occupied different regions of the morphospace, suggesting the absence of a significant convergent signal in the humeral shape of fossorial mammals. This was confirmed by the search.conv() analysis which did not recover evidence of convergence across our morphological sample. This result suggests that fossorial mammals were able to find different solutions to a similar adaptive challenge: the exploitation of the underground. In this case natural selection has produced some of the most divergent and highly derived humeral morphologies among vertebrates (Luo and Wible, 2005;Piras et al., 2015).
Biological structures showing morphological similarity, despite phylogenetic distance, might be used for different purposes. A striking example is represented by the marsupial N. philcreaseri and the placental chrysochlorids, these two taxa share great expansion of the medial epicondyle and the deltopectoral process. In fact, both N. philcreaseri or the chrysochlorids heavily rely on the m. triceps and m. pectoralis action during the power stroke, as evidenced by the enlarged deltopectoral crest and by a similar orientation of these muscles in the these taxa (Jenkins, 1970;Warburton, 2006;Piras et al., 2015;Beck et al., 2016). However, on the medial epicondyle two different muscles groups insert. In the chrysochlorids the m. latissimus dorsi inserts on the medial epicondyle, whereas in the marsupial moles the m. latissimus dorsi inserts on the fascia of the elbow-forearm rather than on the medial epicondyle, where different flexor muscles insert (Warburton, 2006;Beck et al., 2016). This adaptation is unique to the marsupial moles and improves the mechanical advantage of the m. latissimus dorsi for flexion of the shoulder and the effect of the elbow flexion (Warburton, 2006). In this context, it is not surprising that N. philcreaseri, despite the morphological similarities displays different stress distributions to chrysochlorids. Furthermore, our analysis did not identify any morpho-functional convergence between these two clades, suggesting that marsupial-and golden-moles are optimized for different trade-offs.
We found only one instance of morpho-functional convergence accounting for trade-off measures between the Japanese and North-American fossorial moles (Mogera, Scalopus, and Scapanus genera, see Figure 4). This result strongly suggests that fossorial mammals evolved around multiple different optima rather than around a single adaptive peak, as predicted by many-to-one mapping (Losos, 2011). It is likely that natural selection favored the evolution of humeral morphologies characterized by high performance (low stress) rather than promoting the optimization of humeral mobility. However, the trade-off between humeral strength and mobility resulted in different morphologies having similar fitness in the subterranean environment (Marks and Lechowicz, 2006;Losos, 2011).
For example, the lack of significant convergence between European (genus Talpa) and North-American moles may reflect the longstanding separation and phylogeographic histories of these two groups (Ziegler, 1999;Sánchez-Villagra et al., 2006;Schwermann et al., 2019). Scalopine moles went extinct in Europe during the Messinian salinity crisis (Krijgsman et al., 1999), whereas European mole (Talpini) diversity rose after the crisis ended, possibly exploiting the ecological niche left empty by the extinction of the Scalopini in Eurasia (Ziegler, 1999;van den Hoek Ostende and Fejfar, 2006). In this scenario, differences in the initial conditions for the two populations (Scalopini moles survived in North-America only, with one relict species in China; Scapanulus oweni), the phenotypes, genetic variation and the constraints acting on the ancestral populations, may have driven the evolution of trade-off values into different pathways for the two fossorial mole sub-clades .
Yet, the lack of convergent signals between clades evolving under similar selective conditions could be the result of phylogenetic history. In such a scenario, two species might be more similar to each than to their ancestors, but the magnitude of trait changes was not great enough to rule out pre-existing interclade differences (Losos, 2011). In these cases species do not reach "complete convergence" (Herrel et al., 2004) because natural selection might optimize their phenotypes within different selective contexts (environmental conditions might not be identical), or that developmental constraints might preclude two species to evolve around the same adaptive peak.
In the present study this is exemplified by the similar medial and proximal expansion of the humerus in both T. aculeatus and the fossorial talpids (Jenkins, 1970). However, in the former a wide humerus is an adaptation to rotation with the limitation of abduction, providing horizontal placement of the manus to sprawl the soil. On the other hand the expansion of the humerus in Talpidae is an adaptation to its abduction with limitation of rotation providing parasagittal positioning of the hands during burrowing. Hence, the wide humerus of the echidna and the similar humeral expansion in moles evolved pointing to opposite directions and opposite functions (Gambaryan, 2000).
Many examples of convergent evolution are based on simplistic descriptions of both function and phenotypes, but accounting for trade-off measures between different behaviors under similar selective contexts may reveal more complex scenarios. Here we have shown how, contrarily to expectations, fossorial mammals display only a very limited signal of morphofunctional convergence. We have detailed means by which different factors can play key roles in excluding species from evolving toward a single adaptive optimum, as predicted by the many-to-one mapping. On the contrary, our results show that fossorial mammal lineages have found different solutions to a similar adaptive regime evolving around multiple selective peaks defined by different biomechanical trade-offs. Our investigations highlight the fact that natural selection, even operating in similar selective conditions, may fail to produce convergent phenotypes in different evolutionary lineages. In this context, future investigations of morpho-functional convergence should expand on other skeletal elements participating in performing a particular function, this would allow a greater accuracy in defining the extent of convergence.

DATA AVAILABILITY STATEMENT
The data used in all the analyses are available in Supplementary Material. The raw data used to generate the results form this study have been made available at: doi: 10.25952/5e65995b36b6c.

ETHICS STATEMENT
The sample used in this study is entirely composed by specimens conserved in museums, including fossils for which no ethical approval is required.

AUTHOR CONTRIBUTIONS
GS and PR conceived the study. GS, PR, SC, AP, and PP analyzed the dataset. GS, SW, SH, AP, PP, and MA sampled and processed the dataset. GS, PR, SC, and BD wrote the manuscript with contribution form all the co-authors. Botanico Tropical/IICT, Lisboa, Ursula Goelich from Wien Natural History Museum, Ugo Messana and Enrico Francucci from Studio Dentistico Messana. We are grateful to Ronald Eng and Timothy C. Cox from the Bourke Museum, Seattle, US. We are grateful to the two referees for their useful comments.