Coagulation factor V is a T-cell inhibitor expressed by leukocytes in COVID-19

Summary Clotting Factor V (FV) is primarily synthesized in the liver and when cleaved by thrombin forms pro-coagulant Factor Va (FVa). Using whole blood RNAseq and scRNAseq of peripheral blood mononuclear cells, we find that FV mRNA is expressed in leukocytes, and identify neutrophils, monocytes, and T regulatory cells as sources of increased FV in hospitalized patients with COVID-19. Proteomic analysis confirms increased FV in circulating neutrophils in severe COVID-19, and immunofluorescence microscopy identifies FV in lung-infiltrating leukocytes in COVID-19 lung disease. Increased leukocyte FV expression in severe disease correlates with T-cell lymphopenia. Both plasma-derived and a cleavage resistant recombinant FV, but not thrombin cleaved FVa, suppress T-cell proliferation in vitro. Anticoagulants that reduce FV conversion to FVa, including heparin, may have the unintended consequence of suppressing the adaptive immune system.


INTRODUCTION
Dysregulation of both the immune (Zhou et al., 2020) and coagulation systems (Bikdeli et al., 2020) occurs in severe COVID-19 infection. Immunological responses include T-cell lymphopenia, which we have found can persist for months after the initial illness (Bergamaschi et al., 2021). Coagulopathy is an important cause of morbidity and mortality in patients with COVID-19, and a marked increase in circulating Factor V (FV) activity has been reported in patients with severe COVID-19, associated with increased risk of thromboembolism (Stefely et al., 2020). FV is primarily produced by the liver, and circulates as a 330 kD inactive anticoagulant form composed of a heavy chain, B domain, and light chain. The B domain contains basic (aa 963 to 1008) and acidic (aa 1492 to 1538) regions, which are predicted to bind together, exposing the intervening motif on the surface of the folded structure. The B domain is required for the anticoagulant activity of FV (Cramer and Gale, 2012;Thorelli et al., 1998). Thrombin, and in the presence of phospholipids, Factor Xa (Monkovic and Tracy, 1990) sequentially cleave FV at Arg 709 , Arg 1018 and Arg 1545 , releasing a cleaved B domain to create a heavy and light chain dimer, designated as FVa, that can bind factor Xa and act as a pro-coagulant.
Production of FV by T-cells (Shen et al., 1993) and monocytes (Dashty et al., 2012) has been previously reported. We were interested to learn whether leukocyte derived FV may increase in and if so what the functional consequences of this may be. We report here that circulating and lung-infiltrating leukocytes are a source of increased FV in patients with COVID-19. We further show that intact FV, but not FVa, inhibits T-cell proliferation in vitro. We propose that neutrophil-, monocyte-and Treg-derived Factor V may be an important cause of the dysregulated T-cell response to SARS-CoV-2.
Correlation of the FV module eigengene expression with severity of COVID-19 is highly significant ( Figure 1B). FV module expression is similar to healthy controls in asymptomatic or mildly symptomatic individuals in the community. In hospitalized patients with mild disease FV module expression is elevated at presentation, but declines as patients recover ( Figure 1B). In patients with severe disease FV module expression is elevated at presentation and increased levels persists for several weeks. Analysis of the absolute number of the cell types that were found to express FV over time shows that CD4 cells and Tregs are suppressed early in disease, whereas neutrophil counts are elevated and remain high in patients with severe disease ( Figure S2).

