To investigate the mutational space sampled by each germline V gene, we first compared the somatic mutations observed in human antibody repertoires of three healthy donors and three HIV-1-infected donors. Briefly, B cell receptor transcripts from peripheral blood B cells were sequenced using either Roche 454 pyrosequencing or Illumina MiSeq technologies. Starting from hundreds of thousands or millions of reads from each sample (Table S1), we assigned germline V and J genes for each transcript, removed low-quality sequences, identified clonal lineages, and selected one representative sequence per clonal lineage to build GSSPs for each V gene (Figures 1 and 2). We used the program partis (21) to predict novel germline alleles for all quality-filtered repertoires (see Materials and Methods), providing high confidence that all germline gene polymorphisms were excluded from the profiles.

In order to test the robustness of GSSPs to noise in the data, we subsampled the lineages found in each of the three healthy donors and built profiles for common VH genes using 25, 50, 100, 200, or 300 lineages per profile. We then calculated the Jensen–Shannon divergence between GSSPs (see Materials and Methods). We found that the between-donor Jensen–Shannon divergence between these profiles began to converge when 300 lineages were used to create a profile (Figure S1A in Supplementary Material). We therefore used 300 as the minimum number of lineages to build a mutational profile for further analyses except for quantifying the similarities of GSSPs (see next section), in which seven of 69 GSSPs built using 100–300 lineages were included but they did not change our conclusions. We also determined that there were no significant differences between GSSPs from IgM and IgG repertories, and therefore treated all VH data together. High-resolution plots of each profile, as well as the underlying numerical data, can be found at DOIs 10.6084/m9.figshare.3511083 and 10.6084/m9.figshare.3511085.

For each position of each germline V gene (numbered using IMGT scheme), we calculated the rarity of all possible substitutions as (1)rVi,a=1−(mV,i*fVi,a​) where m_V,i is the substitution frequency of position i in germline gene V and f_Vi,a is the substitution bias, or frequency at which a particular non-germline amino acid a was observed in all V gene lineages with substitutions at position i. The rarity is undefined for germline amino acids and equal to 1.0 for substitutions that are never observed in a particular dataset. When we examined the effects of choosing different representative sequences for each lineage (see Materials and Methods), the rarity scores of all substitutions were highly correlated regardless of representative sequence (Figure S1B in Supplementary Material), demonstrating again that the observed GSSPs are not significantly affected by noise in the data.

Figure 2 shows the GSSPs of IGHV1-69 and IGKV3-20 from three healthy donors in which the distributions of somatic mutation levels of lineages for each gene are similar (Figure S2A in Supplementary Material). Overall, the amino acid GSSPs of each V gene are remarkably consistent, similar to previous studies of mutations and selection (12, 18, 19, 21) In contrast, different V genes have noticeably divergent profiles, even within a single donor, as can be seen from the GSSPs of IGHV1-2 and IGHV1-69 (Figures S2B,C in Supplementary Material). The GSSPs depend on two factors, substitution frequency and substitution bias, each of which is dealt with separately in the following sections. We also combine these two to define rarity, shown in Equation 1.

The substitution frequency at each position of each V gene is similar among donors. As expected, the substitution frequency is higher within CDR 1 and CDR2 (Figure 2, as defined by IMGT), consistent with the fact that the CDRs contain higher proportions of AID hotspot motifs (18). For IGHV1-69, six positions in framework region (FWR) 3 exhibited substitution frequencies as high as positions in the CDRs (≥30%), probably because substitutions in the FWRs can play important roles in recognizing antigens and regulating antibody structural stability and conformations (28–31). Moreover, we also observed high frequencies of substitutions at positions adjacent to the CDRs (position 39 of IGHV1-69 and position 66 of IGKV3-20, IMGT numbering). Since these positions connect the CDRs to the β-strands of the FWRs, it is possible that these positions are important for regulating the flexibility and conformations of the CDR loops (28, 29).

For many positions, we observed 1–3 dominant substitutions with much higher frequency than other substitutions, indicating a substantial substitution bias. As expected, a major bias is toward substitutions that require only a single nucleotide change (Figure 3; Figures S3A,B in Supplementary Material). In addition, we found that amino acid substitutions requiring 2 or 3 nucleotide changes are significantly rarer than those which require only a single nucleotide change (Figure S3C in Supplementary Material). However, the exact substitutions that are preferred vary based on context, even for a specific codon (Figure 3). This is consistent with previous findings that mutation frequency varies by position within the V gene, independent of and in addition to variations due to SHM hotspots and coldspots (12, 17).

