Evolution of methicillin-resistant Staphylococcus aureus towards increasing resistance

Abstract

Objectives

To elucidate the evolutionary history of Staphylococcus aureus clonal complex (CC) 8, which encompasses several globally distributed epidemic lineages, including hospital-associated methicillin-resistant S. aureus (MRSA) and the highly prevalent community-associated MRSA clone USA300.

Methods

We reconstructed the phylogeny of S. aureus CC8 by mutation discovery at 112 genetic housekeeping loci from each of 174 isolates, sampled on five continents between 1957 and 2008. The distribution of antimicrobial resistance traits and of diverse mobile genetic elements was investigated in relation to the isolates' phylogeny.

Results

Our analyses revealed the existence of nine phylogenetic clades within CC8. We identified at least eight independent events of methicillin resistance acquisition in CC8 and dated the origin of a methicillin-resistant progenitor of the notorious USA300 clone to the mid-1970s. Of the S. aureus isolates in our collection, 88% carried plasmidic rep gene sequences, with up to five different rep genes in individual isolates and a total of eight rep families. Mapping the plasmid content onto the isolates' phylogeny illustrated the stable carriage over decades of some plasmids and the more volatile nature of others. Strikingly, we observed trends of increasing antibiotic resistance during the evolution of several lineages, including USA300.

Conclusions

We propose a model for the evolution of S. aureus CC8, involving a split into at least nine phylogenetic lineages and a subsequent series of acquisitions and losses of mobile genetic elements that carry diverse virulence and antimicrobial resistance traits. The evolution of MRSA USA300 towards resistance to additional antibiotic classes is of major concern.

MRSA, phylogeny, acquisition of resistance, mobile genetic elements, plasmids

Introduction

The worldwide emergence and spread of methicillin-resistant Staphylococcus aureus (MRSA) over the last 50 years represents one of the most serious challenges to clinical microbiologists worldwide. The first MRSA isolates, which appeared in the UK and Denmark shortly after the introduction of penicillinase-stable methicillin in 1959, were assigned to the so-called Archaic MRSA clone.¹ Methicillin resistance had been acquired through horizontal transfer of the mecA gene from an unknown source into a methicillin-susceptible S. aureus (MSSA). Comparison of early MRSA and contemporary MSSA and MRSA isolates revealed close phenotypic and genotypic similarity, suggesting specific MSSA clones as ancestors of the first MRSA in Europe.¹ These ancestors belonged to clonal complex (CC) 8, which encompasses several MRSA, including the Archaic clone (ST250-MRSA-I, where ST stands for sequence type),^2–4 the Iberian clone (ST247-MRSA-I)³ and other multiresistant, hospital-associated MRSA of more geographically limited importance, such as the Hanover clone (ST254-MRSA-IV; UK EMRSA-10), which was prevalent in Europe during the 1990s. The Brazilian/Hungarian clone (ST239-MRSA-III) and the Portuguese clone (ST239-MRSA III variant) are also related to CC8, as they originated through a recombination event involving CC8 and CC30 S. aureus.⁵ Further, one of the most prevalent community-associated MRSA (CA-MRSA) clones, termed USA300 (ST8-MRSA-IV), belongs to CC8 and its emergence likely represents a recent evolutionary process within this CC.⁶

We have previously demonstrated the utility of a genome-wide mutation discovery strategy to elucidate the population structure of S. aureus.^7–10 In this study we used this approach to reconstruct the evolutionary history of CC8. We included contemporary MRSA and MSSA as well as isolates from the 1950s through the 1960s,^1,3 collected from five continents.

Materials and methods

Bacterial isolates and their characterization

