Proteome and phospholipidome interrelationship of synovial fluid-derived extracellular vesicles in equine osteoarthritis: An exploratory ‘multi-omics’ study to identify composite biomarkers

Osteoarthritis causes progressive joint deterioration, severe morbidity, and reduced mobility in both humans and horses. Currently, osteoarthritis is diagnosed at late stages through clinical examination and radiographic imaging, hence it is challenging to address and provide timely therapeutic interventions to slow disease progression or ameliorate symptoms. Extracellular vesicles are cell-derived vesicles that play a key role in cell-to-cell communication and are potential sources for specific composite biomarker panel discovery. We here used a multi-omics strategy combining proteomics and phospholipidomics in an integral approach to identify composite biomarkers associated to purified extracellular vesicles from synovial fluid of healthy, mildly and severely osteoarthritic equine joints. Although the number of extracellular vesicles was unaffected by osteoarthritis, proteome profiling of extracellular vesicles by mass spectrometry identified 40 differentially expressed proteins (non-adjusted p < 0.05) in osteoarthritic joints associated with 7 significant canonical pathways in osteoarthritis. Moreover, pathway analysis unveiled changes in disease and molecular functions during osteoarthritis development. Phospholipidome profiling by mass spectrometry showed a relative increase in sphingomyelin and a decrease in phosphatidylcholine, phosphatidylinositol, and phosphatidylserine in extracellular vesicles derived from osteoarthritic joints compared to healthy joints. Unsupervised data integration revealed positive correlations between the proteome and the phospholipidome. Comprehensive analysis showed that some phospholipids and their related proteins increased as the severity of osteoarthritis progressed, while others decreased or remained stable. Altogether our data show interrelationships between synovial fluid extracellular vesicle-associated phospholipids and proteins responding to osteoarthritis pathology and which could be explored as potential composite diagnostic biomarkers of disease.


Introduction
Osteoarthritis (OA) is the most prevalent arthritic phenotype and is one of the most important causes of perception of pain and loss of quality of life in the older population [1].OA has often been classified as a chronic degenerative joint disease resulting from a process of wear and tear.However, OA has an important inflammatory component, the mediators of which trigger an aberrant remodelling of joint structures inside the afflicted articulation [2].These may include synovial membrane dysfunction, abnormal bone proliferation, and subchondral bone sclerosis [3].Age, gender, obesity, genetics, inactivity, joint loading, aberrant morphology and alignment, previous injuries, and muscle weakness are the most prevalent risk factors for OA [4].OA in horses is a major cause of lameness, with over 60 % of lameness cases associated with a clinical diagnosis of OA [5].This results in impaired mobility, pain, poor performance, and early retirement, making equine OA a serious welfare issue that also leads to significant economic losses for the equine industry [6].
Previously, it has been shown that human and equine osteoarthritic pathogenesis follows a similar route from initial injury to disease progression and outcome, and as such, the horse is widely regarded as a clinically relevant model for musculoskeletal disease in humans [7].In addition, the horse's articular cartilage biology is anatomically comparable to that of humans with respect to both composition and thickness [8].The horse as a model for disease offers numerous further benefits, including the applicability of advanced diagnostic methodologies, such as magnetic resonance imaging (MRI) and arthroscopy, as well as serial sampling of biological material for analysis making it possible to monitor disease development, disease progression and response to intervention in great detail [7].
Presently, OA pathophysiology is not fully understood.The diagnosis is commonly based on clinical examination and radiographic imaging and, due to the insidious character of the disorder is often made at late stages when cartilage damage is already substantial and far exceeds the tissue's capacity for intrinsic repair [9].Therefore, it is paramount to identify biomarkers of disease that can be used to develop diagnostic tests that are both sensitive and specific for early OA, which could ultimately enable a timelier management of therapeutic interventions and decelerate disease progression.
In recent years, the concept of composite biomarkers has become popular; by definition, they are a non-linear combination of multiple measurements used to diagnose disease or predict outcomes [10].Thus far, they have been used in neurological diseases such as Alzheimer's disease and bipolar disorder [11], often using neuronal networks, artificial intelligence or machine learning algorithms.As such, extracellular vesicles (EVs) can be considered a biological source for composite biomarker discovery.
EVs are nanoscale-sized vesicles with a phospholipid bilayer membrane secreted by cells and specialised in restoring homeostasis or facilitating intercellular communication [12].Furthermore, EVs transport bioactive molecules that can elicit a response in recipient cells, resulting in physiological and phenotypic changes [13,14].They are present in tissues and body fluids, such as blood, urine and synovial fluid (SF) [15][16][17].It has been proposed that EVs may play a vital role in cartilage homeostasis and in the propagation of OA by promoting inflammation and regulating extracellular matrix (ECM) turnover [18][19][20][21].EVs are found in abundance in SF due to its close proximity to EV-secreting sources, such as native cells found within the joint space and periarticular tissues, including but not limited to chondrocytes and synoviocytes [22].For joint disorders such as OA, SF is thus the most appropriate source of biochemical information [20,21,23].
The translation of EV biomarkers to the clinic has been pioneered in the fields of cancer and neurodegenerative diseases [24,25].Nowadays, EVs are increasingly seen as a source for biomarker discovery for various disorders, including joint disease [18,20,26].A comprehensive understanding of the molecular composition of EVs and their role in disease requires the interpretation of molecular intricacy by accounting for multiple biological levels, such as the proteome and phospholipidome [27,28].Such a comprehensive experimental and data analysis approach provides a more thorough understanding of the complete spectrum of molecular changes contributing to cellular response, disease development and pathogenesis and is helpful for the identification of naturally occurring composite biomarkers.Recent studies in ovarian cancer [29] and Alzheimer's disease [30] have identified a relationship between the proteome and phospholipidome of EVs.
The hypothesis of our exploratory study was that the proteome and phospolipidomic cargo changed with OA severity.To investigate this we exploited omics-based technologies to analyse the proteome and phospholipidome of SF-derived EVs (SF-EVs) to 1) enable comprehensive profiling of a healthy state versus clinically diagnosed mild and severe OA in horses and 2) identify candidate composite diagnostic biomarkers of OA.

