Study design

The study was conducted in Lima, Peru, where both DTP and DTaP vaccines are available and there is no preferential recommendation for their use. The study was conducted in accordance with the study protocol, International Council on Harmonization Good Clinical Practices standards, the Declaration of Helsinki, and applicable regulatory requirements. Both the Vanderbilt University Medical Center (VUMC) Human Research Protection Program and the Ethics Review Board of the Instituto de Investigación Nutricional (IIN) approved the study protocol and amendments.

Participants

Following written informed consent by a parent or legal guardian, healthy male and female infants ≥2 months of age were enrolled. A pilot phase of the study (n = 12 participants) was conducted first to optimize workflows and establish sufficiency of blood volume for subsequent analysis. Pilot study participants were included only in analyses from Day 150 (pertussis-protein reactive B-cell frequency and pertussis serology), Day 480 (pertussis serology), and Day 510 (reactive-B cells and pertussis serology) and did not contribute to correlation analyses involving RNA-Seq, RP, or cytokine data. The pilot study was followed by the main study phase (n = 44 participants). Infants from whom blood could not be obtained on Days 2 or 8 were excluded and replaced. Participants included infants receiving care in either the Peruvian Ministry of Health system who received DTP by standard of care or infants receiving care at private pediatric clinics who received DTaP as routine vaccination. Infants were eligible to participate if they were ≥37 weeks’ gestation at birth, 50–89 days of age at the time of enrollment, and considered to be in good general health by medical history and physical examination. Participants were excluded if they were of low birth weight (<2.5 kg); anemic (hemoglobin <7 mg/dL); had evidence of a primary immunodeficiency or HIV; were receiving systemic corticosteroids at enrollment or had received corticosteroids for more than 14 days since birth; or had a history of receiving any blood products since birth. Children with acute illnesses remained eligible, though enrollment was delayed until all symptoms were resolved for >48 h. Participants who received a vaccine that differed from the expected sequence, or participants for whom data were not available due to loss to follow-up/parent withdrawal, were excluded from analyses after the relevant event. Maternal pertussis immunization was not available during this time for pregnant individuals in Peru.

Study procedures

The overall study scheme and consort diagram are illustrated in Fig. 1. Briefly, 12 evaluable participants were enrolled during the pilot phase of the study. Following successful venipuncture, 7 infants received a hexavalent vaccine containing DTaP (Hexaxim, Sanofi Pasteur, containing PT/FHA, diphtheria, tetanus, hepatitis B, HiB, and inactivated poliovirus) and 5 infants received a pentavalent vaccine containing DTP (Serum Institute of India, containing whole-cell pertussis, diphtheria, tetanus, hepatitis B, and Hib), Table 1. Infants returned to the clinic on Day 2, approximately 24 h after vaccination, for physical examination and repeat venipuncture. Blood samples from Days 1 and 2 in the pilot study were used for RNA-Seq, RP, and cytokine measurements. Infants subsequently received doses 2 and 3 of DTP or DTaP per standard of care at 4 and 6 months of age and returned on Day 150 (approximately 1 month following completion of the primary series) for measurement of pertussis antibodies and enumeration of pertussis-specific B-cells. Where possible, the same vaccine formulation was administered for each vaccine dose. At the time of the pertussis vaccine booster (Day 480), pre-immunization blood was obtained for serum pertussis antibodies and on Day 510, approximately 1 month following the booster dose, blood was obtained for measurement of serum pertussis antibodies and enumeration of pertussis-specific reactive B-cells³³.

In the main study phase, infants were randomly allocated at enrollment to either a Day 2 or Day 8 blood sample to ascertain early and late transcriptomic changes in response to vaccine. As with the pilot study, infants in whom paired blood samples pre- and post-vaccination could not be obtained were withdrawn from the study and replaced. The remainder of the study schedule was unchanged during the main study phase from the pilot study.

All participants were provided a thermometer and memory aid to measure solicited adverse events (reactogenicity) for the 7 days following the first dose of vaccine. Unsolicited adverse events, deemed by the investigator to be due to vaccine, were also recorded throughout the study.

Sample preparation, RNA-Seq, and RP assays

