Philippe Pérot

Anne Biton

Jacques Marchetta

Anne-Gaelle Pourcelot

André Nazac

Henri Marret

Thomas Hebert

Delphine Chrétien

Marie-Christine Demazoin

Michael Falguières

Thomas Hébert

Michaël Falguières

Laurence Arowas

Hélène Laude

Isabelle Heard

Marc Eloit

Broad-Range Papillomavirus Transcriptome as a Biomarker of Papillomavirus-Associated Cervical High-Grade Cytology

INTRODUCTION

Human papillomaviruses (HPV) infections are associated with the development of cervical carcinoma, one of the most common cancers among women, and other cancers like anal cancer [START_REF] Lin | Human papillomavirus types from infection to cancer in the anus, according to sex and HIV status: a systematic review and meta-analysis[END_REF] and head and neck cancer [START_REF] Chaturvedi | Human papillomavirus and rising oropharyngeal cancer incidence in the United States[END_REF] . HPV are the etiologic agents responsible for over 99% of all cervical cancers [START_REF] Walboomers | Human papillomavirus is a necessary cause of invasive cervical cancer worldwide[END_REF] . HPV are small, non-enveloped DNA viruses commonly transmitted through sexual contact, which infect basal cells and replicate in the nucleus of squamous epithelial cells. HPV include more than 200 genotypes characterized by their oncogenic potential, with highly oncogenic HPV types (high-risk HPV) having a unique ability to drive cell proliferation [START_REF] Schiffman | Carcinogenic human papillomavirus infection[END_REF] .

The genomic organization of papillomaviruses is divided into functional early and late regions. The model of HPV infection, which is mainly derived from knowledge on HPV16, is that following the infection of basal cells in the cervical epithelium, the early HPV genes (E6, E7, E1, E2, E4 and E5) are expressed and the viral DNA replicates from the episomal form of the viral DNA. As the cells divide, in the upper layers of the epithelium the viral genome is replicated further, and the late genes (L1 and L2) and E4 are expressed. Viral shedding then further initiates new infections [START_REF] Woodman | The natural history of cervical HPV infection: unresolved issues[END_REF] .

HPV infection during the development of cervical cancer is associated with a shift from productive infection (which in most of the cases will be cleared by the immune system), towards non-productive persistent and transforming infection (in a minority of cases) characterized in particular by a high level of E6 and E7 mRNAs and low expression of E2 and late genes such as L1 [START_REF] Doorbar | The biology and lifecycle of human papillomaviruses[END_REF][START_REF] Shulzhenko | Ménage à trois: an evolutionary interplay between human papillomavirus, a tumor, and a woman[END_REF] . High-risk HPV infection may result in low-grade lesion, with highly productive infection and high rate of spontaneous regression. In contrast, high-risk persistent HPV infection is responsible for high-grade lesion, the true precancerous lesion.

Cervical cancer screening allows detection and treatment of precancerous lesions before the development of cervical cancer. Screening is based on different algorithms, some allowing detection of HPV, and others identifying abnormal cells. Despite the role of high-risk HPV in cervical cancer, screening tests of cancer or precancerous lesions remain in many countries mainly based on the Papanicolaou (Pap) cytology test and do not include molecular virology tests [START_REF] Schiffman | Carcinogenic human papillomavirus infection[END_REF] . This is largely due to the low Positive Predictive Value (PPV) of current molecular tests. Indeed, because most of the current molecular diagnostic methods rely on the detection of HPV genome (DNA) and do not address the patterns of viral expression (RNA), they remain weak predictors of the evolution from low-grade squamous intraepithelial lesion (LSIL) to high-grade squamous intraepithelial lesion (HSIL) of the cervix [START_REF] Tornesello | Viral and cellular biomarkers in the diagnosis of cervical intraepithelial neoplasia and cancer[END_REF] . In addition, DNA identification of high-risk HPV is not fully predictive of cancer since only persistence for years of high-risk HPV is associated with an increased risk of cancer development [START_REF] Schiffman | Carcinogenic human papillomavirus infection[END_REF] .

Thus, the use of HPV DNA tests, as a screening assay, is currently increasing worldwide and shows high sensitivity [START_REF] Ogilvie | Effect of Screening With Primary Cervical HPV Testing vs Cytology Testing on High-grade Cervical Intraepithelial Neoplasia at 48 Months: The HPV FOCAL Randomized Clinical Trial[END_REF] but low PPV for HSIL detection [START_REF] Cuzick | Comparing the performance of six human papillomavirus tests in a screening population[END_REF] .

HPV RNA tests and in particular the detection of E6 and E7 mRNAs of high-risk HPV have been proposed as better molecular markers of cancer development, but E6 and E7 are also expressed during HPV transient infection so it remains difficult to define a threshold of expression associated with the persistence and evolution to high-grade lesions and cancer. There is no consensus that HPV RNA tests have a better diagnostic accuracy compared to HPV DNA tests and cytology for the detection of cervical precancerous lesions [START_REF] Virtanen | Performance of mRNA-and DNA-based high-risk human papillomavirus assays in detection of high-grade cervical lesions[END_REF][START_REF] Cook | Aptima HPV Assay versus Hybrid Capture® 2 HPV test for primary cervical cancer screening in the HPV FOCAL trial[END_REF][START_REF] Ge | Aptima Human Papillomavirus E6/E7 mRNA Test Results Strongly Associated With Risk for High-Grade Cervical Lesions in Follow-Up Biopsies[END_REF] . There is therefore a need for a novel generation of molecular diagnostic tests that can not only detect HPV infection, but also have the ability to accurately predict precancerous stages to offer a better and cost saving medical benefit [START_REF] De Thurah | Concordant testing results between various human papillomavirus assays in primary cervical cancer screening: systematic review[END_REF][START_REF] Hawkes | Not all HPV nucleic acid tests are equal: only those calibrated to detect high grade lesions matter for cervical screening[END_REF][START_REF] De Thurah | Not all HPV nucleic acid tests are equal: only those calibrated to detect high grade lesions matter for cervical screening: Response to "Concordant testing results between various human papillomavirus assays in primary cervical cancer screening: systematic review" by de Thurah, Bonde, Uyen, Lam and Rebolj[END_REF] .

We took advantage of Next-Generation Sequencing (NGS) technologies that now make it possible to study populations of transcripts as a whole, instead of focusing only on one or two specific messenger RNA as done with former techniques such as quantitative RT-PCR used in HPV RNA tests. To do so, we have developed a multiplexed amplification system targeting the virus splice junctions coupled with NGS analysis, tentatively named HPV RNA-Seq (based on the AmpliSeq technology), that allows to describe fine equilibrium among transcript species of 13 high-risk HPV (HPV16, 18, 31, 33, 35, 39, 45, 51, 52, 56, 58, 59, 66) plus 3 putative high-risk HPV (HPV68, 73, 82), in a single reaction. This molecular approach makes, in particular, possible to take a snapshot of the early vs late populations of HPV transcripts and to define a model based on a combination of reads that reflects the biology of the virus, which can then be correlated with the evolution of cervical lesions. The ultimate goal is to replace the current combination of cytology (Pap smear) and HPV molecular screening by a single molecular test for both the detection of high-risk or putative high-risk HPV and the triage of women at risk of transforming infection, before colposcopy.

In this proof-of-concept study conducted on 55 patients (27 HSIL, 28 LSIL), starting from cervical smears conserved at room temperature, we have shown that (i) HPV RNA-Seq can detect papillomaviruses with performances comparable to a HPV DNA-based reference diagnostic kit [START_REF] Heard | Clinical and analytical performance of the PapilloCheck HPV-Screening assay using the VALGENT framework[END_REF] , and (ii) a combination of the number of sequencing reads at specific early vs late HPV RNA spliced junctions can be used as a marker of high-grade cytology, with encouraging diagnostic performances as a triage test. Evaluation of transport medium for RNA conservation HPV16-positive cervical squamous cell carcinoma SiHa cells were cultivated and inoculated at a final concentration of 7x10 4 cells/mL in four transport medium: PreservCyt Solution (Hologic, USA), NovaPrep HQ+ Solution (Novaprep, France), RNA Protect Cell Reagent (Qiagen, Germany) and NucliSens Lysis Buffer (BioMérieux, France). The mixtures were aliquoted in 1mL tubes and kept at room temperature for 2 hours (D0), 48 hours (D2), 168 hours (D7), 336 hours (D14) and 504 hours (D21). In parallel, 7x10 4 cells pellets without transport medium were kept frozen -80°C for 2 hours, 48 hours, 168 hours, 336 hours and 504 hours as a control. At D0, D2, D7, D14 and D21, room temperature aliquots were centrifuged, the medium removed, and the pellets were frozen -80°C for a short time (<1h) before proceeding with RNA extraction. In the particular case of the NucliSens Lysis Buffer since the cells were lysed, the entire 1mL aliquot was frozen -80°C for a short time without prior centrifugation. For each sample, RNA was extracted using the PicoPure RNA Isolation kit (Thermo Fisher Scientific, USA), together with the corresponding (time match) frozen control, so that all samples have undergone one freezing cycle. RT-qPCR was performed to quantify the expression of the two human genes G6PD (forward primer: TGCAGATGCTGTGTCTGG; reverse primer: CGTACTGGCCCAGGACC) and GAPDH (forward primer: GAAGGTGAAGGTCGGAGTC; reverse primer: GAAGATGGTGATGGGATTTC) and the expression of the two viral genes HPV16 E6 (forward primer: ATGCACCAAAAGAGAACTGC; reverse primer: TTACAGCTGGGTTTCTCTAC) and E7 (forward primer: GTAACCTTTTGTTGCAAGTGTGACT; reverse primer: GATTATGGTTTCTGAGAACAGATGG) (Supplemental Figure 1). RNA integrity was assessed on a Bioanalyzer instrument (Agilent, USA) (Supplemental Figure 2).