Materials and methods
An extended description of the methodologies used in this study can be found as supplementary information.

Ethical considerations
Equine SF was collected from horses presenting at the EQI VET SERWIS clinic in Buk, Poland, with various disorders of the locomotor system before the intra-articular application of a local analgesic as a standard part of the clinical lameness examination.Sample collection was approved by the University of Liverpool's Veterinary Research Ethics Committee (VREC1180).Ethical approval was not required in Poland, as the procedures were considered non-experimental clinical veterinary practices, in accordance with Polish and EU law (Dz.U. poz.266 and 2010-63-EU directive).

Sample collection
SF was collected via aseptic arthrocentesis from one metacarpophalangeal joint of each biological donor into a plain Eppendorf tube.Samples were spun at 2540×g at 4 • C for 5 min.The supernatant was then transferred to a new Eppendorf tube, snap-frozen in liquid nitrogen and stored at − 80 • C. A description of diagnostic criteria and classification methods can be found in the supplementary information.
The horses with no lesions in the joints (which featured locomotor abnormalities caused by disorders of other, unrelated structures) were classified as horses with healthy joints.Three biological replicates were pooled per sample resulting in 5 mL of SF.Pooled samples came from horses with healthy joints or with the same disease severity.Donors for the pooled samples were randomised with respect to age and sex.A total of 42 donors were used, resulting in 14 pooled samples (Healthy joints n = 7, mild OA n = 4, and severe OA n = 3).

Extracellular vesicle isolation and quantification
EVs from SF were isolated using a published and validated method [16].Samples of pooled cell-free synovial fluid (5 mL) were incubated with Hyase (5 mg/mL; Sigma-Aldrich, St. Louis, MO, USA), followed by centrifugation to remove protein aggregates and debris.Subsequently, the supernatants were processed through ultracentrifugation to isolate extracellular vesicle (EV) pellets, which were then resuspended in μL of PBS with 0.1 % Bovine Serum Albumin (BSA) depleted of EVs prior to use (EV-depleted BSA).Next, sucrose density gradient centrifugation was employed.Additional details of the EV isolation and EV-depleted BSA are provided in the supplementary information section.The EV-containing fractions, validated previously in Refs.[16,21], were pooled based on densities (1.10-1.16g/mL) for further lipidomics and proteomics analyses.Relevant data regarding the experimental details for EV isolation and characterisation have been submitted to the EV-TRACK knowledgebase (EV-TRACK ID: EV230607) [31].
EVs were labelled with PKH67, as previously described [21,32], involving resuspension in PBS+0.1 % EV-depleted BSA, the addition of PKH67 dye, and subsequent density gradient ultracentrifugation.A procedural control sample was analysed by high-resolution flow E. Clarke et al. cytometry to determine the background.A BD Influx jet-in-air flow cytometer optimized for single EV analysis was used for single particle analysis [32,33].The EV concentration was calculated based on fluorescent events in EV-enriched sucrose fractions F7-F10 (densities 1.10 g/mL-1.16g/mL), using BD FACS Software for data collection and FlowJo software for analysis.A comprehensive description is available in the supplementary section.The MIFlowCyt author checklist can be found as Suppl.Table 2 and the MIFlowCyt-EV framework as Suppl.Table 3 [34].

Lipid isolation
Lipids were extracted following the Bligh & Dyer method [35] with slight modifications, as described previously [21].Isolated EV samples were mixed with methanol and chloroform for the separation of hydrophilic and hydrophobic phases.The extracted lipids were dried and stored in a nitrogen atmosphere.A more extensive description is accessible in the supplementary section.During the lipid extraction, one sample (composed of n = 3 individual horses) from healthy joints was lost.Therefore, n = 6 SF-EV samples of the group with healthy joints were used for all lipidomics analyses and subsequent omics integration.

Mass spectrometry lipidomics
Lipid pellets obtained from dried samples were resuspended in chloroform/methanol and analysed using hydrophilic interaction liquid chromatography (HILIC) coupled with a Fusion Orbitrap mass spectrometer as described previously [21,36].A quality control sample, along with a lipid standard, was included for quantification and verification of the mass spectrometry run quality.Thorough details are found in the supplementary materials.

Lipid annotation
RAW data was converted to mzML format, and LC/MS peak-picking and retention time correction were done using XCMS in R to annotate identified peaks based on retention time and exact m/z-ratio.Annotation criteria included presence in at least 3 out of 13 pooled samples, utilizing an in-silico phospholipid database, and adjusting for isotope overlap, with a focus on major adducts to prevent under-quantification of less prominent lipid species.Detailed information can be found in the supplementary section.The RAW and mzML converted mass spectrometry data is deposited in the YODA repository of Utrecht University [37].

SDS PAGE & silver stain
Proteins from EV protein extracts were separated using sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE).The procedure involved adding Novex™ Tris-Glycine SDS Sample Buffer to the sample, denaturing proteins through heating, electrophoresis using a NuPAGE™ gel, and visualization using silver stain according to the manufacturer's guidelines as previously done [23].Additional details are provided in the supplementary information section.

On bead digestion
Hydrophilic and hydrophobic magnetic beads were employed to digest EV proteins, facilitating the removal of the incompatible urea lysis buffer.Following the treatment of lysed equine synovial fluid-derived EVs with reducing and alkylating agents, the EV proteins were extracted using magnetic carboxylate SpeedBeads, subjected to trypsin digestion, and subsequently prepared for mass spectrometry analysis through desalting and normalisation steps [23].A comprehensive explanation is available in the supplementary section.

