Patients and donors
Demographic and clinical characteristics BM-MSCs donors and of the RA patients are shown in Additional file 1: Table S1 and Additional file 2: Table S2, respectively.
Bone marrow mesenchymal stem cells
Mesenchymal stem cells (MSCs) were obtained from bone marrow (BM) aspirates collected from the iliac crests of three donors, following informed consent. We included subjects older than 18 years old, with no previous diagnosis of autoimmune disease or lymphoproliferative/neoplastic conditions. Briefly, aspirates were diluted in an equal volume with saline and centrifuged over a Ficoll layer at 2000×g for 20 min. Cellular fraction recovered was washed two times in Dulbecco’s modified Eagles medium (DMEM) (Lonza). The cell pellet obtained was suspended in 5 ml with complete culture medium (DMEM supplemented with 2 mM glutamine, 0.06% penicillin, 0.02% streptomycin, and 10% FBS). Cultures were incubated at 37 °C in a 5% CO2 humidified atmosphere in 25-cm2 flasks. After several days, non-adherent cells were removed and fresh medium was added. The medium was exchanged every 4 days of culture. When cultured cells reached 80–90% confluence, adherent cells were trypsinized (0.05% trypsin/1.0 mM EDTA), harvested, and expanded in 25-cm2 flasks. MSC characterization was performed according to the minimal criteria recommended by the ISCT (International Society for Cellular Therapy) described by Dominici et al [16]. Cells in the fourth passage were used in the experimental analysis.
Bone marrow aspirates and blood were obtained in accordance with Good Clinical Practices and the principles expressed in the Declaration of Helsinki. The study was approved by our institutional Ethics Committee (Comité Ético de Investigación Clínica Hospital Clínico San Carlos—Madrid).
Isolation of PBMCs from RA patients
Five consecutive patients diagnosed with RA according to the 2010 ACR/EULAR criteria, attending the Hospital Rheumatology Outpatient Clinic, were included in this study. Patients were over the age of 18 at disease onset and had no previous history of any other chronic disease such as diabetes, chronic kidney disease, and/or lymphoproliferative/neoplastic conditions. We excluded from the study those patients receiving more than 10 mg/day of prednisone or equivalent, those who had received any intramuscular dose of corticosteroids in the previous 2 months, or those treated with drugs that can affect lymphocyte activation, such as calcium antagonists or statins. At inclusion, demographic and clinical data were collected, and a fasting venous blood sample was extracted on EDTA as anticoagulant. PBMCs were separated by centrifugation on a Ficoll-Hypaque gradient at 900×g, for 20 min at 25 °C.
In vitro BM-MSC–PBMC co-cultures
BM-MSCs (2 × 105 cells) were cultured alone for 24 h in non-treated Falcon® 6-wells flat bottom plates (Corning) in complete low-glucose DMEM (Lonza Group Ltd., Basel, Switzerland) at 37 °C and 5% CO2. After this time, the medium was replaced by supplemented RPMI. PBMCs (2 × 106 cells) were added to the wells in ratio 1:10 (BM-MSC:PBMCs) mostly based on cellular and cellular membrane size.
In some wells, anti-CD3-/anti-CD28-coated beads (Dynabeads® Human T-Activator; Life Technologies) were also added at a ratio of 1 bead per 4 PBMCs.
Three days after co-culturing, supernatant was collected removing the anti-CD3-/anti-CD28-coated beads as well as PBMCs not attached to the BM-MSCs. Furthermore, those weakly attached to the surface of the BM-MSCs were collected by gently pipetting the bottom of each plate with clean RPMI. CD3/CD28 beads were magnetically removed (following the manufacturer’s protocol), and PBMCs were recovered after centrifugation and counted, and their viability was assessed by Trypan blue staining. BM-MSCs were treated with Trypsin/EDTA and the remaining PBMCs attached were removed by using anti-CD45 antibodies conjugated to paramagnetic microbeads (Miltenyi Biotec, Spain). Purified BM-MSCs were stored in RNA protect solution (Qiagen Iberia, Madrid, Spain) at − 80 °C.
RNA extraction and processing
Total RNA from BM-MSCs was extracted using the RNeasy Mini Kit (Qiagen Iberia, Madrid, Spain), according to the manufacturer’s protocol. RNA concentration was quantified spectrophotometrically (Nanodrop ND-1000, Wilmington, DE), and its quality and integrity assessed by capillary electrophoresis on an Agilent 2100 Bioanalyzer (Agilent Technologies. Spain). Barcoded cDNA libraries were prepared from poly(A) enriched mRNA using NebNext Ultra Directional RNA Library Prep Kit (New England Biolabs, Spain). Pooled libraries were sequenced on an Illumina NextSeq 500 instrument to generate on average 32.9 million single-end reads of 76 bp length. To avoid a batch effect bias, all samples were run simultaneously twice, and results merged.
Assessment of PBMC RNA contamination
In order to assess the degree of RNA contamination from the PBMCs, different approaches were followed (see Additional file 4).
Bioinformatic analyses
Raw sequence quality control was performed using FastQC [17] (Babraham Bioinformatics). The raw sequence reads (FASTQ format) were aligned to the cDNA sequences of the human GRCh37 reference assembly available in UCSC Genome Browser (http://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/hg19.2bit), using the Rsubread Bioconductor package v1.20.3 [18], using the default settings, and reporting only uniquely mapped reads.
Read summarization was performed with featureCounts [19], using the in-built gene annotations from the NCBI RefSeq for Hg19, included in the Rsubread package. Mapped reads for each sample were summarized into the meta-feature “gene,” thus obtaining gene-level expression counts that were used as input for gene expression analysis.
In order to improve the statistical power by decreasing the number of multiple comparisons to adjust for and to reduce the possible bias of very small counts with no biological significance, we removed those genes with low expression according to the number of counts for million mapped reads (CPMs). We set a threshold of at least 1 count per million (CPM), in at least 4 of the 5 active and/or 4 of the 5 resting samples.
To adjust for variable sequencing depths between samples, the raw gene counts were normalized using a weighted trimmed mean of the log expression ratios (Trimmed Mean of M values [TMM] algorithm) as implemented in the edgeR Bioconductor package [20, 21]. A multidimensional scaling (MDS) plot was used as an unsupervised approach to visualize the data structure of the analyzed samples.
Statistical analysis
The edgeR Bioconductor package [20] was used to identify genes that were differentially expressed between BM-MSCs co-cultured with activated or resting PBMCs. Our study design had three experimental factors: the BM-MSC donor (three healthy donors), the PBMCs donors (five RA patients), and the PBMC activation state (two levels: resting or activated). Therefore, this comparison (exposure to resting or activated PBMCs from RA patients) was nested by the RA patient from whom the PBMCs were obtained and, in turn, nested by the BM-MSCs donor (Fig. 1).
Based on this design, we used the following formula:
$$ \sim \mathrm{MSC}\_\mathrm{ID}+\mathrm{MSC}\_\mathrm{ID}:\mathrm{RA}\_\mathrm{ID}+\mathrm{MSC}\_\mathrm{ID}:\mathrm{Activ} $$
Due to the unbalanced design, we manually edited the matrix generated by this formula, removing those columns containing no information. Differential expression was analyzed using the quasi-likelihood F test [20]. Significant differentially expressed genes (DEGs) were defined as those with a log2 fold change ≥ 2 and a false discovery rate (FDR) of ≤ 1% (adjusted with the Benjamini-Hochberg method [22]).
Quality of RNA extraction
The mean (SD) RNA concentration was 255.9 (147.4) ng/μl for BM-MSCs co-cultured with resting PBMCs and 78.7 (28.5) ng/μl for BM-MSCs co-cultured with activated PBMCs. We used 250 ng of RNA from each sample. All samples had a RNA integrity number of 10.
After sequencing, we generated a total number of unique reads mapped to the human genome between 29.6 and 34.9 million, with a mean (SD) across samples of 32.9 (1.7) millions. These reads were summarized into gene-level expression counts, resulting in a mean (SD) of 23.8 (1.5) million successfully assigned reads [a mean (SD) percentage of 72.2 (1.7)] (Table 3). The difference in the number of reads between BM-MSCs co-cultured with activated or resting PBMCs was not statistically significant (Student’s t test, p = 0.82), nor was the difference in the number of reads among BM-MSC samples regarding the donor (ANOVA p-value = 0.66).
Initial summarization of reads into genes revealed 25,702 metafeatures. Count data was filtered based on the number of CPMs in that particular sample. We set a CPM threshold of 1, which represented a minimum gene count between 22 and 26, depending on the library size. After filtering, 12,821 genes remained. We further removed those genes lacking Gene Symbol identifier (n = 219). Therefore, 12,602 genes were analyzed for differential expression.
The MDS plot clearly distinguishes the transcriptional profiles of those BM-MSCs exposed to resting or activated PBMCs, separated along the first dimension. Furthermore, along the second dimension, a difference between the BM-MSCs from the first donor and those from donors 2 and 3 was observed, regardless the BM-MSCs had been exposed to active or resting PBMCs.
Gene Ontology analysis
To understand the biological impact of the gene expression changes, we performed functional enrichment analysis. Considering that for RNA sequencing data, gene length and read count can introduce biases in the gene ontology (GO) enrichment analysis, we used GOseq [23] in order to minimize this bias. We manually introduced the gene length, based on the data from the inbuilt NCBI RefSeq Hg19 annotation from the Rsubread package. We analyzed separately those genes up- or downregulated, considering significantly enriched those terms with an FDR ≤ 1%. Furthermore, since we expected to observe overlapping themes, we collapsed these terms into “supra-categories.”
Network analysis
To further analyze the biological impact of the DEGs, we created protein-protein interaction networks using the STRING database [24]. In order to assess only those interactions more likely to take place, we used a confidence score cutoff of 900, required experimental evidence in order to consider protein-protein interaction, and used only first-order interactions (meaning only molecules directly interacting with our DEG genes).
We also used the InnateDB [25], another curated database of experimentally verified proteins interactions and signaling pathways involved in the innate immune response.
GBP5 analysis
Co-culture experiments
Bone marrow-derived MSCs (BM-MSC) from 2 donors and peripheral blood lymphocytes (PBL) from 4 donors were co-cultured on transwell 12-well culture plates (#3460, Corning® Transwell®) in RMPI 1640 medium (#BE12-167F, Lonza®) supplemented with 10% FBS and antibiotics at 37 °C in a 5% CO2 atmosphere. In each well, 100,000 MSCs were seeded and then added 800,000 PBLs in the upper chamber. Previously, when needed, PBLs were activated with Dynabeads® Human T-Activator CD3/CD28 for T Cell Expansion and Activation (#11131D, Gibco®). Finally, INF-ϒ was neutralized using anti-INF-ϒ antibody (10 μl/ml, B27, BioLegends®)
Five conditions were analyzed: BM-MSCs with anti-IFN-ϒ, BM-MSCs with PBLs, BM-MSCs with PBLs and anti-IFN-ϒ antibody, BM-MSCs with activated PBLs, and BM-MSCs with activated PBLs and anti-IFN-ϒ antibody.
Quantitative PCR
Following 3 days of co-culture, medium and PBL were removed and RNA from MSCs were extracted using a commercial RNA extraction kit (SPEEDTOOLS Total RNA Extraction Kit, #21.212-4210, Biotools®). For cDNA generation, Superscript® VILO® cDNA Synthesis Kit (#11754050, Invitrogen®) was employed. Finally, cDNA was resuspended in 20 μl nuclease-free water and stored at − 20 °C until required for PCR.
The expression of three genes was quantificated with specific TaqMan® assays: 18S (Hs99999901_s1, #4331182), ACTB (Hs99999903_m1, #4453320), and GBP5 (Hs00369472_m1, #4448892). Quantitative PCR was carried out following master mix indications (TaqMan® Fast Advanced Master Mix, #4444557, Applied Biosystems®) on a MasterCycler RealPlex4 PCR System (Eppendorf). A triplicate of each sample was done. Fold change of GBP5 for each condition was calculated via 2−ΔΔCt method.
