Lf, extracted from bovine milk (bLf), possesses a high homology of sequence and similar functions with human Lf (hLf), and it has been approved as a generally recognized as safe (GRAS) compound by the United States Food and Drug Administration (FDA United States) and as a dietary supplement by the European Food Safety Authority.

For in vitro experiments, highly purified bLf was generously given by Armor Proteines Industries (France). BLf was controlled through SDS-PAGE and silver nitrate staining. BLf purity was about 98%, and its concentration was confirmed via UV spectroscopy according to an extinction coefficient of 15.1 (280 nm, 1% solution). The iron saturation of bLf used, determined via optical spectroscopy at 468 nm, was about 7% according to an extinction coefficient of a 1% solution of bLf completely iron saturated corresponding to 0.54. LPS contamination of bLf, assessed via Limulus Amebocyte assay (Pyrochrome kit, PBI International, Italy), was 0.6 ± 0.05 ng/mg of bLf. Before each in vitro assay, bLf solution was sterilized using a 0.2-µm Millex HV filter at low protein retention (Millipore Corp., Bedford, MA, United States).

The African green monkey kidney–derived Vero E6 and human colon carcinoma–derived Caco-2 cells were purchased from American Type Culture Collection (ATCC). Cells were cultivated in high-glucose Dulbecco’s Modified Eagle’s Medium (DMEM) (Euroclone, Milan, Italy) with 10% fetal bovine serum (FBS) (Euroclone, Milan, Italy) at 37°C in humidified incubators with 5% CO₂. SARS-CoV-2 strain, isolated from nasopharyngeal specimen of a positive COVID-19 patient, was propagated in Vero E6 cells. SARS-CoV-2 titers were obtained by 50% tissue culture infectious dose (TCID50) assays in Vero E6 (Spearman-Kärber method) by microscopic scoring. All assays were performed by infecting Vero E6 and Caco-2 cells with SARS-CoV-2 strain in the Department of Molecular Medicine, University of Padua, under Biosafety Level 3 (BSL3) procedures, in agreement with laboratory containment protocols endorsed by the University of Padua.

Cytotoxicity was evaluated by incubating 100 and 500 μg of bLf—the concentrations used for in vitro experiments—in DMEM with 10% of FBS for 72 h at 37°C with Vero E6 and Caco-2 cells in 96-well plates. Cell viability and proliferation were evaluated by MTT assay (Merck, Italy). The MTT assay is colorimetric assay based on the reduction of a tetrazolium salt to formazan by metabolically active cells. The formazan dye was assessed by spectrophotometric absorbance at 600 nm.

For infection assay, Vero E6 cells were seeded in 24-well tissue culture plates at a concentration of 1 × 10⁵ cells/well for 24 h at 37°C in humidified incubators with 5% CO₂, while Caco-2 cells were seeded at a concentration of 2 × 10⁵ cells/well for 48 h at 37°C in humidified incubators with 5% CO₂. In order to evaluate the putative inhibition of SARS-CoV-2 strain infection on Vero E6 monkey cells, 100 μg/ml of bLf was used. Conversely, the supposed antiviral activity against SARS-CoV-2 strain on Caco-2 human cells was investigated using not only 100 but also 500 μg/ml of bLf. In order to investigate the putative interaction of bLf with viral particles and/or host cells, different experimental approaches in both Vero E6 and Caco-2 cells were carried out. To evaluate if bLf can interfere with the viral infectivity rate by binding viral surface components, a multiplicity of infection (MOI) of 0.1 and 0.01 of SARS-CoV-2 was preincubated with bLf for 1 h at 37°C in humidified incubators with 5% CO₂. The cells were then infected with these suspensions for 1 h at 37°C in humidified incubators with 5% CO₂. In order to evaluate if bLf interferes with the viral attachment to host cells, the cells were preincubated in DMEM without FBS with bLf for 1 h at 37°C in humidified incubators with 5% CO₂. The cells were then washed with phosphate buffered saline (PBS) and infected with SARS-CoV-2 at an MOI of 0.1 and 0.01 for 1 h at 37°C in humidified incubators with 5% CO₂. To assess if bLf can interfere with both viral and host cell components, bLf was added together with SARS-CoV-2 at an MOI of 0.1 and 0.01 to cell monolayer for 1 h at 37°C in humidified incubators with 5% CO₂. In addition, the preincubation of SARS-CoV-2 with bLf for 1 h at 37°C was used to infect cell monolayer previously pretreated with bLf for 1 h at 37°C.