Data-dependent acquisition for the generation of an equine SF EV spectral library
Equine SF was pooled using samples from the metacarpophalangeal joint from our equine musculoskeletal biobank (VREC561) and samples collected in previous studies from the carpal and metacarpal joint of healthy horses as well as those with OA, resulting in a total of 11 ml pooled SF (n = 1) for library generation.These samples were analysed as previously described in order to generate the necessary reference library [38].

Data-independent acquisition proteomics (SWATH)
Data-independent acquisition proteomics (SWATH) was employed using a Triple TOF 6600, with a 2-h gradient (as the library fractions [23]) and a comprehensive precursor m/z range of 400-1500.Retention time alignment and quantification were conducted through Data-Independent Acquisition by Neural Networks (DIA-NN), ensuring a precursor false discovery rate (FDR) of 1 %.The mass spectrometry proteomics data were deposited to the ProteomeXchange Consortium via PRIDE proteome exchange [23] (identifier PXD042765).Both proteomics and lipidomics datasets have been submitted to vesiclepedia [39].

Statistical analysis 2.6.1. EV characterisation
Comparison in EV concentration (fluorescent events/mL) between the mean of two groups, healthy and OA, was done by using a Student's t-test.

Proteomics
Statistical analysis of proteomics data was carried out using the R statistical programming environment or Metaboanalyst [40].The data were quality controlled; proteins with complete observations were normalised using probabilistic quotient normalisation (PQN) and log-transformed (base 10) for downstream analysis, where a normal distribution of the data is a prerequisite for functional enrichment analysis using the Ingenuity Pathway Analysis (IPA) software.Unsupervised multivariate analysis in the form of principal component analysis (PCA) was performed, along with heat map analysis using analysis of variance (ANOVA) and Pearson distance.One-way ANOVA with Tukey's post hoc test was attributed to statistical significant proteins in their respective group comparison.Following ANOVA, a fold change analysis was conducted.

Lipidomics and omics data integration
For lipidomics analysis, the data were normalised based on the sum of total lipids per pool samplei.e. each lipid value in a pooled sample was divided by the total sum of lipids in the same pool sample and multiplied by 0.01; thus, the relative abundances sum up to 100.A minimum of three biological-pool replicates were used for statistical analyses.
Data analysis was run with R version 4.1.2[41].Pareto scaling was performed for the PCA, thus dividing each variable by the square root of its standard deviation.Heatmap and cluster analysis was performed on Spearman correlations with a set speed of twoamong the 50 most abundant lipid species in all sample groupsusing the R-package ComplexHeatmap v1.12.0 [42].
Data integration was performed with the R package mixOmics v6.12.2 [43].on lipidomic and proteomics data normalised by the sum E. Clarke et al.
(as described for lipidomics analysis) followed by R scaling and centring, which determines the vector's mean and standard deviation, deducts the mean from the vector and divides it by the standard deviation.An unsupervised sparse Partial Least Squares (a linear, multivariate regression method for data reduction to assess the relationship between independent and dependent variables) was used to integrate the datasets.The relevance network plot was set with a correlation cut off of 0.7 to allow readability of the displayed proteins and phospholipids.Differences between the proposed proteins and phospholipid percentages for the composite biomarker were analysed with the rank-based non-parametric Kruskal-Wallis test (since the data remained in this scenario skewed and thus non-parametric), followed by the multiple pairwise comparisons with Dunn's test.Significance was set at α = 0.05.Statistical tests were done with GraphPad Prism 9.

Functional enrichment analysis
Functional enrichment analysis was performed on proteomic data using Ingenuity Pathway Analysis (IPA; Qiagen, Hilden, The Netherlands).In order to provide functional analyses, networks, canonical pathways, and related molecular and pathological functions the protein p-values obtained through One-way ANOVA with Tukey's post hoc test, and associated log2 fold change, were used.UniProt_Horse accession codes were used as protein identifiers, and the Qiagen Ingenuity Knowledge Base was used as a reference for exploratory pathway analysis.For network generation, default settings were used to identify molecules whose expression was significantly differentially regulated.These molecules were overlaid onto a global molecular network contained in the Ingenuity Knowledge Base.Networks of 'network-eligible molecules' were then algorithmically generated based on their connectivity.The functional analysis identified the biological functions and diseases that were most significant to the data set.A righttailed Fisher's exact test was used to calculate p-values.Canonical pathway analysis identified the pathways from the IPA library that were most significant to the data set.Analysis was performed on all proteomics data, comparing healthy, mild OA and severe OA groups, and those proteins correlated to phospholipids.

Synovial fluid-derived extracellular vesicle numbers do not significantly differ between osteoarthritic and healthy phenotypes
Recently we found that an inflammatory insult in the joint, such as LPS, can strongly affect the quantity of SF-EVs [21].Therefore, we investigated, using the same technology, if the quantity of SF-EVs was altered as a result of OA using samples from equine patients with radiographically diagnosed OA and comparing these with samples from healthy joints.The quantity of EVs was assessed by single-EV fluorescence-based flow cytometric analysis of PKH-labelled EVs [21,33] on 3 representative samples of the group with healthy joints, 2 samples of the mild OA group and 1 of the severe OA group.The PKH + events were measured in individual sucrose fractions ranging from 1.08 to 1.18 g/mL.The peak of fluorescent events was identified in the densities from 1.10 to 1.16 g/mL (Fig. 1A); those were considered the EV-enriched fractions and were used for calculating EV numbers (Fig. 1B).We did not observe statistically significant differences with a p = 0.47 between the numbers of SF-EV from the healthy joints group (where each sample consisted of SF-derived from 3 different horses) (3.1 × 10 8 per mL SF ± 6.6 × 10 7 ; mean ± SD) and the OA group (i.e.mild OA n = 2 and severe OA n = 1, each consisting of SF-derived from 3 different OA-diagnosed horses with the respective severity degree of OA) (4.0 × 10 8 per mL ± 9.4 × 10 7 ; mean ± SD).

