Qualitative fingerprinting of psychoactive pharmaceuticals, illicit drugs, and related human metabolites in wastewater: A year-long study from Riga, Latvia.

The COVID-19 pandemic has become an unprecedented public health emergency causing immense societal and socio-economic consequences. Multiple studies have outlined that interventions to curb the spread of the virus are likely to have an effect on substance use patterns. In this study, we explored the presence of psychoactive pharmaceuticals, illicit drugs and related human metabolites in 24-h composite wastewater samples that were collected weekly in 2021 from the central WWTP of Riga, Latvia. The analysis was performed via suspect screening approach using three separate high-resolution mass spectrometry (HRMS) workflows, which relied on reversed-phase liquid chromatography (RPLC), hydrophilic interaction liquid chromatography (HILIC) and direct infusion HRMS. In total, 39 out of 149 substances were detected throughout the sampling period. These include pharmaceuticals (mainly antiepileptics, antidepressants and antipsychotics), illicit drugs (e.g., MDMA, MDEA, cocaine, etc.) and new psychoactive substances (alpha-PVP). The results were evaluated in relation to COVID-19 incidence rate and the severity of containment and closure policies. For some compounds we observed temporal changes that may be potentially linked to the state of the pandemic. For instance, higher detection rates were observed for several illicit drugs during periods, when restrictions on public events were relaxed. Meanwhile, some psychoactive pharmaceuticals and drugs used to treat upper respiratory tract infections displayed increased prevalence in weeks when the national COVID-19 incidence rates were higher. However, without baseline reference data from previous years, it is difficult to discern how much of the relationships seen are linked to pandemic progression and seasonal variability.

Introduction

The coronavirus disease 2019 (COVID-19) pandemic has brought about various changes in healthcare systems, economies and societies nearly everywhere in the world. To reduce the spread of the SARS-CoV-2 virus, most countries have introduced a multitude of interventions such as nationwide lockdowns, school closures, cancellation of public events, travel restrictions and curfews [1]. However, these policies have had far-reaching consequences for both individuals and societies. Studies have shown that the first wave of the COVID-19 pandemic has been associated with a greater decline in mental health in several populations [2], [3]. Current evidence also suggests that the circumstances of the pandemic have led to changes in sales, purchases and consumption of pharmaceuticals. For instance, a slight increase in psychotropic medication sales was observed in Rome, Italy, from March 2020 to February 2021 compared with the same period in the preceding year, whereby the change of sales seemed to vary according to the pandemic phases and related lockdowns [4]. Furthermore, there are indications that restrictions imposed to curb the pandemic may have affected illicit drug use and the illicit drug supply across the world. In accordance to the European Monitoring Centre for Drugs and Drug Addiction (EMCDDA) trendspotter briefing published in June 2020, an overall decline in drug use, or some forms of drug use, was detected in Europe during the first 3 months of the pandemic [5]. Yet, a more recent EMCDDA trendspotter study (April 2021) reports that, despite some reductions reported during the initial lockdown period, the use of most substances returned close to previous levels as social distancing measures were eased over the summer period [6]. In general, it appears that the most significant changes concerning pharmaceutical and illicit drug use were observed during the first wave of pandemic, which inflicted an unprecedented global impact on our daily lives and caused severe disruptions in supply chain networks. Nevertheless, stringent short-term restrictions (e.g., lockdowns and curfews) during later stages of the pandemic may have a temporal impact on substance use and drug consumption patterns.

Wastewater-based epidemiology (WBE) is an approach based on the surveillance of parent substances, metabolites, pollutants, pathogens or lifestyle biomarkers in untreated wastewater samples to obtain qualitative and/or quantitative data on different aspects of population within a given wastewater catchment [7]. It can provide complementary information by pinpointing quickly developing changes in WW composition. During the COVID-19 pandemic, WBE has been used for a wide range of applications to explore temporal trends of substance use, self-medication practises and other lifestyle aspects, which may be affected depending on the circumstances caused by the pandemic. In this respect, several major topics are often addressed: alcohol consumption [8], [9], [10], use and production of illicit drugs and new psychoactive substances (NPS) as well as human pharmaceutical consumption patterns [10], [11], [12], [13], [14], [15]. In addition, the WBE methodology provides insights beyond questions concerning a given catchment population. A large fraction of compounds analysed in WBE-based studies, especially active pharmaceutical ingredients (APIs), illicit drugs and pollution biomarkers, can be considered environmental contaminants. Therefore, the results of small molecule WBE studies can also provide some information on the discharge of pollutants, as the removal of these compounds in wastewater treatment plants (WWTPs) is often insufficient [16].

The most common analytical technique for small molecule WBE is based on solid-phase extraction (SPE) followed by reversed-phase liquid chromatography (RPLC) coupled with electrospray ionization mass spectrometry (MS) [17]. In terms of the stationary phase, several options may be considered: conventional reversed-phase C18, polar end-capped C18 phases and even highly polar hydrophilic interaction liquid chromatography (HILIC) phases, which may be used alone or in two-dimensional LC application [18], [19]. Similar methodological flexibility exists in the selection of MS platform. Tandem MS (MS/MS) systems are primarily used for targeted quantitative methodologies, whereas high-resolution MS (HRMS), such as quadrupole time-of-flight and Orbitrap MS technology, is applied to augment the methods with non-target and suspect screening capabilities. Given the increasing quality and scope of HRMS data libraries (e.g., MoNA, mzCloud, METLIN, EU-MassBank, etc.), availability of in-silico fragmentation and retention time prediction tools and presence of well-established identification criteria, suspect screening approach has become a reliable complementary technique, which can be used in WBE-based studies [16], [20]. This enables qualitative screening for a wide range of compounds without authentic standards, under assumption that the chosen sample preparation and instrumental analysis method is compatible with the suspect list. In addition, intra-analyte comparison of signal intensities can be useful for estimating possible quantitative changes within the sample batch.

In this study, we applied three separate HRMS methods (i.e., RPLC-HRMS, HILIC-HRMS and DI-HRMS) to perform a wide scope qualitative screening of psychoactive pharmaceuticals, illicit drugs and related human metabolites in 24-h WW samples during the COVID-19 pandemic in Riga, Latvia. The samples (n = 51) were collected from the central WWTP of Riga on a weekly basis from December 22, 2020 to December 14, 2021. There were three main goals of study: (i) thoroughly assess the wastewater composition in terms of what substances from the suspect list can be confidently identified, (ii) explore the identification capabilities of each of the applied HRMS workflows and (iii) investigate the observed detection rates and adjusted signal intensities of the identified substances in relation to the state of the pandemic (i.e., the 14-day notification rate of newly reported COVID-19 cases per 1,00,000 population and the severity of restrictions). This way, we tried to illuminate whether or not the observed variations of certain substances might be potentially linked with the progression of the pandemic, suggesting a change in consumption behavior of the population within the WWTP catchment area.

Materials and methods

Chemicals and reagents