MATERIALS AND METHODS

Ethics

HPV selection and splice sites analysis

HPV reference clones made available by the International Human Papillomavirus Reference Center (Karolinska University, Stockholm, Sweden) served as reference genomes, except for HPV68 which was retrieved from Chen et al. [START_REF] Chen | Evolution and Taxonomic Classification of Alphapapillomavirus 7 Complete Genomes: HPV18, HPV39, HPV45, HPV59, HPV68 and HPV70[END_REF] . Accession numbers used in this study were: K02718 (HPV16), X05015 (HPV18), J04353 (HPV31), M12732 (HPV33), X74477 (HPV35), M62849 (HPV39), X74479 (HPV45), M62877 (HPV51), X74481 (HPV52), X74483 (HPV56), D90400 (HPV58), X77858 (HPV59), U31794 (HPV66), KC470267 (HPV68), X94165 (HPV73) and AB027021 (HPV82). Multiple alignment of HPV genomes was done with ClustalW v2.1 using Geneious 19 v10. Previously known splice donor (SD) and splice acceptor (SA) sites for HPV16 [START_REF] Zheng | Papillomavirus genome structure, expression, and post-transcriptional regulation[END_REF] and HPV18 [START_REF] Wang | Construction of a full transcription map of human papillomavirus type 18 during productive viral infection[END_REF] were reported on the alignment, and predictions of unknown SD and SA sites were done manually for the other genotypes by sequence analogy (Figure 1A and 1B).

HPV RNA-Seq AmpliSeq custom panel A custom AmpliSeq panel was designed to be used on both PGM and Ion Proton instruments (Thermo Fisher Scientific). Five categories of target sequences were defined as follow:

HPV spliced junctions (sp): a set of target sequences, which are specific HPV splice events, involving a pair of splice donor (SD) and splice acceptor (SA) sites. The nomenclature includes a "sp" tag. For example, "31_sp_1296_3295_J43-46" stands for HPV31 (31), splice junction (sp), SD at position 1296 on HPV31 genome, SA at position 3295 on HPV31 genome, and junction (J) at position 43-46 on amplicon. The junction coordinates are given in a 4-bases interval, where the first 2 bases correspond to the donor part (or left part) and the last 2 bases to the acceptor part (or right part) of the sequence.

HPV unspliced junctions (unsp): a set of target sequences which are specific HPV genomic regions spanning either SD or SA sites, in the absence of any splice event. The nomenclature includes an "unsp" tag. For example, "31_unsp_1296_1297_J43-46" stands for HPV31 (31), unspliced (unsp), last base of the left part of the amplicon at position 1296 on HPV31 genome, first base of the right part of the amplicon at position 1297 on HPV31 genome, junction (J) at position 43-46 on amplicon. In this context, the term 'junction' refers to the exon-intron interface (ie the position where a donor or acceptor site would be found in case of a splice event), and the associated junction coordinates are used to characterize unspliced sequences bioinformatically as described in section "Sequencing data processing".

HPV genome away from spliced junctions (gen): a set of target sequences which are specific HPV genomic regions, away from any SD or SA sites. The nomenclature includes a "gen" tag. For example, "45_gen_1664_1794_NoJ" stands for HPV45 (45), HPV genomic region (gen), amplicon coordinates from position 1664 to position 1794 on HPV45 genome.

HPV-human fusion sequences (fus): a set of hypothesis-driven viral-cellular fusion transcripts, based on previous descriptions [START_REF] Wentzensen | Characterization of viral-cellular fusion transcripts in a large series of HPV16 and 18 positive anogenital lesions[END_REF][START_REF] Tang | The landscape of viral expression and host gene fusion and adaptation in human cancer[END_REF][START_REF] Peter | MYC activation associated with the integration of HPV DNA at the MYC locus in genital tumors[END_REF][START_REF] Lu | Multiple-integrations of HPV16 genome and altered transcription of viral oncogenes and cellular genes are associated with the development of cervical cancer[END_REF][START_REF] Kraus | The majority of viral-cellular fusion transcripts in cervical carcinomas cotranscribe cellular sequences of known or predicted genes[END_REF] . For each HPV, 18 fusion sequence candidates involving SA2 or putative breakpoint 1 or 2 (put. bkpt, see Figure 1B) for the viral part, and specific exons from MYC or PVT1 oncogenes for the cellular part, were added to the design. For example, "18_fus_929_MYC_001_exon3_J37-40" stands for HPV18 (18), candidate fusion transcript (fus), break/fusion at position 929 on HPV18 genome, fused with MYC mRNA isoform 001 exon 3, junction (J) at position 37-40 on amplicon.

Human sequences (hg): a set of 30 human sequences used as internal controls retrieved from publically available AmpliSeq projects and representing housekeeping genes (ACTB, B2M, GAPDH, GUSB, RPLP0), epithelial markers (KRT10, KRT14, KRT17), oncogenes, tumor suppressor genes, and direct or indirect downstream effectors of HPV oncoproteins (AKT1, BCL2, BRAF, CDH1, CDKN2A, CDKN2B, ERBB2, FOS, HRAS, KRAS, MET, MKI67, MYC, NOTCH1, PCNA, PTEN, RB1, STAT1, TERT, TOP2A, TP53, WNT1). The nomenclature for these sequences includes an "hg" tag. For example, "hg_TOP2A_E21E22" stands for human topoisomerase 2A mRNA exon 21-22.

In total, 750 target sequences were included into the panel (Table 1) and can be amplified with a pool of 525 unique primers (Supplemental Table S1). The custom panel is registered under number WG_WG00141 (Ion AmpliSeq Designer). The average amplicon size of the panel (primers included) is 141bp (range: 81-204bp). A detailed table including amplicons names and characteristics along with their corresponding primers and amplicon sequences is given in the Supplemental Table S1.

Study participants

Study participants were women aged from 25 to 65 years old referred for colposcopy consultation in French hospitals. The patients were referred for colposcopy in the context of a LSIL or a HSIL result at their cytology test performed in accordance with French recommendations regarding the cervical cancer screening program. Patients provided written informed consent according to French legislation.

Specimen collection

Genital samples were collected just before performing colposcopy using a cervical sampling device, immersed and rinsed in a vial filled with 20mL of PreservCyt Solution (Hologic), and sent at room temperature to the HPV National Reference Center (CNR) at Institut Pasteur, Paris, France. From July 2014 to April 2015, 84 patients were enrolled in the study, coming from 3 different French centers: CHU Angers (n=66); CHU Kremlin-Bicêtre (n=10); CHU Tours (n=6). Samples were removed from the study because of technical reasons (sample leakage, n=1), legal issues (n=7) or because they were used for initial technical tests (RNA conservation, RNA extraction and amplification, n=4). The remaining 72 samples (HSIL=37; LSIL=35) were processed (Supplemental Table S2).

Data collection

The following bio-clinical data were collected: date and results of the cytology test, age at the time of the cytology test, date and results of all available histological results posterior to colposcopy. As colposcopy was performed in the context of routine healthcare, biopsies were not performed in case of normal colposcopy.

HPV DNA detection using the PapilloCheck Test Kit (HPV DNA) Upon reception at CNR, 16 mL of cytological sample were transferred into a 50 mL Falcon tube and centrifuged at 4,500 g for 10 minutes. The supernatant was removed and the pellet washed with 1 mL of PBS. Sample was then centrifuged again at 5,000 g for 10 minutes and the supernatant removed.

The pellet was frozen at -80°C before DNA extraction. Following DNA extraction (Macherey Nagel, Germany), HPV detection was done using the PapilloCheck Test Kit (Greiner Bio-One GmbH, Germany) according to manufacturer instructions (Supplemental Table S2).

RNA extraction and characterization

In parallel to the HPV DNA procedure, 3 x 1 mL aliquots of cytological specimen were centrifuged at 14,000 rpm for 7 minutes, the supernatant was removed and the pellet was washed with 1 mL of PBS.

Sample was then centrifuged again at 14,000 rpm for 7 minutes and the supernatant removed. The pellet was frozen at -80°C before RNA extraction. RNA extractions were done using the PicoPure RNA Isolation kit (Thermo Fisher Scientific), including on-column DNAse treatment, with a final elution volume of 30 µL. Total RNA was quantified on a Nanodrop (Life Technologies, USA) and RNA integrity was evaluated on a Bioanalyzer RNA 6000 pico chip (Agilent) using the RIN (RNA Integrity Number), a quality score ranging from 1 (strongly degraded RNA) to 10 (intact RNA). For each sample, RT-qPCR targeting mRNA from housekeeping genes ACTB (forward primer: CATCGAGCACGGCATCGTCA; reverse primer: TAGCACAGCCTGGATAGCAAC; amplicon size = 210bp), and GAPDH (forward primer: GAAGGTGAAGGTCGGAGTC; reverse primer: GAAGATGGTGATGGGATTTC; amplicon size = 226bp) were done in a SYBR Green format with 45 cycles of amplification. RT-negative (RT-) PCR were run to evaluate the presence of residual DNA after RNA extraction (Supplemental Table S2).

Amplification and sequencing

Starting from RNA, cDNA were generated using the SuperScript III (n=17 samples) or SuperScript IV (n=55 samples) first strand synthesis system (Thermo Fisher Scientific) with random hexamers and a final RNAse H treatment. Libraries were prepared using the Ion AmpliSeq Library Kit 2.0 and AmpliSeq custom panel WG_WG00141, with 21 cycles of amplification before adapter's ligation. Each sample was barcoded individually. Only positive libraries were sequenced (Supplemental Table S2). In total, 55 clinical samples plus 1 cellular model (SiHa) were sequenced on 4 Ion Proton runs. Raw data (.fastq files) are available on the NCBI SRA database under BioProject accession number PRJNA525642.

Sequencing data processing

