Skip to main content
Download PDF

Pulmonary artery embolism: comprehensive transcriptomic analysis in understanding the pathogenic mechanisms of the disease

BMC Genomics volume 24, Article number: 10 (2023) Cite this article

Abstract

Background

Pulmonary embolism (PE) is a severe disease that usually originates from deep vein thrombosis (DVT) of the lower extremities. This study set out to investigate the changes in the transcriptome of the pulmonary artery (PA) in the course of the PE in the porcine model.

Methods

The study was performed on 11 male pigs: a thrombus was formed in each right femoral vein in six animals, and then was released to induce PE, the remaining five animals served as a control group. In the experimental animals total RNA was isolated from the PA where the blood clot lodged, and in the control group, from the corresponding PA segments. High-throughput RNA sequencing was used to analyse the global changes in the transcriptome of PA with induced PE (PA-E).

Results

Applied multistep bioinformatics revealed 473 differentially expressed genes (DEGs): 198 upregulated and 275 downregulated. Functional Gene Ontology annotated 347 DEGs into 27 biological processes, 324 to the 11 cellular components and 346 to the 2 molecular functions categories. In the signaling pathway analysis, KEGG ‘protein processing in endoplasmic reticulum’ was identified for the mRNAs modulated during PE. The same KEGG pathway was also exposed by 8 differentially alternative splicing genes. Within single nucleotide variants, the 61 allele-specific expression variants were localised in the vicinity of the genes that belong to the cellular components of the ‘endoplasmic reticulum’. The discovered allele-specific genes were also classified as signatures of the cardiovascular system.

Conclusions

The findings of this research provide the first thorough investigation of the changes in the gene expression profile of PA affected by an embolus. Evidence from this study suggests that the disturbed homeostasis in the biosynthesis of proteins in the endoplasmic reticulum plays a major role in the pathogenesis of PE.

Peer Review reports

Background

Pulmonary embolism (PE) is a life-threatening disease that is one of the major causes of sudden unexpected deaths [1, 2] Most PEs originate from deep vein thrombosis (DVT), a condition occurring when a blood clot forms in the veins, typically of lower extremities (calf, femoropopliteal, and less frequently, the iliac veins) [3]. The thrombus consists mainly of fibrin, collagen, platelets and erythrocytes, and it reacts directly with the wall of the venous vessel, triggering various inflammatory processes [4]. This process is not unilateral as the endothelium exerts a feedback effect on clot formation and rearrangement [5, 6]. DVT characteristically arises in regions of reduced blood flow and then propagates as a result of hypercoagulability initiated by hypoxia and increased concentration of blood cells [1]. PE begins when the embolus disconnects from the point of origin and lodges in the pulmonary artery (PA), causing various severe consequences, including hemodynamic collapse and death [7]. When the embolic material moves into the PA, it enters a different environment. The question that arises is whether the thrombus formed in the vein will cause additional reactions of the arterial endothelium after reaching the PA, exemplary changes similar to those observed in acute coronary syndrome or ischemic stroke. When DVT and PE affect the same organism, they together constitute the disease known as venous thromboembolism (VTE) [2]. The prognosis in PE differs depending mostly on the degree of vessel obstruction, and it may lead to nonfatal VTE, fatal PE, postthrombotic syndrome, or chronic thromboembolic pulmonary hypertension [8]. However, the prognosis of PE/VTE patients is much worse than that of patients with only DVT. Treatment for both DVT and PE is similar and is targeted at preventing further increase of the clot size and new clots from forming, restoring the perfusion of the plugged vessel [9, 10]. To achieve this goal, anticoagulant therapy is widely used, initially with unfractionated or low molecular weight heparin, then, for at least 3 months, oral anticoagulants as vitamin K antagonists, or currently more preferred rivaroxaban, dabigatran, apixaban or endoxaban are used [9, 11]. In some cases of PE, fibrinolytic therapy can be applied in patients with clinical symptoms of haemodynamic instability [12]. The exact understanding of PE pathophysiology, including the knowledge of an intricate interrelationship between the PA endothelium and the lodged blood clot, helps in risk evaluation and allows for determining the most relevant and effective treatment. In the rat model, IL‑1β, TF, XOD and NF‑κB were found to be significantly differentially expressed, and thus were considered DVT biomarkers [13]. In the rabbit, it has been discovered, that genes related to the inflammation process or immune, pulmonary and cardiovascular diseases, such as IL-8, TNF-α, and CXCL5, may have functional importance for PE [14]. Due to multi-omics studies of the transcriptome and metabolome data of PE patients, LIPC and NAT2 were found to be differentially expressed in PE [15]. Further studies investigating the molecular background of thrombus formation and its effect on the venous and arterial endothelium in course of the DVT and PE are needed. Therefore, in our previous study, we already investigated the transcriptomic profile of femoral veins in DVT, and we discovered that dysregulation of TNF, NF-κB, and apoptosis pathways affect inflammatory responses and the formation, remodeling, and resolution of thrombus in DVT pathogenesis [16]. Whereas, the present research was aimed at analyzing the changes in the gene expression associated with PE in the PA samples (PA-E). To accomplish this goal the experiment was performed in a new pig model, in which the thrombus was first formed in vivo in the femoral vein, and then was released to induce PE in the same animal, mimicking the natural way of its formation, as well as its passage through the vascular bed. As the PA-E samples were collected 3 h after the thrombus release, the changes observed in the investigated tissues most closely reflect these typical for acute PE.

Methods

Surgical procedures in the experimental and control animals

The study was performed in 11 castrated male pigs (24 weeks old, 60 kg body weight, b.w.) of the Polish Landrace breed; 6 animals served as the experimental group (further referred to as pulmonary arteries affected by embolism PA-E); the remaining 5 pigs served as controls (CTR). Several different experimental models of PE involving rats, mice, pigs and rabbits have been described [17,18,19,20,21,22,23]. In the present study experimental model of VTE was based on the one proposed by Robinson et al. [24] as described previously. In the experimental group first, the DVT was developed by forming the thrombus in the right femoral veins using the porcine model described previously [25]. Briefly, before performing any surgical procedures, all the DVT pigs were pretreated with atropine (Atropinum Sulfuricum, Polfa, Warsaw, Poland, 0.05 mg/kg b.w., s.c.) and azaperone (Stresnil, Janssen Pharmaceutica, Beerse, Belgium, 2.5 mg/kg b.w., i.m.). Thirty minutes later, to induce anesthesia, the main anesthetic drug propofol (Propofol-Lipuro, B. Braun Melsungen AG, Melsungen, Germany, 10 mg/kg b.w.) and the main analgesic drug ketamine (Bioketan, Vetoquinol, Poland, 10 mg/kg b.w.) were given intravenously in a slow, fractionated infusion. The depth of anesthesia was monitored by testing the corneal reflex. Once the animals were transported to an operating theater, general anesthesia during the entire procedure was maintained with inhalation of sevoflurane (Sevoflurane, Baxter, Ontario, CA USA; administered at one human MAC end-tidal concentration of 2.0%) under continuous pulse oximetry and heart-rate monitoring. Under sterile conditions, the left femoral veins were carefully exposed in each animal by dissecting the fascia and exposing and preparing the sartorius muscle. Next, the vein was closed proximally and distally at the length of about 30–40 mm with surgical ligatures. To accelerate the thrombus formation, 200 units of thrombin (BioTrombina, Biomed Lublin S.A., Lublin, Poland; the doses of the treatment agents were corresponding to those previously used in other experiments; [24]) were administered into the closed segment of the vein. Immediately afterwards, to prevent the formation of thrombi outside the closed segments, a bolus of unfractionated heparin (Heparinum WZF, Warsaw, Poland) was administered, in a dose of 100 U/kg, through a catheter inserted in the ear vein. Two hours later (an optimal time needed for thrombus formation after the closing of the femoral vein), the thrombus was released from the right femoral vein to induce PE. The thrombus release was facilitated mechanically by bending the right hind limb (thus the movements of muscles were acting as natural pistons). If this procedure was ineffective, 10 mL of saline was injected into the right femoral vein distally from the formed thrombus. Computed tomography pulmonary angiography (CTPA; Aquilion PRIME TSX-303A Toshiba, Tustin, CA, USA), a method of choice serving to assess the pulmonary vessels when PE is suspected, was performed in all experimental pigs after the release of the thrombi from the right femoral veins (Fig. 1).

Fig. 1
figure 1

A representative images of: (A) an embolus collected from the segmental pulmonary artery (B) an embolus, indicated by a white arrow, in the segmental pulmonary artery found at autopsy

Full size image

Using a maximum intensity projection (MIP) technique the thrombi were visualized and measured with the Toshiba Aquilion PRIME TSX-303A as a standard software tool. During the CTPA procedure, the animals were sedated with propofol (Propofol-Lipuro, B. Braun Melsungen AG, 10 mg/kg b.w.). Three hours after the CTPA, all experimental pigs were euthanized with sodium pentobarbital (Euthasol, FATRO, Ozzano dell’Emilia BO, Italy, 140 mg/kg). The entire heart and both lungs were dissected from each animal. After the incision of the right ventricle and the right atrium, the pulmonary trunk was incised immediately above the pulmonary valve. Then, very carefully and delicately, PAs were successively cut open, beginning from the pulmonary trunk, and moving peripherally towards segmental and subsegmental arteries, to determine the size and location of the thrombus released from the vein. After removal of the thrombus with tweezers, the part of the PA where the blood clot was previously located was collected for transcriptome analyses. The control animals were first pretreated with atropine and azaperone, and thirty minutes later, propofol and ketamine were given intravenously to induce anesthesia; next, the control pigs were euthanized with sodium pentobarbital (the drugs, their doses and routes of administration were analogous to those used in the experimental group). The entire heart and both lungs were dissected from each animal. Then, very carefully and delicately, PA was exposed in each control animal and cut out for transcriptome analysis (the samples of the PA were collected in a way that exactly corresponded with that applied in the experimental animals, i.e., both the place of sampling and the length of the artery segment were the same). All the animals were kept under standard laboratory conditions. They were fed standard fodder (Grower Plus, Wipasz, Wadąg, Poland) and had free access to water. The animals were housed and treated according to the guidelines of the local Ethics Committee for Animal Experimentation in Olsztyn (affiliated with the National Ethics Committee for Animal Experimentation, Polish Ministry of Science and Higher Education; decision No. 90/2018 of 13.02.2018). All procedures conformed to the guidelines from Directive 2010/63/EU of the European Parliament on the protection of animals used for scientific purposes or the NIH Guide for the Care and Use of Laboratory Animals.