All the solvents were of HPLC grade and were purchased from Honeywell Fluka™ (France) or Supelco (PA, USA). Ammonium formate (99 %) was obtained from ACROS Organics (NJ, USA). Formic acid (≥98 %) and hydrochloric acid (≥37 %) was supplied by Sigma-Aldrich (Finland) and Honeywell Fluka™ (Austria), respectively. Ammonia hydroxide (25 %) was purchased from Chempur (Poland). Deionised water was obtained via Milli-Q plus system from Millipore (MA, USA). All native analytical standards used in this study were obtained either from Sigma-Aldrich (MO, USA), Fluka (Switzerland), Dr. Ehrenstorfer (Germany) or Biopure™ Romer Labs (Austria). Isotopically labelled standards of salbutamol-d3 and mebendazole-d3 were purchased from Dr. Ehrenstorfer (Augsburg, Germany). Glass fibre filters (0.6 µm) were obtained from Advantec (Japan). Strata-X reversed phase (60 mg/3 mL) and Strata-X-CW weak cation exchange (200 mg/6 mL) SPE cartridges were supplied by Phenomenex (CA, USA), while Oasis HLB (60 mg/3 mL) cartridges by Waters (MA, USA).

Sample collection and treatment

Sampling campaign was carried out from December 2020 to December 2021 in the WWTP “Daugavgriva” of Riga (Latvia), which has a capacity of 200 000 m3/day and serves a population of 697 000 inhabitants. Every week a twenty-four-hour composite influent wastewater sample was collected from Monday morning until Tuesday morning using a refrigerated automatic sampler ASP-Station 2000 RPS20B (Endress+Hauser). The automatic sampler collected the wastewater 24 times (200 mL per sub-sample) over the period of 24 h. An aliquot (1.5 L) of the composite sample was transferred into polyethylene terephthalate bottle, immediately frozen and stored at − 20 °C prior to extraction and analysis. A total of 51 samples were collected.

Sample preparation

The sample preparation protocol was optimised in terms of SPE stationary phase, elution conditions and reconstitution medium. More detailed information about the optimisation steps is given in Section S1. Samples were thawed at room temperature, filtered through a 0.6 µm glass fiber filters and adjusted to pH 3.0 ± 0.5 with 37 % HCl. For the SPE procedure, we used manually packed two-layer columns that were prepared by adding 180 ± 5 mg of Strata-X-CW cation exchange resin/water slurry (1:2, w/w) on top of the Oasis HLB (60 mg/3 mL) cartridge and compressing the top layer with a polyethylene frit. Prior to extraction, the cartridges were conditioned with 5 mL of MeOH and 5 mL Milli-Q grade water (pH 3.0 ± 0.5). Next, 50 mL of filtered sample was loaded into the cartridge at a flow rate of around 5 mL/min. After drying of the SPE cartridge by applying vacuum for 1 h, the analytes were eluted with 9 mL of 5 % NH4OH in MeOH. The elutes were evaporated to dryness under a gentle stream of nitrogen at 35 °C and reconstituted in 100 μL of MeOH:water (1:1, v/v) mixture that contained isotopically labelled syringe standards (salbutamol-d3 and mebendazole-d3) at a concentration of 0.2 ng/uL. The final extract was split into two aliquots - 50 μL were transferred to 350 μL fused insert vial for LC-HRMS analysis, while 25 μL were transferred to 1.7 mL glass vial and diluted with 800 μL of 0.1 % formic acid in MeOH:water (4:1, v/v) for DI-HRMS analysis.

5-Hydroxyindoleacetic acid (5-HIAA) analysis was carried out in accordance to the method reported by Pugajeva et al. via dilute-and-shoot approach [19]. In brief, 2 mL of wastewater sample was diluted with 8 mL of Milli-Q grade deionized water and 10 μL of formic acid was added. Finally, the diluted samples were filtered through Ultrafree-MC centrifugal PVDF filters (0.22 µm, Merck, Germany) and analysed by LC-MS/MS.

Instrumental analysis

Both LC methods (reversed-phase C18 and HILIC) were carried out on an UltiMate 3000 (Dionex, Germany) HPLC system. The column compartment and autosampler temperatures were kept at 35 °C and 14 °C, respectively. Injection volume was set to 4 μL. Chromatographic separation for the HILIC method was achieved on a Kinetex® HILIC (150 × 3 mm, 2.6 µm, Phenomenex, USA) column. Mobile phase A consisted of 10 mM ammonium formate and 0.1 % (v/v) formic acid in MeCN:water (9:1, v/v) while mobile phase B was 10 mM ammonium formate and 0.1 % (v/v) formic acid in MeCN:water (1:9, v/v). Meanwhile, chromatographic separation for reversed-phase C18 method was carried out on a Kinetex® C18 (50 × 3 mm, 1.7 µm, Phenomenex, USA) column using water as mobile phase A and methanol (MeOH) as mobile phase B. Gradient elution program for both LC methods is presented in Table S1 in the Supplementary material. Optimisation of the LC method is discussed in Section S1.

Mass spectra were acquired using a Q-Exactive Orbitrap HRMS system (Thermo Scientific, Germany) equipped with a heated electrospray ionisation source (ESI) in positive ionisation mode. Full scan and data dependant MS/MS (ddMS2) modes were operated at a resolution of 140,000 and 17 500 FWHM (at m/z 200), respectively. ESI source conditions and the main instrumental parameters for both full-scan and ddMS2 modes are summarized at Table S2. Scanning range in full-scan mode was set to 100 – 900 m/z, while for ddMS2 mode an automatic range selection was enabled without a fixed first mass. The ddMS2 precursor selection was based on an inclusion list (only [M+H]+ species). Fragmentation was triggered using 10 ppm mass accuracy window and stepped normalized collision energy (NCE) was applied at 10 %, 25 % and 40 % to generate product-ion spectra. Dynamic exclusion was set to 10 s. Mass accuracy of the Orbitrap HRMS system was ensured via external calibration approach by introducing a mixture containing caffeine, MRFA, Ultramark 1621 and n-butyl amine prior every sequence. Data acquisition was performed on Thermo Scientific Xcalibur (v. 4.1.). Compound identification was conducted in Thermo Scientific TraceFinder™ (v. 4.1) and peak integration was done in Skyline 21.1. software [21].

Direct infusion HRMS analysis was performed on a 7 T SolariXR FT-ICR-MS system (Bruker Daltonics, Germany) equipped with an ESI source. The method was adapted from our previous study [22]. Instrumental parameters were kept the same as in the original study. Details about the method are presented in Section S2 in the Supplementary material. Data acquisition and post-acquisition data treatment was performed using ftmsControl 2.2.0 and Compass DataAnalysis 5.0 from Bruker Daltonics, respectively. Suspect screening workflow was carried out using the R 4.1.1 programming language in RStudio (v. 2021.09) integrated development environment.

Determination of 5-HIAA was achieved on UltiMate 3000 UHPLC system (Dionex, Germany) coupled to the TSQ Altis™ triple quadrupole MS/MS system (Thermo Scientific, Germany) in positive ESI mode. The method was directly adapted from Pugajeva et al. without modifications [19]. Additional information about the method is given in Section S3 in the Supplementary material.

Suspect database

