E6/E7 mRNA Expression of the Most Prevalent High-Risk HPV Genotypes in Cervical Samples from Serbian Women

Cervical cancer caused by persistent infection with HR HPV genotypes is the second leading cause of death in women aged 15 to 44 in Serbia. The expression of the E6 and E7 HPV oncogenes is considered as a promising biomarker in diagnosing high-grade squamous intraepithelial lesions (HSIL). This study aimed to evaluate HPV mRNA and DNA tests, compare the results according to the severity of the lesions, and assess the predictive potential for the diagnosis of HSIL. Cervical specimens were obtained at the Department of Gynecology, Community Health Centre Novi Sad, Serbia, and the Oncology Institute of Vojvodina, Serbia, during 2017–2021. The 365 samples were collected using the ThinPrep Pap test. The cytology slides were evaluated according to the Bethesda 2014 System. Using a real-time PCR test, HPV DNA was detected and genotyped, while the RT-PCR proved the presence of E6 and E7 mRNA. The most common genotypes in Serbian women are HPV 16, 31, 33, and 51. Oncogenic activity was demonstrated in 67% of HPV-positive women. A comparison of the HPV DNA and mRNA tests to assess the progression of cervical intraepithelial lesions indicated that higher specificity (89.1%) and positive predictive value (69.8–78.7%) were expressed by the E6/E7 mRNA test, while higher sensitivity was recorded when using the HPV DNA test (67.6–88%). The results determine the higher probability of detecting HPV infection by 7% provided by the mRNA test. The detected E6/E7 mRNA HR HPVs have a predictive potential in assessing the diagnosis of HSIL. The oncogenic activity of HPV 16 and age were the risk factors with the strongest predictive values for the development of HSIL.


Introduction
It is estimated that approximately every fourth malignancy can be linked to an infectious agent, that is, its contribution to various stages of cancer development (reviewed in [1]). About a third of this contribution is related to the human papillomavirus (HPV) (reviewed in [2]). Today, significant evidence confirms the association of high-risk (HR) HPV as a carcinogen or promoter in developing malignant diseases in different locations: the cervix, vulva, vagina, penis, anus, and certain head and neck regions. In first place are neoplasias of the lower genital tract, such as cervical cancer [3]. According to estimates by the World Health Organization (WHO), that is, by the International Agency for Research on Cancer (IARC), 604,000 new cases and 342,000 deaths were registered around the world in 2020, which makes cervical cancer the fourth most frequently diagnosed cancer in women [4]. In Serbia, organized cervical cancer screening has been conducted since Novi Sad, Serbia, and the Oncology Institute of Vojvodina, Serbia. All of the women in the study did not receive any prior treatment for cervical dysplasia or cancer, and all were unvaccinated against HPV infection. The samples were collected using the ThinPrep Pap test (Hologic Inc., Marlborough, MA, USA) according to the manufacturer's instructions and sent for further analyses to the Center of Virology, Institute of Public Health of Vojvodina, Novi Sad, Republic of Serbia.
The classification of cytological findings was performed according to the criteria of the Bethesda System 2014. It was categorized into negative for intraepithelial lesion or malignancy (NILM), atypical squamous cells of unknown significance (ASCUS), low-grade squamous intraepithelial lesions (LSIL), and high-grade squamous intraepithelial lesions (HSIL). All of the women enrolled in the study were informed about the research objective and signed an informed written consent form. The study protocol was reviewed and approved by the Medical Ethical Committee of the Institute of Public Health of Vojvodina, Novi Sad, Serbia (approval number: 01-252/3).

HR HPV Detection and Genotyping
The ThinPrep cervical smear samples were stored at 4-8 • C for up to 3 days from the sampling day. The 2 mL of collected samples were transferred to nuclease-free tubes and centrifugated at 8000× g for 5 min. The formed pellet was dissolved in 200 µL of nuclease-free water and used for nucleic acid extraction. According to the manufacturer's instructions, DNA extraction was carried out using the SaMag STD DNA Extraction Kit (Sacace Biotechnologies, Como, Italy). The extracted DNA was eluted in 100 µL elution buffer. The detection and genotyping of 12 HR HPV genotypes (16, 18, 31, 33, 35, 39, 45, 51, 52, 56, 58, and 59), marked as the HPV DNA test, were performed using the High Risk Typing Real-TM Kit (Sacace Biotechnologies, Como, Italy) following manufacturer's instructions. The E7 gene of specific HPV genotypes was amplified using primers and TaqMan probes in the multiplex reaction performed in a total of 13 µL. The β globin gene is used as an internal control. Real-time PCR was performed on the SaCycler-96 (Sacace Biotechnologies, Como, Italy). After the initial activation of the DNA polymerase at 95 • C for 15 min, five cycles of amplification were performed under the following conditions: 95 • C/5 s, 60 • C/20 s, and 72 • C/15 s, and 40 amplifications were performed under the following conditions: 95 • C/5 s, 60 • C/30 s (fluorescence detection), and 72 • C/15 s. The kinetics of the detected fluorescence signals were monitored using the SaCycler-96 software package (Sacace Biotechnologies, Como, Italy).