Regarding Vero E6 cells, after each experimental approach, the cells were washed with PBS, covered with DMEM containing 0.75% of carboxymethylcellulose and 2% of FBS and incubated for 48 h at 37°C in humidified incubators with 5% CO₂. After 48 h, the cells were washed, fixed with 5% of formaldehyde for 10 min at room temperature, and stained with crystal violet at 1% for 5 min. The number of plaques was determined after extensive washing.

The other infection experiments were carried out with Caco-2 cells. Significant cell death was not observed until 7 days on Caco-2 cells after SARS-CoV-2 infection at MOI 0.1 (Chu et al., 2020). In this respect, after each experimental procedure, the cell monolayers were replaced with DMEM with 2% of FBS, and after 6, 24, and 48 h postinfection (hpi), the supernatant samples were collected for RNA extraction and quantitative real-time reverse transcription (rRT)-PCR assay of viral particles. Briefly, we lysed 200 μl of supernatant in an equal volume of NUCLISENS easyMAG lysis buffer (Biomerieux, France). Viral RNA detection was assayed by in-house real-time RT-PCR in accordance with the protocol and the primers and probes designed by Corman et al. (Corman et al., 2020), targeting the genes encoding the SARS-CoV-2 envelope (E) (E_Sarbeco_F, E_Sarbeco_R and E_Sarbeco_P1). Quantitative rRT-PCR analyses were executed with 5 μl of purified nucleic acids in a final volume of 25 μl, employing One-Step Real-Time kit (Thermo Fisher Scientific) and run on ABI 7900HT Fast Sequence Detection Systems (Thermo Fisher Scientific). Cycle threshold (Ct) data from rRT-PCR tests were carried out for E genes. Genome equivalent copies per ml were inferred according to linear regression performed on calibration standard curves.

