Genetic evidence of female kin clusters in a continuous population of a solitary carnivore, the Eurasian lynx

Abstract Large terrestrial carnivores can sometimes display strong family bonds affecting the spatial distribution of related individuals. We studied the spatial genetic relatedness and family structure of female Eurasian lynx, continuously distributed in southern Finland. We hypothesized that closely related females form matrilineal assemblages, clustering together with relatives living in the neighboring areas. We evaluated this hypothesis using tissue samples of 133 legally harvested female lynx (from year 2007 to 2015), genotyped with 23 microsatellite markers, and tested for possible spatial genetic family structure using a combination of Bayesian clustering, spatial autocorrelation, and forensic genetic parentage analysis. The study population had three potential family genetic clusters, with a high degree of admixture and geographic overlap, and showed a weak but significant negative relationship between pairwise genetic and geographic distance. Moreover, parentage analysis indicated that 64% of the females had one or more close relatives (sister, mother, or daughter) within the study population. Individuals identified as close kin consistently assigned to the same putative family genetic cluster. They also were sampled closer geographically than females on average, although variation was large. Our results support the possibility that Eurasian lynx forms matrilineal assemblages, and comparisons with males are now required to further assess this hypothesis.

possibilities for assessment of social organization as the cause of local family genetic structuring, independent of the potentially confounding effects of population fragmentation and geographic isolation. In recent decades, several large carnivore species have recolonized parts of their former distribution range and in some areas regained a continuous distribution over large unfragmented landscapes (Chapron et al., 2014).
Due to their different social organization, group-living and solitary large carnivores may differ in the frequency and strength in which they form long-lasting family bonds and, more generally, in the rate in which they interact with other individuals during their lifetime. Little is known about the association between social and family genetic structure, especially in solitary large carnivores. Behavioral studies of leopard (Panthera pardalis), puma (Puma concolor), and tiger (Panthera tigris) have revealed that also solitary species can exhibit kin clusters (Elbroch, Quigley, & Caragiulo, 2015;Fattebert et al., 2016;Goodrich et al., 2010;Logan & Sweanor, 2002), but the family genetic structure is less studied.
The Eurasian lynx (Lynx lynx) is a solitary predator with a social organization based on territoriality. The species has one of the most widespread distributions of the currently living felids (Breitenmoser et al., 2015). The lynx populations inhabiting Europe differ in their population history and degree of habitat discontinuity (von Arx, Breitenmoser-Wuersten, Zimmermann, & Breitenmoser, 2004), but in many parts, their populations are highly fragmented (Kaczensky et al., 2013). Thus, it is often challenging to determine the influence of social, ecological, and evolutionary constraints on genetic relatedness and family genetic structure, as this ideally requires high-resolution genetic data from a continuous and unfragmented population.
For the lynx family, studies on kin clusters and philopatry have given inconclusive results. A link between kin structure and dispersal has been found for bobcat (Lynx rufus) (Croteau, Heist, & Nielsen, 2010;Janečka et al., 2006), but not for Canada lynx (Lynx canadensis) (Campbell & Strobeck, 2006). In Sweden, telemetry studies showed that about one third of Eurasian lynx female offspring remained philopatric (Samelius et al., 2011), indicating the potential for geographic clustering of female relatives. Schmidt, Davoli, Kowalczyk, and Ettore (2016) found an insignificant pairwise geographic genetic distance relationship in a small isolated population of lynx in Białowieza, Poland, but it remains unknown if these results can also be found in other, continuous populations on larger scales. In Latvia, a genetic approach was used to identify the number and location of family groups based on parent-offspring relationships, but details on family genetic structure were not provided (Bagrade et al., 2016).
In this study, we have investigated the genetic relatedness and family structure of Eurasian lynx females in southern Finland. The Eurasian lynx is the only felid species in Finland and the population has shown a substantial population recovery and range expansion during late 1990s and early 21st century (Chapron et al., 2014). The population estimate is based on family group counts, method modified from that in use in Scandinavia (Andrén et al., 2002;Linnell et al., 2007), and has increased from 1,100 to 2,700 adult individuals during years 2007 to 2015 (Luke, 2017), corresponding to an average yearly increase of 12% (Lambda = (2,700/1,100) 1/8 = 1.12).
Distribution area covers whole of Finland, with highest densities in south and central parts of the country (details in : Holmala, 2013)).
In this continuously distributed population, we analyzed the genetic family structure using 23 autosomal microsatellite markers and tested for possible spatial genetic family structure using a combination of Bayesian clustering, spatial autocorrelation, and forensic genetic parentage analysis. This allowed us to study how genetic relatedness and family structure are organized in space in female lynx unconstrained by low population size, isolation or fragmentation.

