Monitoring the SARS-CoV-2 pandemic: screening algorithm with single nucleotide polymorphism detection for the rapid identification of established and emerging variants.

OBJECTIVES: To evaluate a testing algorithm for the rapid identification of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) variants that includes the use of PCR-based targeted single nucleotide polymorphism (SNP) detection assays preceded by a multiplex PCR sensitive to S-Gene Target Failure (SGTF).
METHODS: PCR SNP assays targeting SARS-CoV-2 S-gene mutations ΔH69-V70, L452R, E484K, N501Y, H655Y and P681R using melting curve analysis were performed on 567 samples in which SARS-CoV-2 viral RNA was detected by a multiplex PCR. Viral whole-genome sequencing (WGS) was performed to confirm the presence of SNPs and to identify the Pangolin lineage. Additionally, 1133 SARS-CoV-2 positive samples with SGTF were further assessed by WGS to determine the presence of ΔH69-V70.
RESULTS: The N501Y-specific assay (n = 567) had an overall percentage agreement (OPA) of 98.5%. The ΔH69-V70-specific (n = 178) and E484K-specific (n = 401) assays had OPA of 96.6% and 99.7%, respectively. Assessment of H655Y (n = 139) yielded a 100.0% concordance when applied in the proposed algorithm. The L452R-specific (n = 67) and P681R-specific (n = 62) assays had an OPA of 98.2% and 98.1%, respectively. The proposed algorithm identified six variants of concern/interest (VOC/VOI)-Alpha (n = 149), Beta (n = 65), Gamma (n = 86), Delta (n = 49), Eta (n = 6), Kappa (n = 6)-and 205 non-VOC/VOI strains-including the variants under monitoring B.1.214.2 (n = 43) and B.1.1.318 (n = 18) and Epsilon (n = 1). An excellent concordance was observed for the identification of all SARS-CoV-2 lineages evaluated.
CONCLUSIONS: We present a flexible testing algorithm for the rapid detection of current and emerging SARS-CoV-2 VOC/VOIs, which can be easily adapted based on the local endemicity of specific variants.

Introduction

The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) pandemic is rampant with emerging variants of the virus being reported globally [1,2]. So far, four variants have been classified by the World Health Organization as variants of concern (VOC). Early January 2021, Pango lineage B.1.1.7 (Alpha) became the predominant SARS-CoV-2 variant worldwide [2,3]. Rapidly thereafter B.1.351 (Beta) and P.1 (Gamma) and recently B.1.617.2 (Delta) have spread worldwide [4]. Others have been designated Variants Of Interest (VOI): lineages B.1.525 (Eta), B.1.526 (Iota), B.1.617.1 (Kappa) and C.27 (Lambda) [1,5]. The former VOIs B.1.427/B.1.429 (Epsilon), P.2 (Zeta) and P.3 (Theta) were reclassified because it was demonstrated that they no longer pose a major added risk compared with other circulating variants [5]. VOC/VOIs harbour mutation clusters that might impact virus transmissibility, infectivity, disease severity or immune escape (e.g. N501Y, E484K/Q, L452R), or were associated with diagnostic failure (ΔH69-V70) [[6], [7], [8], [9], [10]]. Whole-genome sequencing (WGS) is the reference standard for variant detection but has some limitations: WGS is time-consuming with a turn-around time of 4–7 days, limiting its impact on rapid healthcare responses [11]. Additionally, WGS is currently not universally available because it is an expensive and labour-intensive technique, and sequence failures are also frequent when viral loads are low [12]. The WGS capacity is limited, although major improvements have been made in high-throughput sequencing capacities and in reducing costs. As several VOC/VOIs can be discriminated based on the presence of key mutations, the European Centre for Disease Prevention and Control (ECDC) described several alternative strategies for rapid detection of VOCs/VOIs. This includes the use of screening single nucleotide polymorphism (SNP) assays by real-time RT-PCR melting curve analysis and the use of a multiplex RT-PCR detecting S-gene target failure (SGTF) as a surveillance strategy [3,11]. We propose a flexible testing algorithm based on the use of a multiplex RT-qPCR able to detect SGTF followed by multiple mutation-specific probe-based PCRs for rapid identification of SARS-CoV-2 variants. WGS was used as reference method to validate the algorithm.

Methods

