Simona Manasra

Andrey V Kajava
email: andrey.kajava@crbm.cnrs.fr

Why does the first protein repeat often become the only one?

Keywords: protein repeats, symmetry, 3D structure, classification, swapping, evolution

Proteins with two similar motifs in tandem are one of the most common cases of tandem repeat proteins. The question arises: why is the first emerged repeat frequently fixed in the process of evolution, despite the ample opportunities to continue its multiplication at the DNA level? To answer this question, we systematically analyzed the structure and function of these proteins. Our analysis showed that, in the vast majority of cases, the structural repetitive units have a two-fold (C2) internal symmetry. These closed structures provide an internal structural limitation for the subsequent growth of the repeat number. Frequently, the units "swap" their secondary structure elements with each other. Moreover, the duplicated domains, in contrast to other tandem repeat proteins, form binding sites for small molecules around the axis of C2 symmetry. Thus, the closure of the C2 structures and the emergence of new functional sites around the axis of C2 symmetry provide plausible explanations for why a repeat, once appeared, becomes fixed in the evolutionary process. We have placed these structures within the general structural classification of tandem repeat proteins, classifying them as either Class IV or V depending on the size of the repetitive unit.

Introduction

Proteins contain a large portion of sequences with arrays of repeats that are directly adjacent to each other [START_REF] Heringa | Detection of internal repeats: How common are they?[END_REF][START_REF] Marcotte | A census of protein repeats[END_REF][START_REF] Kajava | Review: Proteins with repeated sequence -Structural prediction and modeling[END_REF][START_REF] Andrade | Protein repeats: structures, functions, and evolution[END_REF]. These tandem repeats (TRs) are ranging from the repetition of a single amino acid to motifs of 150 or more residues. Approximately half of TR regions are naturally unfolded [START_REF] Tompa | Intrinsically unstructured proteins evolve by repeat expansion[END_REF][START_REF] Simon | Tandem and cryptic amino acid repeats accumulate in disordered regions of proteins[END_REF][START_REF] Jorda | Protein tandem repeats -the more perfect, the less structured[END_REF]. The other half of them folds in a variety of either elongated or ring-like structures having a continuum of shapes, which contrast sharply with globular shapes of non-repetitive proteins. Despite a great diversity of shapes, the structures of TRs can be well-described in a synthetic manner by using a classification based on the length of their repetitive units [START_REF] Kajava | Tandem repeats in proteins: From sequence to structure[END_REF]. In this classification, the protein structures can be subdivided into five classes: (1) crystalline aggregates formed by repeats of 1 or 2 residues, (2) fibrous oligomers stabilized by inter-chain interactions with 3-7 residue repeats, (3) elongated structures with the repeat length of 5-40 residues such as, for example, solenoid folds, (4) closed (ring-like) structures consisting of 30-60 residues repeats, (5) beads on a string structures with a typical size of repeats over 50 residues, which are large enough to fold independently into stable domains. The non-globular 3D structures of TR proteins with the unusually large solvent-accessible surface area allow them to carry out specific biological functions forming, for example, building blocks of supramolecular structures, binders of periodic surfaces and molecules, or bind ligands of different types and sizes [START_REF] Andrade | Protein repeats: structures, functions, and evolution[END_REF][START_REF] Kajava | Tandem repeats in proteins: From sequence to structure[END_REF], Paladin et al.,, 2021). From an evolutionary standpoint, TRs have several specific features. Extension or shortening of the repeating sequences is known to be the result of unequal crossing-over between TR sequences during meiosis occurring between misaligned chromatids on homologous chromosomes [START_REF] Buard | Complex recombination events at the hypermutable minisatellite CEB1 (D2S90)[END_REF][START_REF] Andrade | Protein repeats: structures, functions, and evolution[END_REF]. This error-prone process occurs far more frequently than the background rate of point mutations and represents the basis for quicker evolution of the repetitive sequences in comparison with nonrepetitive ones. As a result, the number of repeats can vary even between orthologues [START_REF] Saupe | A gene responsible for vegetative incompatibility in the fungus Podospora anserina encodes a protein with a GTP-binding motif and Gβ homologous domain[END_REF]. Selective functional advantage of multiple repeats results in these TRs being fixed among populations. The critical event for the initiation of TR sequences is the first intragenic duplication, which sets the stage for unequal crossing-over and further functional specialization [START_REF] Andrade | Protein repeats: structures, functions, and evolution[END_REF]. The emergence of the first repeat is a rare event in comparison with the multiplication/reduction of already existing TRs. Its mechanism is still poorly understood. At the same time, proteins with only two domains in tandem are one of the most common cases representing over 20% among TR proteins with the known 3D structure (see Results). The question arises, why the first emerged repeat, being able to either multiply or regress to a single copy domain, nevertheless often becomes fixed in the process of evolution. Since this phenomenon cannot be explained at the DNA level, in this work, we systematically analyzed the structure and function of the proteins with only two repetitive units in order to answer this question. It is important to note that this significant class of proteins has not yet received sufficient attention from the scientific community. The primary focus of researchers has been on proteins containing three or more repeats, and the majority of bioinformatics tools, analyses, and classification methods have been developed for these multi-repeat proteins [START_REF] Kajava | Review: Proteins with repeated sequence -Structural prediction and modeling[END_REF][START_REF] Paladin | RepeatsDB in 2021: Improved data and extended classification for protein tandem repeat structures[END_REF]. To fill this gap, we developed bioinformatics protocols to select the known 3D structures with only 2 repetitive units in tandem (abbreviated here as 2RUT) and to describe their symmetries. This allowed us to classify these structures and described their typical functions. The results of this work shed light on why the first repeat frequently may become evolutionarily conserved.

Methods

Selection of protein structures containing 2RUT.

