The genetics of immune and infection phenotypes in wild mice, Mus musculus domesticus

Wild animals are under constant threat from a wide range of micro‐ and macroparasites in their environment. Animals make immune responses against parasites, and these are important in affecting the dynamics of parasite populations. Individual animals vary in their anti‐parasite immune responses. Genetic polymorphism of immune‐related loci contributes to inter‐individual differences in immune responses, but most of what we know in this regard comes from studies of humans or laboratory animals; there are very few such studies of wild animals naturally infected with parasites. Here we have investigated the effect of single nucleotide polymorphisms (SNPs) in immune‐related loci (the major histocompatibility complex [MHC], and loci coding for cytokines and Toll‐like receptors) on a wide range of immune and infection phenotypes in UK wild house mice, Mus musculus domesticus. We found strong associations between SNPs in various MHC and cytokine‐coding loci on both immune measures (antibody concentration and cytokine production) and on infection phenotypes (infection with mites, worms and viruses). Our study provides a comprehensive view of how polymorphism of immune‐related loci affects immune and infection phenotypes in naturally infected wild rodent populations.


| 4243
CHEYNEL et al. Penn et al., 2002). Moreover, there is a much wider repertoire of immune-related loci that may also affect individuals' immune responses (Acevedo-Whitehouse & Cunningham, 2006). Animals use Toll-like receptors (TLRs) to detect bacterial and viral pathogenassociated molecular patterns (PAMPs). These TLRs then initiate signalling events, including via cytokines and chemokines, that ultimately activate innate and adaptive immune responses contributing to parasite resistance (Kawai & Akira, 2006). Most mammals have 10-12 TLR coding loci that detect different PAMPs (Roach et al., 2005).
In humans, numerous studies have described associations between polymorphisms in TLR and cytokine-coding loci, and susceptibility or resistance to a range of infections (Mukherjee et al., 2019;Smith & Humphries, 2009). In humans, the cumulative contribution of non-MHC loci, such as those coding for cytokines or PAMP receptors, to the immune phenotype likely exceeds that of the MHC (Jepson et al., 1997;Roederer et al., 2015).
Most of our understanding of the effect of polymorphism in immune-related loci on immune and/or infection phenotypes comes from studies of humans or laboratory animals; few empirical studies have been conducted with wild animals. Major differences exist between the infection and immune states of wild animals and their laboratory-reared counterparts. This has been extensively studied for mice, showing the significantly different immunological effect of wild vs. laboratory gut microbiomes (Abolins et al., 2017;Beura et al., 2016;Rosshart et al., 2017). Outbred lines of laboratory mice are available, though these do not reflect the full genetic diversity of wild populations (Turner & Paterson, 2013). Given the disconnect between the immunological state, infection state and genetic diversity of wild and laboratory animals, the paucity of our understanding of how host genetics affects immune responses and infection phenotypes in wild animals is a notable knowledge gap, and points to the critical need to study these phenomena directly (Turner & Paterson, 2013).
To date, there have been about 30 immunogenetic studies on 16 different genera of wild mammals that sought specific associations between polymorphisms in immune-related loci and measures of immune or infection phenotypes (Table 1). Among these the effect of the MHC has been investigated most frequently and these studies commonly, but not universally, find effects of MHC diversity on micro-and macroparasite infection phenotypes. Studies of polymorphisms in TLR and cytokine-coding loci have also found associations with infection phenotypes. Most of these studies have studied the effects of polymorphism in immune-related loci on infection phenotypes (especially of macroparasite burden), and so there remains major gaps in our understanding of effects on immune responses themselves (but see Charbonnel et al., 2010;Cutrera et al., 2011;Turner et al., 2011;Huang et al., 2022). Overall, polymorphism in rather few immune-related loci has been studied in wild mammals, again in notable contrast to the very large number of loci studied in laboratory animals. It is not yet known whether these laboratory observed effects also occur in the wild. An important remaining question is therefore whether-and if so, how-naturally occurring polymorphism in immune-related loci impacts both immune and infection phenotypes in wild mammals. Wild mice have the potential to be powerful study systems to improve our understanding of these relationships, because they will reveal effects occurring in natural environments, but with the advantage of exceptionally good genetic, physiological and immunological understanding of the study species, which has accumulated from decades of experimental laboratory work (Graham, 2021;Turner & Paterson, 2013).
Here we have investigated the effect of diversity in immunerelated loci on a wide range of immune and infection phenotypes in wild house mice, Mus musculus domesticus. We first assessed the extent of genetic variation in 23 immune-related loci (coding for cytokines, TLRs or the MHC) in our study populations, finding variation in 14 of these loci among 435 mice. We then tested the effect of each genotype found in these 14 loci on 24 discrete immune phenotypes, encompassing both cellular and humoral components of the innate and the adaptive immune response. We also tested the effect of the genotype on infection phenotypes of nine different microand macroparasites (mites, worms, viruses and a bacterium). To look for effects of the genotype on the phenotype, we used a two-step model selection approach, firstly testing the effect of all non-genetic factors (site, sex, age and body mass), then adding each genetic factor to the model in turn to assess its effect on the relevant phenotype. Our study thus provides a uniquely complete picture of the association between genotype and immune phenotype, linking allele in 14 immune-related loci to 24 immune phenotypes and 9 infection phenotypes in wild rodents.

