Sequence identification and expression profile of seven Dermacentor marginatus glutathione S-transferase genes.

Dermacentor marginatus is a widespread tick species and a vector of many pathogens in Eurasia. Due to the medical importance of D. marginatus, control measures are needed for this tick species. Currently tick control approaches rely mostly on acaricide application, whereas wrong and irrational acaricide use may result in drug resistance and residue problems. Vaccination as an alternative approach for tick control has been proven to be effective towards some tick species. However, immunization against D. marginatus has not yet reached satisfactory protection. The effort of in silico based analysis could predict antigenicity and identify candidates for anti-tick vaccine development. We carried out an in silico analysis of D. marginatus glutathione S-transferases (DmGSTs) in order to identify blood-feeding induced GSTs as antigens that can be used in anti-tick vaccine development. Phylogenetic analysis, linear B-cell epitope prediction, homology modeling, and conformational B-cell epitope mapping on the GST models were performed to identify highly antigenic DmGSTs. Relative gene expressions of the seven GSTs were profiled through real-time quantitative PCR (RT-qPCR) to outline GSTs up-regulated during blood feeding. The phylogenetic analysis indicated that the seven GSTs belonged to four classes of GST, including one in epsilon-class, one in zeta-class, one in omega-class, and four in mu-class. Linear B-cell epitope prediction revealed mu-class GSTs share similar conserved antigenic regions. The conformational B-cell epitope mapped on the homology model of the GSTs displayed that GSTs of mu-class showed stronger antigenicity than that of other classes. RT-qPCR revealed DmGSTM1 and DmGSTM2 were positively related to blood feeding. In sum, the data suggest that DmGSTM1 and DmGSTM2 could be tested for potential anti-tick vaccine trials.

Introduction

Ticks are obligate hematophagous ectoparasites of great medical and veterinary importance (de la Fuente 2018). As one of the hard ticks adapted to colder regions, Dermacentor marginatus (Sulzer) is a widespread tick species over northern Eurasian (Zhang et al. 2019), which carries and transmits a variety of tick-borne pathogens to livestock and human, such as Babesia (Hilpertshauser et al. 2006), Theileria (Chisu et al. 2019), Anaplasma (Song et al. 2018), Rickettsia (Li et al. 2019), Coxiella (Bonnet et al. 2013), Brucella (Wang et al. 2018), Borrelia burgdorferi (Chrudimska et al. 2014) and Orthonairovirus (Liu et al. 2019). Due to the medical importance, D. marginatus—as many other tick species—has become an emerging threat to the livestock and public health, which needs control measures to mitigate.

As obligate blood-feeding ectoparasites, ticks need to process a huge volume of blood up to 100-fold their unfed body weight (Magdas et al. 2015). Blood digestion could cause oxidative stress referred to as reactive oxygen species (ROSs), as mentioned in the literature review (Torres-Rivera and Landa 2008). High ROS levels may result in the oxidation of nucleic acids, lipids, and proteins (Torres-Rivera and Landa 2008). Over time, hematophagous arthropods have honed their antioxidant systems to avoid intracellular oxidative stress by Cu, Zn, and Mn superoxide dismutases (SOD), catalase (CAT) and glutathione S-transferase (GST). The former two enzymes work in succession to directly decrease ROS levels by turning superoxide anions into oxygen and water, whereas GST catalyzes the conjugation between glutathione and other toxic molecules, and helps excreting toxic metabolites by increasing their hydrophilicity (Whiten et al. 2017).

GSTs are a widely distributed complex multifunctional phase II metabolite-detoxifying enzyme superfamily involved in cellular detoxification of endogenous and exogenous toxic metabolites (Torres-Rivera and Landa 2008). Due to the importance of GST in ticks, functional characterizations were conducted to study tick GSTs. A gene encoding mu-glass GST from Rhipicephalus (Boophillus) microplus (Canestrini) was isolated from salivary gland cDNA library, and the gene expression profile revealed that BmGST is transcribed in salivary glands and midguts of the adult females but not in the larval stage (Rosa de Lima et al. 2002). Recently, a recombinant GST of Haemaphysalis longicornis Neumann was tested on rabbit against two closely related species Rhipicephalus appendiculatus Neumann and Rhipicephalus sanguineus (Latreille) and cross protection was obtained against adult R. appendiculatus females (Sabadin et al. 2017). Our previous study using recombinant D. marginatus GST mu-1, a homolog of HlGST, as antigen to protect rabbits showed partial protection toward adult D. marginatus infestation (Huercha et al. 2020). A more recent study on a delta-class GST of Ixodes ricinus (L.) revealed that the IrGST is induced by heme in blood meal digestion (Perner et al. 2016). The same team also found that the IrGST is a heme-binding protein which sequesters free heme to ease the capacity pressure on hemosome in the digestive cell of the midgut (Perner et al. 2018).