Lipidomic analysis 3.2.1. Synovial fluid-derived extracellular vesicle phospholipid profiles change during the development of osteoarthritis
Previously we had observed a drastic change in the phospholipidome following an inflammatory stimulus [21]; here we analysed whether the phospholipid profile of the SF-EVs was modified as a result of OA.The phospholipidome profile of the SF-EV from healthy joints (n = 6), mild OA (n = 4), and severe OA (n = 3) equine patients was determined through a bioinformatics analysis that uncovered 280 lipid species after lipid annotation (and background adjustment), isotope and adduct correction and normalisation by the cumulative sum to unity (Supplementary Fig. 1).A PCA, an unsupervised dimensionality reduction Fig. 1.Quantitative flow cytometric analysis of EVs isolated from equine joints with a healthy or osteoarthritic phenotype.A) Single EV-based high-resolution FCM of representative healthy SF-EVs (n = 3) and OA SF-EVs (n = 2) from the mild OA group and n = 1 from the severe OA group.Sucrose density gradient fractions containing EVs labelled with the lipophilic dye PKH67 were measured for 30 s.The majority of EVs floated at densities of 1.16-1.10g/mL.FL -Events: Fluorescent Events.B) EV concentration in SF was calculated as the sum of single fluorescent events measurements (PKH67+ events) in EV-containing sucrose gradient densities (1.16-1.10g/mL).Mean ± SD. p = 0.47.ns: non-significance by Student's t-test.The uppermost point in the OA group reflects the severe OA phenotype.method, revealed a combined explained variance of 69 % with the first and second principal components (Fig. 2A).
A Spearman correlation heatmap showed that the EV populations of the three different clinical groups differed in the distribution of their phospholipid composition (Fig. 2B).The heatmap was split into three sections (slices and clusters) based on Partitioning Around Medoids (PAM) clustering.For the first slice, the predominant lipid classes were phosphatidylserine (PS), ester-linked phosphatidylcholine (PC), esterlinked phosphatidylethanolamine (PE) and ether-linked phosphatidylethanolamine (PE-O), which account for half of the lipid distribution of EVs in the healthy joints but for less in both OA groups.The second slice included other members of the PS and PC classes, and the phosphatidic acid (PA), lysophosphatidylcholine (LysoPC), and phosphatidylinositol (PI) lipid classes.There were no clearly identified clusters in this slice.The third slice consisted solely of sphingomyelin (SM), and the distribution was one-third per group; thus, the OA-derived EVs had a higher presence than the EVs from healthy joints.These results showed a subtle variance among SF-EVs from the healthy and the mild and severe OA phenotypes.