To quantitatively compare different GSSPs, we calculated weighted average of the Jensen–Shannon divergence between homologous positions over the entire V gene (see Materials and Methods). We then used multidimensional scaling (MDS) to visualize these distances for a subset of most frequently used V_H (Figure 4A) and V_κ genes (Figure 4B). (MDS plots of all V_H, V_κ, and V_λ genes are shown in Figures S4A–C in Supplementary Material). These plots confirmed that the GSSPs of each V gene from all donors clustered together and that the GSSPs of each V gene family are more similar than between V gene families. Moreover, we did not observe any differences based on HIV status or the sequencing technology used to obtain data. The GSSPs of a V gene are also similar across longitudinal samples of the same donor (Figures 4C,D; Figures S4D,E in Supplementary Material). Finally, the distributions of substitution frequency and rarity were both highly correlated between donors (Pearson’s r ~0.9 and ~0.8, respectively) (Figures 4E,F; Figures S5A,B in Supplementary Material). In order to increase our sampling depth, we therefore combined all lineages from the three healthy donors and generated a single set of GSSPs. We recalculated substitution frequency and rarity from these profiles, which were used for all further analyses.

The fact that similar preferred substitutions are detected in the repertoires of multiple donors could be explained by a limited range of possibilities. One possibility, though unlikely, is that observed GSSPs are dominated by convergent selection against common antigens. It is also possible that antigen-driven selection, while critical for the development of each individual lineage, is “averaged out” over the entire repertoire such that the observed GSSPs reflect the biased action of the underlying SHM machinery. Finally, it is possible that constraints on structural stability and other sources could limit the substitutional space explored by the antibody repertoire.

In order to distinguish between these possibilities, we first built a GSSP using lineages derived from non-productive rearrangements of IGHV3-23 (32). Because these sequences are derived from the “passenger” allele, they are not subject to selective pressure, even though AID continues to act on them (6, 33). We find that the GSSPs of the functional and non-productive repertoires are highly similar (Figure 5A), although mutations to cysteine observed at several CDR2 and FWR3 positions were suppressed in the GSSP of functional antibodies. We next compared the substitution frequency at each position of the two profiles and showed that they are as similar to each other as between functional antibody profiles of two donors (Pearson’s r = 0.90, p << 0.01) (Figures 5B and 4E). The same was true for rarity (Pearson’s r = 0.80, p << 0.01) (Figures 5C and 4F). While the residual differences between the functional and non-productive GSSPs are likely to reflect the effect of antigen-driven selection, the strong overall correlation suggests that antigen-driven selective effects are mostly averaged out when calculating GSSPs across the entire functional repertoire. Indeed, calculations of rarity based on two single lineages showed much lower correlations to their respective functional repertoires (Pearson’s r = 0.11, p = 0.11 for lineage 08248-00037 and r = 0.26, p = 1e-6 for lineage CAP256-VRC26) (Figure S6 in Supplementary Material), demonstrating that we can observe the modulation of substitution preferences for individual lineages by antigen-driven selection.

Previous work has found that, when site-specific selection coefficients for functional sequences are compared to those derived from non-productive sequences, approximately 30 and 5% of positions in the framework 3 region of human heavy chain genes are under negative and positive selection, respectively (34), which is roughly consistent between donors. McCoy et al further demonstrated that positions under negative selection tend to have less exposed surface area and proposed that the types and frequencies of substitutions at these sites may be restrained by negative selection for structural stability. To understand whether selection for structural stability modulates the similarities of substitution preferences observed in GSSPs, we compared substitution frequencies and rarity scores in the framework 3 region of IGHV3-23 from functional and non-productive repertoires at positions inferred as being under positive, negative, or neutral selection by McCoy et al (Figure S7 in Supplementary Material). The analysis revealed that the frequency of substitutions is less correlated between the functional and non-productive repertoires for positions under either positive or negative selection (Pearson’s r = 0.70, 0.63, 0.30 for neutral, negative, and positive selection respectively), suggesting that selection modulates the substitution frequencies observed in GSSPs. The analysis further showed that sites under positive selection showed reduced but still significant correlation between the functional and non-productive repertoires (Pearson’s r = 0.49 vs. r = 0.67 for neutral selection), suggesting selection modulates substitution bias in GSSPs. However, we were surprised to observe an increased correlation for sites under negative selection (Pearson’s r = 0.78), but the correlation is comparable to the measured similarity of the rarity scores for all V gene positions between the functional and non-productive repertoires (Pearson’s r = 0.79, Figure 5C). We note that only ~38 sites with estimates of selection pressure by McCoy et al are also present in the set of non-productive sequences used in this study, which may allow the introduction of sampling bias. Nonetheless, the sites under neutral selection showed a high correlation coefficient between the functional and non-productive repertoires, suggesting that there are other factors modulating substitution preferences.