Analysis of peripheral blood cells for expression of other coagulation factors showed expression of Factor
XIIIa and low levels of Factor XII but not other coagulation factors (data not shown). scRNAseq was performed on peripheral blood mononuclear cells (PBMCs) derived from 47 individuals recruited in Cambridge UK for whom FV data is available. These volunteers form part of a larger cohort of 130 volunteers in whom scRNAseq was performed (Stephenson et al., 2021). Tregs and monocytes express FV in patients with COVID-19. Tregs (CD4 + , FoxP3+) had the highest expression of FV, with expression also detected in monocytes and lower levels in other CD4 cell subsets ( Figure 1C).
Analysis of the Blueprint (https://www.blueprint-epigenome.eu/) database of hematopoietic genomescale datasets from healthy volunteers shows the highest level of FV expression in neutrophils, eosinophils, Tregs and monocytes ( Figure S3).

FV mRNA expression correlates with protein expression and parameters of disease severity and lymphopenia
To determine whether increased peripheral blood cell FV mRNA levels is associated with FV protein expression we assayed plasma FV levels and performed liquid chromatography -mass spectrometry on neutrophil lysates from healthy controls and patients with severe COVID-19.
There was a modest but statistically significant correlation between FV gene expression and FV protein expression (Figure 2A; p = 0.023, R = 0.17). Proteomic analysis showed significantly higher levels of FV in neutrophil lysates from patients with severe COVID-19 compared to healthy controls ( Figure 2B). We next explored whether FV expression correlated with biomarkers of disease severity. Consistent with our observation that FV gene expression increases in severe disease, we found that FV module gene expression correlates with predictors of disease severity (age, male gender, CRP) and increased plasma levels of IL6, IL10, interferon-g and TNF, whereas FV plasma levels correlate only with fibrinogen and IL6 ( Figures 2C i and 2C.ii). Furthermore, FV module gene expression correlates with suppression of T-cell counts, a marker of disease severity (Bergamaschi et al., 2021), during the first 24 days after symptom onset ( Figure 2C iv), and T and B cell counts after 24 days from the onset of symptoms ( Figure 2C vi). In contrast, there was very little correlation between plasma factor V levels and T and B cell counts during the first 24 days after symptom onset ( Figure 2C iii) or after 24 days from the onset of symptoms ( Figure 2C v).  (B) Weighted gene co-expression network analysis identified a module containing group of genes co-expressed with FV, in which FV is a hub gene and its expression correlates strongly with genes expressed in neutrophils ( Figure S1). Mixedeffects model with quadratic time trend showing the longitudinal expression of the FV module over time, grouped by severity. Gray band indicates the IQR of HCs. A significant effect of time versus severity group interaction term (p = 3.33e-07) indicates that disease severity has a significant effect on longitudinal expression. (C) scRNAseq of PBMCs derived from HC, HCW and patients with COVID-19 showed the highest expression of FV in CD4 + , FoxP3+ Tregs, with expression also detected in monocytes and at lower levels in other CD4 cell subsets. FV expression was increased in Tregs versus CD4 + Naive cells in healthy controls (p = 8.33 3 10 À4 ), and in severe COVID-19 versus healthy controls in CD14 monocytes (p = 0.016) but not other cell subsets. See also Figures S1, S3 and Table S1. Boxes denote IQR with median shown as horizontal bars. Whiskers extend to 1.5x the IQR; outliers are shown as individual points.   Figure 3B). This effect was prevented by thrombin, and enhanced by the thrombin inhibitor hirudin. The cleavage-resistant recombinant FV was a potent inhibitor of Tcon proliferation. Recombinant full length FV and a cleavage-resistant recombinant FV also suppress CD8 + T-cell proliferation, but not B cell proliferation ( Figure 3C).

FV is expressed by lung-infiltrating leukocytes in fatal COVID-19
The significant correlation of lymphopenia with leukocyte expression of FV rather than with plasma levels prompted us to determine if production of the protein at a site of infection outside of the bloodstream might be occurring. To address this, we stained lung tissue from autopsies of four patients who died with COVID-19 lung disease. Strikingly, FV was increased within the lung parenchyma in fatal COVID-19 and was associated primarily with infiltrating neutrophils and monocytes ( Figure 4). Some staining was also seen in a few CD3 + T-cells and in alveolar epithelia. Control lung showed only low levels of FV in epithelial cells.