The suspect database used in this study contained a total of 149 compounds, belonging to the following classes: APIs, illicit drugs, NPS and their respective metabolites. The list of APIs was based on the annual report “2019 Statistics on Medicines Consumption” by the State Agency of Medicines of the Republic of Latvia [23]. Only APIs that matched the Anatomical Therapeutic Chemical (ATC) classification "N - Nervous System" were included in the study. Illicit drugs and associated metabolites were selected from the following studies: Terzic et al., Hahn et al. and van Nuijs et al. that have investigated the occurrence of these compounds in wastewater from different locations in Europe [24], [25], [26]. Finally, selection of NPS class compounds was based on a highly comprehensive review article by Gent and Paul (2021), where over 40 different studies (2013–2020) from all over the world have been complied to summarize the occurrence rate of NPS in wastewater samples. Substances that were reported by at least two different literature sources, were included in the suspect database of this study [18].

Each suspect list entry was supplemented with structural information (e.g., the molecular formula, SMILES and InChIKey) and HRMS full-scan features (accurate mass for two most abundant [M+H]+ ions as well as their theoretical ratio). The latter were generated using the “enviPat” R package [27]. Next, “webchem” R package was used to retrieve computationally predicted octanol-water partition coefficients (XLogP) and topological polar surface area (TPSA) values from PubChem repository [28]. Experimental fragment spectra were obtained from MassBank of North America (MoNa) and the MassBank of Europe. Only ESI-HRMS entries (positive ionisation mode) were included for further evaluation. Each experimental MS2 entry was reprocessed with MetFrag (command line), an in-silico fragmentation algorithm [29], that allowed product ion annotation and m/z value reassignment (mass error threshold, ± 10 ppm; ion intensity threshold, 10 %). The resulting list of fragment ions was combined, grouped by precursor and the five most intense and frequently detected fragments were selected as the experimental MS2 “fingerprint”. Meanwhile, in-silico spectra were generated using Competitive Fragmentation Modeling-ID (CFM-ID) algorithm [30]. Finally, the “ReTip” R package developed by Bonini et al. (2020) was used for the retention time prediction for both reversed-phase LC and HILIC methods [31]. Computation of the model relied on retention time values obtained from the compound list that was used for method development. Training and testing datasets were split 70:30, respectively, and 50 random iterations (seed 1–50) were used to select the model that performed the best in terms root-mean-square error and coefficient of determination. The suspect database is available as a spreadsheet (the first sheet of the dataset S1 in the Supplementary material).

Method development relied on a different set of compounds, as only a limited number of authentic standards were available from the original list of suspect substances. Thus, more than 300 compounds from various classes (e.g., veterinary drugs, human pharmaceuticals, ergot alkaloids, pesticides and pyrrolizidine alkaloids) were considered for the preliminary method development compound list (further referred to as “the development list”). Each compound on the list was supplemented with computationally predicted XLogP, TPSA values as well as hydrogen bond donor and acceptor count. Compounds for which any of the parametric values fell outside the min-max range of the suspect list were discarded. Next, one thousand random datasets were created, each containing 80 different compounds from the development list. The dataset that displayed the least difference in terms of mean and median TPSA and XLogP values was selected as the final development list. The best match is depicted in Fig. S2 in the Supplementary material. Finally, the list was manually supplemented with compounds from the suspect database for which we had reference standards. Several ergot alkaloids and N-oxides were also included to increase the proportion of basic compounds. The final method development database is presented in Table S3.

Data processing and analysis

Two suspect screening workflows were used in this study. Tentative identification of suspects for DI-HRMS methodology was based on the protocol described by Perkons et al. [22]. In line with the findings of the study, DI-HRMS may be prone of false positive hits, if the unknown substance is composed only of C, H, N and/or O atoms. Therefore, the suspect database for DI-HRMS screening was narrowed from 149 to 57 compounds that all contained at least one halogen atom or phosphorus/sulphur. In brief, one full-scan spectrum was acquired by accumulating 32 individual scans. The spectrum was internally calibrated via lock mass (ammonia adducts of polyethylene glycol) and blank subtracted from the procedural blank. Resulting m/z values were exported as comma-separated values and processed in the RStudio environment (v. 2021.0.9) with R (v. 4.1.1). The full-scan HRMS features were extracted and matched against the suspect list using the following criteria: accurate mass error for the two most abundant ions ± 1.25 ppm and the relative error between the theoretical and experimental ion ratios ± 50 %. Next, three MS2 spectra were obtained (collision energy values: 5, 15, and 25 V) for each precursor that was found in full-scan data. The obtained product ions were then matched against the suspect database (mass accuracy ± 5 ppm) and, if at least one experimental fragment was found, the signal was considered to be tentatively identified.

The LC-HRMS data were processed using Thermo Scientific TraceFinder™ (v. 4.1) software. Extracted ion chromatograms were built with the precursor mass tolerance set at ± 5 ppm. Peak picking S/N threshold was set to 5, but peak area cut-off value and retention time window limits were disabled, except for compounds that had authentic standards. In this case, 30 s retention time window was used (predicted retention time indices were used in later steps as an additional verification criterion). Precursor identification criteria were set as follows: isotopic pattern fit: ≥ 80 %, allowed intensity deviation: 10 %, mass accuracy: ± 5 ppm. Product ion confirmation shared the same mass accuracy and at least two (out of five) unique fragments had to be detected for a compound to be considered tentatively identified.

According to the confidence levels proposed by Schymanski et al. [20], level 1 is assigned only when the candidate signal matched the reference standard in all criteria (i.e., MS, MS2 and retention time). It has been stated that level 2a (probable structure) can be assigned if spectrum-structure match between an unknown and experimental library data is unambiguous. Given that only two experimental product ions had to be matched to meet the MS2 criterion, our approach may not lead to a unanimous conclusion regarding the true identity of the unknown signal. Hence, all identified entries were manually checked against mzCloud database to ensure that correct peak assignment had been made. For entries that displayed a similarity score ≥ 70 %, level 2a was assigned. If the score fell between 50 % and 70 %, level 2b (diagnostic evidence) was given. When experimental data about the compound in question were absent from the mzCloud library, MoNa database was surveyed instead. The same criteria were applied (score≥700 - level 2a, 500 ≤score<700 – level 2b). Finally, level 3 (tentative candidate) was given to suspects that could not meet experimental MS2 criteria (or experimental data were unavailable), but satisfied precursor mass accuracy and isotopic pattern requirements, produced at least two matching in-silico fragment ions and displayed a retention time within the 95 % confidence window of the prediction model. For the DI-HRMS data, manual verification of MS2 spectra via mzCloud or MoNA libraries was not conducted. Level 3 was assigned to those substances that passed all full-scan and MS2 criteria. Suspects that failed to meet any of the DI-HRMS identification criteria were discarded without further evaluation. The suspect screening workflow is presented in Fig. 1.

Fig. 1

Scheme of the applied suspect screening workflow and confidence levels of identification (1, 2a, 2b and 3) using LC-HRMS and DI-HRMS approaches.

Quality control and quality assurance (QC/QA)