RNA extraction, library construction and high-throughput sequencing

Total RNA was isolated from the collected parts of the PA of both the experimental and control pigs with the use of the mirVana kit following the manufacturer’s procedure (Thermo Fischer Scientific, Waltham, MA., USA). Extracted RNA was subjected to the quantity and the quality evaluation by the Bioanalyzer 2100 (Agilent Technologies, Waldbronn, Germany). Further, the aliquots with the highest RIN values (RIN > 8) and concentrations (at least 50 ug/ul) were selected for the next steps of sequencing procedures. The sequencing procedure was held by an outsourcing company (Macrogen, South Korea) exploiting Illumina NovaSeq 6000 System (Illumina, San Diego, CA, USA). Briefly, 1 ug of total RNA for each sample was selected for library construction by the Illumina TruSeq mRNA LT Sample Prep kit (Illumina, Inc., San Diego, CA, USA). The laboratory procedure assumed the following steps: the purification of mRNA molecules using poly-T-attached magnetic beads, cutting mRNA into small fragments with divalent cations, amplification of the first-strand cDNA using SuperScript II reverse transcriptase (Invitrogen, Waltham, MA, USA) with random primers and finally the second strand cDNA synthesis with DNA Polymerase I and RNase H. The qPCR Quantification Protocol Guide (KAPA Library Quantification kits for Illumina Sequencing platforms) and TapeStation D1000 ScreenTape system (Agilent Technologies, Waldbronn, Germany) were used to quantify and qualify the RNA-seq libraries. Next, the prepared libraries were transferred to the NovaSeq 6000 platform (Illumina, San Diego, CA, USA), and the sequencing procedure was run with the expected configuration: 150 bp sequence length of paired-end reads and assumed sequencing depth of 40 M reads per sample.

In silico profiling of pulmonary arteries transcriptome affected by embolism

The raw high-throughput sequencing dataset retrieved from NovaSeq 6000 was evaluated according to the quality control standards using FASTQC software version 0.11.7 [26]. The paired-end reads (2 × 150 bp, type stranded) were trimmed after the Illumina adaptors identification within the sequences and also low-quality reads (PHRED cut-off score < 20) were excluded from downstream analysis using Trimmomatic software v. 0.38 [27]. The 120 bp trimmed reads were mapped to the pig reference genome according to ENSEMBL annotation (Sus_scrofa.Sscrofa11.1.99). The mapping results were saved and sorted by coordinate in Binary Alignment Map (BAM) files. Sequences aligned multiple times were not considered for subsequent analysis. The StringTie v. 2.2.1 pipeline [28] was incorporated to re-evaluate the ENSEMBL annotations and BAM files to obtain the novel annotations of intergenic-expressed regions. The integrity of the RNA-seq libraries was validated and clustered with the ggplot2 Bioconductor package v. 3.3.5 of R software v. 4.1.3 (R Core Team; www.R-project.org). Whole-transcriptome high-throughput sequencing (RNA-seq) of PA-E libraries was applied to identify expression profiles of differentially expressed genes (DEGs), other noncoding transcripts, the alternative splicing events and RNA nucleotide polymorphisms (Fig. 2).

Fig. 2
figure 2

Target enrichment workflow of the experiment assigning step-by-step molecular and bioinformatic procedures. ASEs – allele-specific expression variants; DAS – differential alternative splicing; DEGs – differentially expressed genes; DELs – differentially expressed lncRNAs; lncRNAs – long noncoding RNAs; SNVs – single nucleotide variants

Full size image

Differentially expressed genes and differentially expressed long non-coding RNA connected with the pulmonary embolism