E6/E7 mRNA HPV Detection
E6/E7 mRNA of the most prevalent HPVs was tested in the cervical samples positive for the most prevalent HR HPVs DNA and HR HPV DNA negative samples. The HR-HPV-negative samples were included in E6/E7 mRNA testing because the study aimed to determine the mRNA test's clinical characteristics by evaluating and comparing it with the HPV DNA test. Total RNA was extracted from the prepared sample using the miRNeasy Mini Kit and QIAcube robotic workstation (Qiagen, Hilden, Germany) following the manufacturer's instructions. The total RNA was eluted in 50 µL ultrapure water free from nucleases. Following the manufacturer's recommendations, potentially present contaminants were removed using the TURBO DNA-free Kit (Invitrogen/ThermoFisher Scientific, Waltham, MA, USA). The routine procedure for removing contaminants using the kit above included the addition of 5 µL of 10× TURBO DNase Buffer and 1 µL of TURBO DNase enzyme into each sample of extracted total RNA, with incubation for 30 min at a temperature of 37 • C. After the action of the enzyme, 5 µL of inactivation reagent (Dnase Inactivation Reagent) was added, with incubation for 5 min, at room temperature and occasional vortexing. After that, centrifugation was performed (90 s, 10,000× g). The supernatant was carefully transferred to a nuclease-free tube. The real-time reverse transcription PCR (RT-PCR) analysis, marked as the HPV mRNA test, was performed using specific primers and TaqMan probes to detect the E6/E7 mRNA of individual HPV genotypes. The sequences for the primers and probes (Table 1) were adopted from Lindh et al. (2007) [19] and purchased from Life Technologies (Carlsbad, CA, USA). The AgPath-ID One-Step RT-PCR Kit (Applied Biosystems, Waltham, MA, USA) was used for the real-time RT-PCR. A separate reaction mixture was prepared for each set of primers and TaqMan probes. The reaction was prepared to a final volume of 25 µL containing: 12.5 µL 2× RT-PCR Buffer, 1 µL of 25× RT-PCR Enzyme Mix, primers to a final concentration of 300 nM, the probe to a final concentration of 200 nM, 1 µL of RNase Inhibitor reagent (Applied Biosystems, Waltham, MA, USA), 5 µL of isolated total RNA, and DEPC-treated nuclease-free water (Invitrogen, Waltham, MA, USA). Real-time PCR was performed on the Applied Biosystems 7500 Real-Time PCR Systems (ThermoFisher Scientific, Waltham, MA, USA). After the reverse transcription reaction at 48 • C for 30 min, the inactivation of reverse transcriptase and the activation of Taq polymerase were performed at 95 • C for 10 min. After that, 45 cycles of PCR amplification were carried out with denaturation at 95 • C for 15 s and annealing and elongation at 58 • C for 1 min. The data were analyzed with the Applied Biosystems Software v2.0.6 (ThermoFisher Scientific, Waltham, MA, USA) and the GraphPad Prism 8 (GraphPad Software, San Diego, CA, USA).

Statistical Analysis
All of the statistical analyses were performed using SPSS statistics software Version 21.0 (Chicago, IL, USA). Testing the difference in frequencies of attributive features was performed using the Chi-square (χ 2 ) test of independence and quality of the match. The Student's t-test was used to compare values between the two age groups, a numerical characteristic. A one-way analysis of variance (ANOVA) and the Bonferroni post-hoc test were applied to compare values between three or more data groups. Frequencies were used to present the analysis of the oncogenic activity of multiple HR HPV infections. Sensitivity, specificity, positive predictive value (PPV), negative predictive value (NPV), and their 95% confidence intervals (CIs) of HR HPVs DNA and E6/E7 mRNA HPVs detection and cytology test were calculated. To quantify the diagnostic capabilities of the selected test and evaluate its significance, the receiver operating characteristics (ROC) curve was used, which enables testing the significance of differences in the discriminating potential of different variables for the same binary outcome. It is based on a graphical presentation of pairs of sensitivity and specificity that can be obtained by estimating the threshold value for all values of discontinuous variables of the sample. Univariate and multivariate logistic regression were used to analyze the connection between two or more features, generating adequate statistical models. Multivariate logistic regression analysis was applied to all of the analyzed factors to construct a predictive model and named the most relevant predictors for the development of HSIL. A p-value of less than 0.05 defined as statistically significance.

Cervical Cytology
A total of 365 specimens obtained from women in the north part of the Republic of Serbia (Vojvodina) were classified based on the Bethesda System 2014 into four categories by cytological criteria.

HR HPV DNA in Cervical Samples
The cervical samples were analyzed for 12 HR HPVs, where 246 out of 365 (67.4%) had HPV-DNA-positive results, which indicates that the overall prevalence of HPV in the study population was 67.4%. All of the HPV genotypes covered by the HPV DNA test were identified (n = 274) in the study population (246 HR-HPV-positive cervical samples).  Figure 1). The results show that those HR HPVs make up 73% (200/274) of the detected genotypes, including multiple infections, which determined those cervical samples (n = 172) for further examination of oncogenic activity, according to the study's aim. Multiple HPV infections were found in 15.7%. The most common co-infections were those with HPV 16 and 31, found in 7.6% (13/172) of cases with multiple infections (Table 2).
Molecular data were compared with the cytological results and age categories of patients. The distribution of cytological groups was analyzed within the most prevalent HR-HPV-DNA-positive samples (172 cervical samples), including multiple infections. A minority of women, 16.9% (n = 29), had normal results, whereas 83.1% (n = 143) showed different abnormalities. A total of 26.7% (n = 46) of the examined women had ASCUS; in 25.6% (n = 44), LSILs were found, whereas HSILs were detected in 30.8% (n = 53). The mean age of the patients was 36.7 years. Among the specimens, the number of Serbian women who were ≤30 years, 31-44 years, and ≥45 years old accounted for 36.5% (n = 68), 36.0% (n = 62), and 24.4% (n = 42) of the samples, respectively (Table 3). used, which enables testing the significance of differences in the discriminating potential of different variables for the same binary outcome. It is based on a graphical presentation of pairs of sensitivity and specificity that can be obtained by estimating the threshold value for all values of discontinuous variables of the sample. Univariate and multivariate logistic regression were used to analyze the connection between two or more features, generating adequate statistical models. Multivariate logistic regression analysis was applied to all of the analyzed factors to construct a predictive model and named the most relevant predictors for the development of HSIL. A p-value of less than 0.05 defined as statistically significance.

Cervical Cytology
A total of 365 specimens obtained from women in the north part of the Republic of Serbia (Vojvodina) were classified based on the Bethesda System 2014 into four categories by cytological criteria.