Reads were aligned to the reference sequences of the amplicons using STAR 27 v2.5.3a in local alignment mode (parameter --alignEndsType EndToEnd), by only reporting uniquely mapped reads (--outFilterMultimapNmax 1) and turning off splicing alignment (--alignIntronMax 1). The expression of each amplicon was evaluated by the number of sequencing reads uniquely mapping to their respective sequence (read counts). For reference sequences containing a splice junction, only reads mapping at the junction site and encompassing at least 10 bases before and 10 bases after the junction were kept. Read counts for each sequence and each sample are provided in the Supplemental Table S3.

HSIL prediction model

Selection of amplicons

Read counts were normalized by the size of the library (each read count was divided by a ratio of the library size for a given sample to that of the average library size across samples) and the 215 amplicons capturing splice junctions (sp) of the 16 high-risk or putative high-risk HPV were selected.

These amplicons have been annotated with generic names with respect to the type of transcripts they capture, which are shared across HPV species (e.g. "SD1-SA1", see Figure 1B and Supplemental Table S1). Amplicons capturing splice junctions conserved across the 16 HPV species were summed up, leading to the definition of 18 variables used as predictors in the model. 33 out of the 55 clinical samples have been selected as presenting enough coverage of these specific amplicons (20 monoinfected and 13 multi-infected samples). The remaining 22 samples of the dataset were not used in the logistic regression analysis because they had missing or too low expression signal at spliced junctions for the prediction, reflecting for example HPV-negative samples.

Logistic regression model

Calling high-grade cytology Y as taking the value 1 for high-grade (HSIL) and 0 for low-grade (LSIL), and a set of amplicons x, a logistic regression model was used to predict the probability that a given observation belongs to the "1" class versus the probability that it belongs to the "0" class. Logistic regression models the log odds of the event (here the grade of the cytology) as a function of the predictor variables (here the amplicon expression estimated by its read count). Formally, the logistic regression model assumes that the log odds is a linear function of the predictors:

logit (p ) =ln p 1-p ae è ç ö ø ÷ = b 0 + b t x
where p = P( y =1| x) indicates the probability of the event (being of high-grade), β i are the regression coefficients, and x i the explanatory variables, in our case the log2 number of reads mapping to the amplicons.

Solving for π, this gives: p =

1 1+ e -(b 0 +b x)
Implementation of the logistic regression model