Sample collection, extraction and multiplex RT-qPCR detection of SARS-CoV-2 viral RNA

Nasopharyngeal swabs (n = 1700) collected between January and May 2021 were stored immediately after collection in DNA/RNA shield lysis buffer (Zymo Research, Irvine, CA, USA) and heat-inactivated at 71 ± 2°C for 30 minutes. SARS-CoV-2 RNA extraction was performed using the DEXR-15-LM96 extraction Kit (Diagenode, Liege, Belgium) on a KingFisher™ Flex Purification System (Thermofisher Scientific, Waltham, MA, USA) with input and elution volumes of 350 μL and 50 μL, respectively. RT-qPCR detection of the ORF1ab -, N - and S-genes was performed using the TaqPath™ COVID-19 RT-PCR kit on a QuantStudio™ 5 (Thermofisher Scientific). Data analysis was performed using FastFinder Analysis software (v4.3.3, UgenTec, Hasselt, Belgium). SGTF was defined as the lack of signal in the channel detecting the S-gene whereas the N- and ORF1ab-genes were detected in their respective channels. In all, 567 SARS-CoV-2-positive samples with Ct values of the N -gene below 25 were further evaluated by target-specific PCR melting curve analysis and WGS, the latter in the context of the Belgian National Surveillance programme for SARS-CoV-2 variants. Additionally, 1133 positive samples with SGTF in the multiplex PCR were further assessed by WGS to identify the presence of the ΔH69–V70.

SARS-CoV-2 variant typing by target-specific PCR melting curve analysis

RNA extracts were stored at –20°C in the original extraction plates or freshly prepared when the remaining volume was insufficient. Using a one-step protocol, cDNA was prepared from the RNA extracts and amplified on a LightCycler® 480II (Roche Diagnostics, Rotkreuz, Switzerland) applying VirSNiP SARS-CoV-2 Spike N501Y, del H69/V70 (ΔH69-V70), E484K, H655Y, L452R and P681R probe kits (TIB Molbiol, Berlin, Germany) followed by PCR melting curve analysis using the LightCycler 480 SW1.5.1 analysis software. Five-microlitre RNA extracts were processed with 5 μL TaqMan™ Fast Virus 1-Step MasterMix (Thermofisher Scientific) and 0.5 μL specific reagent mix in a total PCR-volume of 20 μL (reverse transcription: 5 min at 55°C; denaturation: 5 min at 95°C and 45 cycles of 5 seconds at 95°C; 15 seconds at 60°C and 15 seconds at 72°C; probe hydrolysis was performed at 40°C, 75°C and 95°C). Melting temperatures (Tm) obtained for the wild-type and SNPs are displayed in Table 1 .

Table 1

Overview of melting temperatures for different PCR melting curve analysis assays obtained with the described protocol

	WILD-TYPE	Variant	N
N501Y	57.3 ± 0.5	62.5 ± 0.4	106
ΔH69–V70	58.5 ± 0.4	64.9 ± 0.5	65
E484K	54.4 ± 1.0	60.1 ± 0.4	89
H655Y	49.2 ± 0.4	56.1 ± 0.5	77
L452R	46.4 ± 0.3	58.2 ± 0.3	27
P681R	52.4 ± 0.4	63.1 ± 0.3	27

Abbreviations: n, number of samples; Variant, variant, mutation present.

Data are presented as mean ± standard deviation.

Viral whole-genome sequencing

Whole-genome sequencing on viral RNA was performed using Nextera XT and CovidSeq protocols (Illumina Inc, San Diego, CA, USA). Briefly, for the in-house protocol (Nextera XT), cDNA synthesis was performed using the LunaScript RT SuperMix Kit (New England Biolabs, Ipswich, MA, USA), followed by SARS-CoV-2 whole-genome amplification with multiplex PCR using ARTIC v3 (nCoV-2019) to generate amplicons with overlaps. Amplicons were purified using AmpureXP beads (Beckman Coulter, High Wycombe, UK) and quantified using the Qubit double-strand DNA (dsDNA) High Sensitivity kit (Life Technologies, Carlsbad, CA, USA). Nextera XT library preparation kits were sequenced on a MiSeq (v2, 500 cycles) (Illumina). The commercial Illumina COVIDSeq test was performed using paired-end sequencing on an Illumina NextSeq (v2, 150 cycles) according to the manufacturer's instructions followed by data analyses with Illumina DRAGEN COVIDSeq Test Pipeline software. Sequence quality was validated using NextStrain. Sequences were uploaded on GISAID (https://www.gisaid.org; strain names starting with hCoV-19/Belgium/UZA-UA) [2].