The SARS-CoV-2 spike glycoprotein structure in prefusion conformation was extracted from a clustering procedure used in a previously published article (Romeo et al., 2020). The 3D structure of the diferric forms of bLf and hLf, refined at 2.8 and 2.2 Å resolution, respectively, were downloaded from the PDB database (PDB IDs: 1BLF (Moore et al., 1997) and 1B0L (Sun et al., 1999)). The protein-protein docking analysis between the modeled SARS-CoV-2 spike glycoprotein (Romeo et al., 2020) and the Lf structures was carried out using the Frodock docking algorithm (Ramírez-Aportela et al., 2016). Frodock’s approach combines the projection of the interaction terms into 3D grid-based potentials and the binding energy upon complex formation, which is approximated as a correlation function composed of van der Waals, electrostatics, and desolvation potential terms. The interaction-energy minima are identified through a fast and exhaustive rotational docking search combined with a simple translational scanning (Garzon et al., 2009). Both docking procedures were performed using Frodock’s (http://frodock.chaconlab.org/) web-server.

Topology and coordinate files of the input structures have been obtained through the tLeap module of the AmberTools 19 package (Salomon-Ferrer et al., 2013). The spike glycoprotein and Lf were parametrized through the ff19SB force field and were inserted into a rectangular box of TIP3P water molecules, imposing a minimum distance of 12.0 Å from the box walls, while the solution was neutralized adding 0.15 mol/L of NaCl ions. To remove steric interactions, all structures underwent four minimization cycles, each composed by 500 steps of steepest descent minimization followed by 1,500 steps of conjugated gradient minimization. An initial restraint of 20.0 kcal mol⁻¹ Å⁻² was imposed on protein atoms and subsequently reduced and removed in the final minimization cycle. Systems were gradually heated from 0 to 300 K in an NVT ensemble over a period of 2.0 ns using the Langevin thermostat, imposing a starting restraint of 0.5 kcal mol⁻¹ Å⁻² on each atom, which was decreased every 500 ps in order to slowly relax the system. The systems were simulated in an isobaric-isothermal (NPT) ensemble for 2.0 ns, fixing a pressure of 1.0 atm using the Langevin barostat and imposing the temperature at 300 K. Covalent bonds involving hydrogen atoms were constrained using the SHAKE algorithm (Ryckaert et al., 1977). 30 ns of production run were performed through the NAMD 2.13 MD package (Phillips et al., 2005), using a time step of 2.0 fs. The PME method was applied to take into account long-range interactions, while a cutoff of 9.0 Å was set for short-range interactions. System coordinates were saved every 1,000 steps.

Distance analysis was carried out using the distance module of the GROMACS 2020 analysis tools (Abraham et al., 2015), while hydrogen bond persistence was evaluated using the hbonds module coupled to in-house written codes. The hydrophobic contacts were identified using the contact_map and contact_frequency routines of the mdtraj Python library (McGibbon et al., 2015). Generalized Born and surface area continuum solvation (MM/GBSA) analyses were carried out using the last 15 ns of the trajectories, through the MMPBSA. py.MPI program as implemented in the AMBER16 software (Case et al., 2016) on two nodes of the ENEA HPC cluster CRESCO6 (Iannone et al., 2019). Pictures of the spike-Lf and spike CTD1-ACE2 complexes were generated using the UCSF Chimera program (Pettersen et al., 2004).

For in vitro experiments, the number of plaque-forming units (pfu)/ml of SARS-CoV-2 on Vero E6 cells and the number of SARS-CoV-2 RNA copies/ml on Caco-2 cells in each experimental approach was compared to the control ones (untreated SARS-CoV-2 and cells) at the same time point in order to assess the statistically significant differences by using unpaired student’s t tests. Results were expressed as the mean values ± standard deviation (SD) of three independent experiments. In each case, a p value ≤0.05 was considered statistically significant.

Preliminarily, the doses of bLf in native form (7% iron saturated) corresponding to 100 μg/ml for Vero E6 cells and 100 and 500 μg/ml for Caco-2 cells were assayed to detect their putative cytotoxicity by measuring cell morphology, proliferation, and viability after 72 h of incubation. Both 100 and 500 μg/ml of bLf did not exert any cytotoxic effect (data not shown).

Then, the efficacy of different concentrations of bLf in the inhibition of SARS-CoV-2 infection was tested on Vero E6 and Caco-2 cells according to different experimental procedures: I) control: untreated SARS-CoV-2 and cells; II) preincubation of bLf with virus inoculum for 1 h at 37°C before cell infection; III) preincubation of bLf with cells for 1 h at 37°C before virus infection; IV) bLf added together with SARS-CoV-2 at the moment of infection; and V) virus and cells separately preincubated with bLf for 1 h at 37°C before infection.

The results obtained with Vero E6 cells are shown in Figure 1A (MOI 0.1) and 1B (MOI 0.01).

Regarding Vero E6 cells, an inhibition of SARS-CoV-2 replication of about one log at MOI 0.1 and about two log at MOI 0.01 was observed when 100 μg/ml of bLf was preincubated for 1 h with virus before infection compared to untreated SARS-CoV-2 infection (p < 0.001 and p < 0.0001, respectively) (Figures 1A,B).

On the contrary, the data illustrated in Figures 1A,B, independently from the MOI used, indicate that bLf, at this concentration, does not block SARS-CoV-2 infection when it is preincubated with Vero E6 cells or when bLf is contemporary added to viral particles and cells at the moment of infection (Figures 1A,B). BLf is also ineffective when it is preincubated for 1 h at 37°C separately with virus and cells before infection (Figures 1A,B).

The efficacy of 100 and 500 μg/ml of bLf against SARS-CoV-2, assayed in Caco-2 cells, is showed in Figures 2A,B (MOI 0.1) and in Figures 2C,D (MOI 0.01), respectively.

Regarding Caco-2 cells, at MOI 0.1, no significant differences were observed in all experimental conditions compared to the control ones when using bLf at 100 μg/ml (Figure 2A). At MOI 0.01, an inhibition of viral load in supernatants was observed at 24 hpi only when 100 μg/ml of bLf was preincubated with the viral inoculum and when the cells were preincubated with 100 μg/ml of bLf compared to the control one (p < 0.05) (Figure 2B). At 48 hpi, an inhibition of viral load was observed only when the cells were preincubated with bLf (p < 0.05) (Figure 2B).

