Antibiotic Resistance Gene Diversity and Virulence Gene Diversity Are Correlated in Human Gut and Environmental Microbiomes

Every year, thousands of tons of antibiotics are used, not only in human and animal health but also as growth promoters in livestock. Consequently, during the last 75 years, antibiotic-resistant bacterial strains have been selected in human and environmental microbial communities. This implies that, even when pathogenic bacteria are the targets of antibiotics, hundreds of nonpathogenic bacterial species are also affected. Here, we performed a comparative study of environmental and human gut microbial communities issuing from different individuals and from distinct human populations across the world. We found that antibiotic resistance and pathogenicity are correlated and speculate that, by selecting for resistant bacteria, we may be selecting for more virulent strains as a side effect of antimicrobial therapy.

including tissues such as skin, mucosa, and, most importantly, the gastrointestinal tract. Thus, even when virulent bacteria are the targets of antibiotics, the administration of these drugs may also affect many nonpathogenic mutualistic or commensal bacterial species present in individuals undergoing treatment (3).
The resistome (the collection of all antibiotic resistance [AR] genes, which exist in both pathogenic and nonpathogenic bacteria [4]) is frequently carried on mobile genetic elements. Similarly, the virulome, the set of genes encoding virulence, can also be carried on the mobilome (5)(6)(7). Therefore, many bacterial virulence factors (VFs) are easily spread in bacterial populations by horizontal gene transfer, converting mutualistic or commensal bacteria into potential opportunistic pathogens. Several examples of virulent (and highly virulent) and multiresistant disseminated clones have already been reported throughout the literature (8-10; for a review, see reference 11).
A given bacterial population may constitute a kind of life insurance to other bacterial populations present in the same microbiome through at least two different mechanisms. First, if antibiotic resistance genes are carried by plasmids or other mobile genetic elements, they may transfer to pathogenic cells and save them from the negative effects of antibiotics. Indeed, these genetic elements can spread into the bacterial community by horizontal gene transfer, even crossing species (12)(13)(14). Curiously, some bacterial strains from different species have been shown to be extremely good donors of certain plasmids. As such, they are able to amplify the number of plasmids in a bacterial community and spread those plasmids to other bacterial cells (15), a phenomenon probably explained by interactions between different plasmids (16,17).
Second, even bacterial cells not coding for antibiotic resistance determinants nor receiving mobile elements may be protected by cells coding for certain drug resistance genes. Indeed, gene products that inactivate antibiotics by degrading or modifying antibiotic molecules are also decreasing their concentration in the local environment. (For a review of mechanisms of possible indirect resistance, see reference 18.) This mechanism of indirect resistance has been shown to occur in different systems and to be pervasive (19)(20)(21)(22)(23).
In this article, we ask whether the abundance of antibiotic resistance protein families correlates with that of virulence among bacterial populations. For this, we have chosen to study metagenomes. There are three main reasons to use metagenomes to address this question. First, it is known that, for millions of years, bacteria had to cope with the presence of many other species, mostly competing with them (24), but also cooperating, both phenomena relying on the full set of genes of the metagenome (24). Second, horizontal gene transfer promotes genetic relationships between species, thus enforcing cooperation (6), avoiding the emergence of cheaters within the microbiome. Third, the study of metagenomes gives us access to the repertoire of genes involved in adaptation to the environment, given that many of these traits are often encoded in the mobilome and thus can be shared by different, eventually unrelated bacteria. In this context, mining for genes coding for antibiotic resistance and virulence traits in metagenomes is a reliable way to access the selective pressures the population is subject to, as well as the cooccurrence of genetic traits of the whole microbiome.
The present study aims at understanding the relationship between antibiotic resistance genes and those coding for virulence. We show that there is indeed a linkage between the dissemination of virulence factors and genes coding for antibiotic resistance, within natural microbiomes, and that this relationship could be influenced by the behaviors of human populations spanning very different geographical locations across the world.