We used the PDB Release of March 2021 for the analysis [START_REF] Berman | The Protein Data Bank[END_REF]. Since the PDB contains the 3D structures of many identical or similar proteins, it was reduced using the CD-HIT program [START_REF] Huang | CD-HIT Suite: A web server for clustering and comparing biological sequences[END_REF] to have only one representative structure for sequences with more than 80% identity. Then, we removed from this dataset structures consisting of less than 60 amino acid residues assuming that most of the repeats with less than 30 amino acid residues are likely to have only one element of the secondary structure, for example, an -helix. As a result, our non-redundant set had 46595 protein structures. To select a set of proteins with two similar structural motifs arranged in tandem, we ran TAPO program [START_REF] Do | TAPO: A combined method for the identification of tandem repeats in protein structures[END_REF] against these 46595 protein structures. The TAPO program searches for repetitive structural motifs in a protein by calculating different features such as RMSDs of its fragments, by the periodicity of strings generated by conformational alphabets, residue contact maps, and arrangements of vectors of secondary structure elements. Then, the dataset of proteins with 2RUTs selected by TAPO was reduced by the CD-HIT to have one representative structure for sequences with more than 60% identity. TAPO found 4717 protein structures with 2RUT. Further, these structures were annotated using Pfam and CATH databases [START_REF] Mistry | Pfam: The protein families database in 2021[END_REF][START_REF] Sillitoe | CATH: Increased structural coverage of functional space[END_REF] and only one representative structure was left if several of them had the same CATH and Pfam ID. The detection of Pfams was made using HMMER search (http://hmmer.org) against the Pfam database. Some of the selected proteins had three and more repeats in tandem of the complete sequences but their 3D structures in PDB had two units because of the shorter fragments used for crystallization or NMR studies. Therefore, we also checked if the proteins with a given Pfam ID have domain architectures with more than two repetitive motifs. For this purpose, we calculated a parameter «Pfam2ratio»*, which is the ratio between the number of times the domain occurs in the sequence more than twice in a row and the number of all sequences where the domain occurs.

Then, we manually examined the structures with values of «Pfam2ratio» more than 10% and removed ones, which contained more than two repetitive units. In addition, our manual analysis of several randomly selected structures showed that a number of them have Rossmann fold with several α/β units. They were erroneously included in our dataset because TAPO detected only two α/β units from several ones in their structures. Therefore, we deleted the structures defined as Rossmann Fold by CATH, which had repeat unit length < 85 residues and repeat cover < 0.5. After the application of all mentioned filters, we obtained 1198 unique protein structures containing 2RUT (Supplementary Data 1).

Methods used to obtain characteristics of the selected structures.

The identified structures have the following information generated by TAPO: PDB code, ID of the chain, positions of the beginning and end of each 3D repeat, the scores evaluated the probability of correct detection of the repeats. We also obtained the secondary structure of the repeats by using the DSSP algorithm from the Bio.PDB module (DSSP of the Biopython package). Internal symmetry of the structures was determined by using CE-Symm program [START_REF] Bliven | Analyzing the symmetrical arrangement of structural repeats in proteins with CE-Symm[END_REF]. In addition, the C2 symmetry was detected by our own algorithm described in the section 3.3. Detection of swapping between the repetitive elements was made by using an algorithm developed in this work (see section 3.5). As a result, all structures from our dataset were characterized by the length of the repeats, secondary structure, fraction of the chain occupied by the repeat, domain classification from CATH and Pfam databases, frequency of TR to have more than two Pfams, type of the internal symmetry, the presence/absence of swapping (Supplementary Data 1). These parameters were further used to classify and cluster structures by the k-means method.

Clustering of the structures

Clustering was carried out separately for proteins with repeat lengths of less than and more than 55 amino acid residues because the 2RUTs with repeat lengths of less than 55 residues frequently form one common structural domain while the ones of more than 55 residues are able to form two structural domains. The clustering was made by using the KMeans module of the scikit-learn library for the Python programming language [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]. The method is based on minimizing the total square deviation of cluster points in n-dimensional (n features) space from cluster centers. The number of clusters was chosen based on the dependence of inertia on the possible number of clusters calculated using the KMeans module. Inertia, also called the within cluster sum of squares, is the principal measure of clustering performance, which measures how far apart the points are in each cluster. We also clustered the structures with Foldseek all-against-all structural alignment method [START_REF] Kempen | Crystal structure of yeast allantoicase reveals a repeated jelly roll motif[END_REF], which is making a sequence from a structure with 3D-interactions (3Di) alphabet for the further alignment.

Results and discussion

Distribution and characteristics of 2RUTs in proteins

To estimate the prevalence of the proteins with 2RUT, first, we analyzed a non-redundant protein set (sequence identity lower than 80%) from the PDB by a pipeline MetaRepeatFinder (MRF) (https://bioinfo. crbm.cnrs.fr/index.php?route=tools&tool=15 [START_REF] Richard | TRDistiller: A rapid filter for enrichment of sequence datasets with proteins containing tandem repeats[END_REF]) developed to detect TRs in protein sequence. Proteins with the repeat length of more than 30 residues were considered. The analysis showed that the cases with 2RUTs are one of the most occurring and second only to TRs with three units (Fig 1A). It is known that the PDB has a bias towards small globular structures amenable for the crystallization and NMR studies in contrast to the large elongated structures of TR-containing proteins. Hence, for a more accurate estimation of the prevalence of 2RUT proteins, it is advisable to employ a protein set derived from wellestablished proteomes. In line with this objective, we further analyzed the human proteome using MRF. The analysis revealed a comparable distribution in the occurrence of proteins with varying numbers of repeats (Fig S1 in Supplementary Data 2). It is known that structure-based methods for detecting tandem repeats generally exhibit higher sensitivity compared to sequence-based methods alone. Therefore, to compare the occurrence of 2RUT structures with those that contain more repetitive elements, we also used the results of the TAPO run against the same non-redundant set from the PDB (Fig 1B). The analysis revealed that the 2RUT structures vastly over numbers the other repetitive structures. This once again proved the importance of studying these structures. [START_REF] Richard | TRDistiller: A rapid filter for enrichment of sequence datasets with proteins containing tandem repeats[END_REF]), and only proteins with the repeat length of more than 30 residues were considered. To avoid redundancy, proteins having a level of identity higher than 80% were clustered by CD-HIT (Y. [START_REF] Huang | CD-HIT Suite: A web server for clustering and comparing biological sequences[END_REF] and represented by only one protein. (B) Occurrence of TRs found in the non-redundant set of the PDB structures by TAPO. Only structures having the repeat length of more than 30 residues were considered. The results show that 2RUT proteins represent 48% of all TRs.

Following the protocol described in Methods, we obtained 1198 structures containing two repetitive units (Supplementary Data 1) with 60% sequence similarity cutoff and unique PFAM and CATH entries. Our survey of these structures suggests that, in general, they can be subdivided into two groups. The structures with the repeat length less than approximately 55 residues fold into one structured domain, and ones with longer repeats usually have two domains each corresponding to one repeat. The analysis of the occurrence of the structures depending on the repeat length shows that 2RUT structures are especially frequent in the range of repeats between 30 and 55 residues (Fig 2). 2RUT structures containing repeats of over 200 residues are rare. These structures may contain several structural domains in one repetitive unit. The analysis of the secondary structures shows the prevalence of α-helices in the structures with repeats less than 55 residues (=52.0%; /=26.7%; =21.0%) and α/β structures (=28.1%; /=53.6%; =18.3%) in cases with repeats more than 55 residues (an insert in Fig. 2). This distinction is due to the frequent presence of four-helix bundles in the 2RUT structures with repeats of less than 55 residues.

Internal symmetry of 2RUT structures

Protein structures built of repetitive structural units frequently have an internal symmetric arrangement, which is associated with a specific biological function and allows us to trace back evolutionary events linked to the formation of this structure [START_REF] Goodsell | Structural symmetry and protein function[END_REF][START_REF] Bliven | Analyzing the symmetrical arrangement of structural repeats in proteins with CE-Symm[END_REF]. The repetitive structures can be subdivided into ones having closed and open symmetries. The closed symmetries of proteins can be further classified as cyclic (Cn), generated by a single n-fold rotational operator; dihedral (Dn), which requires an n-fold rotation and n-perpendicular 2-fold operators [START_REF] Bliven | Analyzing the symmetrical arrangement of structural repeats in proteins with CE-Symm[END_REF]. Typically, repetitive structures with closed symmetry have a strict limit on the number of repeats. Structures with the most common open symmetries have helical (H), translational (R), and superhelical (SH) arrangements of the repetitive units. In principle, these structures can have unlimited number of repeats. Since the structures with one repeat have two repetitive units it was interesting to know whether they are arranged symmetrically or not and if yes, what kind of symmetry they have. Several symmetry detection algorithms have been developed (reviewed in [START_REF] Myers-Turnbull | Systematic detection of internal symmetry in proteins using CEsymm[END_REF]. In this work, we used CE-Symm program (version 2.22), which was tested and performed well on the benchmark sets and can be downloaded locally for large-scale analyses [START_REF] Myers-Turnbull | Systematic detection of internal symmetry in proteins using CEsymm[END_REF][START_REF] Bliven | Analyzing the symmetrical arrangement of structural repeats in proteins with CE-Symm[END_REF]. The CE-Symm was run against our set of 2RUT proteins analyzing only those parts of the structures that contain two repetitive units previously detected by TAPO. Among 1198 protein structures, we identified 583 with C1, 528 with C2, 70 with R, 15 with D2, 1 with H, and 1 with C3 symmetry (Supplementary Data 1). Our manual analysis revealed a larger number of structures containing C2 symmetry, which were not detected by CE-Symm program. Therefore, we complemented the CE-Symm results with ones obtained by a new approach developed in this work (see the next section).

Method for detection of C2 symmetry based on contacts between two repetitive units.

In the C2 arrangement, inter-repeat contacts between structurally equivalent positions exist in pairs n-m' and n'-m (Fig 3 ), where n, m and n', m' are residues at n-and m-positions of the first and second repetitive unit, respectively. Since the repetitive units are often not identical in terms of sequence and structure, the equivalence of the residue positions was assessed using structural superposition of two domains. Our algorithm was based on the analysis of inter-unit contacts and if the analyzed structure had at least one pair n-m'; n'-m, it was considered to have C2 symmetry. The benchmark of this method showed that the best detection of C2 was obtained when we allowed an error ±3 in the assignment of equivalent residue positions. The program was written in Python. Input data includes PDB entry with chain ID and numbers of the first and last residues of the repeats. The repetitive structural units are superposed by using TM-align program [START_REF] Zhang | TM-align: A protein structure alignment algorithm based on the TM-score[END_REF]. The inter-unit contacts were detected by using Contact Map Explorer based on Python library MDTraj (https://contactmap.readthedocs.io/en/latest/index.html). The contact was counted when the distance between two Cβ atoms (or Cα in the case of glycine) was less than 1.0 nm. Our program "C2_contacts" is available at https://github.com/mlkndt/c2 3.4 Prevalence of C2 symmetry in 2RUT structures. CE-Symm program detects the C2 symmetry in 582 structures of 1198 proteins from the dataset (Supplementary Data 1). Of these 582 structures, our algorithm found the C2 symmetry in 557 structures, that is, in 95.7% of cases. Additionally, the new algorithm detected the C2 symmetry in 249 other structures. Our manual examination of these structures showed that indeed the vast majority of them have C2 symmetry. By summing up all structures having C2 in accordance with either CE-Symm program or our algorithm we got 831 structures, which represents 69,4% (Fig 4). The remaining structures were detected as asymmetrical or having another type of symmetry. The prevalence of the 2RUT structures among all TR structures taken together with the dominance of the C2 symmetry in the 2RUT structures agrees well with the previous observation that the C2 symmetry is the most frequent in TR structures [START_REF] Myers-Turnbull | Systematic detection of internal symmetry in proteins using CEsymm[END_REF]. The C2 symmetry prevalence for 2RUT structures (Fig 4) is especially noticeable in the structures with the repeat length less than 55 residues, which usually fold into one structural domain. The preservation of the 2RUT structures can be attributed to the closed C2 symmetry, which imposes structural constraints on the expansion of repeats.

Identification of swapping in 2RUT structures

The 2RUT protein structures with C2 symmetry, like no other structure, creates ideal conditions to exchange parts of the repetitive units with each other, in another word to "swap". Indeed, usually, elements of the secondary structure of a protein domain, being under evolutionary pressure for a long period of time, are well-fitted each other within the 3D structure. When two such units of 2RUT structure interact with each other "face-to-face" in C2 arrangement, a secondary structure element(s) of one unit, which do not fold back to their own domain, can easily find the similar well-fitted nest in the other unit (Fig 5). To a large extent, the tendency of a protein to undergo swapping is governed by thermodynamic factors. In the case of swapped oligomers, the entropic component of the free energy typically tends to disfavor the swapped state in comparison to the monomeric state [START_REF] Schlunegger | Oligomer formation by 3D domain swapping: A model for protein assembly and misassembly[END_REF]. In contrast, the 2RUT structures have a covalent bond linkage between the two repetitive units, resulting in a significant reduction of the entropic effect and promoting swapping. Moreover, this covalent linkage diminishes the influence of protein concentration on swapping, thereby increasing the likelihood that the observed swapping is genuine and not a structural artifact caused by crystallization. Therefore, it was interesting to evaluate the frequency of swapping in these structures.

To identify the "swapping" we developed an algorithm, which uses the residue-residue contact maps (Fig 5). The algorithm was slightly different for the structures with repeat length less and more than 55 residues. In the latter case, of two structural domains, the residues of the swapping part do not have contacts with the residues of their own repeat, but they have contacts with the residues of another repeat. As an example, we show in Figure 5 the residue-residue contact map of the L-arabinose-binding protein structure (pdb code 1abe:A). Residues of the last α-helix (positions 104-124) of the first repeat contact only residues of the second repeat, and vice versa. This situation is clearly visible if to follow a vertical lane at positions 104-124 on the contact map (Fig 5). The lower part of the lane does not have any dots (contacts) while the upper part of the second repeat is full of dots. One can see the opposite pattern in the lane 251-271 corresponding to the equivalent α-helix of the second repeat. Thus, our algorithm detected residue contacts, differentiating between those made within the same repeat and those involving residues from another repeat. A 2RUT structure considered to have swapping if one of its repeats contains more than 15 residues that exclusively form contacts with the other repeat. When the repeat length is less than 55 residues, two repetitive units form a single structural domain and the residues have contacts with both their own and neighboring repeats. As illustrated on the contact map of one such domain formed by two repetitive units (pdb code 1j27:A), the inter-repeat contacts are concentrated in the top left and bottom right areas. The algorithm detects well swapping when more than 10 residues from one repeat have contacts with the other repeat. Our analysis revealed that 2RUT structures with C2 symmetry have a very high percentage of the swapping. Among 832 structures with C2 symmetry, 269 have swapping that represent 32%. This tendency is especially noticeable for two-fold structures with repeats less than 55 residues where swapping is presented in 36% of structures. For comparison, 367 2RUT structures, which do not have C2 symmetry have only 9% of cases with of the swapping.

Clustering of 2RUT structures

Although our dataset was filtered to remove structures with high sequence identity, a significant number of similar structures with the same CATH or Pfam IDs remain. Therefore, for a more synthetic representation of the structures we performed their clustering by using variables from Supplementary Data 1: mean length of repeats, type of internal symmetry, secondary structure, fraction of protein occupied by repetitions, existence of swapping, CAT IDs (topology level).

For the clustering, we used the k-means method, which was implemented in klearn.cluster.Kmeans module from the scikit-learn Python library [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]. The clustering was done separately for two groups of proteins with the repeat length less and more than 55 residues. As a result, all structures were grouped into 100 clusters, 40 and 60 of which represent repeats less and more than 55 residues, respectively. The clusters are shown in Supplementary Data 2 (Fig S2). The code for its implementation, which also includes preprocessing and data normalization is presented at https://github.com/mlkndt/c2. We also performed clustering of our dataset by Foldseek software, which enables fast and sensitive comparisons of large sets of protein structures [START_REF] Kempen | Crystal structure of yeast allantoicase reveals a repeated jelly roll motif[END_REF]. The clustering was made with a threshold of 80% sequence coverage by the alignment and generated 252 clusters with similar structures (Supplementary Data 1). The results of this clustering can be used in a more detailed structural classification of proteins with TRs.

Classification of 2RUT structures

Although 2RUT structures are the most common among the structures with TRs, they are not included in the current classification [START_REF] Kajava | Tandem repeats in proteins: From sequence to structure[END_REF][START_REF] Paladin | RepeatsDB in 2021: Improved data and extended classification for protein tandem repeat structures[END_REF]. The classification is based on the length of the repeats, which is related to the ability of a single repetitive unit to fold into a stable structure. The shorter the repeat length the more of the repeats in the tandem and more chains associated together are needed to stabilize the structure. For example, TRs with the length between 1 and 30 residues require one another to maintain the structure in the form of crystallite (Class I), oligomers (Class II), or elongated structures with the inter-repeat stabilizing interactions (Class III). Although structures from Class IV consist of relatively long repeats (30-60 residues), individually, these repeats are at the borderline in terms of the structure stability and become stable in the context of closed (ring-like) structures. Finally, proteins with the size of repeats over 50 residues (Class V) are large enough to fold independently into stable domains. The question arises: where is the place of the 2RUT structures in this classification? The 2RUT structures have a typical repeat length between 30 and 200 residues. Our analysis of the know 2RUT structures revealed that when their repeat size is between 30 and 55 residues they frequently form one common structural domain with the closed C2 symmetry (Fig. 6). By the repeat length and the necessity of stabilizing interactions to fold into the single structured domains, these 2RUT proteins resemble Class IV structures. The other structures are composed of repeats over 55 residues, which, similarly to the repeats of Class V proteins, are able to form structural domains independently. Thus, we conclude that the structures with units between 30 and 55 residues belong to Class IV structures, while the ones with units of more than 55 residues are a part of Class V structures. Structures with the repeat length over 55 residues: (D) mainly α-helical, conserved histidine α-helical domain (PDB entry 6qv5:A), (E) mainly β-structural, LU domain (Ly-6 antigen/uPAR) of three-finger protein superfamily (PDB entry 6iom:A), (F) α/β, Imp4-like protein (PDB entry 1w94:A), (G) a structure with very long repeat of 480 residues, which contains several structural domains in each repeat, acriflavine resistance protein (PDB entry 4k0e:A), (H) Box repeat structure of heterotrimeric PCNA sliding clamp (PDB entry 2hik:A). In this ring-like structure, three molecules with C2 symmetric beads interact with one another. The position of one of the twofold axes within the bead is depicted by a black, lens-shaped symbol on the left side of the ring.

Functions of 2RUT structures

The 2RUT structures also have their own special functions. Unlike many TR proteins, with extended surfaces well suited for protein-to-protein interactions [START_REF] Andrade | Protein repeats: structures, functions, and evolution[END_REF][START_REF] Peterson | A domain shared by the Polycomb group proteins Scm and ph mediates heterotypic and homotypic interactions[END_REF][START_REF] Bork | The immunoglobulin fold: Structural classification, sequence patterns and common core[END_REF], 2RUT structures usually do not have such a function. In accordance with Pfam annotation, they are often associated with binding to small molecules. Typical examples are ferredoxin, which contains an iron-sulfur cluster and carries out electron transport [START_REF] Otaka | Examination of protein sequence homologies: IV. Twenty-seven bacterial ferredoxins[END_REF], and heavy metal-associated domains, containing two conserved cysteines involved in binding [START_REF] Gitschier | Solution structure of the fourth metal-binding domain from the Menkes copper-transporting ATPase[END_REF]. Besides metals and their transport, structures with two domains can participate in the binding of other small molecules [START_REF] Anantharaman | The SHS2 module is a common structural theme in functionally diverse protein groups, like Rpb7p, FtsA, GyrI, and MTH1598/Tm1083 superfamilies[END_REF], including fatty acid molecules [START_REF] Schulze-Gahmen | Crystal structure of a hypothetical protein, TM841 of Thermotoga maritima, reveals its function as a fatty acid-binding protein[END_REF]. Among the 2RUT structures, there are many enzymes whose active site is located between two repeating domains [START_REF] Akey | Crystal Structures of Dehydratase Domains from the Curacin Polyketide Biosynthetic Pathway[END_REF][START_REF] Chisuga | Structural analysis of the dualfunction thioesterase SAV606 unravels the mechanism of Michael addition of glycine to an α,β-unsaturated thioester[END_REF][START_REF] Haynes | Structures of nitroreductase in three states. Effects of inhibitor binding and reduction[END_REF][START_REF] Hofmann | Structure and mechanism of activity of the cyclic phosphodiesterase of Appr>p, a product of the tRNA splicing reaction[END_REF][START_REF] Huang | Crystal structure of β-ketoacyl-acyl carrier protein synthase II from E.coli reveals the molecular architecture of condensing enzymes[END_REF][START_REF] Iyer | Evolution of the deaminase fold and multiple origins of eukaryotic editing and mutagenic nucleic acid deaminases from bacterial toxin systems[END_REF][START_REF] Lloyd | Refinement of Haemophilus influenzae diaminopimelic acid epimerase (DapF) at 1.75 Å resolution suggests a mechanism for stereocontrol during catalysis[END_REF]. For example, the Allantoicase domain, involved in purine cleavage, contains two jelly-roll motifs, each of which has conserved pockets constituting a common active site (Leulliot et al., 2004). On the other hand, repeats do not always have the same effect on enzyme activity. The M16 family peptidase consists of two structurally related domains, but only one of them is an active peptidase, while the other remains inactive. However, the C2 symmetrical structure may help the enzyme to "clamp" the substrate [START_REF] Taylor | Crystal structures of mitochondrial processing peptidase reveal the mode for specific cleavage of import signal sequences[END_REF]. In the case of phosphofructokinases, involved in the regulation of the glycolytic pathway, one domain of the structure, in addition to the binding site, contains an allosteric site that affects the activity of the enzyme [START_REF] Shirakihara | Crystal structure of the complex of phosphofructokinase from Escherichia coli with its reaction products[END_REF]. Occasionally, 2RUT structures exhibit functions that are closely tied to their shapes. For example, the Tli4 family domain, which is part of the effector-immunity Tle-Tli pair in some bacteria, has two domains that together have shape similar to a crab claw. It is used to grasp the Tle domain to prevent its activation [START_REF] Lu | The structural basis of the Tle4-Tli4 complex reveals the self-protection mechanism of H2-T6SS in Pseudomonas aeruginosa[END_REF]. In the proteins of the TBP family that are required to initiate the transcription process by RNA polymerases the two units form a saddle-like shape, which interacts with transcription initiation factors and regulatory proteins after TATA-box binding [START_REF] Nikolov | Crystal structure of TFIID TATA-box binding protein[END_REF]. [START_REF] Hofmann | Structure and mechanism of activity of the cyclic phosphodiesterase of Appr>p, a product of the tRNA splicing reaction[END_REF]. (B) Diaminopimelate epimerase (PDB entry 2gke:A), active-site cleft is formed from both domains [START_REF] Lloyd | Refinement of Haemophilus influenzae diaminopimelic acid epimerase (DapF) at 1.75 Å resolution suggests a mechanism for stereocontrol during catalysis[END_REF]. ConSurf web server [START_REF] Ashkenazy | ConSurf 2016: an improved methodology to estimate and visualize evolutionary conservation in macromolecules[END_REF] was used for estimation of the evolutionary conservation of residues.

In some cases, C2 structures with two domains can also be incorporated into multi-repeat structures with a "beads-on-a-string" architecture, where the C2 structures represents the beads. The ample binding surfaces of these proteins make them particularly adept at engaging with large macromolecules [START_REF] Whitley | Burkholderia oklahomensis agglutinin is a canonical two-domain OAA-family lectin: Structures, carbohydrate binding and anti-HIV activity[END_REF][START_REF] Arrías | The repetitive structure of DNA clamps: An overlooked protein tandem repeat[END_REF].

Conclusions

To investigate still overlooked protein structures with two repetitive units, we built a nonredundant dataset of the known 3D structures. It allowed us to analyze these structures for the occurrence depending on repeat length, internal symmetry, swapping, annotation in Pfam and CATH databases, and functional roles. The analysis revealed that such structures are the most widespread among the TR proteins emphasizing their particular importance. Depending on the repeat length, all 2RUT structures can be subdivided into two groups: (i), units of less than about 55 residues assemble into one structural domain (Fig. 6A-C) and (ii) the repeats over approximately 55 residues fold in two structural domains (Fig. 6 D-F). The repeats of over 200 residues are rare but if one does occur, it often contains multiple structural domains within that repeat (Fig 6G). The vast majority of these structures exhibit two-fold symmetry (C2) and this finding provides insight into why the first repeat, which could theoretically continue to multiply at the DNA level, becomes fixed in the process of evolution. Indeed, the closed C2 arrangement provides a structural restriction for the subsequent expansion of the repeats. Due to the C2 symmetry arrangement, two domains cannot provide the necessary binding interface for a third domain, rendering it superfluous. Thus, one can suggest an evolutionary scenario where after the emergence of the second structural domain which is able to interact with the first one in the C2 structure, the third such domain will no longer be able to join the first two domains in a similar manner. Further multiplication is possible if the C2 structure with two repetitive units duplicates as a whole generating tandems with an even number of repeats. Several such structures were observed, including Box-repeat structures (Fig 6H) [START_REF] Arrías | The repetitive structure of DNA clamps: An overlooked protein tandem repeat[END_REF]. The other specific characteristic of the 2RUT structures with C2 symmetry is an extremely high percentage of swapping. Swapping enhances the fidelity of the face-to-face interaction of the two repetitive units through a kind of "hugging". The fact that a given symmetrical structure was fixed in evolution does not imply the structural but rather functional reasons when the duplication of a region within a gene imparts a new function, which is beneficial for the organism fitness. Indeed, many 2RUT proteins have very specific functions in comparison to the other TR-containing proteins. The duplicated domains with the C2 symmetry provide structural basis for the emergence of new binding sites for small molecules at and around the axe of two-fold symmetry (Fig. 7). Therefore, these proteins are often associated with binding to small molecules including metals, and their subsequent transport or chemical transformations. In some other cases, one may assume an evolutionary scenario when duplication of an ancient gene leads to a protein with two domains retaining their original function, and then, this protein undertakes a kind of subfunctionalization with one of the domains losing enzymatic activity for a new function. For example, in phosphofructokinases (PDB 3opy:B) [START_REF] Shirakihara | Crystal structure of the complex of phosphofructokinase from Escherichia coli with its reaction products[END_REF] only one domain has an enzyme activity and another one is responsible for allosteric effect of the enzyme.

It is worth mentioning that sometimes the C2 structures with two domains can also be a part of multi-repeat structures with a "beads on a string" architecture where the two domain structures represent the beads. These proteins have extended binding surfaces suitable for interactions with large macromolecules. For example, in Box-repeat proteins, two-three molecules with such C2 symmetrical beads interact with each other to form ring-like bracelets (see [START_REF] Arrías | The repetitive structure of DNA clamps: An overlooked protein tandem repeat[END_REF]. Finally, our analysis led to the classification of these structures within the general structural classification of TR proteins. We suggest that the 2RUT structures with repeats between 30 and 55 residues belong to Class IV structures, while the ones with units over 55 residues are a part of Class V structures. Given that this type of structures is the most abundant among TR proteins, inclusion of them in the RepeatsDB [START_REF] Paladin | RepeatsDB in 2021: Improved data and extended classification for protein tandem repeat structures[END_REF] would enhance the database's comprehensiveness.

Figure 1 .

 1 Figure 1. (A)Occurrence of tandem repeats with a different number of repeats in the non-redundant set of proteins from PDB. The repeats were identified by a pipeline MetaRepeatFinder (MRF) (https://bioinfo.crbm.cnrs.fr/index.php?route=tools&tool=15[START_REF] Richard | TRDistiller: A rapid filter for enrichment of sequence datasets with proteins containing tandem repeats[END_REF]), and only proteins with the repeat length of more than 30 residues were considered. To avoid redundancy, proteins having a level of identity higher than 80% were clustered by CD-HIT (Y.[START_REF] Huang | CD-HIT Suite: A web server for clustering and comparing biological sequences[END_REF] and represented by only one protein. (B) Occurrence of TRs found in the non-redundant set of the PDB

Figure 2 .

 2 Figure 2. Occurrence of structures with 2RUT depending on the repeat length. The repeat length starts from 30 residues. Structures with repeat lengths ranging between 30 and 55 are in blue, while structures with repeat lengths exceeding 55 are represented in orange. An insert with swarm plots shows the distribution of secondary structures depending on the repeat length.

Figure 3 .

 3 Figure 3. (A) A scheme of contacts in the interacting repetitive units having C2 symmetry. (B) The 3D structure of cysteine protease (Peptidase_Prp) domain with C2 symmetry (PDB entry 2idl:A). Pairs of contacts between residues 4-95 and 58-45, which in the .3D alignment have numbers 4-38' and 4'-38, where 4, 38 are positions of these residues in the first unit and 4', 38' in the second one. If protein had at least one such pair of contacts, it was considered to have C2 symmetry.

Figure 4 .

 4 Figure 4. Occurrence of 2RUT structures with C2 in comparison with the other symmetries depending on the repeat length (starting from 30 residues). C2-symmetrical structures are shown in

Figure 5 .

 5 Figure 5. Examples of the structures and their contact maps for swapping identification. Each map is divided into four squares containing intra-and inter-repeat contacts of two repetitive units. Regions with contacts of interest are shown in blue: on the left, these are regions where one repeat has contacts only with another one for the structures with the repeat length over 55 residues (PDB 1abe:A) and, on the right, these are regions where one repeat has numerous contacts with another one for structures with the repeat length less than 55 residues (PDB 1j27:A). The source code of C2_swapping program is available at https://github.com/mlkndt/c2.

Figure 6 .

 6 Figure 6. Representative 2RUT structures having C2 symmetry. Position of the twofold axis is represented by a black lens-shaped symbol. Structures with the repeat length less than 55 residues: (A) mainly α-helical, [2Fe-2S] binding domain (PDB entry 3hrd:D), (B) mainly β-structural, biotin attachment domain (PDB entry 1ghj:A), (C) α/β, heavy-metal-associated domain (PDB entry 1cc7:A). Structures with the repeat length over 55 residues: (D) mainly α-helical, conserved histidine α-helical domain (PDB entry 6qv5:A), (E) mainly β-structural, LU domain (Ly-6 antigen/uPAR) of three-finger protein superfamily (PDB entry 6iom:A), (F) α/β, Imp4-like protein (PDB entry 1w94:A), (G) a structure with very long repeat of 480 residues, which contains several structural domains in each repeat, acriflavine resistance protein (PDB entry 4k0e:A), (H) Box repeat structure of heterotrimeric PCNA sliding clamp (PDB entry 2hik:A). In this ring-like structure, three molecules with C2 symmetric beads

Figure 7 .

 7 Figure 7. Examples of C2 binding sites. (A) Cyclic phosphodiesterase-like protein domain (PDB entry 1jh7:A), where the central water-filled cavity forms the enzyme's active site[START_REF] Hofmann | Structure and mechanism of activity of the cyclic phosphodiesterase of Appr>p, a product of the tRNA splicing reaction[END_REF]. (B) Diaminopimelate epimerase (PDB entry 2gke:A), active-site cleft is formed from both domains[START_REF] Lloyd | Refinement of Haemophilus influenzae diaminopimelic acid epimerase (DapF) at 1.75 Å resolution suggests a mechanism for stereocontrol during catalysis[END_REF]. ConSurf web server[START_REF] Ashkenazy | ConSurf 2016: an improved methodology to estimate and visualize evolutionary conservation in macromolecules[END_REF] was used for estimation of the evolutionary conservation of residues.

Acknowledgments

This work was supported by REFRACT MSCA RISE project within European Union's Horizon 2020 research and innovation programme under grant agreement No 823886, as well as from ML4NGP (CA21160), supported by COST (European Cooperation in Science and Technology) under the EU Framework Programme Horizon Europe.


Reyd A Smith

Jérôme Fort

Pierre Legagneux

Olivier Chastel

Mark L Mallory

Paco Bustamante

Jóhannis Danielsen

Sveinn A Hanssen

Jón Einar Jónsson

Ellen Magnúsdóttir

Børge Moe

Charline Parenteau

Kyle J L Parkinson

Glen J Parsons

Grigori Tertitski

Oliver P Love

Do foraging ecology and contaminants interactively predict parenting hormone levels in common eider?

Keywords: stable isotopes, carbon-13, nitrogen-15, mercury, seabird, Arctic

Blood samples from 106 common eiders from 6 colonies were analyzed (δ 13 C, δ 15 N, THg, prolactin)

• A significant interaction was found between δ 13 C, δ 15 N and Hg on prolactin • In the 3-way interaction, δ 15 N positively correlated with incubating prolactin

• Low δ 13 C and high Hg had a consistent interactive impact in combination with δ 15 N • This complex interaction may have downstream implications for reproductive success

Introduction

Increasing anthropogenic activity has resulted in rapid and wide-ranging abiotic shifts in marine systems associated with climate change, including rising air and ocean temperatures [START_REF] Screen | Increasing fall-winter energy loss from the Arctic Ocean and its role in Arctic temperature amplification[END_REF][START_REF] Najafi | Attribution of Arctic temperature change to greenhouse-gas and aerosol influences[END_REF], changes in ocean circulation and salinity [START_REF] Timmermans | Surface freshening in the Arctic Ocean's Eurasian Basin: An apparent consequence of recent change in the wind-driven circulation[END_REF][START_REF] Zika | Improved estimates of water cycle change from ocean salinity: the key role of ocean warming[END_REF], and dramatic changes in sea ice abundance and distribution within polar systems [START_REF] Johannessen | Arctic climate change: observed and modelled temperature and sea-ice variability[END_REF][START_REF] Hoegh-Guldberg | The impact of climate change on the world's marine ecosystems[END_REF]. Arctic regions disproportionately exhibit these simultaneous, cumulative effects (IPCC 2019; AMAP 2021), resulting in impacts on food web relationships [START_REF] Vasseur | A mechanistic approach for modeling temperaturedependent consumer-resource dynamics[END_REF][START_REF] Frederiksen | From plankton to top predators: bottom-up control of a marine food web across four trophic levels[END_REF][START_REF] Parmesan | Ecological and evolutionary responses to recent climate change[END_REF][START_REF] Pecuchet | Novel feeding interactions amplify the impact of species redistribution on an Arctic food web[END_REF]. Alterations in the abundance and distribution of lower trophic level organisms (i.e., primary producers) have been evident, with organisms such as phytoplankton displaying advanced phenology and higher abundance relative to sea ice algae [START_REF] Frederiksen | From plankton to top predators: bottom-up control of a marine food web across four trophic levels[END_REF][START_REF] Renaut | Northward expansion and intensification of phytoplankton growth during the early ice-free season in Arctic[END_REF][START_REF] Lewis | Changes in phytoplankton concentration now drive increased Arctic Ocean primary production[END_REF]). Additionally, changes in higher trophic level foraging ecology have been associated with climate change, as prey disruptions further disturb the relationships among food webs [START_REF] Moline | High latitude changes in ice dynamics and their impact on polar marine ecosystems[END_REF][START_REF] Watt | A shift in foraging behaviour of beluga whales Delphinapterus leucas from the threatened Cumberland Sound population may reflect a changing Arctic food web[END_REF]). Thus, bottom-up and top-down trophic level alterations can result in modifications to individual and species foraging ecology and, in some cases, sizeable population declines [START_REF] Legagneux | Disentangling trophic relationships in a High Arctic tundra ecosystem through food web modeling[END_REF][START_REF] Perkins | Combined bottom-up and top-down pressures drive catastrophic population declines of Arctic skuas in Scotland[END_REF][START_REF] Durant | Contrasting effects of rising temperatures on trophic interactions in marine ecosystems[END_REF]. Seabirds are particularly susceptible to adverse outcomes resulting from changes in both directions of trophic alterations, which can influence body condition and, importantly, downstream reproductive success [START_REF] Baird | Influence of abiotic factors and prey distribution on diet and reproductive success of three seabird species in Alaska[END_REF][START_REF] Hamer | The influence of food supply on the breeding ecology of kittiwakes Rissa tridactyla in Shetland[END_REF][START_REF] Suryan | Site-specific effects on productivity of an upper trophic-level marine predator: bottom-up, top-down, and mismatch effects on reproduction in a colonial seabird[END_REF][START_REF] Perkins | Combined bottom-up and top-down pressures drive catastrophic population declines of Arctic skuas in Scotland[END_REF][START_REF] Hovinen | Relationships between isotopic ratios, body condition and breeding success in a High Arctic seabird community[END_REF]. Additionally, changing predator-prey dynamics and spatio-temporal prey abundance can alter, either positively or negatively, individual seabird exposure to environmental contaminants such as mercury (Hg) through changes in prey selection and foraging locale [START_REF] Bearhop | Variability in mercury concentrations of great skuas Catharacta skua: the influence of colony, diet and trophic status inferred from stable isotope signatures[END_REF][START_REF] Braune | Changes in food web structure alter trends of mercury uptake at two seabird colonies in the Canadian Arctic[END_REF][START_REF] Fort | Does temporal variation of mercury levels in Arctic seabirds reflect changes in global environmental contamination, or a modification of Arctic marine food web functioning?[END_REF], Tartu et al. 2022).

Arctic marine systems are contaminated with Hg via northward transportation of Hg through air and ocean circulation as well as through Arctic rivers [START_REF] Macdonald | Contaminants in the Canadian Arctic: 5 years of progress in understanding sources, occurrence and pathways[END_REF][START_REF] Sonke | Eurasian river spring flood observations support net Arctic Ocean mercury export to the atmosphere and Atlantic Ocean[END_REF][START_REF] Dimento | The impact of sea ice on the air-sea exchange of mercury in the Arctic Ocean[END_REF], and the melting of the cryosphere which releases stored Hg from years of higher contaminant output back into the system [START_REF] Rydberg | Climate driven release of carbon and mercury from permafrost mires increases mercury loading to sub-arctic lakes[END_REF][START_REF] Schuster | Permafrost stores a globally significant amount of mercury[END_REF][START_REF] Hawkings | Large subglacial source of mercury from the southwestern margin of the Greenland Ice Sheet[END_REF]). These factors have all cumulatively resulted in spatiotemporally heterogeneous environmental Hg levels in Arctic regions that influence organism Hg levels by factors including spatio-temporal location and foraging ecology [START_REF] Braune | Persistent organic pollutants and mercury in marine biota of the Canadian Arctic: an overview of spatial and temporal trends[END_REF]2015;[START_REF] Leu | Consequences of changing sea-ice cover for primary and secondary producers in the European Arctic shelf seas: timing, quantity, and quality[END_REF][START_REF] De | Temporal and spatial trends in marine carbon isotopes in the Arctic Ocean and implications for food web studies[END_REF]Albert et al. 2021(b)). Individuals foraging at higher trophic levels, usually determined using stable isotope analysis of δ 15 N (in relation to 15 N to 14 N ratios), are often correlated with elevated Hg levels resulting from the methylmercury (MeHg) biomagnification process, although individual variation outside of this correlation often exists [START_REF] Atwell | Biomagnification and bioaccumulation of mercury in an arctic marine food web: insights from stable nitrogen isotope analysis[END_REF][START_REF] Bearhop | Variability in mercury concentrations of great skuas Catharacta skua: the influence of colony, diet and trophic status inferred from stable isotope signatures[END_REF][START_REF] Dimento | The impact of sea ice on the air-sea exchange of mercury in the Arctic Ocean[END_REF]. Moreover, the influence of Arctic marine primary producers on the diets of higher order predators, often determined using ratios of carbon isotopes ( 12 C -phytoplankton enriched, 13 C -ice algae enriched; δ 13 C), can offer insight into potential Hg sources based on foraging location and sea ice abundance [START_REF] Hobson | Determination of trophic relationships within a high Arctic marine food using δ 13 C and δ 15 N analysis[END_REF][START_REF] Hobson | Sources of primary production, benthicpelagic coupling, and trophic relationships within the Northeast Water Polynya: insights from δ 13 C and δ 15 N analysis[END_REF][START_REF] Mcmahon | Benthic community response to ice algae and phytoplankton in Ny-Ålesund, Svalbard[END_REF].

MeHg is a biologically converted form of Hg that affects endocrine functioning, behaviour, and reproductive success in birds at even low, sub-lethal effects [START_REF] Wiener | Ecotoxicology of mercury[END_REF][START_REF] Whitney | Impacts of sublethal mercury exposure on birds: a detailed review[END_REF][START_REF] Evers | The effects of methylmercury on wildlife: a comprehensive review and approach for interpretation[END_REF], Chastel et al. 2022). While the Hg levels in most seabirds (95%) are considered lower risk, there have been recent calls to examine the potential impact of Hg in a cumulative stressor framework to help address Article 22 of the Minimata Convention (Chastel et al. 2022). Avian reproduction can be affected by Hg exposure through its endocrine-disrupting properties that can influence female breeding behaviour (Tan et al. 2009;[START_REF] Ackerman | Avian mercury exposure and toxicological risk across western North America: a synthesis[END_REF][START_REF] Whitney | Impacts of sublethal mercury exposure on birds: a detailed review[END_REF]Cristol 2017, Chastel et al. 2022). For example, in the beginning stages of reproduction, sub-lethal Hg has been correlated with a decreased likelihood to breed, attributed to the effect of Hg on luteinizing hormone (LH) and the subsequent release of sex hormones (i.e., estrogen; [START_REF] Tartu | To breed or not to breed: endocrine response to mercury contamination by an Arctic seabird[END_REF]2014). Importantly, exposure to Hg has been linked to a decrease in prolactin (PRL), a pituitary hormone key for parental egg and chick attachment [START_REF] El Halawani | Regulation of prolactin and its role in gallinaceous bird reproduction[END_REF][START_REF] Tartu | Corticosterone, prolactin and egg neglect behavior in relation to mercury and legacy POPs in a long-lived Antarctic bird[END_REF][START_REF] Angelier | Does prolactin mediate parental and life-history decisions in response to environmental conditions in birds? A review[END_REF][START_REF] Tartu | Mercury exposure, stress and prolactin secretion in an Arctic seabird: an experimental study[END_REF]. However, no link between Hg and PRL has been observed in other studies, indicating a potentially complex relationship and role of PRL in avian incubation [START_REF] Gilmour | Mercury as an indicator of foraging ecology but not the breeding hormone prolactin in seabirds[END_REF].