Ethics

After initial SARS-CoV-2 diagnostic testing, residual nasopharyngeal samples allocated for destruction were used for validation of the target-specific PCR melting curve analyses according to presumed consent defined by Belgian legislation.

Statistics

Positive per cent agreement (PPA) was defined as the frequency of SNP detection by PCR versus SNP detection by WGS. Negative per cent agreement (NPA) was defined similarly based on wild-type detection. Overall per cent agreement (OPA) was defined as the number of correctly identified wild-type/SNP by melting curve analysis versus WGS. The degree of agreement was evaluated by Cohen's κ performed by SPSS (v27; IBM, Armonk, NY, USA).

Results

Concordance between PCR and WGS

For the six evaluated targets, the PPA, NPA, OPA and κ-coefficients for PCR melting curve analysis versus WGS are displayed in Table 2 .

Table 2

Overview of SARS-CoV-2 strains evaluated by PCR melting curve analysis and whole genome sequencing

Target	Samples tested	PCR			WGS		NA (%)	PA (%)	Kappa (± SE)g	OPA (%)
		WT	Var.	Ab Tm	WT	Var.
N501Y	567	261a	302b	–	267	300	97.8	99.3	0.979 ± 0.009	98.5
ΔH69-V70	178	13	159	6	13	165	100.0	96.4	0.795 ± 0.081	96.6
E484K	401	204c	190	6d	205	190	99.5	100.0	0.995 ± 0.005	99.7
H655Y	139	63	76	–	63	76	100.0	100.0	1.000 ± 0.000	100.0
L452R	67	11	55e	–	11	56	100.0	98.2	0.948 ± 0.052	98.5
P681R	62	8	53f	–	8	54	100.0	98.1	0.932 ± 0.067	98.3

Abbreviations: Ab Tm, number of isolates with aberrant melting temperatures (Tm); NPA, negative per cent agreement; NI, not interpretable because of analysis failure; OPA, overall per cent agreement; PPA, positive per cent agreement; SE, standard error; SNP, single nucleotide polymorphism; Var, number of isolates with specific SNP defined by the kit insert (501Y, ΔH69-V70, 484K, 655Y, 452R or 681R); WGS, whole-genome sequencing; WT, number of isolates without the assessed mutation, considered as wild-type (N501, absence of ΔH69-V70, E484, H655, L452 or P681).

Four samples were NI. In all of them, the SNP was absent in WGS.

Two false positives (B.1.214.2).

One sample was NI in which the SNP was absent in WGS.

B.1.617.1 (n = 6): Aberrant Tm due to E484Q.

One sample was NI in whichL452R was observed by WGS.

One sample was NI in which P681R was observed by WGS.

p < 0.001 for all assays.

N501Y

The N501Y-specific PCR was performed on 567 RNA extracts and the N501Y mutation was detected in 302 strains assigned to four SARS-CoV-2 lineages: B.1.1.7 (n = 149), B.1.351 (n = 65), P.1 (n = 86) and B.1.214.2 (n = 2). All were confirmed by WGS except for two samples identified as B.1.214.2, which were false positives by PCR. The detection of N501 (wild-type) by PCR was confirmed by WGS in all other extracts (n = 261). Four RNA extracts failed to amplify in the N501Y assay and were not interpretable, sequencing revealed that none of these were VOC/VOI, resulting in an OPA of 98.5%.

SGTF and ΔH69-V70

Several lineages displayed SGTF, which has been linked to ΔH69-V70 [13]. In total, 1133 samples with SGTF were analysed by WGS. The presence of ΔH69-V70 was confirmed in 1131 samples; in two sequences the sequencing coverage was too poor for analysis.