RESULTS
AR gene families in the metagenomes. We started by evaluating the quantity of different antibiotic resistance (AR) proteins present in the metagenomes under study. We wanted to know how does the diversity of antibiotic resistance gene families' homologues (i.e., AR diversity [ARd], the number of different protein families) vary with the metagenome protein family richness (i.e., the total number of protein families in a metagenome)? To answer this question, we used a data set composed of natural metagenomes issuing from diverse ecosystems and biomes, such as oceans, coral atolls, deep oceans, Antarctic aquatic environments and snow, soils, hypersaline sediments, sludges, microbial fuel cell biofilms, and animal microbial populations (25), and which from now on will be referred to here as environmental metagenomes, in contrast to those composed solely of human gut metagenomes. In this data set, there is a broad variation in AR gene diversity (ARd) (Fig. 1A). However, there is a correlation between ARd and the protein family richness of the various metagenomes (overall correlation of 0.754). In Fig. 1A, regression lines for human gut microbiome (in blue, slope of 0.0037) and for other biomes (in black, slope of 0.0052) are shown. Nevertheless, the P value for the difference of these slopes is 0.014, which is significant but not highly significant, hence not allowing us, in our opinion, to draw solid conclusions from this difference.
Another data set was composed of 110 human gut metagenomes sampled from healthy individuals with ages ranging from 0.05 to 53 years, spanning different regions of the world: the United States (USA), Malawi, and Venezuelan Amerindians from Amazonia (26) ( Fig. 2A). Metadata can be retrieved at https://bioinfo.mnhn.fr/abi/ public/escudeiro-nogueira/Escudeiro_ARd_VFd_human_gut.xlsx.
VF gene families in the metagenomes. Bacterial virulence factors are often proteins that enable pathogenic bacteria to parasitize a host, including gene products involved in adhesion and invasion, secretion, toxins, and iron acquisition systems (27).
The environmental microbiomes of our data set reveal a great diversity of virulence factor (VFd) protein families' densities ( Fig. 1B). Since virulence may also be associated with the colonization of different types of biomes (some virulence traits are involved in adaptation to new environments, by allowing bacteria to adhere to and colonize substrates, to access resources such as iron, among other functions) besides the context of infection, one can expect different types of these genes in environmental microbiomes. Moreover, Fig. 1B and 2B show that in human gut microbiomes, VFd versus protein family richness of the metagenomes achieves an overall correlation coefficient of 0.669. Regression lines for human gut microbiome (in blue, slope of 0.0035) and for other biomes (black, slope of 0.0106) are shown. Here the P value for the difference of these slopes is 1.7 ϫ 10 Ϫ8 , indicating that there could be a difference, but the points are widely spread over all types of biome, and we prefer not to draw conclusions about these data (metadata can be retrieved at https://bioinfo.mnhn.fr/abi/public/escudeiro -nogueira/Escudeiro_ARd_VFd_human_gut.xlsx).
AR/VF correlations. The main purpose of the present work is to evaluate the relationship, if any, between ARd and VFd in microbiomes. Therefore, we excluded from our analysis all the gene products that are homologous both to AR and VF determinants. Thus, we avoid the introduction of a bias due to both correlations of AR and VF diversities (ARd and VFd) with metagenome richness. In the analysis, we systematically analyzed the residuals of the regression of ARd on metagenome richness and the residuals of the regression of VFd on metagenome richness. In Fig. 3 one can see that the residuals of ARd and those of VFd are positively correlated. Although there is a large variation of ARd (Fig. 1A), and VFd ( Fig. 1B) in environmental metagenomes, there is a strong correlation between ARd and VFd ( Fig. 3): the overall partial correlation between ARd and VFd given metagenome richness is 0.84 (P ϭ 8 ϫ 10 Ϫ36 [i.e., ϳ0.0]). Indeed, the partial correlation between ARd and VFd given the metagenome richness is the same as the correlation between residuals of ARd and residuals of VFd (the partial correlation coefficient for environmental biomes is 0.86, while the one for the human gut is 0.90). When we draw the regression line for the human gut microbiomes from this data set (human gut), we can see that the ratio of residuals of ARd to residuals of VFd is higher, with a slope of 0.63, than for the case of environmental metagenomes, with a slope of 0.44. The difference in slopes is significant, but not very significant, with a P value of 0.03. Despite the fact that the values for ARd and VFd are lower in gut metagenomes ( Fig. 1), one can witness a higher ARd/VFd ratio (larger slope) of  metagenomes pertaining to the human feces biome relative to the slope calculated for all environmental metagenomes (Fig. 3). This suggests that there could be a greater accumulation of ARd over VFd for this particular biome than in environmental ones. Figure 4 shows residuals of AR diversity versus the residuals of VF diversity in the human gut metagenomes for each population (country) and for the three populations together. The partial Pearson correlation coefficient between ARd and VFd given metagenome richness is lower in the human gut metagenome data set (r ϭ 0.68; P ϭ 3 ϫ 10 Ϫ22 [ϳ0.0] [ Fig. 4]) than in the environmental samples (r ϭ 0.86 as seen above [ Fig. 3]).
We can further distinguish different trends upon geographical localization of the human populations under study. The largest contribution to this graph comes from the North American samples, which account for 66/110 (60%) of the individuals (Fig. 4). The Amerindians (Venezuela) account for 21/110 metagenomes (19%), and the Malawians account for 23/110 (21%). These two non-Western human populations have very contrasting lifestyles and exposure to antibiotics: the uncontacted Amerindians harbor a resistome similar to that of pre-antibiotic era, and the people from Malawi have a high prevalence of infectious diseases and widespread use of unprescribed antibiotics (25,26).  Fig. 1A and B. Each dot represents a metagenome. The blue line represents the regression line of ARd residuals of environmental metagenomes on VFd residuals (with gray shading indicating the 95% confidence interval): slope, 0.438, R 2 ϭ 0.739, and P Ͻ 2.2eϪ16. The red line represents the regression line of ARd residuals of human feces metagenomes on VFd residuals: slope ϭ 0.625, R 2 ϭ 0.8829, and P ϭ 1.753eϪ05. The difference in slopes is significant, but not very significant, with P ϭ 0.03.
The Pearson partial correlation coefficients for metagenomes in the three populations are positive and in the same range (North Americans, 0.66; Malawians, 0.81; Amerindians, 0.70). Regressions of the residuals of ARd on the residuals of VFd (Fig. 4) show a smaller slope for Amerindians (0.20) than for the United States (0.31) and Malawians (0.33), but the differences are not significant: P values range from 0.08 for the difference in Malawian-Amerindian slopes to 0.72 for the difference in Malawian and North American slopes.
Cooccurrence of AR and VFs belonging to the cell envelope. Our results suggest that cooccurrence of AR and VFs might be taking place amid bacterial communities. We wondered, however, which were the genetic traits that were more prone to this effect. We then computed the partial correlations, given the metagenome richness, between each of the AR and VF protein families in the metagenomes under study to generate the correlogram presented in Fig. 5 The vast majority of associations fall into the functional category of multidrug efflux pumps (ARs) associated with either secretion systems, as well with iron uptake and adhesion mechanisms (VFs), respectively (Fig. 5). Among the most representative associations between AR and VF traits are those  Fig. 2A and B. Each dot represents a metagenome. The black line represents the regression of the ARd residuals on VFd residuals for all the metagenomes. The green line represents the linear regression of the ARd residuals on VFd residuals for the USA metagenome subset, the red line is the regression line for the Malawian metagenomes, and the blue line is the regression line for the Venezuelan metagenomes. Black line (all), slope ϭ 0.29015, R 2 ϭ 0.5226, and P Ͻ 2.2eϪ16; green line (USA), slope ϭ 0.3083, R 2 ϭ 0.4332, and P ϭ 1.904eϪ09; red line (Malawi), slope ϭ 0.3324, R 2 ϭ 0.6627, and P ϭ 2.288eϪ06; blue line (Venezuela), slope ϭ 0.2052, R 2 ϭ 0.4882, and P ϭ 0.0004258.
Antibiotic Resistance and Virulence in Microbiomes belonging to the cell envelope and the general secretion mechanisms. On the other hand, the AR protein families for ␤-lactamases and general ␤-lactam resistance mechanisms showed no statistically significant correlations whatsoever.

