Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Evaluating the Stability of RNA-Seq Transcriptome Profiles and Drug-Induced Immune-Related Expression Changes in Whole Blood

  • John F. Bowyer ,

    john.bowyer@fda.hhs.gov

    Affiliation Division of Neurotoxicology, National Center for Toxicological Research, United States Food and Drug Administration, Jefferson, Arkansas, United States of America

  • Karen M. Tranter,

    Affiliation Division of Neurotoxicology, National Center for Toxicological Research, United States Food and Drug Administration, Jefferson, Arkansas, United States of America

  • Joseph P. Hanig,

    Affiliation Center for Drug Evaluation and Research, United States Food and Drug Administration, Silver Spring, Maryland, United States of America

  • Nathaniel M. Crabtree,

    Affiliation Division of Neurotoxicology, National Center for Toxicological Research, United States Food and Drug Administration, Jefferson, Arkansas, United States of America

  • Robert P. Schleimer,

    Affiliation Division of Allergy and Immunology, Northwestern Feinberg School of Medicine, Chicago, Illinois, United States of America

  • Nysia I. George

    Affiliation Division of Bioinformatics and Biostatistics, National Center for Toxicological Research, United States Food and Drug Administration, Jefferson, Arkansas, United States of America

Abstract

Methods were developed to evaluate the stability of rat whole blood expression obtained from RNA sequencing (RNA-seq) and assess changes in whole blood transcriptome profiles in experiments replicated over time. Expression was measured in globin-depleted RNA extracted from the whole blood of Sprague-Dawley rats, given either saline (control) or neurotoxic doses of amphetamine (AMPH). The experiment was repeated four times (paired control and AMPH groups) over a 2-year span. The transcriptome of the control and AMPH-treated groups was evaluated on: 1) transcript levels for ribosomal protein subunits; 2) relative expression of immune-related genes; 3) stability of the control transcriptome over 2 years; and 4) stability of the effects of AMPH on immune-related genes over 2 years. All, except one, of the 70 genes that encode the 80s ribosome had levels that ranked in the top 5% of all mean expression levels. Deviations in sequencing performance led to significant changes in the ribosomal transcripts. The overall expression profile of immune-related genes and genes specific to monocytes, T-cells or B-cells were well represented and consistent within treatment groups. There were no differences between the levels of ribosomal transcripts in time-matched control and AMPH groups but significant differences in the expression of immune-related genes between control and AMPH groups. AMPH significantly increased expression of some genes related to monocytes but down-regulated those specific to T-cells. These changes were partially due to changes in the two types of leukocytes present in blood, which indicate an activation of the innate immune system by AMPH. Thus, the stability of RNA-seq whole blood transcriptome can be verified by assessing ribosomal protein subunits and immune-related gene expression. Such stability enables the pooling of samples from replicate experiments to carry out differential expression analysis with acceptable power.

Introduction

From a clinical standpoint, transcriptomic profiling of human blood (peripheral), or its various cells types, can be a useful tool for identifying biomarkers related to disease, immune response, and exposure to irritants [15]. Transcriptome analysis of blood in animal studies has also shown potential for identifying candidate blood-based biomarkers [6]. Until recently, much of the prior research has been conducted using microarray technologies. However, due to greater precision, increased coverage, and larger dynamic range, RNA sequencing (RNA-seq) platforms offer a competitive alternative to analyze changes in transcriptome expression [7]. Although whole blood analysis faces a unique set of challenges because globin mRNA depletion is necessary to ensure data quality, the purported advantages of RNA-seq technologies in terms of dynamic range and novel differential expression have also been observed in whole blood transcriptome studies [2,8]. These researchers found that despite some of the advantages of RNA-seq over microarrays, RNA sequencing reported greater within-group variability in whole blood. Typically, RNA-seq expression demonstrates higher reproducibility and lower variability in genes with high expression [7,9,10] but globin reduction can affect signal variability in whole blood. Thus, we evaluated the stability of globin-depleted mRNA from whole blood RNA-seq expression over time. Since RNA-seq technologies have been approved for clinical use by the FDA [11], changes in blood transcriptome due to toxicity in laboratory animals will play a large role in translating such changes to human data.

Globin transcripts compose over 70% of the expressed mRNA signal isolated from whole blood, and their removal has been reported to increase the integrity of transcript expression [12,13]. However, it has also been suggested that globin transcript depletion affects RNA yields, which would likely affect expression profiles [13]. As an additional means of assessing the overall stability of globin-depleted whole blood RNA-seq profiles, we have studied two sets of important comparisons. One set of comparisons is based on subunits of the 80s ribosome. In humans and mammals, the 80s ribosome consists of over 70 protein subunits, each of which is encoded by its own specific transcript [14]. These subunits are relatively unaffected by most drug treatments and experimental manipulations and are often used as “housekeeping” genes for RT-PCR [1518]. Because all the subunits are necessary for 80s ribosomal function, their transcripts should be expressed at high levels (relative to the entire transcriptome).

The second comparison profiles transcripts found in the various types of leukocytes (white blood cells). Here, we compared the expression of transcripts specific for monocytes/macrophages, lymphocytes (T-cells and B-cells), and all leukocytes. In rat blood, there are roughly 10-times as many B- and T-cells (approximately equal in number) as monocytes. B-cells and T-cells make up approximately 80% of all leukocytes in rat blood and monocytes comprise only about 3%. Thus, one might suppose that the expression of monocyte specific transcripts to those specific for T-cells or B-cells might be 10-fold less and 30-fold less than transcripts expressed in all leukocytes. However, such an assumption is not correct for the chemokine, C-C motif, receptor 1 (Ccr1) since this gene is expressed at very high levels in monocytes. In contrast to this, neutrophils consist of over 10 to 20% of the white blood cells (WBC) in the circulating blood of rat but have relatively low levels of mRNA present, as is the case with most end stage cells including eosinophils and basophils [19].

The second primary objective of this study was to assess the stability of control and drug-treated transcriptome expression in the whole blood of Sprague-Dawley rats over time. Many laboratories often replicate the same experiment over time to ensure reproducibility of significant changes in expression resulting from drug or toxicant exposure. Generally, data from all time points are then pooled. However, when detecting modest changes (e.g. 1.5- to 3-fold), pooling data can be problematic if the measured endpoints fluctuate by 2-fold or more over time in the control group. In the current study, stability was assessed across four time points (ranging over a 2 year period) for a control group and a treatment group. Stability was evaluated primarily through comparisons of the expression profile of ribosomal subunit proteins and secondarily through comparisons profiling the entire transcriptome.

