Changes to the sebum lipidome upon COVID-19 infection observed via non-invasive and rapid sampling from the skin

Population metadata overview

The study population analysed in this work included 67 participants, comprising 30 participants presenting with COVID-19 clinical symptoms (and an associated positive COVID-19 RT-PCR test) and 37 participants presenting without. A summary of the metadata is shown in Table 1.

There were more male participants in the COVID-19 positive group (M:F ratio of 0.57) compared to the participant population overall (M:F ratio of 0.52); given recruitment took place in a hospital environment, this may reflect increased severity amongst males. ²⁴ Age distributions for COVID-19 positive and negative cohorts were almost identical (mean age of 64.7 years and 65.0 years respectively). Comorbidities are associated with both hospitalisation and more severe outcomes for COVID-19 infection, but will also alter the metabolome of participants, representing both a causative and confounding factor. The impact on classification accuracy of these comorbidities was tested by stratifying participant data by comorbidity to see if separation improved; this process is described in the following sections. In this pilot study, comorbidities were less well represented in the cohort of COVID-19 positive participants than in the cohort of COVID-19 negative participants.

In terms of diagnostic indicators, levels of C-Reactive Protein (CRP) were significantly higher for COVID-19 participants, whilst lymphocyte and eosinophils levels were lower. A two-tailed Mann Whitney U test on the CRP indicator provided a p-value of .031, and on the lymphocytes a p-value of .004. Effect sizes (calculated by Cohen’s D) were 0.56 and 0.85 respectively. COVID-19 positive participants were also markedly more likely to present with bilateral chest X-ray changes (21 out of 30 COVID-19 positive patients, versus just 2 out of 37 COVID-19 negative patients). With regard to outcomes, COVID-19 positive participants experienced higher rates of requiring oxygen CPAP, higher rates of escalation, and lower survival rates. These observations were in agreement with literature descriptions of COVID-19 symptoms and progression. ¹¹

Overview of features identified by Liquid Chromatography Mass Spectrometry (LC-MS)

998 features were identified reproducibly (present in greater than 90% of pooled QC LC-MS injections, with coefficient of variation below 20% across all pooled QCs and signal to noise ratio greater than three) by LC-MS across a series of pooled QCs and these formed the basis of the analysis in this work. Differences between COVID-19 positive and negative participants were observed across a number of lipids and metabolites, with the most consistent difference seen in reduced lipid levels, especially triglycerides (Figure 1).

• Download figure
• Open in new tab

Figure 1: Volcano plot of features for COVID-19 positive (n=30) versus negative (n=37), labelled features validated by MS/MS, points scaled to significance

Aggregate levels of triglycerides identified by MS/MS were depressed for COVID-19 positive participants, and also for ceramides, albeit fewer lipids of the latter class were identified and validated. The distributions of the natural log of aggregated lipid ion counts by class were not characterised as normal by Shapiro-Wilk normality tests. ²⁵ Two-tailed Mann-Whitney U-tests were performed to test the significance of aggregate levels of these lipid classes. These resulted in p-values of 0.022 and 0.015 for triglycerides and ceramides respectively, with effect sizes (calculated by Cohen’s D) of 0.44 and 0.57, indicative of medium effect size. These results are suggestive of dyslipidemia within the stratum corneum due to COVID-19. The alteration in levels of triglycerides between positive and negative cohorts is comparable to that for CRP or for lymphocytes as indicators of COVID-19 status (Figure 2).

• Download figure
• Open in new tab

Figure 2: Boxplots of diagnostic indicators versus triglyceride levels

Other work has found evidence of dyslipidemia in plasma from COVID-19 positive patients, ^16,15,18 although evidence of whether upregulation or downregulation is dominant for these lipid classes is mixed. Plasma triglyceride (TAG) levels have been found to be elevated in blood plasma for mild cases of COVID-19, but TAG levels in plasma may also decline as the severity of COVID-19 increased. ²⁶

It should be remembered, however, that the primary role of skin is barrier function, and lipid expression in the stratum corneum depends on de novo lipogenesis – in fact nonskin sources such as plasma provide only a minor contribution to sebum lipids. ²⁷ To the extent that the virus sequesters lipids for its own reproduction, it is possible that this causes deficiency in the expression of sebum lipids.

Population-level clustering analyses