Procedural blank samples were analysed in each sample batch. Analyte-specific limits of detection were established from procedural blanks to reduce false positive detection rate that can stem from elevated background noise or laboratory contamination. Corresponding blank sample areas within the retention time window were multiplied by a factor of ten and the resulting value was set as the lowest acceptable area for each suspect. Only compounds exceeding this threshold were declared detected. Aliquots of six influent wastewater samples were fortified with a mixture containing amisulpride, azithromycin, carbamazepine, fluoxetine, gabapentin, levetiracetam, o-desmethyl venlafaxine, oxcarbazepine, paroxetine, sertraline, tramadol and venlafaxine at 1 μg/L before the SPE procedure and analysed as QC samples. The mean relative standard deviations of measured peak areas for QC substances were as follows: 30.6 % for HILIC-HRMS method, 28.5 % for DI-HRMS method and 35.4 % for RPLC-HRMS method. Instrumental drift in terms of sensitivity was adjusted by isotopically labelled syringe standards (i.e., salbutamol-d3, mebendazole-d3). Peak areas of each syringe standard were normalized across the sample batch by dividing peak area with mean peak area. Sample-wise means were then calculated from the two syringe standards and peak areas of unknowns were divided by this value. This allowed to compensate for some variability in instrumental response of the HRMS system between injections. In addition, the instrumental sequence was randomized to minimize order effects.

The method for the determination of 5-HIAA has been validated in our previous research study [19]. The linear range of the method was from 1 to 200 μg/L in wastewater matrix (R2 >0.99). The absolute recovery of 5-HIAA ranged from 63 % to 85 %, while the precision varied from 9 % to 10 % (relative standard deviation). Method limit of quantification was 1 µg/L.

Results and discussion

Identification results

In total, we identified 39 out of 149 suspects across all three methodologies in 24-hour composite samples from the influents of WWTP “Daugavgriva” of Riga (Latvia). Within the class of APIs, 29 substances (23 parent compounds and 6 metabolites) were detected, while 10 substances (7 parent compounds and 3 metabolites) were detected for the class of illicit drugs and NPS. HILIC-HRMS approach was able to identify a slightly higher number of suspects (n = 28) compared to RPLC-HRMS (n = 25). For HILIC-HRMS method, nine substances were confirmed with authentic standards and, thus, were eligible for the highest confidence level. The remaining 19 candidates were manually confirmed using mzCloud and MoNA spectral databases. Among them, 18 were successfully matched against the library MS2 spectra (≥ 50 % similarity score) and levels 2a (n = 12) or 2b (n = 6) were assigned depending on the spectral match quality. One candidate did not meet the library match criterion and was therefore downgraded to confidence level 3. Likewise, 25 substances were identified using RPLC-HRMS approach with the following confidence levels: level 1 (n = 13), level 2a (n = 9) and level 3 (n = 3). An overview of the identified compounds by RPLC-HRMS and HILIC-HRMS are presented in Table 1 and Table 2, respectively.

Table 1

List of suspect compounds found using RPLC-HRMS methodology from 24-hour composite samples from the influents of WWTP “Daugavgriva” of Riga (Latvia), with corresponding m/z values, accurate mass errors, observed retention times, MS2 fragmentation pattern match scores and confidence levels in accordance to Schymanski et al. (2014).

Compound Name	Formula	Expected m/z	Mass error, ppm	RT, mina	MoNA or mzCloud match score, %	Number of matched MS2 fragmentsb	Confidence level
Human pharmaceuticals
10,11-dihydroxycarbamazepine	C15H14N2O3	271.1077	-2.45	6.2 (0.4)	80.5	3/5 (1/5)	2a
10-hydroxycarbamazepine	C15H14N2O2	255.1128	-1.94	6.4 (0.1)	84.1	3/5 (2/5)	2a
Amantadine	C10H17N	152.1434	-1.76	6.6 (0.5)	86.4	2/5 (1/5)	2a
Amisulpride	C17H27N3O4S	370.1795	-2.05	5.4 (−0.4)	72.8	2/4 (2/5)	1
Azithromycin	C38H72N2O12	749.5158	-1.45	11.4 (2.4)	42.1	1/5 (0/5)	1
Carbamazepine	C15H12N2O	237.1022	-3.18	7.4 (−0.5)	91.7	1/5 (2/5)	1
Carbamazepine-10,11-epoxide	C15H12N2O2	253.0972	-2.28	6.5 (0)	75.7	2/5 (2/5)	1
Ephedrine	C10H15NO	166.1226	-0.67	4.5 (1.1)	91	4/5 (3/5)	2a
Fluoxetine	C17H18F3NO	310.1413	-3.06	9.5 (−0.6)	34.5	1/5 (1/5)	1
Gabapentin	C9H17NO2	172.1332	-3.28	2.8 (−0.5)	93.1	5/5 (3/5)	1
Lamotrigine	C9H7Cl2N5	256.0151	-1.60	6.8 (0.6)	65.1	0/5 (0/5)	3
Levetiracetam	C8H14N2O2	171.1128	-1.82	4 (0.9)	93.7	3/4 (2/5)	1
O-desmethyl venlafaxine	C16H25NO2	264.1958	-2.80	6.2 (−0.4)	62.4	3/5 (5/5)	1
Oxcarbazepine	C15H12N2O2	253.0972	-3.24	6.8 (0.1)	84.7	3/3 (4/5)	1
Paroxetine	C19H20FNO3	330.1500	-2.79	9.3 (0.5)	68.9	2/5 (1/5)	1
Phenibut	C10H13NO2	180.1019	-1.25	–	82.3	–	2a
Sertraline	C17H17Cl2N	306.0811	-0.78	9.6 (−0.2)	44.1	3/5 (2/5)	1
Sulpiride	C15H23N3O4S	342.1482	-2.78	4.1 (−1.7)	84	3/3 (2/5)	2a
Tramadol	C16H25NO2	264.1958	-2.58	6.4 (−0.4)	75.3	1/5 (0/5)	1
Venlafaxine	C17H27NO2	278.2115	-1.99	8.8 (0.8)	92.5	3/5 (2/5)	1
Illicit drugs
2-Phenethylamine	C8H11N	122.0964	-0.20	5 (1)	90.5	1/3 (4/5)	2a
Benzoylecgonine	C16H19NO4	290.1387	-2.24	5.3 (−3.1)	77.7	2/2 (2/5)	2a
Cocaine	C17H21NO4	304.1543	-2.02	6.7 (−1.6)	83.1	2/5 (2/5)	2a
MDEA	C12H17NO2	208.1332	-1.84	7.6 (2)	55.1	3/5 (2/5)	3
MDMA	C11H15NO2	194.1176	-1.68	5.9 (1)	48.5	4/5 (3/5)	3

The value in parenthesis refers to the difference between the observed and predicted retention time.

Number of fragment ion signals matched per precursor using experimental and predicted (in parentheses) MS2 features.

Table 2

List of suspect compounds found using HILIC-HRMS methodology from 24-hour composite samples from the influents of WWTP “Daugavgriva” of Riga (Latvia), with corresponding m/z values, accurate mass errors, observed retention times, MS2 fragmentation pattern match scores and confidence levels in accordance to Schymanski et al. (2014).