When bLf was used at a concentration of 500 μg/ml, a decrease of viral load up to 48 hpi was observed when the viral inoculum was preincubated with bLf compared to the control group, independently from the MOI used (p < 0.05) (Figures 2C,D). When the cells were preincubated with bLf, a decrease of SARS-CoV-2 load up to 24 hpi was observed compared to the control at MOI 0.1 (p < 0.001 after 6 hpi and p < 0.05 after 24 hpi) (Figure 2C), while at MOI 0.01 the decrease of viral load remained statistically significant up to 48 hpi compared to the control group (p < 0.05) (Figure 2D). When bLf was added together with SARS-CoV-2 during the adsorption step a decrease of viral load up to 24 hpi was observed compared to untreated SARS-CoV-2 infection, independently from the MOI used (p < 0.001 after 6 hpi and p < 0.05 after 24 hpi for MOI 0.1; p < 0.05 after 6 and 24 hpi for MOI 0.01) (Figures 2C,D). When the cells were preincubated with bLf and infected with SARS-CoV-2 previously preincubated with bLf, a decrease of viral load up to 24 hpi was observed for MOI 0.1 compared to untreated SARS-CoV-2 infection (p < 0.001 after 6 hpi and p < 0.05 after 24 hpi for MOI 0.1) (Figure 2C), while at MOI 0.01 the decrease of viral load remained statistically significant up to 48 hpi compared to untreated SARS-CoV-2 infection (p < 0.05) (Figure 2D).

The molecular docking simulation suggests a potential interaction of the bLf structure with the spike glycoprotein CTD1 domain in the up conformation (Figure 3A). The first three solutions obtained by Frodock clustering procedure account for more than 60% of the total generated complexes and are almost completely superimposable to that shown in Figure 3A. Starting from the first Frodock solution, we performed a 30 ns long classical MD simulation in order to verify the stability of the complex and check for the presence of persistent interactions between the two proteins. As shown in Supplementary Figure S1A, the distance between the centers of mass of spike and bLf, calculated as a function of time, oscillates around the value of 4.5 nm, indicating a constant close contact between the two molecules. MM/GBSA analysis confirmed the high affinity of the bLf for the spike CTD1 domain (Supplementary Table S1A), showing interaction energy of −28.02 kcal/mol. In particular, MM/GBSA results underlined that the Van der Waals term mainly contribute to the binding energy (Supplementary Table S1A).

A detailed analysis of the interaction network revealed the presence of 28 different interactions, which persist for more than 25% of the simulation time, in agreement with the high interaction energy calculated. In detail, we found three salt bridges, five hydrogen bonds and 20 residue pairs involved in hydrophobic contacts (Supplementary Table S2 left side).

To check if some of the spike residues targeted by the bLf protein were involved in the binding with ACE2, we compared the average structure extracted from the simulation with the ACE2/CTD1 domain complex structure (PDB ID: 6LZG (Wang et al., 2020)) (Figure 4). Surprisingly, only two spike residues (Gly502 and Tyr505) were shared between the complexes interfaces (Supplementary Table S2 left side), as evaluated from the inspection of the superimposed structures and from the article analysis (Wang et al., 2020). Despite this, Lf holds the same position assumed by the ACE2 enzyme, that is, above the up CTD1 domain.

We performed the same analysis over the evaluated hLf–spike complex, obtaining a binding pose superimposable to that observed for the bovine protein (Figure 3B). Although using the human protein we could still observe a persistent and close contact between the two molecules (Supplementary Figure S1B), the analysis of the interaction network indicated the presence of a larger number of interactions (45 interactions), in agreement with the higher interaction energy revealed by the MM/GBSA approach (−48.25 kcal/mol, Supplementary Table S1B). In detail, we found 12 salt bridges, 10 hydrogen bonds and 23 residue pairs involved in hydrophobic contacts (Supplementary Table S2 right side), in agreement with the presence of a negative electrostatic contribution term (Supplementary Table S1B). Comparing the average structure extracted from the simulation with the ACE2/CTD1 domain complex structure (PDB ID: 6LZG (Wang et al., 2020)) (Supplementary Figure S2), we observed that also for the hLf, only two residues (Thr500 and Tyr505) were shared between the complexes interfaces (Supplementary Table S2 right side).

These results allow us to hypothesize that, in addition to the HSPGs binding (Lang et al., 2011), both bLf and hLf should be able to hinder the spike glycoprotein recognition of the ACE2 receptor, blocking the virus from entering into the cells.