A step-by-step workflow in this study is presented in Fig. 1.
Datasets used in this study
We collected the data of HDL-c, LDL-c, TG and TC from a published genome-wide association study (GWAS)  involving in a total of 312571 individuals (the proportion of European-ancestry: 73.74%). The GWAS summary data of dyslipidemia, defined as the disorders of lipoid metabolism (ICD-9 code: 272), was obtained from the Genetic Epidemiology Research on Adult Health and Aging (GERA) (https://cnsgenomics.com/content/data) with 53991 European individuals .
For assessment of associations with COVID-19 risk, we used the latest COVID-19 GWAS results from GRASP database (https://grasp.nhlbi.nih.gov/Covid19GWASResults.aspx). The phenotype used in this GWAS was case/control for COVID-19 infection containing 1221 positive COVID-19 tests and 4117 negative tests from UK Biobank individuals (released on June 5, 2020). To explore the causal effect of blood lipids on severity of COVID-19, we also accessed GWAS of severe COVID-19 defined as respiratory failure  (http://www.c19-genetics.eu/). This GWAS included 1610 severe COVID-19 patients and 2205 control participants from Italy and Spain. These two GWASs were all conducted with the correction of age, sex and top 10 principal components.
SNPs Filter and Data Standardization
For each exposure GWAS, we filtered single nucleotide polymorphisms (SNPs) and standardized effect size using the following criteria:
1) Remove the SNPs located in major histocompatibility complex (MHC) region.
2) Select the SNPs with a common frequency of the effect allele (EAF) (> 0.01 and < 0.99).
3) Standardize the effect size (β) and standard error (se) for each GWAS data with the function of minor allele frequency and sample size using the following formula :
where z = β/se from the original summary data, p is the minor allele frequency, and n is the total sample size.
Instrumental variables (IVs) selection
We selected independent and genome-wide significant GWAS SNPs of HDL-c, LDL-c, TG, TC and dyslipidemia by use of the clumping algorithm in PLINK (http://pngu.mgh.harvard.edu/purcell/plink/)  at a suggestive threshold (r2 threshold = 0.001, window size = 1 Mb, p-value = 5 × 10-8). The 1000 Genomes Project (http://www.internationalgenome.org/) data were used as the reference for linkage disequilibrium (LD) estimation. For each outcome, we harmonize the data according to the SNPs included in COVID-19 GWAS and their effect allele. After data harmonization, we then removed outlier pleiotropic SNPs using RadialMR  with the p-value threshold of 0.05. RadialMR  identified outlying genetic instruments via heterogeneity test (modified Q-statistics). After the removal of pleiotropy, the remaining exposure related SNPs for each outcome as instrumental variables (IVs) were utilized to perform MR analyses.
MR analyses and pleiotropy assessment
We conducted four complementary two-sample MR methods, including Inverse-Variance Weighted (IVW) method, weighted median method, weighted mode method and MR-Egger method, which make different assumptions about horizontal pleiotropy.
The IVW method assumes balanced pleiotropy . The pleiotropy is assessed via Cochran’s Q statistic and presented as excessive heterogeneity which will inflate the estimate of MR analysis . MR-Egger is based on the assumption which indicate instrument strength independent of the direct eﬀects . It can be evaluated by the regression dilution I2 (GX)  according to the assumption that no measurement error (NOME) in the SNP exposure eﬀects. I2 (GX) is an adaptative I2 statistic which propose to quantify the strength of NOME violation for MR-Egger method. If I2 (GX)  was sufficiently low (I2 (GX) < 0·9), the correction analysis was conducted to assess the causal effect by SIMEX, which can substantially mitigate adverse effects by simulation extrapolation. The intercept term of MR-Egger method can used for evaluating the directional pleiotropic effect . When the intercept is zero or its p-value was not significant (p-value > 0.05) were considered as non-pleiotropy. Moreover, we also used the Rucker’s Qʹ statistic  to measure the heterogeneity for MR-Egger method. If the difference Q − Qʹ is sufficiently extreme with respect to a χ2 distribution with the 1 degree of freedom, we indicated that directional pleiotropy is an important factor and MR-Egger model provides a better fit than the IVW method . All methods of two-sample MR analyses were measured by TwoSampleMR package in R. For various estimates for different measures, we select the main MR method as following rules:
1) If no directional pleiotropy in MR estimates (Q statistic: p-value > 0.05, MR-Egger intercept: intercept = 0 or p-value > 0.05, Q – Q’: p-value > 0.05), IVW method was used.
2) If directional pleiotropy was detected (MR-Egger intercept: intercept ≠ 0 and p-value < 0.05, Q – Q’: p-value < 0.05) and p-value > 0.05 for the test of Q’, MR-Egger method was used.
3) If directional pleiotropy was detected (MR-Egger intercept: intercept ≠ 0 and p-value < 0.05, Q – Q’: p-value < 0.05) and p-value < 0.05 for the test of Q’, weighted median method was used.
Leave-one-out sensitivity analysis was implemented to assess whether any significant results were generated by a specific SNP in IVW models.
MAGMA (https://ctg.cncr.nl/software/magma)  is commonly used for gene and gene-set analyses based on GWAS and genotype data. In order to explore the association of TC and COVID-19, we implemented MAGMA to identify genes and gene sets in which multiple SNPs show moderate association to TC without reaching the stringent genome-wide significance level. The genome-wide gene-based association study (GWGAS) is based on the model of multiple linear principal components regression and calculated the gene p-value using F-test . All 19427 protein-coding genes from the NCBI 37.3 gene definitions were employed for SNPs annotation. We mapped SNPs to genes by a defined window around each gene of 2kb away from the transcription start site (TSS) upstream and 1kb away from the transcription stop site downstream based on human reference assembly (GRCh37 or hg19) . The GWGAS analysis was performed to quantify the degree of association for each gene to TC and to compute the correlations between genes are estimated according to LD statistics. The LD reference was also from Phase 3 of 1000 Genomes.
Multiple testing correction
We employed false discovery rate (FDR) to address multiple comparisons issue and the adjusted p-value < 0.05 was used for judging significance. For MR estimates, we adopted an FDR control procedure for the susceptibility and severity of COVID-19 separately.