To limit overfitting, we used L2-norm (ridge) regularization, which allows shrinking the magnitudes of the regression coefficients such that they will better fit future data. We estimated the logistic model using the R (http://www.r-project.org/ ; last accessed on January 29, 2019) package glmnet [START_REF] Friedman | Regularization Paths for Generalized Linear Models via Coordinate Descent[END_REF] . Leaveone-out (LOO) cross-validation was used to pick the regularization parameter λ, the one that gives minimum mean cross-validated misclassification error was used. Using λ as the regularization parameter, the model output consisted in an estimate of a coefficient value β for each variable in the logistic regression model. This model was then used to predict the grade of the multi-infected observations, by treating each HPV species separately.

Training set and test set

The model was built upon the clinical outcome LSIL or HSIL obtained from the cytological analysis, and estimated on a training set consisting of 20 mono-infected samples (5 LSIL and 15 HSIL) in order to avoid a confusion bias. It is indeed anticipated that, in the case of multi-infected samples, several HPV could contribute differently to the progression of the lesion, or to a mix of several grades within the same sample, because they are engaged in different stages of their cycle. The performance of the model was then evaluated on a test set consisting of 13 multi-infected samples. In this case, the set of amplicons of each HPV species was used separately to classify the multi-infected samples, to get one prediction per HPV, as done for the mono-infected samples. For example if a sample had expression of amplicons from both HPV16 and HPV32, two predictions were given: one using only sequencing reads mapping to HPV16, and one using only sequencing reads mapping to HPV32. Like this it became possible to interpret the results finely from a virological point of view, as we could discriminate which HPV was responsible of the lesion.

RESULTS

Evaluation of transport medium for RNA conservation

The stability of total RNA from cervical cells at room temperature was evaluated in four solutions:

PreservCyt (Hologic), the most widely used solution for gynecological specimen collection; NovaPrep HQ+ Solution (Novaprep), a competitor product used for cells and DNA recovery but never evaluated for RNA conservation; RNA Protect Cell Reagent (Qiagen), a popular solution for RNA stability; and

NucliSens Lysis Buffer (BioMérieux), a lysis buffer part of the NucliSens automated acid nucleic procedure which has been described as a RNA stabilizer (unpublished data). The amount of spiked HPV16-positive cervical squamous cell carcinoma cells (SiHa) was calibrated to be representative of a cervical smear. After 48h at room temperature, RT-qPCR measurement of cellular and viral transcripts showed no or little RNA loss in PreservCyt, only limited RNA degradation (<1 log) in RNA Protect and NucliSens Lysis Buffer, and a marked RNA loss in NovaPrep HQ+ Solution (>2 log) (Supplemental Figure 1). After 7 days and up to 21 days, only the PreservCyt solution provided RNA quality with a limited RNA degradation pattern as indicated by the detection of 18S and 28S rRNA (Supplemental Figure 2).

We therefore decided to use the PreservCyt solution to collect the gynecological specimen of the study.

HPV RNA-Seq AmpliSeq custom panel

Transcriptomic maps known for HPV16 [START_REF] Zheng | Papillomavirus genome structure, expression, and post-transcriptional regulation[END_REF] and HPV18 [START_REF] Wang | Construction of a full transcription map of human papillomavirus type 18 during productive viral infection[END_REF] were used to predict unknown but likely splice donor and splice acceptor sites for HPV31, 33, 35, 39, 45, 51, 52, 56, 58, 59, 66, 68, 73, and 82 (Figures 1A and1B). The resulting reconstructed transcripts, as well as HPV genomic sequences, were used as a template for the design of an Ion AmpliSeq panel targeting 16 high-risk or putative high-risk HPV and named HPV RNA-Seq. Putative breakpoints in HPV genomes, and 30 human cellular genes used as internal controls, were also added to the design. In total, 750 sequences are targeted by a single mix made of 525 unique primers (Table 1 and Supplemental Table S1).

Samples, RNA & sequencing

72 gynecological samples (HSIL=37; LSIL=35) coming from 3 different French centers (Angers, Kremlin-Bicêtre and Tours) and collected in PreservCyt solution were processed with RNA extraction using a method designed to recover total RNA from as little as a single cell (PicoPure RNA Isolation kit, Thermo Fisher Scientific) (Supplemental Table S2). In most of the cases total RNA was measurable using a Note, the presence of residual cellular DNA or HPV DNA in RNA preparation is not a major concern since the AmpliSeq assay can differentiate between HPV transcripts and genomic sequences. AmpliSeq libraries were initiated from total RNA and were positive after 21 cycles of amplification for 55 samples (i.e. detectable on a Bioanalyzer HS DNA chip). Attempts to add one or two amplification cycles did not bring any significant improvement to the results (data not shown).

In total, 55 patients (HSIL=27; LSIL=28), plus SiHa HPV16-positive cells as a control, were sequenced on Ion Proton. The sequencing reads were aligned to the target sequences and read counts were generated (Supplemental Table S3). An average of 2.4 million usable reads per sample was reached (min=0.02M; max=8.36M), among which an average of 2.1 million reads mapped to the human sequences (hg) used as internal controls (min=0.01M; max=8.06M) (Supplemental Table S3). The detection of highly expressed human sequences in all samples, even though inter-sample variations were observed, contributed to validate the sequencing procedure, which is important especially for the interpretation of HPV-negative samples. Rare non-zero values were also observed for some of the numerous HPV-human fusion sequences (fus) that were hypothesized (Supplemental Table S3) but were all false positives, identified as such because only half of the reference sequences were covered by reads.

HPV RNA-Seq used for HPV detection and genotyping

The first application of HPV RNA-Seq is to detect the presence in a given sample of any of the 16 highrisk or putative high-risk HPV targeted by the panel. The number of reads mapping to HPV-specific amplicons (i.e. the sum of categories "sp", "unsp" and "gen") was used to detect the presence of a given HPV genotype. To help determining a threshold for detection, we took as a reference a HPV DNA test validated for clinical use (PapilloCheck, Greiner Bio-One GmbH). The best sensitivity and specificity values between the two tests were obtained for threshold of 100-200 reads (Figure 2). For example, a threshold value of 150 reads resulted in a Sensitivity (Se (HPV-DNA) ) of 97.3%, a Specificity (Sp (HPV-DNA) ) of 83.3%, leading to a Positive Predictive Value (PPV (HPV-DNA) ) of 92.3% and a Negative Predictive Value (NPV (HPV-DNA) ) of 93.8% for detecting high-risk HPV in this population composed of around 50% of HSIL and 50% of LSIL (Table 2, raw data in the Supplemental Table S2 and Supplemental Table S3). A more detailed view of the genotypes identified by both techniques is given in Figure 3. The number of mono-infected, multi-infected, or HPV-negative samples identified by the two tests is summarized in Using a threshold value of 150 reads, HPV RNA-Seq detected two more positive patients than the HPV DNA test (n=39 vs n=37, Table 2). HPV RNA-Seq identified the presence of more than one HPV for three more patients than the HPV DNA test (n=13 vs n=10 multi-infected samples, Table 3). Globally, HPV16 was found at a slightly weaker occurrence by HPV RNA-Seq (n=18 vs n=19) in favor of other genotypes such as HPV31, 33, 45, 52, 56, 58 or 66 which were less commonly found by the HPV DNA test (HPV31 n=5 vs n=4; HPV33 n=3 vs n=1; HPV45 n=3 vs n=2; HPV52 n=5 vs n=3; HPV56 n=4 vs n=2; HPV58 n=5 vs n=4; HPV66 n=2 vs n=1, Figure 3). Apart from HPV16, only HPV51 was less frequently found by HPV RNA-Seq than by HPV DNA (n=2 vs n=3). The cellular model (SiHa) gave only HPV16 signal in both tests, as expected (Supplemental Table S3).

HPV RNA-Seq used as a marker of high-grade cytology

We conducted an exploratory analysis on 20 of the mono-infected samples in which we showed that HPV RNA spliced junctions could be used to predict high-grade cytology. We focused the analysis on amplicons capturing splice junctions (category "sp") to be sure to detect HPV transcripts. However, the number of mono-infected samples (n=20) used as training set was small, in particular the number of samples of LSIL (n=5). In this configuration, we were not able to perform a fully accurate variable selection, i.e. to select the strict minimum number of amplicons that were necessary for HSIL versus LSIL prediction and set the others to zero coefficient. In addition, we were not able to avoid overfitting, as using only 5 LSIL and 15 HSIL samples did not allow capturing the diversity of the whole population. LOO cross-validation was used to pick the lambda giving the minimum cross-validated error using ridge regularization. Lambda = 0.08 gave a mean cross-validated error of 15%. We also computed a 20% prediction error using nested cross-validation. This error rate can be seen as an indicator of how the model could fit future datasets. We used the corresponding parameter to fit a regularized logistic regression model, assigning a coefficient to each amplicon (Table 4) and a probability of being of high-grade to each sample (Table 5). The grade of the 20 mono-infected samples was classified correctly, except for one observation (Table 5). It is interesting to note that this unique misclassified sample (IonXpress_019_2613), which was classified LSIL by the cytological analysis, was further found as containing a mixture of LSIL and HSIL lesions after histological examination performed more than one year after the sampling done for HPV RNA-Seq/cytology.

The estimated model was then used to classify the 13 multi-infected samples, with each HPV species present within one sample being classified individually for its implication in HSIL development. If at least one HPV species gave a HSIL prediction, the sample was considered to be HSIL. We calculated performances for HSIL prediction for all samples, considering as not being of high-grade both the six samples without sufficient coverage of the splice junctions and the 16 HPV-negative samples not exceeding the threshold of HPV detection. The calculated performances for HSIL prediction in comparison to cytology for the 55 patients (mono-infected, multi-infected and HPV-negative) were Se (cyto) =66.7%, Sp (cyto) =85.7%, PPV (cyto) =81.8% and NPV (cyto) =72.7% (Table 6A). The performances were also calculated for the subset of 39 samples having at least one HPV identified by HPV RNA-Seq, giving in this case Se (cyto/HR+) =94.7%, Sp (cyto/HR+) =80.0%, PPV (cyto/HR+) =81.8% and NPV (cyto/HR+) =94.1% (Table 6B).

Note that the ratio HSIL to LSIL remained similar between these two populations (around 1:1), making the comparison of the PPV and the NPV possible. Finally a summary of the results for HPV detection and genotyping (HPV RNA-Seq vs HPV DNA) and high-grade cytology prediction (HPV RNA-Seq vs cytology), including posterior histological data of cervix biopsies when available, is presented in Table 7.

HPV RNA-Seq used as a triage test

The performances of HPV RNA-Seq as a triage test were evaluated using histology as gold standard.

Results from histological examination were, however, not available for all patients. The time interval separating HPV RNA-Seq/cytology tests from histological analysis, varying between 0 and 780 days, was another limitation in this study. To try to overcome these drawbacks, we compared the performances of HPV RNA-Seq vs histology to the performances of cytology vs histology, considering either all available samples (Supplemental Table S4 A), or only samples for which histology was done less than 3 months after HPV RNA-Seq/cytology (Supplemental Table S4 B), or only samples for which histology was done less than 6 months after HPV RNA-Seq/cytology (Supplemental Table S4 C). In addition and for each category, we made the distinction between the performances obtained when HPV RNA-Seq HPV-positive and HPV-negative patients were grouped together (Supplemental Table S4 1&1'), or when only HPV-positive patients were considered (Supplemental Table S4 2&2'). Calculation of the PPV as a function of HSIL prevalence in the population was also done (Supplemental Figure 3 and Supplemental Table S4). All these results are given as Supplemental Data.

DISCUSSION

We have developed a highly-multiplexed RT-PCR assay coupled with Next-Generation Sequencing (HPV RNA-Seq) combining HPV detection and genotyping together with predicting high-grade cytology, starting from cervical specimens conserved at room temperature. A pilot study was conducted on 55 patients.

The performances of HPV RNA-Seq used as a HPV detection and genotyping assay were evaluated in comparison to the HPV DNA PapilloCheck kit (HPV DNA), which is officially approved for clinical use. A good concordance of the results was observed between the two assays (Area Under Curve > 0.95, Figure 2). A positive threshold of 150 reads resulted in high sensitivity and negative predictive value of HPV RNA-Seq (Se (HPV-DNA) =97.3%, NPV (HPV-DNA) =93.8%, Table 2), along with a relatively high but lower specificity and positive predictive value (Sp (HPV-DNA) =83.3%, PPV (HPV-DNA) =92.3%, Table 2) linked to the identification of additional genotypes by HPV RNA-Seq not detected by HPV DNA. Because cervical samples were split before independent extractions of RNA (HPV RNA-Seq) and DNA (HPV DNA), the few differences observed between the two tests can reflect a non-homogeneous distribution of infected cells. It is also important to note that PapilloCheck, like other HPV DNA tests, is not 100% accurate [START_REF] Dalstein | Analytical evaluation of the PapilloCheck test, a new commercial DNA chip for detection and genotyping of human papillomavirus[END_REF][START_REF] Klug | Comparison of the performance of different HPV genotyping methods for detecting genital HPV types[END_REF] , so it remained difficult to identify potential false positive results of HPV RNA-Seq versus better sensibility. For example, three patients were classified as HPV-negative by Papillocheck but not by HPV RNA-Seq. The number of RNA-Seq reads associated to HPV species in these three potential false HPV-positive patients was close to the limit of detection for some of them (<= 400 reads) but not for all (e.g. 39527 reads mapped to HPV58 for sample 2065) (Supplemental Table S3). The calculated sensitivity and specificity may therefore not reflect optimally the added value of HPV RNA-Seq. These limitations are common for any novel diagnostic test when compared to older references.

Effective cervical cancer screening requires high Se and NPV for high-risk HPV detection, as women with a negative HPV test are usually tested again only after several years. The positive threshold for HPV genotyping was set at 150 reads in this study because it optimized both Se and Sp values, but lowering this threshold in order to maximize the sensitivity remains possible. Such adjustments will be possible following the study of larger cohorts.

As a second application of HPV RNA-Seq, as a triage test, a logistic regression model for the prediction of high-grade cytology was built based on a combination of the number of reads captured at specific HPV RNA spliced junctions, using the grade found by cytology as a reference. This evaluation was conducted in a population of women with LSIL or HSIL cytology results. Where at least one HPV was given a high-grade signature, the patient's prediction was set as "HSIL". Conversely, "Not HSIL" was used when either no HPV was detected in the sample (threshold of 150 reads) or when none of the genotypes detected by HPV RNA-Seq were given high-grade prediction (absence of detectable transcripts). We preferred to use the terminology "Not HSIL" rather than "LSIL" because the protocol did not allow the comparison of the HPV-DNA positive samples evaluated as LSIL with the ones evaluated as normal in cytology. Also, because there is a possibility that cervical lesions could in some rare cases originate from causes other than HPV infections, the use of "Not HSIL" instead of "LSIL" in the case of HPV-negative samples seemed more appropriate.

As far as the comparison with cytology could be used as a benchmark, when the 55 patients were considered (including mono-infected, multi-infected and HPV-negative), the number of HSIL predicted by HPV RNA-Seq (n=22) was lower than the number of HSIL identified by cytology (n=27), resulting in Se (cyto) =66.7%, Sp (cyto) =85.7%, PPV (cyto) =81.8% and NPV (cyto) =72.7% (Table 6A). Interestingly, when only HPV RNA-Seq HPV-positive samples were considered, the PPV(cyto/HR+) for the detection of high grade lesions remained unchanged but the Se(cyto/HR+) and the NPV(cyto/HR+) increased to 94.7% (+28.0) and 94.1% (+21.4) respectively, with the number of HSIL predicted by HPV RNA-Seq (n=22) becoming superior to the number of HSIL identified by cytology (n=19) (Table 6B). In this case the only one patient identified HSIL by cytology but not predicted HSIL by HPV RNA-Seq (sample IonXpress_020_3137) was found HSIL by the histological examination done 350 days later, which

argument in favor of a potential added medical value of HPV RNA-Seq over cytology, although studies on larger cohorts are now required. Another observation is that the Se of HPV RNA-Seq vs histology was always higher on the subset of HPV-positive patients (+12.0 to +33.3, Supplemental Table S4), similarly to the evaluation done with cytology taken as reference (Table 6). Lastly, the Sp of HPV RNA-Seq vs histology increased on the subset of patients for whom histology was performed less than 3 months after sampling (+5.0 to +8.3, Supplemental Table S4) but decreased on the subset of patients for whom histology was performed less than 6 months after sampling (-1.7 to -2.3, Supplemental Table S4), which may indicate that some lesions have evolved in the meantime.

Although the minimum number of reads required for the assay was not evaluated, our observations tend to support that 1 million reads or less per sample is enough for performing HPV genotyping, but more would be needed for HPV transcripts detection. We do realize for example that the absence of detectable transcripts for a given HPV was assimilated to the absence of HPV transcripts in the sample, which may not be true if sequencing depth was insufficient. Generally speaking, the questions of the suggesting that the format of the test can also evolve to keep only the most informative target sequences, while potentially reducing the depth of sequencing required for analysis and the associated costs. For example, a reduction in the number of human sequences (hg) used as internal controls could be considered.

HPV RNA-Seq will be further developed and validated as a companion test in HPV DNA-positive patients or when the result of cytology is uncertain, in order to allow focusing the colposcopies to the most relevant patients. It has recently been shown that only one third of women recommended for colposcopy after primary HPV testing (DNA) and cytology had actually HSIL [START_REF] Ogilvie | Effect of Screening With Primary Cervical HPV Testing vs Cytology Testing on High-grade Cervical Intraepithelial Neoplasia at 48 Months: The HPV FOCAL Randomized Clinical Trial[END_REF] . By increasing the positive predictive value in detecting HSIL, HPV RNA-Seq could significantly increase the medical benefit-cost ratio of colposcopies. The case of Atypical Squamous Cells of Undetermined Significance (ASCUS) would also constitute an important patient's category to demonstrate an added value of the assay.

Once the performances and the medical benefit have been evaluated on large cohorts, such broad range genotype papillomavirus transcriptome assay could eventually replace first line cytology and DNA-based tests, by providing in a single procedure both HPV detection and genotyping together with a molecular marker of high-grade lesions. Other diagnostic applications in HPV-associated anogenital or head and neck cancers can also be envisioned.

In conclusion, HPV RNA-Seq can provide a second line test in HPV-positive patients in order to reduce unnecessary colposcopies and even be used as a two-in-one test combining HPV typing with triage capabilities. HPV RNA-Seq is minimally-invasive and is convenient for sample conserved at room temperature. The assay will now require further clinical validation in larger cohorts. Multi-infected samples 13 10 HPV-negative samples 16 18 Number of mono-infected, multi-infected and HPV-negative samples by HPV RNA-Seq at threshold value of 150 reads, vs HPV DNA (PapilloCheck). Data: 55 patients. The first and fourth columns give the id of the splice junction captured by the amplicon, the second column gives the coefficient assigned by the logistic regression and the third column indicates whether the splice junction comes from a "late" or "early" transcript. The first column gives the sample id, the second column gives the probability estimate that the sample is HSIL, the third and fourth columns give the corresponding prediction, the fifth column contains TRUE if the prediction is consistent with the grade evaluated by cytology. HPV genotypes included in the scope of the HPV DNA test (PapilloCheck) but not in HPV RNA-Seq are indicated into brackets. For each genotype identified by HPV RNA-Seq, a classification is given: either "Not enough coverage on splice junctions" (no prediction was possible for the genotype), "Not HSIL"

or "HSIL". When at least one HPV was given high-grade signature, the patient's prediction was set as "HSIL". Conversely, a final "Not HSIL" means that either no HPV was detected in the sample, or that none of the HPV genotypes detected were given HSIL prediction.

Supplemental

  approval and consent to participate: This work was approved by the Comité de Protection des Personnes Ile de France 1 (CPP IDF1) and by the Agence Nationale de Sécurité des Médicaments et des Produits de Santé (ANSM). The data processing was authorized by the Commission Nationale de l'Informatique et des Libertés (CNIL). Patients provided written informed consent to participate in the study.

Nanodrop ( 70 /

 70 72 positive, average on positive RNA eluates = 18 ng/µL) and was detectable on a Bioanalyzer pico RNA chip with a pattern indicating RNA degradation (63/72 positive, average RNA Integrity Number on positive = 2.2). RT-qPCR performed for all samples on ACTB mRNA (amplicon size = 210bp) and GAPDH mRNA (amplicon size = 226bp) indicated that RNA quality was compatible with amplification of 200-250bp size fragments (ACTB mRNA average Ct=27.8; GAPDH mRNA average Ct=30.1). Samples that failed passing this initial RT-PCR quality control were not sequenced. qPCR performed after omitting the reverse transcription step (RT-) were also run and showed in general no or little traces of residual genomic DNA (ACTB DNA average Ct=38.4 ; GAPDH DNA average Ct=35.6).

  format of the test and of the model of use are of importance in the perspective of deploying NGSbased in-vitro diagnostic tests. AmpliSeq, former product by Thermo Fisher Scientific developed for Ion Proton and PGM instruments, has been transferred in 2018 to Illumina and is now fully compatible with Illumina sequencers. In a decentralized laboratory model, it becomes possible that 4-6 samples could run on a benchtop iSeq100 sequencer for a cost per sample around 200 USD, with RNA extraction, quality controls and data analysis included (salaries and equipment excluded). In a more centralized view where all regional samples would converge to one laboratory, the use of productionscale sequencers such as the HiSeq or NovaSeq instruments could allow multiplexing up to 381 samples per run, potentially reducing the cost per sample to around 10-20 USD and thus making NGSbased assays competitive over PCR-based tests. Another point is that not all sequences of HPV RNA-Seq contributed equally to the result, with some of them giving useless or redundant information,

Figure 1 :

 1 FIGURES AND TABLES

Figure 2 :

 2 Figure 2: Receiver Operating Characteristic (ROC) curve. HPV DNA (PapilloCheck) was used as a

Figure 3 :

 3 Figure 3: Comparison of the number of HPV genotypes identified by HPV RNA-Seq and HPV DNA.

Figure 1 :

 1 Evaluation of transport medium for RNA conservation: detection of cellular and viral transcripts by RT-qPCR at D0, D2, D7, D14 and D21. (A) Cellular transcripts GAPDH and G6PD

Table 3 .

 3 Note that, because the HPV DNA test can detect the 16 high-risk or putative high-risk HPV

	captured by HPV RNA-Seq (HPV16, 18, 31, 33, 35, 39, 45, 51, 52, 56, 58, 59, 66, 68, 73, and 82) plus 8
	additional low-risk HPV (HPV6, 11, 40, 42, 43, 44/45, 53 and 70), the comparison was based only on
	the 16 HPV common to both tests.

Table 1 :

 1 HPV RNA-Seq AmpliSeq custom panel contents.

	692
			sp	unsp	gen	fus	hg
		HPV16	14	11	4	18	0
		HPV18	18	12	4	18	0
		HPV31	14	11	4	18	0
		HPV33	13	9	4	18	0
		HPV35	14	10	4	18	0
		HPV39	10	8	4	18	0
		HPV45	14	10	4	18	0
		HPV51	10	9	4	18	0
		HPV52	16	11	4	18	0
		HPV56	16	10	4	18	0
		HPV58	13	8	4	18	0
		HPV59	14	8	4	18	0
		HPV66	15	8	4	18	0
		HPV68* 10	9	4	18	0
		HPV73* 13	10	4	18	0
		HPV82* 11	9	4	18	0
		human	0	0	0	0	30
		TOTAL	215	153	64	288	30	750
	693	The number of target amplicons is indicated for each category (sp, unsp, gen, fus, hg) and for each
	694	viral and cellular origin. Putative high-risk HPV are indicated by a star (*).
	695

Table 2 :

 2 Performances of HPV RNA-Seq for HPV detection.

	696

Table 3 :

 3 Comparison of HPV RNA-Seq and HPV DNA for the classification of samples.

	HPV RNA-Seq HPV DNA

Table 4 :

 4 Coefficients of the (ridge) logistic regression.

	junction	coefficient	name_transcript_category name_transcript_contents
	(Intercept)	0.468298365
	SD2_SA10 -0.693322203	late	L1
	SD3_SA4	0.545728771	early	(E1) E4 E5
	SD1_SA4	0.387642812	early	(E6) E2 E5
	SD2_SA4	-0.262522618	early	(E7) E2 E5

Table 5 :

 5 Classification results of the (ridge) logistic regression.

	sample	prediction_score prediction_class prediction_class prediction_accuracy
	IonXpress_039_115	0.115	-1	LSIL	TRUE
	IonXpress_033_730	0.204	-1	LSIL	TRUE
	IonXpress_038_114	0.259	-1	LSIL	TRUE
	1492	0.425	-1	LSIL	TRUE
	IonXpress_019_2613	0.562	1	LSIL	FALSE
	IonXpress_027_598	0.653	1	HSIL	TRUE
	729	0.716	1	HSIL	TRUE
	567	0.718	1	HSIL	TRUE
	IonXpress_018_2439	0.902	1	HSIL	TRUE

Table 6 :

 6 Performances of HPV RNA-Seq for the prediction of high-grade cytology. Seq vs cytology for HSIL detection (A) for the 55 patients and (B) for the 717 subset of 39 patients having at least one HPV identified by HPV RNA-Seq. Sensitivity (Se), Specificity 718 (Sp), Positive Predictive Value (PPV) and Negative Predictive Value (NPV) are given. "Not HSIL" means 719 that either no HPV was detected in the sample by HPV RNA-Seq or that none of the HPV genotypes 720

	A		Cytology
			HSIL LSIL	Se(cyto)	66.7%
		HSIL	18	4	Sp(cyto)	85.7%
	HPV RNA-Seq
		Not HSIL	9	24	PPV(cyto)	81.8%
					NPV(cyto)	72.7%
	B		Cytology
			HSIL LSIL	Se(cyto/HR+) 94.7%
	HPV RNA-Seq	HSIL	18	4	Sp(cyto/HR+) 80.0%
	HR+	Not HSIL	1	16	PPV(cyto/HR+) 81.8%
					NPV(cyto/HR+) 94.1%

Table 7 :

 7 HPV detection and genotyping and HSIL prediction for the 55 clinical samples. 723

				HPV RNA-Seq
			Genotyping		Marker of HSIL
			Per patient		Per HPV		Per patient
	Sample name	HPV DNA	Detection	Not enough	Not HSIL	HSIL	Prediction	Cytology	Histology	Time (days)
				coverage on						cyto-histo
				splice junctions
	D-15-0041_1066_BC13	16	16			16	HSIL	HSIL	HSIL	55
	D-15-0041_1122_BC14	16	16			16	HSIL	HSIL	HSIL	130
	D-15-0041_1124_BC5	16,39	16,39	39	16		Not HSIL	LSIL	HSIL	[70-434]
	D-15-0041_1490_BC6	16,39	16,35,39		39	16,35	HSIL	LSIL	HSIL	67
	D-15-0041_1492_BC7	16	16		16		Not HSIL	LSIL	LSIL	81
	D-15-0041_151_BC15	16,(53)	16			16	HSIL	LSIL	HSIL	130
	D-15-0041_152_BC16	16,(42)	16,52,82	16,52,82			Not HSIL	LSIL	LSIL	41
	D-15-0041_2209_BC11	16,(42),52	16,39,52	39	16,52		Not HSIL	LSIL	HSIL	n.d.
	D-15-0041_250_BC12	16,39,(42)	16,39		16,39		Not HSIL	LSIL	LSIL	55
	D-15-0041_25_BC4	16	16			16	HSIL	HSIL	HSIL	75
	D-15-0041_2612_BC8	16	16			16	HSIL	HSIL	n.d.	n.d.
	D-15-0041_567_BC9	16	16			16	HSIL	HSIL	HSIL	n.d.
	D-15-0041_610_BC2	16	16			16	HSIL	HSIL	HSIL	113
	D-15-0041_729_BC3	16	16			16	HSIL	HSIL	HSIL	59
	D-15-0041_752_BC10	16	16			16	HSIL	HSIL	HSIL	444
	IonXpress_017_2437	(43),51	51		51		Not HSIL	LSIL	LSIL	195
	IonXpress_017_251	neg	neg				Not HSIL	HSIL	LSIL	85
	IonXpress_018_2439	58	58			58	HSIL	HSIL	LSIL	164
	IonXpress_018_440	neg	neg				Not HSIL	LSIL	LSIL	38
	IonXpress_019_2613	16	16			16	HSIL	LSIL	HSIL	[416-780]
	IonXpress_020_3137	(53)	56	56			Not HSIL	HSIL	HSIL	350
	IonXpress_021_10	56,(44/55)	56		56		Not HSIL	LSIL	n.d.	130
	IonXpress_021_443	58	33,58	33		58	HSIL	HSIL	LSIL	99
	IonXpress_022_23	neg	neg				Not HSIL	HSIL	HSIL	n.d.
	IonXpress_022_444	16,33	16,33		33	16	HSIL	HSIL	HSIL	69
	IonXpress_023_24	(6),(11),(53)	neg				Not HSIL	HSIL	HSIL	[0-13]

Acknowledgments: The authors want to thank the Clinical Core of the Center of Translational Science (CRT-CC) of Institut Pasteur for the management of all legal and ethical aspects of the study, PathoQuest (Paris, France) and the Center for Translational Science (CRT) / Cytometry and Biomarkers Unit of Technology and Service (CB UTechS) at Institut Pasteur for providing access to the sequencers.

This work was funded by Institut Pasteur (grants InnovDiag and Pasteur Innov), Paris, France.

Availability of data and material:

The datasets supporting the conclusions of this article are included within the article and its additional files. The AmpliSeq custom panel is registered under number WG_WG00141. Raw sequencing data are available on the NCBI SRA database under BioProject accession number PRJNA525642.

opens the possibility that this sample might be positive if the patient would be tested again by HPV RNA-Seq at a date closer to the histological examination.

In clinical use, after primary screening for high-risk HPV in the general population, a triage test with high Sp and PPV is needed for the triage of women at risk of transforming infection before colposcopy.

In countries that have adopted HPV DNA as a screening test, cytological analysis can be used for the triage of women at risk because cytology has better Sp and PPV than HPV DNA tests [START_REF] Cuzick | Management of[END_REF]32 . In line with that, the PPV (cyto) =81.8% of HPV RNA-Seq outperformed HPV DNA and other RNA tests, whose PPV as triage assays never exceed 50% in a population of women referred for colposcopy composed with a similar 1:1 ratio of HSIL to LSIL [START_REF] Cuzick | Comparing the performance of six human papillomavirus tests in a screening population[END_REF] .

We thus sought to evaluate the added value of HPV RNA-Seq over cytology for the triage of women at risk of developing cervical cancer. To do so, histology was used as the gold standard for the diagnosis of cervical lesions. However, an inherent limitation of this work was that histology was not concomitant with the sampling carried out for HPV RNA-Seq/cytology tests, which means that by the time histology was performed (between 0 and 780 days after initial sampling), the lesion could have evolved spontaneously in one direction (LSIL to HSIL) or another (HSIL to LSIL). To help clarify this point we compared side to side the performances of HPV RNA-Seq vs histology to the performances of cytology vs histology, considering different categories of samples (Supplemental Table S4).

Remarkably, whatever the category considered, the Sp of HPV RNA-Seq vs histology was always greater than or equal to the Sp of cytology vs histology (+0.0 to +11.1), and the resulting PPV of HPV RNA-Seq vs histology was always greater than the PPV of cytology vs histology (+2.4 to +7.4 in this population reflecting others studies 9 ) (Supplemental Table S4). Calculation of PPV as a function of HSIL prevalence allowed anticipating a delta PPV between HPV RNA-Seq vs histology and cytology vs histology. For example, in the case where the ratio of HSIL to LSIL would tend to 1:2 as seen elsewhere [START_REF] Cuzick | Comparing the performance of six human papillomavirus tests in a screening population[END_REF] , the delta PPV could be up to +10.4 in favor of HPV RNA-Seq vs histology (range: +4.4 to +10.4, (Supplemental Figure 3 and Supplemental Table S4). This observation constitutes a solid S4. (JPEG 2780 Ko)


Noura Slimani
email: noura.slimani@edu.mnhn.fr

David Sánchez-Fernández

Eric Guilbert

Moncef Boumaïza

Simone Guareschi

Jean Thioulouse

Assessing potential surrogates of macroinvertebrate diversity in North-African Mediterranean aquatic ecosystems

Keywords:

The need to use surrogates of biodiversity is quite relevant in threatened habitats harboring high values of biodiversity, such as the Mediterranean aquatic ecosystems. In this study, we assess the performance of eight macroinvertebrate groups (Coleoptera, Heteroptera, Odonata, Trichoptera, Plecoptera, Ephemeroptera, Crustacea, and Mollusca) as surrogates of the whole aquatic macroinvertebrate assemblage in 49 localities from Northern Africa (Tunisia). Specifically, we aimed to test i) the congruence of the patterns of species richness and composition among these eight groups (at species level) in order to propose which groups could be accurate as indicators of diversity of the whole community, and ii) if higher-taxon levels (genera or families) are suitable for predicting overall species richness and composition in these ecosystems. In total, we found 72 families, 157 genera and 280 species. Our results show a high congruence between the patterns of species richness and composition of Ephemeroptera, Coleoptera (even at higher taxonomic levels, especially genus) and the whole community. Thus, we recommend the use of these two groups as surrogates of macroinvertebrate diversity in inland aquatic ecosystem in the study area. They can be used for both i) the rapid and inexpensive monitoring of biodiversity in aquatic ecosystems and ii) conservational studies in order to identify areas with the highest values of freshwater biodiversity in Mediterranean areas. Finally, high values of congruence among taxonomic levels were found suggesting that, in general, higher taxa can be used as biodiversity surrogates for cost-effective practical survey in Mediterranean aquatic ecosystems from Northern Africa.