| Mice
The sampling, processing, determination of infection state and immune phenotyping of the mice have previously been reported (Abolins et al., 2017(Abolins et al., , 2018. In brief, 460 mice (Mus musculus domesticus) were live trapped from 12 sites in the southern United Kingdom (UK) between March 2012 and April 2014. Population genetic analyses show strong genetic structuring of mice among the different sample sites (Abolins et al., 2018). After being humanely killed, mice were sexed, weighed and measured from the tip of the snout to the base of the tail. Female mice that were later found to be pregnant had the mass of the foetuses subtracted from their total mass, and these values were used in all subsequent analyses. To estimate body condition, the Scaled Mass Index (SMI) was calculated as previously described (Peig & Green, 2009). Age was determined using eye lens mass as previously described (Rowe et al., 1985). The distribution of individuals by study site, sex and age class is given in Table S1.

| Immune measures
The immune state of the mice was assessed by measuring (i) 12 immune cell populations, (ii) the concentration of 3 immunoglobulins, and (iii) the concentration of 9 cytokines produced after in vitro TA B L E 1 Studies in wild mammals, testing for associations between genetic variation in specific immune-related loci (i.e. MHC, TLR and/or cytokines) and measures of immune or infection phenotypes. These articles were searched for from the Web of Science using combinations of key words 'immuno-genetic', 'genetic variation ', 'immunity', 'infection' and 'wild' on

TA B L E 1 (Continued)
stimulation of spleen cells with 4 different stimuli (anti-CD3/anti-CD28, CpG, lipopolysaccharide and peptidoglycan). Together, these measures therefore encompass cellular and humoral components of the innate and adaptive immune responses, so providing a broad assessment of humoral and cellular immune responses (as reported in Abolins et al. (2017Abolins et al. ( , 2018).
The 12

Current or prior infection with Corona, Mouse Hepatitis, Sendai,
Minute, Noro and Parvo viruses and to Mycoplasma pulmonis was inferred from serological data, as previously described (Abolins et al., 2017).
The control locus Myo1a was also successfully amplified, sequenced and two SNPs (Myo1a_1, Myo1a_2) were found to be polymorphic in these mice. A summary table of the SNPs identified and variants found is provided in Table S2.
The PCR amplification of these fragments used the primer sequences shown in Table S3.  Table S4 and variants detected are shown in Table S5). Successful PCR products were purified using the GeneJet PCR clean-up kit (Thermo Fisher Scientific), following the manufacturer's instructions and single-read sequencing of the purified products was carried out by Eurofins Genomics, UK.

| Statistical analyses
We tested all possible associations between individual SNPs and immune or infection phenotypes. We used generalized linear models for continuous variables (immune measures and worm number) and ordinal logistic regression for categorial data (i.e. intensity of infection with mites and microparasites). Each immune or infection phenotype was used as the response variable in the model. Raw data were used when the immune parameter had a Gaussian distribution, while traits with a non-Gaussian distribution were log(X + 1) transformed before analysis to improve the normality of the data.
As our dataset comprised concentrations of nine different cytokines produced in vitro after spleen cell stimulation with four different stimuli (i.e. 36 cytokine measures for each mouse), we performed a principal component analysis (PCA) on all mice for which all 36 cytokine measures were available (N = 172) to identify the main axes of variation, for which we used the R package 'ade4' (Dray & Dufour, 2007). We removed data for three mice that contributed abnormally to the PCs (i.e. contributing 8% of the first principal component) possibly because of a severe inflammatory state. PC1 and PC2 represent 34% and 19% (total 53%), respectively, of the covariation among the cytokine measures. PC1 was mostly influenced by IL-1B, IL-10, IL-12p40, IL-12p70, IL-13, IL-4; PC2 was mostly influenced by IL-6, IFNγ and MIP-2α (Table S6). PC1 thus appears to represent overall cytokine concentration, weighted more to specific cytokines and PC2 to Th1/Th2-related cytokines. For cytokine measures, we therefore used PC1 and PC2 as the response variables in the modelling approach below.
Our modelling approach had two steps. Firstly, for each immune or infection phenotype, we built a base model to investigate the effect of potential non-genetic confounders, specifically site effects (nine sites, Figure 1), age (linear function and three classes, Table S1), sex (two classes: male and female), and body condition (SMI, continuous).
We merged data from some closely located sites with small numbers of mice, giving us nine sites for analysis across the whole population ( Figure 1; Table S1). The three age classes were chosen to differentiate individuals ecologically according to their behavioural and lifehistory state, linked to their likelihood of encountering a range of environmental antigens as: 0-6 weeks (immature animals unlikely to be venturing far from their nest site); 7-12 weeks (young animals unlikely to have reproduced); >12 weeks (sexually mature adults). We tested the two-way interactions between site, age and sex. Among female mice, 15% were pregnant. We tested the effect of pregnancy on immune traits, finding that there was only an effect on the cytokine measures, and so we included pregnancy as a random factor in models analysing the cytokine PCAs. Secondly, we tested the effect of the SNPs in the immune-related loci, with the different allele combinations tested as fixed effects (three classes: homozygous 1, homozygous 2, heterozygous). Depending on the variables retained in the base model, two-way interactions between SNP, site, sex and age were also included. Using a model selection procedure based on the Akaike information criterion (AIC, Burnham & Anderson, 2002), we retained the model with the lowest AIC. When the difference in AIC between competing models was less than 2, we retained the model with the fewest parameters to satisfy parsimony rules (Burnham & Anderson, 2002). Residuals plots were checked to ensure the fit of the selected regression models. Finally, to be more selective about the associations identified in our study, the significance of the genetic effects selected with the model testing was also tested using deletion testing and the log-likelihood ratio test (McCullagh & Nelder, 1989).
The p-values obtained from this additional deletion testing step are reported in Table 2 and provided in Table S12. All statistical modelling was performed using R version 4.0.5 (R Core Team, 2021).

| RE SULTS
We genotyped 435 mice for 23 SNPs in 21 immune-related loci; 14 of these 23 SNPs were found to be polymorphic. We successfully genotyped 378 (87%) mice for all 14 SNPs; the remaining 57 (13%) mice were successfully genotyped for at least 11 SNPs. The minor allele frequencies and the observed and expected heterozygosity are shown in Tables S7 and S8, respectively. The distribution of each SNP among sample sites is shown in Table S9.

| Genetic variation in immune-related loci
There was a high level of variation in the MHC loci within our wild mouse population. H2-Ab was the most diverse H-2 locus, with three different genotypes present at similar frequencies and with a heterozygosity of 0.33. For H2-Aa and H2-Eb, one genotype dominated (found in 85% and 70% of mice, respectively) and with lower heterozygosity (0.09 and 0.21, respectively).
For four of the cytokine-coding loci (IL-1a, IL-1b_U, IL-17a_U and IL-17F), there was relatively high level of genetic diversity, with overall heterozygosity between 0.12 and 0.34. In contrast, the other cytokines were each dominated by one very common genotype (89%-95%), and had low heterozygosity (0.01-0.06). A similar situation pertains for the two TLR loci where one genotype dominated (95% and 96% for TLR5 and TLR9, respectively) with a heterozygosity of 0.04.

TA B L E 2
Effect of polymorphism in immune-related loci on immune phenotypes. Results of linear models where first the effect of different age functions (linear, classes), sex (female vs. male), site (9 sites), and scaled mass index were tested to select the base model of non-genetic factors. Secondarily, the effect of each polymorphism was tested and shaded boxes indicate where genetic effects were selected in the model (Table S10).

| Genetic effects on immune phenotypes
We    Table 2). We then tested the effect of the three different combination of alleles of each SNP on the immune phenotype (Tables 2 and 3). We found no relationships between the control locus and the immune and infection phenotypes. Full details of the two-step model selection are reported in Table S10, effect estimates of all the models selected are shown in Table S11, and the p-values of the log-likelihood ratio test after deletion testing are reported in Table S12.

| Antibody concentration
We found strong associations between the genotype of the MHC and cytokine-coding loci and the concentration of antibodies (  (Figure 2b).
Finally, the MHC H2-Ab genotype showed an association with the concentration of IgG (Table 2), where heterozygous mice had lower concentrations than either of the two homozygous genotypes, that is GG and AA genotype (Figure 2b). While these patterns are evident in the whole data set, effects of IL-1a or IL-17a_U genotype on IgE concentration and of H2-Aa polymorphism on IgG concentration also varied among sites and among mouse age classes (Table 3,   Table S11).

| Cytokine production
We also found associations between the genotype of cytokinecoding loci and the PCA of cytokine responses (Table 2, Figure 2c).
Polymorphism in IL1-b_U was associated with PC1 (p < 0.05), which mainly reflects cytokines IL-1B, IL-10, IL-12p40, IL12p70, IL-3 and IL-4. Specifically, heterozygous GT or homozygous GG mice had higher PC1 values than the most common genotype TT (Table 3). The genotype of the MHC locus H2-Eb was associated with PC1 (p < 0.001), with heterozygous AG mice having lower PC1 values than the rare AA genotype or the most common GG genotype (Table 3). The genotype of IL-17a_U was associated with cytokine PC2 (p < 0.05), which mainly reflects cytokines IL-6, IFNγ and MIP-2α; heterozygous mice had higher PC2 values than the homozygous genotypes (Table 3).
Concerning the PCA analyses, we generally note that PC1 has negative loadings for all measures and thus indicates variation in the levels of cytokine production, while PC2 presents a more qualitative separation between IFNγ and IL-6 against IL-13 and IL-4, consistent with representing the Th1/Th2 phenotype.

| Immune cells
We found some associations between the MHC and cytokine loci genotypes and immune cell numbers (Table 2, all p > 0.05; Figure 2a).
The H2-Aa genotype was associated with the number of DCs, CD4 + and CD8 + T cells, and MDSC: mice with the rare AA genotype had more cells than either of the other genotypes (Figure 2a). In contrast, the rare AA H2-Eb genotype was associated with fewer CD4 + T cells compared to the other genotypes (Table 3). Finally, there were some associations between the genotype of IL-17F and IL-17a_U and the number of PMNs (Table 2). Mice heterozygous for IL-17F had more TA B L E 4 Effect of polymorphism in immune-related loci on infection phenotypes. Results of linear models where first the effect of different age functions (linear, classes), sex (female vs. male), site (9 sites) and scaled mass index were tested to select the base model of non-genetic factors. Secondarily, the effect of each polymorphism was tested and shaded boxes indicate where genetic effects were selected in the model (Table S10).

Sendai virus
Site + age (classes)

Mycoplasma pulmonis Site + age (classes)
Note: The stars shows the significance of the selected genetic effect, determined by deletion testing and the log-likelihood ratio test, as: *** for p < 0.001; ** for p < 0.01; * for p < 0.05 and for p > 0.05 (Tables S12). Loci where the minor allele frequency (MAF) is >0.1 are shown in bold (Table S7); loci with a MFA < 0.1 are not bold and these results should be interpreted with more caution.
PMNs than either homozygous genotype (Table 3), whereas mice with the most common GG IL-17a_U genotype had more PMNs than the other rarer genotypes (Table 3). As above, some of these associations varied among sample sites (Table S11).
For some loci, there was too little genetic variation to allow us to draw strong conclusions about their association with immune phenotypes. However, we found a recurring tendency for associations between alleles of IL-6, IL-10, IL-13 and TNF coding loci and many different components of the immune phenotype (Table 2). This interesting trend would benefit from further study. We did not, however, find any striking associations between the two TLR loci, nor the IL-17a_N cytokine locus, and the immune phenotype.

| Genetic effects on infection phenotypes
As above, we first tested the effect of non-genetic factors on infection phenotypes (Table 4). Again, we found strong site and age effects on virus seropositivity, mite and worm burden, and of SMI on mite and worm burden, as in Abolins et al. (2018) (Table 4). We then tested the effect of the three different allele combinations in each SNP on infection phenotypes (Tables 4 and 5). Full details of the two-step model selection are reported in Table S10, effect estimates of all the models selected are shown in Table S11, and the p-values of the log-likelihood ratio test after deletion testing are reported in Table S12.
We found a strong association between MHC H2-Aa genotype and mite burden: mice with the rare homozygous AA genotype had fewer mites than those with the other, more common genotypes (Tables 4 and 5, p < 0.05). MHC H2-Ab genotype was associated with worm burden (Table 4, p < 0.01), with a general tendency for mice with the AA genotype to have higher worm burdens than the two other genotypes. However, this general pattern obscures some local variation between sites (see all effect estimates of the H2-Ab x site interaction in Table S11); for instance, H2-Ab heterozygous mice from site WF + WT have higher worm burdens than the other genotypes.
We also found associations between the genotype of MHC loci and seropositivity for virus infections (Table 4). For example, mice with AA or AG H2-Ab genotypes were less likely to be seropositive for Minute virus compared with mice with the most common GG genotype ( were less likely to be seropositive for norovirus-but more likely to be seropositive for parvovirus-compared with the two other homozygous genotypes (Table 5, p < 0.05 for both).
We found only limited evidence of associations between genotypes of cytokine-coding loci and infection phenotypes (Tables 4 and 5), though mice that were heterozygous at the IL-1a locus were more likely to be seropositive for Minute virus genotypes (p < 0.05).
Finally, we found a recurring association of IL-6 alleles and both macroparasite (worm number, p < 0.01) and microparasite infection markers (MHV, sendai virus and M. pulmonis; Table 4, p < 0.05 for all). We also found an association of the TLR-9 locus and the mite infection phenotype (Table 4, p < 0.01). However, for these loci, the relative abundance of the different genotypes means that these results should be considered with caution and would benefit from further study.

| DISCUSS ION
We tested the effect of the genotype of 14 immune-related loci on a wide range of immune and infection phenotypes in wild UK house TA B L E 5 The effect of different genotypes on infection phenotypes. The parameter estimates for genetic effects for loci where the minor allele frequency > 0.1 are shown. Parameter estimates for other factors (age, site, SMI, etc.), interactions, and other loci are in Table S11. In Variable the genotype in brackets is the genotype to which the comparison is made. Statistical significance is represented by: for p < 0.1, * for p < 0.05, ** for p < 0.01 and *** for p < 0.001. Degrees of freedom is df. Most of the immune-related loci that we investigated were polymorphic, though there were differences in the representation of genotypes at loci among mice from different sites. Overall the observed genetic variation was most notable for the MHC and cytokine-coding loci, but was much less extensive for the TLR coding loci (known to be highly conserved in vertebrates), which may explain why we found comparatively fewer clear effects of polymorphism in TLR coding loci. We did, however, find strong associations between the MHC and cytokine-coding loci on different components of the immune phenotype-principally on antibody concentration and on cytokine production-but also on infection phenotypes. The effects of the MHC loci on immune and infection state are consistent with previous results in natural populations (Table 1). Interestingly, we found that the genotype of the MHC H2-Ab locus was associated with both IgG concentration and intestinal worm burden and there is the possibility that this is causal. Associations between MHC loci and antibody response in wild populations have been described in other mammals such as Soay sheep (Huang et al., 2022; but also see in other species reported in Gaigher et al., 2019). We also found associations between the MHC H2-Aa locus and various populations of DCs, CD4+ and CD8+ T cells, and MDSCs, as well as on ectoparasite burden.
We found major associations between the genotype of cytokinecoding loci and antibody concentration. Specifically, the serum concentration of IgE was associated with genotypes of IL-1a, IL1b_U, IL-6 IL-17a_U and TNF, with different IL-1a, IL1b_U, IL-17a_U genotypes having quite strikingly different concentrations. Genotypes of IL-6 were also associated with the concentration of faecal IgA. The potential mechanism of these effects on antibody concentration is that these polymorphisms affect (i) animals' antibody production per se, (ii) animals' resistance or susceptibility to a wide range of infections, which is manifest as different antibody concentrations or (iii) a combination of these.
Among the cytokine-coding loci that we investigated, we found ing immune responses to a wide range of micro-and macroparasites (Couper et al., 2008). Similarly, polymorphism in IL-6 also appears to be pleiotropic, which is also consist with IL-6's role in integrating different aspects of an immune response, and where in humans suppression of IL-6 function increases the risk of serious infection (Rose-John et al., 2017). However, we need to be cautious in interpreting these data since the measures of microparasite infection are serological, and so may be confounded with individuals' antibody production per se, which we also found to be affected by polymorphism in these and other loci.
While we focus on SNPs within these loci of immunological interest the associations we report are to the SNP and to regions that are linked to it, and so it would be inappropriate to seek causal relationships between the focal SNP and altered protein structure as underlying any direct or indirect immunological effect. While it would be of interest to understand how such genetic variation alters protein structure and function, our study was not designed to study such mechanisms.
Our measures of viral infection were serological, and while this is an easy-to-use approach it is potentially problematic because it is itself an immunological measure. An advantage of serological diagnosis is that it shows evidence of historical infection, though in the mice studied here their median age is 6-7 weeks (maximum 20-39, for male and females, respectively; Abolins et al., 2017), meaning that in many cases these serological diagnoses are often likely of relatively recent viral infection. Future work would benefit by using more specific viral diagnosis, as well as considering a wider range of parasites, but also host pathology. Although it is difficult to compare the relative magnitude of the genetic effects between the diverse immune and infection traits measured with different methods and scales, the strong associations that we have observed between genetic effects and antibody concentrations or cytokine production appear to be especially significant.
For some of the immune-related loci (including IL-6, IL-10, IL-13, TNF, TLR), there was a substantial imbalance in the frequency of different genotypes that limited our analytical power. Previous studies have noted that effective sample sizes can limit what can be studied in wild populations and have suggested minimum sample sizes of at least 200 individuals (Gaigher et al., 2019). It is notable that even with our large sample size of more than 400 mice, the relative under-representation of some genotypes has constrained the power of our analyses. Previous genetic analyses of these same mice using putatively neutral loci showed genetic differentiation among mice from the different sites (Abolins et al., 2018), again consistent with differences in genotype frequencies in the immune-related loci that we have observed. Moreover, previous analyses of the immune distance among these mice found that mice within one sample site were on average more immunologically similar to each other, compared with mice from other sample sites (Abolins et al., 2018). In agreement with both of these analyses, in the present study we also found significant effects of sample site on immune and infection phenotypes. Genetic and immunological differences among mice at different sites are therefore a notable feature of the mouse populations that we have studied. This host genetic and immunological diversity therefore presents a heterogeneous environment in which parasites are selected and are evolving, and where hosts continue to evolve to respond to that threat.

| CON CLUS ION
In conclusion, this work has conducted a large-scale study of the genetics of immune and infection phenotypes in wild mice, which adds to the rather limited set of studies in wild mammals. We have detected a range of genetic associations between multiple components of the immune response and on infection phenotypes. Of particular note, we have found significant associations of polymorphism with cytokine-coding loci on the serum concentration of IgE antibodies, and of MHC with the serum concentration of IgG antibodies. There is considerable genetic diversity in the loci that we have studied, suggestive of balancing selection acting, likely due to shifting patterns of infection and co-infection that occur in wild animal populations.
Results such as those we present here could be investigated further, for example by undertaking detailed structure-function studies in laboratory animal systems, to seek to understand how the genetic effects we describe functionally affect immune responses. However, such an approach needs to be cognisant of the context in which these immune responses are produced. Thus, while parasites do cause harm to their hosts, which immune responses can ameliorate, the costliness and risk of immunopathology means that maximal antiparasite immune response may not always maximize fitness.

AUTH O R CO NTR I B UTI O N S
Conceived and designed the experiments: EMR, MV. Performed the experiments: LL, EMR and MV. Analysed and interpreted the data: LC, LL, EMR and MV. Wrote the paper: LC, EMR and MV.

ACK N O WLE D G E M ENTS
This work was supported by NERC (NE/I022892/1). We would like to thank the farm owners and London Underground for generous access to their properties. We also would like to thank Jane Hurst and Steve Paterson for helpful discussions.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declared that they have no competing interests exist.

DATA AVA I L A B I L I T Y S TAT E M E N T
The entire dataset and the script used for the analyses are deposited