No clustering was identifiable at the total population level by principal component analysis (PCA), i.e. by unsupervised analysis. Orthogonal partial least squares discriminant analysis (OPLS-DA) performed on the same data set still revealed limited separation (Figure 3). A confusion matrix was constructed using a pairwise knock-out approach to establish training models; projecting these models onto the excluded participants to test accuracy (Table 2) showed sensitivity of just 63% and specificity of 70%. Given the wide range of comorbidities and the lack of age-matching, this is not unexpected.

• Download figure
• Open in new tab

Figure 3: OPLS-DA of 67 participants, classified by COVID-19 positive / negative

Investigation of confounding factors

To test the impact of age and diagnostic indicators (CRP, lymphocytes and eosinophils), these variables were pareto-scaled and included in the matrix for OPLS-DA modelling. Variable importance in projection (VIP) scores for lymphocytes, CRP, and eosinophils were 2.02, 1.34 and 1.17 respectively, ranking 16, 243 and 293 out of 1,002 total features. As a single feature, depressed lymphocyte levels show high correlation with COVID-19 positive status, consistent with lymphocyte count being both a diagnostic and prognostic biomarker. ²⁸ Age as a vector had a VIP score of just 0.19 (ranking 849 out of 1,002 total features), indicating that age is a smaller influencer of stratum corneum lipids than other factors. Overall, OPLS-DA separation did not improve by the addition of these age and diagnostic indicators.

To test whether separation would improve in smaller / more homogenous groups, separate OPLS-DA models were built for each split of the population by comorbidity. If model performance improved (measured by goodness of fit - R2Y - and sensitivity and specificity) then this could indicate that sebum lipid profiling would perform better if models were constructed based on stratified and matched datasets. Table 4 shows the results for these metrics across the different modelled subsets.

Separation generally improved as the data were grouped more finely, but for most subpopulations there was no improvement in the modelled predictive power. Four subsets did however show more interesting improvements in model performance. These were the subsets with a specific comorbidity that were being treated by medication (high cholesterol, T2DM and IHD) and the subset undergoing treatment with statins. OPLS-DA modelling of the subset of participants under medication for high cholesterol (Figure 4) showed both good separation (R2Y of 1.00) and better sensitivity and specificity (Table 3). This subgroup was treated with lipid-lowering agents, specifically statins.

• Download figure
• Open in new tab

Figure 4: OPLS-DA of 15 participants treated for high cholesterol, COVID-19 / positive negative

Stratification by other comorbidities also yielded improved results, as did stratification by medication. Summary data for these models is set out below in Table 4. The best sensitivity and specificity results were found in the OPLS-DA model for high cholesterol as shown above, but stratification by ischaemic heart disease, type 2 diabetes mellitus and by medication with statins all yielded an improvement versus the base model.

OPLS-DA modelling of the subset of participants under medication for type-2 diabetes mellitus (T2DM) showed good separation (Figure S1, Supplementary Information), with sensitivity of 78% and specificity of 75%. This subgroup was typically being treated with oral hypoglycaemics, for example metformin, in some cases with insulin and in some instances with diet control only. The subgroup comprising participants undergoing treatment for ischemic heart disease (IHD) also showed much better separation (R2Y of 1.00, Figure S2, Supporting Material), again with better sensitivity and specificity of 75% and 86% respectively. This subgroup received varied medication, but participants presenting with IHD were also being prescribed statins.

As shown in Table 4, model performance also improved for a stratified dataset based on those participants taking statins. Given that statins control cholesterol and lipid levels, this may have provided a more similar “baseline” against which to measure perturbance in the lipidome by COVID-19. Analysing all patients taking statins (which includes both participants treated for high cholesterol and also participants with poor diabetic control or history of ischaemic heart disease, where statins are routinely added prophylactically to improve long-term outcomes) showed improved separation by OPLS-DA modelling (Figure S4, Supplementary Information), with R2Y of 0.74, sensitivity of 71% and specificity of 76%.

Looking across the models, there was limited commonality in the features identified as significant in differentiating between COVID-19 positive and negative. Whilst some features had high VIP scores in all subgroups (Figure 5), many did not, a possible indicator of overfitting due to the smaller groups when stratified. Where overlap does occur between the features, this may reflect the natural overlap between the subset populations, for example the subsets of participants presenting with ischaemic heart disease and with high cholesterol are largely subsets of the participants receiving treatment by statins.

• Download figure
• Open in new tab

Figure 5: Heat map of VIP scores ranked by commonality to different subgroup OPLS-DA models