| Study area
Our study area encompasses two game management districts in the southeastern and eastern parts of Finland, which have a total area of 22,936 km 2 ( Figure 1). About 70% of the terrestrial area is covered by forest, where pine, spruce, and birch are the most common trees.
Roughly 10% of the region is covered by lakes. Together with mires, F I G U R E 1 Location of the study area, locations of the Eurasian lynx female tissue samples from years 2007 to 2015, and their identified clusters in southern Finland. Symbols denote the genetic cluster each sample was assigned to cluster 1 = gray triangle; cluster 2 = black square; cluster 3 = dark gray circle; A (admixed) = hollow diamond. N = 132 farm land, and urban areas, these two regions form a heterogeneous mosaic of different land uses.

| Family genetic structure, spatial autocorrelation, and genetic parentage analysis
We used STRUCTURE v.2.3.4 (Falush, Stephens, & Pritchard, 2003;Pritchard, Stephens, & Donnelly, 2000) to test for indications of family genetic clustering among the females. We set the maximum number of populations to 10 (K = 10) with 10 independent runs for each K and assuming population admixture and correlated allele frequencies.
Burn-in period was 100,000 Markov-Chain-Monte-Carlo (MCMC) iterations, with a subsequent sampling of 1,000,000 MCMC iterations. We processed the results with Structure Harvester (Earl & von Holdt, 2012), which implements the ad hoc approach of Evanno, Regnaut, and Goudet (2005), and determined the number of putative family genetic clusters, using a membership value of q ≥ 0.7 as a threshold value (Kopatz et al., 2014;Pelletier, Obbard, Mills, Howe, & Burrows, 2012). Based on cluster assignment, we performed a factorial correspondence analysis (FCA) among all genotypes, and calculated the number of alleles, observed and expected heterozygosity using GenAlEx 6.5 and inbreeding coefficient using Genetix 4.05.2 (Belkhir, Borsa, Chikhi, Raufaste, & Bonhomme, 1996-2004. We performed a spatial autocorrelation analysis using GenAlEx 6.501 (Peakall & Smouse 2006;Peakall & Smouse 2012) to examine the relationship between genetic and spatial distance among each pair of female lynx. In the analysis, the family genetic structuring indicated by STRUCTURE was taken into account by using a multipopulation approach. To strike the balance between sample size TA B L E 1 23 different STR markers and the XY-test used in the genetic analysis. Markers are ordered according to the multiplex set they were assigned to, including the respective annealing temperatures (AT) and fluorescent label (Flag). Final concentration of the respective markers was the same for both forward and reverse primer inferred if the 95% CI around r (the average relatedness) did not contain 0, and if r exceeded the 95% CI around the null hypothesis of r = 0, that is no spatial structure.
We used the program Familias 3.1.9.5 (Egeland & Mostad, 2000;Kling, Tillmar, & Egeland, 2014) to reconstruct parenthood and sib ship from the microsatellite data by calculating likelihood ratios (LRs) for mother-offspring and sibling relationship. The Familias software is used worldwide by forensic laboratories and has been applied to numerous cases, for example resolving family relations and individual identification after disasters. Familias is programmed to recognize several different relationships and makes a pairwise comparison with all individual genotypes against each other and calculates an LR for each relationship. This makes an objective way to distinguish between the most likely relationships. The LR represents the probability of hypothesis one (candidate female is the true mother) divided by the probability of hypothesis two (candidate female is unrelated to the offspring in question) (Marshall, Slate, Kruuk, & Pemberton, 1998).
Based on the a priori probability of being related or unrelated equals to 0.5, the LR shows which relationship is more likely than others. This means a LR of 20 from the genetic analysis corresponds to a 95% probability for relatedness and would be considered as significant support for the relationship in question, while a LR value above 100 would be considered as highly significant support (99% probability of relatedness). Calculating family relationships requires allele frequency data of the population in question and a kinship correction. These we deduced from the data of this study. The Familias software is used worldwide by forensic laboratories and has been applied to numerous cases, for example resolving family relations and individual identification after disasters. In a second step, we used the pairwise geographic distances between death locations of each pair of individuals to calculate the mean and median distance between female parent-offspring, siblings, and all individuals. We also expected that individuals within pairs identified as close relatives would assign to the same family cluster, which was checked against the STRUCTURE results.

| Microsatellite genotyping
We obtained a full genetic profile of 23 STRs for 132 samples (N = 133). One sample did not contain enough DNA for successful amplification with the employed STR markers. All of the samples genotyped successfully received a unique identity. For 13 samples, analysis of one or more markers had to be repeated due to amplification failure in the first round, when three tissue samples showed no amplification success in one or two markers.

| Genetic diversity
Genetic diversity was relatively high, and expected heterozygosity and observed heterozygosity were almost identical, both averaging 0.66 across all markers (

| Fine-scale genetic family structure
Based on the estimated mean likelihood values and Evanno's ΔK

| Genetic parentage analysis
Parentage analyses using the program Familias suggested close, genetic relationships between many different pairs of female lynx: parent-offspring, sibling or half-sibling, and these pairs of close kin were found both within all three genetic clusters and within the group of admixed individuals (Tables 3-5). Individuals within pairs of close kin consistently assigned to the same genetic cluster, except in some instances where one of the individuals in a pair was admixed. Among individuals having at least 2-year difference in age, 32 pairs (47 individuals) were identified as likely mother and daughter by significant support that is a LR above 20 (Table 3). Of these, 97% showed 99% probability of relatedness. Furthermore, 64 pairs (79 individuals) showed LR >20 for being siblings and of these, 68% showed 99% probability of relatedness. In total, 26 pairs (42 individuals) showed LR >20 in both of the analyses above and thus identified as both likely mother-daughter and likely sisters (Table 5). In all of these cases, LR was higher for motherdaughter than for sister, which was also supported by a relatively large age difference in most cases. Assuming that all of these were mother-daughter relationships, the total number of significant sister relationships would be reduced to 38 (57 individuals; F I G U R E 2 (a) Mean of estimated log-likelihood values, and (b) rate of log-likelihood values for 132 Eurasian lynx females for different number of clusters from the software STRUCTURE (Pritchard et al., 2000) postprocessed with Evanno's approach (Evanno et al., 2005) F Samples are grouped according to their harvest location and sorted by longitude from west (left) to east (right). The y-axis indicates the membership coefficient q, that is the likelihood of belonging to a particular cluster. Each bar represents one individual and the length of each section of one bar corresponds to the q value for the respective cluster Table 4), and the total percentage of individuals with at least one close kin (mother-daughter and sister-sister) would be 63.6% (84 ind.), with 36.4% (48 ind.) being unrelated to all others. The median geographic distance between identified mother-offspring pairs was 36.8 km (mean distance 112.6 km; mean age difference 5.3 years) and between siblings 55.6 km (mean distance 108.6 km; mean age difference 2.5 years). In comparison, the pairwise median geographic distance between all females was 228.6 km (mean distance 262.8 km; mean age difference 3.03 years; Figure 6).

| D ISCUSS I ON
Using three different genetic approaches, we found support for family genetic structure among females of a solitary carnivore, the Eurasian lynx. The potential family groups were represented by spatially overlapping clusters and a flat but negative pairwise geneticgeographic distance relationship. Furthermore, we found evidence that within each of the putative family genetic clusters, there were several different sibling or mother-daughter pairs. Thus, our results suggest that in a spatially unrestricted population of Eurasian lynx, closely related females tend to cluster together geographically, in agreement with the hypothesis that they may form matrilineal as- Also, no effective geographic or anthropogenic-associated barriers for lynx are known from southern Finland. Moreover, the fine-scale genetic clustering found in our study is on a much smaller spatial scale than would normally be considered relevant for studying genetic subpopulations, stretching only to the similar extent as the diameter of several lynx home ranges (K. Holmala unpublished, Linnell et al., 2007), and with siblings or mother-daughter pairs consistently occurring within clusters. Furthermore, Eurasian lynx F I G U R E 4 Visualization of the factorial correspondence analysis (FCA) for female lynx genotypes sampled in southern Finland in the time period from 2007 to 2015. Different colors represent the clusters identified by the Bayesian clustering approach: cluster 1 (black squares), cluster 2 (light gray triangles), cluster 3 (gray diamonds), and admixed individuals with a cluster membership value q < 0.7 (white circles) F I G U R E 5 Results of the spatial autocorrelation analysis, that is combined correlation between genetic and spatial distance with GenAlEx 6.501 of lynx samples from southern Finland in the time period 2005 to 2015. The estimated relatedness coefficient (r; y-axis) for each distance class (x-axis) is given as a solid line. The 95% confidence interval for the Null hypothesis of random distribution is given as a dashed line, the bootstrap errors are displayed as whiskers. N = 132 is capable of dispersing long distances, for example up to 215 km females and to 428 km males (Samelius et al., 2011), thus potentially allowing for strong gene flow across large areas and among subpopulations (Wayne & Koepfli, 1996) when populations persist through the landscape. In support of this, Ratkiewicz et al. (2014) found evidence of active gene transfer between Finnish and Russian lynx populations. Thus, overall it seems likely that the genetic clustering observed in our study is a signature of the many different pairs of close relatives identified within each cluster, but not between them.
The observed significant correlation of genetic and geographic distances, that is the shorter the distance between individuals, the higher the pairwise relatedness, further corroborates this interpretation. Some of the individuals not clearly assigned to clusters and thus characterized as admixed genotypes, had close relatives among those individuals assigned to a genetic cluster. This indicates mixing of the genetic groups and that those individuals had the same mother but probably a different father. Another possible explanation is that the admixed individuals originate from family groups located outside the study area.
The median distance between all females in our study was over four times longer than between siblings (55.6 km) and over six times more than between parent and offspring (36.8 km). There was a large variation in distances especially for siblings, whereas most mother-daughter pairs were relatively close (although some outliers were found). The cost of defending a given territory may increase with population density, resulting in increased home-range overlap (Rodgers et al., 2015). Under these conditions, it may be that related individuals tolerate the costs of sharing resources due to benefits gained from inclusive fitness (Anderson, 1989). Kin clusters that are associated with home-range overlap could potentially support higher local population densities also for Eurasian lynx. However, a strong kin cluster in an area may also potentially hinder immigrating unrelated female lynx from establishing new territories in the vicinity. Immigration by many species of territorial mammals and birds appears to be limited by crowding (Lambin, Aars, & Piertney, 2001).
TA B L E 4 Recognized 38 sibling relationships (60 individuals) in a lynx population in southern Finland during years 2007-2015. The likelihood ratios (LR) above 20 were sorted from the highest to the lowest. Individual identification contains ID-number, year of birth, year of death, and genetic cluster (1, 2, 3, and A = admixed) Inversely density-dependent dispersal, impeding both immigration and emigration, seems to be true for Eurasian lynx in fragmented populations (Zimmermann, 2004;Zimmermann, Breitenmoser-Würsten, & Breitenmoser, 2005, but whether or not this is also the case in unfragmented ones, such as in Finland, remains to be investigated. Our results support the hypothesis that Eurasian lynx form matrilineal assemblages at regional scale. However, Elbroch et al. (2015) TA B L E 5 Recognized 26 pairs of Mother-offspring and siblings in a lynx population in southern Finland during years 2007-2015. The likelihood ratios (LR) above 20 were sorted from the highest to the lowest. Individual identification contains ID-number, year of birth, year of death, and genetic cluster (1, 2, 3, and A = admixed)  (Fattebert et al., 2016).
When studying social organization based on genetic data, special attention should be given in selecting the right combination of methods, taking into account species, population size, dispersal capability, the degree of population fragmentation, and the spatial and temporal extent of the study. These may be the reasons why a previous study on a small, isolated population showed generally lack of a relationship between the spatial distance and relatedness among individuals, but on the other hand, showed the domination of the entire population by a limited number of reproducing individuals, which partly indicates kin clustering (Schmidt et al, 2016).
It could also mean that social structure is flexible and changes with external conditions. Philopatry, however, does not necessarily lead to genetic clustering (Biek et al., 2006). Biek et al. (2006) found that even though female pumas remained philopatric, there was no genetic clustering, although genetic legacy of females with high reproduction success could be traced. They assumed that either the females were not successful in leaving philopatric offspring or the males, which emigrated from more distant populations, brought enough different alleles to outweigh the clustering phenomenon.
It is also worthwhile to ask, is genetic clustering of related individuals always a proof of philopatry? The question relates to the ecologically meaningful spatial scale of the species and the scale used in the study. In species with a social organization determined by family bonds, the spatial distribution of related individuals and the level of genetic clustering we observed might potentially be higher in small, isolated, and fragmented populations as a result of lower dispersal possibilities. Indeed, subadults female lynx may be less prone to cross barriers such as highways and densely populated valleys than males (Zimmermann, Breitenmoser-Würsten, & Breitenmoser, 2005. However, the opposite effect is also not unthinkable based on the inversely density-dependent dispersal behavior of lynx (Zimmermann, Breitenmoser-Würsten, & Breitenmoser, 2005 and, moreover, the possibility that matrilineal structure might disappear at low population densities, as observed for leopard (Fattebert et al., 2016).
The whole Finnish lynx population is of native origin and historically, it experienced periods of population decline. However, it is continuous and well connected to the population in Russia via extensive woodlands (Chapron et al., 2014). As such, our study contributes reference values for genetic parameters from a large lynx population in an almost unfragmented habitat. Robust reference values from large lynx populations are required for the assessment of the genetic status, management, and remedy of the still small reintroduced lynx populations in central Europe (e.g., Zimmermann, 2004;Zimmermann, Breitenmoser-Würsten, & Breitenmoser, 2005;Zimmermann, Breitenmoser-Würsten, & Breitenmoser, 2007). Whether the patterns observed among female lynx in our study represent true matrilineal assemblages will be further clarified by including males in analysis. This will also lead to obtaining estimates for the genetic parameters of the whole population. Studies including also males are now needed to further assess the social organization of this species.

ACK N OWLED G M ENTS
We thank our numerous technicians and volunteers in the field and laboratory for their rigorous work on samples, and Ida Fløystad for establishing the microsatellite method in the laboratory. Work has been partly financed by Raija and Ossi Tuuliainen Foundation.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
KH, SBH, AH, AK, JS, and HGE designed the study. KH was responsible for collected samples. JS generated the raw data in the laboratory. AK, SBH, and JS analyzed the data. All authors contributed to the writing of the manuscript. All authors approved the final version of the manuscript.