Cardiovascular diseases (CVDs) include atherosclerosis, which is an inflammatory disease of large and medium vessels that leads to atherosclerotic plaque formation. Atherosclerosis is characterized by early endothelial cell (EC) dysfunction and altered contractility of vascular smooth muscle cells (VSMCs). Recruitment of blood leukocytes to the injured vascular endothelium characterizes the initiation and progression of atherosclerosis and involves many inflammatory mediators, modulated by the cells of both innate and adaptive immunity (Libby, 2021; Soehnlein and Libby, 2021). The key factors contributing to the early stages of atherosclerosis and plaque development include the pro-inflammatory cytokine interferon (IFN)α and IFNγ and the pattern recognition receptor (PRR) Toll-like receptor 4 (TLR4) (Szelag et al., 2016). Together, they trigger the activation of members of the signal transducer and activator of transcription (STAT) and interferon regulatory factor (IRF) family (Michalska et al., 2018; Antonczyk et al., 2019). STAT activation is mediated by a highly conserved SH2 domain, which interacts with phosphotyrosine (pTyr) motifs for specific STAT–receptor contacts and STAT dimerization. The active dimers induce gene transcription in the nucleus by binding to target genes through a DNA-binding domain (DBD) interacting with gamma-activated sequence (GAS) sites. IRFs possess a conserved DBD and an IRF-association domain (IAD) that participate in the interactions with other members of the IRF family, other transcription factors (i.e., STATs), and co-factors. IRFs bind to the IFN regulatory element (IRE) and IFN-stimulated response element (ISRE) to activate the transcription of IFN and IFN-stimulated genes (ISGs) (Szelag et al., 2016; Antonczyk et al., 2019).

In many immune cells and cells from the vasculature, IFNγ and TLR4 participate in signaling cross-talk through combinatorial actions of distinct and overlapping transcription factors on ISRE, GAS, ISRE/GAS, ISRE/NF-κB, or GAS/NF-κB binding sites (Platanitis and Decker, 2018; Piaszyk-Borychowska et al., 2019). As such, inflammation-induced activation of multiple STATs, IRFs, and NF-κB coordinates robust expression of various chemokines, adhesion molecules, and antiviral and antimicrobial proteins. Thus, signal integration between IFNγ and LPS in vascular cells and atheroma-interacting immune cells modulates important aspects of inflammation, with STATs and IRFs being important mediators. In particular, STAT1, STAT2, and STAT3 (Szelag et al., 2016; Plens-Galaska et al., 2018) and IRF1, IRF4, IRF5, IRF8, and IRF9 (Döring et al., 2012; Zhang et al., 2014a; Guo et al., 2015; Cheng et al., 2017; Liu et al., 2017; Seneviratne et al., 2017; Clément et al., 2018; Antonczyk et al., 2019; Leipner et al., 2021) have recently been recognized as prominent modulators of inflammation, especially in immune and vascular cells during atherosclerosis. Based on this, these proteins represent interesting therapeutic targets, and combined inhibition could be a novel treatment strategy in CVDs (Szelag et al., 2016; Antonczyk et al., 2019).

STAT inhibitory strategies are numerous, and by exploring the pTyr-SH2 interaction area of STAT3, searches for STAT3-targeting compounds are numerous and yielded many small molecules (Szelag et al., 2016). Recently, we developed a pipeline approach that combines comparative in silico docking of multi-million CL and CDL libraries to multiple STAT-SH2 models with in vitro STAT inhibition validation as a novel selection strategy for STAT-targeting inhibitors (Czerwoniec et al., 2015; Szelag et al., 2015). This approach allowed us to identify a new type of multi-STAT inhibitor, C01L_F03, targeting the SH2 domains of STAT1, 2, and 3 with equal affinity. Moreover, we observed a similar STAT cross-binding mechanism for STATTIC, a previously identified STAT3 inhibitor (Plens-Galaska et al., 2018). This novel class of multi-STAT inhibitors was shown to mediate the genome-wide inhibition of pro-atherogenic gene expression directed by the cooperative involvement of STATs with IRFs and/or NF-κB (Plens-Galaska et al., 2018).

To date, indirect IRF modulation has been mainly studied in terms of antiviral response regulation and cancer treatment, using, i.e., antisense oligonucleotides and siRNA knockdown strategies (Antonczyk et al., 2019). However, recently, promising small-molecule-based IRF4 (Agius et al., 2023) and cell-penetrating peptide (CPP)-based IRF5 direct-inhibition strategies were developed in connection with multiple myeloma and SLE, respectively (Banga et al., 2020; Song et al., 2020). Here, we extended our STAT inhibitor pipeline approach with 3D structure models for IRF1, 2, and 8 DBDs (Szelag et al., 2015; 2016). Using comparative in silico docking of these IRF-DBD models on a multi-million natural compound ZINC library (Irwin and Shoichet, 2005), we identified the novel multi-IRF inhibitor, ALEKSIN, which exhibited genome-wide inhibition potential toward IRF-, STAT-, and NF-κB-mediated transcription, similar to STATTIC. Furthermore, we discovered an ALEKSIN and STATTIC commonly inhibited pro-atherogenic gene signature that could help monitor plaque progression during experimental atherosclerosis. Together, this suggests that the application of a multi-IRF/STAT inhibitory strategy could offer great promise for the diagnosis and treatment of CVDs.