Detailed descriptions of laboratory procedures are described elsewhere³⁴ accompanying manuscript, by Leguia et al.). Briefly, whole blood was collected on Day 1 (pre-vaccination) and on either Day 2 or Day 8 and immediately processed to maintain fidelity of transcripts for sequencing. In the pilot study, plasma was separated first, and peripheral blood mononuclear cells (PBMC) were recovered using SepMate and Lymphoprep (STEMCELL Technologies) according to manufacturer’s instructions. Due to insufficient recovery of ribosome footprints for RP in the pilot phase, we modified the approach in the main study and added 100 ug/mL cycloheximide (CHX), which freezes translating ribosomes on their in vivo positions on mRNA, to blood collection tubes prior to obtaining blood samples. As a result, only sera were available for cytokine analyses in the main study. PBMCs were frozen after separation. At the time of analysis, cells were removed from liquid nitrogen, thawed, and resuspended in lysis buffer; 80% of the cell lysate of each sample was used to prepare ribosome profile libraries, as described elsewhere³⁴ accompanying manuscript, by Leguia et al., and 20% of the lysate was used for RNA-Seq, using standard methods³⁵.

Serology

The humoral immune response to vaccine was assessed by measuring serum IgG antibody responses by liquid bead-based assay to pertussis toxin (PT), pertactin (PRN), fimbriae 2/3 (FIM, List Biologicals, Campbell, CA), and filamentous hemagglutinin (FHA, Enzo Biochem, Framingdale, NY)³⁶. Antigens were individually and covalently conjugated to carboxylate modified Luminex beads, each having a distinct spectral signature. The reportable value of the assay is the mean of two tests and expressed as the serum concentration of antigen-specific IgG in IU/mL; the assay is calibrated to the WHO International Standard Pertussis Antiserum WHO 06/140.

Pre- and post-vaccination differences in antibody concentrations for each pertussis antigen and each time point were evaluated for each group using a two-sided Welch’s t test on the log scale adjusting for unequal variance.

Cytokine assay

A commercially available 21-plex inflammatory cytokine panel (Milliplex, Millipore) was used to measure pro-inflammatory cytokines from plasma (pilot study) and sera (main study). Briefly, 25 µl of neat sera was added to 25 µl of assay buffer in a prepared flat bottom 96-well plate with standard and controls. All wells had 25 µl of the bead mixture added and incubated overnight at 4 °C. After washing, 50 µl of detection reagent was added for 1 h followed by 50 µl of streptavidin-phycoerythrin for 30 min. The plate was then acquired on a Luminex FlexMap 3D. Inflammatory cytokine concentrations, as well as their post-vaccination fold changes, were summarized using the median, 95% CI of the median, first and third quartile, minimum and maximum. The 95% CI of the median concentration and median fold change from pre-vaccination were determined using the bootstrap method (1000 replications) and were visualized using time trend plots. To identify inflammatory cytokines that showed a differential response from pre-vaccination (Day 1), a Wilcoxon signed-rank test was carried out for each cytokine, vaccine group, and post-vaccination time point (Days 2 and 8).

Reactive B-cell frequency assay

Methods for the enumeration of antigen-specific B cells have been published previously³⁷ and includes EBV-transformation of cells, formation of lymphoblastoid cell lines (LCL), screening for PT antibody, and calculation of the number of LCL’s in each well. Briefly, previously frozen PBMC were thawed rapidly in a 37 °C water bath and washed with a warm medium. Each well was seeded with 50 μL of PBMC resuspended in a medium with B95.8 cell supernatant at a concentration of 0.5 ×10⁶ cells/mL and CpG oligonucleotides to stimulate transformation. Plates were incubated at 37 °C and checked daily for the formation of LCLs. After 6–12 days, cell supernatants were harvested and screened for antibodies to PT in a plate-based ELISA, resulting in a positive, indeterminate, or negative result. Indeterminate results occur due to a lack of cell viability or failure to form LCLs. The RBF was calculated as the number of positive ELISA wells divided by the number of LCLs. Significant differences between the pertussis-specific B-cell responses were evaluated using a Wilcoxon rank-sum test.

RNA-Seq and ribosome profiling analysis