HR HPV DNA in Cervical Samples
The cervical samples were analyzed for 12 HR HPVs, where 246 out of 365 (67.4%) had HPV-DNA-positive results, which indicates that the overall prevalence of HPV in the study population was 67.4%. All of the HPV genotypes covered by the HPV DNA test were identified (n = 274) in the study population (246 HR-HPV-positive cervical samples).   The distribution of the most frequently detected HPVs concerning cytology is shown in Table 4. The prevalence rates of HR HPV 16 ranged from 44.8% in the group of NILM cytology to 75.5% in the HSIL group. Contrarily, the prevalence of HR HPV 31 is similar in the groups of NILM (37.9%), ASCUS (34.8%), and LSIL (31.8%), while it is lower in the group of HSIL (11.3%). HPV genotypes 33 and 51 are present in all of the cytological groups in less than 21%. The statistically significant difference in the prevalence between the number of positive findings of HPV 16 (χ 2 test; p = 0.035) and HPV 31 (χ 2 test; p = 0.017) was determined, depending on the degree of severity of the cytological findings, which was not determined for HPV 33 and 51 (χ 2 test; p = 0.706, p = 0.790, respectively) ( Table 4).  A statistically significant difference was found in the number of female patients concerning the cytological findings and the age of the patients (χ 2 test; χ 2 = 29.500; p = 0.000) ( Table 5). The statistically significant difference was determined regarding the cytological findings and the age of the patients, where the female patients diagnosed with HSIL were significantly older compared to the other groups (ANOVA; F = 9.321; p < 0.001). The Bonferroni post-hoc test determined that the female patients with HSIL are statistically significantly older than those with ASCUS (p < 0.001), NILM (p < 0.001), and LSIL (p = 0.012) ( Figure 2, Table 5). Female patients with confirmed HPV 31 are statistically significantly younger (33 years) than the other HR-HPV-DNA-positive patients (t = 2.317; p = 0.022). The average age of the female patients with confirmed HPV DNA 16 was 36.9 years; with HPV DNA 33, it was 34.1 years, while the average age of patients with HPV DNA 51 was 40.5 years. The statistical analyses show that the proportion of HPV 16 positivity maintained at the same level as age. Conversely, the proportion of HPV 31 positivity decreases with age. The detection of HR HPV 33 genotypes decreases with age, while HR HPV 51 increases (Table 5).   Figure 3 shows the study design with sample processing to analyze the expression of the E6/E7 mRNA of the most prevalent HR HPV in cervical samples from Serbian women. A total of 291 cervical samples, which include HPV 16-, 31-, 33-and 51-positive (n = 172) and HR-HPV-negative samples (n = 119), were tested by the HPV mRNA test. The E6/E7 mRNA HPV was detected exclusively in HR-HPV-DNA-positive samples ( Table 6). E6 and E7 transcripts of the four most frequent HR HPVs were detected in 57.5% (115/200) of the HR HPV DNA confirmed genotypes. Accordingly, the distribution of E6/E7 mRNA HR HPV 16, 31, 33, and 51 are shown in Table 6. The E6/E7 mRNA HR HPV 16 was the most abundant, which accounted for 25.5% (51/200) of HR HPV genotypes. Next in frequency was the E6/E7 mRNA HR HPV 31 in 16.5% (33/200), while the E6/E7 mRNAs HR HPV 33 and 51 were equally represented in 8% (16/200) and 7.5% (15/200), respectively. Almost every second HPV 16 genotype is oncogenically expressed (48.6%; 51/105), and it was detected in 29.7% (51/172) HPV-DNA-positive samples. The oncogenic activity of HPV 31 was detected in approximately every fifth (19.2%; 33/172) HPV-DNA-positive sample. Regarding the oncogenic activity of the remaining tested genotypes, HPV 33 and HPV 51 were detected in roughly every tenth HPV-DNA-positive sample ( Table 6). The results of expressing E6 and E7 HR HPV oncogenes were expressed through the dispersion of the obtained Ct values. The oncogenic activity of HPV 16 is detected by the lowest registered value (Ct = 16) (Supplementary Materials Figure S1).  Figure 3 shows the study design with sample processing to analyze the expression of the E6/E7 mRNA of the most prevalent HR HPV in cervical samples from Serbian women. A total of 291 cervical samples, which include HPV 16-, 31-, 33-and 51-positive (n = 172) and HR-HPV-negative samples (n = 119), were tested by the HPV mRNA test. The E6/E7 mRNA HPV was detected exclusively in HR-HPV-DNA-positive samples ( Table 6). E6 and E7 transcripts of the four most frequent HR HPVs were detected in 57.5% (115/200) of the HR HPV DNA confirmed genotypes. Accordingly, the distribution of E6/E7 mRNA HR HPV 16, 31, 33, and 51 are shown in Table 6. The E6/E7 mRNA HR HPV 16 was the most abundant, which accounted for 25.5% (51/200) of HR HPV genotypes. Next in frequency was the E6/E7 mRNA HR HPV 31 in 16.5% (33/200), while the E6/E7 mRNAs HR HPV 33 and 51 were equally represented in 8% (16/200) and 7.5% (15/200), respectively. Almost every second HPV 16 genotype is oncogenically expressed (48.6%; 51/105), and it was detected in 29.7% (51/172) HPV-DNA-positive samples. The oncogenic activity of HPV 31 was detected in approximately every fifth (19.2%; 33/172) HPV-DNA-positive sample. Regarding the oncogenic activity of the remaining tested genotypes, HPV 33 and HPV 51 were detected in roughly every tenth HPV-DNA-positive sample ( Table 6). The results of expressing E6 and E7 HR HPV oncogenes were expressed through the dispersion of the obtained Ct values. The oncogenic activity of HPV 16 is detected by the lowest registered value (Ct = 16) (Supplementary Materials Figure S1).