Compound Name	Formula	Expected m/z	Mass error, ppm	RT, mina	MoNA or mzCloud match score, %	Number of matched MS2 fragmentsb	Confidence level
Human pharmaceuticals
Amisulpride	C17H27N3O4S	370.1795	-1.97	6.0 (0.3)	82.5	3/4 (3/5)	1
Amitriptyline	C20H23N	278.1903	-0.52	3.6 (0.6)	60.5	2/5 (1/5)	2b
Azithromycin	C38H72N2O12	749.5158	1.07	9.5 (3.2)	77.1	2/5 (0/5)	1
Clozapine	C18H19ClN4	327.1371	-2.75	4.5 (0.2)	57.2	4/5 (1/5)	2b
Codeine	C18H21NO3	300.1594	-3.14	6.4 (1.3)	73.7	2/5 (0/5)	2a
Ephedrine	C10H15NO	166.1226	-1.22	4.8 (0.2)	92.4	4/5 (3/5)	2a
Gabapentin	C9H17NO2	172.1332	-1.33	6.8 (−0.1)	88.1	5/5 (4/5)	1
Levetiracetam	C8H14N2O2	171.1128	-3.43	2.0 (−1.3)	82.2	4/4 (3/5)	1
Melperone	C16H22FNO	264.1758	-0.38	3.7 (0.2)	84.6	2/3 (2/5)	2a
Mirtazapine	C17H19N3	266.1652	-2.33	5.4 (2.3)	68.4	2/5 (1/5)	2b
O-desmethyl venlafaxine	C16H25NO2	264.1958	-0.50	5.0 (0.5)	89.2	4/5 (2/5)	1
O-desmethyl tramadol	C15H23NO2	250.1802	-4.33	4.6 (1)	85.3	2/5 (2/5)	2a
Paroxetine	C19H20FNO3	330.1500	-3.99	3.3 (−2.4)	84.4	5/5 (3/5)	1
Pregabalin	C8H17NO2	160.1332	-1.41	6.4 (0)	88.1	5/5 (3/5)	2a
Quetiapine	C21H25N3O2S	384.1740	-1.13	4.3 (−0.6)	61.5	3/5 (0/5)	2b
Sertraline	C17H17Cl2N	306.0811	-4.07	3.4 (−0.2)	84.9	3/5 (2/5)	1
Sulpiride	C15H23N3O4S	342.1482	0.07	6.4 (1.8)	51.8	2/3 (2/5)	2b
Tramadol	C16H25NO2	264.1958	-2.23	4.2 (0.2)	89.3	1/5 (0/5)	1
Venlafaxine	C17H27NO2	278.2115	-3.30	4.7 (0.1)	90.9	3/5 (2/5)	1
Illicit drugs
2-Phenethylamine	C8H11N	122.0964	-0.77	4.7 (−0.8)	80.0	2/3 (3/5)	2a
Alpha-PVP	C15H21NO	232.1696	-4.35	4.0 (1.4)	63.0	2/2 (2/5)	2b
Amphetamine	C9H13N	136.1121	-2.39	4.3 (−1)	86.1	2/2 (2/5)	2a
Benzoylecgonine	C16H19NO4	290.1387	-1.71	6.3 (2.8)	87.2	2/2 (4/5)	2a
Cocaine	C17H21NO4	304.1543	-2.32	4.8 (1.1)	83.5	5/5 (3/5)	2a
Ecgonine methyl ester	C10H17NO3	200.1281	-1.18	7.3 (3.1)	74.8	2/2 (3/5)	2a
MDA	C10H13NO2	180.1019	-0.64	4.2 (1.3)	35.2	3/5 (3/5)	3
MDMA	C11H15NO2	194.1176	-1.76	4.5 (1.8)	88.0	4/5 (2/5)	2a
Methamphetamine	C10H15N	150.1277	-1.21	4.5 (0)	88.3	2/2 (3/5)	2a

The value in parenthesis refers to the difference between the observed and predicted retention time.

Number of fragment ion signals matched per precursor using experimental and predicted (in parentheses) MS2 features.

As anticipated, DI-HRMS approach yielded the lowest number of successful identifications (see Table S4 in the Supplementary material). Eight APIs were detected in samples, while azithromycin and fluoxetine were only identified in fortified QC samples. Paliperidone was the only compound that was detected by DI-HRMS but not by chromatographic methods. This is a rather unusual finding hinting a possible false positive identification because it is extremely rarely used in Latvia. In 2020, only 0.04 defined daily doses (DDD) per 1000 inhabitants were sold for this API. Yet, paliperidone, also known as 9-hydroxyrisperidone, is the primary active metabolite of the antipsychotic agent risperidone, which has substantially higher DDD per 1000 inhabitants (1.04). Besides, risperidone is one of the APIs manufactured by JSC "Grindeks", which is a pharmaceutical company located in Riga, Latvia. The factory operates a state-of-the-art WWTP but the treated effluent is discharged into the municipal sewerage system and therefore trace level detection of the manufactured APIs and their transformation products could not be ruled out. In general, lower number of positive identifications can be attributed to the fact that DI-HRMS methodology relied on a condensed database comprised solely of compounds containing at least one halogen atom or phosphorus/sulphur (n = 57). Of these 57 compounds, only 10 were detected in the samples using one of the chromatographic methods. Given that DI-HRMS was able to detect 8 unique compounds (and 2 from QCs), the overall performance of the method can be assessed as sufficient. Only melperone was detected via LC-HRMS but not by DI-HRMS. In total, four compounds were successfully identified by all three methodologies: amisulpride, paroxetine, sertraline and sulpiride. An identification example of supliride via RPLC-HRMS, HILIC-HRMS and DI-HRMS approaches is presented in Fig. 2.

Fig. 2

RPLC-HRMS and HILIC-HRMS chromatograms, DI-HRMS full-scan data and corresponding MS2 spectra for an antipsychotic drug sulpiride measured in wastewater. Diagnostic MS2 fragments from the mzCloud library (ID: 3452) are at m/z 214.0169, m/z 112.1121, m/z 98.0964 and m/z 58.0651.

During manual verification of RPLC-HRMS data for 3,4-methylenedioxyamphetamine (MDA), an intense peak was observed in several samples at 6.55 min that satisfied the monoisotopic mass and isotopic pattern criteria of MDA. However, the mzCloud database query returned an 82.3 % match for an entry (ID: 6077) belonging to phenibut, a nootropic agent used to treat anxiety and insomnia. Phenibut has not been listed in the annual report of API consumption statistics in Latvia [23] and therefore was not included in the initial suspect list that was derived from the annual report. Yet, it is available as a prescription medication in Eastern European countries, including Latvia [32]. Considering that phenibut matches the ATC classification "N - Nervous System" and displayed a high quality MS2 match, it was included in the study with a confidence level of 2a.