Definitive separation has not proved possible in this pilot study, given that too few datapoints were available to rigorously stratify by medication or by comorbidity. Nonetheless, at the aggregate level, participants with a positive clinical COVID-19 diagnosis present with depressed lipid levels (triglycerides and ceramides in particular), with the possibility of consequently reduced barrier function. Furthermore, these findings suggest that better matching of participants could yield a clearer separation of positive and negative COVID-19 participants by their lipidomic profile. Overfitting is a risk in any pilot study with small n; this risk can only be reduced through both a larger training set of data and subsequently testing the models on future validation sets, made possible through cohesive efforts such as the work of the MS Coalition.

Another point to note is a possible lack of confounders in the participant population from seasonal respiratory viruses. Whilst the COVID-negative patients included patients with respiratory illnesses (e.g. COPD, asthma) and COVID-like symptoms, samples were collected between May and July, when the incidence of respiratory viruses is generally low. Both the common cold and influenza have some symptoms overlap with COVID-19 and may possibly lead to alterations to lipid metabolism that could interfere with the identification of features related to COVID-19 infection. Such viruses within the UK are more prevalent in autumn and winter. ²⁹ Whilst it seems unlikely that seasonal respiratory viruses were a major confounding factor in this work, this is a factor that will need to be taken into account in future studies, and may also allow the opportunity to test sebum’s selectivity and specificity with regard to other respiratory viruses.

A final limitation of this study is inconsistency in the timeline between onset of symptoms, hospital admission, PCR test and sebum sampling, which was an inevitable consequence of collecting samples in a pandemic situation. Patients were sampled immediately upon recruitment to the study. This means that the range in time between symptom onset and sebum sampling ranged from 1 day to > 1 month. Future work should explore longitudinal sampling of patients, to establish how quickly (or whether at all) the sebum lipidome returns to normal, and whether it has prognostic power. This will help to inform the practical utility of sebum in clinical or mass testing.

Participant recruitment

The participants included in this study were recruited at NHS Frimley Park NHS Trust, totalling 67 participants. Ethical approval for this project (IRAS project ID 155921) was obtained via the NHS Health Research Authority (REC reference: 14/LO/1221).

Materials and chemicals

The materials and solvents utilised in this study were as follows: gauze swabs (Reliance Medical, UK), 30 mL Sterilin™tubes (Thermo Scientific, UK), 10 mL syringes (Becton Dickinson, Spain), 2 mL microcentrifuge tubes (Eppendorf, UK), 0.2 µm syringe filters (Corning Incorporated, USA), 200 µL micropipette tips (Starlab, UK) and Qsert™clear glass insert LC vials (Supelco, UK). Optima™(LC-MS) grade methanol was used as an extraction solvent, and Optima™(LC-MS) grade methanol, ethanol, acetonitrile and 2-propanol were used to prepare injection solvents and mobile phases. Formic acid was added to the mobile phase solvents at 0.1% (v/v). Solvents were purchased from Fisher Scientific, UK.

Sample collection, inactivation and extraction

Collection of the samples was performed by researchers from the University of Surrey at Frimley Park NHS Foundation Trust hospitals. Participants were identified by clinical staff and were categorised by the hospital as either “query COVID” (meaning there was clinical suspicion of COVID-19 infection) or “COVID positive” (meaning that a positive COVID test result had been recorded during their admission). Each participant was swabbed on the right side of the upper back, using 15 cm by 7.5 cm gauzes that had each been folded twice to create a four-ply swab. The surface area of sampling was approximately 5 cm x 5 cm, pressure was applied uniformly whilst moving the swab across the upper back for ten seconds. The gauzes were placed into Sterilin polystyrene 30 mL universal containers. Samples were transferred from the hospital to the University of Surrey by courier within 4 hours of collection, whereupon the samples were then quarantined at room temperature for seven days. Finally, the vials were transferred to minus 80°C storage until required.

Alongside sebum collection, metadata for all participants was also collected covering inter alia gender, age, comorbidities (based on whether the participant was receiving treatment), the results and dates of COVID PCR (polymerase chain reaction) tests, bilateral chest X-Ray changes, smoking status, and whether the participant presented with clinical symptoms of COVID. Values for lymphocytes, CRP and eosinophils were also taken here the most extreme values during the hospital admission period were recorded. These were not collected concomitantly with the sebum samples.