In birds, PRL secretion commences after the parent physically touches their eggs and increases following the cessation of laying and the onset of incubation behaviour [START_REF] Buntin | Neural and hormonal control of parental behavior in birds[END_REF][START_REF] Sockman | The role of prolactin in the regulation of clutch size and onset of incubation behavior in the American kestrel[END_REF]. In common eider (Somateria mollissima, hereafter eider), PRL increases at the onset of incubation and remains at a consistent level during incubation, followed by a steep increase immediately before hatching [START_REF] Criscuolo | Factors affecting plasma concentrations of prolactin in the common eider Somateria mollissima[END_REF]. Additionally, the role of PRL in clutch size determination has been previously proposed, with increased PRL signalling for a decrease in LH that stops the laying period in anticipation of incubation [START_REF] Sockman | The role of prolactin in the regulation of clutch size and onset of incubation behavior in the American kestrel[END_REF].

Disruption of PRL secretion during incubation has been proposed as a likely mechanistic link that results in lowered quality of incubation behaviour via decreased incubation attentiveness (Angelier and Chastel 2009;Angelier et al. 2009). Incubating females in worsened body condition may have lower PRL levels, thus reducing incubation quality and increasing predation risk [START_REF] Criscuolo | Factors affecting plasma concentrations of prolactin in the common eider Somateria mollissima[END_REF]. Overall, the cumulative and interactive effects of prey consumption and Hg accumulation on PRL can influence key reproductive traits such as incubation behaviour and reproductive success. These linkages are important to examine, especially given the accelerating rate of environmental shifts in Arctic marine systems.