Three-dimensional models of DNA-binding domains of IRF1, IRF2, and IRF8 were prepared based on the existing crystal structures for IRFs deposited in RCSB Protein Data Bank. A detailed description of the utilized homology modeling procedure and optimization of the models is presented in the study by Szelag et al. (2016) and Szelag et al. (2016). For a better understanding of the interaction of IRF DNA-binding domains with their target sequence, 3D models were designed in complex with the IRE DNA (consensus sequence: 5′-GAG​AAG​TGA​AAG​T-3′). To find an “ideal” cavity, the molecular probe of the active site to which ligands are matched, “protomol” was generated (Szelag et al., 2016).

A natural compound library (NCL) containing 131,582 small molecules that are natural metabolites (Irwin and Shoichet, 2005) was selected and downloaded from the ZINC database. The compounds found in the NCL are characterized by low molecular weight, chemical parameters fulfilling the criteria of the Lipinski’s rule of five (Lipinski, 2004), and ready-to-dock parameters of protonation state and partial atomic charges (Irwin and Shoichet, 2005).

Geometries of potential IRF inhibitors used for docking were obtained from the ZINC database (code names presented in Table 2). The structures were also provided in ready-to-dock 3D formats with molecules represented in biologically relevant forms (Irwin and Shoichet, 2005).

A five-step comparative approach, CAVS, developed by our team in order to select STAT-specific compounds (Czerwoniec et al., 2015; Szelag et al., 2016; Plens-Galaska et al., 2018) was utilized to identify potential IRF inhibitors. This novel tool combines comparative in silico docking to the IRF-DBD with the in vitro validation of potential inhibitors. Structural models of IRF-binding domains and a collection of compounds from the NCL were utilized in the virtual screening procedure. As a result, we identified a number of potential IRF inhibitors for further validation and characterization. Compounds were selected using the pscreen algorithm throughout the in silico ligand–protein docking procedure (Szelag et al., 2016; Antonczyk et al., 2019). In other words, a library of small molecules was docked to the binding pockets of IRF1-, IRF2-, and IRF8-DBD of their 3D structural models. The binding affinity of the individual compounds was compared by using the binding score (BS) and comparative binding affinity value (CBAV), which also allowed us to distinguish the inhibitory potential for each of the IRF-DBDs. Based on the BS and CBAV, a list of the most promising potential IRF inhibitors was generated. After confirming the purchasability and availability of these compounds in numerous vendors, we ordered 20 of them for further in vitro testing.

In order to compare ALEKSIN and STATTIC with compounds obtained from the NCL virtual screening, docking simulations of ALEKSIN and STATTIC to IRF1, IRF2, and IRF8 DBDs were performed using the pgeom algorithm implemented in Surflex-Dock 2.6 (Jain, 2003; 2007). For each structure in the predefined area of IRF1, 2, and 8 DNA-binding domain, we obtained 20 binding poses. Then, for each compound, the best of 20 binding poses was filtered out for further analysis. Finally, IRF1 CBAV was determined to compare the binding between IRF1, IRF2, and IRF8 for both compounds.

Human microvascular endothelial cells (HMECs) (Ades et al., 1992) were provided by the Center for Disease Control and Prevention (Atlanta, GA) and cultured in MCDB-131 medium (IITD PAN, Wroclaw, Poland) containing 10% of fetal bovine serum (FBS) (Gibco, Thermo Fisher Scientific), 100 U/mL penicillin, 100 μg/mL streptomycin, 0.01 μg/mL EGF, 0.05 μM hydrocortisone, and 2 mM L-glutamine. STATTIC was purchased from Sigma and ALEKSIN (ZINC9547778, MolPort-004-931-223) from IBS, STOCK6S-36352. Recombinant IFNα and IFNγ were purchased from Merck, while LPS was provided by Sigma-Aldrich. Rabbit polyclonal antibodies against STAT1-pTyr701, tSTAT1, tSTAT2, IRF1, and IRF9 were obtained from Santa Cruz, and STAT2-pTyr689, from Merck. The Tubulin antibody was purchased from Merck, and the anti-rabbit HRP-conjugated antibody, from Sigma-Aldrich.

Depending on the experiment design, the medium was changed from complete to starving medium containing 1% FBS 8–12 h before treatment. Then, the cells were pre-treated with various combinations of inhibitory compounds and IFNα, IFNγ, and LPS. HMECs were treated with either a single stimulus, 200 U/mL (4 h for RNA isolation and 2 h for protein isolation) of IFNα, 10 ng/mL of IFNγ (8 h for RNA isolation and 4 h for protein isolation), or a combination of IFNγ for 8 h of LPS (1 μg/mL) for 4 h. Depending on the experiment, different concentrations of ALEKSIN were administered 24 h before trypsinization, while in the case of STATTIC, 8 h before. The effect of ALEKSIN on cell viability was quantified by comparing cell counts and morphology [using the ZOE Cell Imager (brightfield channel, ×20 objective)] of HMECs before and after treatment with different concentrations of the compound.

