I A Orlov

G L Ataev

B Gourbal

A S Tokmakova

A V Bobrovskaya

E E Prokhorova

The transcriptomic analysis of Planorbarius corneus hemocytes (Gastropoda) naturally infected with Bilharziella polonica (Schistosomatidae)

Keywords:

The study of molluscan innate immunity is essential for understanding the evolution of the immune system. An advance in the knowledge of their immune system can be achieved by increasing the number of model species. Our study focuses on the immunity of Planorbarius corneus, a pulmonate snail widely distributed in Eurasia. These snails are intermediate hosts of many trematodes, including Bilharziella polonica (Schistosomatidae). In this paper we obtained and analyzed transcriptomes of hemocytes of uninfected snails Planorbarius corneus and snails naturally infected with Bilharziella polonica. The transcriptomes were found to contain transcripts encoding all major groups of immune factors previously described for other gastropods. Pathogen-recognition molecules were the most diverse group of immune factors. Comparison of the transcriptomes of the infected and the uninfected molluscs showed that the expression of some genes changed during infection. Our results extend the knowledge of immune responses of pulmonate snails to trematode invasion and promote P. corneus as a new model for the study of molluscan defence reactions KEYWORDS Planorbarius corneus -Bilharziella polonica -Hemocytes -

Introduction

Biomphalaria glabrata is a popular object of immunological research and the best studied gastropod species in this respect. This is mainly due to the fact that these snails are the intermediate hosts of schistosomes, including Schistosoma mansoni, the agent of human schistosomiasis. Most of the information about gastropod immunity has been obtained in studies made on Biomphalaria glabrata. Therefore, it remains unknown how conservative and universal immunityrelated mechanisms and factors reported in B. glabrata in different gastropod species. Moreover, the experiments are usually performed on laboratory molluscs, while disease-related events occur in natural settings, which can differ significantly from the laboratory conditions. This raises a question of whether we can confidently translate the laboratory findings to wild animals.

To answer this question, one should conduct studies on molluscs taken from the nature. Comparison of data obtained on molluscs from natural habitats and molluscs of laboratory strains can provide information about the impact of biological environmental constraints on their immune reactions. For a long time, the research community has lacked established models for studies of gastropod immunity, which hampered the progress in the field. A need for new model organisms was evident, and recently the situation has to change, albeit slowly.

Increasingly more studies have been focussed on the immunity in different gastropod species, such as a transcriptomic analysis of laboratory strains of Lymnaea stagnalis [START_REF] Seppala | Transcriptome profiling of Lymnaea stagnalis (Gastropoda) for ecoimmunological research[END_REF] and molluscs from natural populations of Littorina littorea [START_REF] Gorbushin | Lectin-like molecules in transcriptome of Littorina littorea hemocytes[END_REF], Biomphalaria pfeifferi [START_REF] Buddenborg | Transcriptomic responses of Biomphalaria pfeifferi to Schistosoma mansoni: investigation of a neglected African snail that supports more S. mansoni transmission than any other snail species[END_REF], and Physella acuta [START_REF] Schultz | RNA-seq: the early response of the snail Physella acuta to the digenetic trematode Echinostoma paraensei[END_REF]. In these studies, common groups of immune reaction factors are detected in all infected (more generally speaking, immunologically challenged) snails. However, there are also substantial differences between species, which may be due to physiological features, taxonomic position and the influence of the habitat, including the impact of pathogens that an organism encounters and coevolves with. It has become fairly obvious that in order to identify the general patterns of pulmonate immune reactions and to build a model of the gastropod immune system, it is necessary to expand the range of molluscan species under study and use molluscs from both laboratory lines and natural populations.

Planorbarius corneus is a promising model for the study of gastropod immunity. This species, common across Eurasia, is of particular interest because of its role as an intermediate host for a number of trematodes, including Bilharziella polonica and some other schistosomatids. B. polonica causes cercarial dermatitis, sometimes called the swimmer's itch. Bilharziella polonica -Planorbarius corneus host-parasite system has already been investigated in a number of studies [START_REF] Zbikowska | Thermal preferences of bird schistosome snail hosts increase the risk of swimmer's itch[END_REF] and can be considered, with some reservations, as analogous to Schistosoma mansoni -Biomphalaria glabrata pair.

In our previous studies, we have obtained valuable information about Planorbarius corneus hemolymph, its cell composition, haematopoiesis and immune system reactions to exposure to various immunity-inducing agents [START_REF] Ataev | The influence of trematode infection on the hemocyte composition in Planorbarius corneus (Gastropoda, Pulmonata)[END_REF]. Our data indicate that the immune mechanisms of P. corneus, especially the cell-based immunity, are in general very similar to those of Biomphalaria glabrata. These processes can now be productively explored in detail with the help of well-developed methods of next generation sequencing.

In this study, we performed a transcriptomic analysis of hemocytes, the main effectors of molluscan immune system, in Planorbarius corneus snails naturally infected with Bilharziella polonica. Our aims were to obtain the first ever transcriptomic data on Planorbarius corneus hemocytes, to provide an improved characterization of P. corneus as a model organism for gastropod immunity studies and to compare the processes occurring during natural infection with those recorded in laboratory snail strains. We described the main molecular factors of P. corneus immunity and recorded the response of these molluscs to Bilharziella polonica infection, thus significantly expanding our scarce knowledge about this non-model organism.

Materials and methods

Molluscs, hemolymph collection and RNA extraction