A recent study from our group indicated that eider colonies highly varied in multidimensional niche dynamics: migratory colonies had higher δ 15 N and THg, but lower δ 13 C and sea-ice cover positively correlated with δ 15 N, δ 13 C and THg [START_REF] Smith | Environmental and life-history factors influence inter-colony multidimensional niche metrics of a breeding Arctic marine bird[END_REF]). The study also found that individual eiders with lower δ 13 C and higher trophic positions had higher THg [START_REF] Smith | Environmental and life-history factors influence inter-colony multidimensional niche metrics of a breeding Arctic marine bird[END_REF]). In the current study, we built from these results to examine the interactive effects of foraging ecology, inferred with stable isotope analysis [START_REF] Kelly | Stable isotopes of carbon and nitrogen in the study of avian and mammalian trophic ecology[END_REF], and Hg on PRL levels of incubating female common eiders. Eiders are long-lived, colonial-nesting marine ducks with dispersed breeding populations across the circumpolar and sub-Arctic. As such, they are exposed to diverse environmental conditions, which influence variation in their foraging ecology and Hg exposure (Albert et al. 2021(a); [START_REF] Smith | Environmental and life-history factors influence inter-colony multidimensional niche metrics of a breeding Arctic marine bird[END_REF]. Eider prey includes a diversity of benthic invertebrates (e.g., urchins and mussels) and pelagic macroinvertebrates (e.g., amphipods; [START_REF] Sénéchal | Do purely capital layers exist among flying birds? Evidence of exogenous contribution to arctic-nesting common eider eggs[END_REF]. Eiders forage extensively in the pre-breeding period in various ice conditions depending on the timing of spring melt to obtain an appropriate body condition needed to breed [START_REF] Jean-Gagnon | The impact of sea ice conditions on breeding decisions is modulated by body condition in an arctic partial capital breeder[END_REF]. Importantly, female eiders undergo a 24-26 day uniparental incubation fast, relying solely on their endogenous stores to fuel the incubation period [START_REF] Hanssen | Incubation start and egg size in relation to body reserves in the common eider[END_REF][START_REF] Bottitta | Effects of experimental manipulation of incubation length on behavior and body mass of common eiders in the Canadian Arctic[END_REF][START_REF] Sénéchal | Do purely capital layers exist among flying birds? Evidence of exogenous contribution to arctic-nesting common eider eggs[END_REF]. During this fasting period, contaminants including Hg increase in the blood due to the mobilization of endogenous lipid and protein stores [START_REF] Wayland | Concentrations of cadmium, mercury and selenium in common eider ducks in the eastern Canadian arctic: influence of reproductive stage[END_REF][START_REF] Bustnes | Strongly increasing blood concentrations of lipid-soluble organochlorines in high arctic common eiders during incubation fast[END_REF][START_REF] Mcpartland | Lead isotopic signatures in blood from incubating common eiders (Somateria mollissima) in the central Baltic Sea[END_REF]. These levels can go over the threshold known to be negatively correlated with PRL secretion in other seabird species [START_REF] Tartu | Corticosterone, prolactin and egg neglect behavior in relation to mercury and legacy POPs in a long-lived Antarctic bird[END_REF][START_REF] Wickham | ggplot2: Elegant Graphics for Data Analysis[END_REF]. Consequently, in combination with reduced body condition, elevated blood levels of Hg may interfere with hormones such as PRL, with the potential for downstream influences on incubation behaviour and reproductive success [START_REF] Tartu | Mercury exposure, stress and prolactin secretion in an Arctic seabird: an experimental study[END_REF]. However, intricate relationships between these and other physiological aspects have been demonstrated previously in the literature.

In common eiders, complex and integrated relationships have been shown between clutch size, body condition/mass, incubation duration, PRL and the primary stress hormone in birds, corticosterone (CORT) [START_REF] Erikstad | Clutch-size determination in precocial birds: a study of the common eider[END_REF][START_REF] Criscuolo | Body mass and clutch size may modulate prolactin and corticosterone levels in eiders[END_REF]. Clutch size and body condition are positively related, while body condition can be negatively correlated with incubation duration but positively correlated with incubation attentiveness [START_REF] Erikstad | Clutch-size determination in precocial birds: a study of the common eider[END_REF][START_REF] Erikstad | Does the cost of incubation set limits to clutch size in common eiders Somateria mollissima?[END_REF]. Meanwhile, body mass was positively linked to PRL levels, while baseline levels of CORT were negatively correlated with plasma PRL, altogether leading to reductions in overall reproductive success [START_REF] Criscuolo | Factors affecting plasma concentrations of prolactin in the common eider Somateria mollissima[END_REF]2005;[START_REF] Parmesan | Ecological and evolutionary responses to recent climate change[END_REF]. These relationships provide a background for examining the effects of foraging ecology and Hg on PRL in female eiders.

Our overall hypothesis was that the combined effects of foraging ecology and Hg contamination would affect PRL levels in incubating eiders. Specifically, we predicted that lower levels of PRL would be correlated with lower δ 15 N (lower trophic level foraging, lesser quality prey), higher δ 13 C (sea ice algae-based foraging, more restrictive foraging environment), and elevated Hg levels. Examining these relationships will allow for a more extensive determination of the indirect mechanisms by which climate change affects at-risk Arctic seabird populations.

Methods

Sample collection and physiological analysis

Seven Arctic and sub-Arctic breeding colonies were monitored in 2018 (Table 2,N=136).

Incubating eiders were captured on their nests using noose poles, bownets, or dogs (John's Island location only; [START_REF] Milton | Sex-Specific Survival of Adult Common Eiders in Nova Scotia, Canada[END_REF]. A 200-1,000 µL blood sample was collected from each female's tarsal vein, jugular vein (Grindøya and Kongsfjorden), or brachial vein (Iceland) using a 23G thinwall, 1-inch needle attached to a heparinized 1 mL syringe (10 mL syringe -Grindøya and Kongsfjorden). Samples were centrifuged within 8-12 h of collection at approximately 10,000 rpm for 10 min to obtain red blood cells (RBC) for stable isotope (isotope-ratio mass-spectrometer) and total Hg (THg; Advanced Mercury Analyser; see [START_REF] Smith | Environmental and life-history factors influence inter-colony multidimensional niche metrics of a breeding Arctic marine bird[END_REF] for further details on these methods and results). Plasma was collected from the centrifuged whole blood, frozen for storage in the field, and later sent for PRL analysis to the Centre d'Etudes Biologiques de Chizé (CEBC), France. Plasma samples were analyzed for PRL using radioimmunoassay (RIA) following the procedure detailed in [START_REF] Cherel | Seasonal and fasting-related changes in circulating gonadal steroids and prolactin in king penguins, Aptenodytes patagonicus[END_REF] and validated in common eiders in previous studies [START_REF] Criscuolo | Factors affecting plasma concentrations of prolactin in the common eider Somateria mollissima[END_REF]. The analysis consisted of duplicates with colonies split into two runs: John's Island, Iceland and Tern Island in one run (intra-assay variation = 12.01%), and Grindøya, Faroe Islands, Kongsfjorden, and Solovetskiy Islands in a second run (intra-assay variation = 9.28%). Inter-assay variation between the two runs was 14.75%. PRL was measured in 10µL, and the lowest detectable limit was 0.45 ng/mL.

As specific individual incubation stage was unable to be determined at many sites or hens, we could not correct the potential effect of incubation fasting on THg or δ 15 N levels [START_REF] Hobson | Stable-nitrogen isotope enrichment in avian tissues due to fasting and nutritional stress: implications for isotopic analyses of diet[END_REF][START_REF] Wayland | Concentrations of cadmium, mercury and selenium in common eider ducks in the eastern Canadian arctic: influence of reproductive stage[END_REF][START_REF] Smith | Environmental and life-history factors influence inter-colony multidimensional niche metrics of a breeding Arctic marine bird[END_REF]). However, using RBCs to determine isotopic and THg values minimized this effect as δ 15 N enrichment is lesser in blood cells than in plasma [START_REF] Cherel | Nutrition, physiology, and stable isotopes: new information from fasting and molting penguins[END_REF]. Further, as PRL levels are consistent throughout incubation stages in eiders prior to hatching [START_REF] Criscuolo | Factors affecting plasma concentrations of prolactin in the common eider Somateria mollissima[END_REF], we were able to conduct inter-colony and individual comparisons of PRL without correction for individual or colony incubation stage.

Statistical analysis

To examine whether PRL levels varied between colonies, we first ran an ANOVA followed by a posthoc Tukey's Honest Significant Difference (HSD) test using the R package car v.3.0-10 (Fox and Weisberg 2019). We then used compact letter display, in combination with the R package ggplot2 v.3.3.6, to visually determine significance at an alpha level of 0.05 between the seven colonies using the R package multcompView v.0.1-8 [START_REF] Graves | Package 'multcompView[END_REF][START_REF] Wickham | ggplot2: Elegant Graphics for Data Analysis[END_REF].

Following this, we calculated AICc values for 12 potential models to determine the most parsimonious model using the R package AICcmodavg v.2.3-1 (Mazerolle 2020). The model with the lowest AICc also had a significantly lower log-likelihood value, thus, was a better fit than the other models (Table 1). We then used a general linear mixed model (GLMM) with a normal distribution to determine the interactive effect of δ 15 N, δ 13 C and THg levels on PRL secretion. The dependent variable was PRL, with δ 15 N, δ 13 C and THg values included as fixed independent variables in all single, two-and three-way interactions. Additionally, we included female body mass (g) and clutch size as fixed independent variables in a two-way interaction. All independent variables were standardized and centred for comparison to remove misrepresentation of effects and increase the interpretability of coefficients [START_REF] Schielzeth | Simple means to improve the interpretability of regression coefficients[END_REF]. We included breeding colony in the model as a random variable to control for inter-colony PRL variation, with the objective being to examine individual variation. All statistical analyses were conducted in RStudio version 2022.07.2 (RStudio Team). All model assumptions were met by visually inspecting histograms to confirm a normal data distribution, quantile-quantile plots, and residual plots from the model outputs. We also confirmed the absence of multicollinearity between all variables using correlation coefficients and variance inflation factors (VIF), which were below 0.36 and 3.90 (3-way interaction), respectively. In interpretations of three-way interactions, the interaction among two factors (A * B) differs along with levels of a third factor (C). To assist with the visualization of these complex relationships, significant interactions were graphed using the R package Interactions v.1.1.3 [START_REF] Long | Interactions: Comprehensive, User-Friendly Toolkit for Probing Interactions[END_REF]. This package calculated Johnson-Neyman intervals and simple slopes analysis to determine the significant slopes within the three-way interaction.

Results

Median PRL levels were variable among eider colonies, ranging from 245.18 ng/mL (Tern Island) to 348.99 ng/mL (Kongsfjorden), with the median PRL for all eiders at 287.60 ng/mL (Table 2, Figure 2). At the individual eider level, the lowest (142.72 ng/mL) and highest (489.20 ng/mL) values were found at the Kongsfjorden colony. Mean colony PRL levels significantly varied between colonies, with the Kongsfjorden colony having significantly higher mean PRL levels than every colony except for Grindøya (Figure 2). Conversely, Tern Island had significantly lower mean PRL levels than all colonies but Faroe Islands (Figure 2). Although we report the PRL values for the John's Island colony, researchers at this location did not collect body mass measurements.

Therefore, this colony was omitted from the statistical model (n=18) and any individuals that did not have body mass and/or clutch size data (n=12, total individuals omitted = 30).

The GLMM showed that δ 15 N (t 80 =2.289, p=0.025; Table 3) as well as the three-way interaction between δ 15 N, δ 13 C and THg were significant (t 95 =2.089, p=0.039; Table 3). Johnson-Neyman intervals and simple slopes analysis determined that 6/9 of the slopes in the three-way interaction were significant. In particular, 3/3 slopes were significant at THg +1 standard deviation (SD), 2/3 were significant at mean THg levels, and 1/3 were significant THg -1 SD (Figure 3). At THg at +1 SD, δ 13 C values at -1 SD, mean levels and at +1 SD all resulted in significant δ 15 N slopes on PRL (Figure 3). At mean THg levels, the slopes of δ 15 N at mean δ 13 C and δ 13 C -1 SD were also significant (Figure 3). THg at -1 SD resulted in only the slope of δ 15 N at δ 13 C -1 SD having significance (Figure 3). However, there was no significant relationship between female body mass, clutch size, or a two-way interaction of those variables with PRL (Table 3).

Discussion

We examined the interactive effects of individual foraging ecology and THg contamination on PRL levels in incubating female common eiders across six colonies at a pan-Arctic scale. We predicted an effect of the interaction between δ 15 N, δ 13 C and THg on PRL levels, given the potential effect of these variables on female body condition and hormone regulation. Consistent with our prediction, individuals at lower trophic levels (lower δ 15 N) showed the lowest PRL levels, most commonly in combination with individuals at higher THg levels and those foraging in predominantly phytoplankton-based environments (lower δ 13 C). Our findings suggest a potential cumulative, interactive effect of foraging ecology and elevated THg exposure on PRL levels in wild eider colonies.

Interactive relationship between foraging ecology and PRL levels

We found a significant positive relationship between δ 15 N and PRL levels. Pre-breeding eiders foraging at lower trophic levels may potentially have a relatively lower body condition going into and then throughout incubation. PRL and incubating body condition have been linked in previous studies, where lower-condition females showed reduced PRL levels [START_REF] Criscuolo | Factors affecting plasma concentrations of prolactin in the common eider Somateria mollissima[END_REF].

Although body mass did not correlate with PRL, there is nonetheless potential for the indirect linkages between trophic level (i.e., δ 15 N) and body condition to influence PRL secretion (see section 4.3). Changes in food web structure are apparent with continuing climate change and are likely to vary spatio-temporally between colonies with differing prey availability (see section 1; [START_REF] Smith | Environmental and life-history factors influence inter-colony multidimensional niche metrics of a breeding Arctic marine bird[END_REF]). These effects are important for a species such as eiders since a lower body condition beginning in the pre-breeding stage could reduce the amount of PRL during incubation, with potential implications on nest attentiveness and overall breeding success [START_REF] Erikstad | Clutch-size determination in precocial birds: a study of the common eider[END_REF][START_REF] Criscuolo | Factors affecting plasma concentrations of prolactin in the common eider Somateria mollissima[END_REF]. Eiders forage primarily on amphipods and bivalves, which have demonstrated sensitivity to warming ocean temperatures [START_REF] Węsławski | Climate change effects on Arctic fjord and coastal macrobenthic diversity-observations and predictions[END_REF][START_REF] Dalpadado | Distribution and abundance of euphausiids and pelagic amphipods in Kongsfjorden, Isfjorden and Rijpfjorden (Svalbard) and changes in their relative importance as key prey in a warming marine ecosystem[END_REF][START_REF] Vinarski | Freshwater Mollusca of the Circumpolar Arctic: a review on their taxonomy, diversity and biogeography[END_REF]. Based on our current results, changes in eider trophic niches due to ecosystem shifts (likely to vary widely and depend on spatio-temporal prey availability) may have downstream impacts on body condition and PRL levels. For example, individuals that consume a greater abundance of lower trophic level prey could have comparatively lower body condition with expected negative impacts on PRL levels. This prey may include bivalves, which have lower energy returns due to a low ratio of flesh content to shell size compared to shell-less macroinvertebrates [START_REF] Hamilton | Size-selective predation of blue mussels (Mytilus edulis) by common eiders (Somateria mollissima) under controlled field conditions[END_REF][START_REF] Schummer | Dietary overlap of sympatric diving ducks during winter on northeastern Lake Ontario[END_REF]. In colonies where trophic niches shift to a greater abundance of higher trophic level prey, such as arthropods, the reverse effect is possible where higher trophic level prey may improve reproductive output, as seen in blue-footed boobies [START_REF] González-Medina | High quality diet improves lipid metabolic profile and breeding performance in the blue-footed booby, a long-lived seabird[END_REF]. Interestingly, the colonies with the lowest median PRL levels (Faroe Islands, Tern Island) were also identified previously as having the lowest median δ 15 N out of the colonies studied [START_REF] Smith | Environmental and life-history factors influence inter-colony multidimensional niche metrics of a breeding Arctic marine bird[END_REF], indicating that their lower PRL values may be linked with their utilization of lower trophic level prey such as bivalves. Despite these potential patterns, the eventual outcome of prey abundance shifts on eider hormonal mechanisms at the heart of successful incubation outcomes remains to be seen as climate change effects increase annually. However, our work provides a valuable baseline to examine subsequent future changes.

We also detected a significant relationship between δ 13 C and PRL levels in combination with δ 15 N and THg. Along with the food web changes described with δ 15 N above, the most consistently detected relationship was associated with phytoplankton-based foraging (low δ 13 C) and PRL levels in our three-way interaction. Due to the increasing effects of climate change, phytoplankton blooms are expected to become more intense and occur earlier in the year [START_REF] Kahru | Are phytoplankton blooms occurring earlier in the Arctic?[END_REF][START_REF] Ardyna | Phytoplankton dynamics in a changing Arctic Ocean[END_REF]. Phytoplankton blooms are also increasingly apparent when sea ice is still present, even in the benthos, due to increased light transmission with thinning ice [START_REF] Arrigo | Massive phytoplankton blooms under Arctic sea ice[END_REF][START_REF] Horvat | The frequency and extent of sub-ice phytoplankton blooms in the Arctic Ocean[END_REF][START_REF] Shiozaki | Bottomassociated phytoplankton bloom and its expansion in the Arctic Ocean[END_REF]. In eiders, phytoplankton blooms in the Baltic Sea have been linked to a reduction in nesting eiders [START_REF] Larsson | Effects of an extensive Prymnesium polylepis bloom on breeding eiders in the Baltic Sea[END_REF]. In this study, the relationship was suggested to be connected to an impact of the phytoplankton bloom on blue mussel (Mytilus edulis) availability (either through toxic or non-toxic pathways), followed by a subsequent reduction in eider body condition that impacted their breeding propensity [START_REF] Larsson | Effects of an extensive Prymnesium polylepis bloom on breeding eiders in the Baltic Sea[END_REF]. Our results showing a consistently significant effect of low δ 13 C in our three-way model on PRL may be related to phytoplankton-based environments resulting in less prey availability (such as blue mussels), either quality or quantity, that affects eider body condition and subsequently PRL.

Conversely, we also detected a significant relationship between sea ice algae-based foraging (high δ 13 C) and PRL levels in our three-way interaction, only in high THg individuals.

The relationship between high δ 13 C and PRL, along with elevated THg, may be due to constraints in foraging areas with higher sea ice cover, as seen previously in eiders that rely on access to open water to reach benthic and pelagic areas [START_REF] Jean-Gagnon | The impact of sea ice conditions on breeding decisions is modulated by body condition in an arctic partial capital breeder[END_REF]). As such, constraints in prey availability, combined with a high proportion of ice cover, could result in lower energetic stores (i.e., body condition) available to eiders during incubation, with possible downstream impacts on PRL levels.

Overall, the relationships we have detected between foraging ecology (based on isotopic niche) and PRL provide an initial look into how pre-breeding prey selection may indirectly influence key reproductive hormones in eiders. Future studies may benefit from integrating research into whether pre-breeding prey availability and quality may influence incubating PRL.

Importance of THg levels in relation to incubating PRL

The interactive relationship between δ 13 C and δ 15 N with PRL, as described in section 4.1, was most often significant in combination with high THg. The importance of THg in this relationship could potentially relate to the endocrine-disrupting properties of THg [START_REF] Zhu | The endocrine disruptive effects of mercury[END_REF]Tan et al. 2009). In addition, negative relationships between THg and body condition, as seen in previous studies [START_REF] Provencher | Implications of mercury and lead concentrations on breeding physiology and phenology in an Arctic bird[END_REF], but see [START_REF] Carravieri | Quantitative meta-analysis reveals no association between mercury contamination and body condition in birds[END_REF], may have also played a role in the effect of THg on PRL in combination with δ 13 C and δ 15 N. Importantly, our findings support concern over THg as an important component of cumulative stressors systems in the Arctic, where the combination of stressors such as contaminants, climate change and food web dynamics are expected to have cumulative impacts on birds [START_REF] Chastel | Mercury contamination and potential health risks to Arctic seabirds and shorebirds[END_REF][START_REF] Esparza | Mercury, legacy and emerging POPs, and endocrine-behavioural linkages: Implications of Arctic change in a diving seabird[END_REF][START_REF] Smith | Exposure to cumulative stressors affects the laying phenology and incubation behaviour of an Arctic-breeding marine bird[END_REF]. Overall, we found that THg played an important interactive role along with foraging ecology on PRL levels. However, whether this is related to the direct effects of THg on PRL, or indirectly mediated by foraging ecology or other intermediate physiological parameters is unclear and requires further study.