We next sought to determine whether the substitution frequency and substitution biases observed in GSSPs are modulated by the biased action of SHM machinery. To accomplish this, we conducted simulations and generated a virtual repertoire of IGHV3-23 lineages using S5F (18), an antibody-specific mutation model that estimates a mutability and mutation preference for each nucleotide based on the four surrounding nucleotides (two on either side). The comparisons of the GSSPs showed that the simulated repertoire reproduces many dominant mutations in the functional and passenger allele repertoires, such as G10D and V89I/L (Figure 5A). Indeed, the mutability and rarity of the simulated repertoire was reasonably similar to substitution frequency and rarity of the functional and passenger allele repertoire, respectively (Figures 5B,C). Conversely, simulations using a previous antibody-specific amino acid substitution model (AB) (20) showed a heavy bias toward the use of histidine. The Pearson’s rho of the rarity scores derived from AB compared to those derived from the functional repertoires was very low (Figure S8 in Supplementary Material). This is likely due to a number of causes, including the omission of positional effects (17) and the mixture of sequences from different loci and species used to build the model.

Since the S5F model was designed to reflect intrinsic AID hotspots and coldspots (18), the strong correlation of rarity scores between simulated and actual data provides additional evidence demonstrating that the molecular machinery of SHM plays a key role in producing the observed GSSPs. Because the substitutions observed in the GSSPs recapitulate the majority of the sampled mutational space, we will use the GSSPs to understand the features of mutations available for antigen-driven selection in analyses in the following sections.

We next investigated the distribution of the observed rarity of all possible mutations. Strikingly, over 70% of all possible mutations were never observed in any of our datasets. Moreover, nearly 45% of observed mutations were “extremely rare” (defined as rarity > 0.995) (Figure 6A). Importantly, rarity is defined as a per-lineage score, as only one sequence per lineage is used to construct the GSSPs (Figure S1B in Supplementary Material). Thus, an extremely rare substitution is one that is expected to be sampled only by one in two 100 lineages. However, particular rare mutations can become fixed in highly expanded lineages (Figure S6 in Supplementary Material) and therefore be present in a larger fraction of individual B cells.

We then analyzed the substitution bias f_Vi,a for each substitution type at positions with substitution frequency m_V,i ≥ 0.05 averaged over all germline genes (Figure 6B). We used the within-position bias rather than rarity to control for differences in overall mutabilities. The results demonstrated that conservative substitutions tend to be favored, but we observed many non-conservative substitutions such as R to T and P to S. Interestingly, we also observed many substitutions that were asymmetric. For instance, R to T is more common than T to R; E to D is more common than D to E. To quantify the difference between conservative and non-conservative mutations, we compared the within-position bias for each observed substitution to its substitution score in the BC0030, a substitution matrix that scores amino acid changes based on functional interchangeability (35), which we used as a proxy for physicochemical similarity (Figure 6C). BC0030 only accounted for ~28% of the variation (i.e. Pearson’s r = 0.28) in within-position biases (Figure 6C), suggesting that substitutions in the antibody context are not solely constrained by physicochemical properties. In addition, the observed bias is different from a previous antibody-specific amino acid substitution model (AB model) (20), which showed high propensities for histidine to exchange with nearly all other amino acids.