Differences in lipid class composition of synovial fluid extracellular vesicles are related to osteoarthritis progression
Having established a difference between the SF-EVs from healthy joints (n = 6) compared to mild OA (n = 4) and severe OA SF-EVs (n = 3), we proceeded to analyse in more detail how the lipid classes were distributed in the respective groups (Fig. 2C, Suppl.Fig. 2).The most abundant phospholipid classes in all three clinical groups were SM (20-40 %), PC (25-40 %), PS (12-16 %), PE O-(8-12 %) and PE (5-10 %) (Fig. 2C).However, a relative increase of SM was observed in the OA groups (healthy 19.9 %, mild OA 35.5 % and severe OA 37.5 %), while the amounts of PC, PI and PS relatively decreased in OA groups which was most pronounced in the severe OA group (healthy: PC 38.Despite variations in the total lipid classes with respect to the whole phospholipidome, the individual lipid species contributing to the lipid classes were similarly distributed throughout the clinical groups following normalisation within each class (Suppl.Fig. 2).Thus, the observed shifts in lipid classes cannot be directly attributed to changes in individual lipid species.Overall, these findings demonstrate that the phospholipidome is gradually transformed as OA develops.

Proteomic analysis 3.3.1. Principal component analysis of proteomics demonstrates variable protein distribution according to osteoarthritic phenotype
Unsupervised multivariate analysis using PCA was conducted on the proteome of all samples exploring the variability between SF-EVs derived from healthy joints (n = 7), mild OA (n = 4) and severe OA (n = 3).A total of 5774 unique peptides were identified, translating to 290 proteins with no missing values.Missing values as such were imputed (using impute 2,1,1) using the following method: For the 7 healthy samples, up to 2 missing values were imputed by inserting the mean of the healthy values for that particular protein.Similarly, for the 4 mild OA and 3 severe OA samples, up to 1 missing value was imputed by inserting the mean of the mild OA or severe OA values for that particular protein, resulting in a total of 598 proteins identified and quantified across all samples and used for statistical analysis (Supplementary Fig. 1).The first two components (Fig. 3A) reduce the total variation of all the individual data points by 36.4 %.

Differentially expressed proteins identified across osteoarthritic phenotypes
Using ANOVA, 40 proteins were identified as being significantly differentially expressed (p < 0.05) prior to false discovery rate (FDR) adjustment across all experimental groups (SF-EVs from healthy joints (n = 7), mild OA (n = 4) and severe OA (n = 3)).Uncorrected values were used due to this being an exploratory study, whereby multiple testing correction methods can fail to identify statistically significant values due to stringent thresholds [44].Following Tukey's post hoc analysis, a remaining 37 were significant (p < 0.05).Table 1 demonstrates the top 25 differentially expressed proteins and their respective fold change expression, as well as the specific comparison found to be significant following Tukey's post hoc tests.It was revealed that microtubule-associated protein (ANOVA p = 0.006, Tukey test: severe OA compared to mild OA p = 0.006, and mild OA compared to healthy p = 0.03) was present at higher levels in mild OA compared to the severe form of the disease.Further proteins with an increased expression in severe OA compared with the group with healthy joints and mild OA were fibroblast activation protein alpha (ANOVA p = 0.03, Tukey test: severe OA compared to healthy p = 0.03) and Interleukin 1 receptor accessory protein (ANOVA p = 0.02, Tukey test: severe OA compared to healthy p = 0.02).Conversely, platelet-activating factor acetylhydrolase IB subunit alpha (ANOVA p = 0.004, Tukey test: severe OA compared to mild OA p = 0.003) exhibited increased expression in mild OA but was decreased in severe OA, as shown in Table 1.Other significant (p < 0.05) proteins attributed to EVs that were identified in our dataset included RAB GTPases, such as RAB GDP dissociation inhibitor (ANOVA p = 0.03, Tukey test: severe OA compared to healthy p = 0.02) and RAB8 (ANOVA p = 0.004, Tukey test: severe OA compared to mild OA p = 0.005, and mild OA compared to healthy p = 0.04).Overall, a change in the proteome was observed in response to an altered OA phenotype, with significant proteins attributed to pathways known for propagating OA disease development within the joint.

Functional enrichment analysis of the synovial fluid-derived extracellular vesicles proteome using IPA highlights dysregulation in pathways associated with cartilage homeostasis and an inflammatory phenotype
Functional enrichment analysis was performed in order to provide biological meaning to the identified and quantified proteome.In both mild OA and severe OA groups, the top canonical pathways were identified using the Ingenuity Knowledge Base Library and accounting for protein p-value following Tukey's post hoc analysis and log2 fold change.It was found that signalling by Rho family GTPases (p = 0.0000244), liver x receptor/retinoid x receptor (LXR/RXR) activation (p = 0.0353), complement system activation (p = 0.0107-0.000239),clathrin-mediated endocytosis (p = 0.000329), and macrophage alternative action signalling (p = 0.00604) were all significant to OA pathology when considering EV cargo, as shown in Fig. 4A, B, and C. Additionally, significant diseases and functions in both severe and mild OA included inflammation of an organ (p = 0.04).Molecular functions found to be significant in severe OA compared to mild included injury of joint (p = 0.00552), complement activation (p = 0.0171) and accumulation of macrophages (p = 0.000414).Disease and molecular functions identified in a severe OA compared to healthy included: fibrosis (p = 0.0189), systemic inflammation (p = 0.0211), acute inflammation of tissue (p = 0.0289) and osteoarthritis (p = 0.0316).Finally, mild OA compared to healthy identified significant functions including: complement activation (p = 0.000609), development of articular cartilage (p = 0.00374), injury of joint (p = 0.0118), inflammation of joint (p = 0.0169), osteoarthritis (p = 0.0280) and chronic inflammation (p = 0.0338).A complete list of significant diseases and function can be found in Suppl.Table 4.

Multi-omic integration 3.5.1. Proteomics and lipidomics data integration demonstrates a high correlation between proteins and phospholipids in synovial fluid extracellular vesicles
Integration of the proteome and phospholipidome datasets was performed to determine if biologically feasible correlates could be established; and thus, identify candidate composite protein-lipid biomarkers.An unsupervised approach was selected to integrate the dataset, consisting of a PCA assessment followed by sparse Partial Least Squares (sPLS2) regression which was tuned by cross-validation.
The initial exploratory analysis employing PCA was undertaken to recognise how the individual proteomic and lipidomic datasets behaved under the same normalisation conditions and to determine the optimal data integration model (Suppl.Figs.3A and 3B).The omics datasets were normalised by the summed intensity of the sample, followed by centring and scaling of the data, thus subtracting the mean and dividing by the standard deviation.It was observed that clustering of the samples was comparable to the previous PCA (Figs. 2A and 3A).Subsequently, to integrate the omics data sets, the unsupervised sPLS2 model was constructed separately for the proteomics and lipidomics data (Suppl.Figs.3C and 3D).As an unsupervised analysis, the information about the groups (healthy joints (n = 6), mild OA (n = 4) and severe OA (n = 3)) was not taken into consideration; however, the samples were labelled to understand how they clustered.In Suppl.Figs.3C and 3D, both sPLS2s project the respective data similarly, with the superior subspace primarily composed of SF-EV samples from healthy joints, the inferior one of mild OA and severe OA SF-EV samples and the top left subsection of overlapping samples from all groups.Afterwards, both sPLS2s were averaged (Fig. 5A).The integrated averaged sPLS2 had a similar structure in components as the individual sPLS2.Fig. 5B assesses the degree of agreement between the proteomic and lipidomic datasets by plotting the position of each sample from both sPLS2s in the same space and connecting them with an arrow that indicates at its base the location in the proteomics data set and at the tip the location in the lipidomic data set.Most samples were located relatively close to each other indicating a correlation between the phospholipidome and proteome of SF-EVs.
This correlation was further explored with a Cluster Image Map (CIM) (Fig. 5C) to examine the connection between the features and components in a broad range, drawing attention to the relevant variables that collectively accounted for the covariance between the two datasets.According to the CIM, the phospholipid variables were divided into three slices that were either positively or negatively related to two main protein clusters.The left slice corresponded to 4 SMs (SM 36:0; 2,  SM 43:3; 2, SM 41:3; 2, SM 40:3,2), PC O-32:3, and PC 40:9, which had a positive association with the lower protein cluster.The middle slice, consisting of three PC species (PC 37:3, PC 28:0 and PC 31:1), had an inverse pattern of the cluster, with the upper group depicting the strongest association.Finally, the right slice, comprising the PC 34:4 and the two PI species (PI 32:1 and PI 38:6), had a similar association pattern as the middle one; however, the lower cluster exhibited a negative correlation, while the cluster above was positively correlated to the proteins.

Relevance network for the selection of candidate proteins and phospholipids as composite OA biomarkers
To better comprehend the correlation between the proteins and phospholipids, a relevance network plot was created (Fig. 6).Three substructures could be identified from the network.The larger cluster contained the same lipids as the middle slice from the CIM (Fig. 5C; PC 28:0, PC 31:1 and PC 37:3), with all the correlations depicted being positive.The second substructure consisted of the right-side slice lipids from the CIM (Fig. 5C; PC 34:4, PI 32:1, PI 38:6), with primarily positive correlations to the proteins except to the anion exchange protein.This cluster also overlapped with some of the same proteins as PC 28:0, PC 31:1 and PC 37:3.The third substructure was composed of the lipids from the left-side slice of the CIM (SM 36:0; 2, SM 43:3; 2, SM 41:3; 2, SM 40:3,2, PC O-32:3, PC 40:9).This cluster displayed only positive correlation with the depicted proteins, including the anion exchange protein.Moreover, the proteins that correlated to phospholipids from the relevance network plot (Fig. 6) were found to be associated with pathways such as actin cytoskeleton signalling (p = 5.71 × 10 − 7 ) and  signalling by Rho family GTPases (p = 5.89 × 10 − 8 ), as shown in Table 2.
Since the sPLS2 analysis is an unsupervised approach (i.e., no information regarding the groups is entered in the model), neither the CIM nor the network explained how the SF-EV phospholipids and the correlated proteins relate to the healthy joints, mild OA and severe OA.
To determine differences between the clinical groups, all lipids and proteins with a correlation above 0.754 based on the network (Fig. 6) were assessed with a Kruskal-Wallis test (Fig. 7).A significant decrease in PC 34:4 and PI 38:6, and decline in PI 32:1 and the related proteins showed a similar trend in SF-EVs derived from severe OA compared to healthy joints (Fig. 7).Conversely, SM 36:0; 2, SM 41:3; 2 and PC O-32:3 and the correlated proteins showed a trend to increase with the severity of OA, with significant differences for PC O-32:3, moesin and vasodilator-stimulated phosphoprotein (Fig. 7).A list of potential candidate proteins for composite biomarkers is provided in Suppl.Table 5.Overall, we here show a strategy for composite biomarker discovery based on SF-derived EV-associated phospholipids and proteins and revealed potential candidates that could be explored as composite phospholipid-protein EV biomarkers in OA pathology.

Discussion
In this exploratory study we designed a workflow for a multi-omics approach based on phospholipidomic and proteomic integration to identify composite SF-derived EV-biomarkers for OA, based on the analysis of EVs isolated from SF of horses with clinically defined OA (mild OA and severe OA), or from healthy equine joints.Hereto we investigated the phospholipidome and proteome of purified SF-EVs and designed a strategy for multi-omics data integration and differential expression analysis.To identify genuine composite EV-biomarkers, we used differential centrifugation followed by sucrose density gradient centrifugation to purify EVs from SF by removing most types of the contaminating lipoproteins and protein aggregates [16,21].While the numbers of EVs in SF were unaffected by OA, consistent with other studies [45,46], the proteomic and phospholipidomic profiles of SF-EVs were correlated to the presence of OA.
We found that OA pathology directly impacted the phospholipidome at the lipid class level, showing a relative gradual changes in several lipid classes associated with disease severity.The relative reduction in PC, PS and PI in mild OA and more drastically in severe OA, could be explained by the relative increase in SM.Since PC and SM are primarily located in the outer layer of the plasma membrane, the increase in SM disrupts the balance and reduces the amount of PC [47].Similarly, although PS and PI are predominantly found in the inner leaflet of the lipid bilayer, they can also be affected by an increase in SM.Furthermore, we found relatively higher levels of PC compared to PC O-, while the PE and PE O-classes showed an opposite trend.These findings align with the lipidomics findings in the EV field and highlight the importance of ether lipids, especially PE O-, in EV biology, including membrane trafficking and cholesterol regulation [48].SM, one of the main lipid classes detected in the SF-EVs, plays a crucial role in the plasma membrane composition, cellular proliferation, differentiation, growth, signal transduction, and apoptosis [49].SMs are instrumental in the formation of lipid rafts enabling the selection of membrane proteins involved in signal transduction and intracellular transport [50].The notable relative increase of SMs with OA severity suggests that more lipid raft-like domains may be present in SF-EVs as the OA pathology progresses, facilitating and enhancing the cell-to-cell communication of SF-EVs.
Functional enrichment analysis of the differentially expressed SF-EV proteins identified a range of activated canonical pathways associated with disease phenotype.Specifically, Rho family GTPases, including RAC family small GTPase1 and ezrin were identified as activated in severe OA compared to healthy joints.Dysregulation of Rho GTPases has been implicated in rheumatic disorders in humans like rheumatoid arthritis, OA, and psoriatic arthritis, contributing to hypertrophic changes and cartilaginous matrix destruction [51][52][53].Rac1, a pro-inflammatory factor, stimulates MMP13 production and upregulates markers of chondrocyte hypertrophy, such as COLX and ADAMTS-5 [52].Dysregulated activation of Rho GTPases, particularly CDC42, can lead to the degradation of articular chondrocytes through IL-6/STAT3 signalling [54].The presence of these proteins in SF-EVs from diseased groups suggests their potential role in propagating disease within the joint by carrying cargo that induces phenotypic and metabolic changes.
Functional enrichment analysis also revealed disease and molecular functions related to complement system activation and macrophage alternative action signalling, with proteins such as complement C6 and ATP citrate lyase attributed to such pathways respectively.In fact, the complement system activation has previously been attributed to OA pathology, with its activation implicated in the formation of terminal complement complex (TCC) on chondrocytes, resulting in cell death, or the initiation of the production of matrix degrading enzymes, such as MMP13 [55,56].In previous studies an imbalance of macrophage subtypes (M1 and M2) has been proposed to contribute to the chronic low-grade inflammation associated with OA and to be implicated in OA pain mechanisms [57,58].In addition, in this study, it was found that LXR/RXR activation was implicated in severe OA phenotypes compared with mild OA with proteins such as inter alpha trypsin heavy chain 4. In previous studies it has been shown that LXR/RXR signalling is dysregulated in OA tissue and associated with inflammation [59], and has been identified in early and late stage OA [60].Additionally, a reduction in LXR signalling has been found to contribute to catabolic processes in OA in human articular cartilage [61].Hence, our findings suggest that a significant involvement of the immune system in the later stages of OA pathogenesis is reflected in the SF-EV proteome.
Overall, the observed changes in both phospholipid classes and proteins between SF-EVs derived from healthy joints and OA patients and the gradual changes associated with the severity of OA, suggest that these SF-EV parameters may be used as natural composite biomarkers for OA diagnosis and progression.Our multi-omics integration approach, using unsupervised sPLS2 regression and PCA, indeed revealed a remarkably strong similarity in the space distribution induced by the SF-EV phospholipidome and proteome, indicating a strong interrelationship, which is mainly due to a strong correlation between specific phospholipids with a certain set of proteins.Functional enrichment analysis of the proteins from this correlation network revealed several canonical pathways, such as signalling by Rho Family GTPases as previously identified and actin cytoskeleton signalling.
Integration of data revealed potential composite biomarkers consisting of downregulated and upregulated phospholipids and proteins as OA severity progressed.Downregulated proteins and the respective phospholipids were comprised of phospholipids PC 34:4, PI 32:1 and PI 38:6, and proteins such as heat shock protein 90 (HSP90AA1) and CD163.Interestingly, HSP90AA1 has been demonstrated to be downregulated in blood and cartilage of human patients with OA, and levels correlated with the risk incidence of OA [62], while CD163, a transmembrane protein of M2 macrophages [63], was shown to decline as OA progressed in this study.It has been suggested that the inability of macrophages to transform from M1 to M2 might contribute to the onset and development of OA [64].
Among the upregulated proteins, several structural proteins were detected, including α-2 smooth muscle actin, erythrocyte membrane protein band 4.1-like 2 (EPB41L2), ezrin, and moesin.These proteins likely indicate changes in diseased joint tissues, which were reflected in the structural protein composition of SF-EVs.α-smooth muscle actin is known to be expressed in fibroblast-like synoviocytes (FLSs) undergoing a change to a myofibroblast-like phenotype in the presence of transforming growth factor β (TGFβ), linked to OA pathogenesis [65], as well as to colocalise with fibronectin, which is associated with inflammation in OA [66].Ezrin, moesin, and EPB41L2 activation promote enhanced proliferation and formation of fibrillated OA cartilage by blocking cell-cell contact inhibition in chondrocytes [67].Ezrin has also been connected to the RhoGTPase signalling pathway in OA synovial fluid [67].Additionally, cluster of differentiation 90 (CD90) and CD109 transmembrane proteins, upregulated as OA progresses, regulate the pathological response in rheumatoid arthritis (RA) fibroblast-like synoviocytes, driving inflammation and fibrosis [68,69].The upregulated proteins were associated with phospholipids SM 36:0; 2, SM 41:3; 2, and PC O-32:3.The combinations of these proteins and phospholipids could potentially serve as candidate composite SF-EV biomarkers for OA onset and progression.The inherent constraints of this exploratory study include a relatively small clinical sample size and a large volume of SF required.In future E. Clarke et al. studies, it is important to properly identify composite biomarkers that can be translated to a clinical setting.Moreover, each sample was composed of a pool of three other horses rather than individual donors, which can affect the composition of the SF-EVs.Furthermore, overcoming challenges such as non-conformity in radiological and clinical parameters for OA severity assignment and the lag in developing analytical tools for comparing mass spectrometry proteomics and lipidomics pipelines is crucial in future studies.Nonetheless, the approach of our exploratory study in equine OA highlights the potential for identifying important molecular mechanisms of OA and aims to serve as a framework for the discovery of SF-derived EV-based composite biomarkers having the potential to inform disease severity and enable targeted disease management in the future.

Fig. 2 .
Fig. 2. Lipidomic profile of equine synovial fluid-derived EVs from healthy joints or from mild OA or severe OA patients.Healthy samples (n = 6), mild OA (n = 4), severe OA (n = 3).Each sample is comprised of a pool of three different animals.Lipids were extracted from EVs isolated by differential centrifugation up to 100,000 g, followed by purification with sucrose density gradients.A) Principal component analysis of lipids isolated from the three different clinical groups.The principal components (PC)-1 and − 2 explain 49 % and 20 % of the variance, respectively.Healthy samples (green circle), mild OA (orange circle), and severe OA (red circle).B) Lipid species correlation of SF-EVs.Combined heatmap (cluster dendrogram) of Lipid-Lipid Spearman correlations between the 50 most abundant lipid species in all EV sample groups.Lipid order was based on Partitioning Around Medoids, also known as K-Medoids, a centroid-based clustering algorithm.On top of the figure is the cluster dendrogram.Below is the group distribution, the relative lipid intensity of each species, and the heatmap.Under the heatmap, the degree of saturation, the lipid class of each lipid, and the respective annotation of each lipid species are indicated.C) Changes in EV lipid classes during OA development.Vertical slices plot of SF-EVs showing the relative molar abundances for individual lipid classes.Abbreviations: LysoPC, (lysophosphatidylcholine); LysoPG, (lysophosphatidylglycerol); LysoPI, (lysophosphatidylinositol); LysoPS, (lysophosphatidylserine); PC, (ester-linked phosphatidylcholine); PC O-, (ether-linked phosphatidylcholine); PE, (ester-linked phosphatidylethanolamine); PE O-, (ether-linked phosphatidylethanolamine); PI, (phosphatidylinositol); PS, (phosphatidylserine); SM, (sphingomyelin).(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) 3 %, PI 1.73 % and PS 16.0 %; mild OA: PC 31.8 %, PI 1.27 % and PS 13.35 %; severe OA group: PC 26.3 %, PI 0.59 % and PS 11.8 %).Additionally, compared to healthy SF-EVs, ether-linked phosphatidylcholine (PC O-) and PA classes demonstrated a relative rise in mild OA SF-EVs (healthy: PC O-1.71 % and PA 1.02 % mild OA: PC O-1.90 % and PA 1.34 %).However, the levels declined in severe OA-derived SF-EVs even more than the baseline levels in healthy joint derived EVs (severe OA: PC O-1.19 % and PA 0.55 %).Inversely, both PE types (ester-linked and etherlinked) showed a reduction in the mild OA-derived EVs compared to the healthy joint derived EVs (healthy: PE 8.06 % and PE O-11.8 %; mild OA: PE 5.72 % and PE O-8.17 %), while there was an increment in EVs isolated from the severe OA group (Severe OA: PE 10.2 % and PE O-10.2 %) with the ester-linked PE class level even higher than in EVs derived from healthy joints.