Other physiological and external factors with a potential influence on PRL

We were able to examine the relationship between body mass and PRL levels in this study. However, while previous studies successfully used this metric to determine eider body condition [START_REF] Descamps | Individual optimization of reproduction in a long-lived migratory bird: a test of the condition-dependent model of laying date and clutch size[END_REF], this measure of female condition may have been an incomplete look into the relationship between body condition and PRL in incubating eiders. Body condition, often measured as body mass divided by total head length [START_REF] Jamieson | An evaluation of methods used to estimate carcass composition of common eiders Somateria mollissima[END_REF] or scaled mass index [START_REF] Peig | New perspectives for estimating body condition from mass/length data: the scaled mass index as an alternative method[END_REF], may provide a more in-depth look into female quality at the time of sampling, rather than only body mass which does not adjust for overall body size. Additionally, individual body mass declines as eiders progressively use their endogenous resources during their incubation fast [START_REF] Criscuolo | Factors affecting plasma concentrations of prolactin in the common eider Somateria mollissima[END_REF][START_REF] Hanssen | Incubation start and egg size in relation to body reserves in the common eider[END_REF]. We were unable to control the incubation stage or body size in this study, but this would be ideal for future research on this topic.

It is worth highlighting that the primary glucocorticoid hormone in birds, corticosterone (CORT), is an important factor in eider pre-incubation body condition, where higher CORT corresponds to more intensive foraging [START_REF] Hennin | Pre-breeding energetic management in a mixed-strategy breeder[END_REF][START_REF] Evers | The effects of methylmercury on wildlife: a comprehensive review and approach for interpretation[END_REF]. Baseline CORT may play an essential regulating role in the relationship we have detected between PRL and foraging ecology during pre-breeding periods. For example, during pre-breeding foraging, females with higher baseline CORT may result in downstream elevated PRL levels during incubation due to more intensive foraging, leading to higher body conditions at the start of incubation [START_REF] Criscuolo | Factors affecting plasma concentrations of prolactin in the common eider Somateria mollissima[END_REF][START_REF] Angelier | Does prolactin mediate parental and life-history decisions in response to environmental conditions in birds? A review[END_REF]. Additionally, elevated baseline CORT has also been linked to a greater need for resources during incubation, as seen in incubating macaroni penguins (Eudyptes chrysolophus; [START_REF] Crossin | Corticosterone predicts foraging behavior and parental care in macaroni penguins[END_REF]. Eiders, which fast during one single incubation bout, have demonstrated a negative relationship between CORT and PRL during incubation periods in previous studies that regulates incubation behaviour [START_REF] Criscuolo | Corticosterone alone does not trigger a short term behavioural shift in incubating female common eiders Somateria mollissima, but does modify long term reproductive success[END_REF][START_REF] Parmesan | Ecological and evolutionary responses to recent climate change[END_REF]. The same pattern has been seen in birds that forage during incubation, with stressed individuals leaving more often on foraging trips (Angelier and Chastel 2009;Angelier et al. 2009). In addition, elevated CORT levels may negatively affect female immunity and mass loss during incubation, as seen in previous common eider studies [START_REF] Bourgeon | selectively decreases humoral immunity in female eiders during incubation[END_REF][START_REF] Bourgeon | Proximate and ultimate mechanisms underlying immunosuppression during the incubation fast in female eiders: roles of triiodothyronine and corticosterone[END_REF]). Higher energetic stress could have foraging ecology-mediated links to PRL levels. While high CORT levels may be initially beneficial during the pre-breeding period, these elevated may be detrimental to PRL during incubation, especially in fasting female eiders. Further implications may also arise from the added effects of elevated stress-induced CORT, especially as predator-prey relationships and climate conditions shift in Arctic systems.