While the overall distribution of SHM was similar in all three healthy donors (Figure S2A in Supplementary Material), we found that lineages with higher SHM were more likely to contain extremely rare mutations (Figure S9A in Supplementary Material). Indeed, when we subsampled IGHV1-2 lineages with either low SHM (8 or fewer amino acid substitutions; 854 lineages with a combined 3,193 substitutions from germline) or high SHM (more than 20 amino acid substitutions; 132 lineages with a combined 3,136 substitutions) from the combined dataset of the three healthy donors, we found that they resulted in substantially similar profiles (Figure S9B in Supplementary Material) and highly correlated rarity scores for substitutions that were observed in both sets of lineages (Figure S9C in Supplementary Material), after accounting for the smaller number of lineages in the high SHM set. However, the high SHM lineages contained more unique substitutions that were not observed in the low SHM set (Figure S9D in Supplementary Material). Although this may be expected, as most possible mutations are rare, it confirms that the elicitation of antibodies with specific rare mutations is expected to require higher overall levels of SHM.

Anonymized human PBMCs from normal, healthy donors were obtained through the NIH Clinical Center Department of Transfusion Medicine apheresis program by automated leukapheresis. Signed informed consent from the donors was obtained in accordance with the Declaration of Helsinki, and the study was approved by the National Institute of Allergy and Infectious Diseases (NIAID) Institutional Review Board. PBMCs were prepared by density gradient separation using Ficoll Paque Plus (GE HealthCare Life BioSciences, AB). Cells were then frozen in heat-inactivated fetal calf serum:DMSO (90:10) and stored at −185°C until needed.

Immunoglobulin genes were amplified from PBMC samples from three HIV- and hepatitis C-negative individuals as previously described in Ref. (25, 42–44). Briefly, human PBMCs (6 × 10⁷) were previously obtained from three HIV-1 and hepatitis C-negative individuals (LP32647, LP08248 and LP23810), and the PBMCs were pelleted at 1,200 rpm for 8 min. mRNA was then extracted and eluted in 50 µl elution buffer using μMACS mRNA isolation kit (Miltenyi Biotec) according to the manufacturer’s instructions. The mRNA was then aliquoted for cDNA synthesis, and a 5′RACE approach was used to amplify IgG genes from one aliquot from each donor. The datasets were published and deposited in the NCBI Short Reads Archive (SRA ID: SRP067168) (25).

In this study, we prepared additional IgG libraries from aliquots of the same mRNA using 5′ RACE, as well as an IgG library using 5′ multiplex primers. In addition, we generated libraries for IgM, IgK, and IgL Ig genes using both 5′RACE or 5′multiplex methods. 5′RACE was conducted based on previously described methods (42). Briefly, to synthesize cDNA, 10 µl mRNA was mixed with 1 µl 5′CDS Oligo dT primers (12 µM) and incubated at 70°C for 1 min and then −20°C for 1 min. Then, 1 µl SMARTER Oligo Primer (12 µM) (Clontech), 4 µl 5× RT buffer, 1 µl DTT 20 (20 mM), 1 µl dNTP (10 mM), 1 µl RNAse out, and 1 µl SuperScript II reverse transcriptase (Invitrogen) were added to the reaction. After 2 h of incubation at 42°C, the cDNA products were purified using Nucleospin II kit (Macherey-Nagel) and eluted in 50 µl water. For Ig gene recovery, 10 µl cDNA was used for each PCR reaction. The first PCR amplification was performed with a common 5′ primer II A (Clontech) and Ig constant region-specific 3′ primer (IgG: 5′GGGGAAGACCGATGGGCCCTTGGTGG3′; IgM: 5′GAGGGGGAAAAGGGTTGGGGCGG3′, IgK: 5′GGAAGATGAAGACAGATGGTGCAGCCACAG3′, IgL: 5′CCTTGTTGGCTTGAAGCTCCTCAGAGGAGG3′) using KAPA HIFI qPCR kit (Kapa Biosystems). The PCR products were purified with 2% Size Select Clonewell E-gel (Invitrogen) and Agencourt AMPure XP beads (Beckman Coulter). The second PCR amplification was performed with primers with 454 sequencing adapters (454-RACE-F: 5′CCATCTCATCCCTGCGTGTCTCCGACTCAGAAGCAGTGGTATCAACGCAGAGT3′; 454-IgG-R: 5′CCTATCCCCTGTGTGCCTTGGCAGTCTCAGGGGGAAGACCGATGGGCCCTTGGTGG3′; 454-IgM-R: 5′CCTATCCCCTGTGTGCCTTGGCAGTCTCAGGAGGGGGAAAAGGGTTGGGGCGG3′; 454-IgK-R: 5′CCTATCCCCTGTGTGCCTTGGCAGTCTCAGGGAAGATGAAGACAGATGGTGCAGCCACAG3′; 454-IgL-R: 5′CCTATCCCCTGTGTGCCTTGGCAGTCTCAGCCTTGTTGGCTTGAAGCTCCTCAGAGGAGG3′). The PCR products were again purified with 2% Size Select Clonewell E-gel and Agencourt AMPure XP beads. For 5′ multiplex method, 10 µl mRNA was used for cDNA synthesis as described previously (43). cDNA was eluted in 50 µl water, and 10 µl cDNA was used for each PCR reaction. As described previously (43), IgG and IgM genes were amplified using mixed primers for all VH families, while IgK and IgL genes were amplified using mixed kappa and lambda primers, respectively. PCR products were purified with 2% Size Select Clonewell E-gel and Agencourt AMPure XP beads.

