The SARS-CoV-2 Mu variant should not be left aside: It warrants attention for its immuno-escaping ability.

The COVID-19 pandemic continues to have a threatening impact on a global scale, largely due to the emergence of newly SARS-CoV-2 variants. The Mu (PANGO lineage B.1.621), was first identified in Colombia in January 2021 and was classified as a variant of interest (VOI) in August 2021, due to a constellation of mutations that likely-mediate an unexpectedly enhanced immune resistance to inactivated vaccine-elicited antibodies. Despite recent studies suggesting that the Mu variant appears to have less infectivity than the Delta variant, here we examined the structural effect of the Mu spike protein mutations and predicted the potential impact on infectivity of the Mu variant compared with the Delta and Delta plus spike protein.

INTRODUCTION

The COVID‐19 pandemic continues to have a threatening impact on a global scale, largely due to the newly emerging SARS‐CoV‐2 variants. Among the most recent variants, the Mu (PANGO lineage B.1.621), as designed by WHO, was first identified in Colombia in January 2021 and was classified as a variant of interest (VOI) in August 2021 (https://www.who.int/en/activities/tracking-SARS-CoV-2-variants/).

The Mu variant showed an unexpectedly enhanced immune resistance to inactivated vaccine‐elicited antibodies. On the other side, the Mu variant demonstrated less infectivity than the Delta variant predominant in the world at the time of its emergence. This observation suggests a biological trade‐off between viral transmission and immune escape.

After being identified in Colombia in January 2021, the presence of the SARS‐CoV‐2 Mu variant was sporadically documented, with a prevalence of 0.1% (https://unric.org/en/covid-19-what-is-the-mu-variant/). Although the Mu strain is not yet classified as a “variant of concern” by WHO (which means it is being more closely monitored than a VOI), the mutations found in the B.1.621 strain may have “a synergistic impact on attributes such as reduction of vaccine‐induced protection from severe disease, increased transmission, and disease severity.”

Due to these important epidemiological properties, the Mu variant is now responsible for 39% and 13% of infections in Colombia and Ecuador, respectively. Subsequently, in September 2021, the Mu variant has been detected in more than 20 countries.

Mutations arise as a natural by‐product of viral replication. RNA viruses typically have higher mutation rates than DNA viruses. Coronaviruses, however, make fewer mutations than most RNA viruses because they encode an enzyme that corrects some of the errors made during replication. In most cases, the fate of a newly arising mutation is determined by natural selection. Those that confer a competitive advantage with respect to viral replication, transmission, or escape from immunity will increase in frequency, and those that reduce viral fitness tend to be culled from the population of circulating viruses. However, mutations can also increase and decrease in frequency due to chance events.

The Mu variant spike protein carries eight mutations, including T95I, Y144S, Y145N, R346K, E484K, N501Y, D614G, P681H, and D950N. Several of them are also present in other variants such as the E484K in Beta and Gamma, N501Y in Alpha and Beta, P681H in Alpha, and D950N in Delta. (https://www.ecdc.europa.eu/en/covid-19).

There are currently very few studies looking into the Mu variant. As research on this specific variant is still at a preliminary stage, it is too early to conclude whether the strain is more transmissible or more severe than other variants. However, the spike protein mutations seem to facilitate viral entry of the Mu variant into ACE2‐expressing cells.

We examined the structural effect of the Mu spike protein mutations and predicted the potential impact on infectivity of the Mu variant compared with the Delta and Delta plus spike protein. In fact, some scientists have warned that the Mu strain could become a threat if it were to become the dominant variant. This is due to its ability to spread more easily through cells in the human body, something that is also seen in the Delta variant.

As a matter of fact, at the time this manuscript is being written, the new omicron variant seems to be a candidate to become one of the worldwide dominant variants. Nonetheless, studying the structural properties of the Spike Mu variant can provide a general and deeper understanding of the evolution of the SARS‐CoV‐2 and of the relationships between virus protein structure, infectivity, and interaction with the immune system (as in the case of lambda variant ).

MATERIALS AND METHODS

Spike variant sequences were retrieved from the GISAID translated protein data set. The reference Wuhan Spike sequence is labeled by the RefSeq code yp_009724390. Representative isolates of the Spike sequences of the Mu, Delta, and Delta plus variants are identified by the GISAID codes EPI_ISL_2178393, EPI_ISL_1634920, and GenBank UEN64961, respectively.

Variant spike RBD‐ACE2 and antibody complexes have been modeled by in‐silico mutagenesis using the ad hoc tools in PyMOL. Energy minimization protocol embedded in the molecular graphics program Swiss‐PdbViewer has been applied to remove residue steric overlaps at the interface. The protocol used the GROMOS96 43B1 force field, cutoff 10 Å, and 100 steps of steepest descent minimization followed by 1000 steps of conjugate gradients in vacuo. The minimization was stopped if the energy difference between two consecutive steps was lower than 0.05 kJ/mol. Only residues at the interface have been minimized. This forcefield does not include the parameters to describe glycans which were therefore ignored during minimization. However, the glycans present in the RBD and NTD structures are apparently far from the complex interface and do have not the potential to interfere with it.

The noncovalent interactions taking place at the interface of the predicted complexes were identified with the RING 2.0 web server. The server provides a fast tool to detect intra and interchain interactions including solvent and ions. Protein structural analysis and visualization have been carried out with PyMOL. Sequence alignment and visualization relied on MAFFT and Jalview, respectively.

Computational alanine scanning of the interface residues of the Spike complexes was carried out through the webserver DrugScorePPI. The server provides a fast and accurate tool to predict the binding free energy changes upon alanine mutations at protein‐protein interfaces using a knowledge‐based scoring function.

Molecular dynamics was applied to study the complex RBD‐ACE2. The starting structure of the reference form of the spike RBD bound to ACE2 was taken from the Protein Data Bank (PDB code: 6M0J), while the structures of the Mu, Delta, and Delta plus variants were modeled starting from the reference form as described above. Protonation states of titratable protein residues were determined based on empirical pKa prediction performed with the PropKa program. The position of the proton on the neutral histidine residues was set by visual inspection of the structure. The systems were solvated in a dodecahedron box, imposing a minimum distance of 1.6 nm between the protein and the box, using the TIP3p water model. Na+ and Cl− ions were added to the systems at a concentration of 150 mM with an excess of Na+ ions to compensate for the net negative charge of the solute. Protein residues were described by the Amber03 force field. The simulated systems were composed of about 230 000 atoms each.

Long‐range electrostatic interactions were treated using the particle mesh Ewald (PME) method, with a short‐range cutoff of 1.2 nm and a grid spacing of 0.12 nm. The bond lengths of the hydrogen atoms were constrained to a constant value applying the LINCS algorithm. The temperature was kept constant at 310 K using the Nose–Hoover algorithm with a coupling time constant τ = 0.4 ps. The system was coupled to a pressure bath at 1 bar with τ = 1.0 ps, using the Parrinello–Rahman barostat. A time step of 2 fs was used for numerical integration of the equations of motion.

All equilibrium MD simulations were preceded by energy minimization using the steepest descent algorithm followed by 20 ns of MD simulation at constant volume and temperature (310 K) with harmonic position restraints (force constant 1000 kJ mol−1  nm−2) applied on the heavy atoms of the protein residues present in the original X‐ray structure. Finally, unrestrained MD simulations were carried out for 130 ns for each system. MD simulations were carried out using the GROMACS software package.

RBD‐ACE2 binding energies were estimated as the sum of the short‐ranged Lennard‐Jones (LJ) interactions and short‐range Coulomb interactions between the two groups consisting of the spike RBD and the ACE2 domain as implemented in the function energygrps of the Gromacs package. Binding energies were calculated averaging over 100 ns of unrestrained MD simulations (i.e., neglecting the first 30 ns of each simulation).

RESULTS

RBD‐ACE2 complexes

Interactions occurring at the RBD‐ACE2 interfaces in Wuhan and Mu, Delta, and Delta plus variants were compared. RING 2.0 analysis predicts that mutations characterizing the RBD of the variants alter the interaction patterns. Mu mutation at site 501 of RBD adds a π–π stacking interaction between ACE2 Y41 and RBD Y501 (Figure 1A). The Y501 also interacts with ACE2 K353 through a Van der Waals contact (Table 1). Increased interaction stability induced by the presence of interface Y501 was already described. Mutations in the Delta variant do not directly involve interface interactions; however, it has been suggested that the increased positive surface electrostatic potential induced by the replacing basic residues may stabilize interaction with the negatively charged ACE2 surface. ,  A similar effect on the electrostatic surface potential can be induced by the Mu mutation E484K. The Delta plus variant possesses one more mutation with respect to the parent Delta, namely K417N (Table 1). Apparently, this mutation removes a salt bridge that in the reference Spike connects K417 and ACE2 D30 (Figure 1B).

Figure 1

(A) Complex between RBD and ACE2 from the Protein Data Bank (PDB) structure 6M0J. Green and grey cartoons indicate ACE2 and RBD, respectively. (B) Top: comparison of the interfaces between the reference and the model Mu RBDs and ACE2 at position 501. Middle: possible salt bridge between E75 and K484 in the Mu variant and comparison with the reference domain. Bottom: comparison of the interfaces between the reference and the model Delta plus RBDs and ACE2 at the position 417. Residue side chains and NAG substituent are reported as stick models

Table 1

Predicted molecular interaction differences between the reference and variant RBDs and ACE2

Reference residuea	Mua		Delta		Delta plusa
R346	K	Exposed; not at interface	‐	‐	‐	‐
K417 (0.29)	‐	‐	‐	‐	N	Removal of the salt bridge between K417 and ACE2 D30 (0.16)
L452	‐	‐	R	Exposed: near the interface	R	As in Delta
T478	‐	‐	K	Exposed: near the interface; possible salt bridge with ACE2 E87	K	As in Delta
E484	K	Near the interface; possible salt bridge with ACE2 E75	‐	‐	‐	‐
N501 (0.57)	Y	Van der Waals interaction with ACE2 Y51 and K353 (2.46)	‐	‐	‐	‐

Note: Mutant residues characteristic of each variant are compared.

DrugscorePPI ΔΔGs of interface residues are reported in parenthesis as kcal/mol.

The contribution of the residues to the stability of the RBD‐ACE2 interface has been predicted with the application of DrugScorePPI in‐silico alanine scanning. Results are reported in Table 1. Mu Y501 is an interface hot spot as it contributes 2.46 kcal/mol to the interaction energy while the reference N501 is only 0.57 kcal/mol. Somehow surprisingly, the same analysis predicted that the substitution K417N of the Delta plus variant has a minimal impact on the interface stability.

To estimate the effect of the mutations of the three investigated variants on the interaction energy between RDB and ACE2, molecular dynamics simulations of the RBD‐ACE2 complexes were carried out for the reference RDB as well as the three variants Mu, Delta, and Delta plus. The energy contributions due to the Coulomb and LJ interactions were estimated for the three variants averaging the values along the respective simulated trajectories and compared with the values obtained from the simulation of the reference form. The calculated values are shown in Table 2. Considering the total RDB‐ACE2 interaction energy, given by the sum of the Coulomb and LJ contributions, the lowest value was found for the reference system. Nevertheless, the total interaction energy calculated for the Mu and Delta variants is only 2.04 and 3.64 kcal/mol less stable than the reference system, while the Delta plus variant apparently shows the least stable interaction (ΔΔE = + 22.4 kcal/mol). This behavior is not unexpected if we consider that, among the four simulated systems, the only structure experimentally determined is the reference structure, while the three variants have been modeled from it. This aspect could be particularly relevant in the case of the Mu variant. Indeed, the residue K484 located on a loop of RDB, could interact with the residue E75 of the ACE2 protein, thus potentially creating a salt bridge with a consequent decrease of the Coulombic energy (Figure 1B). Indeed, along 130 ns of simulation, this salt bridge was found to be transiently formed. The formation of the salt bridge could significantly affect the Coulombic interaction energy of the Mu variant, calculated here to be 8.42 kcal/mol, thus less stable than the reference system, and possibly making the total interaction energy of the Mu variant more stable. Similar considerations are valid in the case of Delta and Delta plus variants, where an increased sampling could turn out into a stable interaction between K478 of RDB and E87 of ACE2. However, in the case of Delta plus, the huge destabilization of the RDB‐ACE2 affinity, largely due to the Coulombic interaction (Table 2, +20.24 kcal/mol), could be hardly compensated even by the formation of a salt bridge. The major reason for this destabilization should be the loss of the salt bridge K417‐D30 of the reference RBD‐ACE2 complex even if DrugScorePPI fails to detect any clear perturbation upon mutation of K417 (Table 1).

Table 2

Predicted molecular interaction differences (∆∆E) between the reference and variant RBDs and ACE2 as obtained from MD simulations

	SR Coulomb	SR Lennard‐Jones	Total
Reference	‐	‐	‐
Mu	+8.42	−6.38	+2.04
Delta	−2.58	+6.22	+3.64
Delta Plus	+20.24	+2.16	+22.40

Note: SR Coulomb and Lennard‐Jones contributions are explicitly reported. Values reported as kcal/mol, are averaged over the last 100 ns of the simulated equilibrium trajectories.

Abbreviation: SR, short‐ranged.

RBD‐antibody complexes

Effect of Mu, Delta, and Delta plus mutations on RBD interaction with neutralizing antibodies has been considered in two study cases: interaction with the antibodies COVA2‐04 and COVA2‐39 corresponding to the complexes reported in the PDB coordinate sets 7JMO and 7JMP, respectively. The two antibodies, coded by the gene IGHV3‐53, differ from each other mainly for the length of the CDR H3 (namely, Complementarity‐Determining Region 3 of Heavy chain) that is 10 and 15 residues long in COVA2‐04 and COVA2‐39, respectively. The two antibodies interact with the ACE2‐binding site on RBD using different modes, dominated by heavy chains, characterized by different residue contact patterns.

The Mu, Delta, Delta plus RBD complexes with the two antibodies were modeled via in‐silico mutagenesis starting from the PDB structures 7JMO and 7JMP. Swiss‐PdbViewer energy minimization was applied to correct possible structural inconsistencies as described in the materials and methods section.

In the interaction with COVA2‐04 (reported in the PDB structure 7JMO), the Mu RBD mutation Y501 may form an additional H‐bond interaction with the antibody residue S29 of the Light chain CDR L1. This observation has been corroborated by the DrugScorePPIin‐silico alanine scanning that assigns to Y501 a higher ΔΔG than that of the corresponding reference residue (Table 3). Moreover, the mutation E484K may create a local patch of positive surface electrostatic potential complementary to a negatively charged region in COVA2‐04 (Figure S2). Overall, these results suggest that Mu RBD can bind this antibody with a higher affinity.

Table 3

Predicted differences between variants and reference Spike domain in interaction with antibodies

Variant	Complexed antibody	Mutation		DrugScorePPIΔΔGa (kcal/mol)
Mu	COVA2‐04 (7JMO)	N501Y	H‐bond to S29 of L1 CDR	2.46 (0.57)
Delta plus	COVA2‐04 (7JMO)	K417N	Loss of K417‐E97 salt bridge (CDR H3)	0.52 (1.13)
			Loss of K417‐Y52 H‐bond (CDR H2)	‐
Mu	COVA2‐39 (7JMP)	E484K	Loss of H‐bond E484‐NH‐G54 (CDR H2)	‐
Mu	4A8 (7C2L)	Y145N	Loss of VdW interaction to A101, V102, and Y111 (CDR H3)	0.01 (1.98)

ΔΔG of the reference residue is reported in parentheses.

The Delta variant RBD displays an interaction pattern with COVA2‐04 similar to that of the reference Spike. However, it has been suggested that the mutations occurring in the RBD may increase the positive surface potential that can influence interaction with the negative patches of the antibody (Figure S2).

The K417N mutation that characterizes the Delta plus RBD, removes the reference salt bridge formed between K417 and E97 of the antibody CDR H3 and Y52 of CDR H2. DrugScorePPI confirms that the energy contribution of N417 to the stability of the complex is about half with respect to K417 in the reference RBD (Table 3). Overall, these results suggest a possible decrease of the affinity for the antibody.

Interaction to antibody COVA2‐39 (PDB code 7JMP) displays a similar pattern in the case of reference, Delta, and Delta plus variants. In fact, the antibody interacts with the RBD regions not affected by mutations. On the contrary, the mutation E484K on the Mu RBD induces disruption of the hydrogen bond between the carboxylic group of E484 and the NH of the peptide bond of G54 of antibody CDR H2. In this case, also, DrugScorePPI is not able to detect any significant variation of the contribution of residue 484 to the interface stability. The most evident effect of the E484K mutation is on the electrostatic potential surface that becomes more positive (Figure S2).

NTD‐antibody complexes

Effect of the Mu, Delta, and Delta plus mutations on NTD interaction with antibodies has also been modeled in two study cases: NTD–4A8 and NTD–5‐7 corresponding to the PDB codes 7C2L and 7RW2, respectively. The NTD of Delta and Delta plus variants has been modeled onto the templates of the reference NTDs contained in the selected PDB complexes because of the presence of the deletion of the positions 157–158 (Figure S1). The Mu NTD has been created by in‐silico mutagenesis starting from the reference structures of the selected complexes. Mutations occurring within Mu, Delta, and Delta plus NTDs do not change the residue charge and, as a consequence, modification of the electrostatic surface potential is not expected to be relevant for the interaction with the antibodies (Figure 2).

Figure 2

(A) Complex RBD‐COVA2‐04 reported in the Protein Data Bank (PDB) structure 7JMO. Red and blue cartoons indicate the Ab heavy and light chain, respectively, while the RBD is colored grey. (B) Top: comparison of the interfaces between the reference and the model Mu RBD and the mAb COVA2‐04 at position 501. Bottom: comparison of the interfaces between the reference and the Delta plus RBD and the mAb COVA2‐04 at position 417. Residue side chains and NAG substituent are reported as stick models

Interface interactions in the complexes with 4A8 antibody (PDB 7C2L) display the same pattern between the reference, Delta, and Delta plus NTDs. Indeed, the interacting interfaces are free from mutations. On the contrary, the Mu Y145N mutation in NTD removes of the Van der Waal interaction with antibodies A101, V102, and Y111 located in the CDR H3 (Figure 3). This observation suggests that removal of interfacial atomic contacts may decrease affinity between these two proteins. Details are reported in Table 3. DrugScorePPIin‐silico alanine scanning supports these considerations: contribution to the interface interaction energy is predicted to be 0.01 and 1.98 kcal/mol for N and Y145, respectively.

Figure 3

(A) Complex NTD‐4A8 reported in the Protein Data Bank (PDB) structure 7C2L. Red and blue cartoons indicate the Ab heavy (H) and light chain (L), respectively, while the NTD is colored grey. (B) Comparison of the interfaces between the reference (top) and the model Mu (bottom) NTD and the mAb 4A8

Mu, Delta, and Delta plus mutations do not appear to significantly influence the interaction to 5‐7 antibody. In this case, also, mutations do not alter the NTD interface to the antibody (results not shown).

DISCUSSION

Over the course of the severe acute respiratory syndrome coronavirus 2 pandemic, the clinical, scientific, and public health communities have had to respond to new viral genetic variants. Each one has triggered a flurry of media attention, a range of reactions from the scientific community, and calls from governments to either “stay calm” or pursue immediate countermeasures. While many scientists were initially skeptical about the significance of the D614G alteration, the emergence of the new “UK variant” (lineage B.1.1.7) has raised widespread concern. Understanding which variants are concerning, and why, requires an appreciation of virus evolution and the genomic epidemiology of SARS‐CoV‐2.

Mu was first identified in Colombia in January 2021, and outbreaks have since been “sporadically reported” in South America and Europe, WHO said.  Although the global prevalence of Mu is less than 0.1% among sequenced COVID‐19 cases, its prevalence has “steadily increased” in Colombia and Ecuador, where it now accounts for about 39% and 13% of infections, respectively. Up to now the Mu strain is not yet classified as a “variant of concern” by WHO, the mutations found could have a sort of synergistic impact on characteristics such as reduction of vaccine‐induced protection from severe disease, increased transmission, and disease severity.

Compared to the D614G variant, the Mu variant also displayed a lower viral transmission than the Delta variant, when the effect of the Mu spike protein mutations on cell‐cell fusion was examined. A study by Xie et al. suggested that Mu is likely to be less transmissible than Delta, although more evidence is needed.

Messali et al. carried out a neutralization assay using the isolated virus from a patient soon after confirmation by WGS. Sera samples efficiently neutralized the SARS‐CoV‐2 Mu variant isolate, demonstrating that this variant could be not a concern for vaccine effectiveness. Usually, viruses including SARS‐CoV‐2 change over time. Most changes have little to no impact on the virus's biological properties. However, some changes may affect some critical ones, such as how easily it spreads, disease severity, or the performance of vaccines, therapeutic, and diagnostic tools.

However, there is insufficient evidence to say conclusively whether and to what extent the Mu variant will be able to evade protection from coronavirus vaccines, that is, a principal concern for this variant. The Japan study, using serum samples from people who had received the primary vaccination, and collected one week after their second dose, showed that sera, although reduced in effectiveness by around 70%–80% were still able to neutralize the Mu variant. Rapid genome sequencing in hospitals to identify potential vaccine‐escape SARS‐CoV‐2 variants. However, a study published by Snell et al. suggested that the Mu variant could immuno‐escape from vaccine protection. The authors reported spike mutations within Mu able to reduce neutralization by antibodies. The Mu variant was also found to have the same spike mutation that has been associated with a weakened vaccine response to the Beta variant.

To contribute to the molecular characterization of the Mu variant, we compared the predicted effect of the Mu Spike mutations with those of the Delta and Delta plus variants in relation to the reference Wuhan Spike. In particular, the impact of the mutations of the RBD on the interaction with ACE2 and with two neutralizing antibodies has been considered. Moreover, also the effect of the NTD mutations on the interaction with two different antibodies was compared.

The molecular dynamics experiment suggests that in the interaction with ACE2, Delta relies mainly on electrostatic interaction while Mu prefers Van der Waals stabilization (Table 2) even if the overall effect on interaction stability is similar. This is coherent with the significant contribution of the Mu Y501 mutation to the interface stability as predicted by DrugScorePPI. Interestingly, Delta plus appears to have the least stable interaction with ACE2.

Similar observations can be made about the interaction between RBDs and the two selected antibodies. With the antibody COVA2‐04, Delta appears to have a mode of interaction similar to the reference Spike. In the case of Mu variant, however, the substitution N501Y can add further stabilization resulting in a stronger interaction. In this case also, Delta plus shows a predicted weaker interaction caused mainly by the loss of electrostatic interactions (Table 3).

The examination of the antibody COVA2‐39 suggests that Mu has a potentially weaker interaction to it compared to the reference virus strain Spike. On the contrary, the Delta and Delta plus variants display a reference‐like interaction pattern.

Analysis of the interactions with the Abs involving the Spike NTD suggests that Mu has a potentially weaker interaction with 4A8 compared to the reference, Delta and Delta plus variants. Indeed, Mu mutations remove several Van der Waals bonds (Table 3). Delta and Delta plus NTDs are not hit by any mutations and behave like the reference NTD.

All the variants and the reference Spike NTDs have similar interactions with the 5‐7 antibody.

In conclusion, Mu and Delta Spike are predicted to have a similar affinity for ACE2 although Delta privileges electrostatic binding. Delta plus is the weaker binder. In relation to the interaction with antibodies, Mu has higher or lower predicted affinities for the different antibodies selected for this study. In a sense, this is coherent with the experimental evidence and suggests the presence of potential Mu immuno‐escape ability.

Considering the high immune resistance of Mu to inactivated and mRNA vaccines (reference) global epidemiological monitoring of this variant is required. Anyway, our study strongly calls for emergent (urgent?) evaluation of the protective efficacy of current COVID‐19 vaccines against the Mu variant. Indeed, scientific theories are not perfect replicas of reality, but we have good reason to believe that they capture significant elements of it. And experience reminds us that when we ignore reality, it sooner or later comes back to bite us.

CONFLICT OF INTERESTS

The authors declare no conflict of interest.

AUTHOR CONTRIBUTIONS

Conceptualization: Stefano Pascarella, Antonio Cassone, and Massimo Ciccozzi; Performed experiments: Stefano Pascarella, Martina Bianchi, and Daniele Narzi; Validation: Stefano Pascarella, Martina Bianchi, Marta Giovanetti, Daniele Narzi Roberto Cauda, Antonio Cassone, and Massimo Ciccozzi; Helped with study design and data interpretation: Stefano Pascarella, Martina Bianchi, and Daniele Narzi; Validation: Stefano Pascarella, Martina Bianchi, Marta Giovanetti, Daniele Narzi Roberto Cauda, Antonio Cassone, and Massimo Ciccozzi; Wrote the initial manuscript, which was reviewed by all authors: Stefano Pascarella, Martina Bianchi, and Daniele Narzi; Validation: Stefano Pascarella, Martina Bianchi, Marta Giovanetti, Daniele Narzi Roberto Cauda, Antonio Cassone, and Massimo Ciccozzi.

Supporting information.

Click here for additional data file.