Many eider colonies are facing increasing predation pressure from polar bears (Ursus maritimus), arctic foxes (Vulpes lagopus) and Larus sp. gulls [START_REF] Smith | Shorebird incubation behaviour and its influence on the risk of nest predation[END_REF][START_REF] Hanssen | A natural antipredation experiment: Predator control and reduced sea ice increases colony size in a long-lived duck[END_REF][START_REF] Prop | Climate change and the increasing impact of polar bears on bird populations[END_REF][START_REF] Jagielski | Polar bears are inefficient predators of seabird eggs[END_REF]. Eiders must also contend with elevated thermal environments that can influence female incubation behaviour, especially when combined with Hg exposure [START_REF] Smith | Exposure to cumulative stressors affects the laying phenology and incubation behaviour of an Arctic-breeding marine bird[END_REF]. Interestingly, studies on eiders in the Baltic Sea have determined that those with elevated PRL and greater body condition may have higher nest predation risk [START_REF] Mohring | Parental Investment Under Predation Threat in Incubating Common Eiders (Somateria mollissima): A Hormonal Perspective[END_REF]. However, the correlation between PRL and elevated predation risk was speculated to be attributed to an absence of low-quality individuals with reduced PRL levels at nesting sites with higher predation risk [START_REF] Mohring | Parental Investment Under Predation Threat in Incubating Common Eiders (Somateria mollissima): A Hormonal Perspective[END_REF]. In general, lower PRL may result in a greater risk for egg predation due to lowered incubation consistency. In comparison, higher PRL, thus greater incubation consistency, may lead to higher female predation risk. However, further research is needed to establish these links. Ultimately, the cumulative stressors eiders are now more often exposed to both pre-breeding and during incubation have the potential to greatly decrease nesting success with a 'perfect storm' of stressors influencing reproductive parameters via PRL secretion.