The differentially expressed analyses were performed for the protein coding and long noncoding RNA (lncRNA) transcripts. The expressed transcripts were grouped according to the genomic localization and tagged as transcriptomic active regions (TARs). The differentially expressed analysis was performed by the binominal statistical test implemented in the DESeq2 method [29]. The TARs localizated in the range of protein coding genes were classified and signed as DEGs if their expression fluctuation between PA-E and CTR was statistically significant (adjusted p-value < 0.05). Moreover, lncRNAs with changed expression profiles between compared groups were classified as DELs. Differential expression analysis was conducted on the genomic level. Only long noncoding RNA sequence identification was operated on the transcriptomic level. The research assumed the prediction of all differentially expressed long non-coding RNA (DELs) shinned during the course of PE. The process of lncRNAs identification assumed the prediction of long (> 200 bp) transcripts without coding potential. Noncoding genes with the multi-exonic structures without Rfam and Pfam annotations and characterized by a lack of coding potential, confirmed by three methods—Coding-Potential Assessment Tool (CPAT) [30], Coding Potential Calculator (CPC2) [31] and Flexible Extraction of LncRNAs (FEELnc) [32], were named as lncRNAs. The count gene expression matrix was created from StringTie [28] and Ballgown [33] tools compilation supported by prepDE Python script provided by the StringTie’s authors (https://github.com/gpertea/stringtie/blob/master/prepDE.py; accessed on 21 January 2022). The candidate DEGs and DELs were processed in the MA, volcano, heatmap plots with ggplot2 [34] and circlize R packages v. 0.4.14. Trans–regulatory elements were investigated to uncover the relationship between lncRNAs and protein-coding genes. Based on correlations (Pearson’s coefficient—r) of lncRNA-mRNA expression profiles, the regulatory connections were statistically confirmed (|r|> 0.9) using Hmisc v. 4.7.0 package R software.

Differential alternative splicing genes involved in pulmonary embolism

Alternative splice junction positions of genes were recovered by the rMATS [35]. Briefly, aligned paired-end reads were evaluated according to the difference in transcript alternative splicing usage profile. The splicing patterns of genes were determined by means of reads arranged around splicing junction sites. The percentage of splicing inclusions (PSI) value was used as a basic unit for splicing event descriptions. Differential alternative splicing genes (DASGs) were statistically tested with the following sensitive parameters: adjusted p-value < 0.01 and |ΔPSI|> 0.1. All discovered differential alternative splicing (DAS) events were grouped into five categories: alternative 5’ splice site (A5SS), alternative 3’ splice site (A3SS), mutually exclusive exons (MXE), retained intron (RI) and skipping exon (SE). Significant DAS events were visualized using the maser R package v. 1.12.1 [36, 37] and ggsashimi Python tool v. 0.5.0 [38].

Single nucleotide variants manifest allele frequency bias

SNVs calling results were achieved by compilation of GATK [39], rMATS-DVR [40] and custom R scripts. After recalibration of BAM files by Picard tool v. 2.6.0 (https://broadinstitute.github.io/picard/; accessed on 13 September 2021), the SNVs were recognized and filtered out according to standard GATK parameters. The applied threshold to filter the SNVs was established on minimum coverage (10 reads). The Variant Call Format (VCF) file was created and transferred to the R environment for the subsequent stages of the filtration procedure. The polymorphic sites were positioned on chromosomes and SNVs ranged within bidirectional transcription regions, intronic flanking regions, simple sequence repeat regions, and pseudogenes regions were eliminated from the further statistical comparison of allele frequency. To elucidate the above filtration procedure the GMATo v. 1.2 [41], BLAT v. 36 [42], BEDTools v. 2.30 [43] tools and vcfR v. 1.12 R package [44] were adapted. The significant allele frequency differentiation (|ΔAAF|> 0.3; FDR < 0.001) was retrieved by counting and comparing reference and alternative nucleotides using Samtools’ mpileup function v. 1.6 [45]. The obtained significant results were divided into two groups: localised within genic regions or signed as single nucleotide polymorphism (SNPs) and intergenic or other noncoding SNVs. Both pools of SNVs were tested by VEP [46] to predict the gene structure annotation and effects of polymorphic variants. Based on the smallest standard deviation between biological replicates, the final set of nucleotide substitutions was retrieved and signed as allele-specific expression (ASE). The circos software package [47] was implemented to visualise all variants.

Functional annotations (gene ontology and KEGG enrichment) of differentially expessed genes, differentially alternative splicing genes and allele-specific expressions

Obtained DEGs, DASGs and ASEs were scanned by gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) annotation using g:Profiler software [48]. The essential genes were annotated to ontological terms within three aspects, such as biological processes (BP), cellular components (CC), molecular functions (MF) and also assigned to signalling or metabolic pathways (KEGG database) [49]. The enrichment analysis (Fisher's exact test with Benjamini and Hochberg correction; FDR < 0.05) was applied to uncover ontology and pathway annotations regulated by expressed genes and post-transcription mechanisms, including spliceosome and SNV modifications. To obtain the cardiovascular gene signatures, the DEGs were also scanned according to the Human Phenotype Ontology (HP) database according to the terms representing individual phenotypic anomalies (https://hpo.jax.org/app/; accessed on 20 September 2021). To visualise obtained molecular network interaction, a Pathview v. 1.30.1 R package [36] was applied.

Real-Time PCR validation

The mRNA level of selected genes was determined by Real-Time PCR. The primers for chosen genes were designed using the Primer3Plus software [50] based on the sequences listed in Table S1. The cDNA was obtained with the use of the Applied Biosystems High-Capacity cDNA Reverse Transcription Kit (Thermo Fisher Scientific, Vilnius, Lithuania) according to the manufacturer’s protocol. The Real-Time PCR was performed with the use of Applied Biosystems PowerUp SYBR Green Master Mix (Thermo Fisher Scientific, Vilnius, Lithuania) according to the manufacturer’s protocol on the QuantStudio 3 Real-Time PCR System (Applied Biosystems, Thermo Fisher Scientific Inc., Waltham, MA, USA). In brief, each reaction contained 5 μL of master mix (2X), forward and reverse primers in the amount of 1000 nM of each, 10 ng of cDNA, and a proper volume of nuclease-free water to a final volume of 10 μL. The reactions were performed in four technical replicates for each biological sample. The expression of each gene was calculated using the comparative Pfaffl method [51], where the expression is presented as the fold change relative to the control, as well as normalized to an endogenous reference gene (ACTB, GenBank U07786.1, and GAPDH, GenBank AF017079) (relative quantification RQ = 1). The results were expressed as means of biological replicates ± standard deviations. Statistical analysis was performed using Student's t-test (two-tailed) in Prism 8 software (GraphPad Software Inc., San Diego, CA, USA). P-values were considered statistically significant where 0.0332 (*), 0.0021 (**), 0.0002 (***), and < 0.0001 (****).

Results

Pulmonary arteries transcriptomic statistics and the abundance of expression profiles

After sequencing 395 647 360 raw reads were obtained, globally (Table 1). The filtration procedure removed 32 564 994 reads with a poor quality score (PHRED33 < 20). The surviving 363 082 366 paired-end reads were mapped to the porcine reference genome. Results of the mapping process were applied to the identification of the DEGs and DAS events. Uniquely mapped reads contained an average of 96% out of all processed sequences. According to the gene structure, 46.04% of paired-end reads’ nucleotides were aligned to the coding sequence (CDS), 11.24% to the intronic sequences, 26.58% to the untranslated regions (UTR) and 16.14% to the intergenic localizations. The transcriptomic data from this study have been submitted to the European Nucleotide Archive under accession No. PRJEB43020 (https://www.ebi.ac.uk/ena/browser/view/PRJEB43020).

Table 1 The statistical metrics for the RNA libraries. The sequencing and mapping results for the 11 RNA-seq libraries: CTR (1–5) refers to the controls; PA-E (1–6) refers to the pulmonary arteries affected by embolism. The uniquely-mapped values were related to the reads mapped to only one location of the porcine genome. The multi-mapped values were assigned to reads aligned to more than one locus on the reference genome. The raw, trimmed and mapped values denote in millions
Full size table

The principal component analysis (PCA) revealed high homogeneity of expression profiles within the PA-E group. CTR biological replicates indicated less homogeneity (Fig. 3A). Among 33 270 TARs, 3 330 transcripts were highly expressed (FPKM > 100), according to the FPKM sum of all RNA libraries. The 11 923 TARs indicated medium expression (FPKM > 10) and the rest of 18 017 TARs revealed low expression values in PA-E libraries.

Fig. 3
figure 3

Expression profiles overview of pulmonary embolism. A Graphical representation of the first (PC1) and second (PC2) principal components distributed the sample expression pattern of pulmonary arteries affected by embolism (PA-E; n = 6) and control (CTR; n = 5) libraries. B MA chart with logarithmic values of fold change (logFC; Y axis) plotted against normalised counts (X axis) for PA-E and CTR libraries. C Volcano plot shows log2FC plotted against negative logarithmic adjusted p-values. The dotted horizontal line indicates a negative logarithmic adjusted p-value (0.05) cut-off. (B, C) Red dots represent upregulated differentially expressed transcribed active regions (TARs); green dots depict downregulated TARs; black dots are not significant transcripts, according to DESeq2 method

Full size image

Genes and long noncoding RNAs changed their expression profiles during pulmonary embolism

Screening RNA-seq data for differential gene expression analysis revealed that the PE was associated with 604 differentially expressed TARs (Fig. 3B,C). Among them, 473 TARs were classified as DEGs, which encoded protein sequences (Fig. 4; Table S2). The depth transcriptome analysis also revealed 46 DELs (Fig. 4; Table S3). Within DELs, 4 were annotated in the ENSEMBL database and the other 42 were newly-identified in this research. The trans-correlation of 32 DELs mediated regulation of 57 DEGs. Most trans gene-lncRNA interactions (101) were positively correlated, though only 4 displayed negative expression co-action.

Fig. 4
figure 4

Circular heatmap visualisation of differentially expressed genes (DEGs) and long non-coding RNAs (DELs) in the course of the pulmonary embolism. The 11 upper tracks visualise the normalised (Z-score; red-green scale) expressions for DEGs in each biological replicate (PA-E—pulmonary arteries affected by embolism and CTR-control libraries). The large part of the circle depicts DEGs and the smaller part describes DELs. The inner track shows the correlation links between the co-expressed DEGs and DELs, whereas blue links depict positive (> 0.9) and yellow negative (< -0.9) Euclidean correlation

Full size image

Alternatively spliced transcripts engaged in pulmonary embolism

The reads mapped in the range of acceptor–donor sites were evaluated according to the alternative splicing variants. The splicing detection algorithm (rMATS) allowed the identification of 136 666 alternative splicing events. The 348 splicing events pointed out significant changes in the percentage of exon/intron inclusions in the PA-E vs. CTR comparison. Significant events were incorporated in 229 DASGs. Within detected DASGs, 13 were also classified as DEGs. Among all of the detected DASGs, 43 were classified as A5SS, 56 as A3SS, 27 as MXE, 65 as RI, and 157 as SE (Table S4). Sashimi plots were prepared to visualise exon inclusion levels of the TRAF2 and UBQLN4 (Fig. 5).

Fig. 5
figure 5

Quantitative visualisation (Sashimi plot) of alternatively splicing events statistically significant in changes of percentage splicing inclusion (adjusted p-value < 0.01 and |ΔPSI|> 0.1) between PA-E (pulmonary arteries affected by embolism; red) and CTR (control; violet) samples. The numbers in each upper right corner of the violet and red tracks present the percent splicing inclusion (PSI) values. Splicing differentiations are covered by the number of reads mapped in the range of junction sites (A) TRAF2 (A5SS—alternative 5’ splice site), (B) UBQLN4 (SE – skipping exon). The bottom black tracks show the genomic localisation of splicing events

Full size image

Single nucleotide variants callings with a difference in alternative allele fraction as the signature of pulmonary embolism

In silico investigation revealed the existence of the 456 993 SNVs. To elucidate the SNVs with allele frequency bias between PA-E and CTR samples, the filtration process qualified 70 349 nucleotide variants. The 11 115 ASEs variants revealed statistically significant allelic imbalance (|ΔAAF|> 0.3 and FDR < 0.001) and within them 9 595 were assigned as known porcine SNP candidates (with ENSEMBL ‘rs’ ids). Among SNVs, 68.64% of alternative alleles indicated the highest frequencies in the PA-E libraries. The remaining SNPs were annotated giving 9 215 substitutions localised within the protein-coding genes. In accordance with the VEP annotations, 443 SNVs were assigned to lncRNA regions. Identified SNVs were characterised according to all transcripts main annotation regions: upstream (5.78%), 5’ UTR (2.24%), CDS: missense mutations (4.93%) and synonymous mutations (15.71%), introns (29.25%), 3’ UTR (22.48%) and downstream (19.62%). According to the standard deviation measure between biological replicates, the 536 SNVs were selected for the final set of ASEs. Genes localization consequences of the substitution with allele-specific bias within the transcripts components were summarised in Fig. 6. In the range of ASEs, 792 annotations were described by VEP (Table S5). The 61 ASEs were functionally assigned to ‘endoplasmic reticulum’.

Fig. 6
figure 6

Circos visualisation of single nucleotide variants (SNV) candidates indicated allelic bias between pulmonary arteries affected by embolism (PA-E) and control (CTR) libraries. The external circle illustrates porcine chromosomes and other unclassified scaffolds, where region length is proportional to the number of alternative allele substitutions. The genes belonging to the cellular components of the endoplasmic reticulum were marked on the outer track. The second track depicts canonical (A-G, C-T; black) and other types of nucleotides substitutions (grey). The third track (histogram) depicts differences in the alternative allele fraction (ΔAAF) between PA-E and the CTR libraries. In the eight middle scatter plots, RNA editing substitutions are shown as upstream (blue triangles), 5' UTR (red rectangle), missense (yellow circle), synonymous (purple circle), intron (green squares), splice acceptor (grey circle), 3' UTR (orange triangle) and downstream (grey rectangle). Each inner track shows the coverage of alternative (red histogram) and reference variants (green histogram) in all RNA-seq libraries

Full size image

Gene ontology classification of the transcripts expressed in pulmonary arteries affected by embolism

Active genes involved in PE were classified by GO and 347 of them were assigned to the 27 BP, 324 to the 11 CC and 346 to the 2 MF terms annotations. According to the GO annotation of these genes, the BP enrichment was found for example in the ‘negative regulation of response to stimulus’ and ‘negative regulation of epithelial cell proliferation’ (Figs. 7A and 8).

Fig. 7
figure 7

Selected members of the functional and pathway enrichment annotations. A Circos plot represents the three most important terms of biological processes (BP) and KEGG pathway associated with differentially expressed genes (DEGs) engaged in pulmonary embolism (PE). Gene symbols with logarithmic values (blue-red scale) of fold change (logFC) are located on the left side of the circos. B The outer circle shows the scatter plot for cellular components (CC) and metabolic function (MF) ontology terms enriched with DEGs. Red points represent up-regulated genes and blue points symbolise down-regulated genes. The inner circle contains information on the abundance of genes assigned to the Gene Ontology (GO) term and the overall impact (increasing or decreasing) trend on this functional term. Z-score is calculated using the summarized expression DEGs’ values within particular GO terms

Full size image
Fig. 8
figure 8

UpSet plot merging most important GO terms and associated DEGs. The heights of the bars indicate the number of DEGs annotated to single or multiple GO terms. Joined dots depict connections between particular GO terms and intersected common DEGs

Full size image

In the MF category, the DEGs were involved in only two processes of the general meaning, like ‘binding’ and ‘binding protein’ (Fig. 7B). While the CC category grouping revealed that PA-E mRNAs were involved for instance in ‘endoplasmic reticulum’, ‘extracellular space’ and ‘nucleoplasm’ (Fig. 7B). The comprehensive GO enrichment classification was summarised in Table S6.

Twelve genes (DAD1, DDOST, UBQLN4, PDIA6, DNAJB11, CALR, LMAN2, PDIA4, NFE2L2, MARCHF6, RRBP1, TXNDC5) belonging to ‘protein processing in endoplasmic reticulum pathway’ were also associated with other GO terms. Eight of those genes were associated directly with the ‘endomembrane system’. Other genes like CALR, NFE2L2, UBQLN4 were connected with multiple GO terms, namely: ‘system development’; ‘endomembrane system’; ‘negative regulation of cellular metabolic processes’; ‘negative regulation of nitrogen compound metabolic process’, negative regulation of response to stimulus’ and negative regulation of signal transduction’. Exploring the HP and GO database, DEGs were also grouped according to the potential connection with the development and abnormality of the cardiovascular system. The following DEGs: AMHR2, PTPRB, EPHB2, THBS2, RORA, PTGIS, TMEM204, XDH, EMILIN1, SFRP2, CREB3L1, COL1A2, ACKR3, GLI3, HOXB3, SERPINF1, PIK3C2A, HIPK2, PTPRM, PRICKLE1, COL1A1, GADD45A, FOXC1 were classified as signatures of the ‘cardiovascular system development’ in the course of PE. ‘Abnormality of the cardiovascular system’ was represented by the following DEGs: AAAS, C4A, EPHB2, GREB1L, USP45, LIPC, BNC2, COQ4, RBCK1, PTGIS, WIPI2, GUSB, TTC7A, GUCY1A1, SMARCB1, ACSL4, SERPING1, WT1, ABCC8, CALR, CRLF1, CEP120, CEP164, COL1A2, BIN1, NFE2L2, IMPDH1, GLI3, LIFR, PSMD12, MEOX1, COL6A1, RAD51C, EVC2, PIBF1, CDK4, PIK3C2A, GGCX, COL6A2, SELENOI, IMPG2, ZNF469, PTCH2, EVC, PYCR1, GDF6, COL1A1, FANCF, FKBP14, NSMCE2, TRIM28, SNAI2, FOXC1, CDKN1C, SKI, which may be crucial as genetic biomarkers of the PE etiology. In the signaling pathway analysis, one KEGG pathway was identified for the mRNAs modulated during PE. Within twelve genes enriched in KEGG ‘protein processing in endoplasmic reticulum’ pathway, ten were underexpressed (DAD1, DDOST, UBQLN4, PDIA6, DNAJB11, CALR, LMAN2, PDIA4, RRBP1, TXNDC5) and two others were overexpressed (NFE2L2, MARCHF6) (Figs. 7A, 8 and 9). The same KEGG pathway (ssc04141) was also exposed by DASGs (EIF2AK2, UBQLN4, BAG2, SEC31A, EDEM2, EIF2AK4, TRAF2, MAPK8), with alternative splicing bias between PA-E and CTR samples (Fig. 9; Table S7). Within SNVs with significant allelic imbalance, the 61 polymorphisms were localized in the vicinity of the following genes: ITGA8, CCDC3, ITPR1, BNIP3L, ACSL5, ZDHHC6, SEC23B, HM13, LARS1, PTDSS2, TMEM178A, TMEM214, NBAS, ARL6IP1, SULF1, CERS5, TMTC3, KDELR3, DMPK, TMEM50A, PML, PKD2, SEC31A, TMEM33, BCAP29, ACER3, THY1, that belong to the cellular components of ‘endoplasmic reticulum’ (Fig. 6).

Fig. 9
figure 9

Enrichment Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of differentially expressed (DE) genes and differentially alternatively splicing (AS) genes engaged in the ‘protein processing in endoplasmic reticulum’. Green and red rectangles present downregulated and upregulated genes, respectively. Blue and yellow rectangles reflect DAS genes with increased and decreased PSI values in pulmonary arteries affected by embolism (PA-E) samples. Logarithmic fold change (logFC) red-green scale describes gene expression values. The PSI purple-yellow scale (|ΔPSI|> 0.1) describes percentage inclusion values

Full size image

Comparison of DEGs detected within this research with DEGs from other studies focused on expression analysis in DVT and PE revealed that some of them were overlapping, regardless the species, model, or transcriptomic assay used (Fig. S1). Among 473 identified DEGs, 29 were present also in analyzed other data sets [14, 52,53,54]. Most of the DEGs (IGFBP6, MOXD1, SLC4A11, TTYH3, ATOH8, PDIA6, SLK, VCAN, ME3, COL1A2, PSMD12, CNP, DACT2, ENTPD4, SHMT2) were overlapping between current project and data from Tang et al. [14] and Zagorski et al. [52]. There was one DEG uncovered in all of the analyzed datasets—PLOD2.

Gene expression quantification

The genes for validation were selected following the assessment of the function, expression values and read’s distribution within samples. The overexpression of LIFR and PIK3C2A genes (over 2 arbitrary units compared to the control) was assigned to the abnormalities of the cardiovascular system, making them crucial biomarkers of the PE etiology. The significant overexpression in PA-E compared to the CTR libraries was also noted in the ACSL4 and SELENOI engaged in vascular homeostasis and various membrane-related cellular processes. On the other hand, COL1A1, COL1A2 and COL6A1, which are also crucial for vascular remodelling, were significantly downregulated in PA-E libraries, the same as C4A involved in local inflammatory processes. The RNA-Seq results were validated using Real-Time PCR which revealed overexpression of six genes: ACSL4, ENSSSCT00000024889, LIFR, PIK3C2A, SELENOI and WDR37 whereas underexpression of another six genes: C4A, COL1A1, COL1A2, COL6A1, PTGIS and PYCR1 (Fig. S2).

Discussion

DVT and PE are common life-threatening conditions, therefore the understanding of their molecular background is of great importance, as the identification of genes and pathways is crucial in disease pathomechanism and may reveal potential treatment targets. The current research, which aimed at analyzing the changes in the transcriptome of PA with induced PE, is the continuation of our previous studies devoted to the transcriptomic profiling of femoral veins in DVT. We mainly disclosed that dysregulation of TNF, NF-κB, and apoptosis pathways affect inflammatory responses and the formation, remodeling, and resolution of thrombus in the course of DVT [16]. Whilst, the present project not only uncovered the genetic biomarkers of PE but also unrevealed that the disturbed homeostasis in the biosynthesis of proteins in the ER plays a major role in the pathogenesis of PE.

Genetic biomarkers of the pulmonary embolism (PE)

The present research revealed several DEGs involved in vascular wall integrity and strongly associated with dysfunction of the cardiovascular system, which may be potentially considered genetic biomarkers of PE. In particular, underexpression of COL1A1 (logFC = -2.39), COL1A2 (logFC = -2.14), COL6A2 (logFC = -1.17), COL6A2 (logFC = -1.17), C4A (logFC = -1.56), GLI3 (logFC = -0.89), and overexpression of ACSL4 (logFC = 1.27), PIK3C2A (logFC = 1.02) genes in the PA-E libraries is worth emphasizing. The extracellular matrix, composed of collagens, elastin, and proteoglycans, is an essential component of the PA wall [55]. Our study disclosed four underexpressed DEGs in the PA-E, which were associated with collagen synthesis: Collagen Type I Alpha 1 Chain (COL1A1), Collagen Type I Alpha 2 Chain (COL1A2), Collagen Type VI Alpha 1 Chain (COL6A1) and Collagen Type VI Alpha 2 Chain (COL6A2). COL1A1 and COL1A2 belong to group I collagen (COL1; fibrillar forming collagen), whereas COL6A1 and COL6A2 belong to group VI collagen (COL6), which acts as a cell-binding protein [56]. COL1 is the main type of collagen, present in all three layers of the vessel wall [55], and it not only provides structural support but also is crucial for the growth, homeostasis and repair of the tissue [57]. COL6, which is non-fibrillar collagen, is engaged in the vascular wall remodeling process. It has been found, that COL6 was significantly decreased in the tunica media of rat PA, following the hypoxia [58]. Thus, it can be concluded that the amount of COL6 has a strong impact on the proliferation and migration of vascular smooth muscle cells. Undoubtedly, the accumulation of the clot in the PA creates hypoxic conditions, which may be the reason for COL61 and COL62 underexpression. Intriguingly, Luther et al. [59] revealed that deletion of COL6 in knockout model mice plays a critical protective role following myocardial infarction by limiting infarct size, chronic apoptosis, aberrant remodeling, and fibrosis. Therefore, the underexpression of COL genes observed in the present experiment may potentially reflect protective mechanisms triggered to preserve the functional integrity of the vessel wall. Furthermore, Complement C4A (C4A), GLI Family Zinc Finger 3 (GLI3) and Acyl-CoA Synthetase Long Chain Family Member 4 (ACSL4) genes which were differentially expressed in the investigated samples, have also a strong impact on the vascular integrity. C4A participates in the complement classical activation pathway [60]. Derived from proteolytic degradation of complement C4, C4A anaphylatoxin is a mediator of local inflammatory processes [61]. It induces the contraction of smooth muscles and increases vascular permeability. An increasing body of evidence supports the functional role of complement activation in the pathogenesis of cardiovascular disease through pleiotropic effects on endothelial, hematopoietic cell function, and haemostasis [62]. Therefore, the observed underexpression of C4A may again reflect the protective processes induced to limit the local thromboinflammation. Possibly this effect may be enhanced by the underexpression of GLI3, which acts as a transcription factor that regulates vessel growth, endothelial-cell activity, and apoptosis [63]. It has been found that overexpression of GLI3 upregulates the expression of pro-angiogenic factors in endothelial cells, including interleukins and chemokines [63]. Thus, respectively, its downregulation may mute inflammation. Acyl-CoA Synthetase Long Chain Family Member 4 (ACSL4) expressed in arterial smooth muscle cells may regulate some processes dependent on the release of arachidonic acid mediators [64]. For instance, it has been found that the overexpression of ACSL4 results in increased synthesis of eicosanoids that play important roles in vascular homeostasis. Moreover, ACSL4 modulates prostaglandin E2 release from smooth muscle cells [64]. These observations suggest that the upregulation of ACSL4 disclosed in the PA-E samples, may, again reflect the local inflammation process and elucidates the usefulness and necessity of non-steroidal anti-inflammatory drug application in the PE treatment. Additionally, the present study revealed the overexpression of Phosphatidylinositol-4-Phosphate 3-Kinase Catalytic Subunit (PIK3C2A) in the PA-E. Phosphatidylinositol 3–kinases (PI3K) are recruited during the clot formation after binding collagen to GPVI on platelets. Activation of PI3K causes the release of secondary mediators, such as ADP and thromboxane A2, culminating in platelet aggregation mediated by fibrinogen [65]. Many aspects of cardiovascular physiology and disease have been linked to PI3K isoforms. The key roles of PI3Ks in platelets suggest that they could all represent valid antithrombotic targets, to inhibit thrombosis while having a limited effect on normal haemostasis [66].

Protein processing in the endoplasmic reticulum (ER) pathway

The endoplasmic reticulum (ER) is a large cellular organelle that serves important functions in the cell, including synthesis, folding, modification, and transport of proteins, production of lipids and sterols, as well as storage of free calcium [67]. Numerous cellular stresses can lead to disturbances between the demand for protein folding and the capacity of the ER to efficiently sustain this process, thereby initiating ER stress pathway [68]. This specific pathway is involved in a wide range of pathologies [69], for instance, it may play a crucial role in the development of cardiac hypertrophy and heart failure [70]. ER stress pathway is activated to protect cells, but when stresses are severe, apoptosis is induced to remove damaged units. Current in-depth bioinformatics analysis revealed that the ‘protein processing in endoplasmic reticulum’ pathway (ssc04141; Fig. 8) significantly contributes to the pathogenesis of PE. Our research revealed twelve genes enriched in this KEGG pathway, ten were underexpressed: DAD1 (logFC = -0.60), LMAN2 (logFC = -0.58), RRBP1 (logFC = -1.09), DDOST (logFC = -0.59), PDIA4 (logFC = -0.76), UBQLN4 (logFC = -0.70), DNAJB11 (logFC = -0.71), PDIA6 (logFC = -0.77), CALR (logFC = -0.90), TXNDC5 (logFC = -0.60), and two others were overexpressed: NFE2L2 (logFC = 0.69), MARCHF6 (logFC = 0.44). The obtained results indicate that genes within this pathway were not only differentially expressed but also exposed by DASGs (EDEM2, MAPK8, TRAF2, UBQLN4, EIF2AK2, BAG2, SEC31A, EIF2AK4) and significant allelic imbalance.

Our study reveals, that PE was followed by downregulation of Defender Against Cell Death 1 (DAD1), which is one of the genes that contribute to the synthesis of oligosaccharyltransferase (OST) complex that catalyzes the transfer of a high mannose oligosaccharide onto asparagine residues. It has been found that disruptions in the OST complex induce chronic ER stress by activating one of the pathways in unfolded protein response (UPR), leading to caspase-dependent apoptosis [71]. Thus, most possibly, the underexpression of DAD1 observed in the PA-E libraries represents cell death initiated by severe ER imbalance induced by PE. The downregulation of Lectin Mannose Binding 2 (LMAN2) and Ribosome Binding Protein 1 (RRBP1) disclosed in the present study, also reveal distinct disruption in ER protein processing. LMAN2 is an intracellular protein that shuttles between the Golgi apparatus and ER, responsible for recognizing misfolded proteins that are produced under the ER stress conditions [72]. In turn, RRBP1, underexpressed in PA-E arteries, takes part in polysome assembly and therefore in protein biosynthesis along with the expansion of the Golgi complex. Additionally, RRBP1 may influence protein synthesis, necessary in some unfolded protein response (UPR) pathways [73]. Another underexpressed gene Dolichyl-Diphosphooligosaccharide Protein Glycosyltransferase (DDOST; also known as AGER1) is a receptor of advanced glycation end products (AGEs) and it can decrease prooxidant actions [74]. Accumulation of AGEs, resulting from DDOST overexpression, triggers platelets hyperreactivity and in specific circumstances may promote thrombosis [75]. Therefore, correspondingly, the observed DDOST underexpression most possibly should be considered as a protective mechanism, decreasing platelets activation, thus, leading to blood clot degradation. Disulfide bond formation, a critical step in protein folding which occurs within the ER, is catalyzed by the protein disulfide isomerase family (PDI). One of the PDI members Protein Disulfide Isomerase Family A Member 4 (PDIA4) activates platelets and promotes thrombosis initiation [76]. Accumulation of abnormally folded proteins and increased PDIA4 expression trigger ER response. Additionally, PDIA4 is highly transcribed in the megakaryocytes and effectively joins with the cell surface following platelet activation [77]. Therefore, most possibly underexpression of PDIA4 observed in the PA-E samples, similarly to the underexpression of DDOST reflects the protective mechanism limiting further clot formation in the affected PA.

So far, there are no reports supporting the direct effect of Ubiquilin 4 (UBQLN4) and DnaJ Heat Shock Protein Family (Hsp40) Member B11 (DNAJB11) on the course of PE. Both these genes, which are downregulated in PA-E libraries, are associated with cytosolic chaperones, major regulators of ER protein folding and ER stress [78]. Cytosolic chaperone UBQLN4 belongs to the ubiquitin-binding proteins, which exhibit diverse biological functions, including enhancement of autophagy-mediated degradation and participation in ER-associated protein degradation [79]. DNAJB11 encodes one of the main cofactors of the ER chaperon, a heat-shock protein that performs multifaceted functions in coordinating ER and extracellular proteostasis, required for efficient protein folding and trafficking [80]. Furthermore, Protein Disulfide Isomerase Family A Member 6 (PDIA6) is an ER-resident protein that plays a crucial role in blocking the onset of ER stress, as it suppresses the activation of the UPR [81, 82], therefore the observed UBQLN4, DNAJB11 and PDIA6 underexpression, again reflects the ongoing ER stress response in the PA-E. Another downregulated gene, calreticulin (CALR), is well-known for sequestering calcium within the ER lumen [83]. Mutations within CALR have been shown to occur in 67% of essential thrombocythemia (ET) cases. ET is a myeloproliferative neoplasm characterized by a clonal expansion of megakaryocytes, which can result in both arterial and venous thrombosis [84].

Our research revealed two overexpressed genes involved in the ‘protein processing in endoplasmic reticulum’ pathway—NFE2 Like BZIP Transcription Factor 2 (NFE2L2) and Membrane Associated Ring-CH-Type Finger 6 (MARCHF6). NFE2L2 controls cellular response to oxidative insults, thus preventing damage to components sensitive to redox changes such as proteins, lipids and DNA. NFE2L2 is a putative regulator of the cardiac stress-response [85], and it has been found that NFE2L2 gene polymorphisms are associated with the risk and severity of acute type A aortic dissection [86]. MARCHF6 is a ubiquitin ligase embedded in the membranes of the ER, which participates in ER-associated degradation, including autoubiquitination [87]. Moreover, MARCHF6 is involved in vascular endothelial-cadherin-based adherens junctions function and sprouting angiogenesis [88].

This project examined also the relationship between several genes within the identified ER-KEGG pathway and genes with alternative splicing bias, exposed by DASGs. For instance, a significant allelic imbalance was detected in ER Degradation Enhancing Alpha-Mannosidase Like Protein 2 (EDEM2), Mitogen-Activated Protein Kinase 8 (MAPK8; also known as JNK1 or JNK), TNF Receptor Associated Factor 2 (TRAF2) or in previously described UBQLN4 genes. It has been found, that SNPs with the genome-wide significance of the EDEM2, influence protein C concentration, which is an important endogenous anticoagulant in hemostasis and its deficiencies significantly increase the risk of VTE [89]. MAPK8 is a member of the mitogen-activated protein kinase family (MAPK). MAPK signaling pathway plays a vital role in the pathogenesis of cardiovascular diseases, for example, MAPK8 was identified as the target gene in myocardial infarction [90].

As was already mentioned in the Results section, in the PA-E libraries we have found not only DEGs directly involved in the ‘protein processing in endoplasmic reticulum’ pathway but we have also revealed 61 polymorphisms localized in the vicinity of the genes belonging to the cellular components of ER, for example, ITGA8, CCDC3, ITPR1, PKD2, BNIP3L, SEC23B, CERS5, ACER3, KDELR3, PML, SEC31A, THY1 and BCAP29. These genes, similar to the previously described DEGs, are entangled in the protein biosynthesis process and ER stress-induced apoptosis.

TNFα signaling pathway

A previous analysis of the changes in the transcriptome of the femoral vein caused by DVT in the porcine model based on the formation of the thrombus in vivo, disclosed that TNF, NF-κB and apoptosis pathways are strongly dysregulated during the thrombus formation [16]. The present study also revealed genes with significant allelic imbalance involved in the TNFα signaling pathway, such as the Coiled-coil domain containing 3 (CCDC3), Thy-1 Cell Surface Antigen (Thy1) and TRAF2. CCDC3 is highly expressed in vascular endothelial cells where it is involved in the regulation of TNFα-induced inflammatory response [91]. It has been found, that TNFα, which is responsible for triggering the expression of pro-inflammatory genes by activation of NF-κB, downregulates CCDC3 expression in endothelial cells [92]. Moreover, Azad et al. [91] have proven that the interaction between TNFα and CCDC3 is bidirectional, and an increased CCDC3 level inhibits NF-κB activation and thus decreased TNFα-induced endothelial gene expression. Thy1 is expressed in microvascular endothelial cells during angiogenesis [93]. The evidence suggests that Thy-1 induced on the endothelium regulates vascular permeability at sites of inflammation [94]. The expression on endothelial cells, the inducibility by TNF, and the ability to transmit signals suggest that Thy-1 may again have an important role in inflammatory responses [95]. TRAF2 encodes proteins belonging to the TNF receptor-associated factor (TRAF) protein family, which mediate the transduction of signals from the tumor necrosis factor (TNF) receptor superfamily. This protein is required for TNFα—mediated activation of MAPK8 and nuclear factor κB (NF-κB), for example in response to ER stress. The activation of MAPK8 is essential for the regulation of TNFα induced apoptosis [96]. Activated NF-κB is a factor that modifies the transcription of numerous target genes encoding pro-inflammatory cytokines, chemokines, cell adhesion molecules, and enzymes that produce pro-inflammatory factors, such as nitric oxide or prostaglandins.

The comparison of the present data with the results obtained in other animal and human models

To gain a deeper insight into the molecular background of pathomechanisms that lead to PE currently detected DEGs were compared with the results of four other research performed either in humans or in other animal models (Fig. S1). These studies included: an analysis of the human VTE patients’ microarray dataset [53], an investigation of the differential effect of pulmonary vascular occlusion during mild and severe PE on the lung transcriptome in the rat [52], data mining of public gene expression datasets from VTE and other cardiovascular diseases (i.e. myocardial infarction, peripheral arterial occlusive disease and stroke; [54]) and genome-wide gene expression profiling in rabbit autologous thrombus PE model [14]. Intriguingly, we have found twenty-nine overlapping genes which were also differentially expressed in the present study. Four of these DEGs (ACSL4, FKBP14, PDIA6, COL1A2), were already mentioned in the Discussion section. However, the most striking observation to emerge from the data comparison was that Procollagen-Lysine,2-Oxoglutarate 5-Dioxygenase 2 (PLOD2) gene was uncovered in all of the analyzed datasets, which may indicate its crucial role in PE pathogenesis. PLOD2 protein is involved in the stability of collagen fibrils, that in turn plays a key role in the maintenance of integrity and elasticity of normal as well as atherosclerotic vessel walls and it stimulates the thrombus formation processes [97]. Although compared DEGs were associated with distinct pathogenic mechanisms obtained data suggest that the mechanisms responsible for cardiovascular diseases such as PE and VTE are at least partially overlapping. Thus, such interspecies comparison, including results from various methodological approaches may constitute the fundamental framework for future clinical research on human cardiovascular diseases, hopefully leading to the identification of potential biomarkers for prognosis or therapeutic targets.

Potential limitations

Like all experimental models, ours has some limitations that must be recognized. Unlike humans, animal species do not develop spontaneous VTE, therefore despite the fact that pigs share many similarities with humans, the creation of thrombus in their veins requires the application of external factors. Thus, it is difficult to entirely relate also PE pathomechanism in pig to man. However, it must be stressed that this problem refers to all VTE animal models.

Conclusions

This research provides the first thorough investigation of the changes in the genes expression profile and other posttranscriptional modifications affected by the blood clot lodged in the PA in the course of PE. The current data revealed that dysregulation of protein biosynthesis in ER plays a crucial role in PE pathogenesis. The differentially expressed genes associated with the examined ER pathway may exert a twofold impact on the course of PE, resulting from both their main function and being a link between disturbed ‘protein processing in endoplasmic reticulum’ pathway and ER stress response induced by the clot. Downregulated genes were mostly involved in ER primary function, that is protein synthesis and maturation (including their cleavage, glycosylation, folding and assembly). Whereas, the upregulated genes were involved in the regulation of the redox status, sterol, and lipid metabolism. Furthermore, the present study revealed DEGs classified as signatures of the dysfunction and abnormality of the cardiovascular system. The findings of this research provide insights for understanding PE pathogenesis and its complications. Therefore, we hope that our research will lay the groundwork for future clinical trials investigating new therapeutic agents.

Availability of data and materials

Data relating to this article can be found in the article itself or online supplementary material. Raw data were submitted to the European Nucleotide Archive (ENA; (https://www.ebi.ac.uk/ena/browser/view/PRJEB43020) under accession No. PRJEB43020 and used for transcriptome profiling.

References

  1. Byrnes JR, Wolberg AS. New findings on venous thrombogenesis. Hamostaseologie. 2017;37:25–35.

    Article  Google Scholar 

  2. Aggarwal A, Fullam L, Brownstein AP, Maynard GA, Ansell J, Varga EA, et al. Deep Vein Thrombosis (DVT) and Pulmonary Embolism (PE): Awareness and Prophylaxis Practices Reported by Patients with Cancer. Cancer Invest. 2015;33:405–10.

    Article  Google Scholar 

  3. Ouriel K, Green RM, Greenberg RK, Clair DG. The anatomy of deep venous thrombosis of the lower extremity. J Vasc Surg. 2000;31:895–900.

    Article  CAS  Google Scholar 

  4. Furie B, Furie BC. Mechanisms of thrombus formation. Mechanisms of Disease. N Engl J Med. 2008;359:938–49.

    Article  CAS  Google Scholar 

  5. van Hinsbergh VWM. Endothelium - role in regulation of coagulation and inflammation. Semin Immunopathol. 2012;34:93–106.

    Article  Google Scholar 

  6. Peng Z, Shu B, Zhang Y, Wang M. Endothelial response to pathophysiological stress. Arterioscler Thromb Vasc Biol. 2019;39:e233–43.

    Article  CAS  Google Scholar 

  7. Turetz M, Sideris A, Friedman O, Triphathi N, Horowitz J. Epidemiology, pathophysiology, and natural history of pulmonary embolism. Semin Intervent Radiol. 2018;35:92–8.

    Article  Google Scholar 

  8. Douketis JD. Prognosis in pulmonary embolism. Curr Opin Pulm Med. 2001;7:354–9.

    Article  CAS  Google Scholar 

  9. Jaff MR, McMurtry MS, Archer SL, Cushman M, Goldenberg N, Goldhaber SZ, et al. Management of massive and submassive pulmonary embolism, iliofemoral deep vein thrombosis, and chronic thromboembolic pulmonary hypertension. Circulation. 2011;123:1788–830.

    Article  Google Scholar 

  10. Martinez Licha CR, McCurdy CM, Maldonado SM, Lee LS. Current management of acute pulmonary embolism. Ann Thorac Cardiovasc Surg. 2020;26:65–71.

    Article  Google Scholar 

  11. Weitz JI, Prandoni P, Verhamme P. Anticoagulation for patients with venous thromboembolism: when is extended treatment required? TH Open. 2020;4:e446–56.

    Article  Google Scholar 

  12. Konstantinides S. Should thrombolytic therapy be used in patients with pulmonary embolism? Am J Cardiovasc Drugs. 2004;4:69–74.

    Article  Google Scholar 

  13. Pai R, Fang Q, Tian G, Zhu B, Ge X. Expression and role of interleukin-1β and associated biomarkers in deep vein thrombosis. Exp Ther Med. 2021;22:1–8.

    Article  Google Scholar 

  14. Tang Z, Wang X, Huang J, Zhou X, Xie H, Zhu Q, et al. Gene expression profiling of pulmonary artery in a rabbit model of pulmonary thromboembolism. PLoS ONE. 2016;11:1–16.

    Article  Google Scholar 

  15. Feng R, Lu M, Xu J, Zhang F, Yang M, Luo P, et al. Pulmonary embolism and 529 human blood metabolites: genetic correlation and two-sample Mendelian randomization study. BMC Genom Data. 2022;23:1–9.

    Article  Google Scholar 

  16. Gromadziński L, Paukszto Ł, Skowrońska A, Holak P, Smoliński M, Łopieńska-Biernat E, et al. Transcriptomic profiling of femoral veins in deep vein thrombosis in a porcine model. Cells. 2021;10:1–19.

    Article  Google Scholar 

  17. Diaz JA, Obi AT, Myers DD, Wrobleski SK, Henke PK, Mackman N, et al. Critical Review of Mouse Models of Venous Thrombosis. Arterioscler Thromb Vasc Biol. 2012;32:556–62.

    Article  CAS  Google Scholar 

  18. Shi WY, Wu S, Hu LY, Liu CJ, Gu JP. Swine model of thrombotic caval occlusion created by autologous thrombus injection with assistance of intra-caval net knitting. Sci Rep. 2015;5:1–9.

    Article  Google Scholar 

  19. Albadawi H, Witting AA, Pershad Y, Wallace A, Fleck AR, Hoang P, et al. Animal models of venous thrombosis. Cardiovasc Diagn Ther. 2017;7(Suppl 3):S197-206.

    Article  Google Scholar 

  20. Lin PH, Chen C, Surowiec SM, Conklin B, Bush RL, Lumsden AB. Evaluation of thrombolysis in a porcine model of chronic deep venous thrombosis: an endovascular model. J Vasc Surg. 2001;33:621–7.

    Article  CAS  Google Scholar 

  21. Rectenwald JE, Deatrick KB, Sukheepod P, Lynch EM, Moore AJ, Moaveni DM, et al. Experimental pulmonary embolism: effects of the thrombus and attenuation of pulmonary artery injury by low-molecular-weight heparin. J Vasc Surg. 2006;43:800–8.

    Article  Google Scholar 

  22. Barbash IM, Schenke WH, Halabi M, Ratnayaka K, Faranesh AZ, Kocaturk O, et al. Experimental model of large pulmonary embolism employing controlled release of subacute caval thrombus in swine. J Vasc Interv Radiol. 2011;22:1471–7.

    Article  Google Scholar 

  23. Frisbie JH. An animal model for venous thrombosis and spontaneous pulmonary embolism. Spinal Cord. 2005;43:635–9.

    Article  CAS  Google Scholar 

  24. Robinson VJB, Pineda GE, Salah AK, Pipkin WL, Corley JH, Jonah MH, et al. Latex d-dimer signal in in situ femoral vein thrombus in swine and effect of minidose exogenous tissue plasminogen activator bolus. Chest. 2005;127:622–9.

    Article  CAS  Google Scholar 

  25. Gromadziński L, Skowrońska A, Holak P, Smoliński M, Lepiarczyk E, Żurada A, et al. A new experimental porcine model of venous thromboembolism. J Clin Med. 2021;10:1862.

    Article  Google Scholar 

  26. Andrews S. FastQC: A quality control tool for high throughput sequence data. Babraham Bioinformatics. 2010;:www.bioinformatics.babraham.ac.uk/projects/fastqc.

  27. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.

    Article  CAS  Google Scholar 

  28. Pertea M, Pertea GM, Antonescu CM, Chang T-C, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33:290–5.

    Article  CAS  Google Scholar 

  29. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:1–21.

    Article  Google Scholar 

  30. Wang L, Park HJ, Dasari S, Wang S, Kocher J-P, Li W. CPAT: Coding-Potential Assessment Tool using an alignment-free logistic regression model. Nucleic Acids Res. 2013;41:e74–e74.

    Article  CAS  Google Scholar 

  31. Kang YJ, Yang DC, Kong L, Hou M, Meng YQ, Wei L, et al. CPC2: A fast and accurate coding potential calculator based on sequence intrinsic features. Nucleic Acids Res. 2017;45:W12–6.

    Article  CAS  Google Scholar 

  32. Wucher V, Legeai F, Hédan B, Rizk G, Lagoutte L, Leeb T, et al. FEELnc: a tool for long non-coding RNA annotation and its application to the dog transcriptome. Nucleic Acids Res. 2017;45:8.

  33. Frazee AC, Pertea G, Jaffe AE, Langmead B, Salzberg SL, Leek JT. Flexible analysis of transcriptome assemblies with Ballgown. bioRxiv. 2014;003665.

  34. Wickham H. ggplot2: Elegant Graphics for Data Analysis. New York: Springer-Verlag; 2016.

    Book  Google Scholar 

  35. Shen S, Park JW, Lu Z-X, Lin L, Henry MD, Wu YN, et al. rMATS: Robust and flexible detection of differential alternative splicing from replicate RNA-Seq data. Proc Natl Acad Sci. 2014;111:E5593–601.

    Article  CAS  Google Scholar 

  36. Luo W, Brouwer C. Pathview: an R/Bioconductor package for pathway-based data integration and visualization. Bioinformatics. 2013;29:1830–1.

    Article  CAS  Google Scholar 

  37. Veiga D. maser: Mapping Alternative Splicing Events to pRoteins. R package version 1.0.0,. 2018. https://github.com/DiogoVeiga/maser.

  38. Garrido-Martín D, Palumbo E, Guigó R, Breschi A. ggsashimi: Sashimi plot revised for browser- and annotation-independent splicing visualization. PLoS Comput Biol. 2018;14: e1006360.

    Article  Google Scholar 

  39. McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, et al. The genome analysis toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 2010;20(9):1297–303. https://doi.org/10.1101/gr.107524.110.

    Article  CAS  Google Scholar 

  40. Wang J, Pan Y, Shen S, Lin L, Xing Y. rMATS-DVR: rMATS discovery of differential variants in RNA. Bioinformatics. 2017;33:2216–7.

    Article  CAS  Google Scholar 

  41. Wang X, Lu P, Luo Z. GMATo: A novel tool for the identification and analysis of microsatellites in large genomes. Bioinformation. 2013;9:541–4.

    Article  Google Scholar 

  42. Kent WJ. BLAT —The BLAST -Like Alignment Tool. Genome Res. 2002;12:656–64.

    CAS  Google Scholar 

  43. Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–2.

    Article  CAS  Google Scholar 

  44. Knaus BJ, Grünwald NJ. vcfr: a package to manipulate and visualize variant call format data in R. Mol Ecol Resour. 2017;17:44–53.

  45. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics. 2009;25:2078–9.

    Article  Google Scholar 

  46. McLaren W, Gil L, Hunt SE, Riat HS, Ritchie GRS, Thormann A, et al. The Ensembl variant effect predictor. Genome Biol. 2016;17:122.

    Article  Google Scholar 

  47. Krzywinski M, Schein J, Birol I, Connors J, Gascoyne R, Horsman D, et al. Circos: an information aesthetic for comparative genomics. Genome Res. 2009;19:1639–45.

    Article  CAS  Google Scholar 

  48. Reimand J, Arak T, Adler P, Kolberg L, Reisberg S, Peterson H, et al. g:Profiler—a web server for functional interpretation of gene lists. Nucleic Acids Res. 2016;44:W83–9.

    Article  CAS  Google Scholar 

  49. Kanehisa M, Furumichi M, Tanabe M, Sato Y, Morishima K. KEGG: New perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res. 2017;45:D353–61.

    Article  CAS  Google Scholar 

  50. Untergasser A, Cutcutache I, Koressaar T, Ye J, Faircloth BC, Remm M, et al. Primer3—new capabilities and interfaces. Nucleic Acids Res. 2012;40:e115–e115.

    Article  CAS  Google Scholar 

  51. Pfaffl MW. A new mathematical model for relative quantification in real-time RT-PCR. Nucleic Acids Res. 2001;29:45e–45.

    Article  Google Scholar 

  52. Zagorski J, Kline JA. Differential effect of mild and severe pulmonary embolism on the rat lung transcriptome. Respir Res. 2016;17:1–13.

    Article  Google Scholar 

  53. Su Y, Li Q, Zheng Z, Wei X, Hou P. Identification of genes, pathways and transcription factor-miRNA-target gene networks and experimental verification in venous thromboembolism. Sci Rep. 2021;11:1–16.

    Article  CAS  Google Scholar 

  54. Hounkpe BW, de Oliveira Benatti R, de Sá Carvalho B, de Paula EV. Identification of common and divergent gene expression signatures in patients with venous and arterial thrombosis using data from public repositories. PLoS One. 2020;15(8):e0235501.

    Article  CAS  Google Scholar 

  55. Plenz GAM, Deng MC, Robenek H, Völker W. Vascular collagens: spotlight on the role of type VIII collagen in atherogenesis. Atherosclerosis. 2003;166:1–11.

    Article  CAS  Google Scholar 

  56. Kuivaniemi H, Tromp G, Prockop DJ. Mutations in fibrillar collagens (types I, II, III, and XI), fibril- associated collagen (type IX), and network-forming collagen (type X) cause a spectrum of disease of bone, cartilage, and blood vessels. Hum Mutat. 1997;9:300–15.

    Article  CAS  Google Scholar 

  57. Vuust J, Sobel ME, Martin GR. Regulation of type I collagen synthesis. total pro alpha 1(I) and pro alpha 2(I) mRNAs are maintained in a 2:1 ratio under varying rates of collagen synthesis. Eur J Biochem. 1985;151:449–53.

    Article  CAS  Google Scholar 

  58. Žaloudíková M, Eckhardt A, Vytášek R, Uhlík J, Novotný T, Bačáková L, et al. Decreased collagen VI in the tunica media of pulmonary vessels during exposure to hypoxia: a novel step in pulmonary arterial remodeling. Pulm Circ. 2019;9:1–10.

    Article  Google Scholar 

  59. Luther DJ, Thodeti CK, Shamhart PE, Adapala RK, Hodnichak C, Weihrauch D, et al. Absence of type vi collagen paradoxically improves cardiac function, structure, and remodeling after myocardial infarction. Circ Res. 2012;110:851–6.

    Article  CAS  Google Scholar 

  60. Lintner KE, Wu YL, Yang Y, Spencer CH, Hauptmann G, Hebert LA, et al. Early Components of the Complement Classical Activation Pathway in Human Systemic Autoimmune Diseases. Front Immunol. 2016;7:36.

  61. Ekdahl KN, Teramura Y, Hamad OA, Asif S, Duehrkop C, Fromell K, et al. Dangerous liaisons: complement, coagulation, and kallikrein/kinin cross-talk act as a linchpin in the events leading to thromboinflammation. Immunol Rev. 2016;274:245–69.

    Article  CAS  Google Scholar 

  62. Carter AM. Complement activation: an emerging player in the pathogenesis of cardiovascular disease. Scientifica (Cairo). 2012;2012:1–14.

    Article  Google Scholar 

  63. Matissek SJ, Elsawa SF. GLI3: a mediator of genetic diseases, development and cancer. Cell Commun Signal. 2020;18:54.

    Article  CAS  Google Scholar 

  64. Golej DL, Askari B, Kramer F, Barnhart S, Vivekanandan-Giri A, Pennathur S, et al. Long-chain acyl-CoA synthetase 4 modulates prostaglandin E2 release from human arterial smooth muscle cells. J Lipid Res. 2011;52:782–93.

    Article  CAS  Google Scholar 

  65. Gibbins JM, Briddon S, Shutes A, Van Vugt MJ, Van De Winkel JGJ, Saito T, et al. The p85 Subunit of Phosphatidylinositol 3-Kinase Associates with the Fc Receptor γ-Chain and Linker for Activitor of T Cells (LAT) in Platelets Stimulated by Collagen and Convulxin *. J Biol Chem. 1998;273:34437–43.

    Article  CAS  Google Scholar 

  66. Durrant TN, Hers I. PI3K inhibitors in thrombosis and cardiovascular disease. Clin Transl Med. 2020;9:8.

    Article  Google Scholar 

  67. Schwarz DS, Blower MD. The endoplasmic reticulum: Structure, function and response to cellular signaling. Cell Mol Life Sci. 2016;73:79–94.

    Article  CAS  Google Scholar 

  68. Lin JH, Walter P, Yen TSB. Endoplasmic reticulum stress in disease pathogenesis. Annu Rev Pathol. 2008;3:399–425.

    Article  CAS  Google Scholar 

  69. Lemmer IL, Willemsen N, Hilal N, Bartelt A. A guide to understanding endoplasmic reticulum stress in metabolic disorders. Mol Metab. 2021;47: 101169.

    Article  CAS  Google Scholar 

  70. Sabirli R, Koseler A, Mansur N, Zeytunluoglu A, Sabirli GT, Turkcuer I, et al. Predictive value of endoplasmic reticulum stress markers in low ejection fractional heart failure. In Vivo (Brooklyn). 2019;33:1581–92.

    Article  CAS  Google Scholar 

  71. Zhang Y, Cui C, Lai ZC. The defender against apoptotic cell death 1 gene is required for tissue growth and efficient N-glycosylation in Drosophila melanogaster. Dev Biol. 2016;420:186–95.

    Article  CAS  Google Scholar 

  72. Yu B, Xu C, Tang X, Liu Z, Lin X, Meng H, et al. Endoplasmic reticulum stress-related secretory proteins as biomarkers of early myocardial ischemia-induced sudden cardiac deaths. Int J Legal Med. 2022;136:159–68.

    Article  Google Scholar 

  73. Ortega A, Roselló-Lletí E, Tarazón E, Molina-Navarro MM, Martínez-Dolz L, González-Juanatey JR, et al. Endoplasmic reticulum stress induces different molecular structural alterations in human dilated and ischemic cardiomyopathy. PLoS ONE. 2014;9: e107635.

    Article  Google Scholar 

  74. Liu X, Zhang R, Huang L, Zheng Z, Vlassara H, Striker G, et al. Excessive oxidative stress contributes to increased acute ER stress kidney injury in aged mice. Oxid Med Cell Longev. 2019;2019:2746521.

    Google Scholar 

  75. Zhu W, Li W, Silverstein RL. Advanced glycation end products induce a prothrombotic phenotype in mice via interaction with platelet CD36. Blood. 2012;119:6136–44.

    Article  CAS  Google Scholar 

  76. Wang Z, Zhang H, Cheng Q. PDIA4: the basic characteristics, functions and its potential connection with cancer. Biomed Pharmacother. 2019;2020(122): 109688.

    Google Scholar 

  77. Holbrook LM, Watkins NA, Simmonds AD, Jones CI, Ouwehand WH, Gibbins JM. Platelets release novel thiol isomerase enzymes which are recruited to the cell surface following activation. Br J Haematol. 2010;148:627–37.

    Article  CAS  Google Scholar 

  78. Gidalevitz T, Stevens F, Argon Y. Orchestration of secretory protein folding by ER chaperones. Biochem Biophys Acta. 2013;1833:2410–24.

    Article  CAS  Google Scholar 

  79. Yun Lee D, Arnott D, Brown EJ. Ubiquilin4 is an adaptor protein that recruits Ubiquilin1 to the autophagy machinery. EMBO Rep. 2013;14:373–81.

    Article  Google Scholar 

  80. Huynh VT, Audrézet MP, Sayer JA, Ong AC, Lefevre S, Le Brun V, et al. Clinical spectrum, prognosis and estimated prevalence of DNAJB11-kidney disease. Kidney Int. 2020;98:476–87.

    Article  CAS  Google Scholar 

  81. Kiouptsi K, Finger S, Garlapati VS, Knorr M, Brandt M, Walter U, et al. Hypoxia evokes increased PDI and PDIA6 expression in the infarcted myocardium of ex-germ-free and conventionally raised mice. Biol Open. 2019;8:1–7.

    Google Scholar 

  82. Souza-Neto FV, Jiménez-González S, Delgado-Valero B, Jurado-López R, Genty M, Romero-Miranda A, et al. The interplay of mitochondrial oxidative stress and endoplasmic reticulum stress in cardiovascular fibrosis in obese rats. Antioxidants. 2021;10:1274.

    Article  CAS  Google Scholar 

  83. Biwer LA, Good ME, Hong K, Patel RK, Agrawal N, Looft-Wilson R, et al. Non-endoplasmic reticulum-based calr (Calreticulin) can coordinate heterocellular calcium signaling and vascular function. Arterioscler Thromb Vasc Biol. 2018;38:120–30.

    Article  CAS  Google Scholar 

  84. Nazha B, Garcia G, Kandov R, Odaimi M. Calreticulin Mutated Essential Thrombocythemia Presenting as Acute Coronary Syndrome. Case Rep Hematol. 2015;2015:161764.

  85. Quiles JM, Pepin ME, Sunny S, Shelar SB, Challa AK, Dalley B, et al. Identification of Nrf2-responsive microRNA networks as putative mediators of myocardial reductive stress. Sci Rep. 2021;11:1–12.

    Article  Google Scholar 

  86. Zhang Y, Zheng Q, Chen R, Dai X, Zhu Y, Ma L. Association of NFE2L2 Gene Polymorphisms with Risk and Clinical Characteristics of Acute Type A Aortic Dissection in Han Chinese Population. Oxid Med Cell Longev. 2021;2021:5173190.

  87. Scott NA, Sharpe LJ, Brown AJ. BBA - Molecular and Cell Biology of Lipids The E3 ubiquitin ligase MARCHF6 as a metabolic integrator in cholesterol synthesis and beyond. BBA - Molecular and Cell Biology of Lipids. 2021;1866: 158837.

    Article  CAS  Google Scholar 

  88. Tan JME, van der Stoel MM, van den Berg M, van Loon NM, Moeton M, Scholl E, et al. The MARCH6-SQLE Axis Controls Endothelial Cholesterol Homeostasis and Angiogenic Sprouting. Cell Rep. 2020;32(5):107944.

  89. Cooper DN, Tuddenham EGD. Molecular genetics of familial venous thrombosis. Br Med Bull. 1994;50:833–50.

    Article  CAS  Google Scholar 

  90. Wang Q, Liu B, Wang Y, Bai B, Yu T, Chu X. The biomarkers of key miRNAs and target genes associated with acute myocardial infarction. PeerJ. 2020;8: e9129.

    Article  Google Scholar 

  91. Azad AK, Chakrabarti S, Xu Z, Davidge ST, Fu YX. Coiled-coil domain containing 3 (CCDC3) represses tumor necrosis factor-α/nuclear factor κB-induced endothelial inflammation. Cell Signal. 2014;26:2793–800.

    Article  CAS  Google Scholar 

  92. Kobayashi S, Fukuhara A, Taguchi T, Matsuda M, Tochino Y, Otsuki M, et al. Identification of a new secretory factor, CCDC3/Favine, in adipocytes and endothelial cells. Biochem Biophys Res Commun. 2010;392:29–35.

    Article  CAS  Google Scholar 

  93. Lee W Sen, Jain MK, Arkonac BM, Zhang D, Shaw SY, Kashiki S, et al. Thy-1, a novel marker for angiogenesis upregulated by inflammatory cytokines. Circ Res. 1998;82:845–51.

  94. Ishizu A, Ishikura H, Nakamaru Y, Takeuchi E, Kimura C, Koike T, et al. Thy-1 induced on rat endothelium regulates vascular permeability at sites of inflammation. Int Immunol. 1995;7:1939–47.

    Article  CAS  Google Scholar 

  95. Mason JC, Yarwood H, Tárnok A, Sugars K, Harrison AA, Robinson PJ, et al. Human Thy-1 is cytokine-inducible on vascular endothelial cells and is a signaling molecule regulated by protein kinase C. J Immunol. 1996;157:874–83.

    Article  CAS  Google Scholar 

  96. Urano F, Wang XZ, Bertolotti A, Zhang Y, Chung P, Harding HP, et al. Coupling of stress in the ER to activation of JNK protein kinases by transmembrane protein kinase IRE1. Science. 1979;2000(287):664–6.

    Google Scholar 

  97. Nukala SB, Tura-Ceide O, Aldini G, Smolders VFED, Blanco I, Peinado VI, et al. Protein network analyses of pulmonary endothelial cells in chronic thromboembolic pulmonary hypertension. Sci Rep. 2021;11:5583.

Download references

Funding

Leszek Gromadziński is the recipient of the Statutory Fund of the School of Medicine, Collegium Medicum (61.610.001–300), University of Warmia and Mazury in Olsztyn.

Author information

Authors and Affiliations

  1. Department of Cardiology and Internal Medicine, School of Medicine, Collegium Medicum, University of Warmia and Mazury in Olsztyn, Warszawska Str 30, 10-082, Olsztyn, Poland

    Leszek Gromadziński

  2. Department of Botany and Nature Protection, University of Warmia and Mazury in Olsztyn, Plac Łódzki 1, 10-727, Olsztyn, Poland

    Łukasz Paukszto

  3. Department of Human Physiology and Pathophysiology, School of Medicine, Collegium Medicum, University of Warmia and Mazury in Olsztyn, Warszawska Str 30, 10-082, Olsztyn, Poland

    Ewa Lepiarczyk, Agnieszka Skowrońska & Marta Majewska

  4. Department of Gynecology, and Obstetrics, School of Medicine, Collegium Medicum, University of Warmia and Mazury in Olsztyn, Żołnierska Str 18, 10-561, Olsztyn, Poland

    Aleksandra Lipka

  5. Department of Animal Anatomy and Physiology, Faculty of Biology and Biotechnology, University of Warmia and Mazury in Olsztyn, Oczapowskiego 1A, 10-719, Olsztyn, Poland

    Karol G. Makowczenko

  6. Department of Biochemistry, Faculty of Biology and Biotechnology, University of Warmia and Mazury in Olsztyn, Oczapowskiego Str 1A, 10-719, Olsztyn, Poland

    Elżbieta Łopieńska-Biernat

  7. Department of Plant Physiology, Genetics and Biotechnology, Faculty of Biology and Biotechnology, University of Warmia and Mazury in Olsztyn, Oczapowskiego Str 1A, 10-719, Olsztyn-Kortowo, Poland

    Jan P. Jastrzębski

  8. Department of Surgery and Radiology With Clinic, Faculty of Veterinary Medicine, University of Warmia and Mazury in Olsztyn, Oczapowskiego Str 14, 10-719, Olsztyn, Poland

    Piotr Holak

  9. Clinic of Cardiology and Internal Diseases, University Clinical Hospital in Olsztyn, Warszawska Str 30, 10-082, Olsztyn, Poland

    Michał Smoliński

Authors
  1. Leszek Gromadziński
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Łukasz Paukszto
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Ewa Lepiarczyk
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Agnieszka Skowrońska
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Aleksandra Lipka
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Karol G. Makowczenko
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Elżbieta Łopieńska-Biernat
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Jan P. Jastrzębski
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Piotr Holak
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Michał Smoliński
    View author publications

    You can also search for this author in PubMed Google Scholar

  11. Marta Majewska
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

Conceptualization: L.G. and M.M. methodology: L.G., A.L., P.H., M.S. and M.M. formal analysis: Ł.P., K.G.M. validation: E.Ł.-B. investigation: L.G., A.S., E.L., A.L. and M.M. resources: L.G. writing—original draft preparation: L.G., Ł.P., E.L., A.S., A.L. and M.M. writing—review and editing: L.G., Ł.P., E.L., A.S., A.L. and M.M. visualization: J.P.J. supervision: L.G. and M.M. project administration: L.G. and M.M. funding acquisition: L.G. The author(s) read and approved the final manuscript.

Corresponding author

Correspondence to Marta Majewska.

Ethics declarations

Ethics approval and consent to participate

Animal welfare and experimental procedures were performed in accordance with the guidelines of the local Ethics Committee for Animal Experimentation in Olsztyn (affiliated with the National Ethics Committee for Animal Experimentation, Polish Ministry of Science and Higher Education; decision No. 90/2018 of 13.02.2018). The study was reported in accordance with ARRIVE guidelines (https://arriveguidelines.org). All procedures conformed to the guidelines from Directive 2010/63/EU of the European Parliament on the protection of animals used for scientific purposes or the NIH Guide for the Care and Use of Laboratory Animals.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

Sequences of primers used in Real-Time PCR validation for differentially expressed protein-coding genes and reference genes.

Additional file 2: Table S2.

Differentially expressed protein-coding genes identified in the porcine pulmonary arteries affected by embolism.

Additional file 3: Table S3.

Differentially expressed long noncoding RNAs identified in the porcine pulmonary arteries affected by embolism.

Additional file 4: Table S4.

Differentially alternative splicing events discovered by rMATS tools. Five alternative splicing event types show an imbalance in the percent splicing inclusion (PSI) values between CTR and PA-E libraries.

Additional file 5: Table S5.

Differences in alternative allele fraction for significant allele-specific expression sites identified in the porcine pulmonary arteries affected by embolism.

Additional file 6: Table S6.

Differentially expressed genes enriched in the Gene Ontology annotations.

Additional file 7: Table S7.

Differentially alternative splicing genes classified to the Gene Ontology terms.

Additional file 8: Figure S1.

Venn diagram of DEGs comparing the obtained results with other datasets focusing on expression analysis of the DVT and PE. The numbers within the intersected ovals indicate the count of DEGs common within experiments. The outer ovals depict DEGs which are unique in the particular research.

Additional file 9: Figure S2.

The mRNA expression of selected genes obtained using Real-Time PCR. The expression of endogenous controls is shown as normalised to a value 1, and the samples indicate the changes relative to the controls. All the relative expression units [a. u.] were presented as logarithmic values. The exact values of expression are above the bars. P-values were considered statistically significant where 0.0332 (*), 0.0021 (**) and <0.0001 (****); ns - not significant, but tended to have expression changes.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Gromadziński, L., Paukszto, Ł., Lepiarczyk, E. et al. Pulmonary artery embolism: comprehensive transcriptomic analysis in understanding the pathogenic mechanisms of the disease. BMC Genomics 24, 10 (2023). https://doi.org/10.1186/s12864-023-09110-0

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-023-09110-0

Keywords

  • Pulmonary embolism
  • RNA-seq
  • Molecular pathways
  • Endoplasmic reticulum

BMC Genomics

ISSN: 1471-2164

Contact us