Epidemiological evaluation of human papillomavirus genotypes and their associations in multiple infections

The aim of this study was to determine the frequency of multiple type human papillomavirus (HPV) infections, and whether any types are involved in multiple HPV-type infections (mHPV) more or less frequently than expected. From January 2012 to February 2018, 2848 cervico-vaginal swabs were analysed in the UOC Microbiology and Virology of Policlinico of Bari, Italy. HPV DNA detection was performed using initially nested-polymerase chain reaction (PCR) and subsequently multiplex real-time PCR assay. 1357/2848 samples (47.65%) were HPV DNA positive and 694/1357 (51.14%) showed mHPVs. The median number of mHPVs was 2 (interquartile range: 2–3). HPV-types more frequently detected were 42 (9.97%), 16 (8.92%), 53 (7.23%) and 31 (7.16%). Each detected HPV-type was involved in mHPVs in more than 50% of cases. Statistical analysis showed significant associations for all HPV-types except for 33, 43, 51, 58 and 82 HPV-types. The major number of significant pairwise associations were detected for the types 42 and 70. Only positive associations were detected. Further data are necessary to evaluate the clinical impact of the single combinations.


Introduction
Human papillomavirus (HPV) infection contributes to the burden of diseases as genital warts and cancer, in both men and women [1][2][3]. In particular, nowadays, different studies have demonstrated that there is a high incidence of cervical cancer (CC) by a co-infection of different HPV-types. Thus, HPV co-infections might be considered a risk factor for the carcinogenesis, although different studies reached different conclusions, in which a correlation between multiple HPV-type (mHPV) infections and the carcinogenesis is not always demonstrated [4,5]. Moreover, it is not clear if the risk for the progression to cancer is correlated to some HPV-type specific combinations or if it is common for all the genotypes [6][7][8][9][10]. Some infections, however, may persist for many years, causing cellular alterations, which, if untreated, have the potential to develop into cancer of ano-genital area (carcinoma of the cervix, vagina, vulva, anal, oral and penile areas) [11,12].
Recently, it has been reported that mHPV infections are apparently more widespread than expected [13,14]. Although the genital infections with mHPVs are common in either men or women [1,[15][16][17], the relative rates might not be representative for the prevalence due to the limited number of patients studied [13,14,18]. Nevertheless, since when the HPV-types have been performed by polymerase chain reaction (PCR) assays, mHPV infections have been more consistently described [19]. Among the studies, relevant differences in mHPV prevalences have been described. For example, Dickson et al. have reported that in large studies performed in Costa Rica and Italy, the mHPV infection prevalence range from 24.8% to 52.6% among all HPV-positive tests [5]. Such differences may be explained by the genotyping techniques used and also by the type of and behavioural characteristics of the populations analysed. Furthermore, epidemiological studies have found that patients infected with mHPVs require a management to prevent the progression for cancer, underlining that the monitoring of the genotypes and of the viral load might provide some critical information in the diagnosis of infections either for future management or prognosis [20][21][22]. Nevertheless, the relationship between mHPVs and CC is still debatable.
In the past several studies had attempted to assess patterns of mHPVs in women [1,5,20]. However, different statistical models had been used to analyse the data, obtaining results either different or hardly comparable.
Moreover, mHPVs were seldom detected using earliest diagnostic methods and only the shift towards more sensitive assays allowed confirming that mHPVs are either highly widespread or very common among sexually active women [19]. Moreover, some authors suggested that the apparent clustering of mHPVs may be an artefact due to the measurement process (i.e. diagnostic assay) raising the question of the cross-reactions as confounding factor [1].
The aim of this study was to determine the prevalence of mHPVs in a cohort of women and also evaluate the association of pairwise HPV-types accounting for the analytical genotyping method.

Study population
From January 2012 to February 2018, 2848 cervico-vaginal swabs from 2848 women were consecutively collected. In the case of multiple samples from the same patient, only the first was retained for the analysis. Median age of women analysed was 37 years. The specimens were transferred to the laboratory of Molecular Biology, UOC Microbiology and Virology, Azienda Universitaria-Ospedaliera, Policlinico of Bari, where they were analysed.
All procedures performed in studies involving human participants were in accordance with the 305 ethical standards of the institutional and/or national research committee and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards. Sample information (date of sampling, ward, type of specimen, testing results) together with the data of patients for whom molecular testing was performed (i.e. age and sex) were recorded in an anonymous database by changing sensitive data into alphanumeric codes. No clinical data associated with these specimens were available. The study was approved by Ethical Committee of Policlinico of Bari (No. 0022174|12|03|2018).

Treatment of specimens
In total 2 ml of phosphate-buffered saline (pH 7.4) (Sigma-Aldrich, Milano, Italy) were added to the cervico-vaginal swabs, collected by a rigid cotton-tipped swab applicator (Nuova Aptaca, Cannelli, Italy), and vortexed. Subsequently, the samples were transferred to microcentrifuge tubes (2 ml) and they were stored at −20°C until processing.
At processing moment, microcentrifuge tubes were centrifuged at rcf = 15 700 g for 15 min at 7°C. The majority of supernatant was discarded but 200 µl of supernatant was retained to resuspend the pellet.

DNA isolation and nested-PCR
From January 2012, DNA extraction was performed by automated QIAcube System (Qiagen, Germany) following the manufacturer's protocols. Extracted DNA samples were submitted to a nested PCR amplification, using an Ampliquality HPV-HS Bio Kit (AB Analitica, Padova, Italy) following the manufacturer's instructions.
The method provides for a first amplification of the viral genome L1 region, followed by a second amplification phase with biotinylated primers. PCR products were analysed using 3% agarose gel electrophoresis with ethidium bromide to display DNA under ultraviolet light. Subsequently, PCR products were typed by using a Reverse Line Blot Hybridisation Ampliquality HPV-Type Kit (AB Analitica, Italy). The Ampliquality HPV-Type Kit allows to detect 18 HR-HPVs types (types 16 To assess the suitability of extracted DNA, thiosulfate sulfurtransferase (TST) gene region (202 bp) was amplified at the first amplification. A negative result TST amplification underlines the presence of PCR inhibitors or DNA degraded.
Moreover, the Ampliquality HPV-HS Bio Kit provides either one negative control or one positive control represented by plasmid DNA of HPV-54.

DNA isolation and multiplex real-time PCR
From January 2014, viral nucleic acids were extracted from the resuspended pellet using the automated MagNa Pure 96 system™ (Roche Diagnostics GmbH, Mannheim, Germany) according to the manufacturer's instructions.
DNA extraction was carried out starting from 200 µl of sample. DNA was eluated in a final volume of 100 µl. Extracted DNA samples were submitted to multiplex real-time PCR (mRT-PCR) by an Anyplex™ II HPV 28

Quality controls measures
Currently, in our laboratory, either intra-laboratory quality controls or inter-laboratory quality controls (VEQ) have been defined and activated. In the case of intra-laboratory quality controls, previously analysed samples are used to in order to evaluate the results obtained for a long time. With reference to the interlaboratory quality controls, we currently perform a commercially available VEQ program (DicoCARE VEQ).

Statistical analysis
The association of the DNA extraction methods, PCR assays, sample types, countries of women with either HPV prevalence or mHPV infection was evaluated by the χ 2 test, as appropriate.
To assess the association among each pairing of HPV types, multivariate logistic regression was used. In particular, predictor variables (viral types, PCR assay, extraction method and nationality) were decided a priori. Missing ages were imputed by multiple imputation by fully conditional specification implemented in the Multivariate Imputation by Chained Equations (MICE) package implemented in environment R. Predictive mean matching imputation model was specified on the assumption of missing at random ages and the number of iterations was set to 20. In particular, 50 imputed datasets were generated. For each dataset, a 2 Raffaele Del Prete et al.
logistic regression model was generated and the 50 models were pooled together by the function pool of the package mice. Globally, 26 logistic regression models were evaluated. To simplify the final logistic regression models, nested models with and without extraction method and nationality variables were compared by likelihood ratio tests and the collected P-values were corrected by Benjamini and Hochberg's (BH) procedure with false discovery rate (FDR) <1% [23].
Finally, all P-values collected from the final simplified logistic regression models were corrected for multiple comparison by BH procedure with FDR < 1% .
Calculations of all statistical tests were performed using the open source environment R [24].
The comparison of the 26 logistic regression models containing the variables age, nationality (Italian vs. non-Italian), PCR assay (mRT-PCR vs. nested PCR) and extraction method (Qiagen vs. MagnaPure 96) with the reduced models without the extraction method and nationality variables after BHs correction showed that all P-values were non-significant. For this reason the extraction method and nationality variables were not inserted in the 26 final models (data not shown).
The logistic regression analysis highlighted several statistically associated HPV combinations. Odds ratios are reported in Table 2 and the statistically significant associations after BHs correction have been highlighted in yellow colour. In particular, viral types 42 and 70 were involved in the major number of multiple positive associations. On the contrary, types 33, 43, 51, 58 and 82 were not involved in any statistically significant associations while 11, 16, 18, 40, 56 and 61 were involved in only one positive association. Negative associations were not detected by the regression analysis ( Table 2).

Discussion
Although HPV is considered the major causative agent for CC, little is known about the additional factors as the association among the pattern of HPV-types co-infections and the CC. To date, only few studies have reported a relationship between mHPV infections and CC, but the results obtained are not conclusive. Some studies have suggested a possible role of mHPVs  Raffaele Del Prete et al. in the development of CC [22,[25][26][27][28], while others works have reported that the risk of CC in women infected with mHPVs was no greater than that with single-type infections [29,30].
Our study investigates the prevalence of HPV-types in multiple infections. In many works, the differences observed when the prevalences of HPV infections had been studied, were due to several factors as geographic parameters, age, sexual behaviour, variables affecting immune response and HPV detection methods [29,31]. Similarly, the prevalence of mHPV infections also ranged significantly from 1% to 20%, due to the HPV detection methods used [29,31,32]. On basis of different extraction and amplification assays performed, we did not found significant differences in HPV prevalences (44.69% vs. 50.43% vs. 47.80%).
In agreement with Dickson et al. [33], our study suggests that mHPV infections frequently occurred (51.14%). In particular, the most common HR-HPV-type that we have detected was HPV-16. According to Oliveira et al. [34] and Chatuverdi et al. [22], HPV-16 is involved less frequently in single infections rather than in mHPV infections. Other HR-HPV-types identified more frequently in single infections are 53 (7.23%) and 31 (7.16%). Our results are in accordance with other epidemiological studies [35][36][37][38]. Moreover, only our study shows that 53 and HPV-31 (HR-HPV) together with 42 (LR-HPV) are the most representative HPV-types involved in co-infections. In addition, we have detected the 42 as the most prevalent HPV-type in either single or mHPV infections.
Moreover, Mori et al. [40] studying on the replication interference between HPV-16 and 18 have shown that these genotypes, when co-infected alone in absence of other HPV-types, might decrease the risk of CC, explaining their antagonistic action. However, in our study, we did not found any associations among HPV-16 and HPV-18, in agreement with Dickson et al. [33].
In addition, our study shows that all phylogenetic groups are involved in positively associated combinations, but only the HPV-types belonging to α-9, α-6 and α-7 species are the most prevalent in two-type combinations. On the contrary, these results are not supported by the studies of Trottier et al. [32] and Dickson et al. [33], which have described the α-9 species as the only involved in mHPV infections.
There are several limitations of this study. First, the reported data are not representative of the overall Apulia population so the results are not generalisable. A major limitation is the absence of the cytological data on the enrolled individuals. So, not having the data of histology of cervical lesions, it was not possible to correlate the presence of HPV with the cytological and/or histological findings and it was not possible to argue regarding the role of the mHPVs in the progression of the lesions. Moreover, the lack of behavioural information on the enrolled individuals did not permit to better characterise the analysed population for the roles of such exposure variables in the diffusion of HPV infections, therefore, we have not be able to analyse the association between these aetiological factors and HPV prevalence. Moreover, because of the presence of the same acquisition way, the effect estimate of measures of association between sexually transmitted infections may be confounded by common risk factors.
On the other hand, this study accounted for the analytical methods used to extract and genotype HPV. In fact, the two  Raffaele Del Prete et al. Table 3. Association between HPV phylogenetic groups by multivariate logistic regression models after P-values correction by BHs procedure Strength of associations is expressed by the odds ratios and 95% confidence interval. The multivariate logistic regression models are reported on the 10 columns, while the covariants on the 13 lines. Statistically significant associations are highlighted in yellow. Age: linear variable. PCR Assay: mRT-PCR vs. nested-PCR.
DNA extraction methods were not found associated with any HPV-types while the pairwise mHPVs detected were controlled for the two genotyping assays. To our knowledge, the evaluation of the analytical methods as confounding factor for mHPV associations has been object of very few studies [1].
In conclusion, the impact of HPV co-infections on the clinical outcome is still not clear. However, it is likely that only specific HPV combinations may be associated with a significant clinical impact while other combinations may simply be correlated because of the common infection route or to the diagnostic method used. In particular, some specific combinations could synergistically or antagonistically interact and subsequently affect the risk of development of CC. Therefore, further clinical studies are needed either to determine the mechanism of these mHPV infections or the clinical outcome related to particular combinations of co-infected HPV-types.