Only one API (lamotrigine) and two illicit drugs (MDA and 3,4-Methylenedioxy-N-ethylamphetamine, MDEA) were tentatively identified with a confidence level 3 and could not pass the library match criterion by any of the methods used. For instance, several diagnostic MS2 fragments (i.e., m/z 79.054, m/z 105.070, m/z 135.044 and m/z 163.075) were detected for a peak at RT 4.20 min via HILIC-HRMS that matched both the mzCloud record of MDA (ID: 3061) and literature data [33]. Yet, as seen in Fig. S4, the fragment spectra suffered from interferences (e.g., m/z 163.148, C12H19 + and m/z 134.096, C9H12N+), resulting in an unsatisfactory mzCloud match score (35 %). A co-eluting matrix component of high intensity was discovered at m/z 179.1540 (C11H20N2, −1.48 ppm), which was concurrently subjected to the same ddMS2 scan due the isolation window setting ( ± 1.5 m/z) and prevented the acquisition of clean MS2 spectra. A different issue was encountered when identifying MDEA. In various WW samples during the summer season a cluster of three peaks from 7.05 to 7.55 min were detected using RPLC-HRMS with identical mass accuracy and isotopic pattern features matching those of MDEA. Upon closer investigation of MS2 spectra (Fig. S5), fragment ion at m/z 163.075 was detected for the latest eluting peak at 7.57 min, whereas fragments at m/z 135.044 and m/z 105.070 were detected in all three of them. Querying the mzCloud database revealed that several close matches exist in terms of the similarity score. The following candidates yielded similarity scores ≥ 50 %: 3,4-methylenedioxy-N,N-dimethylamphetamine (MDDM, 62.7 %), 4-methoxy-N,N-dimethylcathinone (55.1 %), MDEA (54.6 %), MDA 2-amido analogue (MMDPPA, 51.5 %), MDMA methylene homologue (52.5 %) and 1,3-benzodioxolyl-N-methylbutanamine (MBDB, 50.7 %). The latter two could be ruled out since they normally yield a diagnostic fragment at m/z 177.091, which was absent from our measurements. When MoNA database was used, MDEA was the top-ranking candidate (55.1 %, ID: EQ360601) compared to alternatives retrieved from mzCloud survey. However, in accordance to the mzCloud data, MDDM (ID: 5150), methoxy-N,N-dimethylcathinone (ID: 5309) and MMDPPA (ID: 6329) seem to be practically indistinguishable from MDEA without the help of reference standards. Based on the global monitoring data reported by Gent and Paul (2021) [18], MDEA seems to be the most likely candidate but the analytical features could not provide a non-ambiguous answer. Hence, it was decided to lower the confidence level of MDEA. Finally, we did not assign level 2b to lamotrigine because none of the fragment ions were matched using the top 5 MS2 feature approach. Since we used this filter for promoting candidates to higher confidence levels, this ultimately led to level 3 despite acceptable mzCloud match (65.1 %).

Particular attention was paid to the identification of synthetic cathinone α-pyrrolidinopentiophenone (alpha-PVP) as it was detected in several WW samples via HILIC-HRMS. Only one peak at 3.95 min (m/z 232.1691, −2.2 ppm) produced the characteristic MS2 features at m/z 91.054, m/z 105.033, m/z 126.128 and m/z 161.096 (see Fig. S6) but the library search ranked alpha-PVP only in the third position with 63.0 % similarity score (ID: 4956). Two candidates ranked higher. 4'-Methyl-α-pyrrolidinobutiophenone (MPBP, ID: 2396) and 4-methyl-α-pyrrolidinobutiophenone (4-methyl PBP, ID: 5806) displayed similarity score of 64.3 % and 64.1 %, respectively. However, the mzCloud database records suggest that both MPBP and 4-methyl PBP yield diagnostic product ions at m/z 133.101 and m/z 112.112, while the presence of butylidenepyrrolidinium ion at m/z 126.128 is not anticipated. Considering that the former two fragments were absent and the prominent ion at m/z 126.128 was detected, it was decided to reject alternatives leaving alpha-PVP as the most plausible candidate (level 2b). We were unable to identify other synthetic cathinones above level 3. Several samples from early autumn of 2021 displayed a peak at 4.44 min that met precursor criteria for mephedrone and produced prominent fragment ions at m/z 119.085, m/z 145.088, m/z 160.112. However, the mzCloud similarity score against mephedrone record (or its positional isomers i.e., 2- and 3- methylmethcathinone) remained below 50 %. In addition, the procedural blank displayed an exceptionally high noise level within the 5-ppm isolation window of mephedrone. Although none of the diagnostic MS2 fragments were observed in the blank sample, confidence level assignment was not made to avoid potential misidentification. High blank levels were also observed for methadone, leading to its removal from the list of successfully identified suspects.

Our identification strategy relied heavily on experimental database matching. Nevertheless, for augmentation purposes we applied two complementary in silico tools for prediction of the expected retention times and generation of fragment ion features. The CFM-ID was used to create a computational reference MS2 library containing the five most frequently generated fragment ions for each suspect compound. Overall, only 14 % of all identified suspects did not contain any of the five CFM-ID generated fragment ions. Among the successfully identified suspects via RPLC-HRMS and HILIC-HRMS, the mean detection rate for in silico and library MS2 fragments was 41 % (109 out of 265) and 62 % (142 out of 228), respectively. For DI-HRMS approach the outcome was comparable - 38 % for in silico and 64 % for library.

The “ReTip” R package developed by Bonini et al. (2020) was used for retention time prediction [31]. Initially, three different models (Xgboost, RandomForest and BRNN) were investigated using RPLC data obtained from compounds included in the method development set (n = 91). The accuracy of each model was assessed over the 50 random train/test splits (train to test ratio: 70/30). Lowest root-mean-square error (RMSE) and 95 % confidence intervals (95 % CI) were found for Xgboost followed by RandomForest and BRNN algorithms (see Fig. S7). The training set that achieved the highest accuracy for the Xgboost model was selected for retention time prediction of suspects (RMSE: 0.73 min, 95 % CI: 1.04 min, depicted in Fig. S8). For the HILIC-HRMS approach other models were not tested and the Xgboost model was trained and tested using the same design. The accuracy of the trained model under HILIC conditions was lower compared to RPLC (RMSE: 1.27 min, 95 % CI: 2.02 min, depicted in Fig. S9). The suspect screening results indicated a similar tendency. Under reversed-phase conditions the RMSE for the identified compounds was 1.10 min and 28 % of the suspects exceeded the ± 1.04 min 95 % CI window. Under HILIC conditions RMSE was 1.29 min and the exceedance rate of 95 % CI was 16 %. The determination coefficient R² between experimental and predicted RT values for HILIC and RPLC was 0.52 and 0.70, respectively. The prediction accuracy of the HILIC model did not comply with the quality criteria established by Alexander et al. (2015) [34]. Yet, it should not have a significant impact on the identification results because predicted RT values were used to qualify only one suspect (MDA) under HILIC conditions.

Detection frequencies of the identified compounds

