Epidemiologic Evaluation of Human Papillomavirus Type Competition and the Potential for Type Replacement Post-Vaccination

Background Millions of women have been vaccinated with one of two first-generation human papillomavirus (HPV) vaccines. Both vaccines remain in use and target two oncogenic types (HPVs 16 and 18); however, if these types naturally compete with others that are not targeted, type replacement may occur following reductions in the circulating prevalence of targeted types. To explore the potential for type replacement, we evaluated natural HPV type competition in unvaccinated females. Methods Valid HPV DNA typing information was available from five epidemiological studies conducted in Canada and Brazil (n = 14,685; enrollment across studies took place between1993 and 2010), which used similar consensus-primer PCR assays, capable of detecting up to 40 HPV types. A total of 38,088 cervicovaginal specimens were available for inclusion in our analyses evaluating HPV type-type interactions involving vaccine-targeted types (6, 11, 16, and 18), and infection with each of the other HPV types. Results Across the studies, the average age of participants ranged from 21.0 to 43.7 years. HPV16 was the most common type (prevalence range: 1.0% to 13.8%), and in general HPV types were more likely to be detected as part of a multiple infection than as single infections. In our analyses focusing on each of the vaccine-targeted HPV types separately, many significant positive associations were observed (particularly involving HPV16); however, we did not observe any statistically significant negative associations. Conclusions Our findings suggest that natural HPV type competition does not exist, and that type replacement is unlikely to occur in vaccinated populations.


Introduction
Infection with high-oncogenic risk human papillomavirus (HR-HPV) is a necessary cause of cervical cancer in women [1] and an important cause of other anogenital cancers in both genders [2]. In addition, some low-oncogenic risk (LR) HPV infections may cause benign lesions known as acuminate condylomata (genital warts), as well as low grade squamous intraepithelial cervical lesions. Two highly effective HPV vaccines have been administered to millions of women around the world (Merck's Gardasil© and GlaxoSmithKline's Cervarix©) [3,4], offering protection against two HR-HPV types (16 and 18)-responsible for approximately 70% of cervical cancer cases. Only Gardasil protects against additional LR-HPV types (6 and 11) that cause approximately 90% of genital warts cases [5][6][7]. Although HPV vaccination is eventually expected to reduce the burden of disease attributable to these HPV types, there is concern that it may lead to "type replacement" [8], i.e., an increase in the prevalence of other non-vaccine HPV types following the reduction of vaccine-targeted types [9,10].
For type replacement to occur, a biological prerequisite is that different HPV types must compete with one another for niche occupation during natural infection [9][10][11]. We recently described different epidemiological approaches to evaluate HPV type competition in order to gain insight regarding the likelihood of type replacement [10]. The two main approaches include construction of Kaplan-Meier curves and Cox models to evaluate sequential acquisition and clearance of HPV types according to HPV status with vaccine-targeted types; and construction of logistic regression models for each vaccine-targeted type to explore whether infection with these types may be associated with infection by other HPV types. A number of cohort studies evaluating the natural history of HPV infections among females have suggested that those infected with HPV (any type) are generally at higher risk of acquiring other types [12][13][14][15], or at about equal risk of acquiring and clearing existing infections [12][13][14][15][16][17]. Similarly, other recent cross-sectional studies that have investigated clustering patterns of different HPV types have found that females infected with HPV (vaccine or other types) are more likely to be infected with additional HPV types [18][19][20][21][22][23][24][25][26]. These previous studies reported very few negative associations, therefore providing some reassurance that type competition does not exist and that replacement is unlikely. Despite the large sample size of some of these studies, few or no co-infections were observed for rare HPV types, leading to non-positivity or low precision for some comparisons. In addition, evaluation of pairwise interactions in these studies did not account for presence of other HPV types, which may have introduced some confounding [10].
To evaluate HPV type competition in the current study, we applied a hierarchical (Bayesian) regression approach that employs shrinkage and adjustment for confounders, as well as other HPV types. Data were available from five pre-vaccination studies conducted among females in Canada and Brazil. . Protocols for each of the five studies have been described in detail elsewhere [27][28][29][30][31]. Briefly, the three cohort studies (Ludwig-McGill, HITCH, and McGill-Concordia) were designed to evaluate the natural history of HPV infection among females, and transmission of HPV among heterosexual couples (male data from the HITCH study was not included in the current analysis). BCCR is a case-control study that was originally designed to evaluate the role of biomarkers in the etiology of cervical precancer and cancer, and CCCaST was the first North American randomized controlled trial to compare Pap cytology versus HPV testing in screening for cervical cancer. Subjects completed questionnaires to collect information on important demographic and lifestyle variables; and provided cervical samples (self or provider collected) for HPV testing at each of their clinic visits. All participants provided written informed consent and each study was approved by review boards or ethical committees at McGill University and other participating institutions.