DISCUSSION
Antibiotic resistance among environmental and human gut metagenomes. As expected, we found different homologues for antibiotic resistance (AR) gene products belonging to different protein families (AR diversity [ARd]) among environmental  (Fig. 1A). We have shown that, in the environmental samples, ARd varies a lot from metagenome to metagenome. This may result from the differential microbial community compositions of the metagenomes, whose genetic diversity can be grouped according to the adaptation to the environment in question (25), but also from the fact that the selective pressure for the maintenance of antibiotic resistance genes in the environmental microbiome varies widely from environment to environment (28)(29)(30). Antibiotic resistance genes can easily travel and spread in the environment due to wind and runoff waters, but also due to human activities and contact with wild animals, in particular, migratory ones (31). The rise of local temperatures due to climate change also has an effect on increasing antibiotic resistance in human bacterial pathogens (32). Natural selective pressures, together with those resulting from such anthropogenic activities as medicine and agriculture, improve the dissemination of resistance genes throughout different environments.
Human gut metagenomes, on the contrary, have a less diversified repertoire of AR determinants (lower ARd) than that of their environmental counterparts.
We found a very strong correlation between the diversity of AR determinants and the metagenomes' protein family richness independently of the geographic origin of the human populations. These similar densities of ARd can indicate that, in human gut microbiomes, the number of different AR genes is not influenced by the human lifestyle, such as diet, medical care, access to antibiotics, or other cultural habits (33), and/or that the adaptation to the intestinal tract shapes microbial AR diversity as well. This is a surprising result given that antibiotic use can vary from country to country (Center for Disease Dynamics, Economics & Policy [CDDEP]; https://cddep.org/), and as a consequence, individuals from different human populations are under different antibiotic exposures. Forslund and coworkers have demonstrated that there are robust differences in the antibiotic resistance arsenals between countries and that these differences follow veterinary and human antibiotic use (34).