The ΔH69-V70-specific PCR was applied to 178 viral RNA extracts (165 with and 13 without SGTF in the RT-qPCR). In 159 strains with SGTF, ΔH69-V70 was detected by PCR and confirmed by WGS. In six SGTF-positive strains, ΔH69-V70 could not be detected by PCR. All were allocated to the B.1.525 lineage and displayed an aberrant Tm (mean ± SD) at 60.6 ± 0.2°C (n = 6) whereas the wild-type had a Tm of 58.5 ± 0.4°C and ΔH69-V70 had a Tm of 64.9 ± 0.5°C (n = 65; Table 1). This resulted in an OPA of 96.6% between PCR and WGS.

E484K, H655Y, L452R and P681R

E484K-, H655Y-, L452R- and P681R-specific PCRs were applied according to the testing algorithm presented in Fig. 1 . E484 and 484K were detected in 204 and 190 RNA extracts, respectively. In one sample harbouring E484, PCR was not interpretable. In six samples, the E484K-specific assay displayed an aberrant Tm at 53.1 ± 0.5°C (versus 54.4 ± 1.0°C and 60.1 ± 0.4°C for E484 and 484K, respectively; n = 89; Table 1), and they were all found to be B.1.617.1. The OPA of the E484K-specific assay was 99.7%. The H655Y mutation was only evaluated when both N501Y and E484K were present and had a 100% concordance with the presence/absence of H655Y in WGS analysis. When both N501Y and E484K were absent in samples without SGTF, L452R and P681R were assessed. These assays had an OPA of 98.5% and 98.3%, respectively, as one sample (identified as B.1.617.1 by WGS) failed to amplify in both PCR assays.

Fig. 1

Testing algorithm for rapid variant detection. The S-gene target failure (SGTF) in a multiplex RT-qPCR assay, indicating the presence of ΔH69-V70 was used as a first discriminator in the testing strategy. Follow-up single nucleotide polymorphism (SNP)-specific PCR melting curve analyses were performed in parallel (blue) or in reflex (yellow) when specific mutations were present. SNPs are displayed in red, the corresponding wild-type in green. Aberrant melting temperatures (grey) in certain SNP-specific assays might be used to identify specific variants.

Rapid detection algorithm

Fig. 1 proposes the testing algorithm based on initial SARS-CoV-2 detection (Ct values < 25.0) with a multiplex RT-qPCR sensitive to SGTF followed by a panel of six probe-based PCR melting curve analysis kits. The order of applying follow-up PCRs depends on the initial results of the multiplex RT-qPCR with two main arms based on the presence/absence of SGTF. The decision to perform follow-up assays depends on the relative frequency of variants in the population. An overview of the observed lineages by WGS and the observed mutations by PCR is shown in Table 3 .

Table 3

Pango lineages identified by WGS and lineage-specific mutations detected by PCR melting curve analysis

WHO label	Pango lineage	Total evaluated	SGTF∗	501Y	ΔH69-V70	484K	655Y	452R	681R
Alpha	B.1.1.7	149	149	149	144a	NA	NA	NA	NA
Beta	B.1.351	65	0	65	NA	65	0	NA	NA
Gamma	P.1	86	0	86	NA	86	76b	NA	NA
Delta	B.1.617.2	49	0	0	NA	0	NA	49	48c
Eta	B.1.525	6	6	0	0d	6	NA	NA	NA
Epsilon	B.1.427/B.1.429	1	0	0	NA	0	NA	1	0
Kappa	B.1.617.1	6	0	0	NA	0e	0	5f	5f
–	B.1.214.2	43	0	2g	NA	0	NA	0	NA
–	B.1.221	37	0	0h	0	0	NA	NA‡	NA‡
–	B.1.1.318	18	0	0	0	18	NA	NA‡	NA‡
–	B.1.258 and sublin.†	13	13	0	13	0	NA	NA	NA
–	B.1.160	13	0	0	0	1	NA	NA	NA
	B.1.619	11	0	0	0	11	NA	NA	NA
–	Other	70	2	0i	2	3j	0	0	0

Abbreviations: NA, not assessed; NI, not interpretable because of analysis failure; SGTF, S-gene target failure; WHO, World Health Organization.

∗Only samples that were assessed by melting curve analysis were included in this table.

†Ten samples were allocated to the B.1.258 lineage, one to the B.1.258.11 sublineage and three to the B.1.258.17 sublineage.

‡These samples were evaluated based on an algorithm that did not yet include follow-up testing for L452R and P681R mutations as neither the Delta, nor the Epsilon variant was endemic in Belgium at the time of testing.