HPV DNA detection and genotyping
In the three cohort studies, cervical specimens were collected and tested for HPV at each clinic visit (every four months during the first year of follow-up/twice annually in subsequent years of follow-up in the Ludwig-McGill and HITCH studies; and twice annually in the McGill-Concordia study). Subjects from the Ludwig-McGill, HITCH, and McGill-Concordia studies contributed an average of 9.0, 4.4, and 4.2 cervical specimens for HPV testing, respectively; whereas subjects from the BCCR and CCCaST studies contributed only one specimen for HPV testing.
Details regarding specific sample collection and HPV testing protocols for each study have been described in detail elsewhere [27][28][29][30][31]. Briefly, all studies employed consensus primer PCR assays (L1 PGMY or MY09/11 and hybridization with oligonucleotide probes and restriction fragment length polymorphism analysis, line blot assay, or linear array), which are capable of detecting between 27 and 40 different HPV types. The MY09/11 and PGMY09/11 protocols are both very sensitive with good overall agreement (kappa range = 0.68-0.83) [32-34] and modifications to the MY09/11 protocol (leading to the PGMY09/11 protocol) has resulted in even greater test sensitivity [32]. Although the genotyping procedure in the Ludwig-McGill study (hybridization with individual oligonucleotide probes and restriction fragment-length polymorphism analysis) did not allow us to distinguish between vaccine-targeted HPV types 6 and 11, these are two of the most closely related HPV types (with similar biological and pathological properties) [35], therefore grouping them was not viewed as a major limitation. Nonetheless, we evaluated HPVs 6 and 11 together, as well as separately in the other four studies. Since types that are phylogenetically related (i.e., from the same species) share a large proportion of their nucleotide sequence (!60%) and display similar biological properties, we suspected that types from the same species would be more likely to compete [35,36]. HPV types belonging to the same species as HPV6/11 (α-10) include 13, 44, and 74; as HPV16 (α-9) include 31, 33, 35, 52, 58, and 67; and as HPV18 (α-7) include 39, 45, 59, 68, and 70.

Statistical analysis
We investigated the association between infection with the vaccine preventable types and infection with each of the other HPV types using pooled data from the five studies. Bayesian hierarchical regression models were constructed for vaccine preventable types 6, 11 (6/11 combined), 16, and 18. Age and lifetime number of sex partners were chosen as covariates a priori, since they are strong predictors of HPV infection [2]. Thus the primary analyses excluded a portion of CCCaST participants who were missing baseline data on lifetime number sex partners. Models for 6/11 combined, 16, and 18 included data from all five studies. Models for 6 and 11 separately excluded the Ludwig-McGill study, as explained above. Secondary analyses included the CCCaST participants with missing information on lifetime number of sex partners by excluding it as a covariate. We also conducted analyses for each study separately.
Specifically, the probability of infection with the vaccine preventable type was modeled in a 2-tier hierarchical model, where subjects' study visits were nested within subjects in order to account for subject-level clustering. At the visit level, a logistic model was fitted with infection with the vaccine preventable type as the outcome and every other HPV type and age at the time of the visit as predictors. At the subject level, the subject-specific intercepts were modeled by accounting for lifetime number of sex partners at baseline, as well as the study that the subject came from for the pooled models. Thus, the odds ratio (OR) estimate for each HPV type represents the odds of detection of the vaccine-preventable type in the presence of that HPV type compared to the odds of detection of the vaccine-preventable type in the absence of that particular HPV type, adjusted for all other HPV types, age at visit, lifetime number of sex partners at baseline, and study.
In order to improve the precision of the estimates for the effect of the presence of other HPV types on the presence of the vaccine preventable type, the logistic regression parameters for all the other HPV types were assumed to be normally distributed around an overall mean effect of co-infection. Diffuse or wide prior distributions were used for all other parameters. All analyses were conducted using WinBUGS software version 1.4.3 (MRC Biostatistics Unit, Cambridge).
The additional hierarchical component on the coefficients of other HPV types produces a shrinkage effect, whereby unstable estimates with large variances are drawn closer to the mean.
The assumption introduces a bias in favour of reducing variance and potentially reducing mean squared error [37]. To explore the possible effect of this bias, we also compared our results with estimates for HPV type associations calculated using the maximum likelihood method.