Raw sequence reads determined using RNA-Seq and RP were pre-processed by removing adapters, low-quality reads, and reads mapping to known rRNA and tRNA genes. For each sample, sequence reads were mapped to the latest human reference transcriptome/genome using HISAT2 (Version 2.2.0). Gene expression quantification was carried out using Subread (Version 2.0.1) against the latest Ensembl gene model reference collection (Version 100; April 2020). Systematic differences in sequence coverage between samples were normalized using the TMM method as implemented in the edgeR R package (Version 3.18.1). Post normalization, genes located on the X or Y chromosomes were excluded to avoid sex-specific effects. Moderated log counts per million (LCPM) for the remaining genes were determined using edgeR. LCPM was used to identify technical bias (e.g., batch effects or GC-content bias) and outlying samples (e.g., sample labeling error). More specifically, principal component analysis, non-metric multidimensional scaling, and correlation heatmaps based on LCPM were used to identify potential global outliers and systematic effects. Samples with outlying profiles based on both principal component analysis and non-metric multidimensional scaling (4 samples for RNA-Seq, 1 sample for RP) were considered global outliers. Reverse cumulative distribution function plots of average LCPM across samples were used to identify a suitable minimum average LCPM filter cut off to exclude lowly expressed genes.

Negative binomial models as implemented in edgeR were used to determine DE genes for each assay and vaccine group after removal of lowly expressed genes and outliers under the assumption that counts are negative binomial distributed. Models included fixed effects for subject, study visit day (pre-vaccination or post-vaccination day). The statistical significance of the study day effect (post – pre-vaccination) was evaluated for each gene using a likelihood ratio test. To compensate for testing multiple genes, p-values were adjusted by calculating false-discovery rates (FDR) based on the Benjamini–Hochberg procedure. Genes with an FDR-value < 0.05 and a fold change of ≥1.5-fold (up or down compared to pre-vaccination) were deemed to be significant DE genes. The following number of participants were included in each DE gene evaluation: n = 10 for RNA-Seq, Day 2, DTP, n = 9 for RNA-Seq, Day 2, DTaP, n = 9 for RNA-Seq, Day 8, DTP, n = 7 for RNA-Seq, Day 8, DTaP, n = 9 for RP, Day 2, DTP, n = 9 for RP, Day 2, DTaP, n = 7 for RP, Day 8, DTP, n = 7 for RP, Day 8, DTaP. Results for each DE gene including FDR-adjusted p-values are detailed in Supplementary Tables 9–15.

Pathway enrichment analysis was carried out separately for each differentially transcribed or translated gene set identified using published gene set information obtained from KEGG (KEGG Pathways Version 98.0; April 2021 and KEGG Modules Version 98.0; April 2021). The analysis was conducted using the goseq R package (Version 1.28) which adjusts for gene length bias when modeling the null distribution. To compensate for testing multiple gene sets, p-values were adjusted using the Benjamini–Hochberg procedure to control the FDR. Gene sets with an FDR-adjusted p-value < 0.1 were considered to be significantly enriched.

To estimate changes in cell type proportions that were present in bulk PBMC RNA-seq data, in silico deconvolution was performed using the R package granulator (Version 1.6.0). Deconvolution algorithms used for this analysis were based on the following methods: ordinary least squares, non-negative least squares, quadratic programming with non-negativity and sum-to-one constraint, quadratic programming without constraints, re-weighted least squares, support vector regression, and linear mixing model, as implemented via the granulator R package³⁸. Cell proportions were then estimated by combining the results from each deconvolution method using the median.

Regularized linear regression analysis in combination with bootstrapping was performed to gene expression fold change responses that predicted antibody response measurements using the glmnet R package (Version 2.0-13). Separate models were fit for each day and assay. To avoid overfitting (n«p and collinearity among gene responses) and to facilitate variable selection, an elastic net regularization step (combination of L1 Lasso and L2 ridge penalization, a = 0.5) was included as part of the fitting procedure. Cross validation was used to determine the optimum regularization parameter. The predictor gene variable set was based on log₂ fold change in LCPM and included genes with an average absolute baseline fold change of ≥1.5 in either treatment group at either of the two post-vaccination visits (Day 2 or 8). The analysis was performed for each combination of assay (RP and RNA-Seq), treatment group, day, and antibody response measurement (PT and FHA).

Regularized logistic regression models were fit to determine gene expression fold change responses that predicted reactogenicity group membership. Analysis was conducted as described for regularized linear regression analysis aside from (1) the absence of bootstrapping and (2) the predictor gene variable set included genes with an average absolute baseline fold change of ≥1.5 in either reactogenicity group at either of the two post-vaccination visits (Day 2 or 8).