454 pyrosequencing for all libraries was performed as described previously (45), the datasets were deposited to NCBI SRA database (Bioprojects: PRJNA336331). The next-generation sequencing data for the antibody repertoires of the HIV infected donors NIH45, CAP256, and Z258 were published in Ref. (44) (SRA ID: SRP052625) (43, 46), (SRA ID: SRP034555 and SRP017087), and (Huang et al. Immunity 2016, Accepted) (SRA ID: SRR4417615-SRR4417632), respectively. The non-productive sequences derived from IGHV3-23 were retrieved from the European Nucleotide Archive (accession numbers AM076988-AM083316) (32).

In order to ensure that our calculations were not affected by sequencing error from the 454 platform, we used a stringent quality control protocol as described in the following sections. We used SONAR² to assign germline V and J genes for each transcript (27). Sequences for which assignments could not be made were removed. In addition, transcripts contained stop codons or out-of-frame junctions were also removed, as these are likely to be the result of sequencing error. Because members of an expanded clonal lineage are likely to share substitutions which did not arise independently, we included only one transcript from each lineage. Lineages were found using SONAR and defined as having the same V and J gene assignments and CDR3 sequences of the same length and at least 90% nucleotide identity. Lineages containing only a single transcript were discarded, as it is impossible to verify the quality of these sequences. For all lineages containing multiple sequences, the transcript closest to the consensus sequence of the lineage was chosen as the representative, which helps minimize the effects of sequencing error (44). In addition, because the primary mode of 454 error is indel, the remaining transcripts were each individually aligned to the assigned germline V gene with CLUSTALW to check for frameshift error. Any sequences with non-codon-length indels were discarded. Finally, each transcript was translated and compared to the amino acid sequence of the assigned germline V gene, and any transcripts with no amino acid changes were discarded, as they do not contribute any information to the GSSP. Paired-end reads obtained from donor Z258 using the Illumina MiSeq platform were merged into a full-length amplicons using USEARCH (47), discarding sequences containing > 25 mismatches in the overlap region. Merged reads with two or more expected errors (based on PHRED scores) were also discarded, and the remaining reads were processed in the fashion described earlier. Overall, we processed nearly 38 million NGS reads across 54 samples from six donors, which resulted in 14.6 million high-quality reads and over 200,000 lineages which we used to construct GSSPs (Table S1).

Donor LP08248 had 284 VH3-23-derived lineages with at least 10 member sequences. We conducted 100 trials in which a random sequence was chosen as the representative of each lineage, and a mutability profile was constructed as described in the following sections. For each of these profiles, we calculated a set of rarity scores. We then calculated the Pearson’s correlation coefficient between the rarity scores of each pair of profiles; the distribution of these coefficients is shown in Figure S1B.

We retrieved 348 VH3-30-derived sequences from the CAP256-VRC26 broadly HIV-1-neutralizing antibody lineage (46), including 33 experimentally isolated monoclonal antibodies. We constructed a GSSP from these sequences and compared the rarity scores derived from this profile to rarity scores derived from a VH3-30 profile constructed from the merged lineages of all three HIV⁻ donors (Figures S6C,D in Supplementary Material). For comparison, we also constructed a profile from 378 VH3-23 derived sequences from a single lineage found in donor LP08248 (Figures S6A,B in Supplementary Material).