Introduction

One of the main concerns for conservation biologists is trying to reduce the high rates of biodiversity loss due to human pressures [START_REF] Kerr | Effects of human activity on global extinction risk[END_REF]. However, most of the species on Earth have not been described so far and major gaps exist in knowledge about their distributions [START_REF] Brown | Biogeography[END_REF]. These taxonomic and biogeographical gaps, known as Linnean and Wallacean shortfalls respectively [START_REF] Lomolino | Conservation biogeography[END_REF], place serious limitations on the ability to conserve biodiversity in the face of the ongoing extinction crisis [START_REF] Hortal | Limitations of biodiversity databases: case study on seed-plant diversity in Tenerife, Canary Islands[END_REF]. Conservation biologists and environmental managers are indeed striving tofind suitable surrogates for mapping and predicting biodiversity as an effective way to overcome this limitations [START_REF] Humphries | Measuring biodiversity value for conservation[END_REF][START_REF] Caro | On the use of surrogate species in conservation biology[END_REF], especially in those countries with high biodiversity levels and scarce naturalist tradition. Biodiversity surrogates are groups of organisms with a sound taxonomy that have been well surveyed in a region, and whose patterns of assemblage structure (species richness, endemism, rarity, etc..) or composition, (seeCorte et al., 2017) are assumed to be indicative of similar patterns of unsurveyed taxa in the same region [START_REF] Pearson | Treatment techniques for adult female survivors of childhood sexual abuse[END_REF]. However, the validity of this assumption is rarely evaluated. In this context, researchers have traditionally used plants and/or vertebrates, especially birds, whereas arthropods have received less attention in this kind of conservation studies [START_REF] Posadas | Using phylogenetic diversity measures to set priorities in conservation: an example from southern South America[END_REF][START_REF] Cardoso | The seven impediments in invertebrate conservation and how to overcome them[END_REF], despite the fact that they represent around 95% of all known animal species [START_REF] Barnes | The Invertebrates: A Synthesis[END_REF][START_REF] Wilson | Half-earth: Our Planet's Fight for Life[END_REF]. An interesting approach to include invertebrate species in the biodiversity assessments could be to consider the number of higher taxonomic groups as a surrogate of the number of local species within the same clade (or other taxonomic groups). The advantage of this approach is that the number of families or genera can be documented more rapidly than the number of species of hyperdiverse groups like invertebrates [START_REF] Williams | Measuring more of biodiversity: can higher taxon richness predict wholesale species richness?[END_REF][START_REF] Caro | On the use of surrogate species in conservation biology[END_REF][START_REF] Baldi | Using higher taxa as surrogates of species richness: a study based on 3700 Coleoptera, Diptera, and Acari species in Central-Hungarian reserves[END_REF][START_REF] Villaseñor | Higher taxa as surrogates of plant biodiversity in a megadiverse country[END_REF]. The need to use surrogates of biodiversity to urgently identify areas of high biodiversity is especially relevant in aquatic ecosystem. These environments are highly diverse, since they constitute only the 0.01% of the world's water, equivalent to only 0.8% of the Earth's surface area, and support at least 100,000 species, i.e., approximately 6% of the estimated 1.8 million described species [START_REF] Dudgeon | Freshwater biodiversity: importance, threats, status and conservation challenges[END_REF][START_REF] Balian | The freshwater animal diversity assessment: an overview of the results[END_REF][START_REF] Heino | Climate change and freshwater biodiversity: detected patterns, future trends and adaptations in northern regions[END_REF]. Besides, freshwater ecosystems, especially those located in the Mediterranean basin, are subjected to a high human pressure (e.g., overexploitation; water pollution;flow modification; degradation of habitat; and invasion by exotic species: [START_REF] Dudgeon | Freshwater biodiversity: importance, threats, status and conservation challenges[END_REF][START_REF] Janssen | Exploring, exploiting and evolving diversity of aquatic ecosystem models: a community perspective[END_REF], and climate change [START_REF] Heino | Climate change and freshwater biodiversity: detected patterns, future trends and adaptations in northern regions[END_REF].