Many of the confirmed suspect APIs were almost ubiquitous in the analysed 24-h WW composite samples. Out of 24 APIs, 15 were found in all samples. Among these substances, antiepileptics (n = 5), antidepressants (n = 3) and antipsychotics (n = 2) were the most prevalent API classes. Similarly, all five identified API metabolites (three carbamazepine metabolites and O-desmethylated metabolites of venlafaxine and tramadol) were detected in all samples. Overall, the detection rates are in line with the consumption data in Latvia and dosage amounts established by WHO Collaborating Centre for Drug Statistics Methodology. To demonstrate this, we calculated the consumption estimate of each suspect API by multiplying the annual DDD per 1000 inhabitants per day in Latvia with the defined daily dose (mg, oral) and combined this value with the detection status. From Fig. S10 it can be seen that APIs with higher doses and consumption rates are also the ones most frequently detected.

Alpha-PVP was the only NPS identified in the samples. It was detected in 18 samples from April 2021 to October 2021. The highest monthly detection frequency was observed during the summer period, especially in June 2021, when traces of alpha-PVP were found in all samples. Traces of 2-phenethylamine (2-PEA), a compound that has been occasionally classified as NPS and reported in WW samples from Europe, were observed in the majority of samples [35], [36]. However, 2-PEA is naturally occurring in humans and from these data alone it is impossible to differentiate between an illicit intake and endogenous origin. MDMA, methamphetamine and MDA were detected in all samples. The latter may originate either from MDA use or from the metabolism of MDMA. MDEA (n = 20) and amphetamine (n = 7) were detected only in a limited number of samples and exhibited pronounced seasonality with higher detection rates in summer 2021. The occurrence of amphetamine can also be linked to the fact that it is one of the metabolites of methamphetamine. Finally, cocaine and its two primary urinary metabolites benzoylecgonine and ecgonine methyl ester were detected in 100 % and 96 % of the samples, respectively. The final detection rates of identified suspects are presented in Table 3. RPLC results for O-desmethyl venlafaxine and tramadol were removed from the table due to unacceptable chromatographic separation.

Table 3

Detection frequency of identified compounds across RPLC-HRMS, HILIC-HRMS and DI-HRMS approaches. The confidence levels are given in the parenthesis. The total number of samples was 51.

Name	RPLC-HRMS	HILIC-HRMS	DI-HRMS	Substance type (class)
Human pharmaceuticals
10,11-dihydroxycarbamazepine	51 (2a)		NI	metabolite of carbamazepine
10-hydroxycarbamazepine	51 (2a)		NI	metabolite of carbamazepine
Amantadine	51 (2a)		NI	parent (dopaminergic agents)
Amisulpride	51 (1)	51 (1)	49 (3)	parent (antipsychotics)
Amitriptyline		51 (2b)	NI	parent (antidepressants)
Azithromycin	42 (1)	51 (1)	QC only (3)	parent (antibiotics)
Carbamazepine	51 (1)		NI	parent (antiepileptics)
Carbamazepine-10,11-epoxide	51 (1)		NI	metabolite of carbamazepine
Clozapine		50 (3)	27 (3)	parent (antipsychotics)
Codeine		51 (2a)	NI	parent (opiates)
Ephedrine	51 (2a)	51 (2a)	NI	parent (sympathomimetics)
Fluoxetine	40 (1)		QC only (3)	parent (antidepressants)
Gabapentin	51 (1)	51 (1)	NI	parent (antiepileptics)
Lamotrigine	51 (3)		40 (3)	parent (antiepileptics)
Levetiracetam	19 (1)	6 (1)	NI	parent (antiepileptics)
Melperone		51 (2a)	ND	parent (antipsychotics)
Mirtazapine		51 (2b)	NI	parent (antidepressants)
O-desmethyl venlafaxine		51 (1)	NI	metabolite of venlafaxine
O-desmethyl tramadol		51 (2a)	NI	metabolite of tramadol
Oxcarbazepine	51 (1)		NI	parent (antiepileptics)
Paliperidone			7 (3)	parent and metabolite of risperidone (antipsychotics)
Paroxetine	22 (1)	12 (1)	1 (3)	parent (antidepressants)
Phenibut	44 (2a)		NI	parent (psychostimulants, agents used for ADHD and nootropics)
Pregabalin		49 (2a)	NI	parent (antiepileptics)
Quetiapine		50 (2b)	3 (3)	parent (antipsychotics)
Sertraline	47 (1)	29 (1)	3 (3)	parent (antidepressants)
Sulpiride	51 (2a)	51 (2b)	49 (3)	parent (antipsychotics)
Tramadol		51 (1)	NI	parent (opioids)
Venlafaxine	51 (1)	51 (1)	NI	parent (antidepressants)
Illicit drugs
2-Phenethylamine	19 (2a)	39 (2a)	NI	parent (phenethylamines)
alpha-PVP		18 (2b)	NI	parent (synthetic cathinones)
Amphetamine		7 (2a)	NI	parent (phenethylamines)
Benzoylecgonine	51 (2a)	51 (2a)	NI	metabolite of cocaine
Cocaine	51 (2a)	51 (2a)	NI	parent (classic illicit drugs)
Ecgonine methyl ester		49 (2a)	NI	metabolite of cocaine
MDA		51 (3)	NI	parent and metabolite of MDMA (phenethylamines)
MDEA	20 (3)		NI	parent (phenethylamines)
MDMA	51 (3)	51 (2a)	NI	parent (phenethylamines)
Methamphetamine	51 (2a)	51 (2a)	NI	parent (phenethylamines)

NI – Not included in the suspect list of DI-HRMS; QC only – detected and identified only in QC samples; ND – included in the suspect list of DI-HRMS but not detected in samples.

Potential temporal trends of the identified compounds

The study was limited to qualitative analysis and no stability assessment of the identified compounds was carried out within the framework of this study. However, the storage of the samples followed the recommendations of McCall et al. , which states that samples should be stored at − 20 °C without additional adjustments [37]. Sample sequence order was randomized and isotopically labelled syringe standards were used to reduce the impact of instrumental drift during the analysis. In addition, we quantitively determined 5-HIAA in each sample, which has been recognized as relatively stable endogenous biomarker for the estimation of the catchment population size [38].

The peak areas of the identified compounds were adjusted with respect to the syringe standard and 5-HIAA followed by compound-wise normalization across all samples via z-score approach (see Fig. S11). First, we explored correlation between different separation techniques for suspects that were identified by more than one methodology. The Pearson product-moment correlation between the HILIC and RPLC data was calculated using the normalized peak areas. Among the 17 compounds that could be detected by both methods, only one (sertraline) did not show a strong correlation, with R< 0.5. The average correlation coefficient between the two methods was 0.78 (see Table S5). Overall, the results show that there is a high level of agreement between the RPLC and HILIC methods. For DI-HRMS the mean correlation coefficient was 0.69. Only one substance displayed R< 0.5 (lamotrigine, R=0.44). The findings suggested that in most cases the DI-HRMS data were generally comparable to those obtained using chromatographic separation methods.

