Main

As of early June 2021, COVID-19 has infected more than 1.7 billion people worldwide, with global deaths attributed to COVID-19 exceeding 3.7 million according to the World Health Organization (WHO) COVID-19 Dashboard (https://covid19.who.int/). Although COVID-19 is primarily an infectious disease caused by the SARS-CoV-2, the infection can trigger a myriad of metabolic perturbations that may alter the trajectory of the course of the disease and the final outcome. The temporal course of COVID-19 pathogenesis is marked by direct viral-mediated effects on host metabolism, as well as collateral damage, which results from the response of the host to the pathogen4. Originally perceived as a means for cells to excrete unwanted materials5, later studies have implicated exosomes in intercellular communication, particularly with respect to viral transmission6,7 and immune responses8. As lipid-bilayer-enclosed structures, EVs exhibit varying lipid membrane compositions depending on the nature and physiological status of their parent cells9. Severity-associated alterations in membrane lipid composition were reported for exosome-enriched EVs in COVID-19 (ref. 10). As exosomes originate primarily from within the MVBs that concentrate more than half of the membrane cholesterol along the endocytic pathway11, cholesterol homeostasis elicits important effects on exosome biology and functions12. Cholesterol homeostasis is identified as a key host pathway that mediates SARS-CoV-2 infection of human cells13,14. Cholesterol with sphingomyelin (SM) also forms ordered raft membrane microdomains that fine-tune membrane fluidity and alter protein localization15. For example, γ-secretase (GS), which cleaves NOTCH-1 to form Notch intracellular domains16 and β-amyloid precursor protein (APP) to form β-amyloid peptides (Aβ)17, displays preferential localization and enhanced activity in lipid rafts18. The lipid composition of EV membranes therefore can potentially alter the localization and activity of GS, modulating downstream pathways pathologically relevant to COVID-19.

We allocated plasma samples to four temporal stages on the basis of the sample collection date relative to the dates of hospital admission and discharge, SARS-CoV-2 RNA test results and associated laboratory findings (Table 1 and Supplementary Table 1; see Methods section ‘Study design and participants’ for details). No significant differences in baseline characteristics (age, sex and body-mass index) were observed across the four disease stages (Table 1). Slight coagulopathy was observed in the symptomatic, hyperinflammatory (S2) stage, indicated by a diminished platelet count and elevated levels of plasma lactate dehydrogenase. These observations corroborate previous clinical findings of appreciable coagulopathy in patients with COVID-19 early in the course of disease19,20,21. Consistent with findings from previous studies22, immune-cell profiles along the disease trajectory indicated a delayed and suppressed innate immune response to infection with SARS-CoV-2. Lymphocyte and T-cell counts were not reduced until the S2 stage, with accompanying hyperinflammation as evident from the sharp increase in C-reactive protein (Table 1). Clinical findings associated with patient blood samples supported our classification based on the disease trajectory of survivors of COVID-19 (Fig. 1a) into an immunosuppressed, presymptomatic (S1) stage followed by a symptomatic, hyperinflammatory (S2) stage to fight endogenous viral loads. The next parts of the disease trajectory comprising the resolution (S3) and convalescent (S4) stages become largely host driven (that is, SARS-CoV-2 RNA negative) and are characterized by a gradual resolution of inflammation (including the restoration of lymphocyte counts and inflammatory markers) (Table 1).

Table 1 Demographics and baseline characteristics of patients with COVID-19
Fig. 1: Investigating the EV composition along the temporal course of COVID-19.
figure 1

a, The temporal trajectory of COVID-19 is generally divided into four stages in survivors, which includes the presymptomatic (S1), hyperinflammatory (S2), resolution (S3) and convalescent (S4) phases. The first half of the trajectory (that is, SARS-CoV-2 RNA positive) is largely virus driven, characterized by immunosuppression and the response of the host to the invading pathogen, whereas the second half comprises mainly host-driven metabolism to resolve the hyperinflammation and restore physiological homeostasis. Exosome-enriched EVs were isolated from fasting plasma, and high-coverage analyses of exosomal lipidome with a focus on sterols and oxysterols were conducted. b, The particle size distribution of isolated EVs. n = 4 independent experiments. c, Particle sizes of isolated EVs analysed by nanoscale flow cytometry. Representative plot from n = 4 independent experiments.

Source data

We used differential centrifugation followed by polymer-based precipitation to isolate exosome-enriched EVs from 78 temporal plasma samples of patients with COVID-19. The size distribution of the isolated EVs was verified by nanoscale flow cytometry analysis and was in agreement with the theoretical size of EVs (<200 nm), in which >85% of EVs were particles < 180 nm and >98% of EVs were particles < 300 nm (Fig. 1b,c). Isolated EVs were enriched in exosome-specific markers23, including TSG101, CD9 and CD63, relative to the supernatant and plasma (Extended Data Fig. 1a). Marker proteins for possible plasma contaminants (annexin A1 as a microvesicle marker; ApoB and ApoE as lipoprotein markers) were not appreciably detected compared with whole plasma samples with equal amounts of total protein loaded (Extended Data Fig. 1b,c). Transmission electron microscopy images of isolated EVs showed that the vesicles exhibited a cup-shaped vesicle morphology that is characteristic of exosomes with sizes around 100 nm (Extended Data Fig. 1d,e). These observations support that the isolated EVs were exosome enriched.

Using high-performance liquid chromatography–tandem mass spectrometry (HPLC–MS/MS), we quantitatively profiled the lipidome of exosome-enriched EVs, with emphasis on various oxidized derivatives of cholesterols and phytosterols known to elicit immuno-modulatory24 and antiviral25 effects. Quantified lipid levels were normalized to the sum of total membrane lipids (that is, molar fractions to total membrane lipids, MFM). A global clustering heat map was constructed using log2-transformed fold-changes and line plots were drawn on the basis of z scores that illustrate the temporal patterns of changes in the EV lipidome during the course of COVID-19 (Fig. 2). Glycosylated sphingolipids such as globotriaosylceramides with very long-chain fatty acyls (that is, C22–C24), which have been reported to increase the susceptibility of CD4+ T cells to HIV infection26, were specifically reduced in the convalescent S4 stage compared with the presymptomatic S1 stage (Fig. 2, box (1)). Numerous plasmalogen phosphatidylethanolamines (PEps) comprising polyunsaturated fatty acyls were specifically elevated in the presymptomatic S1 stage (Fig. 2, box (2)). PEps were found to enrich in HIV particles27 and attenuate cholesterol biosynthesis by facilitating squalene monooxygenase degradation28. Impeded cholesterol biosynthesis compromises the defence of the host against SARS-CoV-2 infection13. Indeed, the Reactome pathway ‘cholesterol biosynthesis’ was significantly downregulated in HepG2 cells treated with S1 EVs relative to EVs from healthy control individuals (Supplementary Table 2; described further below). Sulfatides, which are major constituents of myelin sheath and pancreas β-cells29, were reduced in EVs from the symptomatic S2 and S3 stages, but (partially) restored in the convalescent S4 stage (Fig. 2, box (3) and Extended Data Fig. 2a). Sulfatides regulate autoimmune processes29 and modulate insulin release from pancreatic islets to prevent β-cell exhaustion30. A previous case report and in vitro studies have suggested that SARS-CoV-2 could trigger insulin-dependent diabetes and has the potential to infect pancreatic β-cells31,32. Distinct increases in EV phosphatidylcholines (PCs) (Figs. 2, box (4), 3a–c) and phosphatidylserines (Figs. 2, box (5), 3a–c) were observed in the resolution S3 and hyperinflammatory S2 stages, respectively. In a proinflammatory milieu, the phosphatidylserine content of EVs from vascular cells increases to promote calcification and coagulation33. This agrees with the reduction in platelet count (that is, coagulopathy) that was found in S2 (Table 1).

Fig. 2: Distinct temporal patterns of change in the lipidome of EVs in COVID-19.
figure 2

Hierarchical clustering was performed to aggregate lipids that exhibited similar patterns of change along the progressive course of COVID-19. n = 16, 25, 19 and 18 independent patient samples for S1, S2, S3 and S4 stages, respectively. Clustering patterns were visually examined and six clusters ((1)–(6)) of EV membrane lipids were expanded for further discussion. Left, for the global heat map, log2-transformed fold changes were plotted. Middle, for the line plots, z scores of mean of lipid levels expressed in MFM were plotted individually for lipids within each of the six clusters. Right, zoomed-in images of the six lipid clusters boxed in the global heat map. Statistical comparisons across temporal stages of COVID-19 were conducted using two-sided Kruskal–Wallis rank-sum tests and two-sided Dunn’s tests were used for pairwise comparisons. Only fold changes with P < 0.05 for each pairwise comparison are illustrated (that is, coloured blue or red); red indicates increases and blue indicates decreases in the respective phases analysed (for example, phase 2 relative to phase 1 (S2/S1)). All fold changes with P > 0.05 were set to 1 (shown in white). Gb3, globotriaosylceramide; PS, phosphatidylserine; SL, sulfatide; BMP, bis(monoacylglycero)phosphate; Cer, ceramide; LacCer, lactosylceramide.

Source data

Fig. 3: Perturbed lipid homeostasis along the disease trajectory of COVID-19.
figure 3

a,b, The temporal changes in major lipid classes (a) and individual oxysterols and sterols (b) in plasma EVs from four phases of COVID-19. Statistical significance from two-sided Kruskal–Wallis rank-sum tests is indicated with an asterisk next to the lipid name; statistical significance from pairwise comparisons is indicated using a letter-based representation in which two groups sharing a common letter are not statistically different at P < 0.05. c, The top lipids that were most significantly different in each pairwise comparison between two consecutive temporal stages defined based on both magnitudes of the P value and fold changes (FC). ac, Two-sided Dunn’s tests were used for pairwise comparisons; P < 0.05 was considered statistically significant. d, The lipid correlations between major EV lipid classes. Data were log-transformed and correlations between lipids were calculated using Spearman’s correlations analysis for each stage of COVID-19. The width of the band indicates the number of significant correlations and the colour indicates the direction of correlation. The correlation coefficient cut-off was set to ≥0.7 and P-value cut-off was set to P < 0.05. Blue shade indicates positive correlations and red shade negative correlation between two connecting lipid classes. ad, n = 16, 25, 19 and 18 independent patient samples for S1, S2, S3 and S4 stages, respectively. LPC, lysophosphatidylcholine; LPA, lysophosphatidic acid; LPI, lysophosphatidylinositol; GM3, monosialodihexosyl ganglioside; PCp, plasmalogen phosphatidylcholine; GluCer, glucosylceramide; FF-MAS; follicular fluid meiosis-activating sterol; Sph, sphingosine; Odd-PC, phosphatidylcholine with odd-chain fatty acyls; PG, phosphatidylglycerol; S1P, sphingosine-1-phosphate.

Source data

Several oxysterols and sterols were elevated in both S2 and S4 EVs (Fig. 2, box (6)), indicating dysregulated cholesterol homeostasis on a systems level. Free cholesterol, several sterols as well as oxidized cholesterol derivatives were markedly reduced in S1 compared with S2 (Fig. 3a,b and Extended Data Fig. 3a). Broadly perceived as antiviral compounds that impede viral replication34, S1-specific reductions in oxysterols corroborate a clinical status of ongoing viral replication and a blunted immune response. For example, 27-hydroxycholesterol, which was recently shown to inhibit the replication of two human coronaviruses (SARS-CoV-2 and HCoV-OC43) in vitro25, was particularly reduced in S1 EVs. Notably, 6-keto-5α-hydroxycholesterol and 3β-hydoxy-5-cholestenoic acid represent the only cholesterol derivatives analysed in this study that were significantly increased in S1 compared with S2 (Fig. 3b). 6-keto-5α-hydroxycholesterol, which is also known as 6-oxo-cholestan-3β,5α-diol, is involved in the regulation of inflammation and the immune response24,35. 6-keto-5α-hydroxycholesterol inhibits natural killer (NK) cell-mediated cytotoxicity36 and impedes cytolytic activity of CD8+ T lymphocytes37. As important components of innate immunity22, disrupting the cytolytic activities of these cells is expected to contribute to the delayed immune response and impeded virus clearance in S1. Furthermore, 3β-hydroxy-5-cholestenoic acid is a GS modulator38. Relative to non-steroidal GS modulators, GS modulators with steroid backbones possess appreciable inhibitory activity to Notch substrates in addition to inhibiting the cleavage of APP39. As 3β-hydroxy-5-cholestnoic acid was specifically reduced in the S2 phase (Fig. 3b), we postulate that GS-mediated activity may vary during the course of COVID-19. Marked increases in various oxysterols were evident in S2 EVs compared with S1, including 5α,6α-epoxy-cholesterol, 6α-hydroxy-5α-cholestanol, 7α-hydroxycholestenone, 7-keto-cholesterol, cholesta-4,6-dien-3-one, 24S,25-epoxy-cholesterol and dihydro-testis-meiosis activating sterol (dihydro-T-MAS) (Fig. 3b). General increases in oxysterols could have resulted from the overall state of oxidative stress and hyperinflammation that characterized S2. Lathosterol, a conventional indicator of cholesterol biosynthesis in the circulation40, was distinctly increased in S4 (Fig. 3b), which suggests active cholesterol biosynthesis. Sitosterol, a phytosterol that exhibits moderate effects against viral replication41 and is an indicator of intestinal cholesterol absorption40, was also elevated in S4. Thus, the convalescent S4 stage had an EV sterol profile that is indicative of active cholesterol biosynthesis and uptake.

In addition to sterols, we noted appreciable increases in triacylglycerols (TAGs) in the convalescent S4 stage (Fig. 3c). Chord diagrams illustrating the lipid correlations in each temporal stage revealed a particularly strong positive correlation between EV TAGs and major membrane lipid classes (Fig. 3d), including PC, phosphatidylethanolamine (PE), PEp, phosphatidylglycerol (PG) and phosphatidylinositol (PI), in S1. Such correlations were lost in S2, and were progressively resumed in S3 and S4. The positive coregulation between TAGs and membrane phospholipids in S4, however, was still weaker compared to S1, which was evident in the reduced thickness of connecting blue lines between these lipid classes. These correlations suggest that the mobilization of storage TAGs that drives membrane phospholipid biogenesis was intense in S1, abrogated in S2, and then partially restored in S3 and S4. The ratio of storage TAGs to their lipolytic substrates, including diacylglycerols (DAGs) and free fatty acids (FFAs) (that is, TAG/(DAG + FFA)), was significantly highest in S4 (Extended Data Fig. 3b), indicating suppressed lipolysis in the convalescent phase. In our subsequent EV treatment of A549 cells, both Gene Ontology: Biological Process (GO:BP) ‘GO:BP positive regulation of triglyceride lipase activity’ and Reactome ‘glycerophosholipid biosynthesis’ pathways were elevated in cells treated with S1 EVs compared with S4 EVs (Supplementary Table 2).

Principal component analysis of EV proteomes comprising 142 identified proteins showed that the S1 and S2 stages were well segregated (Fig. 4a,b). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses using differentially expressed proteins (DEPs) showed complement and coagulation cascades and platelet activation as the top dysregulated pathways during the course of COVID-19 (Fig. 4c). The complement system, as an integral component of the innate response that fights viruses primarily through the induction of proinflammatory responses42, was implicated in the pathogenesis of COVID-19 (ref. 43). Pharmacological agents that impede the complement cascade have been proposed as potential therapeutics for COVID-19 (refs. 44,45). We observed distinct increases in complement C1r and C1s (Fig. 4e), but not in complement C1q (Extended Data Fig. 4a,b), in S1 EVs. Complement C1q, C1r and C1s are subcomponents of the C1 complex that constitutes the classical pathway of complement activation46. Localization of classic complement proteins in plasma EVs has been reported previously47. Complement C8 γ chain, an essential component of the membrane attack complex46, was also specifically elevated in S1 EVs (Fig. 4e). By contrast, the EV proteome of S2 was marked by host-driven responses to curb viral spreading. Increase in integrin β3 (Fig. 4e), which has been reported to suppress the intracellular entry of SARS-CoV-2 by blocking viral binding to the ACE2 receptor48, was observed in S2 EVs. SARS-CoV-2 has uniquely acquired an integrin-binding Arg-Gly-Asp (RGD) motif compared with other coronaviruses49, and circulating integrins can possibly bind to and occupy the receptor-binding domains of the spike protein of SARS-CoV-2 to sterically hinder virus–receptor binding48. S2-specific increase in the human leukocyte antigen class I histocompatibility antigen B-α chain (HLA-B) was also observed (Fig. 4e). HLA-B presents viral epitopes to cytolytic CD8+ T lymphocytes and triggers the elimination of infected cells50. S2 EVs also showed marked increases in proteins implicated in hypercoagulopathy, including fibrinogen-γ chain, fibrinogen-α chain and coagulation factor V, which are part of the innate immune response (GO:0045087) and interferon-γ-mediated signalling (GO:0060333) (Fig. 4d). Induction of T-cell-mediated defence and thrombus formation can constrain the continual spreading of SARS-CoV-2. EV-localized increases in fibrinogen-γ-chain proteins, which are specifically enriched in liver hepatocytes (The Human Protein Atlas), suggest that the liver may increase EV secretion to promote coagulopathy on a systems level (described further below).

Fig. 4: Temporal changes in EV proteome along four stages of COVID-19.
figure 4

a, The segregation of the EV proteome analysed by principal component analysis. Bracketed numbers indicate the percentage of total variance explained by each component. Dim., dimension. b, A library of 142 identified proteins from plasma EVs. c, Top KEGG pathways based on pathway enrichment analysis of DEPs. d, Top terms from GO functional enrichment analyses based on DEPs. e, Changes in selected DEPs along the top dysregulated KEGG pathways. ce, DEPs were defined on the basis of statistical significance (P < 0.05) using two-sided Games–Howell’s tests. n = 3 independently pooled patient samples for each clinical phase. ER, endoplasmic reticulum; HTLV-1, human T-cell leukemia virus type 1.

Source data

The transition from the hyperinflammatory to resolution stage denotes a critical turning point of the disease course. We observed elevated proportion of raft lipids (that is, (SM + cholesterol)/PC) (Fig. 5a) and further validated the increased EV membrane anisotropy (Fig. 5b) in the hyperinflammatory compared with the resolution phase. The increases in raft lipids and EV lipid membrane anisotropy, coupled with a reduction in the steroidal GS modulator 3β-hydroxy-5-cholestenoic acid in S2, prompted us to look into the level of GS proteins localized in exosome-enriched EVs. We found that—in agreement with the endocytic localization of GS function51—the catalytic subunit of GS (that is, PS-1) was enriched in exosome-enriched EVs relative to plasma (Fig. 5c) and PS-1 was specifically elevated in S2 compared with the other temporal stages (Fig. 5c,d). Retention of cholesterol in the endocytic compartments induces PS-1 accumulation and increases the production of Aβ40 and Aβ42 (ref. 52). Accordingly, we found that relative to S3, the level of Aβ40 was significantly higher in S2 (Extended Data Fig. 5a). To determine the contribution of EVs from lung epithelial cells during COVID-19, we analysed the content of pulmonary surfactant-associated protein C (SP-C) in the isolated EVs. SP-C, which is exclusively expressed in lung tissues, is a predominantly hydrophobic protein that associates with lipids and is primarily synthesized in type II alveolar epithelial cells53. We found that mature SP-C (Fig. 5c,d) and its precursor protein (Extended Data Fig. 4c) were increased in EVs from S2, S3 and S4 stages compared with S1. EV release from the lungs therefore probably increases after the emergence of symptoms. S2 EVs also exhibited increased levels of interleukin-6 (IL-6) and cleaved caspase 3 (Fig. 5c,d and Extended Data Fig. 4d). IL-6 is a key player in the systemic hyperinflammation that is triggered by SARS-CoV-2 infection54. NOTCH-1 activation induces IL-6 production in a positive feedback loop that triggers hyperinflammatory responses16. Activated caspase 3 induces apoptosis and triggers the NLRP3 inflammasome to generate proinflammatory signals55.

Fig. 5: Lipid membrane anisotropy alters the composition and function of EVs.
figure 5

a, Changes in the proportion of raft lipids ((SM + cholesterol)/PC) across temporal stages of COVID-19. n = 16, 25, 19 and 18 independent patient samples for S1, S2, S3 and S4 stages, respectively. P values from two-sided Dunn’s tests are indicated. Statistical significance is indicated using a letter-based representation in which two groups sharing a common letter are not statistically different at P < 0.05. b, Membrane anisotropy (γ) in S2 and S3 EVs measured by incorporation of the fluorescent dye 1,6-diphenyl-1,3,5-hexatriene. n = 4 independent experiments. P value from a two-sided Wilcoxon rank–sum test is indicated. c, Immunoblot analysis of the PS-1 C-terminal fragment (PS-1(CTF))—the catalytic subunit of GS—indicates that PS-1 is enriched in exosomal fractions (EXO) relative to the supernatant (SUP) and plasma. Immunoblot analysis of PS-1(CTF), precursor pulmonary surfactant protein-C (pre-SP-C), mature SP-C, IL-6, caspase-3 and cleaved caspase 3 (activated) in EVs across four stages of COVID-19. Equal amounts of total proteins (30 µg) were loaded in each lane. d, Quantified levels of PS-1(CTF), SP-C, IL-6 and cleaved caspase 3 for the four temporal stages. P values from two-sided Dunn’s tests are shown. c,d, n = 3, 4, 3 and 3 independent patient samples for S1, S2, S3 and S4 stages, respectively. e, Confocal images of the uptake of PKH67-stained EVs (green) into cultured HepG2 cells after 4 h of incubation at 37 °C compared with the negative control in which PBS was added. Cell nuclei were stained with DAPI (blue). Representative images are from three independent experiments. Scale bars, 20 μm. f,g, Top 20 pathways from the Reactome pathway database (f) and KEGG pathway database (g) based on nominally significantly enriched gene sets obtained using gene set enrichment analysis (GSEA). n = 3 independent cultures for each treatment condition. h, Transition from S2 hyperinflammatory to S3 resolution phase is marked by a proportional reduction in EV membrane raft lipids (SM + cholesterol). The resultant decrease in membrane anisotropy alters the localization of PS-1 and other protein cargoes, providing EVs from distinct stages of infection with different biological properties. ECM, extracellular matrix; HC, healthy control individual; HDACs, histone deacetylases; eIFs, eukaryotic initiation factors.

Source data

We next treated A549 alveolar epithelial cells (as a cell type representing the primary infection site) and HepG2 cells (hepatocytes, representing a secondary, distal site of infection) with plasma EVs isolated from healthy control individuals and different temporal stages of COVID-19 (Fig. 5e–g). Hepatocytes were chosen for investigation because liver dysfunction—as indicated by abnormal levels of alanine aminotransferase and aspartate aminotransferase—was observed in 14–53% of patients with COVID-19 over the course of disease progression56. Visible uptake of EVs stained using the PKH67 dye into HepG2 cells was already observed at 4 h after incubation at 37 °C (Fig. 5e). To choose an optimized EV dosage for treatment, metabolic analysis using the Seahorse XFe96 extracellular flux analyser was conducted to measure the oxygen consumption rate of A549 and HepG2 cells that were incubated with increasing concentrations of EV proteins (Extended Data Fig. 6). The optimal dosage was determined by the minimum EV protein concentration that triggered significant reductions in the oxygen consumption rate and spare respiratory capacity in each cell type, which indicate changes in cell fitness and viability (which was 3.2 µg per 2 × 104 A549 cells and 6.4 µg per 2 × 104 HepG2 cells).

We focused on pathways that were differentially regulated (that is, from upregulated to downregulated or vice versa) by EVs isolated from successive phases of COVID-19 (Fig. 5f,g). In A549 cells, the eukaryotic translation initiation factor kinase 4 (EIF2AK4) pathway was suppressed by S2 EVs, but upregulated by S3 and S4 EVs (Fig. 5f). Response of EIF2AK4 (GCN2) to amino acid deficiency was also among the top altered pathways after SARS-CoV-2 infection of A549 cells57. Recessive mutations in EIF2AK4 are associated with pulmonary veno-occlusive disease58. Pathways governing eukaryotic translation, such as Reactome ‘eukaryotic translation initiation’ and KEGG ‘ribosome’ pathways, were suppressed by S2 EVs, but upregulated by EVs from the resolution and convalescent stages (Fig. 5f,g). SARS-CoV-2 interferes with host protein translation and protein trafficking to suppress the defences of the host against infection59, which is mediated by the viral non-structural protein 1 (NSP1)60. Indeed, the Reactome ‘translation of SARS-CoV-1 structural proteins’ pathway was upregulated in A549 cells treated with S2 EVs compared with S3 EVs (Supplementary Table 2). As we did not detect SARS-CoV-2 NSP1 proteins in both immunoblot (Supplementary Figure 1) and proteomic database searches against the deposited NSP1 sequence, we think that these effects may be mediated by EV viral RNA cargoes that were not characterized in the present study. A distinct set of differentially regulated pathways emerged in EV-treated HepG2 cells (Fig. 5f,g). The Reactome ‘pre-Notch expression and processing’ pathway was increased in HepG2 cells treated with S2 EVs. Quantitative reverse-transcription PCR (RT–qPCR) analyses of HES1 and HES5 (Extended Data Fig. 5b), the downstream transcriptional targets of Notch activation in the liver61,62, also support that Notch signalling was increased in HepG2 cells treated with EVs from the hyperinflammatory compared with the resolution phase. The Reactome ‘haemostasis’ pathway was upregulated by S2 EVs and downregulated by S3 EVs. The increased haemostasis triggered by S2 EVs in HepG2 cells agrees with our lipidomic and proteomic observations that these EVs were enriched in lipids and proteins that promote hypercoagulopathy.

Previous studies have proposed that host lipid metabolism is a primary target of positive single-strand RNA viruses such as SARS-CoV-2, which hijack the metabolism of the host to establish a favourable micromilieu for productive infection63. Using high-coverage lipidomics and sterolomics, we discovered that the temporal trajectory of COVID-19 is marked by dysregulated cholesterol metabolism. In particular, cholesterol and oxysterol were reduced in the presymptomatic S1 phase, and were markedly increased in the hyperinflammatory S2 stage together with the manifestation of clinical symptoms. Cholesterol homeostasis and the retromer complex that regulates endosomal trafficking were among the key host pathways mediating SARS-CoV-2 infection of human cells13,14. A combinatorial approach of lipidomics and proteomics for the characterization of exosome-enriched EVs uncovered several notable compositional attributes of EVs, which may either reflect disease-stage-specific metabolism of the parent cells from which they originate, or confer distinct biological functions that render concerted metabolic response on a systems scale. In this study, we found that increases in raft lipids and membrane anisotropy of EVs during the hyperinflammatory S2 stage were associated with enhanced localization and activity of GS in these EVs. GS-mediated NOTCH-1 cleavage may activate the host inflammatory response to fight SARS-CoV-2 (ref. 16). The reduction in raft lipids and PS-1 localization in EVs during the resolution phase may facilitate the cessation of Notch signalling and attenuate systemic inflammation (Fig. 5h). Increased Aβ peptide production indicates a possibly detrimental metabolic consequence of enhanced GS activity. The role of Aβ peptides in triggering cognitive impairment has been extensively reported64. Neuropsychiatric deficits such as cognitive impairment and brain fog are well-recognized aspects of post-acute COVID-19 syndrome65. Nevertheless, Aβ peptides also possess antiviral and antimicrobial properties66,67, for example, through the induction of the innate immune response mediated by the NALP3 inflammasome68. With regard to the liver, activation of Notch signalling in hepatocytes was previously shown to induce liver fibrosis in nonalcoholic fatty liver disease62. Our RNA-sequencing (RNA-seq) analyses indicated that in contrast to S3 EVs, S2 EVs impose metabolic stress on HepG2 cells. The KEGG ‘oxidative phosphorylation’ pathway was downregulated in HepG2 cells treated with S2 EVs (Fig. 5g). Compared to S2, the pathways ‘WP resistin as a regulator of inflammation’ and KEGG ‘type II diabetes mellitus’ were both downregulated when cells were treated with S3 EVs (Supplementary Table 2). Resistin is an intrahepatic cytokine expressed in response to injury and proinflammatory stimuli69, and hepatic resistin was reported to predict the severity of liver fibrosis in humans70.

In summary, we present a temporal interrogation of the compositional changes in the lipidome and proteome of exosome-enriched EVs, and provide omics-based evidence of the underlying metabolic consequences during the clinical trajectory of COVID-19. Changing EV lipid membrane fluidity can alter the localization of protein cargoes, with EVs from different stages of COVID-19 having distinct biological attributes. In particular, we found that exosome-enriched EVs from the resolution phase trigger opposing effects on cellular metabolism in vitro, and thus may be able to offset cellular responses mounted by S2 EVs. Although different proteins and RNA cargoes not characterized in this study may also contribute to the defined EV functions, we showed that enhanced PS-1 localization to exosome-enriched EVs in the hyperinflammatory S2 phase compared to the resolution S3 phase contributes to their distinct cellular responses on recipient cells. These data also serve as a useful repository of information for future research to elucidate the mechanistic clues that underlie the metabolic dysregulation and consequences in COVID-19 pathogenesis.

Methods

Study design and participants

This retrospective cohort study included two cohorts of inpatients from the Fifth Medical Center of Chinese PLA General Hospital and the Beijing Ditan Hospital Capital Medical University, which were designated hospitals for patients with COVID-19 in Beijing, China. Patients were diagnosed according to the WHO interim guidance and only survivors of COVID-19 were included in this study. We investigated four stages (S1–S4) across the clinical course of SARS-CoV-2 infection (Fig. 1a). Patient blood samples were assigned to four different stages based on clinical records of SARS-CoV-2 RNA results and laboratory findings on the date of blood collection, as well as hospital admission and discharge dates relative to the blood sampling date (Supplementary Table 1). S1 refers to presymptomatic individuals with positive results for SARS-CoV-2 RNA but who do not yet to exhibit any associated symptoms (that is, dry cough, fever and so on). S1 represents the incubation phase of COVID-19 with active SARS-CoV-2 replication and blunted innate immune response22. S2 comprises symptomatic patient samples collected within the first week of symptom onset with positive SARS-CoV-2 RNA results, a period for which the viral load of patients with COVID-19 is expected to peak followed closely by hyperinflammation71,72. S3 consists of post-hyperinflammation samples with negative SARS-CoV-2 results. S4 refers to convalescent patient samples collected within 3 d or on the day of hospital discharge, when clinical symptoms were reasonably resolved. Baseline characteristics and laboratory findings of patients with COVID-19 were summarized in Table 1. The clinical data were collected and checked independently by two physicians. This study was approved by the ethics committee of the Fifth Medical Center of the PLA General Hospital and written informed consent was obtained from all participants enrolled in this study.

Plasma collection and EV isolation

All blood samples were collected after an overnight fast in BD Vacutainer (BD 367525). Plasma was separated by centrifugation at 733g for 10 min. EVs were isolated from plasma using the Total Exosome Isolation Kit from Plasma (Thermo Fisher Scientific) according to the manufacturer’s instructions, with slight modifications. Plasma (200 µl) was first centrifuged at 2,000g for 20 min to remove cells and debris. The clear supernatant was transferred to a fresh tube, and centrifuged again at 10,000g for 20 min to further remove finer debris and large microvesicles. The resultant supernatant was then diluted with 100 µl of ice-cold 1× PBS and mixed thoroughly by vortexing. To minimize artificial oxidation, 2 µg of butylated hydroxytoluene (BHT) was immediately added to the diluted plasma, which was then incubated at 37 °C for 10 min with 10 µl of proteinase K to eliminate non-exosomal plasma proteins. At the end of the proteinase K treatment, 60 µl of Exosome Precipitation Reagent was added and the samples were incubated at 4 °C for 30 min for EV precipitation. After the incubation, the samples were centrifuged at 10, 000g for 5 min. Isolated EVs were collected as a pellet, while the supernatant comprising plasma contents without EVs were transferred to fresh tube and stored at −80 °C. Isolated EVs were resuspended in 100 µl of ice-cold 1× PBS, and immediately inactivated in 650 µl ice-cold chloroform:methanol (1:2, v/v) containing 1% BHT (w/v) for oxysterols and lipid extraction (see the next section).

Lipid extraction from EVs

Lipids were extracted following a modified Bligh and Dyer’s protocol as previously described10. EVs were incubated in extraction solvent comprising chloroform:methanol (1:2, v/v) containing 1% BHT at 4 °C for 1 h at 1,500 r.p.m. After the incubation, 350 µl ice-cold Milli-Q water and 250 µl ice-cold chloroform were added to induce phase separation. The samples were centrifuged at 16,260g at 4 °C for 5 min. The lower organic phase containing lipids and oxysterols was transferred to a new tube. The extraction was repeated once through the addition of another 450 µl ice-cold chloroform to the remaining aqueous phase, and the extractions were pooled and dried in a SpeedVac under organic mode. The lipid extract was resuspended in 200 µl chloroform:methanol (1:1, v/v). One-tenth of the lipid extract from EVs was used for targeted lipidomics analysis, and ninth-tenth of the extract was dried and used for oxysterol extraction and derivatization (see the next section). The upper aqueous–methanol phase together with the remnant EV pellets was dried in a SpeedVac under H2O mode, and the dried pellets were channelled for protein extraction and analyses.

Oxysterol extraction and derivatization

Free oxysterols were analysed in this study. In brief, the lipid extract was resuspended in 500 µl of ethanol containing 5 µg BHT. An internal standard cocktail (50 µl) comprising d7-24-hydroxycholesterol, d7-7β-hydroxycholesterol, d6-25-hydroxycholesterol, d7-27-hydroxycholesterol, d7-7α,24-dihydroxycholestenone, d7-7-keto-cholesterol, d7-7α-hydroxycholestenone, d6-lanosterol, d7-lathosterol, d6-T-MAS, d7-7-dehydrocholesterol and d6-β-sitosterol (Avanti Polar Lipids) was added to the samples. The samples were incubated at 1,200 r.p.m. for 15 min at 4 °C. At the end of the incubation, 250 µl of Milli-Q water and 1 ml of n-hexane was added. The samples were mixed thoroughly by vortexing, then centrifuged at 16,260g for 5 min 4 °C. The clear upper phase containing the oxysterols and sterols in hexane was transferred to new tube. The extraction was repeated once with another 1 ml of n-hexane. The pooled extract was dried in a SpeedVac under organic mode. Oxysterols and sterols were derivatized to obtain their picolinic acid esters before HPLC-MS/MS analysis as previously described73,74.

Lipidomics and sterolomics analysis

The EV lipidome was quantified using a high-coverage targeted lipidomic approach as described previously10, but with further refinement in the selection of internal standards used for quantification. Analyses of phospholipids and sphingolipids were conducted on a system comprising an Exion-UPLC coupled with a 6500 Plus QTRAP that runs Analyst v.1.6.3 (Sciex); neutral lipid analyses were carried out on an Agilent 1260 HPLC (Agilent Technologies) connected to a 5500 QTRAP that runs Analyst v.1.6.3 (Sciex). The internal standard cocktail was adjusted according to the endogenous abundances of various lipid classes in isolated EVs, which comprised d9-PC32:0(16:0/16:0), d9-PC36:1p(18:0p/18:1), d7-PE33:1(15:0/18:1), d9-PE36:1p(18:0p/18:1), d31-phosphatidylserine (16:0/18:1), d7-PG33:1(15:0/18:1), d7-PI33:1(15:0/18:1), d7PA33:1(15:0/18:1), C14-BMP, d8-SM d18:1/18:1, Cer d18:1/d7-15:0, GluCer d18:1/8:0, d3-LacCer d18:1/16:0, GM3 d18:1/18:0-d3, Gb3-d18:1/17:0, sulfatide-d18:1/17:0, d7-LPC 18:1, d7-LPE 18:1, LPA-C17:0, LPI-C17:1, LPS-C17:1, LPG-C17:1, DAG(16:0/16:0)-d5, DAG(18:1/18:1)-d5, d7-cholesteryl sulfate, S1P-d17:1 and Sph-d17:1 obtained from Avanti Polar Lipids, TAG(14:0)3-d5, TAG(16:0)3-d5, TAG(18:0)3-d5, d6-CE18:0 and d6-cholesterol purchased from CDN Isotopes, and d3-16:0-carnitine from Cambridge Isotope Laboratories. d31-FFA-16:0 from Sigma-Aldrich and d8-FFA-20:4 from Cayman Chemicals were used for quantification of saturated/monounsaturated fatty acids and polyunsaturated fatty acids, respectively. Oxysterols and sterols were separated on a 1.9 µm Hypersil Gold C18-selectivity column (50×2.1 mm) (Thermo Fisher Scientific) and analysed on an U3000-DGLC coupled with a QTRAP 6500 Plus (Sciex) in their derivatized forms, as described previously74. A gradient system comprising acetonitrile:methanol:water (45:45:40, v/v/v) containing 0.1% acetic acid as mobile phase A and acetonitrile:methanol:water (45:45:10, v/v/v) containing 0.1% acetic acid as mobile phase B was used at a flow rate of 550 µl min1 and column oven temperature at 35 °C. The gradient program adopted was as follows: 0–2 min, 5% B; 2–3 min, 5–25% B; 3–6 min, 25–65% B; 6–8 min, 65–100% B; 8–12 min, 100% B; 12–13 min, 100–5% B; 13–15 min: 5% B. The identities of individual oxysterols were validated based on both the multiple reaction monitoring transitions and retention time relative to their authentic reference compounds. Quality-control samples were inserted every ten samples throughout the batch analysis, and the performance of quality-control samples is illustrated in Supplementary Figs. 2, 3. PC, phosphatidylcholine; PCp, plasmalogen phosphatidylcholine; PE, phosphatidylethanolamine; PEp, plasmalogen phosphatidylethanolamine; PG, phosphatidylglycerol; PI, phosphatidylinositol; PA, phosphatidic acid; BMP, bis(monoacylglycero)phosphate; SM, sphingomyelin; Cer, ceramide; GluCer, glucosylceramide; LacCer, lactosylceramide; GM3, monosialodihexosyl ganglioside; Gb3, globotriaosylceramide; LPC;lysophosphatidylcholine; LPE, lysophosphatidylethanolamine; LPA, lysophosphatidic acid; LPI, lysophosphatidylinositol; LPS, lysophosphatidylserine; LPG, lysophosphatidylglycerol; DAG, diacylglycerol; S1P, sphingosine-1-phosphate; Sph, sphingosine; TAG, triacylglycerol; CE, cholesteryl ester; FFA, free fatty acid.

Proteomic analysis

Approximately 150 µg of EV protein was obtained from each sample after lipid and oxysterol extraction. Proteomic analyses were conducted using three pooled samples from each temporal stage. Each sample was pooled from two randomly selected samples in the same clinical group. Each pooled sample (50 µg protein) was reconstituted using 8 M urea in 50 mM NH4HCO3. The proteins were then digested with trypsin as previously described75. In brief, the proteins were reduced with 10 mM DTT at 37 °C for 1 h and then alkylated with 25 mM iodoacetamide at room temperature for 1 h in the dark. After reduction and alkylation, the lysates were diluted with 6× the volume of 50 mM NH4HCO3, and digested with trypsin at 37°C for 18 h using an enzyme:substrate ratio of 1:50. All tryptic peptides were desalted by StageTip and dried with a SpeedVac before use. The resuspended peptides were analysed on a LTQ Orbitrap Elite mass spectrometer running Xcalibur v.4.1 (Thermo Fisher Scientific) coupled online to an Easy-nLC 1000 (Thermo Fisher Scientific) in a data-dependent mode. The peptides were separated by reverse-phase LC with a 150 μm (inner diameter) × 250 mm (length) analytical column packed with C18 particles of 1.9 µm diameter. The mobile phases comprise buffer A (100% H2O, 0.1% formic acid) and buffer B (100% acetonitrile, 0.1% formic acid), and a 180-min non-linear gradient was used for chromatographic separation. All MS measurements were performed in positive ion mode. Precursor ions were measured in the Orbitrap analyser at 240,000 resolution and a target value of 106 ions was used for each MS scan. The twenty most intense ions from each MS scan were isolated, fragmented and measured in the linear ion trap. The collision-induced dissociation normalized collision energy was set to 35. The data were analysed using a pre-release version of Thermo Scientific Proteome Discoverer software v.2.3. The Homo sapiens proteome sequence database containing 20,417 entries downloaded from UniProt was used for the database searching. Normalization was performed against the total peptide intensity of each sample.

Nanoscale flow cytometry analysis

EVs were analysed using an A50 Micro Nanoscale Flow Cytometer (Apogee FlowSystems) equipped with a blue (488 nm) laser76. Before analysis, all samples were diluted tenfold to avoid a swarm effect. Parameters in the control panel were set as follows, sheath pressure = 150 mbar, number of flush cycles = 2. Calibration of flow cytometer performance is conducted using a standard reference bead mix (ApogeeMix, Apogee Flow Systems) composed of a mixture of silica nanoparticles with diameters of 180 nm, 240 nm, 300 nm, 590 nm, 880 nm and 1,300 nm, and green fluorescent polystyrene nanoparticles of diameters 110 nm and 500 nm. The threshold of medium angle light scatter (ALS) is set at 22 arbitrary units to eliminate optical and electronic background noise without losing particles of interest. Photomultiplier tube voltages were set to 324 V for small ALS, 350 V for medium ALS and 291 V for large ALS for size analysis. All measurements were performed in log mode. Noise levels in the photomultiplier tube panel were kept below 0.5. All data were analysed using the Apogee Histogram Software.

Transmission electron microscopy

Transmission electron microscopy images of EVs were taken on the basis of a previous protocol77. In brief, isolated EVs were reconstituted in 200 µl of 1× PBS and diluted 10× and 500×. EV samples at two different dilutions were mixed with an equal volume of 4% paraformaldehyde. Drops of fixed EVs (5 µl) were deposited on clean parafilm, and a formvar-carbon-coated grid was floated on each drop of fixed EVs with the coated side facing down. The samples were left to dry at room temperature for 20 min and then washed with PBS. Next, samples were fixed in 1% glutaraldehyde for 5 min, washed with deionized water eight times, and stained with 3% uranyl acetate for 5 min. Excess liquid was removed and the stained EVs were observed using a transmission electron microscope (JEM 1400) at 80 kV.

EV membrane fluidity assay

EV membrane fluidity was measured using the fluorescent dye 1,6-diphenyl-1,3,5-hexatriene (Sigma-Aldrich, D208000) as previously described78. In brief, 10-nmol-equivalent PCs from EVs, quantified on the LC–MS using d9-PC32:0(16:0/16:0) as an internal standard, were resuspended in 1 ml of HEPES buffer containing protease inhibitors at pH 7. The optical density of the mixture was checked to be <0.08 at 348 nm to ensure subsequent fluorescent measurements were not optically obscured. Labelling of EV membranes with 1,6-diphenyl-1,3,5-hexatriene probes was carried out at a molecular ratio of probe:PC of 1:360 for 30 min at room temperature with gentle agitation in the dark. The labelled samples were excited at 35 °C (to mimic body temperature) with vertically polarized light in a Hitachi F-7000 fluorescence spectrophotometer. Emission intensities were measured parallel (IVV) and perpendicular (IVH) to the vertical plane of polarization. To eliminate the effect of instrument polarization response, a grating correction factor (G) was taken into consideration, calculated by exciting the samples with horizontally polarized light and measuring the parallel and perpendicular components, where G = IHV/IHH. Anisotropy (γ) was calculated according to the equation (1),

$$\gamma = ({\mathrm{IVV}} - G{\mathrm{IVH}})/({\mathrm{IVV}} + 2G{\mathrm{IVH}})$$
(1)

Four replicates of EV samples isolated from pooled plasma of S2 and S3 states were measured.

Cell culture

HepG2 and A549 cells were obtained from the American Type Culture Collection. HepG2 cells were maintained in DMEM (high-glucose) medium (Gibco) supplemented with 10% FBS, and A549 were maintained in RPMI 1640 medium (Gibco) supplemented with 10% FBS. Cells were cultured at 37 °C in an atmosphere of 5% CO2.

Confocal imaging of EV uptake by HepG2 cells

EVs isolated from the plasma of patients with COVID-19 were stained with PHK67 dye (Sigma-Aldrich, MINI67) according to the manufacturer’s instructions. In brief, 100 µl of resuspended PKH67 was added to either 100 µl EV suspension or PBS (negative control). EVs were evenly mixed with the dye solution and left to stand at room temperature for 5 min. Next, 200 µl of 1% bovine serum albumin was added to bind excess dye, and the samples were transferred to 100-kDa Vivaspin filters and centrifuged at 16,260g for 3 min. The upper phase containing the stained EVs was washed twice with 1× PBS, and incubated with HepG2 cells for 4 h at 37 °C. At the end of the incubation, the cells were stained with DAPI and observed using a confocal laser scanning microscope (Observer Z1, Carl Zeiss).

Treatment of HepG2 and A549 cells with EVs

EVs were isolated from 1 ml of pooled plasma from each of the five clinical stages, which included healthy control individuals, as well as individuals from the S1, S2, S3 and S4 phases of COVID-19. Cells (106 per well) were incubated with culture medium containing the appropriate concentrations (measured in terms of quantity of EV proteins) of EVs (that is, 160 µg per well for A549 cells; 320 µg per well for HepG2 cells) in 12-well plates for 12 h at 37 °C. At the end of the incubation, the cell culture medium was aspirated. The cells were washed three times with ice-cold 1× PBS. Cells designated for transcriptome analyses were immediately inactivated through the addition of TRIzol reagent (Invitrogen, 15596-026). Three separate wells for each clinical group were collected for transcriptome analyses.

RNA analysis using RT–qPCR

Total RNA was extracted from cells inactivated with 1 ml of TRIzol reagent (Invitrogen, 15596-026). For each phase separation, 0.2 ml of chloroform was added and the aqueous phase was collected in fresh tubes and mixed with 0.5 ml of isopropyl alcohol to pellet the RNA. RNA pellets were washed with 75% ethanol and digested with DNase for 30 min to remove DNA contamination, then pelleted through the addition of NaOAc and LiCl. RNA pellets were then washed with 75% ethanol and reconstituted in nuclease-free water. cDNA was synthesized using the iScript cDNA Synthesis kit (Bio-Rad, 1708891). RT–qPCR was performed using the SYBR Green PCR kit (Bio-Rad, 1725120) on a Bio-Rad CFX Connect 96-Real-Time PCR Detection machine. Relative expression values of HES1 and HES5 were normalized to those of GAPDH. Gene-specific primers used for amplification are listed in Supplementary Table 3.

PromethION library preparation and sequencing

First-strand cDNA was synthesized in a 20 μl reaction mixture with 13 μl of purified nucleic acids from each sample and 100 pmol of oligo(dT)15 primers. For the double-stranded cDNA synthesis, 100 pmol of random primers (Promega, C118A) was added. PromethION library preparation was performed according to the manufacturer’s instructions for barcoding cDNA (SQK-LSK109, EXP-NBD104 and EXP-NBD114). When multiplexing, all of the samples were pooled together. ONT MinKNOW software (v.19.10.1) was used to collect raw sequencing data, and Guppy (v.3.2.4) was used for real-time base calling. The PromethION was run for up to 48 h. Adaptor and low-quality sequences of raw data were trimmed using Qcat with the default parameters. Trimmed reads were aligned to the human genome (version GRCH38) using mimimap2 (ref. 79). The number of reads aligned to human genes was calculated using BEDtools multicov with the human gene annotation file GFF (GENCODE v.37).

Western blot analysis

Proteins were extracted using RIPA lysis buffer with protease inhibitor cocktail (Sigma-Aldrich). The protein content was determined using a Pierce BCA protein assay kit (Thermo Fisher Scientific) according to the manufacturer’s instructions. After boiling, lysates containing 30 µg protein were subjected to SDS–PAGE and blotted onto nitrocellulose membranes (Pall). Primary antibodies included: anti-CD9 (Santa Cruz, sc-13118, 1:500), anti-CD63 (Santa Cruz, sc-5275, 1:500), anti-TSG101 (Abcam, ab125011, 1:2,500), anti-PS-1 (Abcam, ab76083, 1:2,000), anti-caspase-3 (CST, 9662s, 1:2,000), anti-SP-C (Santa-Cruz, sc-13979, 1:500), anti-IL-6 (Abcam, ab233706, 1:2,000), anti-annexin A1 (CST, 32934, 1:1,000), anti-ApoB (Merck Millipore, AB742, 1:3000) and anti-ApoE (Merck Millipore, AB947, 1:1,000) antibodies. Secondary antibodies included: goat anti-rabbit HRP (ZSGB-BIO, ZB-2305, 1:5,000), goat anti-mouse HRP (ZSGB-BIO, ZB-2301, 1:5,000) and rabbit anti-goat HRP (Biosharp, BL004A, 1:2,500) antibodies.

Statistical analysis

To examine the changes in EV membrane lipid composition, for all statistical analyses, lipids and sterols were expressed in molar fractions normalized to total membrane lipids (that is, MFM), which comprises the summed levels of all phospholipids, sphingolipids and free cholesterol quantified in µmol. Statistical comparisons between different temporal stages of COVID-19 for lipids were conducted using Kruskal–Wallis rank-sum tests and Dunn’s tests for pairwise comparisons. P < 0.05 was considered to be statistically significant. Statistical significance of all pairwise comparisons was indicated using a letter-based representation, in which two groups that share a common letter are not statistically different at the P < 0.05 level. Fold change was calculated as the ratio of means of two groups. Hierarchical clustering using the Ward method (ward.D2 in R hclust) was performed on log2-transformed fold-change values, which were set to 1 (no change) if corresponding to P > 0.05. Heat maps were drawn using ComplexHeatmap v.2.7.10.9001. Clustering was visually evaluated and patterns of interests were manually selected. For lipid species within each cluster, standardized group means (mean = 0 and s.d. = 1) of each group are presented as line plots. Standardized group means of lipid class sums and sterols are presented as a heat map, where P < 0.05 in a Kruskal–Wallis test is indicated with an asterisk after the names, and pairwise comparison results from Dunn’s tests are indicated in each cell using the letter-based representation. Fold change and −log10[P] values are presented as volcano plots. Spearman correlation coefficients were calculated between lipid species within samples in each phase using Hmisc v.4.4-2. Correlations with corresponding P < 0.05 were visualized using chord diagrams drawn with chorddiag v.0.1.2 and circlize v.0.4.12, in which the band width indicates the number of correlations and the colour indicates the direction of correlation. The correlation coefficient cut-off was set to ≥0.7 and the P-value cut-off was set to P < 0.05. The blue shade indicates positive correlations and the red shade indicates negative correlation between lipids. Statistical analyses were performed using R v.4.0.2.

For analysis of proteomic data, proteins with missing intensity values in more than one sample were removed and not analysed. Remaining missing values were replaced by zeroes. The intensities of individual proteins were normalized to the sum of total intensities in each sample before statistical analyses. Statistical significance for pairwise comparisons between temporal stages was assessed using Games–Howell’s test. Pathway analyses were conducted using DEPs. DEPs were assigned to GO local database annotations in Perl scripts using Fisher’s exact tests, for which proteins were categorized into three groups, namely biological process, cellular component and molecular function. Pathway enrichment analysis was carried out using the KEGG database (http://www.genome.jp/kegg/pathway.html) using Fisher’s exact tests in Perl scripts.

For analysis of RNA-seq data, differentially expressed genes between the various clinical groups were obtained using the R package DESeq2 (see source data associated with Fig. 5). The Wald statistics from each of the comparisons were then used for pre-ranked GSEA analysis (10,000 permutations) using the fgsea package in R v.3.6.2, ClusterProfiler v.3.14.3 and Bioconductor v.3.1. Gene symbols were mapped to Entrez IDs using org.hs.eg.db v.3.10.0. Gene sets originating from KEGG, Reactome, GO and WikiPathways were obtained from MsigDB v.7.3 and used in the analysis. Summary plots of the nominally significantly enriched gene sets were plotted using ggplot2.

Reporting Summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.