A total of 174 CC8 isolates from 33 different countries on five continents, sampled between 1957 and 2008, were included in this study (Table S1, available as Supplementary data at JAC Online). Isolates to be included were chosen to maximize geographical and temporal coverage. All isolates were subjected to antimicrobial susceptibility testing by the use of microbroth dilution according to EUCAST guidelines (www.eucast.org). The following clinically or epidemiologically relevant antibiotics were tested: penicillin, oxacillin (antibiotic class: β-lactams), gentamicin (aminoglycosides), linezolid (oxazolidinones), erythromycin, clindamycin (macrolides–lincosamides–streptogramins B, MLS_B), tetracycline (tetracyclines), vancomycin, teicoplanin (glycopeptides), ciprofloxacin, moxifloxacin (quinolones), daptomycin (lipopeptides), mupirocin, fosfomycin, rifampicin, fusidic acid and trimethoprim/sulfamethoxazole. Genomic DNA was isolated from overnight cultures with the DNeasy Tissue Kit (Qiagen, Hilden, Germany) using lysostaphin (100 mg/L, Sigma, Taufkirchen, Germany) to achieve bacterial lysis. spa typing and multilocus sequence typing (MLST) were performed as described previously.^11–13 SCCmec elements were typed by a PCR approach including a combination of different primer sets as described before¹⁴ and supplemented for the detection of ccrC (Table S2, available as Supplementary data at JAC Online). LukSF-PV, arcA (ACME) and msrA were detected by PCR using the primers and PCR conditions summarized in Table S2, available as Supplementary data at JAC Online. All PCR reactions were carried out using Fermentas PCR Mastermix (ThermoScientific, Dreieich, Germany). Further characterization of lukSF-PV type was done as reported by O'Hara et al.¹⁵ Plasmid typing was performed as described recently.¹⁶

Mutation discovery

A total of 112 housekeeping loci were investigated for the occurrence of sequence polymorphisms in comparison with S. aureus COL (GenBank accession number NC_002951). All loci were ∼400–500 bp in length (Table S3, available as Supplementary data at JAC Online), thus covering a total of ∼1.5% of the genome from each of the 174 isolates. They were amplified by PCR and amplicons were further analysed by denaturing HPLC (dHPLC-WAVE; Wave^® Nucleic Acid Fragment Analysis System, Transgenomic, Glasgow, UK) as described previously.⁸ All mutations identified were verified by capillary Sanger sequencing. All loci investigated, PCR primers and annealling temperatures are provided in Table S3, available as Supplementary data at JAC Online.

Data analysis