Results
Subject characteristics stratified by study population are listed in Table 1. The average age of participants at enrollment across the five studies ranged from 21.0 (HITCH study) to 43.7 years (CCCaST study). Given that they were studies of young adult women, HITCH and McGill-Concordia studies included few females that were married/common-law (14.1% and 18.0%, respectively) or that had ever been pregnant (9.8% and 16.2%, respectively). Compared with subjects from the four Canadian studies, Brazilian Ludwig-McGill study participants reported fewer lifetime sexual partners (87% had less than five partners) and the majority rarely used condoms (less than 4% used condoms regularly). Most subjects in the McGill-Concordia, HITCH and BCCR studies indicated that they were never smokers (62.7%, 62.3% and 50.0%, respectively); whereas the majority of Ludwig-McGill and CCCaST participants reported that they were current/former smokers (52.5% and 79.8%, respectively). Across all studies, HPV16 was the most common type detected among cervicovaginal specimens: Ludwig-McGill (n = 546, 2.5%), McGill-Concordia (n = 220, 8.2%), HITCH (n = 305, 13.8%), CCCaST (n = 105, 1.0%), and BCCR (n = 47, 4.8%) (Figs 1, 2, 3, 4 and 5). Although the ranking of other common HPV types varied across the studies, the majority were detected as part of a multiple infection (rather than as single infections), except in the Ludwig-McGill study. Subject characteristics that were commonly associated with multiple HPV infection included younger age and higher number of sexual partners ( Table 2). CCCaST participants who reported condom use ("ever" versus "never") and who were widowed/divorced were at higher risk of being infected with multiple HPV types, whereas subjects from the BCCR study who were married/common-law were at significantly lower risk compared with single individuals. Former smoking status was also associated with greater risk of multiple infections in HITCH and CCCaST studies, but not in the others. Figs 6 to 10 display results from the logistic regression models. Each of the graphs present OR estimates for type-type associations on the natural log scale; therefore, (log)OR estimates greater than zero correspond to ORs greater than one (i.e., positive associations between HPV types), and the opposite for (log)OR estimates below zero. In our pooled regression analyses (including data from all five studies), no statistically significant negative associations were observed between vaccine-targeted HPV types (HPVs 6, 11, 16, and 18) and any other types (Figs 6,7,8,9 and 10). In fact, the only point estimate indicating a negative association observed was between HPV18 and 89 (OR = 0.92, 95%CI: 0.49-1.52); however, there was insufficient precision to reject the null hypothesis of no association. These analyses included adjustment for other HPV types, age and lifetime number of sexual partners, but excluded over half of CCCaST study participants (n = 5754) due to missing sexual history information from St. John's study site participants. In our analyses adjusted for other HPV types and age only (including all CCCaST subjects), results were similar, i.e., no negative associations were observed, and OR estimates were generally higher (S1 Fig).
Across the studies with individual typing information for HPVs 6 and 11 (i.e., all other than Ludwig-McGill study), HPV11 was detected in only 23 of 16027 specimens. In our analyses of