DISCUSSION
Our data show that circulating neutrophils, monocytes and Tregs are a source of increased FV in patients with severe COVID-19, both within the bloodstream and within lung parenchyma, and that FV but not thrombin activated FVa decreases T-cell proliferation in vitro. In patients with COVID-19 increased circulating blood cell FV expression correlates with T-cell lymphopenia, which can persist for at least ten weeks after infection.
There was a modest but statistically significant correlation between FV gene expression in circulating blood cells and plasma FV levels, but a role for leukocyte derived FV in regulating the immune response to SARS-CoV-2 is likely to occur following migration of leukocytes into the tissues, including secondary lymphoid organs from which circulating FV is normally excluded.
Broncho-alveolar lavage fluid from patients with severe COVID-19 has shown a predominance of neutrophils, monocytes/macrophages, and eosinophils with few lymphocytes, and in particular CD4 and CD8 T-cell lymphopenia, but an increase in the proportion of Tregs (Ronit et al., 2020). Furthermore lymphocytes were either absent or sparse in areas of SARs-CoV-2 pneumonia infiltrated by macrophages or neutrophils, or containing neutrophil extracellular traps (NETs) composed of neutrophil DNA, histones and granule-derived enzymes (Radermecker et al., 2020;Veras et al., 2020). Thus leukocyte derived FV may suppress local T-cell proliferation at sites of infection. In support of this we found a correlation between plasma FV levels and fibrinogen, biomarkers of hemostasis, whereas T-cell counts correlated with FV module gene expression in circulating leukocytes. We also found expression of FV in lung cells expressing cytokeratin, which in some areas appeared to be forming syncytia. Syncytial formation has been reported in lung tissue of patients with severe COVID-19 (Braga et al., 2021;Bussani et al., 2020), and the presence of syncytia correlates with lymphopenia . Overall these observations would be consistent with the liver being the predominant source of plasma FV, whereas leukocytes are an important source of FV in tissues. Local production of FV by leukocytes may also contribute to thrombosis, and iScience Article neutrophil-rich macrothrombi have been found in heart autopsy tissue of patients who died from COVID-19 (Johnson et al., 2022).
A number of variants in the FV gene are associated with abnormalities of coagulation. The commonest is FV Leiden, a G1691A variant resulting in a R506Q amino acid change, which reduces the inactivation of FVa by activated protein C (Van Cott et al., 2016). Heterozygous FV Leiden is found in about 5% of Caucasians, and homozygosity in about one in 5,000. It is unclear whether FV Leiden alters the immune inhibitory properties of FV, but FV Leiden heterozygosity may reduce the risk of developing sepsis from infection (Yan and Nelson, 2004).
East Texas bleeding disorder is caused by a rare variant at A2440G mutation in exon 13, predicting a S756G mutation in the B domain of FV, which causes upregulation of an alternatively spliced F5 transcript that results in a 250kDa isoform known as FV-short (Vincent et al., 2013). FV-short is thought to inhibit coagulation by forming a high-affinity complex with the coagulation inhibitor tissue factor pathway inhibitor-a (TFPIa). It is unclear whether FV-short has immunomodulatory properties, but we found TFPIa is increased in severe COVID-19 in concordance with FV (data not shown). Immunostaining of COVID-19 lung tissue with anti-CD45 and anti-FV antibodies showed co-staining of CD45 + cells (red) with FV (green, vertical arrows on merged image), but also staining of some CD45 À cells with FV (horizontal arrows on merged image). Co staining with FV (green) was also seen in cells expressing elastase (neutrophil marker, red), CD68 (monocyte lineage marker, red) and some CD3 + cells (red). Cells showing co-expression are identified with vertical arrows on the merged image. Co-staining for FV (green) was also seen in cells expressing cytokeratin (red) in COVID-19 lung tissue, where they form syncytia with bi-nucleate cells (arrow), and to a lesser extent in normal lung tissue. Scale bar 15 mm.