We collected snails Planorbarius corneus (n = 150) with the shell size ranging from 23 to 27 mm from Lake Finskoe in the outskirts of St. Petersburg (60.082,30.285) in the end of July in 2018. To define the schistosome infection status, the snails were kept separately in glasses with artificial spring water for three days, and individuals emitting cercariae were considered as infected. Only snails shedding larvae characteristic of Bilharziella polonica were involved in the experimental group (first confirmation of infection). B. polonica infection was additionally confirmed in each case by morphological examination of sporocysts during dissection (second confirmation of infection) and by rDNA PCR. In addition, we checked the snails for infection with other parasites and assessed their general state. Only B. polonica-infected snails and those exhibiting typical P. corneus characteristics were incorporated into our experiment. Fragments of rDNA were amplified by polymerase chain reaction (for the description of the PCR parameters, see [START_REF] Prokhorova | Analysis of ITS1 and ITS2 of ribosomal DNA in populations of Planorbarius corneus snails (Gastropoda) from the Leningrad and Kaliningrad oblasts of Russia[END_REF], sequenced and annotated in GenBank (MK264353.1).

Before all manipulations, the snails were thoroughly washed with distilled water and then rubbed with 70% ethanol. Hemolymph was collected from a sinus in the head area using a sterile Pasteur pipette and immediately centrifuged at 1000 g for 10 min at 4 • C to sediment the cells. The hemolymph from each snail (approx. 500 μL) was collected and processed individually, separately from the other samples. For each snail the number of cells in the hemolymph was determined using the Goryaev chamber. The snails with a very low or a very high hemocyte count were not included in the experiment. Therefore, for both groups of snails (infected and uninfected) we selected individuals with an average concentration of hemocytes in the hemolymph (400-470 cells per 1 μL). After the second confirmation of infection with B. polonica and the hemocyte count, we formed two samples, each consisting of 10 pooled hemocyte specimens (with an average concentration of hemocytes) from the corresponding group: uninfected (control) and infected animals. RNA was extracted from pooled samples using the Trizol reagent (Invitrogen) and its quality was controlled by capillary electrophoresis on the Agilent Bioanalyzer.

Library preparation, sequencing and read preprocessing

mRNA was enriched and converted to cDNA by RT with polyT primers (Mint-2 kit, Evrogen). Sequencing library with mean fragment size of approximately 450 bp was prepared using QIAGEN QIAseq FX DNA Library Kit, 100bp sequencing in paired-end mode performed on Illumina HiSeq 2500 system. At first, we removed from the obtained reads library-specific adapters and polyA tails, trimmed N's and bases with phred score <5 from both ends, unconditionally clipped 8 bases from 5 ′ and 5 bases from 3 ′ as we noticed a bias in nucleotide distribution there and discarded all reads shorter than 20 bp using cutadapt [START_REF] Marcel | Cutadapt removes adapter sequences from high-throughput sequencing reads[END_REF]. We mapped reads resulting from the previous step to chimeric rRNA-coding cluster with hisat2 [START_REF] Kim | HISAT: a fast spliced aligner with low memory requirements[END_REF] using standard parameters and discarded all the aligned reads (the entire pair was excluded even if only one read aligned). For this purpose, we constructed the rRNA-coding cluster mentioned from the Biomphalaria glabrata sequence replacing ITS1, ITS2 with the Planorbarius corneus sequences previously obtained in our group [START_REF] Prokhorova | Analysis of ITS1 and ITS2 of ribosomal DNA in populations of Planorbarius corneus snails (Gastropoda) from the Leningrad and Kaliningrad oblasts of Russia[END_REF].

Contamination clearing

All the pre-processed reads were classified with kraken2 software [START_REF] Wood | Improved metagenomic analysis with Kraken 2[END_REF] using NCBI NT database (https://www.ncbi.nlm.nih .gov/nuccore/) with confidence parameter equal to 0.6. We used BBduk (https://sourceforge.net/projects/bbmap/) to filter out reads containing 27-mers from specially constructed database. This database consisted of sequences (genomic, mtDNA, rRNA) from the observed contaminating organisms or species closely related to them (supplementary file 1).

De-novo transcriptome assembly and annotation

Reads from both samples (control and infected) were merged and used for a de novo assembly of a reference transcriptome using Trinity [START_REF] Grabherr | Full-length transcriptome assembly from RNA-seq data without a reference genome[END_REF]. Assembly quality evaluation incorporated assessment of technical metrics such as ones obtained by TransRate [START_REF] Smith-Unna | TransRate: reference-free quality assessment of de novo transcriptome assemblies[END_REF] and biologicallyrelevant number of universal single-copy orthologs (USCOs) detected by BUSCO [START_REF] Seppey | BUSCO: assessing genome assembly and annotation completeness[END_REF]used odb10). Transcripts longer than 200 bp were clustered based on 95% sequence identity using CD-hitEST [START_REF] Fu | CD-HIT: accelerated for clustering the next-generation sequencing data[END_REF]. ORFs were predicted and extracted with TransDecoder (https://github.com/T ransDecoder/TransDecoder/wiki). Encoded protein sequences were then clustered based on 99% identity and subsequently annotated with BlastP, e-value threshold equal to 1e-5, against NCBI NR database (ftp://ftp.ncbi.nlm.nih.gov/blast/db/, accessed on August 02, 2021) and InterProScan, with default settings, utilized to predict functional domains. In this annotation we searched for proteins containing immune-related domains or having similarities to previously identified polypeptides with immune function with the Python script. To perform a pathway-based functional annotation we exploited the KEGG database with its related tool BlastKOALA on a genus level. Assigned K numbers were mapped to KEGG modules using KEGG Mapper -Reconstruct tool (https://www.genome.jp/kegg/mapper/reconstruct.html).

Differential expression analysis

To calculate read counts of each transcript we utilized Salmon [START_REF] Patro | Salmon provides fast and bias-aware quantification of transcript expression[END_REF]. Obtained values alongside the offset matrix were imported into EdgeR [START_REF] Robinson | edgeR: a Bioconductor package for differential expression analysis of digital gene expression data[END_REF] package for R with the help of another package tximport [START_REF] Soneson | Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences [version 1; peer review: 2 approved[END_REF]. To conduct differential gene expression analysis with likelihood ratio tests (glmLRT), a conservative value of biological coefficient of variation BCV = 0.4 was specified as suggested in the package manual (https://bioconductor. org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsers Guide.pdf). Transcripts satisfying the following conditions were considered differentially expressed (DE): |logFC| (logarithm of fold-change) >2 and FDR (false discovery rate) < 0.05. Hemocytes from uninfected animals were considered as a reference for comparison, e.g., upregulated DEGs have a higher expression in the infected group compared to the uninfected molluscs. For DE transcripts we performed annotation with BlastX against NCBI NR database (e-value threshold equal to 1e-5).

Other methods

Alignment of identified FBG-containing proteins from P. corneus transcriptome to B. glabrata FRePs and calculation of corresponding percent identity matrices were performed with Clustal Omega [START_REF] Sievers | Fast, scalable generation of high-quality protein multiple sequence alignments using Clustal Omega[END_REF]. Alignments were manually inspected using MEGA 7 software [START_REF] Kumar | MEGA7: molecular evolutionary genetics analysis version 7.0 for bigger datasets[END_REF].

Results and discussion

Planorbarius corneus is a promising alternative model for studying the immune response of pulmonates infected in nature or in the laboratory-controlled environment. As confirmed in numerous studies, the main mechanisms of the immune response in P. corneus do not differ from those in Biomphalaria spp. In particular, the main effector cells of the immune reactions in both these snails are circulating hemolymph cells, the hemocytes [START_REF] Prokhorova | Analysis of ITS1 and ITS2 of ribosomal DNA in populations of Planorbarius corneus snails (Gastropoda) from the Leningrad and Kaliningrad oblasts of Russia[END_REF][START_REF] Prokhorova | Hemocytes of mollusc Biomphalaria glabrata (Gastropoda, Pulmonata)[END_REF]. Trematode invasion is known to affect the ratio of hemocyte types [START_REF] Ataev | The influence of trematode infection on the hemocyte composition in Planorbarius corneus (Gastropoda, Pulmonata)[END_REF]. Hemocytes are involved in cellular reactions such as phagocytosis of bacteria and encapsulation of metazoan parasites, especially sporocysts and metacercariae of trematodes. At the same time, they produce humoral immune factors which mark and opsonize any alien organism.

Gene expression analysis is a key tool for detailed functional activity studies of Planorbarius corneus hemolymph cells. We have previously analyzed the expression of several lectins in the hemocytes of infected and uninfected molluscs [START_REF] Prokhorova | Expression of genes encoding defence factors in the snail Planorbarius corneus (Gastropoda, Pulmonata) infested with trematodes[END_REF] and observed the influence of infection on the expression of some genes in hemocytes and hepatopancreas cells. However, to obtain a more comprehensive picture of the role of hemocytes in the immune challenge and to compare our results with those obtained on other species, a transcriptomic analysis is indispensable. In this work, we sequenced and analyzed the transcriptome of P. corneus hemocytes for the first time.

To study the transcriptome of hemocytes of P. corneus from natural populations under conditions of the background level of immunity challenge and under conditions of trematode invasion, we collected molluscs and verified their infection status. Infection with B. polonica was ascertained by an experienced morphologist, who examined shed cercariae and intramolluscan sporocysts. These two stages of Bilharziella polonica are clearly distinguishable from those of other species on the basis of morphological characters [START_REF] Combes | M´emoires du Mus´eum National d'Histoire Naturelle[END_REF]. B. polonica usually starts to infect the snails actively in late May-early June, and we collected the snails in the end of July. Shedding cercariae was a necessary condition for the inclusion of the snail into the infected group. Therefore, we may be sure that all molluscs in the infected group had a well-established infection.

Pre-processing, decontamination and assembly

Raw and pre-processed library sizes are given in Table 1. All the steps of read pre-processing and transcriptome assembly were performed using merged reads, unless indicated otherwise. After the first technical step, we additionally filtered out probable rRNA contamination, as it can be present even in polyA-enriched RNA preparations. For this purpose, we constructed chimeric rRNA-coding cluster from rDNA cluster of Biomphalaria glabrata with ITS replaced by those of Planorbarius corneus [START_REF] Prokhorova | Analysis of ITS1 and ITS2 of ribosomal DNA in populations of Planorbarius corneus snails (Gastropoda) from the Leningrad and Kaliningrad oblasts of Russia[END_REF]. We found that approximately 4% of pre-processed reads had rRNA origin, which is a modest value.

As already noted, before the sequencing experiment all the molluscs were dissected and scrupulously examined for the presence of visually detectable parasites and commensals. Nevertheless, since we worked with animals from a natural population, we additionally checked the experimental material for any possible contamination with the help of a bioinformatics method. For this, we separately classified pre-processed reads from both samples by kraken2 [START_REF] Wood | Improved metagenomic analysis with Kraken 2[END_REF] using NCBI NT database. As expected, most of them were classified as having molluscan origin. Schistosomal reads (about 50% of all contaminating sequences in the sample from infected animals) could be detected only in the library from infected snails. Traces of bacterial and plant sequences were present in both datasets. While the presence of parasite DNA in the infected sample could reflect the phagocytic activity of hemocytes, we also observed a significant and an almost equal portion of Annelida-originating sequences (approximately 20% of all contaminating sequences) in the infected and the control specimens. Annelid worms are well known to inhabit the mantle cavity of freshwater molluscs [START_REF] Callisto | Malacological assessment and natural infestation of Biomphalaria straminea (Dunker, 1848) by Schistosoma mansoni (Sambon, 1907) and Chaetogaster limnaei (K. Von Baer, 1827) in an urban eutrophic watershed[END_REF]. Therefore, we decided that the contamination occurred during hemolymph collection. Proportions of non-molluscan reads are shown in Fig. 1.

To eliminate contamination by schistosomes and annelids, we initially mapped the analyzed reads to 3 schistosome (Schistosoma mansoni, S. japonicum, S. haematobium) and 3 annelid (Eisenia fetida, Helobdella robusta, Hirudo medicinalis) genomes, discarding all the aligned sequences. Subsequent Kraken2 classification still showed a substantial level of contamination, even after modification of default parameters of splice-aware hisat2 [START_REF] Kim | HISAT: a fast spliced aligner with low memory requirements[END_REF]. To maximally filter out reads from contaminating species we utilized kmer based approach applying BBDuk and custom sequence database. This database included genomic and/or most common sequences (mtDNA, rRNA) from contaminating organisms themselves (defined by Kraken2) or from species evolutionarily close to them in case the data for the original species were lacking. The resulting library consisted of 44744691 paired reads of 88bp. Due to a long evolutionary coexistence and a relatively small phylogenetic distance between gastropods and schistosomes and, especially, between gastropods and annelids, a substantial similarity of some of their genomic sequences may be expected. A double-edged filtration method relying on k-mer occurrence inevitably leads to the loss of some information due to the presence of evolutionary conserved sequences. Nevertheless, we chose this method, and minimized the number of constructed chimeric or non-molluscan transcripts by sacrificing a portion of "true" data.

From the plethora of available transcriptome assemblers, we tried Trinity [START_REF] Grabherr | Full-length transcriptome assembly from RNA-seq data without a reference genome[END_REF], rnaSPAdes and Oases [START_REF] Schulz | Oases: robust de novo RNA-seq assembly across the dynamic range of expression levels[END_REF]. Unexpectedly, obsolete and non-supported Oases worked best for our data in terms of some technical metrics (average and median length) though not in terms of the fraction of mapped back reads (Table 2). Taking into account the modest amount of sequencing data, we finally chose Trinity-generated assembly as produced by the most conservative and reliable assembler, which is additionally supported by the highest number of mapped back reads.

Basic assembly parameters after 95% identity clusterization of transcripts can be found in Table 2. Our assembly included more than 90% of eukaryotic and metazoan (Table 2) universal single copy orthologs (USCOs) according to BUSCO [START_REF] Seppey | BUSCO: assessing genome assembly and annotation completeness[END_REF]. The fact that there were only 55% molluscan-specific USCOs could be explained by the design of our experiment where we sequenced only one specialised tissue which is unlikely to express all the organism-level USCOs. However, the low proportion of molluscan-specific USCOs could also indicate that the assembly was incomplete. We detected a high duplication rate of USCOs despite 95% sequence identity clusterization, which is rather common for transcriptomic assemblies.

There were 60630 ORFs the assembled transcripts (28744 complete and the others partial or internal), which resulted in 53611 polypeptide products after 99% protein sequence identity clusterization.

Transcriptome annotation

Transcripts corresponding to all the main groups of immunity factors previously described in other gastropods were identified in hemocytes of Planorbarius corneus (Table 3). We detected transcripts encoding more than 160 different immune-relevant proteins from various functional classes, such as molecules involved in pathogen recognition, adhesion and intercellular interactions, immune effectors and immune regulators.

Pattern recognition receptors (PRRs) were one of the best represented groups in our data. The entire diversity of PRRs was identified. In particular, we found transcripts encoding ten members of Toll-like receptors (TLR) family and five members of scavenger receptors (SR) family. The diversity of TLRs defines their role as the main receptors recognizing pathogen-associated molecular patterns (PAMPs) in the theory of innate immunity. In addition, the presence of a certain set of TLR genes is considered as a possible mechanism of specific compatibility in the trematode-mollusc system [START_REF] Pila | A novel Toll-like receptor (TLR) influences compatibility between the gastropod Biomphalaria glabrata, and the digenean trematode Schistosoma mansoni[END_REF].

In the widely studied snail species Biomphalaria glabrata, 56 genes encoding TLR of two different types (V-TLR and P-TLR) were identified [START_REF] Adema | Whole genome analysis of a schistosomiasis-transmitting freshwater snail[END_REF]. Analysis of hemocyte transcriptomes showed a higher level of constitutive TLR (BgTLR) gene activity in B. glabrata molluscs resistant to infection with Schistosoma mansoni (BS-90 line) compared with susceptible ones (M line). S. mansoni invasion leads to a 27-fold increase in the expression of this gene in the hemocytes of resistant molluscs, while the expression of this gene in the hemocytes of susceptible molluscs does not change [START_REF] Pila | A novel Toll-like receptor (TLR) influences compatibility between the gastropod Biomphalaria glabrata, and the digenean trematode Schistosoma mansoni[END_REF]. In our work, we did not detect any differences in TLR expression between infected and uninfected snails. It is important to note that the molluscs in our study were taken from nature, and so we do not know their genetic characteristics such as the degree of susceptibility to trematode invasion. However, we studied molluscs with mature infection, which indicates their susceptibility.

Transcripts corresponding to other PRR groups, as well as adhesion molecules, were also abundantly represented in the transcriptome of Planorbarius corneus by scavenger receptors, lectins of different groups, and adhesion molecules, including integrins, selectins, dermatopontins, etc. (Table 4). Besides, transcripts of all major groups of cytotoxicity factors of innate immunity cells were identified in hemocytes: antimicrobial proteins, homologues of complement system components, enzymes involved in the production of reactive oxygen and nitrogen species, enzymes of the phenoloxidase system, etc. These results confirm that hemocytes are universal cells of the innate immunity system of molluscs. Similarly to TLRs, there were no significant differences in expression between infected and uninfected molluscs for these transcripts. Although there is evidence that these factors are upregulated during experimental immunization of molluscs [START_REF] Wang | Pathogen-derived carbohydrate recognition in molluscs immune defense[END_REF], the increase in the expression has a clear temporal dynamic [START_REF] Chen | A novel toll-like receptor from Crassostrea gigas is involved in innate immune response to Vibrio alginolyticus[END_REF]. The main results regarding differential gene expression analysis are discussed in the corresponding section below.

The presence of PRRs and other molecules involved in innate immunity reactions is the basic characteristic of an immunocyte (an immunocompetent cell). The immune molecules mediating defense reactions have to be present in immunocytes in the amounts sufficient for an immediate execution of their function. An increase in the amount of immune factors in circulation is not necessarily due only to an increase in the number of transcripts in hemocytes. Their concentration in the hemolymph may be also elevated as a result of the increase in the number of hemocytes producing them during invasion.

We performed a large-scale functional annotation to reveal metabolic pathways whose elements are encoded by the assembled transcriptome with the help of KEGG tool BlastKOALA. We found 77 modules (metabolic pathways including specialised molecular complexes) that were complete and 33 modules that lacked one only block. Among the fully represented pathways there were some of the carbohydrate and energy metabolism including glycolysis, TCA cycle, pentose phosphate pathway and oxidative phosphorylation. Triacylglycerol, phosphatidylcholine, ceramide, sphingosine and fatty acid biosynthetic reactions enzymes as well as enzymes catalysing beta-oxidation of the latter were observed. There were also all the components of adenine, guanidine and pyrimidine ribonucleotides synthesis pathways and complete sets of elements mediating a number of amino acid metabolism reactions including their synthesis, interconversion and degradation.

Proteins with immune-related domains

According to the InterProScan, we found 5385 immunity-related domains in the transcriptome-encoded proteome and 4405 unique proteins were predicted (Table 4). Slightly less than 30% of them were proteins with IgSF or different lectin-containing domains.

We did not detect any sequences bearing IgSF and fibrinogen domains simultaneously. This feature is a signature of fibrinogen-related proteins (FRePs), which are well-known and well-characterized effectors of immunity in pulmonate molluscs, especially Biomphalaria glabrata [START_REF] Adema | Digenean-gastropod host associations inform on aspects of specific immunity in snails[END_REF]. Indeed, not all molluscan species rely on FRePs. Numerous other molluscs elaborated diversified lectins containing a single FBG domain named FReDs [START_REF] Adema | Digenean-gastropod host associations inform on aspects of specific immunity in snails[END_REF]. On the other hand, we have previously observed the expression of FReP-like proteins in this species under the burden of trematode invasion [START_REF] Prokhorova | Expression of genes encoding defence factors in the snail Planorbarius corneus (Gastropoda, Pulmonata) infested with trematodes[END_REF]. Taking into account both these facts, it is reasonable to suggest that while Planorbarius corneus does possess some classical FRePs in its arsenal, their amount might be lower than in case of the well-studied Biomphalaria glabrata. The main defence factors in Planorbarius corneus are probably FReDs and other types of lectins, which are present in abundance. In particular, we found three polypeptides incorporating C-type lectin or fucolectin domain alongside the motif with an immunoglobulin-like fold, which makes it possible to characterize them as CRePs [START_REF] Gorbushin | Derivatives of the lectin complement pathway in Lophotrochozoa[END_REF].

Taking into account the potential importance of FRePs for snail immunity, we took a closer look on 176 polypeptides with annotated FBG domain found in our transcriptomic data. Out of these proteins, 43 originated from complete, 37 from internal and 96 from partial predicted ORFs. Approximately in 15 percent of polypeptides with FBG domain it was adjoined towards the N-terminus with epidermal growth factor (EGF) or Ca-binding EGF-like domains. Another common neighbour was laminin domain. Some of these proteins included a motif indicating a possibility for aspartate or asparagine posttranslational modification by hydroxylation, which is often found in domains homologous to EGF. Other less frequently observed options were growth factor receptor domain, PAN/Apple domain, Kringle-like motif, and carbohydrate binding signatures such as galactose-binding and WSC domains. A couple of domains characteristic of complement system proteins, such as Sushi and Complement control module/SCR domain, were observed in conjunction with FBG domain as well.

To check if any of FBG-containing proteins appear to be similar to already published FRePs we downloaded previously submitted protein sequences of B. glabrata FRePs from NCBI protein database (174 entries, supplementary file 2). In order to keep only molecules with both domains we annotated protein domains in the downloaded dataset with IPS de novo. After 99% identity clusterization we aligned each of our 176 predicted FBG-containing proteins to the chosen set of B. glabrata FRePs using Clustal omega, calculated percent identity matrices [START_REF] Sievers | Fast, scalable generation of high-quality protein multiple sequence alignments using Clustal Omega[END_REF]. The same approach was utilized but using only IgSF domains from this set of FRePs.

After manual curation of calculated percentage identity matrices and the obtained multiple sequence alignments we inferred that the absolute majority of sequences shared a high similarity with B. glabrata FREPs only in FBG domain. We could find two notable sequences, which were isoforms of one Trinity gene (not to be confused with a gene in its general meaning), originating from predicted complete ORFs. Each of two contained a signal peptide at the N-terminus, and showed high similarity to IgSF domains of 13.1 and several K-type fibrinogen-related proteins of B. glabrata (alignments are in supplementary file 3 and 4).

Though FREPs are described in many molluscs, most of the characterized FREPs belong to one molluscan species, Biomphalaria glabrata. In a close species, B. pfeifferi, 249 sequences fitting the description of FREPs have been revealed. Twelve of them are full-sized homologues of B. glabrata FREPs [START_REF] Buddenborg | Transcriptomic responses of Biomphalaria pfeifferi to Schistosoma mansoni: investigation of a neglected African snail that supports more S. mansoni transmission than any other snail species[END_REF]. Infection of B. pfeifferi with Schistosoma mansoni miracidia significantly changes the expression profile of FREP: expression of FREP2 increases gradually, that of FREP7 and 14 increases significantly on the 3rd day post infection. At the same time, the expression of FREP3, which strongly reacts to schistosomal infection in Biomphalaria glabrata, increases in B. pfeifferi slightly only in the beginning of the invasion. Conversely, FREP5 is expressed much more strongly in B. pfeifferi than in B. glabrata [START_REF] Buddenborg | Transcriptomic responses of Biomphalaria pfeifferi to Schistosoma mansoni: investigation of a neglected African snail that supports more S. mansoni transmission than any other snail species[END_REF]. Thus, the change in FREP expression has a clear temporal dynamic.

The survey of FREPs in other gastropods has not shown such a broad variety of these proteins as in Biomphalaria spp. Only four FREPs genes have been identified in the genome of Physella acuta [START_REF] Schultz | RNA-seq: the early response of the snail Physella acuta to the digenetic trematode Echinostoma paraensei[END_REF]. Infection with Echinostoma paraensei trematodes and bacteria does not upregulate FREP expression but some other lectins respond positively to immunization. Similar results have been obtained for Lymnaea stagnalis, in which only two FREPS were detected [START_REF] Seppala | Transcriptome profiling of Lymnaea stagnalis (Gastropoda) for ecoimmunological research[END_REF]. One protein homologous to FREP3 was identified in Helisoma trivolvis [START_REF] Zhang | Representation of an immune responsive gene family encoding fibrinogen-related proteins in the freshwater mollusk Biomphalaria glabrata, an intermediate host for Schistosoma mansoni[END_REF]. Recently, genes encoding FREPs have been identified in deep-sea gastropods living in thermal springs: Gigantopelta aegis (18 genes) and Chrysomallon squamiferum (14 genes). The maximum expression of FREPs genes in these snails was shown not in hemocytes but in barrier structures: ctenidia and the mantle [START_REF] Lan | Hologenome analysis reveals dual symbiosis in the deep-sea hydrothermal vent snail Gigantopelta aegis[END_REF]. Five genes encoding FREPs were identified in the genome of Littorina littorea, and two FREP genes were identified in the gastropod mollusc Aplysia californica. However, their structure differs significantly from the Biomphalaria FREPs [START_REF] Gorbushin | Lectin-like molecules in transcriptome of Littorina littorea hemocytes[END_REF].

Thus, though FREPs are common in molluscs, their diversity in B. glabrata exceeds that in any other species. It has been suggested that the diversity of FREPs in Biomphalaria could have been formed as a result of a long-term coevolution with Schistosoma mansoni [START_REF] Bayne | Successful parasitism of vector snail Biomphalaria glabrata by the human blood fluke (trematode) Schistosoma mansoni: a 2009 assessment[END_REF][START_REF] Mitta | Compatibility polymorphism in snail/schistosome interactions: from field to theory to molecular mechanisms[END_REF][START_REF] Mitta | The compatibility between Biomphalaria glabrata snails and Schistosoma mansoni: an increasingly complex puzzle[END_REF][START_REF] Adema | Digenean-gastropod host associations inform on aspects of specific immunity in snails[END_REF]. The parasite's attempts to "deceive" the immune system of the host with the help of an increased polymorphism and mimicry led to the development of mechanisms ensuring more effective recognition and destruction of the parasite by the immune system of the mollusc. Production of FREPs was one of these mechanisms. In other molluscan species interacting with other trematode species or other pathogens, the mechanisms of their recognition and elimination could be improved with the help of different PRRs and other components of immunity.

Thus, two explanations of the absence of FREPs in the assembled transcriptome of Planorbarius corneus can be offered. The first is that the mollusc can utilize other molecules for neutralization of the parasite. The second explanation is a constitutively low level of FREPs expression, which made them impossible to detect in our experiment. In order to ascertain whether FREPs are absent or present in P. corneus, it is necessary to analyze its genome and to study its transcriptome at different periods of infection.

Differentially expressed genes

With the help of EdgeR software [START_REF] Robinson | edgeR: a Bioconductor package for differential expression analysis of digital gene expression data[END_REF] we identified 189 differentially expressed genes (DEGs) under the condition of B. polonica invasion: 155 were upregulated and 34 were downregulated. After manual curation 60 and 11 DEGs were found to have sensible Blast annotation (not uncharacterized, putative unknown proteins etc.) for upregulated and downregulated groups, respectively (see supplementary file 5 for the full list of DEGs with sensible Blast annotation). It is important to note that we also included in the analysis sequences present only in one of the conditions as it seems to be highly probable that some of the genes are completely inactive in uninfected animals or, on the opposite, strongly inhibited by the parasite in the infected animals. Since there were no experimental replicates, we used the value of 0.4 for biological coefficient of variation parameter during differential expression analysis. This value was considered conservative enough to minimize expected false-positive results as advised by EdgeR manual. Keeping in mind all the details, we consider our analysis as qualitative but by no means quantitative.

Molecules of viral origin were detected in both down-and upregulated groups of transcripts (Table 5). Two downregulated transcripts were related to retrotransposon sequences coding endonuclease-reverse transcriptase and pol polyprotein. Two more transcripts, one with an increased presence and one with a decreased presence, resembled the replication-associated protein of CRESS ssDNA viruses infecting snails and insects. The level of one more transcript, encoding snail nidovirus orf1ab polyprotein, showed an increase. It has been noted that retrotransposons are somehow involved into the relationship between trematodes and molluscs. It has even been suggested that they play a crucial role in the outcome of parasitic invasion or in the possible horizontal gene transfer between the host and the parasite [START_REF] Smith | PIWI silencing mechanism involving the retrotransposon nimbus orchestrates resistance to infection with Schistosoma mansoni in the snail vector, Biomphalaria glabrata[END_REF]. Our data obtained on P. corneus supports the idea that retrotransposons are involved in trematode-snail interaction, though they do not clarify the nature of this involvement.

Sixteen upregulated DEGs showed a high similarity to molluscan proteins or proteins from closely related species. There seem to be among them collagen alpha and cathepsin B-like polypeptide, which are known to be the effectors of anti-trematode snail immunity [START_REF] Lockyer | Biomphalaria glabrata transcriptome: cDNA microarray profiling identifies resistant-and susceptiblespecific gene expression in haemocytes from snail strains exposed to Schistosoma mansoni[END_REF][START_REF] Myers | Differences in cysteine protease activity in Schistosoma mansoni-resistant and -susceptible Biomphalaria glabrata and characterization of the hepatopancreas cathepsin B Full-length cDNA[END_REF][START_REF] Tetreau | Specific pathogen recognition by multiple innate immune sensors in an invertebrate[END_REF]. The former is suggested to participate in hemolymph coagulation and sporocyst encapsulation while the latter, being an endosomal cysteine protease, recycles heterogeneous material acquired by the cell.

Three transcripts encoding myosin heavy chain-like proteins (II-and 6/7-like) had an increased expression in the infected animals. The times when myosin was considered a muscle-specific polypeptide are long gone and it is now generally accepted that members of the large myosin superfamily can be present in almost every compartment of various cells, mediating the processes of cell adhesion, migration, endocytosis, phagocytosis, chromatin remodelling and others [START_REF] Vicente-Manzanares | Non-muscle myosin II takes centre stage in cell adhesion and migration[END_REF]. Under conditions of endogenous parasite development, hemocytes perform encapsulation and phagocytosis, and undergo proliferation, accompanied with a change in gene expression patterns, probably exploiting myosin molecules at each of these steps. The role of myosin in immunity is still obscure, as is that of many other immune challenge-altered factors. There is evidence that myosin II heavy chain protein was less abundant in sporocyst-encapsulating hemocytes from the resistant Biomphalaria glabrata line compared to non-encapsulating cells [START_REF] Dinguirard | Proteomic analysis of Biomphalaria glabrata hemocytes during in vitro encapsulation of Schistosoma mansoni sporocysts[END_REF]. It has also been reported that two myosin II heavy chain genes were expressed to a lesser degree in the hemocytes of susceptible molluscs exposed to miracidia compared to the hemocytes of the resistant strain [START_REF] Lockyer | Early differential gene expression in haemocytes from resistant and susceptible Biomphalaria glabrata strains in response to Schistosoma mansoni[END_REF]. Upregulation of transcripts for translation elongation factor 1-α and polyubiquitin-B in the infected group can be considered as a sign of more active protein synthesis and degradation, which is to be expected in hemocytes under immunological challenge. An increased level of the same elongation factor-coding sequence has been described in hemocytes of trematode-resistant B. glabrata strain [START_REF] Lockyer | Biomphalaria glabrata transcriptome: cDNA microarray profiling identifies resistant-and susceptiblespecific gene expression in haemocytes from snail strains exposed to Schistosoma mansoni[END_REF].

Different toxic proteins like biomphalysins are exploited by the mollusc during infection to fight the parasite. Though we did not detect any increase in the expression of transcripts encoding biomphalysin, we found another potentially toxic protein resembling the so-called turripeptide. Turripeptide neurotoxins discovered in the marine snail from the superfamily Conoidea can inhibit acetylcholine receptors, harming the parasite in this way [START_REF] Hernandez-Samano | A turripeptide from Polystira nobilis venom inhibits human α3β2 and α7 nicotinic acetylcholine receptors[END_REF]. Moreover, the turripeptide-like protein found in our study contains kazal-type domain according to IPS, thus probably possessing serine protease inhibiting activity. Serine protease inhibitors are well established effectors of snail immunity helping to block the digestive enzymes of the invaders. Another member of the turripeptide-like family has been reported to be overexpressed in non-resistant snails compared to the resistant line [START_REF] Lu | An overview of transcriptional responses of Schistosome-susceptible (M line) or -resistant (BS-90) Biomphalaria glabrata exposed or not to Schistosoma mansoni infection[END_REF]. The authors suggest that it is a susceptibility factor, which mediates tolerance to infection or perhaps contributes to animal survival. The same paper attributes a similar, tolerancemediating role to temptin, a pheromone driving mating behaviour, which had an increased expression in hemocytes of infected molluscs in our research [START_REF] Lu | An overview of transcriptional responses of Schistosome-susceptible (M line) or -resistant (BS-90) Biomphalaria glabrata exposed or not to Schistosoma mansoni infection[END_REF]. At the same time, temptin expression was upregulated in Mytilus galloprovincialis hemocytes stimulated with lipopolysaccharides and in Biomphalaria pfeifferi 3 days after infection with Schistosoma mansoni [START_REF] Buddenborg | Transcriptomic responses of Biomphalaria pfeifferi to Schistosoma mansoni: investigation of a neglected African snail that supports more S. mansoni transmission than any other snail species[END_REF][START_REF] Moreira | Stimulation of Mytilus galloprovincialis hemocytes with different immune challenges induces differential transcriptomic, miRNomic, and functional responses[END_REF].

Apolipophorins binding various lipids demonstrate immune-related properties in insects. Different members of this family are involved in pathogen pattern recognition, toxin detoxification and antimicrobial defence. The role of this family in molluscs is far from being understood but it has been shown that apolipophorin level is decreased in hemocytes of the Biomphalaria glabrata strain susceptible to trematode invasion [START_REF] Dinguirard | Proteomic analysis of Biomphalaria glabrata hemocytes during in vitro encapsulation of Schistosoma mansoni sporocysts[END_REF]. In our study we reported an increase in the expression of apolipophorin-like protein-coding transcript, which is an evidence that it has some immunity-related functions in snails as well as in arthropods.

The secreted vertebrate protein brorin acts as the antagonist of BMP activity, thus participating in neurogenesis and neural tissue development. We observed an upregulation of brorin-like protein gene expression in hemocytes of infected Planorbarius corneus. A vWFC related protein (VREP) similar to brorin has been found to be abundant in coelomic fluid of starfish Asterias rubens several hours after an induced blood loss [START_REF] Shabelnikov | Injury affects coelomic fluid proteome of the common starfish, Asterias rubens[END_REF]. The brorin-like protein detected in Planorbarius corneus might also mediate tissue repair processes aimed at minimizing the consequences of the ongoing infestation, similarly to what has been reported for the starfish.

Despite meticulous data pre-processing after the final assembly, we found in the hemocytes of the infected molluscs some transcripts encoding proteins with a very high similarity to schistosomal ones. We suggest that these sequences have a phagocytotic origin (see also above). It is likely that certain proteins of different trematodes share a high similarity, which is accompanied by a substantial nucleotide sequence divergence of the transcripts encoding them. This seems to be the reason behind such a high proportion of schistosome-related polypeptides in our assembled transcriptome even after a thorough read filtration. In principle, we could also apply a translational-based procedure during pre-processing, e.g. BlastX or DIAMOND. However, this approach would result in lots of false hits mainly due to conserved protein domains. For this reason, we eventually decided against using method and simply manually discarded schistosome proteins that could still be detected after the previous pre-processing steps.

Interestingly, we found a protein from Hsp70 family among other trematode housekeeping polypeptides such as major egg antigen, tubulin or COX III subunit. These proteins are known to localize in cercarial acetabular glands and play a role during the infection of the definitive host. Notably, Hsp70 with molluscan descent according to blast was found in the upregulated group as well. This protein in snail is related to the outcome of miracidium invasion and the course of the disease, but its role is yet unclear and probably depends on the stage of infection, the age of the mollusc, its genetic background and other factors [START_REF] Ittiprasert | Reversing the resistance phenotype of the Biomphalaria glabrata snail host Schistosoma mansoni infection by temperature modulation[END_REF]. The data on Hsp70 expression in hemocytes are also controversial. While a decreased expression in both resistant and susceptible molluscs has been reported in some studies, with the amplitude being lower for resistant molluscs [START_REF] Zahoor | Larval excretory-secretory products from the parasite Schistosoma mansoni modulate HSP70 protein expression in defence cells of its snail host, Biomphalaria glabrata[END_REF], a higher gene activity has been found in resistant snails in other studies, suggesting that Hsp70 exhibits protective properties [START_REF] Lockyer | Biomphalaria glabrata transcriptome: cDNA microarray profiling identifies resistant-and susceptiblespecific gene expression in haemocytes from snail strains exposed to Schistosoma mansoni[END_REF]. As indicated earlier, in our dataset Hsp70 demonstrated an elevated expression. Considering that we also observed some bacterial and protozoan transcripts in the hemocytes of the infected snails, which is an indirect evidence of compromised immunity, these finding suggest that an elevated snail Hsp70, despite its stress-alleviating properties, does not ultimately define the invasion failure.

The coding-transcript of chitin deacetylase, the protein of aminosugar metabolism, was downregulated. Hemolymph cells of some invertebrates contain chitin or molecules that bind to chitin-binding protein but their functionality is unclear [START_REF] Heath-Heckman | The occurrence of chitin in the hemocytes of invertebrates[END_REF]. It is not the first time when proteins related to synthesis or catabolism of chitin are found to change their expression during trematode invasion. For instance, chitinase-like proteins have been found to exhibit a complex time-dependent behaviour during Biomphalaria pfeifferi infection with Schistosoma mansoni [START_REF] Buddenborg | Transcriptomic responses of Biomphalaria pfeifferi to Schistosoma mansoni: investigation of a neglected African snail that supports more S. mansoni transmission than any other snail species[END_REF]. Another underrepresented transcript was that encoding a protein resembling ankyrin repeat domain-containing protein 17, which regulates the cell cycle and participates in antibacterial immune response in human cell lines [START_REF] Menning | A role for the Ankyrin repeat containing protein Ankrd17 in Nod1-and Nod2-mediated inflammatory responses[END_REF]. Downregulation of its transcript in Planorbarius corneus may be indicative of Bilharziella polonica strategy to attenuate snail immunity. However, the role of both these transcripts remains unknown, and so does that of the downregulated transcripts that share some identity with craniofacial development protein 2, filaggrin 2, human pG1 protein.

To sum up, the set of differentially expressed genes reported in our study represents a somewhat mixed snail phenotype. While some upregulated immunerelated genes seem to reflect ongoing attempts of the mollusc to supress the parasite and its reproduction, even at the cercariae-producing stage, others seem to minimize the damage and promote survival of the mollusc with an established infection. The overall picture results from the superposition of the mollusc's efforts to heal its lesions and facilitate immune response and the parasite's efforts to inhibit and overcome them.

The heterogeneity of the population is crucial for its stability under conditions of the changing environment. Some evidence indicates an animal population maintained in the laboratory loses much of its diversity already after several reproduction cycles (unpublished data). Thus, some findings made on laboratory animals may be due to a genetic drift rather than some generalised processes occurring in a free-living population. This is a serious concern to keep in mind. Comparison of our results obtained on free-living animals with the results of experiments on laboratory animals are necessary and planned by our group. It would help to determine the factors defining the success or the failure of the invasion, choose the markers for infection monitoring and assess the extent to which laboratory lines differ from the snails in the nature.

In this study we analyzed and published the first transcriptome of hemocytes of a widespread freshwater pulmonate snail Planorbarius corneus, with an emphasis on their immunity function. We provided an improved characterization of a new organism for gastropod immunity studies, making a step towards going beyond the Biomphalaria glabrata model. Our findings, alongside with the deposited raw data, provide a basis for further research in the field of parasitology and comparative immunology. Reasons to discard sequences and the proportions of reads dropped for a particular reason from the total number of dropped reads. These reasons include resemblance to material from contaminating organisms inferred by Kraken2 (Platyhelminthes, Annelida, plants), rRNA origin or low quality of the data. Likeness of reads to the genetic information from contaminating species was defined as presence of 27-mers from the database made up of probable contaminating sequences (accession numbers can be found in supplementary file 1). To reveal reads of rRNA origin, the reads were mapped to chimeric rRNA-coding cluster (see Methods).

Table 1

Sizes of initial and final libraries. Final library was obtained after standard pre-processing steps (e.g., discarding reads with too many N's, trimming of low-quality bases), removal of reads mapped to chimeric rRNA-coding cluster (see Methods) and decontamination of the data. The latter was achieved by filtering out all the reads sharing 27-mers with sequences from contaminating organisms (accession numbers can be found in supplementary file 1) as detected by Kraken2 or species closely related to them.

Control

Infected Total

Raw library size, 101 read pairs 24225942 27633188 51859130

Final library size, 88bp read pairs 21317065 23427626 44744691

Table 2

Technical and biological metrics of assemblies obtained using Trinity, SPAdes or Oases. Fraction of mapped reads was defined by TransRate which uses snap aligner [START_REF] Smith-Unna | TransRate: reference-free quality assessment of de novo transcriptome assemblies[END_REF]. BUSCO [START_REF] Seppey | BUSCO: assessing genome assembly and annotation completeness[END_REF] Functional groups of immune-related proteins encoded by identified in Planorbarius corneus hemocytes. Transcripts were checked for ORF presence, and one best ORF chosen for each transcript (as predicted by TransDecoder, https://github.com/TransDecoder/TransDecoder/wiki) was translated. The resulting proteome was annotated using BlastP against NR database (e-value<1e-5). In these results we searched for immunity-related proteins utilizing database comprising names of polypeptides with an immune function previously reported in the literature. We detected transcripts encoding more than 160 different immune-relevant proteins from various functional classes, such as molecules involved in pathogen recognition, adhesion and intercellular interactions, immune effectors and immune regulators.

Table 4

List of immunity-related domains found in the proteome encoded by the assembled transcriptome of Planorbarius corneus hemocytes and number of hits for the indicated categories detected by BlastP. IPS with default setting was exploited to predict protein domains using the same proteome as for previously described BlastP annotation. With the help of database comprising names of immune-related proteins and folds we eventually found 5385 immunity-related domains predicted by IPS in 4405 unique polypeptides. Most of these domains were represented by immunoglobulin superfamily and lectin folds.

Table 5

BlastX (e-value < 1e-5) annotation of differentially expressed transcripts determined by EdgeR, transcripts with microorganismal, trematodal or plant origin (according to BlastX) are not shown, for the full list see supplementary file 5. Due to the absence of replicates the analysis is considered to be qualitative, and the data are presented as up-or downregulated without a concrete fold change for every transcript (though Log|FC|>2 was necessary to consider the transcript to be differentially expressed). EdgeR was provided with the conservative numerical value of 0.4 for biological coefficient of variation to conduct the analysis.

Upregulated transcripts Downregulated transcripts

replication-associated protein orf1ab polyprotein collagen alpha-1(I) chain-like collagen alpha-

  Funding informationThe reported study was funded by RFBR according to the research project No. 20-54-15003; CNRS according to the research project No. 232465 and by the grant of the President of Russian Federation for state support of scholarly research by young scholars (project No.МК-1015.2021.1.4).

Fig. 1 .

 1 Fig. 1. Composition of final data used for Planorbarius corneus transcriptome assembly. (A) Proportion of removed sequences in initial data. (B)Reasons to discard sequences and the proportions of reads dropped for a particular reason from the total number of dropped reads. These reasons include resemblance to material from contaminating organisms inferred by Kraken2 (Platyhelminthes, Annelida, plants), rRNA origin or low quality of the data. Likeness of reads to the genetic information from contaminating species was defined as presence of 27-mers from the database made up of probable contaminating sequences (accession numbers can be found in supplementary file 1). To reveal reads of rRNA origin, the reads were mapped to chimeric rRNA-coding cluster (see Methods).

  with odb10 was implemented to assess the presence of Universal Single-Copy Orthologs (USCOs) in the assembled data.

	Parameter	Trinity	SPAdes	Oases
	Number of transcripts	548525	489507	488655
	Number of genes	478273	474345	-
	GC %	35.21	35.05	33.24
	Transcripts over 1 k	60550	73350	130406
	Transcripts over 10 k	6	3	7
	N50	665	831	1206
	Average length, bp	538.07	609.91	825.77
	Median length, bp	354	386	554
	Fraction of mapped reads	0.63	0.61	0.35
	Eukaryotic USCOs present	90.6%	92.2%	90.6%
	Metazoan USCOs present	92.1%	92.7%	92.2%
	Molluscan USCOs present	55.3%	56.0%	55.1%

Table 3

 3

Acknowledgements

The authors are grateful to Natalia Lentsman for her help with the preparation of the English version of the manuscript.

Appendix A. Supplementary data

Availability of data and material

Raw data generated in the current study is available for download from SRA: https://www.ncbi.nlm.nih.gov/sra/PRJNA865135. Data availability Data will be made available on request.

Ethics approval

Not applicable.

Declaration of competing interest

The authors declare that they have no competing interests. Supplementary data to this article can be found online at https://doi.org/10.1016/j.dci.2022.104607.


Joao Cavalcanti Santos

Lénaïc Cuau

Philippe Poignet

Nabil Zemiti

Decoupled Model Predictive Control for Path Following on Complex Surfaces

Decoupled Model Predictive Control for Path Following on Complex Surfaces

João Cavalcanti Santos 1 , Lenaïc Cuau 1 , Philippe Poignet 1 and Nabil Zemiti 1

Abstract-The present letter proposes a predictive path following control (PPFC) that controls 5 degrees of freedom (DoFs) of the end-effector while the remaining (decoupled) translational DoF should be controlled by an external controller. This PPFC is particularly useful for the path following on surfaces with geometric uncertainties such that the external controller can be independently designed to manage the interaction between the tool and the surface. Therefore, the proposed strategy turns out to be a versatile control scheme that can be integrated with external controllers designed for applications such as robotic surface finishing, welding and 3D printing on complex surfaces. The corresponding optimal control problem (OCP) considers mainly the positioning and orientation errors, tangential velocity and control input amplitudes. The proposed PPFC is validated experimentally in the context of robotic 3D bio-printing. A 7-DoF redundant manipulator equipped with a distance sensor is used to handle a print head through the desired print path. The distance measurements are used by an external controller to correct the printed layer height. The obtained accuracy is consistent with the repeatability of the used manipulator and computation time is compatible with high frequency controllers.

Index Terms-Optimization and Optimal Control; Medical Robots and Systems; Motion Control

I. INTRODUCTION

T HE motion control schemes taking as reference a curve within a robot workspace can be broadly divided into two groups: position tracking and path following. The former takes a time parameterized reference, i.e. a specific desired robot state is determined at each sampling time. This methodology is extensively used for the motion control of industrial manipulators [START_REF] Siciliano | Robotics: Modelling, Planning and Control[END_REF], [START_REF] Lewis | Robot manipulator control: theory and practice[END_REF]. In contrast, path following control strategies do not follow a temporal law. Typically, the control scheme aims at converging to the desired path (which is a subset of the robot workspace) along with a desired velocity profile. This kind of strategy is largely used in the control of mobile robots [START_REF] Sujit | Unmanned Aerial Vehicle Path Following[END_REF], [START_REF] Xu | Stabilization and path following of a single wheel robot[END_REF]. In applications in which time parametrization is not necessary, the advantages of path following are well-known. Since this approach drops the temporal law, path following controllers are, in general, less likely to lead to saturated control inputs and result in smoother convergence [START_REF] Lapierre | Robust nonlinear path-following control of an AUV[END_REF]- [START_REF] Matschek | Nonlinear Predictive Control for Trajectory Tracking and Path Following: An Introduction and Perspective[END_REF].

Model predictive control (MPC) is a suitable approach for the design of path following controllers [START_REF] Ritschel | Nonlinear model predictive path-following control for highly automated driving[END_REF]- [START_REF] Yang | Learning-based predictive path following control for nonlinear systems under uncertain disturbances[END_REF]. Thanks to its predictive nature, a MPC strategy is able to find an optimal set of control inputs anticipating changes within its predictive horizon [START_REF] Faulwasser | Implementation of Nonlinear Model Predictive Path-Following Control for an Industrial Robot[END_REF], [START_REF] Yang | Learning-based predictive path following control for nonlinear systems under uncertain disturbances[END_REF]. Moreover, system constraints Manuscript received: September, 16, 2022; Revised December, 15, 2022; Accepted January, 16, 2023. This paper was recommended for publication by Editor Lucia Pallottino upon evaluation of the Associate Editor and Reviewers' comments. This work was supported by the French project ANR-21-ASM2-0001 -Bloc-Print 2 1 The authors are with the LIRMM, Univ Montpellier, CNRS, France can be explicitly considered in the formulation of the optimal control problem (OCP). This characteristic is of particular interest since the best performance is often obtained in the limits of the system capabilities [START_REF] Maciejowski | Predictive control: with constraints[END_REF]. Path following strategies taking advantages of these characteristics are referred to as predictive path following control (PPFC), e.g. [START_REF] Faulwasser | Nonlinear Model Predictive Control for Constrained Output Path Following[END_REF]- [START_REF] Yang | Learning-based predictive path following control for nonlinear systems under uncertain disturbances[END_REF]. Among all the practical cases in which a robotic manipulator should follow a given trajectory, numerous applications involve a desired path projected on a specific surface, e.g. robotic welding, surface finishing, milling and 3D printing. Typically, the robot degrees of freedom (DoFs) related to the interaction between the tool and the surface should respond in a significantly different manner than the remaining DoFs. For instance, in the case of surface finishing, one may constrain the applied contact force perpendicular to the polished surface while the remaining DoFs should follow a given path and velocity profile. Solutions to deal with such an issue include the classic hybrid force/position control [1, Section 9.7], adaptive variable impedance [START_REF] Duan | Adaptive variable impedance control for dynamic contact force tracking in uncertain environment[END_REF] and neuro-adaptive control [START_REF] Karayiannidis | Force/position tracking for a robotic manipulator in compliant contact with a surface using neuro-adaptive control[END_REF].

Hybrid force/position control, as in the early work [START_REF] Raibert | Hybrid Position/Force Control of Manipulators[END_REF], decouples the DoFs in which position tracking and force control should be applied. This approach has the advantage that the control scheme responsible for tracking the desired path can be designed independently from the controller managing the interaction with the surface (henceforth referred to as external controller), as illustrated in Figure 1 for the case studied in this letter. Nevertheless, most of the studies addressing the integration of different control schemes along different DoFs apply position tracking approaches [START_REF] Karayiannidis | Force/position tracking for a robotic manipulator in compliant contact with a surface using neuro-adaptive control[END_REF], [START_REF] Pliego-Jiménez | Adaptive position/force control for robot manipulators in contact with a rigid surface with uncertain parameters[END_REF].

The design of path following control strategies for robotic manipulators capable of integrating different control schemes along different DoFs is scarcely addressed in the state of the art [START_REF] Dahroug | Visual servoing controller for time-invariant 3D path following with remote centre of motion constraint[END_REF], [START_REF] Wen | A novel 3D path following control framework for robots performing surface finishing tasks[END_REF]. Dahroug et al. propose in [START_REF] Dahroug | Visual servoing controller for time-invariant 3D path following with remote centre of motion constraint[END_REF] a visual servoing controller able to follow a desired path keeping a constant remote center of motion. These different aspects are integrated thanks to the application of a task priority strategy. Similarly, Wen and Pagilla introduce in [START_REF] Wen | A novel 3D path following control framework for robots performing surface finishing tasks[END_REF] a procedure in which an optimal path is determined based on a sequence of check points. The obtained path is used in an overall scheme integrating path following and force control. The introduced strategy is used for surface polishing. Both studies [START_REF] Dahroug | Visual servoing controller for time-invariant 3D path following with remote centre of motion constraint[END_REF], [START_REF] Wen | A novel 3D path following control framework for robots performing surface finishing tasks[END_REF] were focused in their corresponding applications. An extension of these schemes to different cases applying a general external controller is not addressed. Most importantly, they have in common that their path following controllers use Frenet-Serret frames [START_REF] Lapierre | Robust nonlinear path-following control of an AUV[END_REF] in order to compute tangential and transversal strategies but without using MPC. This leaves room for the aforementioned improvements obtained with predictive controllers. To the best of the authors' knowledge, the existing PPFC strategies do not integrate an external controller along a particular set of DoFs.

The main contribution of the present letter is the formulation of a PPFC strategy that can be integrated to an external controller managing independently the displacements along the direction normal to a given surface of interest. The proposed PPFC controls 5 tool DoFs while the remaining translational DoF normal to the surface (referred to as decoupled DoF) is controlled by the external controller. As a result, the obtained overall controller can perform path following on a surface with geometric uncertainties in which the behavior of the decoupled DoF controls the interaction (such as distance or force) between the tool and the surface. Since the external control can be designed independently, the interaction between the tool and the surface can be adapted to different applications, e.g. robotic surface finishing, welding and 3D printing. In order to perform the path following, the proposed OCP considers the tangential velocity along with the positioning and orientation errors with respect to the reference curve. The proposed OCP has the particularity that the positioning error with respect to a direction normal to the surface of interest is disregarded, since this DoF is controlled by the external controller.

A numerical optimization algorithm is proposed in order to meet the real-time constraints on the solution of the OCP. Applying such an algorithm, the strategy is validated experimentally in the context of 3D bio-printing. The printing path is generated on an initially unknown skin surface using a RGB-D camera. In this application, the surface geometric uncertainties are mainly caused by the errors related to the calibration between the camera coordinate frame and the robot base frame. As a consequence, the corresponding errors on the distance between the extruder and the skin surface for the theoretical printing path would hinder the appropriate deposition of the material. Similarly to conventional fused deposition modeling (FDM) 3D printing, the distance between the extruder and the surface of interest plays a particularly important role on the final result. Accordingly, the distance between the extruder and the skin surface are corrected with an external controller using a dedicated distance sensor. Satisfactory positioning error and computation time are obtained.

The remainder of this letter is organized as follows. After formulating the addressed problem in Section II, Section III introduces the PPFC strategy, detailing the main aspects considered in the corresponding OCP. The numerical solution of this optimization problem is discussed in Section IV. The experimental validation is detailed in Section V and conclusions are drawn in Section VI.

II. PRELIMINARIES

As depicted in Figure 1, this letter addresses a motion control problem in which the desired path is projected on a complex surface. The desired path consists of a sequence of desired tool positions and a sequence of desired tool orientations. These sequences are denoted as

P d = {p d,1 , p d,2 , . . . , p d,N } and N d = {n d,1 , n d,2 , . . . , n d,N }.
The strictly positive integer N represents the number of reference points while vectors p d,k , n d,k ∈ R 3 determine the desired tool pose and orientation, respectively.

Typically, n d,k is a unit vector normal to the surface at point p d,k . As a result, one may note that a tuple {p d,k , n d,k } constrains 5 DoFs of the tool. The remaining sixth DoF (the tool rotation around each n d,k ) is considered as a redundant DoF. This is the case, for instance, for robotic welding, sanding, surface finishing and 3D printing.

Consider a n-DoF manipulator with joint positions given by q ∈ R n , joint velocities q ∈ R n and n ⩾ 6. The proposed control strategy takes as state vector x = q T qT T and aims at managing these states in order to follow the desired path {P d , N d }.

III. DECOUPLED PPFC

After introducing the system model in Section III-A, the PPFC algorithm is discussed in Section III-B.

A. Model

In order to reduce the real-time computational burden, the control input is considered as the robot joint accelerations u = q. Similar approaches are commonly used in the state of the art, e.g. [START_REF] Yang | Learning-based predictive path following control for nonlinear systems under uncertain disturbances[END_REF]. Accordingly, it is considered that an internal controller is able to track constant joint accelerations for sufficiently small sampling periods ∆t. Therefore, the resulting discrete-time dynamic system determines the (i + 1) th state x i+1 based on the i th state x i and control input u i+1 ∈ R n according to

x i+1 = q i+1 qi+1 = q i + qi ∆t + u i+1 ∆t 2 / 2 qi + u i+1 ∆t = A x i + B u i+1 (1)
with constant matrices A ∈ R 2n×2n and B ∈ R 2n×n . The considered model is hence a linear time-invariant system. Accordingly, for a given initial state x 0 = q T 0 qT 0 T , the i th state vector after the application of a sequence of control inputs u 1 , . . . , u i can be written as Consider also that the kinematic model of the manipulator is known, such that the tool position p ∈ R 3 and orientation n ∈ R 3 can be computed applying the forward kinematics for a given q. As discussed in Section II, note that the unit vector n represents uniquely the tool axis, i.e. the DoF corresponding to the rotation around n is considered as redundant in the present model. The tool linear ṗ and angular ω velocities are computed based on the differential kinematics through ṗT ω T T = J(q) q, where J ∈ R 6×n is the manipulator jacobian matrix. This matrix is computed as a function of the joint positions q for each sampling period.

x i = A i x 0 + A i-1 B A i-2 B . . . B    u 1 . . . u i    (2)

B. Control Scheme

The proposed path following control is formulated as an MPC in which, for each controller cycle, an OCP is solved. Considering an actual measured state x 0 , the next states x 1 , . . . , x hp can be predicted using (2) for a given sequence

of control inputs µ = u T 1 . . . u T hp T
. The strictly positive integer h p is called prediction horizon. The solution of the OCP leads to an optimal µ * that minimizes a cost function J hp . This function should summarize the desired behavior of the system. The design of J hp is discussed in the remainder of this section.

The OCP cost function J hp (µ, x 0 , u 0 ) is a sum of stage costs ℓ(i, µ, x 0 , u 0 ) considering the control inputs and the i th predicted state throughout the prediction horizon:

J hp (µ, x 0 , u 0 ) = hp i=1 ℓ(i, µ, x 0 , u 0 ), (3)
where the stage cost ℓ(i, µ, x 0 , u 0 ) is given by

ℓ(i, µ, x 0 , u 0 ) = ℓ g (q i ) + ℓ v (q i , qi )+ k u ∥u i ∥ 2 + k du ∥u i -u i-1 ∥ 2 + k q∥ qi ∥ 2 , (4)
with positive scalars k u , k d,u and k q. The last three terms in (4) aim at smoothing the trajectory by minimizing the joint velocities and their derivatives. Vectors q i and qi in (4) are obtained from (2) for given x 0 and µ. The terms ℓ g (q i ) and ℓ v (q i , qi ) in ( 4) deal with the positioning error and the tool velocity, respectively. As illustrated in Figure 1, the PPFC does not control the translation along the tool axis n. Therefore, regarding the tool translation, the OCP should consider the positioning error of p with respect to the two axes perpendicular to n. Additionally, the alignment error ϕ related to the orientation of n itself also should be considered. These two aspects are considered in ℓ g (q i ) and illustrated in Figure 2.

The desired path is considered as a linear interpolation between consecutive p d,k and p d,k+1 . For sufficiently dense P d and N d (as defined in Section II), the minimization of the joint velocities and their derivatives in (4) leads to a smooth trajectory taking this linearly interpolated path.

In order to compute the error on the tool position p, one should first determine which segment {p d, k, p d, k+1 } should be considered. To this end, the closest segment to p is obtained using the Algorithm 1, cf. Figure 3.

Algorithm 1 Computation of k(p, k 0 ) and λ(p, k 0 ) Input: k 0 , p, p d,j for j ∈ {k 0 , k 0 + 1, . . . } Output: k, λ 1: k ← k 0 2: λ ← e T k (p -p d,k ) 3: while λ ⩾ ∆p d,k do 4: k ← k + 1 5: λ ← e T k (p -p d,k ) 6:
if λ < 0 then end if 10: end while 11: k ← k For a given positive integer k 0 , consider that the tracking of p d,k0 has already been performed. Thus, the previous segments {p d,k , p d,k+1 } : k < k 0 should not be considered in the algorithm. Let e k be the unit vector

e k = p d,k+1 -p d,k ∥p d,k+1 -p d,k ∥ = p d,k+1 + p d,k ∆p d,k , (5)
and

∆p d,k = ∥p d,k+1 -p d,k ∥. The distance between p d,k0
and the projection of p on the straight line defined by p d,k0

and p d,k0+1 can be written as e T k0 (p -p d,k0 ). If the length e T k0 (p -p d,k0 ) is greater than ∆p d,k0 , the next segment is considered. This procedure is repeated until e T k (p -p d,k ) < ∆p d,k . It is possible to prove that the Algorithm 1 necessarily converges to a finite k under the assumption that there is a positive scalar ε such that every ∆p d,k > ε. Additionally, since the sequence k = {k 0 , k 0 +1, . . . } is considered with an increasing order, the proposed algorithm is able to track selfintersecting curves. This ability is illustrated in the attached video. The convergence proof of Algorithm 1 and its capability to track self-intersecting curves are detailed in [START_REF] Santos | Supplementary Materials -Decoupled Model Predictive Control for Path Following on Complex Surfaces[END_REF].

Similarly to the desired path, the desired orientations are interpolated between two consecutive n d,k and n d,k+1 . For this purpose, let the unit vector r k ∈ R 3 and the positive scalar θ k satisfy

n d,k+1 = R(r k , θ k ) n d,k , (6)
where R(r k , θ k ) represents the rotation matrix defined by the axis of rotation r k and angle θ k . The desired orientation n d,t for a given tool position p, is defined as

n d,t (p) = R r k, c k λ n d, k, (7)
with k, λ computed using Algorithm 1 and

c k = θ k /∆p d,k . Note that n d,t = n d, k if λ = 0, and n d,t = n d, k+1 if λ = ∆p d, k.
Ideally, since the PPFC does not control the translation along the tool axis n, the tool position p should follow the blue surface Ψ k depicted in Figure 2. This surface contains the straight line segment defined by p d,k , p d,k+1 with interpolated orientation according to [START_REF] Matschek | Nonlinear Predictive Control for Trajectory Tracking and Path Following: An Introduction and Perspective[END_REF] for 0 ⩽ λ ⩽ ∆p d,k . This surface can be written as

Ψ k = z ∈ R 3 | ∃ λ, ρ ∈ R : z = p d,k + λ e k + ρ R(r k , c k λ) n d,k . (8)
Nevertheless, the distance between p and this complex surface cannot be obtained analytically and its computation (based on numerical optimization) would significantly increase the computational burden related to the evaluation of J hp . Accordingly, Ψ k is approximated by the plane Π k illustrated in Figure 2. This plane contains the segment p d,k , p d,k+1 and is oriented according to an intermediary vector n m,k defined as

n m,k = (I -e k e T k ) R(r k , θ k /2) n d,k (I -e k e T k ) R(r k , θ k /2) n d,k , (9)
where the term (I -e k e T k ) is used to insure that n m,k is perpendicular to e k .

Therefore, denoting E k = e k n m,k , the projection of p in the plane Π k can be computed as

p Π,k = p d,k + E k E T k (p -p d,k ) (10)
and the distance d k in Figure 2 is written as

d k = (I -E k E T k ) (p -p d, k) . (11)
Reminding that ℓ g (q i ) should quantify the positioning error of p and the alignment error of n, the distance d k and the

u fb q' + + f q, q q f q s ec q d (b) (c) (d)
Fig. 4. Overall control scheme.

scalar product n T n d,t are used for these respective goals so that

ℓ g (q) = k d (I -E k E T k ) (p -p d, k) 2 - k a (n T n d,t ) 2 , (12)
with positive scalar gains k d and k a . Finally, the partial stage cost ℓ v (q, q) should penalize the difference between the velocities along e k and the desired velocity ṗd . This is performed with

ℓ v (q, q) = k v e T k ṗ -ṗd 2 , (13)
with a positive scalar k v . In summary, the proposed PPFC takes the following OCP:

µ * = arg min µ J hp (µ, x 0 , u 0 ) (14a)
s. t. q min ⩽ q i ⩽ q max , i = 1, . . . , h p (14b)

-qmax ⩽ qi ⩽ qmax , i = 1, . . . , h p (14c)
with J hp defined in (3) and q min , q max , qmax the minimum joint position, maximum joint position and maximum joint velocity, respectively. The feedback output from the control scheme is u f b (x 0 , u 0 ) = u 1,µ * , i.e. the first n elements of µ * . Matrices J p , J r ∈ R 3×n are defined as submatrices of J such as J = J T p J T r T and ṗ = J p q, ω = J r q.

For a given actual state q T 0 qT 0 T , the desired joint accelerations u f b computed through ( 14) can be used within a tracking control at joint level. Typically, desired joint velocities can be sent to the internal robot driver. Figure 4 depicts the application of the proposed PPFC in an overall control scheme. The PPFC output of (a) is integrated over time in (b), leading to q′ f . One may note that partial stage costs ℓ g and ℓ v present constant values for positions and velocities varying along n m, k. Therefore, the minimization of ℓ g and ℓ v should not lead to displacements along the direction normal to the surface. Nevertheless, the minimization of q and its derivatives in (4) may induce non-null velocities along n, i.e. n T J p q′ f ̸ = 0. For this reason, block (c) projects q′ f within a subspace of R n , such that n T J p q′ f = 0, insuring that the tool velocity along n generated by qf is null. This procedure is equivalent to the application of a selection matrix (as in [START_REF] Wen | A novel 3D path following control framework for robots performing surface finishing tasks[END_REF], for instance).

Regarding the control of the motion along n, an external controller (d) computes qec based on a given feedback signal s. The measurements s may be obtained, for instance, with force or distance sensors. Note that qec should generate uniquely tool velocities along n, i.e. n n T J p qec = J p qec and J r qec = 0.

IV. NUMERICAL OPTIMIZATION

The applicability of a MPC scheme in a real-time robotic system depends on the computation time necessary to solve the corresponding OCP. One may note that the OCP ( 14) is nonlinear and may have a significant amount of arguments for a typical prediction horizon, e.g. 56 = n × h p for n = 7, h p = 8 as in Section V. Therefore, a detailed discussion on the numerical solution of ( 14) is necessary. The present section proposes an optimization algorithm specifically designed for this goal.

The proposed algorithm described in Section IV-A is derived from the one developed in [START_REF] Santos | A Nonlinear Model Predictive Control for the Position Tracking of Cable-Driven Parallel Robots[END_REF], in which the position tracking of cable-driven parallel robots is addressed. Since a gradientbased solver is used, the computation of the derivatives of J hp is necessary. The numerical estimation of these derivatives requires the computation of the overall cost function (3) at least n×h p +1 times. This may significantly impair the efficiency of the controller. As a solution to this issue, Section IV-B details the analytical expressions for the derivatives of J hp .

A. Overall Algorithm

Consider given x 0 and u 0 , such that the sought solution is µ * minimizing J hp (µ, x 0 , u 0 ), as in [START_REF] Karayiannidis | Force/position tracking for a robotic manipulator in compliant contact with a surface using neuro-adaptive control[END_REF]. For a given µ, (2) leads to predicted tool positions p i for i = 1, . . . , h p and corresponding sequence of integers K = { k1 , . . . , khp } computed according to Algorithm 1. A particular numerical issue related to the computation of µ * is that variations of µ would lead to variations of the integers ki , which results in discontinuous J hp . This would impair the efficiency of purely gradient based solvers. For that reason, the algorithm described in Algorithm 2 is used. end if 10: end loop In short, the Algorithm 2 iteratively solves simplified versions of ( 14) in which the sequence K is taken as constant. More precisely, the steps 3-10 of Algorithm 1 are bypassed in order to prevent the step k ← k + 1, avoiding the variations on K. Accordingly, the optimization problem in Algorithm 2 (step 4) can be solved efficiently with, for instance, sequential quadratic programming (SQP) algorithms. It is worth noting that the convergence of Algorithm 2 is more reliable for reference points P d that are sufficiently close to each other, reducing the variation of the cost function (14a) for successive iterations of the algorithm. The results presented in Section V were obtained with a SQP adapted from the solver used in [START_REF] Santos | A Nonlinear Model Predictive Control for the Position Tracking of Cable-Driven Parallel Robots[END_REF]. This numerical algorithm applies an altered version of the nonlinear gradient project method described in [21, Section 18.6], cf. [START_REF] Santos | Model Predictive Tracking Control of Cable-Driven Parallel Robots : From Concept to Real-Time Validation[END_REF]Chapter 4]. Basic matrix operations were performed using tools provided by the Eigen C++ library [START_REF] Guennebaud | Eigen v3[END_REF].

B. Derivatives of J hp

This section presents the analytical derivatives of J hp for constant K. The analytical derivatives are obtained straightforwardly based on the stage cost described in Section III-B. The complete deduction of these expressions are omitted for the sake of conciseness.

Consider the functions f : R m1 → R and g : R m1 → R m2 for any strictly positive integers m 1 and m 2 . For given input variable y ∈ R m1 , the derivatives of such functions are denoted as

∂ ∂ y f (y) = ∂ ∂ y1 f (y) . . . ∂ ∂ ym 1 f (y) and ∂ ∂ y g(y) =     ∂ ∂ y1 g 1 (y) . . . ∂ ∂ ym 1 g 1 (y) . . . . . . . . . ∂ ∂ y1 g m2 (y) . . . ∂ ∂ ym 1 g m2 (y)     . (15)
The gradient of the cost function in ( 14) is computed through

∂ ∂ µ J hp (µ, x 0 , u 0 ) = hp i=1 ∂ ∂ µ ℓ(i, µ, x 0 , u 0 ). (16)
Therefore, the gradient ( 16) can be computed by means of the derivatives of each ℓ(i, µ, x 0 , u 0 ) with respect to µ. Denoting

ℓ u,i = ℓ u (i, µ, u 0 ) = k u ∥u i ∥ 2 + k du ∥u i -u i-1 ∥ 2 ,
the following expression can be deduced for i > 1:

∂ ℓ u,i ∂ µ =2 0 T (i-2)×n -k du (u i -u i-1 ) T . . . k u u T i + k du (u i -u i-1 ) T 0 T (hp-i-1)×n with 0 m = 0 . . . 0 T ∈ R m . Similarly, for i = 1, ∂ ℓ u,i ∂ µ = 2 k u u T 1 + k du (u 1 -u 0 ) T 0 T (hp-1)×n .
The elements of ℓ(i, µ, x 0 , u 0 ) that depend on the state variables are first derived with respect to q i and qi . Clearly,

∂ ∂ q k q∥ q∥ 2 = 2 k q qT . (17)
Likewise, the derivatives of ℓ g,i = ℓ g (q i ) are given by

∂ ℓ g,i ∂ q i =2 k d d ki w T ki J p (q i )- 2 k a (n T n d,t ) n T ∂ n d,t ∂ q i + n T d,t ∂ n ∂ q i , (18)
with w ki = p i -p Π, ki /∥p -p Π, ki ∥ a unitary vector perpendicular to the plane Π ki . Integer ki is computed as a function of p i using Algorithm 1 and n d,t is obtained with [START_REF] Matschek | Nonlinear Predictive Control for Trajectory Tracking and Path Following: An Introduction and Perspective[END_REF]. The derivatives of n d,t and n can be written as

∂ ∂ q n d,t = θ ∆p d (r × n d,t ) e T J p (q) ∂ ∂ q n = -n × J r (q).
The derivatives of ℓ v,i = ℓ v (q i , qi ) are given by

∂ ℓ v,i ∂ q = 2 k v (e T ṗ -ṗd ) e T J p (19a) ∂ ℓ v,i ∂ q = 2 k v (e T ṗ -ṗd ) e T n j=1 ∂ j p,j ∂ q qj (19b)
where j p,j ∈ R 3 is the j th column vector of J p . Its derivatives

∂ j p,j ∂ q = ∂ jp,j ∂ q1 . . . ∂ jp,j ∂ qn
have vector columns computed as ∂ j p,j ∂ q l = j r,l × j p,j , if j ⩾ l j r,j × j p,l , if j < l where j r,j is the j th column vector of J r . For a given v ∈ R n representing a gradient with respect to q i , the multiplication v T ∂ q i /∂ µ considering the prediction (2) can be written as

v T ∂ q i ∂ µ = (∆t 2 ) 1 2 + (i -1) v T . . . 1 2 + (i -2) v T . . . 1 2 v T 0 T (hp-i)×n (20)
Likewise, the multiplication v T ∂ qi /∂ µ is computed using

v T ∂ qi ∂ µ = ∆t v T . . . v T 0 T hp-i . (21)
Replacing v T in ( 20)-( 21) by the derivatives in ( 17)-( 19), the gradient of the stage cost is given by

∂ ∂ µ ℓ(i, µ, x 0 , u 0 ) = ∂ ℓ u,i ∂ µ + ∂ ℓ g,i ∂ q i + ∂ ℓ v,i ∂ q i ∂ q i ∂ µ + 2 k q qi T + ∂ ℓ v,i ∂ qi ∂ qi ∂ µ .

V. EXPERIMENTAL RESULTS

The experimental validation presented in this section is inspired by the application of the proposed PPFC for 3D robotic bio-printing used in the treatment of serious skin wounds. The goal is to deposit a skin substitute over the wound surface using a print head embedded on a robotic arm endeffector. The used set-up is depicted in Figure 5.

A Revopoint Acusense RGB-D camera is used to determine the three-dimensional geometry of the surface on which the skin substitute should be applied. Based on this threedimensional geometry, a print path covering the whole wound can be generated. The print path consists of a sequence of desired tool positions p d,k and the corresponding unit vectors n d,k normal to the surface at position p d,k . The printing is performed using a 7-DoF Franka Emika Panda robotic arm to position a ViscoTec vipro-HEAD 5 print head along the desired path. The substitute skin is simulated using Pluronic F-127 + HBSS solution mixture.

The definition of the vectors p d,k and n d,k relies on the calibration of the camera coordinate system with respect to the robot base coordinate system. Using the materials described in this letter, such a procedure is subject to errors of ≈ 2 mm, cf. [START_REF] De Oliveira | A hand-eye calibration method for augmented reality applied to computer-assisted orthopedic surgery[END_REF]. Indeed, errors with this magnitude affecting the distance between the extruder and the skin surface would play a critical role.

Using the set-up presented in Figure 5, this distance should be equal to 0.8 mm. Similarly to fused deposition modeling (FDM) 3D printing, variations on this distance affect signifi-cantly the quality of the printed surface and errors of 2 mm would lead to unacceptable issues. Accordingly, as depicted in Figure 5, an Acuity AR-100 distance sensor is used to monitor the distance between the extruder and the skin. This measurement is used by an external controller to move the extruder along its axis n in order to correct the distance to the skin surface. More precisely, the applied external control strategy consists in a simple proportional correction, with a tool translational velocity given by

ṗec = n k lh (h m -h d ), (22)
where h m , h d and k lh are the measured layer height, desired layer height and a scalar positive gain, respectively. The joint velocities qec can be obtained with a classic redundancy resolution method [START_REF] Siciliano | Kinematic control of redundant robot manipulators: A tutorial[END_REF]. One may note that the external controller [START_REF] Santos | Model Predictive Tracking Control of Cable-Driven Parallel Robots : From Concept to Real-Time Validation[END_REF] does not have any feed-forward term and is based uniquely on the tracking error h m -h d . In spite of the simplicity of this controller, h m converges asymptotically to a constant h d if the positioning errors related to n normal to the surface are negligible. A feed-forward term is not necessary to obtain stability because a precise positioning of n results in translational velocities computed by the PPFC scheme that are tangential to the surface. This matter is discussed in detail in [START_REF] Santos | Supplementary Materials -Decoupled Model Predictive Control for Path Following on Complex Surfaces[END_REF].

In addition to the external controller [START_REF] Santos | Model Predictive Tracking Control of Cable-Driven Parallel Robots : From Concept to Real-Time Validation[END_REF], the remaining tool DoFs are controlled with the predictive control scheme proposed in Section III-B. As a result, the desired path is followed with the introduced PPFC, while the external controller corrects the distance between the extruder and the skin surface. Since the controlled robotic arm has n = 7 DoFs and the positioning of the extruder given by the tuple {p, n} represents 5 DoFs, it is worth noting that the proposed control strategy implicitly performs the corresponding redundancy resolution. Table I summarizes the used controller parameters in accordance with the notation introduced in Section III-B. The algorithm programmed in C++ was ran with an Intel(R) Core(TM) i7-7820HQ CPU @ 2.90GHz.

The desired print path is determined based on the three dimensional geometry of the skin obtained with the RGB-D camera. The printed surface covers ≈ 800 mm 2 , and leads to normal vectors n d,k with maximal variation greater than 30 • . Since the desired path consisting of points P d has significant errors due to the camera calibration, the comparison between the original desired path and the printed one P would be meaningless. A printing performed with null error with respect to P d would fail to deposit the substitute skin due to the inappropriate distance between the extruder and the skin. Therefore, in order to compare the desired and the printed path, a rigid transformation d T r is applied to the point cloud determined by P d . The rigid transformation d T r is obtained applying iterative closest point (ICP, [START_REF] Besl | A method for registration of 3-D shapes[END_REF]) in order to register the point clouds P and P d . The registered desired path P r and the measured printed path P are depicted in Figure 6. Regarding the orientation of the extruder, Figure 7 shows the desired N d and measured N orientations along with the corresponding angular error. A summary of the data related to the errors is presented in Table II. The sets P and N are obtained based on the forward kinematics. It is important to highlight that the errors presented in Table II are consistent with the positioning precision obtained with a Panda Franka Emika robotic arm. This manipulator has a pose repeatability of 0.1 mm and path deviation of 1.25 mm (based on the ISO 9283 standard). Additionally, the final printed surface depicted in Figure 8 successfully covers the skin. Regarding the computational burden, the average computing time involved in the solution of ( 14) is 1.05 ms with a standard deviation of 1.06 ms. Recalling that this solution involves the numerical optimization of the nonlinear problem ( 14) taking as arguments 56 variables (h p × n = 8 × 7), the computational efficiency of the proposed scheme is considered satisfactory. It is worth noting that the application of the analytical derivatives described in Section IV-B reduces ≈ 70% of the computation time necessary to solve [START_REF] Karayiannidis | Force/position tracking for a robotic manipulator in compliant contact with a surface using neuro-adaptive control[END_REF].

VI. CONCLUSIONS

The present letter introduced a PPFC scheme with a decoupled translational DoF aligned with the tool axis. The proposed strategy was conceived for the path following on surfaces with geometric uncertainties. Since the design of the external controller responsible for the management of the decoupled DoF is independent of the proposed PPFC, the behavior of the decoupled DoF can be adapted for different applications, e.g. robotic surface finishing and 3D printing on complex surfaces. A numerical algorithm was proposed for the solution of the corresponding optimal control problem. Experimental results in the context of robotic 3D bio-printing led to satisfying precision and computation time. Future works should analyze the stability of the obtained closed-loop system considering suboptimal control inputs and different optimization algorithms.

Fig. 1 .

 1 Fig. 1. Illustration of an application of the proposed strategy: 3D bio-printing in which the PPFC controls 5 DoFs of the tool while an external controller manages the distance between the extruder and the skin surface.

Fig. 2 .

 2 Fig. 2. Main notation used in the design of the PPFC.

Fig. 3 .

 3 Fig. 3. Illustration of the steps involved in Algorithm 1.

Algorithm 2 4 : 5 :

 245 Overall numerical algorithm to solve (14) Input: x 0 , u 0 and initial guess µ p Output: µ * 1: µ * ← µ p 2: loop 3:Update K for µ * using (2) and Algorithm 1 Compute µ * as the solution of (14) for constant K if ∥µ * -µ p ∥ < ϵ then

Fig. 5 .

 5 Fig. 5. Experimental set-up used to test robotic 3D bio-printing.

10 10Fig. 6 .

 106 Fig. 6. Experimental results: comparison between desired and printed path.

Fig. 7 .

 7 Fig. 7. Experimental results: comparison between desired and actual tool orientation.

Fig. 8 .

 8 Fig. 8. Experimental results: final printed surface.