Based on all single-nucleotide polymorphisms detected, we constructed a minimum spanning tree using BioNumerics 6.5 (Applied Maths, Sint-Martens-Latem, Belgium). Additionally, we concatenated the sequences of all loci investigated for each of the 174 isolates, resulting in a sequence alignment 41 679 bp in length. Based on this alignment we constructed a maximum likelihood tree using the PhyML v3.0.1 algorithm included in SeaView 4.2.11 and applying the HKY model of DNA substitution.^17,18 This tree was rooted by including the corresponding sequence from the distantly related S. aureus N315 (NCBI GenBank accession number NC_002745). Visualization and annotation of the tree with external datasets was performed using iTOL (http://itol.embl.de).¹⁹ Root-to-tip distances in the maximum likelihood tree were determined with Path-o-gen software (http://tree.bio.ed.ac.uk/software/pathogen/). DnaSP v5 was used to estimate the nucleotide diversity (π).^20,21 PAUP 4.0 (http://paup.csit.fsu.edu/) was used to calculate the homoplasy index and likelihood scores for maximum-likelihood trees, with and without a molecular clock enforced, respectively. For a likelihood ratio test, the statistical significance of differences between likelihood scores was assessed by assuming a χ² distribution and s − 2 degrees of freedom, where s was the number of sequences included.²²

Evolutionary rates and divergence times were calculated with the BEAST software, 1.7.5 (http://beast.bio.ed.ac.uk/),²³ using the HKY model of nucleotide substitution and strict as well as relaxed clock models based on the concatenated protein coding sequences dated with the year of isolate sampling and with 10⁸ iterations after a burn-in phase of 10⁷ iterations. Data from three independent analyses were averaged to estimate posteriors, resulting in effective sample size values close to or greater than 1000 for all parameters. The resulting set of sampled trees was used to estimate the strength of phylogeny–trait associations for geographical origin, spa types, SCCmec types, virulence gene and plasmid type content, applying Bayesian tip-association significance testing to calculate the maximum monophyletic clade (MC) statistic (BaTS, version 1.0²⁴). For these analyses, indistinguishable isolates were removed from the dataset and 100 trees were sampled from BEAST results.¹⁰ The null hypothesis of no association between trait and phylogeny was estimated on the basis of 100 randomizations of isolate–character associations, and the null hypothesis was rejected if P < 0.05.¹⁰

Results and discussion

Population structure, phylogeography and evolutionary rates

We performed mutation discovery on 1.5% (41 679 bp) of the genomes from each of 174 S. aureus isolates by using denaturing HPLC.⁸ Ascertained polymorphisms and their properties are summarized in Table S4 (available as Supplementary data at JAC Online); they comprised 199 bi-allelic single-nucleotide polymorphisms (SNPs) resulting in 103 different haplotypes. Sixty-two SNPs were parsimony informative, as their derived alleles were found in more than one haplotype. The nucleotide diversity π (the average number of nucleotide differences per site between sequences from two isolates) was 0.00019 ± 0.00001, which is approximately twice the nucleotide diversity found in S. aureus lineage ST5.⁸ A very low level of homoplasy (homoplasy index, 0.03) indicated that CC8 is a highly clonal population, similar to other lineages of S. aureus.^8,10,25

Our maximum likelihood phylogenetic analyses revealed nine clades within CC8 (Figure 1), and a minimum spanning tree displayed an identical clade distribution (Figure S1, available as Supplementary data at JAC Online). Each of these clades can be identified on the basis of specific SNPs (Table S5, available as Supplementary data at JAC Online). Of these nine clades (CC8-A to CC8-I), two were identified by conventional MLST: clade CC8-D, which corresponds to ST254 (Figure S1, available as Supplementary data at JAC Online), an ST that was prevalent among multiresistant, hospital-associated MRSA in Europe during the 1990s; and clade CC8-E, which includes all ST250 and ST247 isolates (Figure S1, available as Supplementary data at JAC Online), representing the Archaic and Iberian MRSA clones. Since these clones represent early epidemics of infections with MRSA,²⁶ it is not surprising that clades CC8-D and CC8-E contain the majority of historic isolates in our collection, which had been collected between the 1950s and 1980s (Figure S2, available as Supplementary data at JAC Online). The remaining clades were dominated by ST8 (Figure S1, available as Supplementary data at JAC Online). Clade CC8-A contains MSSA and MRSA isolates from the 1990s and 2000s from four different continents; this clade includes USA300, an MRSA frequently found in community-associated infections in the USA and beyond, and the direct ancestors of this clone. Isolates in clade CC8-B were collected in Europe during the 1990s and 2000s (Figure S2, available as Supplementary data at JAC Online). This clade includes a number of Danish isolates with spa type t024, which is the predominant healthcare-associated MRSA (HA-MRSA) in Copenhagen. This clone, which was first identified in 2003, accounts for 28% of all MRSA cases from 2003 to 2008. Isolates in clade CC8-C were collected from five continents; most of these isolates were obtained during a large international multicentre study of S. aureus in hospitals.²⁷ Clades CC8-F to CC8-I represent minor clades with limited numbers of isolates.

Figure 1.

Phylogenetic relationships of 174 CC8 isolates from five different continents. Maximum-likelihood phylogenetic tree based on 199 SNPs ascertained in 112 genetic loci (41 679 bp total). The tree is annotated with (I) SCCmec types, (II) possession of genetic determinants lukSF-PV, arcA and msrA, and (III) detection of plasmid replication (rep) genes rep7, rep10, rep13, rep15, rep16, rep20, rep21 and rep22.

Compared with previous findings for various clonal lineages,^8,25,28 we found relatively little association between phylogeny and geographical origin. Accordingly, Bayesian tip-association significance testing (using BaTS analysis²⁴) indicated that isolates from most countries were not more strongly associated with any specific clades in the phylogenetic tree than would be expected from a random distribution (i.e. P > 0.05; Table S6, available as Supplementary data at JAC Online). In clades CC8-A, CC8-C, CC8-D and CC8-E, isolates from different continents cluster together in identical haplotypes. Clade CC8-B consists of isolates from European countries exclusively, however, and the existence of additional continent-specific clades might be masked by the predominance of European isolates in our collection. Thus the absence of a phylogeographical association may in part be due to the small sample size in relation to the large geographical area and time span covered.

Based on the sequence alignment of protein coding sequences from isolates that had been sampled over 51 years (Table S1, available as Supplementary data at JAC Online), we calculated an average nucleotide substitution rate for CC8 of ∼1.8 × 10⁻⁶ substitutions per nucleotide site per year (relaxed clock model, 95% CI, 1.2 × 10⁻⁶ to 2.4 × 10⁻⁶). This average rate is similar to rates previously reported for other clonal lineages of S. aureus.^7,25,28,29 However, substitution rates varied significantly among clades, as indicated by the result of a likelihood ratio test (P < 1 × 10⁻⁹),²² justifying the use of a relaxed molecular clock model.³⁰ When using a relaxed clock (uncorrelated lognormal model), the mean coefficient of variation of the substitution rate was 1.3 (95% CI, 0.8–1.9), supporting the notion of rate heterogeneity among clades.³⁰ This is remarkable as it had not been observed in previously analysed datasets on other CCs of S. aureus,^7,10,25 suggesting that CC8 shows unusually strong variation of evolutionary rates. These preliminary results warrant a systematic investigation of substitution rate variation in CC8, which will, however, require a larger dataset including additional isolates and/or longer sequences (e.g. complete genome sequences).

Evolution of antibiotic resistance

When we tested in vitro susceptibilities to 18 antimicrobial substances encompassing 13 antibiotic classes, CC8 isolates displayed diverse resistance patterns (Table S1, available as Supplementary data at JAC Online), ranging from full susceptibility in some isolates to broad resistance to up to nine different antibiotic classes in others (Table S1 and Figure S3, available as Supplementary data at JAC Online). The majority of broadly resistant isolates were found in clade CC8-E, encompassing the multiresistant hospital-associated Iberian clone (ST247; Table S1 and Figure S3, available as Supplementary data at JAC Online).³ However, isolates resistant to five to seven antibiotic classes were found in all CC8 lineages (Figure S3, available as Supplementary data at JAC Online). Remarkably, for clades CC8-A and CC8-E we found significant positive correlations (P < 0.01 and P < 0.001, respectively) between root-to-tip distances in a maximum-likelihood phylogenetic tree and the number of antibiotic classes to which the corresponding isolates were resistant (Figure 2 and Figure S3, available as Supplementary data at JAC Online). This result indicates a continuous gain of antibiotic resistances during the evolution of these lineages, at least for the limited number of antibiotics tested in this study. The increase in resistance in CC8-E displays a steeper slope than that in CC8-A, which may reflect differential selection pressures in mostly hospital-associated MRSA (CC8-E) and community-associated MRSA (i.e. USA300 in CC8-A). Isolates in clade CC8-E have acquired a higher number of resistance traits (median, 6; maximum, 9) than those in CC8-A (median, 2; maximum, 5), which may be a result of the older age of CC8-E, reflected by longer branch lengths in the phylogenetic tree (Figure 2). In any case, these results indicate that CC8-A/USA300 is on an evolutionary trajectory towards increasing resistance, which is corroborated by recent reports on emerging resistance in USA300 MRSA,^31,32 and which may predict that this strain is likely to gain resistance to additional antimicrobials in the future. Similarly, an increase in antibiotic resistance over time has previously been observed in hospital-associated MRSA in CC22.²⁵ In contrast, no significant correlation was found for clades CC8-B (P = 0.7; Figure 2), CC8-C (P = 0.7) and CC8-D (P = 0.5). For clades CC8-C and CC8-D this may be due to insufficient sampling size and temporal coverage. Clade B contains the Danish t024 clone (haplotype CC8-B-3; Figure S1, available as Supplementary data at JAC Online) and its progeny at the top of the clade and almost all of the isolates clustering here were collected in Northern Europe. In contrast to USA300, the Danish t024 clone has had a very slow expansion, with only 308 cases identified between 2003 and 2008. Due to the Danish MRSA ‘search and destroy’ policy, this is a close estimate of the true number of cases in this time period. Overall, our sample of isolates is limited. Moreover, our collection of isolates is heavily skewed towards those identified in healthcare settings; most probably we would see a different development in lineages containing isolates that are not put under antibiotic-driven selection pressure.

Figure 2.

Number of antibiotic classes affected by resistance plotted against the root-to-tip distance for individual isolates of S. aureus CC8 (CC8-A, red; CC8-B, blue; CC8-E, orange). The graph includes regression lines with 95% CI for each group. Sizes of circles correspond to the number of isolates.

By mapping structurally different SCCmec elements onto the SNP-based phylogeny we identified eight independent SCCmec uptake events in CC8 (Figure 1). The predominance of SCCmecI in clade CC8-E reflects the evolution of ST250-MRSA-I and its most successful descendant, ST247-MRSA-I, from the early 1960s until now (Figure 1 and Figure S1, available as Supplementary data at JAC Online). Haplotypes close to the root of clade CC8-E contain MSSA and MRSA from several decades, confirming previous suggestions concerning the evolution of the first MRSA in Europe.^2–4 Among MRSA in all other clades within CC8, SCCmecIV is the predominant type of SCCmec. In addition, a few isolates with diverse phylogenetic affiliations carry SCCmecV (n = 5; geographical origin: Denmark, Nigeria, Singapore) or SCCmecVI (n = 1; Japan).

Plasmid content

We used a recently described approach for classification of plasmids in S. aureus,¹⁶ which assigns plasmids to rep families based on PCR detection of specific rep gene sequences. In our collection of CC8 isolates, we detected eight of the previously described rep families, including rep7 (n = 29), rep10 (n = 21), rep13 (n = 8), rep15 (n = 6), rep16 (n = 23), rep20 (n = 88), rep21 (n = 71) and rep22 (n = 42) (Figure 1). We found rep genes in 153 isolates (88%) from our collection. Individual isolates carried up to five different rep genes (Figure 1 and Table S1, available as Supplementary data at JAC Online), reflecting the carriage of multiple different plasmids per isolate and/or the presence of several rep genes on some of the plasmids.^16,33 While some rep families were associated with specific clades in the phylogenetic tree (as statistically supported by BaTS analysis at P < 0.01; Table S6, available as Supplementary data at JAC Online; BaTS results), plasmids with rep genes from other families were more widely distributed across sub-lineages within CC8 (Figure 1 and Table S6, available as Supplementary data at JAC Online; BaTS). The association of rep families rep16 and rep21 with clade CC8-A is in agreement with the previously described, uniform presence of a rep16 plasmid of 27 kbp, carrying multiple resistance genes including msrA, and a small cryptic rep21 plasmid of 2.5 kbp in epidemic USA300 isolates from the USA.^33,34 Our data confirm the presence of both these plasmids together with the msrA gene in 19 USA300 isolates from four continents, with a few exceptions where either one of the two plasmids had been lost (Figure 1). Several USA300 isolates in our collection carried additional rep genes, including one isolate with rep15, which corresponds to a conjugative plasmid encoding the tra complex and the resistance genes erm(C) and ileS.³⁴ As a result, this isolate (09-00124) is resistant to five classes of antibiotics (Table S1 and Figure S3, available as Supplementary data at JAC Online), indicating the important role of plasmid transfer in the evolution of antibiotic resistance. Of note, MSSA and MRSA isolates in our collection that represent the predecessors of the USA300 strain, branching off deeply within clade CC8-A, carry a different plasmid with a rep20 gene, which later was lost and replaced with rep16 and rep21 plasmids.

spa typing

spa typing revealed 26 different spa types among the 174 isolates, with four spa types (t008, t051, t024 and t064) covering ∼80% of all isolates (Table S1 and Figure S4, available as Supplementary data at JAC Online). spa type t008 is presumed to be the ancestral spa type within CC8;¹² it is non-randomly distributed in the SNP-based phylogenetic tree, as it is rare in both lineage CC8-C and lineage CC8-E (P < 0.05; Table S6 and Figure S4, available as Supplementary data at JAC Online). spa type t064 was associated with clade CC8-C (P < 0.05; Table S6 and Figure S4, available as Supplementary data at JAC Online). In contrast, the distribution of spa types t024 and t211 on the phylogenetic tree suggests convergent evolution of spa, which has been previously described for other clonal lineages^8,28,35 and which restricts the utility of spa typing for spatial epidemiology.⁹

Distribution of lukSF-PV and arcA

Twenty-four isolates were positive for lukSF-PV (Table S1, available as Supplementary data at JAC Online and Figure 1) and they all share the same PVL-ST ‘R’.¹⁵ All these isolates are affiliated to clade CC8-A, suggesting lukSF-PV was acquired only once in CC8. This is in contrast to other CCs, where multiple acquisitions of lukSF-PV-carrying phages were evident.^8,10,25,36 Clade CC8-A contains only three isolates devoid of lukSF-PV, which has been discussed extensively for its pathogenic relevance in community-acquired S. aureus infection and which is generally accepted as an epidemiological marker associated with community-acquired MRSA.^37,38 Interestingly, the lukSF-PV-negative isolates are located closest to the root of clade CC8-A (Figure 1) and represent isolates with the oldest sampling date in this clade (1996; Figure S2, available as Supplementary data at JAC Online), indicating that the acquisition of lukSF-PV occurred early in the evolution of CC8-A.

arcA, the second epidemiological marker previously connected primarily to CA-MRSA USA300, was detected in clades CC8-A and CC8-B, indicating at least two independent acquisition events for CC8 (Table S1, available as Supplementary data at JAC Online and Figure 1). The detection of lukSF-PV-positive, arcA-negative isolates and their position in clade CC8-A indicates the sequential uptake of the corresponding mobile genetic elements during the course of evolution (Figure 1), eventually leading to the emergence of the USA300 clone.⁶ In contrast, lukSF-PV-negative, arcA-positive isolates that occur in clade CC8-B clearly are unrelated (Figure 1). Such isolates have previously been described for CC8 and beyond;^38–40 they may carry arcA on ACME elements that are structurally distinct from those elements in USA300.^41,42 Six out of eight arcA-positive isolates in clade CC8-B share the same spa type t024 and seven originate from Denmark, indicating the existence of a locally expanding t024 MRSA that is unrelated to t024 isolates from other countries. These isolates are characterized by variations in the J3 region of SCCmecIVa, often including the arcA genes.^40,41,43 Whether this clone originally evolved in Denmark has previously been unknown. Our analyses show that three of the arcA-positive t024 isolates from Denmark are affiliated to SNP-based haplotype CC8-B-3 (Figure S1 and Table S1, available as Supplementary data at JAC Online), which also includes an MSSA isolate from Copenhagen/Denmark from 1967 (09-00010), indicating this haplotype's presence in the region for several decades and suggesting the local evolution of t024 MRSA. The persistence over several years of locally endemic MRSA clones has recently been reported also for other regions.⁴⁴

Stepwise evolution of USA300

Based on the distribution of marker genes on the SNP-based phylogeny (Figure 1), we suggest a model for the stepwise evolution of epidemic USA300 from a methicillin-susceptible ancestor, involving a series of acquisitions and losses of mobile genetic elements as follows: (i) acquisition of rep20 plasmid by a methicillin-susceptible progenitor; (ii) subsequent acquisition of SCCmecIV element and lukSF-PV prophage; Bayesian coalescence analysis dated the emergence of this methicillin-resistant progenitor of USA300 to approximately 1976 (95% CI, 1960–1989); (iii) subsequent replacement of the rep20 plasmid by rep21 and rep16 plasmids (encoding msrA), and acquisition of an ACME element (encoding arcA) occurred; later, occasional acquisition of additional plasmids and loss of individual plasmids. Additional, more subtle changes, including point mutations and small deletions, have evidently occurred in the remainder of the genome.⁶ During its evolutionary history, USA300 has accumulated a growing number of antibiotic resistance traits; if this trend cannot be stemmed, USA300 is likely to become a multiresistant pathogen in the future.

Funding

This study was partially supported by the Robert Koch Institute.

Transparency declarations

None to declare.

Acknowledgements

We thank Mike Henkel, Annette Weller, Franziska Erdmann, Christa Cuny, Birgit Pasemann and the staff at the RKI central sequencing laboratory for excellent technical assistance. We thank Angela Kearns, HPA London, for contributing data concerning the isolates from the HARMONY collection. We are grateful to multiple colleagues who provided their S. aureus strains. We especially acknowledge the contribution of well-characterized strains and helpful comments on the manuscript from Robert Skov and Anders Rhod Larsen (Statens Serums Institute, Copenhagen, Denmark), Teruyo Ito (Juntendo University, Tokyo, Japan), Fred Tenover (Cepheid, Sunnyvale, CA, USA), Karl Gustaf Kristinsson and Barbara Holzknecht (Reykjavik, Iceland/Hvidovre, Denmark) and Tse Hsien Koh (Singapore General Hospital, Singapore, Singapore). Furthermore, we acknowledge the help of Catarina Milheiriço in the selection of the isolates originally characterized at The Rockefeller University, New York, NY, USA and at Instituto de Tecnologia Química e Biológica (ITQB), Universidade Nova de Lisboa, Oeiras, Portugal. Two anonymous reviewers helped to improve the manuscript by providing valuable comments on an earlier version.