OPEN ACCESS
iScience 25, 103971, March 18, 2022 7 iScience Article COVID-19 is associated with an increased risk of thromboembolism, but the role of anti-thrombotic agents remains unclear. In critically ill patients receiving organ support therapeutic dose heparin does not improve clinical outcomes or mortality, and may even cause harm (Goligher et al., 2021;Leentjens et al., 2021). Heparin binds to and activates anti-thrombin (Olson et al., 2010). LMWHs retain some antithrombin activity, but most of their activity is thought to result from inhibition of Factor Xa (Nutescu et al., 2016). Alternatives for VTE prophylaxis include warfarin and direct oral anticoagulants (DOACs). Warfarin depletes the reduced form of vitamin K that acts as a cofactor for gamma carboxylation of prothrombin, VII, IX and X, rendering them inactive (Wallin and Hutson, 2004). DOACs inhibit factor Xa or thrombin. The effect of anticoagulants on circulating factor V levels is unknown, but by inhibiting thrombin heparin reduces factor V activation, and could potentiate the T-cell suppressive effect of FV. In support of this there is evidence that heparin can directly suppress T-cell responses (Gorski et al., 1991). Increased FV expression by cells of the innate and adaptive immune systems may explain the lymphopenia seen in patients with severe Covid19. Heparin may potentiate suppression of the adaptive immune response by reducing FV activation.