Methods

Animal Dosing and Sacrifice

This study was carried out in accordance with the declaration of Helsinki and the Guide for the Care and Use of Laboratory Animals as adopted and promulgated by the National Institutes of Health. The use of animal testing in this study was done under protocols E7295 and E7519 (issued to John Bowyer) that were approved by the NCTR institutional animal care and use committee (IACUC) which is fully accredited (Food and Drug Administration—National Center for Toxicological Research Accreditation #A4310-01) by NIH-OLAW. Nine- to ten-week-old (≈ 65 days) male Sprague-Dawley rats, ninety five total, were obtained from the Charles River Laboratories [Crl:CD(SD)]. Upon arrival at NCTR they received tail tattoos for identification. Prior to testing, rats were housed 2 per cage with food and water available ad libitum. Rats were housed on a daily 12 hr light cycle with lights on at 6:00 am and off at 6:00 pm. During housing, the temperature (23°C) and humidity (53%) were controlled. The rats were tested between 12 and 13 weeks of age (≈ 87 days old).

The data used in this study (a cohort of a larger study contained in GSE62368 and GSE64778) is comprised of data from n = 22 saline dosed “control” rats and n = 21 rats given a neurotoxic exposure to AMPH. The larger study was carried out to identify mRNA biomarkers and to determine transcriptional immune-related responses present in circulating blood after neurotoxic exposures to hyperthermia, neurotoxic AMPH exposure, and non-neurotoxic AMPH exposure. Herein, we focus solely on the control and neurotoxic AMPH treatment groups. Multiple replicates of the same experiment were conducted over a 2 year period. This was done to ensure that any significant effects seen in transcript expression would be robust and reproducible over time and hopefully subsequently observable in other laboratories. The sample size of the control group at each time period was: Con1: n = 3, Con2: n = 7, Con3: n = 6, and Con4: n = 6. The size of the AMPH groups were: AMPH1: n = 5, AMPH2: n = 3, AMPH3: n = 6, and AMPH4: n = 7. The saline control and AMPH groups were sacrificed on: Con1 and AMPH1, September 2012; Con2 and AMPH2, June 2013; Con3 and AMPH3, Aug 2013; Con4 and AMPH4, July 2014.

Dosing commenced at 8:00 a.m. and ended at 2:00 p.m. AMPH-exposed animals were given 4 doses of amphetamine comprised of a sequential exposure to 5, 7.5, 10, and 10 mg/kg AMPH subcutaneously with 2 hr between each dose. The d-amphetamine-sulfate (Sigma-Aldrich, St. Louis, MO) dose was dissolved in normal saline (1 ml/ kg injected). This AMPH dosing paradigm was modified slightly from an original acute methamphetamine dosing paradigm that was shown to produce neurotoxicity [20,21]. The saline animals received 4 injections of 1 ml/kg s.c. In all groups, the behavior and body temperature of each rat were monitored hourly during saline or AMPH exposure and until at least 3 hr after the last dose (time of sacrifice). The lethal effects of hyperthermia and hyperpyrexia in the AMPH groups when body temperatures exceeded 41.4°C were prevented by placing the animals unrestrained on crushed ice for 15 to 30 min in a clean, wood chip free cage to allow their temperatures to drop below 40.0°C.

Collection and processing of cardiac blood

All rats were sacrificed at 4 to 5 pm with an overdose of 300 mg/kg body weight of pentobarbital, resulting in deep anesthesia in less than 3 min. At that point, 3 to 5 ml of cardiac blood was withdrawn from the heart using an 18 gauge needle attached to a 5 ml syringe, and the rats were then euthanatized with decapitation. Approximately 1 ml of the blood was immediately separated into approximately three 300 μl aliquots, frozen on dry ice in cylindrical foil capsules, and then stored at -70°C for later RNA-seq analysis. The remainder of the blood was injected into a 7.0 ml BD (Franklin Lanes, NJ) Vacutainer containing 12 mg EDTA (data from these aliquots are to be reported at later date). One ml was used to count the total RBC and WBC numbers as well as the different types of WBC present. Within 12 h of collection of the cardiac blood in EDTA coated vacutainers, one ml was used to determine total RBC, WBC counts, and the numbers of lymphocytes (T- and B-cells), monocytes, neutrophils, basophils and eosinophils present. Complete blood counts were determined on an ABX Pentra 60 C+ analyzer (ABX, Irvine CA). Maintenance and calibration was done according to the manufacturer’s recommendations. Three levels of assayed controls were included in daily analyses as internal controls.

RNA isolation and RNA-seq processing

RNA isolation was performed using RNAzol BD (Molecular Research Center, Inc., www.mrcgene.com) [9] with modified procedures for frozen blood. Just prior to RNA isolation, the frozen 300–400 μl RBC or whole blood aliquots were placed on dry ice and cut into approximately 100 to 150 μl sections. Two or three of these smaller sections were immediately (were not allowed to thaw) homogenized in 1.6 ml of RNAzol BD containing 45μl of 5 N acetic acid using a 2 ml Teflon pestle homogenizer. Approximately 300 to 400 μl of frozen whole blood was used for RNA isolation. Separation of the upper phase containing the extracted RNA with lower phases containing protein and DNA was accomplished by adding 350 μl of 4-bromoanisole and centrifugations at 12,000 x g for 10 min. To precipitate the RNA, 600 μl 2-propanol was added to 600 to 700 μl the extracted upper phase. After 15 min, centrifugation at 12,000 x g was performed. The RNA ppt was washed once with 70% ethanol and water and a second wash of 100% ethanol. After ≈ 1 min air drying the RNA was dissolved in 25 to 35 μl of RNase-free water. The final total cellular RNA recovered ranged from 5 to 25 μg and was stored at -70°C until shipment for data processing. Two to four μl of each aliquot was set aside for purity analysis at NCTR using an Agilent 2200 TapeStation System (Agilent Technologies, Palo Alto, CA).

The isolated total RNA was then shipped overnight on dry ice to Expression Analysis Inc. [EA; Durham, NC] for globin transcript removal and sequencing. Alpha and Beta Globin mRNA were substantially depleted from total RNA samples using the GlobinClear-Mouse/Rat Kit (Life Technologies # AM1981), essentially as described by the vendor. Briefly, 1.25 μg of total RNA isolated from whole blood was combined with biotinylated capture oligonucleotides complementary to globin mRNAs and the mixture incubated at 50°C for 15 minutes to allow duplex formation. Streptavidin magnetic beads were added to each specimen, and the resulting mixture was incubated for an additional 30 minutes at 50°C to allow binding of the biotin moieties by Streptavidin. These complexes, comprising Streptavidin magnetic beads bound to biotinylated oligonucleotides that are specifically hybridized to globin mRNAs, were then captured using a magnet. The globin-depleted supernatant is transferred to a new container and further purified using RNA binding beads. The final globin mRNA-depleted RNA samples are quantitated by spectrophotometry using a NanoDrop ND-8000 spectrophotometer.

RNA-seq Expression Profiling

EA created the library and performed 50bp paired-end and strand-specific sequencing using an Illumina platform. A total of 25–35 million reads were generated per sample. Sequencing by synthesis methods, as implemented via Illumina technology, were used to generate the RNA-seq data. Globin mRNA-depleted RNA samples were converted into cDNA libraries using the TruSeq Stranded mRNA Sample Prep Kit (Illumina, #RS-122-2103). Starting with 300 ng of globin mRNA-depleted RNA, polyadenylated RNA (primarily mRNA) was selected and purified using oligo-dT conjugated magnetic beads. This mRNA was chemically fragmented and converted into single-stranded cDNA using reverse transcriptase and random hexamer primers, with the addition of Actinomycin D to suppress DNA-dependent synthesis of the second strand. Double-stranded cDNA was created by removing the RNA template and synthesizing the second strand in the presence of dUTP in place of dTTP. A single A base was added to the 3’ end to facilitate ligation of sequencing adapters, which contain a single T base overhang. Adapter-ligated cDNA was amplified by polymerase chain reaction to increase the amount of sequence-ready library. During this amplification the polymerase stalls when it encounters a U base, rendering the second strand a poor template. Accordingly, amplified material used the first strand as a template, thereby preserving the strand information. Final cDNA libraries were analyzed for size distribution using an Agilent 2200 TapeStation (D1000 Screentape, Agilent # 5067–5582), quantitated by qPCR (KAPA Library Quant Kit, KAPA Biosystems # KK4824), then normalized to 2 nM in preparation for sequencing.

The 2nM normalized samples were pooled. cDNA templates were denatured using fresh 0.1N NaOH, and then diluted to a final loading concentration of 13pM. Cluster generation was performed on an Illumina cBOT (v1.5.12.0) using an Illumina TruSeq Paired-End Cluster Kit v3 (Illumina # PE-401-3001) to cluster an Illumina Paired End Flow Cell. Templates were attached to the flowcell via a dense lawn of oligonucleotides that bind to the sequencing adapters added during sample preparation, which are extended and then denatured. The flowcell was then sequenced through 51 bases, paired end, with an 8 base index cycle on an Illumina HiSeq 2000 (HiSeq Control Software v 1.5.15.1). During sequencing cycles, fluorescent reversible terminator dNTPs were added to the clusters, with only a single base per target being incorporated. Following imaging of the clusters, the terminator and fluorescent tag were cleaved so that the next base could be incorporated. Quality control files contained read length and depth results (before and after clipping), presence of artifact/ duplicate sequences, distribution of base quality and base frequency by sample. Also, flow cell total yield, PF reads, barcode quality, alignment summaries and number of genes detected were determined. The basecall files were converted to fastq files using CASAVA 1.8.2. The fastq files were clipped using fastq-mcf with the parameters “—max-ns 4—qual-mean 25-H-p 5-q 7-l 25” [22]. The fastq files were aligned to the rat Ensembl release 70 transcriptome (rn5) using bowtie 0.12.9 with the parameters “-e 500 –m 100 –chunkmbs 256”. The alignments were quantified using RSEM v 1.2.0 with no special parameters.

RNA-seq Data Analysis

RNA-seq by Expectation-Maximization (RSEM) [23] data generated by EA was rounded to produce digital expression values. In total, 16,881 transcripts were mapped to the reference rat genome in all samples. Transcript-level outliers within each experimental group were identified by the iLOO method [24] and subsequently replaced with the trimmed mean. All reports of average transcript expression are summarized using DESeq2 normalized counts [25]. DESeq2 was also used to carry out differential expression analysis. Significance was assessed using p-values adjusted by the Benjamini-Hochberg method for multiple comparison testing. Transcripts with adjusted p-values lower than 0.05 and absolute log2(FC) values exceeding specific cutoffs were declared significant. All reported values of correlation reflect Pearson’s correlation coefficient and are significant at α = 0.05.

A number of criteria were used to identify the most stable transcripts for RNA-seq data. Ideally, a stable transcript would show minimal variability, be consistently expressed in the high region of expression, and demonstrate consistent expression across treatment groups. Thus, candidates for stable expression were transcripts with no detected outliers, no measurements in the low expressed region [26], and no significant difference between the control and treatment group.

All analyses were carried out in R (http://www.r-project.org).

Results

The RNA integrity number (RIN) for the control whole blood samples (n = 22) ranged from 7.0 to 9.4, averaging 8.3 ± 0.15 (mean ± SEM). The RIN for the AMPH group ranged from 6.9 to 9.3, averaging 8.1 ± 0.15 (mean ± SEM), and did not significantly differ from control. Although globin depletion is reported to reduce expression [2], hemoglobin alpha (Hba) remained the highest expressed transcript. Also, as would be expected, the expression of olfactory receptor genes (approximately 1,200 unique transcripts [10]) was virtually absent in blood with only 14 genes having counts greater than 1. In addition, we assessed whether high levels of myoglobin transcript (Mb) were present in the blood samples. Since cardiac muscle contains very high levels of myoglobin [27], high expression of the myoglobin transcript could indicate that cardiac muscle was “contaminating” the samples as a result of the cardiac punctures to obtain whole blood. Two of the control animals from the Con4 group had very high levels (>5000 counts) for Mb. They were removed from further analysis resulting in n = 20 for the control group. None of the animals in the AMPH group had to be removed from the study due to high Mb levels. Mb expression for the remaining control and AMPH animals was relatively low with mean expressions of 29.8 ± 3.4 and 35.8 ± 1.8, respectively.

Average expression for the 70+ genes necessary to encode the 80s ribosome are shown in Table 1. The large 70s subunit transcripts are labeled with an NCBI gene symbol beginning with “Rpl,” while the small 30s subunit transcripts begin with “Rps.” Average expression for most of the transcripts was ranked in the top 5th percentile of expression levels (i.e. rank of 844 or higher). All of the ribosomal subunit genes, except Rpl22, had at least one transcript (splice variant) in the top 5th percentile; Rpl22 ranked 1010 in whole blood. RNA-seq expression of the ribosomal proteins for control animals were highly correlated (r≥0.98 for log2 expression data) across the first 3 time points (Fig 1).

thumbnail
Fig 1. Comparison of the Ribosomal Subunit Expression among Control Groups.

Pearson correlation and scatterplot matrix of log2 normalized expression of the ribosomal subunit proteins in control groups, Con1-Con4B. The red solid line is the identity line. Note that the aberrant/anomalous results of Con4A are due to less than optimal amplification prior to sequencing.

https://doi.org/10.1371/journal.pone.0133315.g001

thumbnail
Table 1. Expression of Rat Ribosomal Subunit Transcripts in the Whole Blood of Control and Amphetamine Treated Rats.

https://doi.org/10.1371/journal.pone.0133315.t001

Initial sequencing of the Con4 group (denoted by Con4A) yielded poorer correlation with the other groups (r≥0.74 for log2 expression data) (Fig 1). To some extent, diminished correlation between Con4A and Con1 – Con3 was also evident when assessing profiles of the entire transcriptome (Fig 2). Interestingly, this was the case despite the fact that total clusters, genes detected, % genes detected, % transcriptome mapped, skewness, median coverage, and % of all genes mapped did not differ between Con1, Con2, Con3, and Con4A. One noticeable fluctuation was the slight decrease in %(G+C), which ranged from 43% to 48% for Con4A, compared to the usual 48% to 53% range for the other groups. The less than optimal sequencing results shown with Con4A were due to an Illumina library preparation kit that contained a less than optimal PCR Master Mix reagent.

thumbnail
Fig 2. Comparison of All Transcripts among the Control Groups.

Pearson correlation and scatterplot matrix of log2 normalized expression of all transcripts in control groups, Con1-Con4B. The red solid line is the identity line. Again, the aberrant results of Con4A are due to less than optimal amplification prior to sequencing.

https://doi.org/10.1371/journal.pone.0133315.g002

Subsequently, aliquots of the same RNA from the 4 animals in the Con4A group were re-sequenced (denoted by Con4B), yielding a more accurate expression of the ribosomal proteins (Fig 1) and comparable %(G+C) values. Differential expression (DE) analysis of pairwise comparisons of the four control groups indicated that error rates (i.e. DE transcripts) were well below the nominal level. For example, DE analysis comparing Con1 and Con4B, the most distant pairwise time points, returned only 3 significant transcripts. Thus, although the groups were dosed and sacrificed two years apart, stability of transcript expression in control animals was maintained. When profiles of the ribosomal subunit expression for the four AMPH groups (AMPH1, AMPH2, AMPH3 and AMPH4) were compared to each other, pairwise correlation values were comparable to those observed for the control groups (correlation values among the four pairs ranged between 0.97 and 0.98). Somewhat surprisingly, correlation between each control group and the time-matched AMPH group was at least 0.97 (Fig 3).

thumbnail
Fig 3. Expression of Ribosomal Subunits in the Control versus Amphetamine Groups.

Pearson correlation and scatterplot matrix of log2 normalized expression of the ribosomal subunit proteins comparing each control group to its time-matched AMPH group. The red solid line is the identity line.

https://doi.org/10.1371/journal.pone.0133315.g003

To further characterize repeatability of the control expression profiles, 29 genes specifically related to the immune system and found predominately on leukocytes were evaluated. The mean (and standard error) of each gene is reported in Tables 2 and 3 for both groups (all four subgroups were included in the calculations). Most of the selected genes are classified as cluster of differentiation antigens (i.e. classification determinant) and are located on the outer membranes of leukocytes. The seven transcripts primarily expressed on monocytes/macrophages (Ccr1, Cd14, Cd163, Cd68, Itgax, Nos2 and Arg1) had a 1000-fold range if indeed interleukin 1 β is exclusively produced in monocytes (Il1b counts reached 18,000). Ten genes preferentially expressed in either T-cells (Cd2, Cd3e, Cd8a, Cd28, Cd247) or B-cells (Cd19, Cd22, Cd79b, Cd180 and Ly86) had a median expression of 1921 with only slightly more than a 10-fold range. The median expression of 6 selected transcripts found in almost all leukocytes (Cd300a, Cd44, Cd48, Cd84, Cd97 and Sell) was 6342 but ranged over 70-fold. The neutrophil related genes (Camp, Elane, Mpo, Mme, Prtn3 and Ctsg) had expression levels below 100 except for Camp.

thumbnail
Table 2. Expression of Selected Immune-Related Transcripts in Control and Amphetamine Treated Rats.

https://doi.org/10.1371/journal.pone.0133315.t002

thumbnail
Table 3. Citations for Gene Expression within Leukocyte Cell Types.

https://doi.org/10.1371/journal.pone.0133315.t003

The expression of these immune-related genes had near-1 correlation (r≥0.97) among the four control groups (Con1-Con3 and Con4B) (Fig 4). However, as expected, the correlation between each control group and its paired time-matched AMPH group was not as strong because of the known effect of the drug on the immune system (Fig 5). Fig 6 presents the log2 fold change of AMPH to control expression for the selected immune-related genes, except for the three neutrophil-related genes with very low expression. Relative expression is presented for all four time-matched groups and also for all of the data (denoted by a purple diamond). AMPH has a significant effect on 15 of the 26 genes. Relative to control, all of the genes primarily expressed in T-cells decreased approximately 1.5- to 3.0-fold, while 4 of the 7 monocyte-specific transcripts increased between 2.0 and 6-fold. Expression of 3 of the 5 genes relatively specific for B-cells decreased around 1.5-fold, and genes that are normally expressed in most leukocytes were not significantly different (less 1.5-fold), with only two of the six genes (Cd300a and Cd84) reaching significance between AMPH and control.

thumbnail
Fig 4. Comparison of Immune-Related Transcripts among the Control Groups.

Pearson correlation and scatterplot matrix of log2 normalized expression of the selected immune-related genes in control groups, Con1-Con4B. The red solid line is the identity line.

https://doi.org/10.1371/journal.pone.0133315.g004

thumbnail
Fig 5. Expression of Immune-Related Genes in the Control versus Amphetamine Groups.

Pearson correlation and scatterplot matrix of log2 normalized expression of the selected immune-related genes comparing each control group to its time-matched AMPH group. The red solid line is the identity line.

https://doi.org/10.1371/journal.pone.0133315.g005

thumbnail
Fig 6. Expression Changes in Genes Primarily Expressed in Monocytes, B-cells or T-cells, and all Leukocytes as a Result of Neurotoxic AMPH Exposure.

Scatterplot of log2 fold change (AMPH to control) expression of the selected immune-related genes. For each gene, values of log2FC are presented for each time-matched subgroup comparison using a different colored dot for each subgroup (AMPH1/Con1, black; AMPH2/Con2, red; AMPH3/Con3, blue and AMPH4/Con4B green). The pooled expression of AMPH to control over all subgroups (AMPH/Con) is represented by a purple diamond. Differential expression was assessed for the latter comparison. Log2FC expression for genes with adjusted p-values smaller than 0.05 and 0.58 < |log2FC| < 1 is marked by one asterisk. Likewise, expression for genes with adjusted p-values smaller than 0.05 and |log2FC|≥ 1 is marked with two asterisks.

https://doi.org/10.1371/journal.pone.0133315.g006

The levels of the major types of leukocytes in blood for control and AMPH are presented in Table 4. The levels of the lymphocytes (72.3%), neutrophils (23.6%), monocytes (3.5%), eosinophils (0.4%) and basophils (0.2%) were as would be expected based on standard reported levels Monocytes were significantly upregulated by AMPH; the levels of neutrophils were also upregulated by AMPH and marginally significant. Lymphocytes were significantly down-regulated by AMPH.

thumbnail
Table 4. Analysis of Leukocyte Cell Types Present in the Whole Blood of Control and Amphetamine Treated Rats.

https://doi.org/10.1371/journal.pone.0133315.t004

Discussion

The results of our study demonstrate the stability of whole blood RNA-seq transcriptome profiles after globin transcript depletion. Globin depletion was implemented since it has been observed to improve the transcriptome profile in whole blood in previous studies using RNA-seq methods [2]. Any loss of RNA due to globin removal was reasonably consistent since comparisons of transcriptome profiles over a period of two years were highly correlated for both the control and treatment group. The level of agreement is also attributed to consistent expression of each experimental group across time. Our study provides evidence that the consistency and repeatability of RNA-seq of whole blood can be validated by profiling expression of the ribosomal subunits. However, further studies will be necessary to determine whether this approach is more sensitive than methods such as the detection of aberrant G-C profiles [28]. The utility of RNA-seq data from whole blood was further demonstrated by the completeness of the expression profile of immune-related genes as well as the consistent effect that drug treatment had on the expression of these genes, at least in the case of a toxic exposure to AMPH [29].

The 80s ribosome is a group of more than 70 transcripts, of which are all expressed at relatively very high levels [14]. Note that there are 10 transcripts with lower expression; however, the gene for each of these transcripts produced another transcript (splice variant) that measured at much higher levels. As seen in this study, the one anomalous group of RNA-seq data (Con4A) that was less than optimal showed large variability in these transcript levels (Fig 1). This was solely due to a very rare instance of amplification kit variability. Although the observed variability in Con4A expression was somewhat apparent in the entire transcriptome (Fig 2), the deviant variation was readily apparent when profiling the ribosomal subunits. As an example, Pearson correlation values for Con4A were at least 0.90 for the whole genome but were as low as 0.74 when comparing expression levels of the ribosomal transcripts.

An additional benefit of using ribosomal transcripts rather than the entire genome is that the majority of ribosomal transcripts do not change much with AMPH treatment. Only two of the ribosomal transcripts were significantly differentially expressed between control and AMPH (p<0.05 and FC>2). However, it is not yet known how other types of drug/toxicant treatments will affect subunit expression relative to control. Several ribosomal proteins have been used as “housekeeping” genes to normalize RT-PCR expression across individual RNA samples, although they may not be optimal in all tissues [1114]. As an example, in our study, ribosomal protein L7, Rpl7, was identified as the highest consistently expressed transcript across all 41 samples. Finally, a relatively high expression for at least one transcript of each gene in the ribosomal subunit is almost as important as consistent expression. This is because expression of all of the genes is necessary to form one of the most numerous protein complexes in an organism. Our data showed that at least one transcript for each ribosomal gene ranked in the upper 5th percentile of expression. One could argue that the use of a set of 10 or so genes commonly used to normalize RT-PCR data would provide a better way to determine the potential for significant variation in whole blood RNA-seq. However, such analysis would be limited in some respects because the specific genes in a set of housekeeping genes for RT-PCR are not linked by a common process or protein complex. Thus, unlike the set of ribosomal transcripts, it is much more difficult to determine whether the relative levels of the various genes previously used to normalize RT-PCR are as they should be (in relationship to each other).

In this study, we profiled the expression of immune-related genes (mostly cell determinant genes that are expressed on leukocytes). The stability and consistency of the expression of these immune-related genes are also important as these are the genes in whole blood, along with cytokines and chemokines protein levels in blood, that are used to evaluate how the immune system is responding to a given disease, treatment or toxic insult [3037]. These genes are major effectors in the chemokine and cytokine signaling system of leukocytes. There have been previous reports on how gene expression, as determined by microarray or RNA-seq, is affected in human blood in anorexic patients or after acute ethanol exposure [38,39]. However, these studies focused primarily on genes not necessary directly related to cytokines, chemokines, or receptors that are exclusively or predominately localized to leukocytes. One very recent RNA-seq study investigated the effects of immunosuppression on the expression of genes related to cytokine/chemokine signaling in mononuclear cells isolated from peripheral blood in patients involved in [40]. A pronounced 10-fold or more down regulation of these genes was observed due to the immunosuppression. Our data clearly indicates that RNA-seq methods are also a sensitive method for detecting changes in gene expression related to cytokine and chemokine signaling resulting from toxic exposures, in this case AMPH.

We observed that the levels of transcripts relatively specific for monocytes, T-cells, B-cells, neutrophils or present in all leukocytes were roughly as expected in the circulating blood of rat. Transcript levels relatively specific for monocytes, the least abundant cell type (3%), were readily detectable in the control and amphetamine groups, and some such as Ccr1 had very high expression. All but one of the genes specific for neutrophils, which are 8 to 15% of the all WBCs in rat, were at low levels. This would be expected due to the nature of mature neutrophils present in circulating blood. However, some genes, mostly those also present in other types of WBCs, can be induced [19]. Therefore, some of the expression that we observed in blood of interleukin 1b (Il1b) and neutrophil cytosolic factor 1 (Ncf1), which were at very high levels (data not shown), could have been due to expression in neutrophils as well as monocytes. Thus, the overall profile of immune expression was well represented.

There was a significant difference between the control and amphetamine groups for several of the selected immune-related genes. Specifically, 10 of the 26 immune-related genes had 2-fold or greater change in expression and 5 additional genes demonstrated at least 1.5-fold change in expression. Two genes, Cd14 and Arg1, increased expression more than 4-fold as a result of amphetamine exposure. Importantly, the response (change in transcript levels) of all these genes, except Cd163, was very consistent across all four treatment groups. In general, many of the changes observed in these immune related genes might be expected to be due to AMPH exposure, including the increase in β-adrenergic stimulation produced by AMPH releasing norepinephrine [4144]. Finding treatment-induced changes in immune-related transcripts is of the utmost importance since previous, and likely subsequent, animal studies involving gene expression in blood or serum have focused on changes related to immune system function.

Our data on cell counts of monocytes and lymphocytes (T-cells and B-cells) in circulating blood indicate that the transcript expression changes may be partially explained by increased or decreased numbers of specific types of leukocytes present in circulating blood after AMPH. The monocyte count increased slightly over two fold while the lymphocyte levels decreased approximately 30%. In regards to the monocytes, the increase in their numbers could partially explain the increase in expression of genes specifically found in them. Nonetheless, the changes in numbers/counts cannot explain all the gene expression changes seen. The expression increases in monocyte-specific genes range from a less than 1.5-fold increase (Ccr1 and Itgax) to almost a 6-fold increase for Cd14, which would strongly argue that there is likely a relative change in the relative expression of the genes present within the circulating monocytes. Likewise, the 30% decreases in the number of circulating T-cells and B-cells may not explain the much greater fold change decreases, particularly in T-cell-related genes. Overall, the increases in the circulating blood expression profile and monocyte and neutrophil cells numbers indicate an activation of the innate immune system by AMPH [4548].

The stability of transcript expression in four control groups of animals (Sprague-Dawley rats in this study) over a two-year timeframe was strong. However, more importantly, the changes in gene expression of the animals given a neurotoxic dosage of AMPH relative to the controls were very consistent across all four time points. Thus, the work presented here demonstrates that under these conditions control and treatment-specific animals can be pooled for analysis. The increase in sample size will subsequently increase power to detect differentially expressed transcripts. The ability to show reproducibility of treatment effect on the transcriptome over time lends to the validation of the expression changes that are reported.

Author Contributions

Conceived and designed the experiments: RPS JFB NIG NMC JPH KMT. Performed the experiments: JFB KMT NMC. Analyzed the data: JFB NIG RPS NMC JPH KMT. Contributed reagents/materials/analysis tools: JFB KMT NIG NMC. Wrote the paper: JFB NIG RPS JPH NMC.

References

  1. 1. Mohr S, Liew CC (2007) The peripheral-blood transcriptome: new insights into disease and risk assessment. Trends Mol Med 13: 422–432. pmid:17919976
  2. 2. Raghavachari N, Xu X, Munson PJ, Gladwin MT (2009) Characterization of whole blood gene expression profiles as a sequel to globin mRNA reduction in patients with sickle cell disease. PLoS One 4: e6484. pmid:19649296
  3. 3. Kam SH, Singh A, He JQ, Ruan J, Gauvreau GM, O'Byrne PM, et al. (2012) Peripheral blood gene expression changes during allergen inhalation challenge in atopic asthmatic individuals. J Asthma 49: 219–226. pmid:22316092
  4. 4. Brouard S, Mansfield E, Braud C, Li L, Giral M, Hsieh SC, et al. (2007) Identification of a peripheral blood transcriptional biomarker panel associated with operational renal allograft tolerance. Proc Natl Acad Sci U S A 104: 15448–15453. pmid:17873064
  5. 5. Broderick G, Ben-Hamo R, Vashishtha S, Efroni S, Nathanson L, Barnes Z, et al. (2013) Altered immune pathway activity under exercise challenge in Gulf War Illness: an exploratory analysis. Brain Behav Immun 28: 159–169. pmid:23201588
  6. 6. O'Loughlin A, Lynn DJ, McGee M, Doyle S, McCabe M, Earley B (2012) Transcriptomic analysis of the stress response to weaning at housing in bovine leukocytes using RNA-seq technology. BMC Genomics 13: 250. pmid:22708644
  7. 7. Marioni JC, Mason CE, Mane SM, Stephens M, Gilad Y (2008) RNA-seq: an assessment of technical reproducibility and comparison with gene expression arrays. Genome Res 18: 1509–1517. pmid:18550803
  8. 8. Raghavachari N, Barb J, Yang Y, Liu P, Woodhouse K, Levy D, et al. (2012) A systematic comparison and evaluation of high density exon arrays and RNA-seq technology used to unravel the peripheral blood transcriptome of sickle cell disease. BMC Med Genomics 5: 28. pmid:22747986
  9. 9. Li S, Tighe SW, Nicolet CM, Grove D, Levy S, Farmerie W, et al. (2014) Multi-platform assessment of transcriptome profiling using RNA-seq in the ABRF next-generation sequencing study. Nat Biotechnol 32: 915–925. pmid:25150835
  10. 10. Zhao S, Fung-Leung WP, Bittner A, Ngo K, Liu X (2014) Comparison of RNA-Seq and microarray in transcriptome profiling of activated T cells. PLoS One 9: e78644. pmid:24454679
  11. 11. Sheridan C (2014) Milestone approval lifts Illumina's NGS from research into clinic. Nat Biotechnol 32: 111–112. pmid:24509734
  12. 12. Shin H, Shannon CP, Fishbane N, Ruan J, Zhou M, Balshaw R, et al. (2014) Variation in RNA-Seq transcriptome profiles of peripheral whole blood from healthy individuals with and without globin depletion. PLoS One 9: e91041. pmid:24608128
  13. 13. Mastrokolias A, den Dunnen JT, van Ommen GB, t Hoen PA, van Roon-Mom WM (2012) Increased sensitivity of next generation sequencing-based expression profiling after globin reduction in human blood RNA. BMC Genomics 13: 28. pmid:22257641
  14. 14. Yusupova G, Yusupov M (2014) High-Resolution Structure of the Eukaryotic 80S Ribosome. Annu Rev Biochem 83: 467–486. pmid:24580643
  15. 15. Cai JH, Deng S, Kumpf SW, Lee PA, Zagouras P, Ryan A, et al. (2007) Validation of rat reference genes for improved quantitative gene expression analysis using low density arrays. Biotechniques 42: 503–512. pmid:17489238
  16. 16. Huggett J, Dheda K, Bustin S, Zumla A (2005) Real-time RT-PCR normalisation; strategies and considerations. Genes Immun 6: 279–284. pmid:15815687
  17. 17. Martinez-Beamonte R, Navarro MA, Larraga A, Strunk M, Barranquero C, Acin S, et al. (2011) Selection of reference genes for gene expression studies in rats. J Biotechnol 151: 325–334. pmid:21219943
  18. 18. Vandesompele J, De Preter K, Pattyn F, Poppe B, Van Roy N, De Paepe A, et al. (2002) Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol 3: RESEARCH0034. pmid:12184808
  19. 19. Subrahmanyam YV, Yamaga S, Prashar Y, Lee HH, Hoe NP, Kluger Y, et al. (2001) RNA expression patterns change dramatically in human neutrophils exposed to bacteria. Blood 97: 2457–2468. pmid:11290611
  20. 20. Bowyer JF, Thomas MT, Schmued LC, Ali SF (2008) Brain region-specific neurodegenerative profiles showing the relative importance of amphetamine dose, hyperthermia, seizures and the blood-brain barrier. Ann NY Acad Sci 1139: 127–139. pmid:18991857
  21. 21. Sonsalla PK, Nicklas WJ, Heikkila RE (1989) Role for excitatory amino acids in methamphetamine-induced nigrostriatal dopaminergic toxicity. Science 243: 398–400. pmid:2563176
  22. 22. Aronesty E (2013) Comparison of sequencing utility programs. Open Bioinform J 7: 1–8.
  23. 23. Li B, Dewey CN (2011) RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics 12: 323. pmid:21816040
  24. 24. George NI, Bowyer JF, Crabtree NM, Chang CW (2015) An Iterative Leave-One-Out Approach to Outlier Detection in RNA-Seq Data. PLoS One 10: e0125224. pmid:26039068
  25. 25. Love M, Huber H, Anders S (2014) Moderated estimation of fold change and dispersion for RNA-Seq data with DESeq2. bioRxiv. Preprint Published online February 19: 2014.
  26. 26. George NI, Chang CW (2014) DAFS: a data-adaptive flag method for RNA-sequencing data to differentiate genes with low and high expression. BMC Bioinformatics 15: 92. pmid:24685233
  27. 27. Hendgen-Cotta UB, Flogel U, Kelm M, Rassaf T (2010) Unmasking the Janus face of myoglobin in health and disease. J Exp Biol 213: 2734–2740. pmid:20675542
  28. 28. Lehning EJ, Balaban CD, Ross JF, Reid MA, LoPachin RM (2002) Acrylamide neuropathy. I. Spatiotemporal characteristics of nerve cell damage in rat cerebellum. Neurotoxicology 23: 397–414. pmid:12387366
  29. 29. Bowyer JF, Patterson TA, Saini UT, Hanig JP, Thomas M, Camacho L, et al. (2013) Comparison of the global gene expression of choroid plexus and meninges and associated vasculature under control conditions and after pronounced hyperthermia or amphetamine toxicity. BMC Genomics 14: 147. pmid:23497014
  30. 30. Mirantes C, Passegue E, Pietras EM (2014) Pro-inflammatory cytokines: Emerging players regulating HSC function in normal and diseased hematopoiesis. Exp Cell Res.
  31. 31. Rabinovitch M, Guignabert C, Humbert M, Nicolls MR (2014) Inflammation and immunity in the pathogenesis of pulmonary arterial hypertension. Circ Res 115: 165–175. pmid:24951765
  32. 32. Wu D, Zhou J, Bi H, Li L, Gao W, Huang M, et al. (2014) CCL11 as a potential diagnostic marker for asthma? J Asthma 51: 847–854. pmid:24796647
  33. 33. Menard-Katcher C, Furuta GT (2014) Non- and semi-invasive methods of monitoring eosinophilic esophagitis. Dig Dis 32: 102–106. pmid:24603390
  34. 34. Locke BA, Dasu T, Verbsky JW (2014) Laboratory diagnosis of primary immunodeficiencies. Clin Rev Allergy Immunol 46: 154–168. pmid:24569953
  35. 35. Signorelli SS, Fiore V, Malaponte G (2014) Inflammation and peripheral arterial disease: the value of circulating biomarkers (Review). Int J Mol Med 33: 777–783. pmid:24535646
  36. 36. Visconti L, Nelissen K, Deckx L, van den Akker M, Adriaensen W, Daniels L, et al. (2014) Prognostic value of circulating cytokines on overall survival and disease-free survival in cancer patients. Biomark Med 8: 297–306. pmid:24521026
  37. 37. Wiseman S, Marlborough F, Doubal F, Webb DJ, Wardlaw J (2014) Blood markers of coagulation, fibrinolysis, endothelial dysfunction and inflammation in lacunar stroke versus non-lacunar stroke and non-stroke: systematic review and meta-analysis. Cerebrovasc Dis 37: 64–75. pmid:24401164
  38. 38. Kim Y, Trace SE, Crowley JJ, Brownley KA, Hamer RM, Pisetsky DS, et al. (2013) Assessment of gene expression in peripheral blood using RNAseq before and after weight restoration in anorexia nervosa. Psychiatry Res 210: 287–293. pmid:23778302
  39. 39. Kupfer DM, White VL, Strayer DL, Crouch DJ, Burian D (2013) Microarray characterization of gene expression changes in blood during acute ethanol exposure. BMC Med Genomics 6: 26. pmid:23883607
  40. 40. Dorr C, Wu B, Guan W, Muthusamy A, Sanghavi K, Schladt DP, et al. (2015) Differentially expressed gene transcripts using RNA sequencing from the blood of immunosuppressed kidney allograft recipients. PLoS One 10: e0125045. pmid:25946140
  41. 41. Bellinger DL, Lorton D (2014) Autonomic regulation of cellular immune function. Auton Neurosci 182: 15–41. pmid:24685093
  42. 42. Felten DL, Felten SY, Bellinger DL, Carlson SL, Ackerman KD, Madden KS, et al. (1987) Noradrenergic sympathetic neural interactions with the immune system: structure and function. Immunol Rev 100: 225–260. pmid:3326822
  43. 43. Nance DM, Sanders VM (2007) Autonomic innervation and regulation of the immune system (1987–2007). Brain Behav Immun 21: 736–745. pmid:17467231
  44. 44. Mahajan SD, Hu Z, Reynolds JL, Aalinkeel R, Schwartz SA, Nair MP (2006) Methamphetamine modulates gene expression patterns in monocyte derived mature dendritic cells: implications for HIV-1 pathogenesis. Mol Diagn Ther 10: 257–269. pmid:16884330
  45. 45. Akira S, Misawa T, Satoh T, Saitoh T (2013) Macrophages control innate inflammation. Diabetes Obes Metab 15 Suppl 3: 10–18. pmid:24003916
  46. 46. Di Santo JP (2014) Staying innate: transcription factor maintenance of innate lymphoid cell identity. Immunol Rev 261: 169–176. pmid:25123284
  47. 47. Lauvau G, Chorro L, Spaulding E, Soudja SM (2014) Inflammatory monocyte effector mechanisms. Cell Immunol 291: 32–40. pmid:25205002
  48. 48. Amulic B, Cazalet C, Hayes GL, Metzler KD, Zychlinsky A (2012) Neutrophil function: from mechanisms to disease. Annu Rev Immunol 30: 459–489. pmid:22224774
  49. 49. Cook DN, Beck MA, Coffman TM, Kirby SL, Sheridan JF, Pragnell IB, et al. (1995) Requirement of MIP-1 alpha for an inflammatory response to viral infection. Science 269: 1583–1585. pmid:7667639
  50. 50. Geissmann F, Jung S, Littman DR (2003) Blood monocytes consist of two principal subsets with distinct migratory properties. Immunity 19: 71–82. pmid:12871640
  51. 51. Biswas SK, Mantovani A (2010) Macrophage plasticity and interaction with lymphocyte subsets: cancer as a paradigm. Nat Immunol 11: 889–896. pmid:20856220
  52. 52. Bierer BE, Burakoff SJ (1988) T cell adhesion molecules. FASEB J 2: 2584–2590. pmid:2838364
  53. 53. Painter MW, Davis S, Hardy RR, Mathis D, Benoist C, Immunological Genome Project Consortium (2011) Transcriptomes of the B and T lineages compared by multiplatform microarray profiling. J Immunol 186: 3047–3057. pmid:21307297
  54. 54. Kioussis D, Ellmeier W (2002) Chromatin and CD4, CD8A and CD8B gene expression during thymic differentiation. Nat Rev Immunol 2: 909–919. pmid:12461564
  55. 55. Moretta A, Ciccone E, Pantaleo G, Tambussi G, Bottino C, Melioli G, et al. (1989) Surface molecules involved in the activation and regulation of T or natural killer lymphocytes in humans. Immunol Rev 111: 145–175. pmid:2697680
  56. 56. Moretta A, Bottino C, Vitale M, Pende D, Cantoni C, Mingari MC, et al. (2001) Activating receptors and coreceptors involved in human natural killer cell-mediated cytolysis. Annu Rev Immunol 19: 197–223. pmid:11244035
  57. 57. Horcher M, Souabni A, Busslinger M (2001) Pax5/BSAP maintains the identity of B cells in late B lymphopoiesis. Immunity 14: 779–790. pmid:11420047
  58. 58. Mansson R, Welinder E, Ahsberg J, Lin YC, Benner C, Glass CK, et al. (2012) Positive intergenic feedback circuitry, involving EBF1 and FOXO1, orchestrates B-cell fate. Proc Natl Acad Sci U S A 109: 21028–21033. pmid:23213261
  59. 59. Miyake K, Yamashita Y, Ogata M, Sudo T, Kimoto M (1995) RP105, a novel B cell surface molecule implicated in B cell activation, is a member of the leucine-rich repeat protein family. J Immunol 154: 3333–3340. pmid:7897216
  60. 60. Sigal LH (2012) Basic science for the clinician 59: polymorphonuclear cells: mechanisms in human defense and in the pathogenesis of autoimmune disease. J Clin Rheumatol 18: 443–449. pmid:23211587
  61. 61. Borrego F (2013) The CD300 molecules: an emerging family of regulators of the immune system. Blood 121: 1951–1960. pmid:23293083
  62. 62. Merzaban JS, Burdick MM, Gadhoum SZ, Dagia NM, Chu JT, Fuhlbrigge RC, et al. (2011) Analysis of glycoprotein E-selectin ligands on human and mouse marrow cells enriched for hematopoietic stem/progenitor cells. Blood 118: 1774–1783. pmid:21659548
  63. 63. Boles KS, Stepp SE, Bennett M, Kumar V, Mathew PA (2001) 2B4 (CD244) and CS1: novel members of the CD2 subset of the immunoglobulin superfamily molecules expressed on natural killer cells and other leukocytes. Immunol Rev 181: 234–249. pmid:11513145
  64. 64. Krause SW, Rehli M, Heinz S, Ebner R, Andreesen R (2000) Characterization of MAX.3 antigen, a glycoprotein expressed on mature macrophages, dendritic cells and blood platelets: identity with CD84. Biochem J 346 Pt 3: 729–736. pmid:10698700
  65. 65. Safaee M, Clark AJ, Ivan ME, Oh MC, Bloch O, Sun MZ, et al. (2013) CD97 is a multifunctional leukocyte receptor with distinct roles in human cancers (Review). Int J Oncol 43: 1343–1350. pmid:23969601
  66. 66. Tedder TF, Steeber DA, Chen A, Engel P (1995) The selectins: vascular adhesion molecules. FASEB J 9: 866–873. pmid:7542213