Conclusion

The interactive effects of foraging ecology and THg exposure on PRL levels were examined in incubating female common eiders. Our results show that individuals with high THg that were interactively foraging at lower trophic levels (low δ 15 N) and in areas with greater phytoplankton abundance (low δ 13 C) had the most significant relationships resulting in lowered PRL levels.

Overall, this complex interactive system, based on foraging ecology and Hg exposure, may have downstream implications for reproductive success through hormone alteration and, therefore, population demography. These effects may be strengthened within the context of continually changing environmental conditions on a global scale. However, further studies are required to fully understand these complex relationships and their reproductive and demographic consequences in a changing Arctic.

Acknowledgments.

  Thank you to the Isotopic Spectrometry Facility at LIENSs (La Rochelle University -CNRS, La Rochelle) and the Service d'Analyses Biologiques at CEBC for their help with sample analysis. We also thank Environment and Climate Change Canada (ECCC), the Natural Sciences and Engineering Research Council of Canada, Canada Research Chairs Program, the Nunavut Wildlife Management Board, Northern Scientific Training Program, Polar Continental Shelf Project, POLAR Canada, the Canadian Network of Centres of Excellence (ArcticNet), the Government of Ontario, the University of Windsor and the Arctic Institute of North America for logistical support, research and/or personal funding. The IUF (Institut Universitaire de France) is acknowledged for its support to P. Bustamante as a Senior Member.This work also contributes to the ARCTOX initiative (https://arctox.cnrs.fr) and the ARCTIC-STRESSORS project (ANR-20-CE34-0006-01).