Limitations of the study
Our study is limited by recruitment of the cohorts from a single geographical region in the UK.
Sample size is critical when studying a heterogeneous disease, and COVID-19 falls into this category. A larger sample would allow more detailed analysis of the effect of demographic or disease related factors, and allow study of the effect of pro-thrombotic variants in the FV gene. In addition, our patients were recruited during the first wave of the pandemic, and a study examining infection with new SARS-CoV-2 strains with different virulence, and in vaccinated as well as unvaccinated patients would be interesting. Finally, it is unclear from this study whether increased leukocyte FV is specific for SARS-CoV-2, or a more general effect that could contribute to lymphopenia in viral infections. The demographics, clinical and laboratory assessments of the participants have been previously described (Bergamaschi et al., 2021). 18 asymptomatic healthcare workers (group A) were 22.2% male and had a mean (SD) age of 32.9 (12.7); 40 symptomatic healthcare workers (group B) were 22.5% male and had a mean (SD) age of 36.0 (11.8); 46 hospitalised patients with mild disease (group C) were 54,3% male and had a mean (SD) age of 58.0 (16.9); 37 hospitalised patients requiring oxygen (group D) were 64.9% male and had a mean (

Peripheral blood cell isolation and culture
PBMCs were isolated from leukapheresis samples and whole blood by polysucrose density gradient centrifugation (Ficoll-Paque, GE Healthcare Life Sciences, UK). Leukapheresis samples are de-identified and age and gender of donor is not known. Cells were grown at 37 C in RPMI media supplemented with 10% human AB serum.
Neutrophils were isolated from blood by dextran sedimentation and discontinuous Percoll gradients (Reyes et al., 2021). Up to 80 mL of whole blood was collected into citrate tubes and centrifuged at 300 3 g (acceleration 5, deceleration 5) for 20 min and the platelet-rich plasma layer removed. Erythrocyte sedimentation and leukocyte-rich plasma were obtained by incubating the remaining contents in the tube with 6 mL of 6% Dextran 500 in saline and final volume adjusted to 50 mL with 0.9% NaCl for at least 20 min at room temperature. The leukocyte-rich portion was centrifuged at 350 3 g (acceleration 5, deceleration 5) for 6 min, with the pellet resuspended in 3 mL of 49.5% Percoll (GE Healthcare) and overlayed onto 61.2% Percoll and 72.9% Percoll. Gradients were centrifuged at 720 3 g (acceleration 1, deceleration 0) for 20 min to obtain neutrophil and PBMC layers.
To prepare proteomic samples neutrophils were centrifuged at 300 3 g for 5minat 4 C and resuspended in 7 mL of 0.2% NaCl (w/v in H 2 O) for 5minat room temperature and topped up with 7 mL of 1.6% NaCl (w/v in H 2 O). Cells were washed twice in Dulbecco's phosphate-buffered saline (DPBS; Thermo Fisher), pelleted at 300 3 g for 5minat 4 C and resuspended in 372 mL of freshly made 5% sodium dodecyl sulfate (SDS, BioRad) lysis buffer and vortexed. Samples were then heat denatured in a heat block for 5minat 100 C and stored at À80 C. Cell lysates were thawed and tris(2-carboxyethyl) phosphine hydrochloride (TCEP) and triethylammonium bicarbonate (TEAB) were added to a final concentration of 10 and 50 mM, respectively. Lysates were shaken at 500rpmat 22 C for 5 min before being incubated at 98 C for 5 min. Samples were allowed to cool and were then sonicated with a BioRuptor (30 cycles: 30 s on and 30 s off). Tubes were centrifuged at 17,000 3 g to collect the cell lysate and 1 mL of benzonase (27.8 units) was added to each sample and samples incubated at 37 C for 15 min. Samples were then alkylated with addition of 20 mM iodoacetamide for 1hat 22 C in the dark. Protein lysates were processed for mass spectrometry using ProTifi s-trap spin columns following the manufacturer's instructions. Lysates were digested with Trypsin at a ratio 1:20 (protein:enzyme) in 50 mM ammonium bicarbonate. Peptides were eluted from s-trap columns by ll

OPEN ACCESS
Poor quality cells were removed based on an excess of mitochondrial expressed genes, defined as > 7% of all UMIs in each single-cell. Sparsely-sequenced cells were also removed where the total UMI count for a cell was <1000. Deconvolution normalisation factors were then estimated for all cells across all samples combined, prior to log10 transformation with a pseudocount (+1), as implemented in scran . The probability of being a doublet was estimated for each cell per sample using the ''doubletCells'' function in scran based on highly variable genes (HVGs). Highly variable genes (HVG) were defined across all single-cells by modelling the mean-variance relationship across genes using a loess fit, as implemented in the modelGeneVar function in scran. Next, we used ''cluster_walktrap''https://arxiv.org/abs/physics/ 0512106) on the shared nearest-neighbour (SNN)-graph that was computed on HVGs to form highly resolved clusters per sample. Per-sample clusters with either a median doublet score greater than the median + 2.5 x MAD or clusters containing more than the median + 2.5 MAD genotype doublets were tagged as doublets. This was followed by a second round of highly-resolved clustering across the whole data set, in which again cells belonging to clusters with a high proportion (>60%) of cells previously labelled as doublets were also defined as doublets.
Collection and analysis of post mortem tissue De-paraffinized sections were exposed to high-pressure antigen retrieval before incubation with anti-FV and cell specific markers overnight (4 C). This was followed by specified-specific secondary antibody-conjugated to Alexa flour-488 or Northern lightÀ557 plus Hoechst 333342 (1 mg/mL) for nuclei detection for 1 h (1:100), mounted in VectaShield (Vector Laboratories) before viewing on a TCS-SPE CLSM (Leica Microsystems).Isotype-specific sera was used as a negative control. Image for each fluorophore was acquired sequentially using the same constant acquisition time and settings rather than simultaneously to avoid crosstalk between channels.

Expression of FV constructs
Three constructs comprising (1) FV B domain ((aa710-1545), (2) full length FV (aa 1 -2224) and (3) [R709A, R1018A, R1545A]FV aa1-2224 were sub-cloned into a proprietary vector for the HEK293-6E system by Peak Proteins. All sequences contained a C-terminal 6His tag to facilitate purification. Cells were transfected at a 500 mL scale for each construct, media harvested 5-6 days after transfection and protein purified using a combination of Ni affinity and size exclusion chromatography and if required ion exchange. Purified proteins were analysed by reducing and non-reducing SDS-PAGE, A280 to determine concentration, size exclusion and mass spectrometry to confirm identity. iScience Article Unless otherwise specified, longitudinally collected data was grouped by bins of 7 or 12 days. Pairwise statistical comparison of absolute cell counts, CRP or serum measures between individuals in a given severity group at a given time bin and HCs, or between severity groups, was conducted by Wilcoxon test unless otherwise specified. For analyses involving repeated-measures, false discovery rate corrected (Benjamini & Hochberg) p value were reported. For individuals sampled more than once within a given time bin, data from the earliest blood collection was used.

QUANTIFICATION AND STATISTICAL ANALYSIS
Weighted gene co-expression network and gene enrichment analysis The weighted gene co-expression network analysis (WGCNA) package in R overcomes the problem of multiple testing by grouping co-correlated genes into modules and relating them to clinic traits. Modules are not comprised ofa priori defined gene sets but rather are generated from unsupervised clustering. The eigengene of the module is then correlated with the sample traits and significance determined. A signed adjacency matrix was generated, and a soft thresholding power chosen to impose approximate scale-free topology. Modules identified from the topological overlap matrix had a specified minimum module size of 30. Significance of correlation between a clinical trait and a modular eigengene was assessed using linear regression with Bonferroni adjustment to correct for multiple testing. Modules were annotated using Enrichr and Genemania. Genes with high connectivity termed ''hub genes'' were identified based on a module membership of 0.8 or above and were selected to have a correlation with the trait of interest R0.8.
Gene enrichment analysis was performed using the Enrichr and GO enrichment analysis tools.

Correlation
The relationships between multiple features were quantified using Pearson's correlation (Hmisc package) and visualized with corrplot.

Mixed effects model
Longitudinal mixed modelling of factor V module eigenvalue changes over time ðy ij Þ was conducted using the nlme package in R, including time ðt ij Þ with a quadratic trend and disease severity ðX j Þ as fixed effects, and sampled individuals as random effects ðu j Þ:

Proteomic data analysis
Data were analysed with Spectronaut 14 using the direct data-independent acquisition option 27 (Skyline, MacCoss Lab Software is a freely available alternative). Cleavage Rules were set to Trypsin/P, Peptide maximum length was set to 52 amino acids, Peptide minimum length was set to 7 amino acids and Missed Cleavages set to 2. Calibration Mode was set to Automatic. Search criteria included carbamidomethylation of cysteine as a fixed modification, as well as oxidation of methionine, deamidation of asparagine and glutamine and acetylation (protein N-terminus) as variable modifications. The false discovery rate threshold was set to 1% Q-value at both the Precursor and Protein level. The single hit definition was to Stripped sequence. Data were searched against the human SwissProt database (July 2020) and included isoforms. The Major Group Quantity was set to the Sum of peptide quantity and the Minor Group Quantity was set to the Sum of the precursor quantity; Cross Run Normalization was disabled. Fold changes and p values were calculated in R utilising the bioconductor package LIMMA version 3.7 28. The Q-values provided were generated in R using the ''qvalue'' package version 2.10.0. Estimates of protein copy numbers per cell were calculated using the histone ruler method 29. The mass of individual proteins was estimated using the following formula: CN 3 MW/NA = protein mass (g cell À1), where CN is the protein copy number, MW is the protein molecular weight (in Da) and NA is Avogadro's Constant.

Single-cell RNA-sequencing clustering and annotation
Highly variable genes (HVG) were defined across all single-cells that passed QC and that were successfully assigned to a donor individual, by modelling the mean-variance relationship across genes using a loess fit, as implemented in the modelGeneVar function in scran. HVGs (FDR 1%) were used as input to estimate the first 50 principal components using the IRLBA implementation in the R package irlba. Single-cell clusters were computed by first constructing a k-nearest neighbour graph (k = 20). Cells were broadly grouped into discrete clusters based on the Walktrap community detection algorithm (https://arxiv.org/abs/ physics/0512106), using the kNN-graph as input.

ll
OPEN ACCESS