A number of studies have focused on the assessment of different taxonomic groups as biodiversity surrogates. Some of them have criticized the use offlagship species and other surrogate concepts in conservation, especially when the chosen areas are not consistently associated with high local biodiversity over space and time [START_REF] Roberge | Usefulness of the umbrella species concept as a conservation tool[END_REF]. In an exhaustive review on freshwater ecosystems, [START_REF] Heino | Are indicator groups and cross-taxon congruence useful for predicting biodiversity in aquatic ecosystems?[END_REF]found that indicator groups and, more generally, crosstaxon congruence do not appear to be particularly relevant for conservation in the freshwater realm. [START_REF] Guareschi | Cross-taxon congruence in wetlands: assessing the value of waterbirds as surrogates of macroinvertebrate biodiversity in Mediterranean Ramsar sites[END_REF]found a limited concordance between assemblage patterns of macroinvertebrates and water birds in Iberian wetlands. However, other studies point out to high congruence among groups, suggesting that water beetles in Mediterranean semi-arid regions [START_REF] Sánchez-Fernández | Are water beetles good indicators of biodiversity in Mediterranean aquatic ecosystems? The case of the Segura river basin (SE Spain)[END_REF] or snails in wetlands of Northeastern China [START_REF] Guan | Snails (Mollusca: Gastropoda) as potential surrogates of overall aquatic invertebrate assemblage in wetlands of Northeastern China[END_REF] could be useful as effective indicator of freshwater biodiversity. In this study, we assess for thefirst time in Africa the performance of eight macroinvertebrate groups as potential surrogates of aquatic (macroinvertebrate) biodiversity. Specifically, we aimed to test i) the congruence of the patterns of species richness and composition among these eight groups (at species level) in order to identify which groups can be used more accurately as biodiversity indicators in inland aquatic ecosystem in Northern Tunisia, and ii) if higher-taxon richness are suitable for predicting overall species richness and composition in these ecosystems.