Figure 1 .

 1 Figure 1. Locations of the seven pan-Arctic and subarctic common eider colonies used in this

Figure 2 .

 2 Figure 2. Boxplot depicting whether the mean PRL levels at each common eider colony vary

Figure 3 .

 3 Figure 3. Interactive effect of δ 15 N, δ 13 C and THg (visualized using standard deviation from the

Table 1 .

 1 Results of model selection on plasma PRL levels of 111 female common eiders from 6 breeding colonies. The number of parameters (k), corrected Akaike information criterion (AICc), change in AICc (ΔAICc), AICc weight, and log-likelihood are provided for each competing model.Colony was included as a random variable in all models, however, when clutch size was included as a random variable it is denoted with R .

	Model	k	AICc	ΔAICc	AICc Weight	Log-likelihood
	THg*δ 15 N*δ 13 C+Mass*Clutch 13	1106.61	0.00	0.80	-538.33
	THg*δ 15 N*δ 13 C+Mass+Clutch	12	1109.70	3.09	0.17	-541.17
	THg*δ 15 N+THg*δ 13 C+ δ 15 N*δ 13 C +Mass*Clutch	12	1113.91	7.30	0.02	-543.28
	THg*δ 15 N*δ 13 C+Mass+Clutch R 12	1115.73	9.12	0.01	-544.19
	THg*δ 13 C+δ 15 N+Mass+Clutch	9	1123.24	16.63	0.00	-551.68
	δ 13 C*δ 15 N+Mass+Clutch	8	1126.41	19.80	0.00	-554.46
	THg+δ 15 N+δ 13 C+Mass+Clutch	8	1126.57	19.96	0.00	-554.54
	δ 15 N+Mass+Clutch	6	1136.47	29.86	0.00	-561.81
	δ 13 C+Mass+Clutch	6	1139.02	32.40	0.00	-563.08
	THg+Mass+Clutch	6	1139.53	32.91	0.00	-563.34
	THg*δ 15 N*δ 13 C+Mass	11	1173.26	66.64	0.00	-574.30
	THg*δ 15 N*δ 13 C+Clutch	11	1191.18	84.57	0.00	-583.28
	THg*δ 15 N*δ 13 C	10	1457.05	350.44	0.00	-717.64

Table 3 .

 3 Results of a general linear mixed effects model examining the interactive effects of THg (µg/g dry weight), δ 15 N and δ 13 C (‰) as well as an interaction between body mass (g) and clutch size on prolactin (PRL) levels (ng/mL). The breeding colony (n=6) of the 106 female arcticbreeding common eiders throughout Canada and Eurasia was a random variable in the model, with the John's Island eiders and 12 additional individuals from other colonies removed from the model due to the absence of morphometrics (significant p-values are bolded).

		Estimate	Standard Error	df	t-value	p-value
	Intercept	306.746	19.540	4.526	15.599	<0.001
	Body Mass	-9.472	7.375	93.941	-1.284	0.202
	Clutch Size	4.628	6.389	94.890	0.724	0.471
	Mass*Clutch	-0.36	6.872	92.220	-0.053	0.958
	THg	-7.108	9.509	93.910	-0.747	0.457
	δ 13 C	-17.543	11.313	77.172	-1.551	0.125
	δ 15 N	25.427	11.107	80.279	2.289	0.025
	THg*δ 13 C	7.484	8.609	94.777	0.869	0.387
	THg*δ 15 N	6.437	5.918	94.889	1.088	0.280
	δ 15 N*δ 13 C	-13.376	10.882	92.772	-1.229	0.222
	THg*δ 15 N*δ 13 C	14.138	6.767	94.561	2.089	0.039