Prevalence of HR HPV Based on E6/E7 mRNA HPV Expression in Different Cytological Groups
The expression of the E6 and E7 genes as indicators of the oncogenic activity of HR HPV 16, 31, 33, and 51 was analyzed concerning cytological findings. The oncogenic activity of the tested genotypes increases with the severity of the cervical intraepithelial lesion. A statistically significant difference in E6/E7 mRNA HPV expression among the various cytological groups was observed (χ 2 test; χ 2 = 108.623; p < 0.001). E6/E7 mRNA HPVs are the most prevalent in patients with HSIL cytological findings (88.9%). In the group of patients with LSIL cytological findings, it was demonstrated in a lower percentage (60%). A two-fold lower prevalence is observed in patients with ASCUS (29.4%). Oncogene activity in women with normal cytological findings is present in 10.9% of samples (Table 7). Subsequently, the E6/E7 mRNA distribution of HR-HPV-DNA-positive samples according to genotype and cytological groups was analyzed. E6/E7 mRNA HR HPV 16 is the most represented in patients with HSIL cytological findings (64.2%), while in the other groups of cytological findings, it was demonstrated in a lower percentage (3.4-22.7%). A statistically significant difference was found in the number of positive findings of E6/E7 mRNA HR HPV 16 concerning the cytological status (χ 2 test; χ 2 = 46.881; p < 0.001). This result singled out the HR HPV 16 genotype for further analyses. The presence of E6/E7 mRNA HR HPV 31 was the least detected in HSIL (7.5%) compared to other cytological groups (22.7-26.1%). The distribution of the oncogenic activity of the remaining genotypes (HR HPV 33 and 51) is approximately the same across different cytological groups and remains at a low level (2.2-11.4%) ( Table 7).   The analyses of the oncogenic activity of multiple HR HPV infections are presented by frequencies. The overall oncogenic activity, including both single and multiple HR HPV infections detected using the E6/E7 mRNA HR HPV test, increases with the degree of cervical lesion severity (60-100%). The oncogenic activity detected in single-genotype infections (20.0-85.7%) is higher compared to the oncogenic activity of multiple genotypes (0-40%) ( Table 8).

Prevalence of E6/E7 mRNA HR HPV Expression According to Age
The results were categorized according to age categories (≤30, 31-44, ≥45 years) to analyze the prevalence of E6/E7 mRNA HPV in the specific age groups. The lowest percent of E6/E7 mRNA HR HPV 16 was detected in the younger group and the highest percent was detected in patients over 44 years. A statistically significant difference was observed in patients with E6/E7 mRNA HR HPV 16 expression (χ 2 test; χ 2 = 7.331; p = 0.026), wherein individuals with positive E6/E7 mRNA HPV results were older than the others. The detection of E6/E7 mRNA HR HPV 31 and 33 decreases with the increasing age of the patient, while E6/E7 mRNA HR HPV 51 increases (Table 9).