Methods

Study area

The study was performed in northern Tunisia (Fig. 1), a Mediterranean region located in the Maghreb (North Africa), encompassing four watersheds (Medjerda; Northwestern; Ichkeul organized within two large sub-basins: Sejnane and Joumine; and northeastern watershed. Tunisia is located in a contact zone marking the transition from the temperate humid Mediterranean climate to the dry Saharan climate [START_REF] Zielhofer | Mid-and Late Holocenefluvial chronology of Tunisia[END_REF]. The northern area of Tunisia is characterized by a typical Mediterranean climate, ranging from humid (Mogods-Kroumirie region) to sub-humid (Bizerte region). A total of forty-nine non-impacted sites were selected to include a quantitative representation of the main aquatic habitats (seeArrignon, 1976), in the study area (Fig. 1and Table S1 in Supplementary Material): lotic freshwater (20 sites), lentic freshwater (7 sites), lotic saline (17 sites) andlentic saline waters (5 sites).

Biological data

For each one of the 49 selected sampling sites, faunistic compositions of eight taxonomic groups of macroinvertebrates (Crustacea, Mollusca and 6 groups of insects: Coleoptera, Heteroptera, Odonata, Trichoptera, Plecoptera, Ephemeroptera) were obtained from standardizedfieldwork. All localities were sampled in 2013, conducting the same sampling effort in each one of them. In each site, macroinvertebrates were sampled using both Surber nets (300 μm) and kick net (filet Troubleau) during twenty-five minutes across the entire habitat heterogeneity. Samples were preserved in 70% ethanol and taken to the laboratory for identification at species level. The applied methodology is common in studies on freshwater macroinvertebrates (e.g., [START_REF] Picazo | Water beetle biodiversity in Mediterranean standing waters: assemblage composition, environmental drivers and nestedness patterns[END_REF].

Statistical analysis

Richness patterns

Spearman rank correlations were used to evaluate the relationship among the species richness patterns of the different groups of macroinvertebrates and the Total Species Richness (TSR). To avoid giving higher weight in the correlation of the groups with a greater number of species, for each taxonomic group we also calculated the Remaining Richness value (RR), defined as the total number of species at a site (of all eight groups considered) minus the number of species belonging to the considered indicator group. Finally, we also conducted Spearman correlations to explore whether the higher taxon richness (family, genus) is correlated with the species richness (and with the RR values).

Community composition patterns

We used non-metric multidimensional scaling (NMDS,'metaMDS' function in vegan) to summarize composition patterns in the aquatic macroinvertebrates dataset. We repeated the NMDS using different taxonomic levels for the entire community (family, genus and species) and within each one of the eight invertebrate groups. The analyses were performed on a biological matrix based on presence-absence data and using the Jaccard index. As some taxonomic groups do not occur in all sites, a dummy taxon was added (when necessary), to allow the comparison among groups. A Procrustean analysis was applied to evaluate the degree and significance of community concordance among the NMDS ordinations of the different taxonomic levels and different taxonomic groups. Procrustean rotation analysis is regarded as a robust method for concordance analysis (Peres-Neto and Jackson, 2001) being frequently used to study aquatic communities (e.g., [START_REF] Virtanen | Community concordance between bryophyte and insect assemblages in boreal springs: a broad-scale study in isolated habitats[END_REF][START_REF] Guareschi | Cross-taxon congruence in wetlands: assessing the value of waterbirds as surrogates of macroinvertebrate biodiversity in Mediterranean Ramsar sites[END_REF][START_REF] Valente-Neto | Selecting indicators based on biodiversity surrogacy and environmental response in a riverine network: Bringing operationality to biomonitoring[END_REF]. Three dimensional NMDS ordinations were compared with the function'protest' (vegan package) where the Procrustean rotation analysis was accompanied by a permutation test (n = 9999). The statistic obtained is a Procrustes correlation r derived from the symmetric Procrustes residual m2(r =√ 1 -m2). Congruence analyses were performed in three parts: i) comparing each order at family-genus and species level with the entire community at species levels; ii) comparing inside each order the concordance among its different taxonomic levels (only where the Order is present); and iii) using the entire community and testing the concordance among taxonomic levels. With thefirst two parts we aim to understand which group may be surrogates of others composition patters, while the last step provides useful information about taxonomic resolution useful in aquatic biodiversity assessment. All statistical analyses were performed with R software (R Core Team, 2014), using ade4 and vegan package (http://pbil.univ-lyon1.fr/ ade4/, [START_REF] Thioulouse | Interactive multivariate data analysis in R with the ade and ade4TkGUIpackages[END_REF][START_REF] Oksanen | and environment for statistical computing[END_REF].

Results

Richness patterns

A total of 72 families, 157 genera and 280 species were found in the 49 localities studied (Table 1). The group with highest richness was Coleoptera (117 species) followed by Mollusca (40) and Heteroptera (36). Coleoptera and Heteroptera were the most widespread groups in the study area, being present in all sampling sites (Table 2). Ephemeroptera, Mollusca, Odonata and Crustacea were absent from one type of habitat (lentic or lotic saline waters). Trichoptera species were found in 10 sites, including two types of habitat (both lotic and lentic freshwaters), while Plecoptera species were confined to freshwater lotic ecosystems appearing only in two sites (Table 2). The richness patterns of all groups except Plecoptera were significantly correlated with the pattern of total richness (TSR) with Coleoptera showing the highest values, followed by Ephemeroptera (Table 3). In the same way, the species richness patterns of all groups but Plecoptera and Trichoptera were significantly correlated with their respective RR values, with Ephemeroptera showing the highest values (Table 3).

Among groups, the strongest correlation was found between Ephemeroptera and Mollusca, followed by Mollusca and Odonata (Table 3). The groups with highest number of significant correlations were Ephemeroptera (with all groups) and Trichoptera (with 4). However, Plecoptera were just correlated with one group, and Coleoptera with two groups (p < 0.05; seeTable 3). On the other hand, the richness of families and genera considering the whole community were significantly correlated with TSR (r = 0.845 and 0.977; p < 0.001, respectively). In the same way, the richness of families and genera of any taxonomic group, were significantly correlated with TSR, their respective species richness and RR values, showing in general high Spearman correlation coefficients, especially with their respective species richness values (seeTable 4).

Community composition: patterns and congruence

The ordination space of thefirst three axes of the NMDS including all groups presented afinal stress value of 0.15 (for both family, genus and species level). Focusing on each order separately, Heteroptera and Coleoptera families showed the highest levels of concordance with the entire macroinvertebrate community at species level (r = 0.55; p < 0.001;Table 5). On the other hand, Plecoptera showed the lowest level of concordance with no significant p-value. In the case of Coleoptera and Heteroptera the same results were always obtained with or without using dummy taxa. At genus level Coleoptera and Heteroptera presentedthe highest level of concordance with the entire community at species level (r = 0.69 and r = 0.61 respectively) while Coleoptera and Ephemeroptera obtained the best concordance at species level (r = 0.88 and r = 0.56, respectively; seeTable 5). Focusing inside each order, the analyses were performed using only the sites where the order appears and it was not possible to proceed with Plecoptera data (just two sites with presence). In these cases Crustacea at family level presented the highest community concordance with their respective species composition (r = 0.90; p < 0,001) followed by Heteroptera (r = 0.80; p < 0,001;Table 5). Overall, every order at genus level acts like a good surrogate of its community composition at species levels (minor value were obtained for Mollusca r = 0.67). Finally, according to protest analysis, the community variation was, in all cases, significantly concordant across different taxonomic levels. However, the concordance between genus and species community assemblages was higher (r = 0.95; p < 0.001) than the concordance between family and species community (r = 0.69; p < 0.001).

Discussion

This is thefirst study addressing the potential use of biodiversity surrogates in North-African inland aquatic ecosystems. One of the main strengths of this study is the wide taxonomic spectrum considered. Unlike other studies (e.g. [START_REF] Heino | Among-taxon congruence in biodiversity patterns: can stream insect diversity be predicted using single taxonomic groups? Can[END_REF][START_REF] Bilton | How wide to cast the net? Cross-taxon congruence of species richness, community similarity and indicator taxa in ponds[END_REF][START_REF] Guareschi | Cross-taxon congruence in wetlands: assessing the value of waterbirds as surrogates of macroinvertebrate biodiversity in Mediterranean Ramsar sites[END_REF], we consider here the main macroinvertebrate groups (excepting Diptera) in Mediterranean freshwater ecosystems, both in terms of species richness and abundance (seeTierno de [START_REF] Tierno De Figueroa | Freshwater biodiversity in the rivers of the Mediterranean Basin[END_REF]. The high correlation values found between the species richness of the whole community and the species richness of most of the considered taxonomic groups, especially Coleoptera (r = 0.847) and Ephemeroptera (r = 0.730) suggest that at least these groups can be used as effective surrogates of the whole macroinvertebrate community, being congruent with the results found inSánchez-Fernández et al. ( 2006) and [START_REF] Bilton | How wide to cast the net? Cross-taxon congruence of species richness, community similarity and indicator taxa in ponds[END_REF]. Coleoptera and, even if in a less way Ephemeroptera, have been already stressed as surrogate of aquatic macroinvertebrate richness in freshwater ecosystems belonging to the Spanish National Parks [START_REF] Guareschi | Aquatic macroinvertebrate biodiversity: patterns and surrogates in mountainous Spanish national parks[END_REF]. Nevertheless, examining the correlation between species richness patterns is only one of the possible ways to evaluate biodiversity indicators [START_REF] Kati | Testing the value of six taxonomic groups as biodiversity indicators at a local scale[END_REF], and these results only partially match with those obtained with the community composition assessment. Previous studies suggest that strong concordance between multiple organism groups should be indicated by r-values > 0.7 (e.g., [START_REF] Heino | Are indicator groups and cross-taxon congruence useful for predicting biodiversity in aquatic ecosystems?[END_REF]. In this case, just Coleoptera at species level showed a r-value higher than 0.7, being Coleoptera, Heteroptera and Ephemeroptera (at any taxonomic resolution) the taxonomic groups whose community composition showed the highest levels of concordance with the entire community at species level. Thus, it seems that Coleoptera at species level stands out in both richness and (even if in less minor) composition assessments. The use of Coleoptera as surrogate shows some important advantages, as it encompasses the gradient of environmental heterogeneity in semi-arid Mediterranean regions , i.e., they are present in all localities irrespective of the type of habitat, showing similar results in each type of habitat separately (see Table S1 in Suplementary Material). However, the correlation values when excluding beetles species from the complete pool (RR), although significant, falls from 0.847 to 0.399 (Table 3), suggesting that the high species richness of this group is biasing the obtained results with the whole community. In the same way, the results on congruence of community composition were conducted only considering the whole community and not separately with RR due to methodological problems to compare values when some species groups are missing. Besides, excepting Ephemeroptera, the congruence among richness patterns of the different groups was weak in most of the cases. This result was similar to those found byHeino et al. ( 2003)across 110 headwater streams in Finland. Patterns of Coleoptera species richness just correlated significantly with patterns of species richness of two groups (Heteroptera and Ephemeroptera). These issues call in question in some degree the performance of water beetles as unique surrogate of biodiversity in Mediterranean freshwater ecosystems (seeHeino, 2010), and suggest its use in combination with others like Ephemeroptera and Heteroptera as a combined effective strategy. On the other hand, the correlation values between family and genus levels with TSR were high, suggesting that, in general, the higher taxa (both genera and families) could be used as biodiversity surrogates for cost-effective practical survey. Besides, we observed that assemblage composition patterns among taxonomic levels are concordant and significant, with relevant values of significance (minor value of r = 0.69). In this context, the community composition even at family level presented a significant level of concordance with the community composition at species level. The r-values obtained between groups at family level and the whole community at species level (0.23-0.55) may be considered quite similar to those obtained byDolph et al. (2011) studying concordance between fish and macroinvertebrate 2015) studying concordance among waterbirds and macroinvertebrates in Mediterranean wetlands. However, the highest level obtained between assemblages at genus and species levels, seems to stress genus level as the ideal compromise between classification effort and gathered information. The identification of individuals at genus level, generally possible at their last larvae stage or with adults, can be conducted using exclusively morphological characters (see e.g. [START_REF] Tachet | Invertébrés d'Eau Douce (2nd corrected impression)[END_REF]. Nevertheless, among adult insects, only the majority of Coleoptera and Heteroptera are strictly aquatic and they generally represent the most frequently recorded taxonomic groups in lentic systems (e.g., [START_REF] Nicolet | The wetland plant and macroinvertebrate assemblages of temporary ponds in England and Wales[END_REF][START_REF] Guareschi | Cross-taxon congruence in wetlands: assessing the value of waterbirds as surrogates of macroinvertebrate biodiversity in Mediterranean Ramsar sites[END_REF]. [START_REF] Valente-Neto | Selecting indicators based on biodiversity surrogacy and environmental response in a riverine network: Bringing operationality to biomonitoring[END_REF]also recommend EPT taxa (Ephemeroptera, Plecoptera and Trichoptera) at genus level as the best indicator group for monitoring the effects of riparian vegetation loss on aquatic biodiversity in Neotropical streams. Thus, considering all results together, we propose the use of Ephemeroptera in combination with Coleoptera both at genus level as the best surrogate of biodiversity macroinvertebrates in inland aquatic ecosystems of Tunisia. Besides, it is important to note that any potential surrogate of biodiversity must cover a reasonably wide geographic range and occur in a broad range of habitat types [START_REF] Caro | On the use of surrogate species in conservation biology[END_REF], and this is the case of both Ephemeroptera and Coleoptera [START_REF] Tachet | Invertébrés d'Eau Douce (2nd corrected impression)[END_REF]. In sum, we recommend the use of higher taxa as surrogates for the rapid assessment of aquatic biodiversity only when accurate information on species level (rare, endemic or endangered species) is not available (most of the cases working with invertebrates, especially in Africa). This proposal of biodiversity surrogates can be used i) for the rapid and inexpensive monitoring of biodiversity in aquatic ecosystems of Mediterranean areas and ii) for conservational studies in order to identify areas with the highest values of biodiversity in which the aquatic environment has been or will be threatened [START_REF] Dudgeon | Freshwater biodiversity: importance, threats, status and conservation challenges[END_REF][START_REF] Vörösmarty | Global threats to human water security and river biodiversity[END_REF][START_REF] Chadwick | Silent streams[END_REF][START_REF] Tierno De Figueroa | Freshwater biodiversity in the rivers of the Mediterranean Basin[END_REF]. Indeed, thefivefirst sites ordered by the criteria of complementarity using the data of Coleoptera and Ephemeroptera at genus level, encompassed 172 species (62% of the total), which is interesting for conservational purposes and could be useful for future biomonitoring programmes and biodiversity research. Ourfindings provide a relevant contribution to the task of searching areas of high aquatic biodiversity in poorly sampled areas such as North Africa, as well as to explore congruencies among taxonomic groups in inland aquatic ecosystems from this area. However further studies in other areas of Tunisia and other North African states are recommended to improve and update the distribution mainly on water beetles and mayflies (also from the whole macroinvertebrate community in some localities) to complement and validate our main results. Thus, different stakeholders (e.g. conservationists, environmental consultants, natural resource managers, universities) may directly benefit from thesefindings establishing thefirst practical information on biodiversity surrogate selection in these aquatic systems.

Fig. 1 .

 1 Fig. 1.Study area. Red dots indicate sampling localities in North Tunisia.

Table 1

 1

	Number of families, genera and species of the eight taxonomic groups
	recorded
	in the study area.
	```
	Families Genera Species
	```
	```
	Coleoptera (Col) 12 51 117
	Heteroptera (Het) 13 19 36
	Ephemeroptera (Eph) 8 19 24
	Plecoptera (Ple) 6 10 16
	Odonata (Odo) 7 14 19
	Trichoptera (Tri) 6 9 18
	Mollusca (Mol) 13 27 40
	Crustacea (Cru) 7 8 10
	Total 72 157 280
	```

Table 2

 2 Number of sites with the presence of each group in the different types of habitat.

	```
	Lotic-freshwater
	(n = 20)
	```
	```
	Lentic-freshwater
	(n = 7)
	```
	```
	Lotic-saline
	(n = 17)
	```
	```
	Lentic-saline
	(n = 5)
	```
	```
	All sites

Acknowledgements

DS-F was supported by a post-doctoral contract funded by the Universidad de Castilla-La Mancha and the University of Murcia.

Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi.org/10.1016/j.ecolind.2019.01.017.