Fig. 3 .
Fig. 3. Proteomic profile of equine SF-EVs derived from healthy joints and from joints with mild and severe OA A) Unsupervised multivariate analysis using principal component analysis.The first two principal components were plotted, accounting for ~36.4 % of the variance.SF-EV samples were plotted based on acquired SWATH-MS data, after PQN normalisation and log transformation.Each plotted point represents a pooled SF-EV sample comprised of three biological replicates, which are colour-coded by OA severity, with severe OA in red (n = 3), mild OA in orange (n = 4), and healthy in green (n = 7).B) Heatmap demonstrating average protein intensities between SF-EV healthy (green (n = 7)), mild OA (orange (n = 4)) and severe OA (red (n = 3)) phenotypes.Protein intensities were transformed and are displayed as colours ranging from red to blue.Both rows and columns are clustered using the Ward method, and distance was calculated using Pearson Distance.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 4 .
Fig. 4. Ingenuity Pathway Analysis networks providing an overview of related molecular mechanisms.Analysis was conducted using the ingenuity knowledge base library and accounting for protein p-value and log2 fold following ANOVA and Tukey's post hoc test analysis.Node colour indicates up-regulated genes (from light pink: low upregulation to red: high upregulation), and solid lines represent direct interactions between pathways according to the Ingenuity knowledge base information.A) Mild OA compared to healthy (16 significant proteins), B) Severe OA compared to healthy (20 significant proteins), C) Severe OA compared with mild OA (9 proteins).(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 5 .
Fig. 5. Unsupervised proteomic and lipidomic data integration.Proteomic and lipidomic datasets from SF-EVs derived from healthy joints, mild OA and severe OA were normalised by the sum.A) Sparse Partial Least Squares-2 regression (sPLS2) of SF-EV samples projected into the area covered by the averaged components of both datasets.Healthy SF-EVs (green triangle (n = 6)), mild OA SF-EVs (orange cross (n = 4)), and severe OA SF-EVs (Sev.OA; red circle (n = 3)).B) Unsupervised multivariate sPLS2 arrow plot from the integration of proteomic and lipidomic data.The base of the arrow shows where a specific sample is in relation to the components of the proteomics dataset, and the tip of the arrow shows where the same sample is located concerning the components of the phospholipidomics dataset.Healthy SF-EVs (green circle), mild OA SF-EVs (orange circle), and severe OA SF-EVs (Sev.OA; red circle).The boxes zoom in on certain samples to better show the arrow direction C) Clustered Image Map from the sPLS2 data integration performed on the SF-EV omic datasets.The graphic shows the degree of similarity between the proteomic and lipidomic variables clustered over two dimensions and grouped using the Euclidean distance approach.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 6 .
Fig. 6.Network representation derived from the sPLS2 analysis of the proteomics and lipidomics integrated data.A relevance network plot with a correlation cutoff of 0.7 was created.Hence, only the variables with a correlation above 0.7 or below − 0.7 are shown.The networks are bipartite, and each edge connects a protein (rectangle) to a phospholipid (circle) node based on a similarity matrix.The colour of the connecting lines represents positive (red) or negative (green) correlations.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 7 .
Fig. 7. Discovery of potential composite biomarkers for OA.Phospholipids and proteins were normalised by the sum of the total amount of material (i.e., lipid or protein).Healthy (green (n = 6)), mild OA (orange (n = 4)) and severe OA (red (n = 3)).*p < 0.05, Kruskal-Wallis with Dunn's post hoc test A) Phospholipids and proteins that decreased with OA severity, B) phospholipids and proteins that increased with OA severity.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Table 1
Top 25differentially expressed (p < 0.05) proteins across SF-EV samples derived from healthy joints and joints with mild OA and severe OA following analysis of variance (ANOVA) and Tukey's post hoc test analysis, identifying significant experimental group comparisons and heatmap analysis.
E.Clarke et al.

Table 2
Top 5 canonical pathways identified using Ingenuity Pathway Analysis, following input of proteins correlated to lipids.