Total RNA was extracted using the GeneMATRIX Universal RNA Purification Kit (EurX, Gdansk, Poland, E3598-02) according to the manufacturer’s instructions. A measure of 500 ng of purified total RNA was then reverse-transcribed using Thermo Fisher Scientific reagents (K1622). The transcripts were quantified via qPCR with Maxima SYBR Green/ROX qPCR Master Mix (K0223, TFS) using the CFX Connect Thermal Cycler System (Bio-Rad Laboratories, Hercules, CA, United States). The target gene levels were normalized to glyceraldehyde-3-phosphate dehydrogenase (GAPDH). The PCR primers used are listed in Table 1.

Table 1. List of PCR primers.

Western blot analysis was essentially performed according to Plens-Galaska et al. (2018). HMECs were washed with phosphate buffered saline (PBS) and lysed using radio-immune precipitation assay (RIPA) lysis buffer (50 mM Tris-HCl, pH = 8.0, 150 mM NaCl, 1% Nonidet-40, 0.5% sodium deoxycholate, 0.1% SDS, 1% protein inhibitor cocktail, 1% EDTA, and 0.1% PMSF) and stored at −80°C. Lysates were quantified using a bicinchoninic acid (BCA) kit (Pierce). A measure of 30 µg of protein was loaded on Blot 4%–12% Bis-Tris Plus Gels, electrophoresed, and transferred to PVDF membranes (Santa Cruz). All Western blot analyses were performed using the SNAP i.d. system (Merck). Membranes were blocked in 0.125% non-fat dry milk or 1% BSA in TBS-Tween (TBS-T) and incubated with primary antibodies (1:500 IRF1, 1:500 IRF9, 1:1,000 pSTAT1, 1:500 tSTAT1, 1:500 pSTAT2, 1:500 tSTAT2, and 1:2,000 tubulin) and then with the secondary anti-rabbit HRP-conjugated antibody (1:20,000). Immunoreactive bands were visualized by enhanced chemiluminescence (ECL) using the Luminata Forte HRP Substrate (Merck) and detected using the G:Box System (Syngene). After detection, the membranes were stripped with a buffer containing 25 mM glycine and 1% SDS, pH = 2.0, and re-probed. ImageJ software (https://imagej.net/ij/) was used for Western blot quantification, with α-tubulin as the reference protein.

The scratch assay was performed according to Plens-Galaska et al. (2018), with minor changes. HMECs were seeded at a density of 400.000/mL on 6-well plates and cultured until they reached around 85%–90% confluency. The cells were pre-treated with ALEKSIN for 24 h and STATTIC for 8 h. After 12 h of treatment with compounds, scratches with a diameter of approximately 9.6 mm were introduced and subsequently treated with or without 10 ng/mL of IFNγ and 1 μg/mL of LPS. At the same time, reference points were generated, and the first image was taken. The second image was taken after 12 h using the AxioObserver Z1 Microscope (Zeiss). The images acquired for each sample from two independent repeats were further analyzed quantitatively using ImageJ to determine the % wound coverage and the rate of cell migration (Rm) in μm/h (Liang et al., 2007; Schneider et al., 2012).

1. The animals were divided into two groups (×2 n = 2). Cross sections of the mouse aorta (left part of the aortic arch) were formalin-fixed, dehydrated with ethanol and xylene, and paraffin-embedded. The paraffin sections were stained with hematoxylin and eosin (H&E) using the Leica Autostainer XL H&E Slide Stainer. The sections were deparaffinized in xylene and rehydrated with distilled water. H&E staining was performed by 3 min of submersion in hematoxylin, followed by a washing step in tap water and 30-s submersion in eosin, after which the slides were dehydrated and covered with a cover slip. Histological staining and immunohistochemistry were assessed under a Nikon Eclipse Ti microscope. Images were taken and processed using NIS-Elements software (Nikon). The slides were observed under a microscope and scored for the presence of early or progressive lesions, as described by Shibata et al. (2017).

2. The animals were divided into two groups (×2 n = 2). Oil Red O staining of the whole aorta was performed as described earlier (Centa et al., 2019), with minor changes. In brief, aortas were dissected, pinned on a hard surface, stained with 2 mg/mL Oil Red O dilution for 20 min, and rinsed twice in isopropanol. Images were taken using an Olympus ×10 microscope, and the lesion area was assessed as the percentage of total area of the aorta using ImageJ software.

The RNA-seq library was essentially prepared according to Sekrecka et al. (2023) and Sekrecka et al. (2023). RNA from HMECs and mouse aortas was quantified using a Qubit RNA BR assay kit (Q10210; Thermo Fisher Scientific), and the quality was assessed via the Agilent 2100 Bioanalyzer using the RNA 6000 Nano kit (5067–1511, Agilent Technologies, Santa Clara, CA, United States), according to the protocols provided by manufacturers. RNA degradation was assessed using the RNA integrity number (RIN), and samples with a RIN higher than 9 were then used for further analysis. RNA libraries were prepared in three biological repeats from 1 μg of total RNA using the NEBNext® Ultra™ or Ultra™ II RNA Library Prep Kit for Illumina® (New England Biolabs [NEB], Ipswich, CA, United States) together with the NEBNext Poly(A) mRNA Magnetic Isolation Module (NEB) and NEBNext® Multiplex Oligos for Illumina® (NEB), according to the manufacturer’s protocol. The quality and fragment distribution of the prepared libraries were estimated using the Agilent High-Sensitivity DNA kit (5067–4626, Agilent Technologies), and the quantity was assessed using the Qubit dsDNA HS assay kit (Q32851, Thermo Fisher Scientific).

Sequencing was performed on the Illumina NextSeq500 platform (75 bp, single-end). After quality control, adapter trimming, and quality filtering using fastp (v.0.22.0) (Chen S. et al., 2018), the reads were aligned to the human GRCh38/hg38 genome using STAR 2.6.1d (Dobin et al., 2013). Raw per-gene counts were calculated using featureCounts (v.1.6.2) (Liao et al., 2014). Before differential expression testing, the genes were prefiltered to include only genes that had a minimum of 10 raw reads in at least one sample. Count normalization and differential gene expression (DEG) analysis were performed using the edgeR package (v.3.30.3) (Robinson et al., 2009) in R version 3.6.3 (R Core Team (2021) R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna; https://www.R-project.org). Tests were corrected for multiple comparisons using the method of Benjamini and Hochberg (BH). Genes with a log2 fold-change ≥1 and adjusted p-value ≤0.05 were considered differentially expressed between conditions.

Libraries were sequenced on the Illumina HiSeq X Ten platform (150 bp, paired-end). Sequence reads were trimmed to remove possible adapter sequences and nucleotides with poor quality using fastp (v.0.22.0) with the default parameters and mapped to the M. musculus reference genome (GRCm38/mm10) using STAR aligner (2.6.1d). Unique per-gene counts were calculated using STAR. Before differential expression testing, genes were prefiltered to include only genes that had a minimum of 10 raw reads in at least one sample. Counts were normalized as logarithmic counts per million. Data analyses were carried out using R (v.4.2.2) packages. The linear regression model package limma (v.3.52.4) (Ritchie et al., 2015) was used to identify DEGs between the normal-diet group and the high-fat diet group. Tests were corrected for multiple comparisons using the BH method. Genes with a log2 fold-change >0.5 and adjusted p-value <0.05 were considered differentially expressed between conditions. All raw and processed sequencing data are accessible via the NCBI Gene Expression Omnibus (GEO) under the accession numbers GSE270277 (human HMEC RNA-seq data) and GSE270260 (mouse ApoEKO RNA-seq data).

GO Enrichment Analysis was performed with the “enrichGO” function of the clusterProfiler package (v.4.6.0) (Wu et al., 2021), using the biological process ontology. Enrichment p-values were corrected for multiple comparisons using the BH method (cutoff for p-value = 0.01 and for qvalue = 0.05). All DEGs were used as background. The “simplify” function was applied to reduce the redundancy of the enriched GO terms, with default parameters and cutoff 0.6. A total of 15 terms with the highest statistical significance were used for visualization as bar plots (ggplot2 3.5.0) (Wickham, 2011). Enrichment was defined as −log10 (adjusted p-value).

For promoter analysis, active promoters for genes differentially expressed between conditions were selected using the proActiv R package (v.1.8.0) (Demircioğlu et al., 2019) with default parameters. In brief, counts and normalized promoter activity estimates for each annotated promoter were generated based on junction files obtained from STAR during alignment. Only promoters with high and medium average activities for each gene (categorized as major and minor, respectively) across the samples were used in further analysis. Promoters with average activities less than 0.25 were considered inactive. proActiv does not provide activity estimates for promoters that are not uniquely identifiable from splice junctions (single-exon transcripts and promoters that overlap with internal exons). In order to obtain these promoters, we used cap analysis of gene expression (CAGE) data tag sequencing of the 5′ end of transcripts from the FANTOM5 project (hg38_fair + new_CAGE_peaks_phase1and2. bed.gz, mm10_fair+new_CAGE_peaks_phase1and2.bed.gz). Only promoters with the highest activity (maximum peak score value per gene) were used.

Subsequently, combined lists of promoters for all DEGs for HMECs and ApoEKO data were prepared using the custom R script and used for the identification of enriched transcription factor-binding sites in HOMER v.4.11 (annotatePeaks.pl) (Heinz et al., 2010). For the identification of GAS and ISRE-binding sites, the set of selected matrices (four for GAS and three for ISRE) from the study by Sekrecka et al. (2023) was used. A similar strategy was applied for NFkB matrix selection and optimization, using two matrices from the HOMER database [Homer Motif 235 (GSE23622) and Homer Motif 268 (GSE19485)] and one motif from the de novo analysis. Known motif search was performed on sequences −850/+150 around the transcription start site (TSS) with gene annotation files from GENCODE (v.43 and v.M25 for human and mouse data, respectively).

Graphical representations of the results were generated using VennDiagram v.1.7.3 (Chen and Boutros, 2011), pheatmap v1.0.12 (Kolde, 2019), and ComplexHeatmap v.2.6.2 (Gu et al., 2016) packages.

Heatmaps for selected genes were created using normalized counts. Hierarchical clustering (by row, with default method: complete, Euclidean distance) was used to select groups of IFNγ+LPS target genes inhibited by ALEKSIN, STATTIC, or both in HMECs. For plotting, row scaling with Z-scores was performed. The color scale indicates the expression change over time for each sample compared to the expression of the control. Colors represent high (red) and low (blue) normalized intensity, respectively.

Single-cell RNA-seq data from mice with deficient low-density lipoprotein receptor and expressing only apolipoprotein B100 (Ldlr−/−/Apob100/100) under a 3-month high-fat diet were obtained from the study by Örd et al. (2023). The data were then processed and analyzed, as previously described (Boroujeni et al., 2024). The processed data were used to generate heatmaps using ComplexHeatmap package version 2.20.0 (Gu et al., 2016). For clustering the rows, the Euclidean method was implemented.

To identify novel multi IRF-DBD inhibitory compounds, an NCL from the ZINC database was screened using the pre-screen algorithm (see Materials and Methods). Compounds were docked to IRF1-, IRF2-, and IRF8-DBD. The compounds were first selected according to their IRF1-, IRF2-, and IRF8-BS and subsequently filtered by IRF1-CBAV(IRF2)∼0 and IRF1-CBAV(IRF8)∼0 for comparable binding affinities. Accordingly, 20 compounds with the highest BS (named pI01–pI20) were purchased (Table 2). To test the inhibitory capacity of these compounds toward IRF target gene expression in vitro, HMECs were treated with 50 μM of each compound for 24 h and 10 ng/mL of IFNɣ for 8 h. Except for pI05, pI11, and pI13 (Figure 1A), none of these compounds inhibited IRF target gene expression (not shown). Interestingly, pI05, pI11, and pI13 exhibited different inhibition patterns, from full (pI13) or partial (pI05 and pI11) inhibition of Irf1 to full (pI11 and pI13) or partial (pI05) IRF1 target gene expression (Ifit3 and Isg15) (Figure 1A). Next, we examined the in silico binding of pI05, pI11, and pI13 to the DBD of IRF1, IRF2, and IRF8. The top-scored binding conformations were visualized using PyMOL (see Materials and Methods), showing similar IRF1-, IRF2-, and IRF8-BS (Figure 1B). Together, this suggested that pI05, pI11, and pI13 inhibit IFNγ-induced IRF target gene expression by targeting the DBD of multiple IRFs. Because of problems with the stability of pI05 and purchasability of pI13, we only continued with the further characterization of pI11 and named it ALEKSIN.

Table 2. List of 20 compounds from the screening of the Natural Compound Library chosen for in vitro testing.

Figure 1. (A) pI05, pI11, and pI13 inhibit IFNɣ-induced expression of IRF1, IFIT1, and ISG15. HMECs were treated with 50 μM of each compound for 24 h and 10 ng/mL of IFNɣ for 8 h. RNA was isolated and analyzed by qPCR. Experiments were performed in two individual repeats, which were compared by a t-test. *p < 0.05, **p < 0.01, and ***p < 0.001. (B) Top scored binding conformations of pI05, pI11, and pI13 in the DBD of IRF1, IRF2, and IRF8. Graphical representation describes in detail the binding mode of the top scored conformation of the inhibitor in the active pocket of the IRF-DBD. Models presented as human IRF-DBD/IRE complexes are shown in cartoon representation. Proteins are presented with a visible secondary structure, alpha helices, and beta sheets, with hot pink denoting IRF1-DBD, tv-blue denoting IRF2, and aquamarine cyan denoting IRF8. The dsDNA fragment of the respective IRF-DBD/IRE complexes implicates the position of the selected target cavity for the inhibitory compound. The target pocket for virtual screening with an idealized active-site ligand (protomol) in the DBD of hIRF1, hIRF2, and hIRF8 is shown in a transparent surface representation in cyan blue. Protomol is based on the interaction plane between DNA and amino acid residues of the respective hIRF-DBD/IRE complexes. Docking of pI05, pI11, and pI13 to the DBD of IRF1, IRF2, and IRF8. Top-scored binding conformations (with the highest BS) of pI05, pI11, and pI13 in the IRF1, 2, and 8 DBD shown in stick representation and colored according to the atomic structure. Docking simulations were performed using Surflex-Dock 2.6 and were visualized in PyMOL.

To study the IRF inhibitory characteristics of ALEKSIN in more detail, we tested it first on IFNγ-treated HMECs at different concentrations and time points. Apparently, treatment for 12 h with ALEKSIN at concentrations varying from 5 to 20 μM did not result in the inhibition of IFNγ-induced IRF target gene expression (not shown). On the other hand, exposure of cells to ALEKSIN for 24 h at 10 or 20 μM completely inhibited IIFNα and IFNγ-induced expression of STAT1, IFIT3, and ISG15, whereas the effect on IRF1 expression was only partial (Figures 2A, B).

Figure 2. (A, B) ALEKSIN inhibits IFNα and IFNy-induced gene expression of IRF1 and IRF1 target genes Stat1, Ifit3, and Isg15 in a concentration-dependent manner. HMECs were pre-treated with 10 μM or 20 μM of ALEKSIN for 24 h and (A) 200 U/mL of IFNα for 4 h or (B) 10 ng/mL of IFNy for 8 h. RNA was isolated and subjected to qPCR analysis. Experiments were performed in two individual repeats and compared by a t-test. *p < 0.05, **p < 0.01, and ***p < 0.001. (C, D) ALEKSIN inhibits IFNα and IFNy-induced protein expression of STAT1, STAT2, IRF1, and IRF9 and partial phosphorylation of STAT1 and STAT2. HMECs were treated with 10 μM of ALEKSIN for 24 h and (C) 200 U/mL of IFNα or (D) 10 ng/mL of IFNy for 0, 2, 4, 8, and 24 h. Protein extracts were collected, and levels of IRF1, IRF9, pSTAT1, pSTAT2, tSTAT1, tSTAT2, and α-tubulin were assessed by Western blotting.

We also tested the effect of ALEKSIN on STAT phosphorylation and STAT and IRF protein expression. Interestingly, pre-treatment of HMECs with ALEKSIN, followed by IIFNα exposure, resulted in the partial inhibition of STAT1 and STAT2 phosphorylation and expression of IRF1 at early time points (2 and 4 h). IRF9, STAT1, and STAT2 expression was also inhibited but with a more delayed pattern (between 4 and 24 h of IIFNα treatment) (Figures 2C; Supplementary Figures S1A, S3). The same was true for IFNγ-treated cells, with comparable effects of ALEKSIN on expression of IRF1, IRF9, STAT1, and STAT2 in a time-dependent manner. The inhibition of STAT1 phosphorylation, on the other hand, was not so pronounced (Figures 2D; Supplementary Figures S1B, S3). However, using ALEKSIN at a higher concentration (25 μM) also resulted in the partial inhibition of IFNγ-induced STAT1 phosphorylation (Figure 3A; Supplementary Figures S2A, S3). Compared with ALEKSIN (10 μM), under these conditions, STATTIC (10 μM) completely inhibited STAT1 phosphorylation (Figure 3B; Supplementary Figures S2B, S3), which is in agreement with its recent identification as a potent multi-STAT inhibitor (Plens-Galaska et al., 2018). Moreover, ALEKSIN showed no toxicity under these conditions based on RNA (Supplementary Figure S4A) and protein (Supplementary Figure S4B) concentration stability, and cell viability (Supplementary Figure S4C) and morphology (Supplementary Figure S4D).

Figure 3. (A, B) ALEKSIN and STATTIC inhibit IFNy-induced expression of IRF1 and IRF9 and phosphorylation of STAT1. HMECs were pre-treated with 20 μM (A) or 10 μM (B) of ALEKSIN for 24 h and 10 μM (B) of STATTIC for 8 h and 10 ng/mL of IFNγ for 2 h. Protein extracts were collected, and levels of IRF1, IRF9, pSTAT1, and α-tubulin were assessed by Western blotting. (C) Docking of ALEKSIN to the SH2 domain of STAT1, STAT2, and STAT3. Graphical representation describes in detail the binding mode of top scored conformation of the inhibitor in the active pocket of the STATs-SH2 domain. Models are presented in cartoon representation with a visible secondary structure, alpha helices, and beta sheets, in yellow for STAT1, tv-orange for STAT2, and deep purple for STAT3. The target pocket for virtual screening with an idealized active-site ligand (protomol) in the SH2 domain of STAT1, STAT2, and STAT3 is shown in transparent surface representation in cyan blue. Top-scored binding conformations of ALEKSIN in the STAT1, 2, and 3-SH2 domains are shown in stick representation (with the highest BS) and colored according to the atomic structure. Docking simulations were performed using Surflex-Dock 2.6 program and were visualized in PyMOL.

At the same time, we docked ALEKSIN to the SH2 domain of STAT1, STAT2, and STAT3 and observed that it interacted with the selected binding cavity but with a much lower binding affinity than STATTIC (Figure 3C). This low-affinity binding of ALEKSIN to STAT-SH2, in addition to high-affinity IRF-DBD binding (Figure 1B), could reflect the presence of IRF-independent effects and aligns with the observed partial inhibition mediated by ALEKSIN toward IFN-induced STAT phosphorylation.

Subsequently, we further investigated the ability of ALEKSIN and STATTIC to inhibit pro-inflammatory and pro-atherogenic signaling communicated by IFNγ and LPS cross-talk. As shown in Figure 4, pre-treatment of HMECs with ALEKSIN or STATTIC resulted in the inhibition of IFNγ+LPS-induced gene expression of IRF1, STAT1, IFIT3, ISG15, and MX1. In general, STATTIC was slightly more potent than ALEKSIN. These data suggested that ALEKSIN and STATTIC may commonly block STAT, IRF, and NF-κB cooperative promotor activation mediated by IFNγ and LPS in human microvascular endothelial cells. To provide further evidence for this, we studied the genome-wide effects of ALEKSIN and STATTIC on IFNγ+LPS-mediated vascular inflammation. For this, we performed RNA-seq on RNA isolated from HMECs treated with IFNγ+LPS in the presence or absence of 10 µM of ALEKSIN or 10 µM of STATTIC (GEO accession: GSE270277). IFNγ+LPS increased the expression of 537 genes by at least two-fold or higher than untreated cells, of which the top 25 are shown in Table 3. These included many recognized IFNγ and LPS target genes associated with chemotaxis/migration (CXCL9, CXCL10, CXCL11, CX3CL1, CCL8, and CCL20) and immune response (GBP4, GBP5, GBP7, and IDO1).

Figure 4. ALEKSIN and STATTIC inhibit IFNγ+LPS-induced gene expression of IRF1, STAT1, IFIT3, ISG15, and MX1. HMECs were treated with (A) 10 μM of ALEKSIN for 24 h or (B) 10 μM of STATTIC for 8 h and for 8 h with IFNγ +4 h with LPS. RNA was isolated and subjected to qPCR analysis. Experiments were performed in two individual repeats, which were compared by the t-test. *p < 0.05, **p < 0.01, and ***p < 0.001.

Table 3. Top 25 upregulated genes in HMECs in response to IFNγ and LPS treatment.

Next, we identified 405 IFNγ+LPS target genes that were commonly inhibited by ALEKSIN and STATTIC (Figure 5A), with the inhibition pattern of the top 25 genes shown in Table 3. Among them, we could recognize a variety of STAT and IRF target genes. The inhibition ratio for ALEKSIN was determined by dividing [FC IFNγ/LPS vs. UN] over [FC IFNγ/LPS + ALEKSIN vs. UN], and for STATTIC, it is [FC IFNγ/LPS vs. UN] over [FC IFNγ/LPS + ALEKSIN vs. UN] (Table 3). From this inhibition ratio, it can be concluded that STATTIC is more potent than ALEKSIN (Table 3). In addition, we also recognized 54 genes that were predominantly inhibited by ALEKSIN, whereas for another 58 genes, STATTIC was the more dominant inhibitor (Figure 5A). The complete list of upregulated genes in response to IFNγ+LPS in the presence or absence of ALEKSIN or STATTIC is shown in Supplementary Table S1.

Figure 5. ALEKSIN and STATTIC commonly inhibit the cross-talk between IFNγ and LPS-induced genes. Genome-wide analysis of ALEKSIN and STATTIC commonly and specifically inhibited IFNγ+LPS target genes. (A) ALEKSIN and STATTIC commonly (dark green) and specifically (light green: ALEKSIN; blue: STATTIC) inhibited IFNγ+LPS target genes were identified using hierarchical clustering. (B) Gene clusters identified under (A) were subjected to GO analysis using clusterProfiler. A total of 15 terms with the highest statistical significance were used for visualization as bar plots. Enrichment was defined as −log10 (adjusted p-value).

GO analysis of the 405 commonly inhibited genes further revealed the enrichment of general biological terms connected to inflammation and atherogenesis, including cytokine-mediated signaling pathway, defense response and immune system process, regulation of cytokine production, inflammatory response, innate immune response, adaptive immune response, T-cell activation, regulation of leukocyte proliferation, leukocyte cell–cell adhesion, neutrophil chemotaxis, cellular response to lipids, and response to the tumor necrosis factor (Figure 5B). GO analysis of the ALEKSIN- or STATTIC-specific gene clusters, on the other hand, recognized similar but more restricted biological terms. For example, ALEKSIN-specific genes were associated with inflammatory response, regulation of cytokine production, fat-cell differentiation, lymphocyte and leukocyte differentiation, response to lipids, regulation of cell–cell adhesion, and response to lipopolysaccharide (Figure 5B). In contrast, STATTIC-specific genes were connected to the cytokine-mediated signaling pathway, defense response and immune system process, regulation of cytokine production, inflammatory response, innate immune response, adaptive immune response, leukocyte proliferation, leukocyte cell–cell adhesion, and myeloid cell differentiation (Figure 5B). This suggests that, in general, functional overlap exists between ALEKSIN and STATTIC common and specific genes.

We subsequently performed promoter analysis on the ALEKSIN and STATTIC common and specific gene clusters. In the region −850 to +150 bp around active promoter sites detected using the proActiv R package (v.1.8.0), we identified transcription factor-binding sites using the HOMER v.4.11 tool and a set of selected and optimized matrices for GAS, ISRE, and NF-kB-binding site motifs (for details see Material and Methods).Figure 6A shows the predicted representation of individual or combined ISRE, GAS, or NF-κB-binding sites in the proximal promoters of ALEKSIN and STATTIC commonly inhibited genes. Accordingly, the majority of these genes contained single ISRE (16.5%) or GAS (19.4%) sites or combinations of ISRE + GAS (24.4%), ISRE + NF-κB (5.29%), GAS + NF-κB (13.5%), or ISRE + GAS + NF-κB (15.3%). In general, under these conditions, ISRE motifs correspond to the potential binding of multiple STATs (STAT1 and STAT2), IRFs (IRF1, IRF7, IRF8, and IRF9), and GAS motifs to that of multiple STATs (STAT1 and STAT3). Surprisingly, 19 genes (5.59%) were assigned to the group with only an NF-κB-binding site in their proximal promoter. However, these included genes like BCL3, CACNA1A, CX3CL1, CXCL1, CXCL2, CXCL3, CXCL6, GBP2, GBP7, IL7R, KRT17, NOD2, NUB1, and OPTN, which contained putative GAS and/or ISRE sequences outside the 850-bp selected promoter area (not shown).

Figure 6. ALEKSIN and STATTIC common and specific genes share the presence of IRF, STAT, and/or NF-κB-binding sites. Gene clusters of ALEKSIN and STATTIC commonly and specifically inhibited IFNγ+LPS target genes were subjected to promoter analysis using HOMER, and the distribution of GAS, ISRE, and NF-kB sites for ALEKSIN and STATTIC common (A), ALEKSIN-specific (B), and STATTIC-specific (C) is shown in individual Venn diagrams.

Analyzing the promoters of the 58 STATTIC-specific (Figure 6B) and 54 ALEKSIN-specific (Figure 6C) genes displayed a similar distribution for single ISRE (12.5% vs. 13.5%), GAS (29.2% vs. 29,7%), or NF-κB (8.33% vs. 10.8) sites, or combinations of ISRE + GAS (10.4% vs. 21.6%), GAS + NF-κB (18.8% vs. 10.8%), or ISRE + GAS + NF-κB (12.5% vs. 13.5%). On the other hand, ISRE + NF-κB sites could only be recognized among STATTIC-specific genes (8.33%).

Nevertheless, in general, these results strongly suggest that ALEKSIN and STATTIC common and specific genes share the presence of IRF, STAT, and/or NF-κB-binding sites and predict that ALEKSIN and STATTIC commonly inhibit pro-inflammatory and pro-atherogenic gene expression directed by the cooperative involvement of STATs, IRFs, and/or NF-κB.

In addition, we aimed at providing evidence that an IRF/STAT-dependent inhibitory strategy could be used to block IFNγ+LPS-induced vascular inflammation. Thus, we performed a wound healing assay to examine the effect of ALEKSIN on IFNγ+LPS-induced EC migration compared to STATTIC (Figure 7). Cells stimulated with IFNγ+LPS displayed increased capacity of migration, resulting in >80% of wound coverage after 12 h of treatment (Figure 7). In contrast, HMECs treated additionally with ALEKSIN or STATTIC demonstrated a drastic reduction in migratory activity. Both inhibitors caused a decrease in the IFNγ+LPS-induced wound healing capacity to less than 10% (Figure 7) compared to 25% in the absence of IFNγ+LPS (Figure 7).

Figure 7. ALEKSIN inhibits IFNγ+LPS-induced EC migration similar to STATTIC. Wound healing assay performed on HMECs treated with ALEKSIN or STATTIC with or without IFNɣ+LPS (A). Border lines (blue) determine scratch borders at the beginning (left) and end (right) of the experiment. Statistical evaluation of wound healing assay (B). Graph shows the percentage of healed wound compared to 0-h control (left) and the rate of cell migration (Rm) in μm/h (right). The experiment was performed in two individual repeats, which were compared by the two-way ANOVA test and Bonferroni correction. *p < 0.05, **p < 0.01, and ***p < 0.001.

To characterize the expression of ALEKSIN and STATTIC commonly inhibited genes during atherosclerotic plaque formation, we used the ApoEKO HFD mouse model (Sanz-Garcia et al., 2017). As shown in Figure 8, a HFD for 12 weeks resulted in the development of aortic atherosclerotic plaques (aortic arch: Figure 8B; whole aorta: Supplementary Figure S5) compared to LFD (aortic arch: Figure 8A; whole aorta: Supplementary Figure S5), which correlated with increased levels of total and LDL cholesterol (Figures 8C, D). RNA-seq on aortic arch RNA isolated from ApoEKO mice on HFD (n = 8) vs. LFD (n = 8) (GEO accession: GSE270260) identified 763 HFD-upregulated genes (log2FC > 0.5; Supplementary Table S2). GO analysis of these genes revealed the enrichment of biological functions mainly involved in cytokine production, the cytokine-mediated signaling pathway, immune response-regulated signaling pathway, leukocyte-mediated immunity, leukocyte proliferation, leukocyte activation, leukocyte migration, leukocyte cell–cell adhesion, T-cell activation, and cell killing (Figure 8E). Interestingly, a comparison with the enriched biological terms from IFNγ+LPS-treated HMEC revealed a strong overlap (Figure 5B). We also performed promoter analysis, with Figure 8F showing the predicted representation of individual or combined ISRE, STAT, or NF-κB binding sites in their proximal promoters (−850 to +150). The majority of these genes contained single ISRE (15.6%), GAS (36.7%), or NF-κB (7.72%) sites, or combinations of ISRE + GAS (20.3%), ISRE + NF-κB (3.01%), GAS + NF-κB (9.6%), or ISRE + GAS + NF-κB (6.97%). This is similar to the binding site distribution seen for IFNγ+LPS-induced genes in HMECs (Figure 5B) and predicts a combined role of IRFs with STATs and/or NF-κB in experimental atherosclerosis as well. The complete list of HFD-upregulated genes is shown in Supplementary Table S2.