Virulence factors among environmental and human gut metagenomes.
In what concerns virulence, we can assert that there are a wide diversity and density differences of VFs in environmental metagenomes, which poses as evidence of the plasticity portrayed by environmental bacteria to adapt to different hosts and niches (Fig. 1B). On the other hand, human gut microbiomes harbor a less diverse VF repertoire, especially in the U.S. samples, which seems to indicate an evolution toward adaptation to the human gut or less contact with pathogens, eventually due to vaccinations and sanitation.
Association of AR and VF. In the present study, we intended to evaluate whether the ability that a bacterial community as a whole has to colonize a host (to encode virulent factors) is linked to its survival under different antibiotic pressures (to encode antibiotic resistance). That is, whether the acquisition of new virulence traits is followed by the acquisition of new resistance traits, or vice versa. However, we are unable to determine the route of dissemination of these characteristics into and within the microbiome. It seems conceivable that the AR and VF genes could travel together in bacterial clones arriving the microbiome and that they could eventually spread inside the microbiome through horizontal gene transfer, but in this study, we have no means to say anything about the mechanisms involved in this spreading.
According to Yatsunenko and colleagues, once the human gut metagenome is established at the age of 3 years, it does not diverge much from individual to individual in terms of both phylogenetic diversity and functional richness (26).
The most relevant results presented here are those that show, when corrected for the metagenomes' protein family richness, ARd and VFd are strongly correlated both in environmental samples (Fig. 3) and in human gut samples (Fig. 4), with special emphasis on the human gut.
The North American intestinal samples (Fig. 4) show a wide variety of associations between ARd and VFd, always presenting a statistically significant correlation between these genetic traits. This result, in itself, reinforces our hypothesis that antibiotic resistance and virulence are in fact coassociated in human gut microbiomes. The U.S. population, like those of other industrialized countries, is culturally exposed to antibiotics from health care facilities such as hospitals, antibiotic therapy, and the use of antibiotics in agriculture and livestock.
Malawians and Amerindians from Venezuela share, phylogenetically, more similar human gut microbiomes than the North Americans (26), which belong mainly to a different enterotype (26,35). Yet, in contrast to the Venezuelan individuals, there is a very strong correlation among the Malawian ones. The uncontacted Amerindian microbiome represents a "frozen" relic of a pre-antibiotic era of the human resistome, while the Malawian gut microbiome is much more exposed, both to antibiotics and to colonization by pathogens (34,36,37).
Amerindians have no known access to pharmaceutical drugs, as they usually make use of the traditional indigenous medicine. Malawi is one of the poorest countries in Africa, where most people live on less than one dollar a day, many people cope with AIDS and bacterial infections, and many children suffer from severe malnutrition (36). Nutrition has been reported to have a big impact on both the human gut microbiome composition and resistome (33,38,39). In Malawi, UNICEF has been implementing a program of ready-to-use therapeutic food (RUTF) to reduce mortality among children. RUTF often contains antibiotics such as co-trimoxazole. It has been questioned, however, whether the success of this therapy is due to renutrition or to the combination with antibiotics (40,41) that could have an effect on the microbial gut composition. It has been reported that there is widespread resistance to almost all of the antibiotics that are empirically used in Malawi due to the lack of routine microbiological culture and sensitivity testing (41), but also due to self-medication. Human gut samples from the Malawi people of the data set used in the present study are described as having a high overall resistance potential, with an overrepresentation of cephalosporin and tetracycline resistance genes, which may suggest extensive use of old, broad-spectrum antibiotics, a known problem in many developing countries (34).
The corepresentation of AR and VF targeting the cell envelope. Between all the possible statistical AR and VF protein families' correlations, the strongest ones are among proteins belonging to the bacterial cell envelope (Fig. 5). This result is not surprising, as many proteins that belong to the secretory system or to the secretome itself, and thus target the bacterial envelope, are frequently encoded on mobile genetic elements (6). As such, this cooccurrence may also be regarded as a direct consequence of the dynamics of these genetic elements coding for both types of determinants.
In human gut, the strongest correlations involve efflux pumps that can extrude antibiotics nonspecifically and adhesion and iron-scavenging mechanisms of pathogenicity. One possible explanation for the fact that the most frequent AR and VF associations in human gut involve efflux pumps could be that they allow a fast and efficient response to new man-made antibiotic molecules, while in environmental genomes, specific resistance mechanisms targeted to specific antibiotic molecules may have had more time to evolve to specific strategies.
Another possible explanation is that extrusion can also be implicated in pathogenic mechanisms of human-associated pathogens as well as bacterium-host interactions. For example, biofilm formation in a Staphylococcus aureus methicillin-resistant (MRSA) strain has been known to be essentially reliant on the activity of fibronectin binding proteins (42), and multidrug efflux pumps have direct implications for the formation and maintenance of such biofilms (43). Furthermore, quinolone-resistant S. aureus strains upregulate the production of fibronectin binding proteins when subjected to a sublethal dosage of ciprofloxacin (44). It has also been acknowledged that physiological levels of some cations present within the host promote the upregulation of genes encoding putative efflux transporters, oxidoreductases, and mechanisms of iron uptake in Acinetobacter baumannii (45) as in Burkholderia cenocepacia (46), which could explain the coassociation of iron acquisition systems with those of multidrug efflux pumps.
We show here that, even after correcting for protein family richness, there is a correlation between the diversity of antibiotic resistance and virulence traits in human gut microbiomes and in environmental microbiomes. Interestingly, this correlation varies from one human population to another. Human antibiotic exposure, due either to therapy or to the environment and food, can have an effect on selecting for potentially pathogenic bacteria in the human gut microbiomes. It can also drive and shape changes in the gene pool of microbiomes. This means that by selecting for resistant bacteria, we might also be selecting for more virulent strains as a side effect of antimicrobial therapy.