Comparison of Tests for the Detection of HR HPV Genotypes and Their Oncogenic Activity
The comparison of tests for the detection of the most prevalent HR HPV genotypes and their oncogenic activity revealed that the presence of the HR HPV genotypes is higher than the presence of the oncogenic activity of the genotypes in younger women (≤30 years), similar in middle-aged women (31-44 years), and lower in women 45 years and older. The prevalence of oncogenic activity of the HPV genotypes increases with the severity of the cervical intraepithelial lesion. Compared to the prevalence of the examined HR HPV genotypes, the same parameter is lower in women with normal and undefined cytological findings, approximately the same in low-grade lesions, and significantly higher in high-grade lesions ( Figure 4). The comparison of tests for the detection of the most prevalent HR HPV genotypes and their oncogenic activity revealed that the presence of the HR HPV genotypes is higher than the presence of the oncogenic activity of the genotypes in younger women (≤30 years), similar in middle-aged women (31-44 years), and lower in women 45 years and older. The prevalence of oncogenic activity of the HPV genotypes increases with the severity of the cervical intraepithelial lesion. Compared to the prevalence of the examined HR HPV genotypes, the same parameter is lower in women with normal and undefined cytological findings, approximately the same in low-grade lesions, and significantly higher in high-grade lesions ( Figure 4). The calculated clinical characteristics of the DNA and E6/E7 mRNA HR HPV tests are shown in Table 10. The sensitivity and NPV of both tests were increased with the severity of the cervical intraepithelial lesion, with the HR HPV DNA test showing a statistically significantly higher level (67.6-98.2%). The specificity of the E6/E7 mRNA HR HPV test (89.1%) The calculated clinical characteristics of the DNA and E6/E7 mRNA HR HPV tests are shown in Table 10. The sensitivity and NPV of both tests were increased with the severity of the cervical intraepithelial lesion, with the HR HPV DNA test showing a statistically significantly higher level (67.6-98.2%). The specificity of the E6/E7 mRNA HR HPV test (89.1%) is statistically significantly higher than the HR HPV DNA test (75.6%). The PPV of the HR HPV DNA test is approximately the same for all types of cytological findings (60.3-64.6%), while for the mRNA HR HPV test, it increases with the degree of the cervical intraepithelial lesion (69.8-78.7%) and it is statistically significantly higher (Table 10). To quantify the diagnostic capabilities of the E6/E7 mRNA HPV test and evaluate its significance, an ROC curve was used to assess the assays for detecting HSIL. The area determined by the ROC curve (AUC) of E6/E7 mRNA HR HPV is 0.812 (CI (95%): 0.752-0.871), while the area under the ROC curve formed by the parameters of the HR HPV DNA test is 0.740 (CI (95%): 0.680-0.799) (Table 11, Figure 5).  Figure 5).  The relationships between the oncogenic activity of all of the tested HPVs and previously singled out HPV 16 vs. the cytological results were analyzed using Spearman's (ρ) correlation ( Figure 6). There is a statistically significant positive moderate correlation between the presence of HPV 16 oncogenic activity and cytological status (Spearman's cor- The relationships between the oncogenic activity of all of the tested HPVs and previously singled out HPV 16 vs. the cytological results were analyzed using Spearman's (ρ) correlation ( Figure 6). There is a statistically significant positive moderate correlation between the presence of HPV 16 oncogenic activity and cytological status (Spearman's correlation; ρ = 0.494; p < 0.001) ( Figure 6A). The oncogenic activity of the tested HPV genotypes is associated with a strong statistically significant positive correlation with the degree of cervical intraepithelial lesion severity (Spearman's correlation; ρ = 0.594; p < 0.001) ( Figure 6B). Univariate multinomial logistic regression examined the individual factors (the HR HPV DNA 16, the oncogenic activity of all of the tested HR HPVs, the oncogenic activity of the HR HPV 16 genotype, and the age category) that indicated an increased probability of developing HSIL (Supplementary Materials Tables S1-S4).
For the influence of HR HPV 16 on the probability of HSIL, a statistically significant predictive value was determined in the NILM and LSIL cytological groups. Patients with the confirmed presence of the HPV 16 genotype through a DNA test have a higher probability of being diagnosed with HSIL than NILM (3.8-fold), while they will have a 2.6-fold higher probability of being diagnosed with HSIL compared to the probability of being diagnosed with LSIL (Supplementary Materials Table S1).
The influence of the oncogenic activity of all tested of the HR HPV genotypes (E6/E7 mRNA HR HPVs) on diagnosing high-grade lesions of the cervical epithelium has a statistically significant prediction found in all cytological groups, NILM, ASCUS, and LSIL. If patients have confirmed indicators of oncogenic activity, they will have an almost sevenfold higher probability of being diagnosed with HSIL compared to the probability of being diagnosed with NILM. The same patients have a 19-fold higher probability of being diagnosed with HSIL compared to the probability of detecting ASCUS and a 5-fold higher probability of detecting HSIL compared to LSIL (Supplementary Materials Table S2).
E6/E7 mRNA HR HPV 16 is an indicator for diagnosing HSIL, and a statistically significant predictive value was determined for all types of cytological findings, NILM, AS-CUS, and LSIL. In a patient with confirmed HPV 16 oncogenic transcripts, the probability of diagnosing HSIL is 50-fold higher than the probability of diagnosing NILM. Patients with a positive result of E6/E7 mRNAs HR HPV 16 are 12-fold more likely to be diagnosed with HSIL compared to the detection of ASCUS, while in patients with the same result, the probability of detection of HSIL is 6-fold higher than the probability of detection of LSIL (Supplementary Materials Table S3).
A statistically significant predictive value for the age category was determined in all cytological groups, NILM, ASCUS, and LSIL. The probability of detecting HSIL concerning normal results in HR-HPV-positive patients is 6.3-fold higher if they are ≥ 45 years of age than women under 30. A similar prediction for the detection of HSIL was shown concerning ASCUS (6.5-fold). Patients from the oldest age category have a 3.4-fold higher probability of detecting HSIL concerning LSIL than the youngest (Supplementary Mate- Univariate multinomial logistic regression examined the individual factors (the HR HPV DNA 16, the oncogenic activity of all of the tested HR HPVs, the oncogenic activity of the HR HPV 16 genotype, and the age category) that indicated an increased probability of developing HSIL (Supplementary Materials Tables S1-S4).
For the influence of HR HPV 16 on the probability of HSIL, a statistically significant predictive value was determined in the NILM and LSIL cytological groups. Patients with the confirmed presence of the HPV 16 genotype through a DNA test have a higher probability of being diagnosed with HSIL than NILM (3.8-fold), while they will have a 2.6-fold higher probability of being diagnosed with HSIL compared to the probability of being diagnosed with LSIL (Supplementary Materials Table S1).
The influence of the oncogenic activity of all tested of the HR HPV genotypes (E6/E7 mRNA HR HPVs) on diagnosing high-grade lesions of the cervical epithelium has a statistically significant prediction found in all cytological groups, NILM, ASCUS, and LSIL. If patients have confirmed indicators of oncogenic activity, they will have an almost seven-fold higher probability of being diagnosed with HSIL compared to the probability of being diagnosed with NILM. The same patients have a 19-fold higher probability of being diagnosed with HSIL compared to the probability of detecting ASCUS and a 5-fold higher probability of detecting HSIL compared to LSIL (Supplementary Materials Table S2).
E6/E7 mRNA HR HPV 16 is an indicator for diagnosing HSIL, and a statistically significant predictive value was determined for all types of cytological findings, NILM, ASCUS, and LSIL. In a patient with confirmed HPV 16 oncogenic transcripts, the probability of diagnosing HSIL is 50-fold higher than the probability of diagnosing NILM. Patients with a positive result of E6/E7 mRNAs HR HPV 16 are 12-fold more likely to be diagnosed with HSIL compared to the detection of ASCUS, while in patients with the same result, the probability of detection of HSIL is 6-fold higher than the probability of detection of LSIL (Supplementary Materials Table S3).
A statistically significant predictive value for the age category was determined in all cytological groups, NILM, ASCUS, and LSIL. The probability of detecting HSIL concerning normal results in HR-HPV-positive patients is 6.3-fold higher if they are ≥45 years of age than women under 30. A similar prediction for the detection of HSIL was shown concerning ASCUS (6.5-fold). Patients from the oldest age category have a 3.4-fold higher probability of detecting HSIL concerning LSIL than the youngest (Supplementary Materials Table S4).
The predictive model contains four independent variables: HR HPV DNA 16, E6/E7 mRNA HR HPVs, E6/E7 mRNA HR HPV 16, and age category. The HSIL lesion, as a representative of a high degree of cervical atypia, represented a dependent variable concerning all of the analyzed relevant factors. Multivariate logistic regression analysis was applied to all of the analyzed factors to construct a predictive model and named the most relevant predictors for the development of HSIL (Table 12). Analyzing the mutual influence of the examined relevant factors for diagnosing HSIL compared to the probability of diagnosing a normal cytological finding, the strongest statistically significant predictor was determined to be the oncogenic activity of the HPV 16 genotype. If it is detected, the probability of diagnosing HSIL concerning the probability of diagnosing normal results increases 19-fold (OR = 19.10; CI (95%): 1.54-236.98; p = 0.022). A statistically significant but almost three-fold weaker predictor for the detection of HSIL compared to NILM is the patient belonging to the oldest age category (OR = 6.65; CI (95%): 1.66-26.60; p = 0.007), while female patients from the age category 31-44 years have a slightly lower statistically significant probability (OR = 5.38; CI (95%): 1.36-21.30; p = 0.016) for the detection of the same lesion. The age category was determined as the strongest statistically significant predictor by analyzing the mutual influence of relevant factors for diagnosing HSIL concerning the probability of detecting ASCUS. HR HPV DNA-positive patients belonging to the oldest age category (≥45 years) have an 8.7-fold higher probability of being diagnosed with HSIL compared to the probability of being diagnosed with ASCUS compared to patients belonging to the younger age category (OR = 8.74; CI (95%): 2.15-35.57; p = 0.002). Female patients with the proven oncogenic activity of HPV 16 have a probability of detecting HSIL 6.4-fold higher compared to the probability of being diagnosed with ASCUS (OR = 6.38; CI (95%): 1.22-33.54; p = 0.029). By analyzing the mutual influence of relevant factors for diagnosing HSIL concerning the probability of detection of LSIL, the strongest statistically significant predictor was determined to be the oncogenic activity of HPV 16 (OR = 5.10; CI (95%): 1.09-23.83; p = 0.038), while weaker statistically significant predictability is shown by the patient's belonging to a specific age category. Female patients belonging to the oldest age category (≥45 years) have a 3.7-fold greater probability of detecting an HSIL change compared to the detection of LSIL compared to the younger age category (≤30 years) (OR = 3.72; CI (95 %): 1.16-11.92; p = 0.027) ( Table 12).

Discussion
Persistent infection caused by HR HPV is the leading risk factor for developing cervical intraepithelial lesions and cervical cancer (reviewed in [20]). Many countries have introduced screening programs based on the cytomorphological examination of cervical samples using the PAP test during the last 60 years to reduce the morbidity and mortality caused by cervical cancer. However, this procedure has shown less than optimal sensitivity (50%) and high inter-and intra-individual variability [21,22]. HPV DNA detection and genotyping provide an efficient screening method and enable risk stratification. However, just proving the presence of the HPV genome in a cervical epithelium does not provide insight into the type of infection. It does not answer whether there is a transient or persistent infection in which the virus is actively multiplying and in which there is a high risk of cancer developing (reviewed in [23]). It was established that cervical carcinogenesis is strongly associated with the HPV-caused infection in which the transcription of E6 and E7 HR HPV oncogenes occurs, with the consequent increase of their mRNA and protein levels. For this reason, the detection of E6 and E7 mRNA of HR HPVs can serve as a promising biomarker of their persistence and oncogenic activity (reviewed in [23]), which could enable a better assessment of the progression to high-grade cervical intraepithelial lesions and, in this regard, significantly influence the algorithm for monitoring patients (reviewed in [18,23,24]).
In our study, 365 female patient samples from Serbia were tested for the twelve HR HPV genotypes. It should be emphasized that the study has limitations. Used HPV tests are non-validated according to European Guidelines for use in primary screening (Meijer HPV test criteria) or according to the VALidation of HPV GENoyping Tests (VALGENT) protocols [25]. The commercial molecular HPV test was chosen for research according to Serbia's general requirements for molecular diagnostics and the limited research expenses.
From the total number of tested samples, 246 (67.4%) samples were positive for at least one of the 12 tested HR genotypes. These results agree with previous studies in the same area and European countries, where a high prevalence of HPV was registered. Studies in Serbia have reported that the presence of HR HPV infections range from 50% to 79% of women [26,27]. In the countries of our region, such as Croatia, a similar prevalence was shown (59%) [28] and in Bulgaria (61%) [29], while a slightly lower prevalence was registered in Italy (53%) [30]. Contrarily, a low prevalence of HPV infection was registered in Western and North European countries. Some of them are Great Britain (20.6%), Sweden (9.7%), and the Netherlands (3.8%) [7].
HPV genotype 16 (38%) emerged as the most frequently detected genotype in the study group, in concordance with previous reports [26,27]. The presence of HPV 16 is more often present (76%) in high-grade lesions compared to lower-grade lesions of the cervical epithelium. Studies on the HPV 16 genotype's prevalence concerning cytological findings confirm these results (reviewed in [6]). The following frequencies also represented the Alpha-9 genotype, HPV 31 (17%) and HPV 33 (9%). A meta-analysis study on five continents shows that HPV 31 is especially frequent in Europe [31]. The results of this research indicate that the frequency of the HPV 31 genotype statistically significantly decreases with the progression of the intraepithelial lesion. The decreasing trend of the same genotype's presence according to the degree of severity of the cytological lesion, more precisely, a lower prevalence in cancers compared to precancerous lesions, is observed in research from other studies [32][33][34]. Like HPV 31, a prevalence decrease of HPV 33 in HSIL compared to normal cytology was also registered in this research. Globally, HPV 33 ranks fourth in frequency and is responsible for 4.2% of all registered cervical cancers [6]. The surprising fact is a markedly high prevalence (9%) of the Alpha-5 genotype HPV 51 in our research. In the context of the causative agents of cervical cancer, HPV 51 is not included in the top ten most frequent HPV genotypes registered worldwide [6]. However, the detection of this genotype within this research, as well as previous studies from our region [26] and certain European countries [35][36][37][38][39][40][41], places it among the first four most prevalent genotypes detected in precancerous lesions of the cervix. Since HPV vaccines do not protect against all oncogenic HPVs, such as HPV 51, a complete understanding of its oncogenic activity is particularly significant [42]. The remaining tested genotypes (HR HPV 52,56,45,18,59,58,39,and 35) were present in a low percentage which is in concordance with previous reports [26]. Although HPV 18 is considered to be responsible for 15% of invasive cervical cancers, it is essential to note that its prevalence is similar to some studies from neighboring countries [29,43,44]; we found that in our area, HPV 18 was present at a lower percentage. Its frequency (3.6%) was 10-fold lower than that of HPV 16 ( Figure 1). The HPV vaccine was introduced in over half of the WHO member countries in 2020 [5]. There is scientific data from numerous countries that have implemented HPV vaccines in their routine immunization programs on decreasing the burden of cervical HPV infections and precancers [45]. According to our data, using the nine-valent vaccine could prevent more than 80% of the cervical precancerous lesions identified in this study.
The presence of HR HPVs determined further examination of their oncogenic potential. In our study, the expression of E6/E7 mRNA HR HPV was identified in 67% of the HPVpositive samples ( Table 6). The percentage of expression of the E6 and E7 HPV-examined HR HPVs is proportional to the degree of severity of the cervical lesion (Table 7). It can be observed that approximately every tenth HPV-infected woman (11%) with normal cytological findings is infected with an oncogenically expressed HPV. At the same time, in HSIL, this relationship is the opposite. Namely, the absence of indicators of HPV oncogenic activity is detected in approximately every tenth woman with HSIL status (11.1%). An undoubted trend in E6/E7 mRNA HR HPV positivity with increasing cytology severity has been observed in the data of studies conducted in different regions of the world [46][47][48][49][50]. In agreement with the report of Argyri et al. (2013) and according to the mRNA test, E6/E7 expression was prevalent in 9.1% of women with normal cytology, similar to our study [51]. However, other previously published data indicated the prevalence of those transcripts in a lower proportion of women with normal cytology findings than our study (0%) [52].
In our study, HPV 16 constituted 29.7%, followed by genotype 31 (19.2%), genotype 33 (9.3%), and genotype 51 (8.7%). The results showed that approximately every second HPV 16 genotype is oncogenically expressed (48.6%). Our data agree with the results reported in other studies. Rossi et al. (2017) observed a positivity rate for E6/E7 mRNA HR HPV, ranging from 58% to 77% in HPV-DNA-positive women [50]. A study by Tüney et al. (2017) in Turkey registered 55.6% E6/E7 mRNA HR HPVs, where HPV genotype 16 constituted 57.8% [24]; similar to that, Bruno et al. (2018) found that the HPV 16 genotype was the most oncogenetically active [46]. As confirmed in this research, it is also stated by several other authors that the transcription product HPV 16 is most often detected, which indicates that the tendency of expression of the essential oncogenes of this genotype is significantly higher compared to the other examined genotypes [24,53,54], which gives it the status of the HPV with the most carcinogenic potential [55]. The oncogenic activity of HPV 31 was detected in approximately every fifth (19.2%) HPV-DNA-positive sample. Contrary to HPV 16, the opposite trend is observed with the degree of cervical lesion severity. It can be assumed that the negative transcriptional status of E6 and E7 oncogenes is more prevalent due to the presence of an episomal form of the virus or an established transcriptional control that enables the spontaneous elimination of infection. Within this research, the results of E6/E7 mRNA HPV 51 detection are approximately the same as those of HPV 33. The distribution of oncogenic activity of these two genotypes is approximately the same across cytological groups and remains at a low level (2.2-11.4%) ( Table 7). Other studies have registered different percentage ranges of transcriptional detection of oncogenes of a particular genotype (HPV 33), from its absence [52] to complete expression of 100% [55]. The reason for those differences can be explained by the variations in their number in the total sample [55]. Furthermore, the total oncogenic activity increases with the degree of cervical lesion severity (60-100%). The oncogenic activity of individual genotypes (20-86%) is higher than that of multiple genotypes (0-40%) in the same type of lesion and increases with the degree of its severity (Table 8). These data are supported by the literature. Each of the detected genotypes is considered to have an independent mechanism of action in oncogenesis [56][57][58][59], which supports the hypothesis that different cells being infected with different viral genotypes rather than their intracellular coexistence is possible [60].
To analyze the oncogenic activity of the four most frequently diagnosed HR HPVs related to the age categories, a statistically significantly higher prevalence of positive E6/E7 mRNA HPV 16 findings was observed in the older age groups (Table 9). In agreement with this result, studies have shown that E6 and E7 mRNA HR HPV detection is significantly higher after 30 years [48,61]. The same result is supported by population-based cohort studies, which state that the majority of young women have an asymptomatic HPV infection which, thanks to the immune response, acquires transitory status [49,62].
The growing interest in molecular diagnostics methods has led various authors to compare the characteristics of the HPV DNA test with the mRNA test, evaluating the diagnostic accuracy in identifying high-grade cervical lesions. Our data showed that the specificity (89%) and PPV (70-79%) of the mRNA test are statistically significantly higher, while the same test has a statistically significantly lower sensitivity than the HPV DNA test. These results agree with the statements of different authors, who emphasized slightly lower values of the sensitivity of the mRNA test compared to the sensitivity of the HPV DNA test. At the same time, the specificity is expressed at a higher level [48,49] (reviewed in [63]). In contrast, others suggest that the sensitivity is similar and that its application would significantly improve the assay's specificity characteristics [49]. Studies evaluating the reliability of the mRNA assay showed heterogeneous findings (reviewed in [32]). The meta-analysis results indicated that the obtained sensitivity value ranged from 41% to 95%, while the registered specificity was 42-97% (reviewed in [64]). The observed difference is explained by the heterogeneous participation of cervical pathology [64] and methodological quality in different studies, which highlights the limitations in the general interpretation of these test characteristics [32]. Although the results are presented broadly, they maintain a unique trend and suggest that the mRNA test over the HPV DNA test improves specificity [64]. By comparing the clinical characteristics of the test, it can be concluded that the detection of E6 and E7 mRNA HR HPV compared to HPV DNA represents a much better marker for more accurate screening of high-grade cellular atypia of the cervix (reviewed in [32]), which make it an appropriate tool for the secondary screening of cervical cancer [47,64].
In our study, the results of the ROC curve analysis indicated that the probability of detecting HPV infection with the mRNA test (81%, AUC = 0.812) in patients with highgrade lesions is higher than the possibility of diagnosing them with the HPV DNA test (74%, AUC = 0.740). The obtained results are in line with the other studies emphasizing the potential usefulness of this test. Sun et al. (2021) also found higher AUC values for the mRNA assay compared to the DNA HPV (0.929 vs. 0.833) [65], while according to Yao et al. (2017), the clinically relevant portion of the AUC of mRNA was 0.721 [66].
According to the previously established degree of influence on diagnosing HSIL, relevant factors were selected in our research. Factors that represent an increased risk for more severe cervical changes were gradually examined (the HPV 16 genotype, the total oncogenic activity of all of the tested HPV genotypes, the HPV 16 oncogenic activity, and the age category of the patient).
Firstly, the HPV-DNA-16-positive results have a statistically significant predictive value for diagnosing HSIL (Table S1). In the same context, the study of Bruno et al. (2018) stated that the presence of the HPV 16 genotype was associated with a five-fold higher risk of developing a high-grade lesion compared to women with the presence of another HPV genotype [46]. Similarly, the data from a recent study indicate that type-specific HPV persistence predicted high-grade lesions, with HPV 16 being the most common type [67].
Secondly, the oncogenic activity of all of the tested HR HPVs has a statistically significant predictive value for diagnosing HSIL (Table S2). The obtained results support the statements related to the research on the importance of oncogenic activity, which show that the detection of E6 and E7 mRNA HR HPV could have a prognostic value in monitoring the development of carcinogenesis [68,69]. The results of the study by Fontecha et al. (2016) [55] showed that in the highest percentage of E6 and E7 mRNA HPV-positive women, the progression of the lesion is diagnosed over time (53%), followed by the persistence of abnormal cytological findings (42%), while regression was recorded in 10-fold lower cases (4%).
Furthermore, in patients with positive results of indicators of the oncogenic activity of HPV 16, a statistically significantly higher probability of diagnosing a high-grade lesion was determined in all types of cytological groups (Table S3). According to the previous insights into the published statements, one of the few prospective follow-up studies that dealt with the predictive value of E6 and E7 mRNA HPV 16 indicated that through the detection of this biomarker, it is possible to identify 87.5% of the HPV infections that progressed. In this case, the risk of progression of negative cytology and low-grade cervical lesions was 10-fold higher than that detected in mRNA-negative women during a follow-up period of 35 months [70]. This is supported by the results of Johansson et al. (2015), which indicated that the absence of E6/E7 mRNA HPV demonstrated a high negative predictive value for the future development of high-grade lesions of the cervix among HR-HPV-DNA-positive women with ASCUS and LSIL [68].
Within the results of this research, a statistically significant final model with all of the predictors was constructed, which shows the strength of the potential of the examined factors, that is, biomarkers for diagnosing precancerous lesions in women with HR HPV infection. Looking at each cytological group individually, two independent variables made a statistically significant contribution to the model and were thus named as the strongest predictors. These are the oncogenic activities of HR HPV 16 and the age category (≥45 years). It is known that the presence of HPV infection is confirmed in all age categories. However, belonging to a particular age group is a determinant significantly associated with the risk of acquiring this infection, depicting the peak in prevalence, which generally takes place around 20-25 years of age. For lesions to progress to more severe forms of cervical disease, a period of 5-14 years is necessary, during which the infection persists and the process of oncogenesis takes place [71]; the results of this research confirm the stated findings. In the examined women, the age category (≥45 years) is a statistically significant prognostic factor for diagnosing HSIL in all of the cytological groups (Tables 12 and S4). Loopik et al. (2020) [58] indicated that the risk of progression of an existing high-grade lesion increases with age, i.e., in women over 50, the risk of developing cervical cancer increases by seven fold. In addition, the risk of developing cervical abnormalities and the need to use an mRNA test in the diagnostic protocol of HPV-DNA-positive postmenopausal women with normal cytology is emphasized by Asciutto et al. (2020) [72].

Conclusions
In summary, this study describes the detection rates of the most common HR HPVs (16, 31, 33, and 51) and E6/E7 mRNA HR HPV expression in 365 Serbian women who showed normal and abnormal cytological findings. Those HR HPV genotypes are oncogenically active in more than half of the examined cases, and the detected oncogenic activity has predictive potential in diagnosing high-grade cervical intraepithelial lesions. According to the constructed predictive model, the oncogenic activity of HPV 16 and age are risk factors with the strongest predictive values for developing those lesions. Thus, our data indicate that mRNA testing may be more relevant than HPV DNA for assessing lesion grade and diagnosing and monitoring women at risk of progressive cervical disease. This way, the mRNA test as a tool for better risk stratification of HPV infection could overcome unnecessary examinations, increased costs, and patient anxiety. However, further follow-up studies are needed to determine the clinical utility of the mRNA HR HPV test.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/diagnostics13050917/s1. Figure S1: Dispersion of Ct values of E6/E7 mRNA HR HPVs. Table S1: Analysis of HR HPV DNA 16's influence on the diagnosis of HSIL; Table S2: Analysis of total E6/E7 mRNA HR HPV's influence on the diagnosis of HSIL; Table S3: Analysis of E6/E7 mRNA HR HPV 16's influence on the diagnosis of HSIL; Table S4: Analysis of the influence of age on the diagnosis of HSIL.  Informed Consent Statement: All of the women provided written consent for the use of the biological specimens for research purposes.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.