a–j Discordances: a In five samples ΔH69-V70 was not assessed, all were SGTF. b In ten samples H655Y was not assessed, 484K detected. c In one sample P681R was not assessed, 452R detected. d ΔH69-V70 was not detected because of an aberrant Tm. e An aberrant Tm was observed in the E484K assay (n = 6). f In one sample, L452R and P681R were NI. g Two samples were falsely positive for 501Y. h In three samples N501Y was NI. i In one sample of the B.1.1.777 lineage, N501Y was NI. j In one sample of the B.1.525 lineage, E484K was NI.

Performance of the testing algorithm

The B.1.1.7 lineage (Alpha) was assigned based on two key mutations: N501Y and ΔH69-V70, the latter causing SGTF. Combined detection of the N501Y and ΔH69-V70 mutations had a 100.0% concordance with WGS for detection of the B.1.1.7 lineage in 144 samples and all were SGTF.

The allocation of the B.1.351 lineage (Beta) and P.1 lineage (Gamma) based on the combination of N501Y-, E484K- and H655Y-specific assays resulted in a 100.0% concordance with WGS. In all samples of the Beta variant (n = 65), the N501Y and E484K mutations were observed in the absence of H655Y and in all strains of the Gamma variant (n = 86), the N501Y, E484K and H655Y mutations were detected.

Identification of the B.1.427/1.429 (Epsilon), B.1.617.1 (Kappa) and B.1.617.2 (Delta) lineages was based upon the assessment of the L452R, P681R and E484Q mutations in the absence of other mutations assessed. The concordance of SARS-CoV-2 strains harbouring L452R and P681R in the absence of N501Y, E484K or H655Y with the B.1.617 lineage (n = 54) was 98.1%. Six samples had an aberrant Tm in the E484K-specific assay and were identified as Kappa. All others (n = 48) were identified as Delta and had a Tm identical to the wild-type (E484). We identified the B.1.427/1.429 lineage in one sample that harboured only the L452R mutation in the absence of the other mutations assessed by PCR.

The Eta variant (B.1.525 lineage) was identified in six samples, displaying 484K on top of SGTF. All these samples displayed an aberrant Tm in the ΔH69-V70-specific assay, which resulted in discordance between the presence of SGTF and the ΔH69-V70-specific PCR. Nonetheless, the presence of ΔH69-V70 was confirmed by WGS.

Discussion

With emerging variants spreading quickly around the globe, rapid detection of VOC/VOIs is vital for monitoring and relevant for decisive management (contact tracing, prolonged isolation) of the SARS-CoV-2 pandemic. In the currently proposed testing algorithm (Fig. 1), a multiplex RT-qPCR that detects the ORF1ab -, N - and S-genes was used as first discriminator to guide follow-up testing based on the presence/absence of an SGTF. We confirmed the excellent concordance between SGTF and ΔH69-V70 that was previously reported [13]. This makes them interchangeable in the testing algorithm when SGTF-sensitive diagnostics are not available. Importantly, B.1.1.7 is not the only lineage harbouring the ΔH69-V70, responsible for SGTF (e.g. B.1.525, B.1.620 and B.1.258). However, the presence of the 501Y mutation in SGTF samples resulted in a 100% concordance rate for detecting the B.1.1.7 lineage.

Further assessment of SGTFs by melting curve analysis (Fig. 1) relies on the endemicity of other variants. Regions where the B.1.1.7 + E484K variant, the B.1.525 and/or B.1.620 lineage is circulating, might need to perform follow-up detection of the E484K mutation systematically after N501Y assessment to obtain sufficient discrimination (B.1.1.7 versus B.1.1.7 + E484K and B.1.525 versus B.1.620). Interestingly, the Eta variant displayed an aberrant Tm at 60.6 ± 0.2°C (n = 6) in the ΔH69-V70 assay, presumably due to the A67V mutation allocated in the probe-binding site. This might be a signature for this lineage (Fig. 1), although it needs to be confirmed in future studies.

In the absence of SGTF and the presence of 501Y and 484K mutations, the H655Y-specific assay proved to be an excellent discriminator between the Beta and Gamma variants (100% concordance for both). This supports the observations of Matic et al., who used a similar set of target-specific RT-PCRs to identify the Alpha, Beta and Gamma variants [14].