The analysis of the obtained samples was adapted from Sinclair et al. ³⁰ To extract analytes from the sample gauzes, the Sterilin vials and contents were allowed to equilibrate to room temperature after which 9 mL methanol was added, followed by vortex-mixing for 10 sec. The solution with gauze was then sonicated for 30 min at ambient temperature. The metabolite-rich methanol was then filtered through a 0.2 μm filter to yield three equal aliquots of 2 mL fractions in 2 mL Eppendorf tubes, and a 0.2 mL aliquot reserved to create a pooled QC in a separate 10 mL scintillation vial. Each 2 mL sample was then dried under nitrogen for 3 hours, leaving a lipid pellet, and frozen at minus 80 °C until the day of analysis.

To reconstitute the samples, the Eppendorf tubes and contents were allowed to equilibrate to room temperature and the dried pellet was dissolved in 200 μL of methanol:ethanol (v/v, 50:50). The reconstituted solution was vortex mixed for 20 seconds followed by sonication (5 min) and centrifugation at 12 000 g for 10 minutes. 150 μL of supernatant was extracted and transferred to a LC-MS insert vial for analysis. Samples were analysed over a period of five days. Each day consisted of a run incorporating solvent blank injections (n=5), pooled QC injections (n=3), followed by 16 participant samples (triplicate injections of each) with a single pooled QC injection every six injections. Each day’s run was completed with pooled QC injections (n=2) and solvent blanks (n=3). A triplicate injection of a field blank was also obtained.

Instrumentation and software

Analysis of samples was carried out using a Dionex Ultimate 3000 HPLC module equipped with a binary solvent manager, column compartment and autosampler, coupled to a Orbitrap Q-Exactive Plus mass spectrometer (Thermo Fisher Scientific, UK) at the University of Surrey’s Ion Beam Centre. Chromatographic separation was performed on a Waters ACQUITY UPLC BEH C18 column (1.7 µm, 2.1 mm x 100 mm) operated at 55 ºC with a flow rate of 0.3 ml min^-1.

The mobile phases were as follows: mobile phase A was acetonitrile:water (v/v 60:40) with 0.1% formic acid, whilst mobile phase B was 2-propanol:acetonitrile (v/v, 90:10) with 0.1% formic acid (v/v). An injection volume of 5 µL was used. The initial solvent mixture was 40% B, increasing to 50% B over 1 minute, then to 69% B at 3.6 minutes, with a final ramp to 88% B at 12 minutes. The gradient was reduced back to 40% B and held for 2 minutes to allow for column equilibration. Analysis on the Q-Exactive Plus mass spectrometer was performed in split-scan mode with an overall scan range of 150 m/z to 2 000 m/z, and 5 ppm mass accuracy. Split scan was chosen to extend the m/z range from 150 to 2 000 m/z whilst maximising the number of features identified. ^31,32 MS/MS validation of features was carried out on Pooled QC samples using data dependent acquisition mode. Operating conditions are summarised in Table 5.

Data processing

LC-MS outputs (.raw files) were pre-processed for alignment, normalisation and peak identification using Progenesis QI (Non-Linear Dynamics, Waters, Wilmslow, UK), a platform-independent small molecule discovery analysis software for LC-MS data. Peak picking (mass tolerance ±5 ppm), alignment (RT window ±15 s) and area normalisation was carried out with reference to the pooled QC samples. Features identified in MS were initially annotated using accurate mass match with Lipid Blast in Progenesis QI, whilst validation was performed using data dependent MS/MS analysis using LipidSearch (Thermo Fisher Scientific, UK) and Compound Discoverer (Thermo Fisher Scientific, UK). This process yielded an initial peak table with 14,160 features. All those features with a coefficient of variation across all pooled QCs above 20% were removed, as were those that were not present in at least 90% of pooled QC injections. These features were then field blank adjusted: all those features with a signal to noise ratio below 3x were also rejected. The remaining set of 998 features were deemed to be robust, reproducible and suitably distinct from those found in the field blank.

Inclusion criteria were also applied to participant data, requiring both full completion of metadata and also agreement between the result of the PCR COVID-19 test (Y/N) and the clinical diagnosis for COVID-19 (Y/N). Whilst these inclusion criteria reduced the total number of participants from n=87 to n=67, this was considered worthwhile given the potential for misdiagnosis to confound the development of statistical models.

Statistical Analysis

Data processing and analysis of the peak:area matrix was conducted through user-written scripts in the statistical programming language R, using the RStudio graphical user interface package. Both PCA and OPLS-DA were performed for classification and prediction of data. Leave-pair-out cross-validation was used for OPLS-DA model validation. The data were pareto-scaled in RStudio as part of all statistical analyses, without replacement of missing values.