The antibody produced by the vaccinated host animal were able to neutralize tick protective antigen and/or pathogens (de la Fuente et al. 2016). Identifying epitopes in tick protective antigens is of great interest in designing epitope-based vaccines (Sanchez-Trincado et al. 2017), as in silico prediction methods could dramatically reduce the time and effort spent in potential epitope mapping by experimental testing (Potocnakova et al. 2016). B-cell epitope prediction study on tick GSTs have revealed important epitopes that share cross-reactivity among multiple tick species which was demonstrated through a peptide-based serological experiment (Ndawula et al. 2020). B-cell epitope-based vaccines could induce a high specific immune response and induce longer immune protection (Ndawula et al. 2020). Similar in silico analysis on tick-specific aquaporin-1 was carried out in the hope of finding specific epitopes vital to tick survival (Ndekezi et al. 2019). The aim of the study was to obtain GST sequence information and identify highly antigenic GST expressed during blood feeding of female D. marginatus hoping to identify potential effective antigen which could be further tested for anti-tick vaccine trials.

Materials and methods

Rabbits and tick specimens

Two 4-month-old female New Zealand rabbits (ca. 2 kg each) used in the experiment were kept in the veterinary facility, subordinate to Xinjiang Agricultural University, in isolated cages under natural conditions and were provided with clean water and commercial feed. The ticks used in the experiment were kept and reared in College of Veterinary Medicine, Xinjiang Agricultural University, under controlled conditions at 24 ± 1 °C, 90 ± 5% RH and L14:D10 photoperiod. Briefly, ca. 600 larvae and 150 nymphs were fed on one rabbit separately in each ear-bag made of dense cloth and adhered to the rabbit’s ears. Adult ticks (20 female ticks and 30 male ticks in each ear-bag) were fed on the other rabbit using the above-mentioned approach. Egg samples were collected from a single egg clutch laid by an engorged female tick. Pooled specimens included eggs (10 µg), unfed larvae (n ≈ 50), fed larvae (n = 10), unfed nymphs (n = 20), fed nymphs (n = 3), whereas single specimens included adult ticks and organs collected from female adult ticks. The experimental procedures related to animals are conducted under supervision of the Animal Welfare Committee of Xinjiang Agricultural University (permission nr 2019018).

Total RNA extraction and cDNA preparation