Discussion
Assessment of pre-vaccine epidemiological data can provide insights concerning natural HPV type competition and the potential for type replacement [10]. HPV types that naturally compete with HPVs 6, 11, 16, and/or 18 may be more likely to fill the ecological niches vacated by these vaccine-target types. The US Food and Drug Administration and Health Canada recently approved Merck's new HPV vaccine (Gardasil 9©) that protects against the same four HPV types as the original Gardasil vaccine (6,11,16, and 18), plus additional oncogenic HPV types 31, 33, 45, 52, and 58 [38]. However, despite the availability of this new nonavalent vaccine, concern about type replacement remains important. Millions of women have already been vaccinated using either the bivalent or quadrivalent formulations, and both first-generation vaccines continue to be administered in many countries.
In general, our results support previous studies, which mainly reported null or positive associations between different HPV types [18][19][20][21][22][23][24][25][26]. Recently, Vaccarella and colleagues used a number of large data sets to evaluate clustering patterns between HPV types (via hierarchical regression models with woman-level random effects), identifying few negative associations and some positive associations, which they generally attributed to diagnostic artifacts [21][22][23][24]. Similarly, Chaturvedi and colleagues reported very few negative or positive associations in examining HPV co-infection patterns among women from the Costa Rica Vaccine Trial, concluding that HPV infections seemed to occur independently in this population [19]. Furthermore, in a recent pooled analysis, including information from three diverse study populations in the Netherlands, Mollers and colleagues also reported no significant pairwise interactions, but did suggest that clustering patterns differed across risk groups and across types, particularly between low-and high-risk HPV types [25]. In general, phylogenetic relatedness did not strongly influence clustering patterns in these prior studies; whereas in our study, HPV16 (α-9) was positivity associated with all related types, and HPVs 6/11 and 18 were positively associated with related types 44 (α-10) and 59 (α-7), respectively.
Across the five studies, there were more than 38,000 cervical specimens with valid HPV testing results, which makes the current pooled analysis one of the largest studies on this topic to date. As a result, we were able to evaluate associations between vaccine-targeted HPV types with all others, including rare types. The application of Bayesian methods incorporating shrinkage further improved our precision, and still allowed us to adjust for all relevant covariates and presence of other HPV types in our models. However, any improvement in precision resulting from shrinkage comes at the expense of introducing some bias [37]. To explore if our results may have been meaningfully different according to traditional analytic methods (i.e., without this bias/precision trade-off), we performed sensitivity analyses using maximum likelihood estimation. As expected, this approach led to wider confidence intervals, but importantly it did not lead to any statistically significant ORs less than one (data not shown). Although we did not observe any statistically significant negative associations in our study, we did observe a large number of positive associations. We suspect that most significant positive associations may be attributed to residual confounding, i.e., due to our inability to control for all risk factors of multiple-type HPV infections (e.g., host susceptibility, immunological differences, or other unmeasured behaviour risk factors). For example, in our analyses including all CCCaST specimens (i.e., unadjusted for sexual history), confounding may explain the higher OR estimates and greater number of HPV types found to be positively associated with HPVs 6, 11, 16       [39][40][41]. However, albeit insightful, this approach leads to a form of selection bias, referred to as collider stratification bias [42], and was therefore not applied in the current study. Despite variation in key demographic and behavioural risk factors across individual studies, results were generally consistent after adjustment for age, lifetime number of sexual partners, and other HPV types. Nonetheless, it is important to consider how differences in important HPV risk factors may have impacted our results and ability to pool information. For example, in the HITCH and McGill-Concordia studies, participants were younger than those in the other studies and therefore we may suspect that infections in these two studies are more likely to represent incident or recently acquired infections rather than persistent infections. This may have important implications since oncogenic vaccine-targeted types, such as HPVs 16 and 18, are more likely to be persistent and detected with other HPV types, leading to higher OR estimates in the current study. In a separate recent analysis conducted to evaluate incidence and clearance of individual HPV types according to infection with vaccine-targeted HPV types, we observed similar two-year incidence rates (any infection) in the Ludwig-McGill, McGill-Concordia, and HITCH cohort studies (23.6%, 27.0%, and 18.3%, respectively) [15]. Also, compared with their younger counterparts in the McGill-Concordia and HITCH studies, Ludwig-McGill participants were more likely to clear their existing infections within two years [15]. These results suggest that infections observed across studies in the current analysis may represent a similar proportion of incident or recently acquired and persistent infections. Importantly, results from this cohort analysis also did not provide any evidence of HPV type competition, i.e., individuals with vaccine-targeted HPV types were not less likely to acquire other types or more likely to clear their existing infections [15]. The five studies from which specimens were collected all utilized broad spectrum PCR assays to test for the presence of HPV. Although these assays are able to amplify and detect a large number of HPV genotypes and may detect as few as 10 copies of viral DNA for most common genital HPV types [32, 43,44], previously we discussed concerns regarding the sensitivity of consensus assays in the context of type replacement evaluation, particularly in situations where specimens are coinfected with multiple HPV types [10]. In a recent analysis conducted to evaluate possible "masking" of HPV52 in the presence of HPV16, we observed a significant positive association between HPV16 viral load and masking of HPV52 [45]. Other PCR assays have been developed with reported high sensitivity for detection of multiple HPV types from coinfected specimens, e.g., using array primer extension (APEX) for typing [46]; however, these methods remain less common. In addition, there is also the possibility that assay specificity may be reduced as a consequence of probe cross-reactivity [47], which may explain the tendency for some phylogenetically related types to cluster together. However, considering that most HPV types from the α-9 species are also classified as definite carcinogens by the International Agency for Research on Cancer (all except for HPV67) [48], they are also more likely to persist (than low-risk types) and therefore are more likely to be detected together with other types [25]. The observation that certain HPV types (e.g., HPV16) were consistently observed more frequently than others across individual studies suggests that a competitive advantage exists for some HPV types.
Previous cross-sectional and cohort studies focusing on different populations and employing unique analytic/genotyping methods have failed to provide consistent or strong evidence that negative pairwise HPV interactions exist [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26]. The current study adds to this literature by providing additional reassurance that-owing to the lack of HPV type competition-type replacement appears unlikely. Since we did not include females who received prophylactic HPV vaccines for comparison in this study, we must assume that no major differences in acquiring other types exist among females who are naturally uninfected with vaccine-target types. Eventually, a definitive answer to this question of whether HPV type replacement has occurred will come from long-term surveillance studies which compare pre-and post-vaccination type-specific HPV prevalence rates, and which properly account for possible diagnostic artifacts [10,45].