MATERIALS AND METHODS
Metagenomic data sets. Our human gut query cohort included 110 previously studied and publicly available metagenomes pertaining to individuals from different regions of Venezuela, Malawi, and the United States, as well as a broad age span (0.05 to 53 years) (see the data set of the project mgp98 "Human gut microbiome viewed across age and geography [WGS]" at the MG-RAST metagenomic analysis server [26]). All of the metagenome files were generated by the same team and project and by using the same bioinformatics pipeline. Our environmental cohort comprised 64 previously selected and publicly available environmental metagenomes belonging to 12 different biomes (see the data set collection of the metagenomes referred to in reference 25). Although Delmont's team report using a data set comprised of 77 metagenomes, there are only 70 MG-RAST accession numbers present in the article's appendix, of which only 64 are publicly available. Both project's metagenomes were downloaded from MG-RAST on 29 January 2016 under the FASTA format, making use of successive calls to its application programming interface (API) (47) using the respective MG-RAST accession numbers available in the aforementioned reference. Each FASTA file comprised clustered protein-coding sequences retrieved from MG-RAST's file-formatting pipeline (550.cluster.aa90.faa files). The protein sequences enclosed in these files are clustered at 90% identity, containing nonredundant translated sequences. These protein-coding formatted FASTA files contain the translation of one representative from each cluster. Thus, the numbers of protein sequences for a given metagenome used here represent its richness.
BLAST, VFDB, and resfams. For every metagenome present in our query, a blastp (50) search was performed against the VFDB database for bacterial virulence factor families (see the virulence factors of pathogenic bacteria VF.pfasta data set of the Institute of Pathogen Biology, China) (27) and the Resfams AR Proteins database of bacterial antibiotic resistance protein families (the Resfams AR Proteins v1.2 database of the Dantas Laboratory at the Washington University School of Medicine, St. Louis, MO) (48). An alternative approach made use of the Antibiotic Resistance Genes Database (ARDB) (reference antibiotic resistance genes data set resisGenes.pfasta of the Center for Bioinformatics and Computational Biology, University of Maryland) (49), but several hindrances concerning its subclassification by functional antibiotic resistance protein families made us discard the possibility of using such a database to address AR and VF protein family cooccurrence in metagenomes.
The BLASTϩ executables package was downloaded on 17 November 2015 (ncbi-blast-2.2.31ϩ version). The VFDB database was downloaded on 11 November 2013 (31 classified FASTA files of bacterial virulence factor subfamilies), and the Resfams database was downloaded on 29 January 2016 (123 classified FASTA files of bacterial antibiotic resistance protein subfamilies). Every protein included in each of the addressed metagenomes was used as a query in order to search for similarities to either AR or VF protein-coding traits. Hence, we aimed at retrieving the best hit (best-scored alignment) that enabled us to assign an AR or VF function to each of the aforementioned metagenomic proteins. Every blastp search was performed with a very stringent E value cutoff of 10 Ϫ15 (10 orders of magnitude lower than the conservative commonly used E value of 10 Ϫ5 ). Next, we filtered the resulting output files so as to retrieve only the alignment hits with Ͼ60% coverage of the query size and whose query and subject length ratio are between 0.75 and 1.5. When using a combination of this very stringent E value of 10 Ϫ15 together with coverage and length ratio filters, we expected to only retrieve true homologues, avoiding false positives. Furthermore, from all the generated alignments between our metagenomic query cohort and the preceding databases, we computed the representative counts for the different gene families that coded for either AR or VF traits. Thus, the amounts of different classes (gene families) that are present in a given metagenome represent its diversity in terms of AR or VF traits.
We have also removed hits for proteins that aligned with both antibiotic resistance and virulence factor proteins (25.9% of hits against the VFDB and 29.9% of the hits against the Resfams database). The above-mentioned filters and algorithms were implemented by making use of Unix scripting languages (GNU Awk version 4.0.1 and Z-Shell version 5.0.2) in a Linux environment.
Statistical analysis. To test for relationships in metagenomes between the presence of both antibiotic resistance and virulence factor traits, we proceeded as follows. Our expectation is that the diversity number for homologues of antibiotic resistance protein families (henceforth denoted as ARd) and the diversity number for homologues of virulence factor protein families (denoted as VFd) in each metagenome (clustered at 90% identity) increases with the protein family richness of the latter (in that the protein family richness here means the total number of cluster representative proteins). Given the potential diversity of these protein families, one does not expect ARd and VFd to level off with a metagenome's protein family richness. Therefore, we assume a linear relationship between ARd and metagenome's protein family richness, with a fixed 0 intercept. The same is assumed for VFd. Our results will show that these assumptions are reasonable. Thus, to avoid spurious correlation, we corrected the diversity of ARd and VFd in a metagenome by taking into account its protein family richness: we have computed the partial correlations between ARd and VFd given the metagenome richness, and equivalently, when plotting ARd versus VFd, we have instead taken residuals of regression of ARd on metagenome richness and residuals of regression of VFd on metagenome richness. Recall that the partial correlation measures the degree of association between two random variables, with the effect of a set of controlling random variables removed. For example, the partial correlation between A and B given a third variable, C, correlated with A and B is the correlation between the residuals e a and e b resulting from the linear regression of A with C and of B with C. If r AB , r AC , and r BC are the correlation coefficients between A and B, A and C, and B and C, the partial correlation between A and B given C, r AB.C , can be computed as

ACKNOWLEDGMENTS
We thank Eduardo P. C. Rocha for critical reading of the manuscript and constant interest and suggestions since the beginning of this work, Guillaume Sapriel for suggestions to improve the methods, Pedro David, Eric Duchaud, Guillaume Achaz, and all the members of our teams for fruitful discussions, and Octávio Paulo, Francisco Pina-Martins, and Sophie Brouillet for technical help and granting access to the servers. We also thank the anonymous referees who helped to improve this manuscript.