The recent spread of the Kappa and Epsilon variants also requires their rapid detection. Therefore, follow-up assessment of the L452R and P681R mutations can be performed in the absence of N501Y and E484K. The B.1.617 sublineages harbour both L452R and P681R whereas B.1.427/B.1.429 harbours only L452R. Based on local frequency of these variants, parallel or sequential (P681R when L452R present) evaluation can be performed. In the E484K assay, the E484Q mutation (e.g. Kappa variant) displayed an aberrant Tm (n = 6, 100%). When confirmed in larger populations, this finding can become helpful for rapid discrimination between the B.1.617 sublineages. Similarly, the Lambda variant might display an aberrant Tm in the L452R-specific assay as it harbours the L452Q mutation, but this was not assessed in the current manuscript.

Even though SARS-CoV-2 displays only one or two nucleotide changes/month/lineage, high incidence rates and selection pressure due to increasing population immunity result in rapid occurrence of new variants [9,15]. Therefore, tailor-made variant detection strategies to contain the SARS-CoV-2 pandemic are warranted. The proposed algorithm implements different mutation-specific PCR melting curve analyses, which allows an easily tailored, rapid detection strategy that can be performed daily in routine analysis. In a parallel set-up, identification is sped up, but more non-discriminating tests are performed. The conditional assessment of specific mutations reduces the number of non-discriminating PCRs and the risk for false positives but prolongs time to identification. Nonetheless, the latter is relative as the turn-around-time of melting curve analysis is approximately 2.5 hours per run. When sample load is high, high-throughput set-ups can be designed employing automated liquid handlers and 384-well PCR plates. Alternatively, fixed multiplex panels for mutation detection—including most of the assessed mutations—recently became available, enabling high-throughput testing [16]. However, the high flexibility of individual target-specific PCRs is a major advantage compared with multiplex panels in terms of quick adaptation guided by local epidemiology. The proposed algorithm uses SGTF as first discriminator, which could be considered a limitation but alternative strategies including the E484K/Q, L452R, P681R/H and H655Y or K417N/T mutations could be used to discriminate between variants without the requirement of the SGTF. The lack of discrimination between respectively the Iota and Zeta variants (former VOI) and the Beta and Theta variants (former VOI) is another limitation. The latter can be overcome by inclusion of a K417N/T-specific assay.

A rapid variant detection algorithm can be a valuable add-on to surveillance sequencing but WGS remains crucial in the monitoring of SARS-CoV-2 variants, identifying new/emerging variants and providing the required data for a tailored rapid variant detection strategy. The data obtained by surveillance sequencing during this study showed increasing occurrence of the Delta variant, leading to the incorporation of the P681R and L452R PCRs in the testing algorithm, demonstrating its flexibility. Importantly, only a fraction of SARS-CoV-2 positive samples can be sequenced because of limited sequencing capacity—especially in periods with high incidence rates. Furthermore, not all laboratories have access to the sequencing technology. SNP variant detection might be able to fill this gap because of higher PCR capacities. When central reference laboratories perform surveillance sequencing, satellite laboratories could monitor variants by a rapid detection strategy. Additionally, determining the signature of variants involved in outbreaks and breakthrough infections after vaccination by rapid variant detection, can reduce the required number of samples in a cluster to be sequenced, thereby reducing the cost and labour significantly.

In conclusion, the current data illustrate that rapid detection of variants by PCR melting curve analysis based on the proposed algorithm is a quick and robust add-on to WGS for mapping SARS-CoV-2 variants. Furthermore, the algorithm can be easily tailored, based on local epidemicity of specific variants, demonstrating its flexibility which is crucial in the constantly evolving pandemic.

Author contributions

VM and JC contributed to study design; JC, MLM, KL, BBX and CL contributed to assay validation and implementation; JC, MLM, KL and BBX contributed to data acquisition and analysis; and JM performed the data processing and statistical assessment. JM, JC and VM drafted the manuscript and all authors contributed to critical revision of the manuscript for important intellectual content and final approval.

Transparency declaration

No relevant conflicts of interest have to be disclosed by any of the authors. This manuscript is part of the output from RECOVER (Rapid European COVID-19 Emergency research Response), which has received funding from the EU Horizon 2020 research and innovation programme (grant agreement number 101003589).