Next, we tested whether the normalized peak areas showed any association with the respective COVID-19 incidence rates. Data on 14-day notification rate of newly reported COVID-19 cases per 1,00,000 population were retrieved from the Centre for Disease Prevention and Control of Latvia. Prior to this analysis three samples were excluded. Sample 37 (October 7, 2021) and sample 41 (November 5, 2021) were discarded due to missing HIAA values (39]. Meanwhile, the positive correlation observed for ephedrine is likely associated with its diastereomer pseudoephedrine. It is found in several over-the-counter preparations in Latvia and used as nasal decongestant medication to relieve cold and allergy symptoms. It is impossible to distinguish ephedrine from pseudoephedrine using standard LC-MS without enantiomeric profiling. Considering that the 14-day notification rate of newly reported COVID-19 cases was lower during the summer, it remains unclear whether the positive correlations for ephedrine (pseudoephedrine) are directly linked with the pandemic or rather to the regular seasonality of upper respiratory tract infections which are known to have lower incidence during the summer months. This also applies to codeine, a prescription opioid which is used as a painkiller and cough suppressant. Meanwhile, moderate negative correlation was found for two antiepileptics levetiracetam (R=−0.42, p = 0.003) and gabapentin (R=−0.32, p = 0.027). This could be attributed to higher incidence of epileptic seizures during summer since they are more likely to occur on bright sunny days [40]. As seen in Fig. 3, the normalized peak areas of both APIs peak at July 2021. On the contrary, two antidepressants mirtazapine (R=0.35, p = 0.016) and amitriptyline (R=0.33, p = 0.021) as well as O-desmethyl venlafaxine (R=0.38, p = 0.008) showed moderate positive correlation suggesting higher consumption during the months when COVID-19 incidence rates were higher. In accordance to Montgomery et al. (2021) the average consumption of antidepressants increased during the first wave of the COVID-19 pandemic in the USA [13]. Yet, such trend was not observed in a study conducted in Innsbruck, Austria, in which WBE provided no evidence for increased use of nervous system drugs during the early weeks of lockdown [15]. Since we do not have a baseline data regarding the consumption of antidepressants prior to the pandemic, a direct link cannot be established.

Fig. 3

Normalized peak areas of selected APIs that were identified during the study from December 22, 2020 till December 12, 2021. The dashed line represents the 14-day notification rate of newly reported COVID-19 cases per 1,00,000 population in Latvia. Pearson product-moment correlation is depicted at the left upper corner of each facet. Asterisks indicate statistical significance (*p < 0.05, **p < 0.01 and ***p < 0.001). The area highlighted in red indicates period from October 21 to November 14 when curfew was imposed as one of the confinement measures to curb the pandemic.

For the class of illicit drugs, only cocaine (R=0.52) and its metabolite ecgonine methyl ester (R=0.50) showed strong positive correlation. However, these data have to be interpreted with caution because the other cocaine metabolite benzoylecgonine did not show any statistically significant relationship (R=0.23, p = 0.17). Cocaine and ecgonine methyl ester exhibit lower stability in WW compared to benzoylecgonine [41]. The two former are therefore more likely to be lost during sample storage, yielding lower concentrations in samples from the earlier months of the sampling campaign. Moderately negative correlation was displayed by MDEA (R=−0.45, p < 0.01). A sudden spike in terms of MDEA peak intensities was observed starting from June 2021 till September 2021. This matches the timeframe when restrictions regarding public events were partially lifted. To assess the impact of containment and closure policies we used the stringency index, which is one of the indicators reported by the Oxford COVID-19 Government Response Tracker [42]. It reflects the extent of government interventions to contain the spread of the virus (e.g., school closings, bans on public gatherings, curfews, etc.). As seen in Fig. 4, normalized peak areas of MDEA start to rise immediately after the stringency index drops as the restrictions imposed on public events and entertainment venues are progressively lifted. Less pronounced, but similar trend can also be seen for MDMA, amphetamine and alpha-PVP.

Fig. 4

Normalized peak areas of selected illicit drugs and NPS that were identified during the study from December 22, 2020 till December 14, 2021. The dashed line represents the stringency index used by the Oxford COVID-19 Government Response Tracker. The area highlighted in red indicates period from October 21 to November 14 when curfew was imposed as one of the confinement measures to curb the pandemic.

Conclusions

In this study, a HRMS-based suspect screening approach was applied to identify the presence of psychoactive pharmaceuticals, illicit drugs and related human metabolites in 24-h WW composite samples that were collected over the period of twelve months (December 22, 2020 to December 14, 2021) from the central WWTP of Riga, Latvia. The analytical workflow relied on three separate methods: RPLC-Q-Exactive-Orbitrap-MS, HILIC-Q-Exactive-Orbitrap-MS and DI-FT-ICR-MS.

A total of 39 out of 149 substances (26 %) were identified with confidence levels between 1 and 3. Both RPLC-HRMS and HILIC-HRMS showed comparable performance in terms of identification capability. Meanwhile, DI-HRMS exhibited inferior sensitivity and lower detection rates but was able to successfully identify all but one of the suspects that were confirmed by chromatographic methods. Overall, the most frequently detected APIs were antiepileptics followed by antidepressants and antipsychotics. 10 compounds were identified for the class of illicit drugs, including alpha-PVP which was the only NPS that was regularly detected during the sampling period.

In 2021, Latvia faced two waves of the COVID-19 pandemic. As a result, a series of containment measures were implemented to limit the spread of the virus, which may have influenced the behavioural patterns associated with the consumption of illicit drugs and medicinal products. In this study, we attempted to illuminate possible temporal trends by exploring the changes in the signal intensities and detection frequencies of the identified suspects. Strong correlation between the observed intensities and the 14-day notification rate of newly reported COVID-19 cases was observed for several APIs, including azithromycin (R=0.77), codeine (R=0.55), ephedrine (R=0.60) and phenibut (R=−0.59). To our knowledge, this is the first time phenibut has been detected in WW samples. Correlation analysis also suggested a potential consumption decline for two antiepileptics (levetiracetam and gabapentin) in months when COVID-19 incidence rates were lower. The opposite trend was noticed for three antidepressants (mirtazapine, amitriptyline and O-desmethyl venlafaxine). The correlation analysis did not show significant relationship between illicit drug results and COVID-19 incidence rates. However, the occurrence of some illicit drugs (MDEA, MDMA, alpha-PVP and amphetamine) seemed to be more pronounced in WW samples collected when the restrictions on public events and entertainment venues were partially lifted (i.e., during the summer period).

Overall, the results demonstrate that the detection rates and analytical responses of the identified compounds are not uniform across 2021 and that some of the variability may be linked to the status of the COVID-19 pandemic and the ensuing consequences. However, given the lack of reference data from previous periods and the qualitative nature of this methodology, the obtained data are subject to a high degree of uncertainty, which may have an impact on the results presented in Section 3.3. Thus, subsequent quantitative study and continuous monitoring is advised to verify the findings.

CRediT authorship contribution statement

Ingus Perkons: Methodology, Formal analysis, Data curation, Visualization, Writing – original draft. Laura Elina Ikkere: Methodology, Formal analysis, Data curation, Writing – review & editing. Veronika Sukajeva: Investigation, Data curation. Romans Neilands: Resources, Conceptualization. Kristina Kokina: Resources, Methodology. Iveta Pugajeva: Methodology, Supervision, Project administration, Funding acquisition.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.