Total RNA was extracted from engorged adult female ticks and eggs using RNAiso Plus trizol extraction (Takara, Japan) and reverse-transcribed to synthesize first strand cDNA templates with a Fastking RT cDNA synthesis kit (Tiangen Biotech, China) following the manufacturer’s instructions. By referring to transcriptome data (not published) and GST sequences deposited in NCBI (https://www.ncbi.nlm.nih.gov/), specific primers were designed for amplifying CDs of the GSTs. The PCR amplifications were performed according to the procedures recorded in Table S1. The PCR amplicons were separated by electrophoresis on 1.5% agarose gels. The obtained PCR product were segregated and purified using the E.Z.N.A. Gel Extraction kit (OMEGA, USA). Subsequently the purified sequence fragments were ligated into the pMD19-T cloning vector (Takara, Japan), and transformed into DH5α strain of Escherichia coli, of which the insert was sequenced from a recombinant plasmid purified from a single colony.

In silico analyses

Amino acid sequences deduced from nucleotide sequences produced in transcriptome data of D. marginatus were predicted using ORFfinder server of NCBI (https://www.ncbi.nlm.nih.gov/orffinder/). Multiple sequence alignments were performed using ClustalW (https://www.genome.jp/tools-bin/clustalw) algorithm and displayed through BoxShade server (https://embnet.vital-it.ch/software/BOX_form.html). Phylogenetic analyses were conducted using PhyML software v.3.0 with a configuration optimized by Prottest software v.3.4. The obtained phylogenetic tree was displayed through iTOL server (https://itol.embl.de/tree/). Linear B-cell epitopes were predicted by Bepipred Linear Epitope Prediction v.2.0 (http://tools.iedb.org/bcell/) with a threshold of 0.5. The structure-based predictions were performed using DiscoTope v.2.0 sever (https://services.healthtech.dtu.dk/) under default setting (−3.7). The predicted molecular weight (Mw) and isoelectric point (pI) were determined using the ExPASy server (http://web.expasy.org/peptide_mass/). Molecular models were constructed using SWISS-MODEL server (https://www.swissmodel.expasy.org/) (Waterhouse et al. 2018), and the predicted models were analyzed using EzMol server (http://www.sbg.bio.ic.ac.uk/~ezmol/) (Reynolds et al. 2018).

Real-time quantitative PCR analysis

The single-strand cDNA of D. marginatus were used in real-time quantitative PCR (qPCR) for GST gene expression profiles. The cDNA samples originated from different developmental stages of D. marginatus including egg, larva (unfed and engorged), nymph (unfed and engorged), and adults of both sex (unfed). The cDNA samples regarding blood feeding state were originated from 1- to 4-day fed, semi-engorged, and engorged female adult ticks. In addition, midguts, salivary glands, Malpighian tubes, and ovaries of engorged and semi-engorged adult female ticks were harvested for synthesis of cDNA. The transcription profile was conducted through qPCR using QuantiNova SYBR Green PCR Kit (Qiagen, Germany) with an Applied Biosystems 7500 Fast Real-Time PCR System (Applied Biosystems, USA). Specific primers for target genes and an internal control are shown in Table S1. The procedure running qPCR was at 95 °C for 2 min initial denaturation, followed by 40 cycles of a denaturation step at 95 °C for 15 s and an annealing/extension step at 60 °C for 35 s. Dermagentor marginatus elongation factor 1 alpha (ef-1α) gene was used as a reference to normalize the gene expression data as in tick spices the ef-1α is expressed comparatively stable in various life stages (Nijhof et al. 2009). The data were analyzed with Applied Biosystems 7500 system SDS software. All the templates used in qPCR were fourfold diluted using double distilled H2O and the qPCR was performed in triplicate.

Results

Cloning of the GST sequences

Full coding regions corresponding to D. marginatus GSTs were obtained by RT-qPCR at annealing temperatures and extension durations according to amplification conditions. The sequencing yielded full coding regions of the seven GST contigs ranging from 654 to 717 bp, and the deduced protein sequences were all aligned with potential tick GST homologs in Non-Redundant database of NCBI with identity of at least 60% or higher (Table S2).

Phylogenetic analysis and amino acid sequences of the GSTs

The sequenced DmGST genes demonstrated that these encoded proteins with amino acid residues from 216 to 238. To classify the GST homologs, sequences of Ixodes scapularis Say, Aedes aegypti L. and Drosophila melanogaster Meigen were aligned for GST classification (Reddy et al. 2011). The best amino acid substitution model for building the phylogenetic tree was LG + I + G using the Prottest software. Subsequently, the phylogenetic tree was built by PhyML, and phylogenetic analysis demonstrated the grouping of the seven GSTs in various clades, respectively (Fig. 1). Among the GSTs, four of them were grouped into mu-class as amino acid sequence identity revealed relatively high sequence similarity within the class. A zeta GST was found in D. marginatus and grouped in a branch with zeta-class GSTs of I. scapularis. One of the GST sequences was identified as an epsilon-class GST, and one GST was identified as an omega-class GST according to their sequence identity between classes (Strange et al. 2001).

Fig. 1

Unrooted phylogenetic tree showing the relationship between Dermacentor marginatus GST orthologs. A phylogeny using protein sequences was constructed for classification of DmGSTs. Sequences were aligned using the ClustalW server. Model selection was performed with Prottest. The tree was built by PhyML using Maximum likelihood (ML) method. GST sequences obtained in the study are marked with dots. Bootstrap value was set to 1000. (For interpretation of the references to colors in this figure, refer to the web version of this article.)

Sequence alignments were performed among and between GST classes (Fig. 2). Alignment of mu-class GSTs showed that they all shared conserved peptide sequences. Some of the conserved regions revealed identical predicted B-cell epitopes using BepiPred (Jespersen et al. 2017) server 2.0 prediction (Fig. 2a). Comparation between sequences among GST classed using DmGSTM1 as a template indicated that the predicted antigenic regions were largely independent, and only limited sequences shared conserved regions that display predicted antigenic activity (Fig. 2b).

Fig. 2

Amino acid sequence alignment of Dermacentor marginatus (Dm) GSTs. The identical residues are depicted in a black background, whereas similar residues are shown in gray shading. The sequence alignments were done using ClustalW, and Boxshade was used to obtain the shaded schematic representation. a Mu-class GST sequences are aligned, and highlighted shade in green represents predicted B-cell linear epitopes. b GST sequences are aligned among different classes (using DmGSTM1 to represent mu-class), and highlighted shade in violet represents predicted B-cell linear epitopes. The aligned sequences and GenBank accession numbers: DmGSTM1 (MN699679), DmGSTM2 (MN699680), DmGSTM3 (MN699681), DmGSTM4 (MT665976), DmGSTE1 (MT665977), DmGSTZ1 (MT665978), DmGSTO1 (MT665979). (For interpretation of the references to colors in this figure caption, refer to the web version of this article.)

GST homodimer modeling and B-cell epitope mapping

To obtain the GST homodimer models, protein structures of the GSTs were predicted in SWISS-MODEL server using existing GST 3D models as templates (Table S3). All templates used for predicting the GST structures have at least sequence coverage over 94%. All the predicted GST structures showed over 40% amino-acid sequence identities to the template models except for DmGSTE1 (34.95%). The DmGSTZ1 showed highest QMEAN score (−0.48) among the predicted homodimer models, whereas the lowest QMEAN score belongs to DmGSTM4 as −2.51 indicating the predicted model are good for mapping the structural epitope of the seven GSTs.

The predicted conformational B-cell epitopes from each GST were mapped and highlighted on the 3D protein structure of homodimer GST models (Fig. 3). The GST models are fixed at the same template to present the same topology and are displayed with one symmetric side of the homodimer facing forwards and rotate 180° right to reveal the other symmetric side. All the GST models exhibited structural epitopes. Generally, the mu-class GSTs displayed larger areas of epitopes than other three groups of GSTs referred in the study, and the DmGSTM3 had the most epitopes predicted compared to the other GSTs. The DmGSTE1 was predicted with few epitopes (Fig. 3). Parallel comparation of predicted linear epitopes and structural epitopes on each GST revealed a same trend that the two predictions were more likely to overlap in the front half of the protein sequence in which the glutathione binding site of the GST locates (Table S4).

Fig. 3

Structural illustration of the predicted conserved B-cell epitopes. The structures for GSTs were built based on existing 3D GST models (PDB: 1b4p. 1 for mu-class, 3mak.1 for epsilon-class, 1fw1.1 for zeta-class, 4is0.1 for omega-class). The homodimer models of the seven GSTs are shown in two colors: blue represents chain A and beige represents chain B, whereas red stands for predicted conformational epitopes. Epitopes were mapped according to the prediction results by DiscoTope (Table S4). The GST models are fixed at the same template to exhibit one symmetric side of the homodimer forwards and the other symmetric side by rotating 180° right. M1–M4 represents DmGSTM1–4; E1, DmGSTE1; Z1, DmGSTZ1; O1, DmGSTO1. (For interpretation of the references to colors in this figure caption, refer to the web version of this article.)

GST gene expression profiles

The expression of the seven GST genes were measured by real-time RT-PCR. The GST genes were found expressed in most of the developmental stages and organs except that DmGSTM4 was not detected in unfed larva (Fig. 4). The DmGstm1 gene was expressed during all developmental stages. Blood feeding increased expression of DmGstm1 in female adult tick at 96 h after attachment and semi-engorged state over 2500-fold compared with its unfed state, however the expression of the GST was relatively low in eggs (Fig. 4a). The expression of DmGstm2 increased after blood feeding and exhibited up-regulation in most of the developmental stage and organs specially in midgut both in semi-engorged and engorged state. However, the expression was relatively low in salivary glands of female adults (Fig. 4b). The DmGstm3 showed down-regulation during the first 4 days after attachment on the host, known as slow feeding period (Franta et al. 2010), but, expression increased significantly on semi-engorged females and in most of their organs except midgut. Relatively high expression levels were observed in larval and nymphal forms of the tick as well. (Fig. 4c). The expression of DmGstm4 was up-regulated in malpighian tube and ovary of semi-engorged adult female and in larval stages over tenfold (Fig. 4d). The DmGste1 exhibited significant up-regulation in semi-engorged female ticks, larvae, and nymphs (Fig. 4e). The DmGstz1 showed relatively high expressions in malpighian tube and ovary in semi-engorged female ticks (Fig. 4f), and the expression profile of the DmGsto1 was characterized in relatively low expression in male adults and high expression in ovary (Fig. 4g).

Fig. 4

Differentially expressed GST genes in Dermacentor marginatus during blood feeding and in developmental stages. The X axis shows various states of the tick. The Y axis shows gene expression fold-changes. Unfed female was used as the reference. Dermacentor marginatus elongation factor 1 alpha was used to normalize the expression profiles. adult female tick; adult male tick; SE semi-engorged; E engorged; SG salivary gland; MG midgut; OV ovary; Mal Malpighian tube

Discussion

Vaccinating animals against ticks is considered one of the most promising tick control alternative approaches to reduce the use of acaricides, but no antigen has so far exceeded the efficacy of Bm86 vaccine against R. microplus (de la Fuente et al. 2016). As mentioned in the literature review that a number of tick protective antigens have been found to be useful in tick control (de la Fuente 2018). The tick GSTs were also investigated as potential candidates for vaccination against many tick species, and reached a protective effect by reducing repletion rate and egg production of female adult ticks, and egg hatchability (Huercha et al. 2020; Ndawula et al. 2019; Oldiges et al. 2016; Sabadin et al. 2017).

The GSTs are a superfamily present in eukaryotes and prokaryotes that are generally involved in the detoxification of endogenous and xenobiotic compounds (Ketterman et al. 2011). As phase II detoxification enzymes, GSTs play a very important part in tick survival in on-host and off-host periods. Feeding large amounts of host blood and the digestion of the blood meal could greatly increase redox pressure in ticks (Graca-Souza et al. 2006). Genomic analysis revealed 35 GSTs of epsilon, delta, mu, omega, zeta, kappa, and microsomal GSTs in I. scapularis (Reddy et al. 2011). In our work, there were seven GST genes with complete coding regions that encode GSTs which belong to four classes. The sequence alignments and phylogenetic analysis indicated these D. marginatus GST genes were homologs of other tick GSTs. All the obtained GST sequences were named according to the suggestion of nomenclature for insect GSTs (Ketterman et al. 2011).

The vaccine approach to control ticks is based on the fact that an antibody from the host animal could traverse midgut cells into hemolymph and reach other organs to carry out its biological function (Ackerman et al. 1981). Prior studies have noted the importance of B-cell epitope prediction, and revealed that tick GSTs share conserved epitopes in different tick species (Ndawula et al. 2020). In the current study, we found there were linear B-cell epitopes in different D. marginatus GSTs shared similar prediction patterns suggesting that they might have cross immunity within GST classes. Conformational epitopes are formed by residues apart in the protein sequence that are orderly folded in space together (Jespersen et al. 2017). As expected, conformation epitope prediction did not display similar mapping patterns among GSTs. Moreover, the conformation epitope prediction results indicated that mu-class GST may generate stronger B-cell immune response than other three classes of GSTs in this study. A previous study of a GST-cocktail vaccine against five species of ticks has obtained 35% efficacy, suggesting that the mu-class GST is a potential target of broad-spectrum antigens against ticks (Ndawula et al. 2019).

The RT-qPCR was designed to determine the expression profile of the obtained seven GST sequences. The current study found that all the GSTs were expressed during various tick-developmental stages. Most of the genes were significantly up-regulated after blood feeding. This result may be explained by the fact that digestion of a blood meal in tick midguts caused a surge in redox pressure, as large amounts of heme and bioavailable iron are released and need to be sequestered and transferred out of the cells to maintain homeostasis during endocytosis and degradation of hemoglobin in addition to the host transferrin in the midgut epithelial cells (Galay et al. 2015; Perner et al. 2016). Thus, the GSTs are up-regulated to counteract the increased redox pressure, whereas high level of redox could be beneficial for limiting the proliferation of certain pathogenic microorganisms in tick species (Budachetri et al. 2018).

Genomic studies of I. scapularis revealed 14 out of 32 cytosolic GSTs were from mu-class (Reddy et al. 2011). Mu-class GSTs were mainly involved in reducing intracellular redox pressure during blood feeding of tick species (Hernandez et al. 2018). In this study, we found that the up-regulated expression of these mu-class GSTs were positively correlated with blood feeding and digestion. Prior studies that have noted the importance of mu-class GST, and characterization of several GSTs were conducted (Hernandez et al. 2018; Oldiges et al. 2016; Shahein et al. 2008). The current study revealed that DmGstm1 and DmGstm2 were highly up-regulated after blood feeding in female adult ticks. These GSTs are likely to be related to detoxification of toxic metabolites produced in the processes of blood meal digestion and embryogenesis, as the GSTs expression in midgut and ovary increased significantly over time.

Besides mu-class GST, epsilon, delta, omega, zeta, kappa, microsomal GST comprised the rest of the tick GST superfamily (Reddy et al. 2011). The epsilon-class GST that was overexpressed at the presence of insecticide was first found in Anopheles gambiae Giles (Ranson et al. 2001). In the current study, the expression profile indicated that the DmGste1 was up-regulated in most of the life stages specifically in the Malpighian tube and ovary. The role of epsilon-class GST in D. marginatus may display the same function as found in mosquitoes through detoxification of insecticides (Saruta et al. 2019; Tchouakui et al. 2018). As mentioned in the literature review, the zeta-class GSTs have a widespread distribution in plants, fungi, and mammals. The GSTZ in animals and plants may play a role in the endobiotic metabolism (Ketterman et al. 2011). The alignment of the DmGSTZ1 showed high sequence identity with maleylacetoacetate isomerase, which corroborates earlier findings that maleylacetoacetate isomerase and GSTZ1 are identical (Blackburn et al. 1998; Fernández-Cañón et al. 2002). Omega-class GST that is a most recently defined GST class appears to be widespread in nature and have been identified in bacteria, insects, yeast, mammals, and plants. The GSTOs were characterized in insects, and are overexpressed under stress response and at the presence of insecticides (Balakrishnan et al. 2018; Yamamoto et al. 2011).

A prior study has noted the importance of a GST that showed heme-binding ability (Perner et al. 2018). This GST was able to reduce the pressure caused by extra free heme in the tick digestive cell by binding to heme at the cost of losing its ability of detoxification. The same team first found that the GST was the only GST gene which significantly up-regulated in response of heme presence by transcriptome analysis of I. ricinus (Perner et al. 2016). However, recombinant IrGST1 did not reveal anti-tick effect (Perner et al. 2018), although it might have a potential to be used as a vaccine candidate to target tick heme metabolic pathway by co-immunization with other tick protective antigens. A vaccine study on a hook worm, Necator americanus, using one of the heme-binding GSTs as antigen to immunization elicited 32 and 39% reductions in adult hookworm burdens in dog and hamster vaccine trials, respectively (Zhan et al. 2010). Nevertheless, there remains a question that some tick GSTs showed anti-tick effect, but the metabolic process and the endogenous substrate of the GSTs were not clear yet. Limitations of this study are that multiple sequence alignment of homologous genes between tick species was not performed, which is important for revealing cross-reaction epitopes between GSTs of different tick species. As further investigation on tick GST increases, further functional characterization would deepen our understanding of tick GSTs crucial for maintaining homeostasis that detoxify toxic endogenous metabolic substrates. In addition, the potential cross antigenicity of the DmGSTs needs to be demonstrated experimentally.

Conclusions

The present study adds to the growing body of research that indicates tick GSTs as a source of potential tick antigen. The following conclusions can be drawn from the present study that among the seven DmGSTs, mu-class GST exhibited better antigenicity in silico. Gene expression profiles supported that DmGSTM1 and DmGSTM2 are of potential value to be developed as anti-tick vaccine. Despite its predictive nature, this study offers some insight into finding tick GST as protective antigen; nevertheless, further studies need to be carried out in order to validate the antigenic effect of the DmGSTs.

Below is the link to the electronic supplementary material.

Electronic supplementary material 1 (DOCX 17 kb)

Electronic supplementary material 2 (DOCX 16 kb)

Electronic supplementary material 3 (DOCX 18 kb)

Electronic supplementary material 4 (DOCX 17 kb)