Identifying mutations requires an accurate knowledge of the germline sequence. Undocumented polymorphisms can appear as donor-specific, high-frequency mutations, skewing the resulting GSSPs and exaggerating differences between donors. We therefore used partis (21) to infer novel germline alleles from each non-redundant dataset. (Allele finding is considered a beta feature, available from https://github.com/psathyrella/partis. We used commit 610ae46 from June 18^th, 2016.) Across all 6 donors, we found 33 unique novel heavy chain alleles, 23 novel kappa chain alleles, and 9 novel lambda chain alleles. Of these, 3, 5, and 4, respectively, were observed in two or more of these donors. The alleles have been deposited in the IgPdb³.

For each non-redundant dataset, we used CLUSTALW2 to align translated transcripts from each V gene to all alleles of that gene, including novel alleles detected by partis. For each position, we calculated the substitution frequency and rarity as described in the text. For polymorphic positions, all possible polymorphisms were considered germline residues, without confirming the presence of each allele in each donor. Gaps and undefined amino acids were excluded from counting. GSSPs were only calculated if there were at least 100 sequences for a particular V gene in the given dataset. When building GSSPs, IMGT alleles IGHV4-4*07 and IGHV4-4*08 were considered as alleles of the IGHV4-59 gene, due to sequence similarity. Logo plots of GSSPs were shown using a customized version of WebLogo (48).

To compare two GSSPs, we calculated the Jensen–Shannon divergence between the distributions of mutations or substitutions in the two profiles at each position in the V gene. This is defined as J(PA,i,PB,i)=H(PA,i+PB,i2)−H(PA,i)+H(PB,i)2 where P_A,i and P_B,i are vectors of length 20 describing the probability of observing each possible non-germline residue at position i for profiles A and B, respectively. H is the Shannon entropy, defined as H(Pi)=−∑j={aa}Pijlog2Pij where j can be any of the 20 amino acids. These position-wise Jensen–Shannon divergences were then averaged, with the divergence of each position weighted by the mean substitution frequency for that position in the two profiles being compared: 〈JA,B〉=∑imA,i+mB,i2J(PA,i,PB,i)∑imA,i+mB,i2 where m_A,i and m_B,i are the substitution frequency (observed frequency of mutation) at position i for profiles A and B, respectively. For comparisons of GSSPs from different V genes, only homologous positions were used. We generated a matrix of divergences among the datasets studied and used the cmdscale command in R to do multidimensional scaling and generate coordinates for plotting.

The non-productive sequences were manually aligned to germline gene IGHV3-23 to remove indels in order to produce a meaningful GSSP. In cases where the boundaries of an indel were ambiguous, the boundary bases were replaced with Ns to avoid potential bias.

Starting from the nucleotide sequences of IGHV3-23*01, the S5F model was used to generate an artificial repertoire of 420 lineages, each with an independently simulated set of mutations. We produced 15 sequences each with 1–28 nucleotide changes; after accounting for the possibilities of silent mutations or multiple changes in a single codon, this produced a distribution of amino acid changes that roughly matched observed distributions of SHM in the functional repertoire (Figure S5C in Supplementary Material). For each sequence, we calculated the mutability of each nucleotide position under the S5F model and randomly selected a site for mutation based on those probabilities. The resulting substitution was then randomly chosen based on the substitution biases indicated for the target base by the S5F model. We excluded mutations that resulted in stop codons, but did not exclude re-mutation at a previously mutated site. For simulations with the AB model, positions to be mutated were selected using our own observations of mutation frequencies as calculated for the GSSP of IGHV3-23 and the substitution was then chosen under the AB model.

The GSSPs of 69 V genes (26 heavy chain, 20 kappa chain, and 23 lambda chain) were used to calculate the substitution frequencies. Because the substitution frequency at positions with low substitution frequency may be undersampled, we excluded positions being mutated in less than 5% of the V gene-specific clones. We also excluded positions containing amino acid changes among alleles. For each selected position of a profile, we normalized the substitution frequency to 1. Finally, all selected positions were grouped based on germline amino acid type. For each of the 20 amino acids, the substitution frequency to each of the other amino acids was averaged over all positions within a group.

The Mann–Whitney U test was used to evaluate the similarities of features between different groups in this study. Pearson’s correlations and statistical significance were calculated using cor.test in R.