\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[english]{babel}



\usepackage{float}







\iflatexml
% \documentclass[cmbright]{WileySTAT-V1} % Not (yet) supported by LaTeXML.
\fi
 %% [cmbright,doublespace], [demo] for draft mode
% \usepackage{moreverb} % Not (yet) supported by LaTeXML.

% \usepackage{gensymb} % Not (yet) supported by LaTeXML.

\usepackage{longtable}

\usepackage{array}

\usepackage{lineno}
% \usepackage[style=apa, backend=biber]{biblatex} % Not (yet) supported by LaTeXML.
\usepackage{xcolor}

% \usepackage{bookmark} % Not (yet) supported by LaTeXML.

% First macro on next line not (yet) supported by LaTeXML: 
% \addbibresource{kinship.bib}
\IfFileExists{upquote.sty}{% \usepackage{upquote} % Not (yet) supported by LaTeXML.
}{}
\def\volumeyear{*}

\def\volumenumber{*}

\def\DOI{*}
% First macro on next line not (yet) supported by LaTeXML: 
% \runninghead{McMaster et al.}{Kinship estimation in non-model species}
\linenumbers


\begin{document}

\title{Evaluating kinship estimation methods for reduced-representation SNP
data in non-model species}



\author[1]{Eilish McMaster}%
\author[2]{Patricia Lu-Irving}%
\author[3]{Marlien van der Merwe}%
\author[4]{Simon Ho}%
\author[5]{Maurizio Rossetto}%
\affil[1]{The University of Sydney School of Life and Environmental Sciences}%
\affil[2]{Botanic Gardens of Sydney}%
\affil[3]{Royal Botanic Gardens and Domain Trust}%
\affil[4]{University of Sydney}%
\affil[5]{The Royal Botanic Gardens Sydney}%


\vspace{-1em}



  
  \date{May 20, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
\internallinenumbers
Kinship estimation is widely used in ecological and evolutionary research, particularly in studies of human genealogy and genome-wide associations. In conservation, restoration, agriculture, and forestry, identifying relationships between individuals can be crucial for successful population management and can provide insight into inheritance patterns. Kinship estimation methods are typically designed for large datasets with hundreds of thousands of single-nucleotide polymorphisms. However, studies of non-model species often use much smaller datasets obtained using reduced-representation sequencing. To evaluate the performance of kinship estimation methods under these circumstances, we applied six algorithms to datasets from six non-model Australian flowering plant species (\textit{Acacia terminalis}, \textit{Acacia suaveolens}, \textit{Banksia serrata}, \textit{Banksia aemula}, \textit{Hakea sericea}, and \textit{Hakea teretifolia}), encompassing 3,390 individuals and 369 families. Our results show different performances of kinship methods on reduced-representation sequence data compared with prior evaluations. PC-Relate, RelateAdmix, and Goudet's beta dosage exhibited limited precision, KING Homo and KING Robust demonstrated high precision with limited sensitivity, while PLINK displayed variable sensitivity and precision. The sensitivity and precision of the methods were affected in various ways by filtering parameters; each method showed its best performance under different thresholds for minor allele frequency and locus missingness. We also present a case study that illustrates a practical application of the methods, demonstrating how estimates of kinship can inform management of seed production areas of the broadleaf hopbush (\textit{Dodonaea viscosa}). Based on our findings, we offer specific recommendations for utilizing kinship estimation methods in studies of reduced-representation sequence data from non-model species.%
\end{abstract}%



\sloppy


\address{\affilnum{a}{School of Life and Environmental Sciences, University of Sydney, 
    New South Wales,
    Australia} \\
\affilnum{b}{Research Centre for Ecosystem Resilience,
    Botanic Gardens of Sydney,
    New South Wales,
    Australia}}

\corremail{eilish.mcmaster@botanicgardens.nsw.gov.au; emcm4052@uni.sydney.edu.au}

\section*{Introduction}\label{introduction}

Understanding the relationships among individuals can provide valuable insights into inheritance patterns, genetic diversity, and evolutionary processes within populations. Kinship estimation algorithms play a crucial role in this context, by assessing relatedness based on genetic data. These estimators find widespread application in human contexts, such as genealogy for identifying unknown relatives and in genome-wide association studies to mitigate family bias \parencite{chang_second-generation_2015}. Applications extend to non-model organisms where there is no prior knowledge of relationships or population history. This is particularly relevant to conservation \parencite{bragg_optimizing_2020} restoration \parencite{breed_potential_2019,rossetto_restore_2019}, agriculture \parencite{garavito_identification_2016,baloch_whole_2017, ndjiondjop_genetic_2017, edet_dartseq-based_2018, odoi_genetic_2023}, and forestry \parencite{hall_parentage_2020,nantongo_silicodart_2022}, where knowledge of relationships between individuals is important for \emph{in situ} management of populations and breeding programs. However, kinship estimation methods have primarily been tested on large datasets from humans \parencite{manichaikul_robust_2010, ramstetter_benchmarking_2017}, well-studied animals \parencite{goudet_how_2018}, or simulations \parencite{moltke_relateadmix_2014}, which are not representative of the data typically used in studies of non-model species.

In population genomic analyses of non-model organisms, resource constraints often limit the accessibility of high-fidelity genomic methods such as resequencing or custom SNP arrays. The emphasis shifts towards adequately representing the study species using cost-effective methods such as reduced-representation sequencing, including RADseq \parencite{davey_radseq_2010, davey_genome-wide_2011}, ddRADseq \parencite{peterson_double_2012}, and DARTseq \parencite{sansaloni_high-density_2010}, to generate data from many samples rather than obtaining higher-quality data from few individuals. However, these methods tend to produce smaller datasets with thousands to tens of thousands of SNPs \parencite{baloch_whole_2017, alam_ultra-high-throughput_2018, rutherford_evolutionary_2022, lu-irving_capturing_2023} rather than hundreds of thousands or millions of SNPs \parencite{gibbs_international_2003, ha_coverage_2014, gudmundsson_variant_2022}. Reduced-representation sequencing methods are also more error prone \parencite{roques_microsatellites_2019, crotti_population_2020} than resequencing \parencite{stoler_sequencing_2021} or SNP arrays \parencite{ceballos_assessing_2018}. Reduced-representation sequencing approaches can offer high quality insights into population structure and gene flow, but there is less certainty about the impact of having fewer SNPs and higher error rate on more intensive analyses, such as kinship determination. Kinship estimation methods are usually designed with large SNP datasets in mind \parencite{ramstetter_benchmarking_2017}, so they might not perform as well in analyses of smaller datasets \parencite{goudet_how_2018}.

The biological characteristics of non-model species, particularly plants, mean that population structure and levels of inbreeding can vary dramatically even within species. For example, monoecious plants, which aren not constrained by single-sex boundaries, can experience extreme levels of inbreeding by biparental, cross-generational, and self fertilization, which is not often observed in animals \parencite{jain_evolution_1976}. Moreover, certain species demonstrate a high tolerance or preference for inbreeding, while others are more inclined towards outcrossing or hybridization \parencite{whitehead_plant_2018}. Even within species, mating system preference can vary between seasons, subpopulations \parencite{whitehead_plant_2018}, and maternal lines \parencite{lu-irving_capturing_2023}. Clonality due to apomixis or asexual budding can also be prevalent in some species. Importantly, these traits explicitly violate the assumptions of homogeneity and Hardy-Weinberg equilibrium often relied upon by kinship estimation methods.

Relationships between individuals are commonly estimated through the calculation of the kinship coefficient (\(\phi\)\textsubscript{ij}, between individuals \emph{i} and \emph{j} ), which relies upon the fundamental concepts of identity by state (IBS: whether two individuals share the same alleles at a locus) and identity by descent (IBD: whether two individuals share the same alleles at a locus because they share a common ancestor) \parencite{powell_reconciling_2010}. The kinship coefficient represents the probability that alleles randomly sampled from the same locus in two individuals are IBD. The degree of relationship between individuals can be directly inferred from the kinship coefficient due to the constraints of inheritance (i.e., a first-degree relative should be kin\textasciitilde0.25, second degree \textasciitilde0.125) \parencite{manichaikul_robust_2010}.

Two of the most common assumptions in kinship estimation models are population allele frequencies and Hardy-Weinberg equilibrium. Ideally, the allele frequencies of the ancestral founder population would be employed; however, in the absence of known pedigrees, kinship estimation methods commonly use the allele frequencies of the current population as a proxy when calculating the probability of IBD \parencite{powell_reconciling_2010}. A notable challenge arises from the finite nature of datasets, where assumed allele frequencies might not be representative of the larger population \parencite{purcell_plink_2007}. Furthermore, allele frequencies can differ between groups in a structured population, confounding IBD probability estimation. This can result in the overestimation of kinship, especially in smaller structure groups with distinct allele frequencies \parencite{weir_unified_2017}. Additionally, many kinship estimators assume Hardy-Weinberg equilibrium to derive IBD probabilities. However, the requirements of Hardy-Weinberg equilibrium, including a homogeneous population with random mating, are seldom met in natural populations due to factors like mate choice, inbreeding, or limitations in population connectivity. Various kinship estimation methods have been developed to address, correct, or circumvent these problematic assumptions, aiming for more robust and accurate kinship inference.
To evaluate the performance of kinship estimation methods in analyses of non-standard data, we conducted extensive testing using DArTseq data from six species of Australian flowering plants commonly used in habitat restoration practices: \emph{Acacia terminalis} , \emph{Acacia suaveolens} , \emph{Banksia serrata} , \emph{Banksia aemula} , \emph{Hakea sericea} , and \emph{Hakea teretifolia} . These data include known mother-offspring pairs, and comprise a total of 3,390 individuals and 369 families. Our study aims to: (i) determine which methods are most effective in detecting relationships while minimizing the occurrence of false positives; (ii) assess the impact of different filtering parameters on the sensitivity and precision of kinship estimation methods and identify optimal filtering parameters; and (iii) investigate the impacts of individual species biology on the performance of kinship estimation methods.

We apply our analytical insights to a case study of the broadleaf hopbush (\emph{Dodonaea viscosa} ), a cosmopolitan species found across Australia, New Zealand, Asia, Africa, and the Americas \parencite{harrington_species_2009}. The dataset comprises wild individuals from New South Wales, and three seed production areas (SPAs) used to make seed for landscape restoration. We analyse kinship in these SPAs to help manage genetic diversity and minimize inbreeding in the progeny. This case demonstrates how kinship estimation enhances species management by effectively regulating kin relationships, thereby preserving genetic diversity in offspring.

\section*{Methods}\label{methods}

\subsection*{\texorpdfstring{Data acquisition }{Data acquisition }}\label{data-acquisition}

To assess kinship estimation methods, we utilized data from six widely distributed plant species: \emph{Acacia terminalis, Acacia suaveolens, Banksia serrata, Banksia aemula, Hakea sericea,} and \emph{Hakea teretifolia.} These species belong to the angiosperm families Proteaceae (\emph{Banksia} and \emph{Hakea} ) and Mimosaceae (\emph{Acacia} ) and are indigenous to Australia.

Target species were sampled along the east coast of Australia, primarily in New South Wales, with some extending to Queensland (\emph{B. serrata} ) and Victoria and Tasmania (\emph{A. terminalis} ) (Figures S1-S18). At the time of field sampling, approximately 3 g of leaf tissue, and seed were collected. The geographic coordinates (latitude and longitude) of each plant were recorded. Seed was germinated \emph{ex situ} at the Australian PlantBank, and leaf tissue was taken from these offspring at 4--6 weeks. The leaf samples were stored at -80 \degree C for a minimum of 12 h, then freeze dried and stored in silica gel to ensure preservation until DNA extraction. The \emph{Hakea} datasets were described in a previous publication \parencite{lu-irving_capturing_2023}.

Genotyping was conducted using DArTseq, a reduced-representation sequencing method implemented by Diversity Arrays Technology Australia Pty Ltd (Canberra, Australia). The DArTseq method involves a genome restriction digest followed by sequencing of the digested products, with single-nucleotide polymorphisms (SNPs) called using proprietary analytical pipelines by DArT Pty Ltd. Genotype data were provided as SNP matrices along with quality informing statistics for each bi-allelic marker and sample \parencite{andrews_harnessing_2016}.

\subsection*{Preliminary assessment of population structure}\label{preliminary-assessment-of-population-structure}

We assessed the population structure of each species by analysing genotype data through PCA, $F_{ST}$, and STRUCTURE analyses. This is important as population structure considerably influences the performance of kinship estimation methods. We imported the DArT SNP data for each species into RStudio version 2023.9.1.494 \parencite{rstudio_team_rstudio_2020}, running R version 4.3.1 \parencite{r_core_team_r_2021}. The data were then filtered using RRtools v0.1 \parencite{bragg_optimizing_2020} to exclude loci with \textless{} 96\% reproducibility, \textgreater{} 10\% missing data (locus missingness), and minor allele frequency \textless{} 5\%. Additionally, we removed fixed loci and retained only one SNP per DArT tag. Samples with \textgreater{} 30\% missing data (sample missingness) were also excluded.

In R, we conducted PCA using the \texttt{glPca} function from adegenet v2.1.20 \parencite{jombart_adegenet_2008}, and $F_{ST}$ analysis used \texttt{snpgdsFst} from SNPRelate v1.36 \parencite{zheng_high-performance_2012}. The \texttt{struct} file for STRUCTURE was generated by \texttt{dart2struct} from RRtools. For each species, STRUCTURE v2.3.4 \parencite{pritchard_inference_2000} ran with 10,000 burn-in steps and 20,000 MCMC steps for ancestral populations (\emph{k} ) from 1 to 10, with 5 replicates per \emph{k} . The Evanno method in StructureHarvester v0.6.63 provided summary statistics for each \emph{k} . Results were visualized in RStudio using pophelper v2.3.1 \parencite{francis_pophelper_2017}, ggplot2 v3.4.4 \parencite{wickham_ggplot2_2016}, ggthemes v4.2.4, scatterpie v0.2.1, and ozmaps v0.4.5 (Figures S1-S18).\selectlanguage{english}
\begin{table}[ht]
\begin{small}\selectlanguage{english}
\begin{longtable}[]{@{}
  >{\raggedright\arraybackslash}p{(\columnwidth - 22\tabcolsep) * \real{0.0968}}
  >{\raggedright\arraybackslash}p{(\columnwidth - 22\tabcolsep) * \real{0.0993}}
  >{\raggedright\arraybackslash}p{(\columnwidth - 22\tabcolsep) * \real{0.1753}}
  >{\raggedright\arraybackslash}p{(\columnwidth - 22\tabcolsep) * \real{0.0514}}
  >{\raggedright\arraybackslash}p{(\columnwidth - 22\tabcolsep) * \real{0.0591}}
  >{\raggedright\arraybackslash}p{(\columnwidth - 22\tabcolsep) * \real{0.0612}}
  >{\raggedright\arraybackslash}p{(\columnwidth - 22\tabcolsep) * \real{0.0661}}
  >{\raggedright\arraybackslash}p{(\columnwidth - 22\tabcolsep) * \real{0.0689}}
  >{\raggedright\arraybackslash}p{(\columnwidth - 22\tabcolsep) * \real{0.0590}}
  >{\raggedright\arraybackslash}p{(\columnwidth - 22\tabcolsep) * \real{0.0926}}
  >{\raggedright\arraybackslash}p{(\columnwidth - 22\tabcolsep) * \real{0.0221}}
  >{\raggedright\arraybackslash}p{(\columnwidth - 22\tabcolsep) * \real{0.1482}}@{}}
\caption{{% \internallinenumbers
Overview of genomic datasets and population structure parameters under minimal filtering parameters (maximum locus missingness of 0.4, maximum sample missingness of 0.3, and minimum minor allele frequency of 0.00001).
}}\label{table:species_summary}\\
\toprule\noalign{}
\begin{minipage}[b]{\linewidth}\raggedright
Species
\end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright
Code
\end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright
Common name
\end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright
Loci$^a$
\end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright
Families$^b$
\end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright
Mothers
\end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright
Seedlings
\end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright
Unrelated$^c$
\end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright
Total sites
\end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright
Sites with families$^d$
\end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright
\emph{k}$^e$
\end{minipage} & \begin{minipage}[b]{\linewidth}\raggedright
DNA 1C (pg)$^f$
\end{minipage} \\
\midrule\noalign{}
\endhead
\bottomrule\noalign{}
\endlastfoot
\emph{Hakea teretifolia} & HakeTere & Dagger hakea, needlebush & 46,761 & 66 & 65 & 302 & 169 & 32 & 7 & 2 & 0.93 \parencite{morgan_relationship_2005} \\
\emph{Hakea sericea} & HakeSeri & Bushy needlewood, silky hakea & 27,627 & 61 & 60 & 310 & 152 & 29 & 6 & 5 & - \\
\emph{Acacia suaveolens} & AcacSuav & Sweet wattle & 28,269 & 40 & 30 & 157 & 210 & 47 & 14 & 3 & 0.86 \parencite{gallagher_invasiveness_2011} \\
\emph{Acacia terminalis} & AcacTerm & Sunshine wattle & 18,782 & 76 & 64 & 339 & 404 & 88 & 13 & 2 & 0.86 \parencite{gallagher_invasiveness_2011} \\
\emph{Banksia serrata} & BankSerr & Saw banksia, old man banksia, wiriyagan & 6,780 & 90 & 87 & 477 & 176 & 44 & 14 & 2 & 0.92 \parencite{jordan_environmental_2015} \\
\emph{Banksia aemula} & BankAemu & Wallum banksia & 7,583 & 36 & 32 & 281 & 75 & 14 & 6 & 2 & - \\
\hline
\end{longtable}
\footnotesize{$^a$total number of polymorphic loci, $^b$groups of related individuals, $^c$wild individuals with sequenced tissue but no collected seed, $^d$sites with associated ex situ germinated seedlings, $^e$optimal number of ancestral populations, $^f$known DNA 1C values for study species are also listed \parencite{leitch_angiosperm_2019}. All species with known C-values are confirmed diploid.}\\
\end{small}

\end{table}

\subsection*{Data filtering}\label{data-filtering}

An aim of this study was to assess the impact of various filtering parameters on the effectiveness of different kinship estimation methods, as filtering can have a large impact on the outcomes of SNP analyses \parencite{shafer_bioinformatic_2017}. To accomplish this, we filtered the SNP data of each species by employing combinations of maximum locus missingness values (0.00001, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, and 0.4), maximum sample missingness values (0.1, 0.2, and 0.3), and minimum minor allele frequency values (0.00001, 0.01, 0.03, 0.05, 0.1, 0.15, and 0.2). The resulting filtered SNP datasets were converted into PLINK bed files using the \texttt{snpgdsGDS2BED} function and GDS files using the \texttt{snpgdsCreateGeno} function from SNPrelate.

\subsection*{Kinship methods}\label{kinship-methods}

We compared six established methods for estimating kinship, each using different approaches. We applied these methods to all filtered datasets. These methods focus on SNP-based, not segment-based IBD estimation, as our data lacks reference genomes for mapping. We chose them for their diverse computational approaches in estimating kinship and handling population structure.

While PC-Relate and RelateAdmix showed superior performance in human datasets according to previous benchmarks \parencite{ramstetter_benchmarking_2017}, many non-model species analyses still prefer older PLINK and KING methods. Hence, we aimed to assess the suitability of PC-Relate and RelateAdmix for smaller datasets. Finally, Goudet's method is more recent and producing promising results \parencite{ochoa_estimating_2021}, but its performance on real non-model species data remains uncertain.

\subsubsection*{PLINK}\label{plink}

PLINK employs a method-of-moments strategy to estimate the likelihood of sharing 0, 1, or 2 alleles (\emph{k} \textsubscript{0}, \emph{k} \textsubscript{1}, \emph{k} \textsubscript{2}) at any given locus due to IBD between individuals \parencite{purcell_plink_2007,chang_second-generation_2015}. This estimation relies on population allele frequencies derived from the input data, assuming Hardy-Weinberg equilibrium. The resulting IBD probabilities for pairs of individuals can then be used to calculate the pairwise kinship coefficient \((\phi = \frac{k1}{4} + \frac{k2}{2})\). We executed PLINK v1.9 on each bed file using the following command: \texttt{{[}plink -\/-bfile filtered\_bed\_file -\/-genome -\/-out output\_file{]}}.

\subsubsection*{KING Homo and KING Robust}\label{king-homo-and-king-robust}

The KING Homo method employs a unique approach to kinship estimation by bypassing the computation of IBD probabilities \parencite{manichaikul_robust_2010}. In contrast to PLINK, which relies on IBD probabilities for kinship derivation, KING Homo estimates the kinship coefficient (\(\phi\)) based on the genetic distance between individuals, determined by the number of shared heterozygous and homozygous loci. This genetic distance is modelled under the assumption of Hardy-Weinberg equilibrium. KING Homo maintains the assumption of population homogeneity. KING Robust applies an additional correction within families identified in the \texttt{fam} file to account for changes in allele frequencies within families.

We ran KING v2.3.2 on each bed file with the command (applied identically for KING Homo and KING Robust): \texttt{{[}king -b filtered\_bed\_file -\/-kinship -\/-prefix output\_file{]}}. In the case of KING Homo, the corresponding \texttt{fam} file featured unique families for each individual. Conversely, the \texttt{fam} file for KING Robust assigned individuals to families based on site information.

\subsubsection*{RelateAdmix}\label{relateadmix}

RelateAdmix uses maximum likelihood to estimate IBD probabilities while accounting for population structure \parencite{moltke_relateadmix_2014}. To do this, it uses ancestry estimation provided by ADMIXTURE software \parencite{alexander_fast_2009} to readjust the expectations of individual-specific allele frequencies based upon the ancestral background. The individual allele frequencies are then used to calculate the IBD probabilities. In the case of no admixture, the method is equivalent to that of \parencite{thompson_estimation_1975, milligan_maximum-likelihood_2003}.

To create the input ancestry estimation files for RelateAdmix, we ran ADMIXTURE linux-v1.3.0 \parencite{alexander_fast_2009} on each filtered \texttt{bed} file with the optimal \emph{k} value for each species (as determined by STRUCTURE using the Evanno method): \texttt{{[}admixture filtered\_bed\_file k{]}}. The ancestral population frequencies and admixture proportion files were then used to run RelateAdmix: \texttt{{[}relateAdmix -plink filtered\_bed\_file -f admixture.P -q admixture.Q -o output\_file{]}}.

\subsubsection*{PC-relate}\label{pc-relate}

PC-relate addresses challenges related to kinship calculation in the presence of population structure by adjusting individual allele frequency expectations \parencite{conomos_model-free_2016}, akin to RelateAdmix. Unlike RelateAdmix, which relies on ADMIXTURE input for defining structure, PC-relate utilizes the results of a principal components analysis in related samples (PC-AiR) to characterize structure. PC-AiR differs from conventional PCA methods in that it incorporates genotype and KING kinship data. In the PC-relate approach, the principal components of a PC-AiR derived from genotype and KING input are used to adjust individual-specific allele frequency expectations.

For each filtered dataset, we conducted a PC-AiR using the \texttt{pcair} function from the GENESIS v2.20.1 R package \parencite{gogarten_genetic_2019} using the GDS and KING Homo kinship files as input. We then ran the \texttt{pcrelate} function on the GDS object with the first two PCs of the PC-AiR to obtain kinship values.

\subsubsection*{Goudet's beta dosage}\label{goudets-beta-dosage}

Goudet et al. (2017, 2018) introduced a novel kinship estimator (denoted as \(r_{ij}^{\beta}\) or beta dosage) designed to address challenges like inbreeding, founder effects, and population structure, which are particularly relevant for species undergoing conservation management. Unlike methods relying on population allele frequencies for kinship calculation, this approach considers the kinship relative to all pairs in the sample.

Because of the calculation's simplicity, we reproduce it here. The kinship estimator is formulated as \(r_{ij}^{\beta} = \frac{M_{ij} - M_{S}}{1 - M_{S}}\), where \(M_{ij}\) represents the proportion of alleles in individuals \emph{i} and \emph{j} that are identical by state, and \(M_{S}\) stands for the average proportion of matching alleles for all pairs of individuals. Notably, the relative nature of this estimator implies that the average kinship in the population is 0. Individuals sharing fewer alleles than the average exhibit negative kinship values. We calculated Goudet's kinship using the \texttt{beta.dosage} function from hierfstat v0.5-11 \parencite{goudet_hierfstat_2005, goudet_how_2018} on the filtered GDS data.

\subsection*{Evaluation of sensitivity and srecision in kinship estimation}\label{evaluation-of-sensitivity-and-precision-in-kinship-estimation}

We classified predicted degrees of relationship between individuals by the pairwise kinship results for all methods using the theoretical thresholds established by \parencite{manichaikul_robust_2010} (Table S1). To evaluate the performance of the kinship estimation methods, we assessed the sensitivity and precision of each run. The calculation of these metrics involves categorizing kinship estimates as true positives (TP), false positives (FP), and false negatives (FN). However, our kinship dataset differs from others in that it lacks complete pedigree data. Having only one known parent for the seedlings introduces complexities such as potential shared parents, selfing, or apomixis. To tackle this issue, we classified relationships among seedlings that share the same mother as  true positives if they were second-, first-, or zero-degrees by kinship. Likewise, for parent-offspring pairs, those categorized as first- or zero-degree were regarded as true positives. Lower degrees of relationship inferred for known relatives were classified as false negatives.

To categorize inferred relationships as false positives, we had to make assumptions about which pairwise relationships were improbable in our data. Pollen and seed are generally expected to disperse short distances most of the time, and long-distance dispersal is rare \parencite{webb_selection_1998, nathan_spatial_2000, wessinger_pollen_2021}. Published evidence from empirical studies of pollination and seed dispersal in the study species and congeners are consistent with this expectation, with pollen and seed dispersing \textless{} 3 km \parencite{hackett_pollination_2001, whitney_dispersal_2002, he_long-distance_2009, krauss_contrasting_2009, groom_implications_2010, mast_time-calibrated_2012, millar_extensive_2014, palfi_seed_2017}. Thus, we deemed inferred zero-, first-, or second-degree relationships between individuals \textgreater{} 10km apart to be false positives. All other inferred relationships, such as those within sites without known parentage or between individuals \textless{} 10 km apart, were disregarded for sensitivity and precision calculations due to the uncertainty surrounding their true relatedness.

We then calculated different sensitivity and precision metrics outlined below for each kinship output. We chose sensitivity and precision because they are widely used to describe the efficacy of tests and are intuitive. Sensitivity is the capacity of a test or method to accurately detect true positives, indicating its effectiveness in identifying the existence of a condition when it is truly present. For example, in a kinship analysis, sensitivity indicates how accurately the method identifies real family relationships among individuals. Precision, on the other hand, measures the proportion of positive results that are correctly determined. Precision tells us how confident we can be that a positive result truly indicates the presence of the condition being tested for. In kinship analysis, precision would indicate the reliability of the method in identifying at correctly identifying family relationships among those that it flags as related.

In an empirical scenario, the desired balance between sensitivity and precision depends on the specific goals and constraints of the study. If the primary concern is to avoid missing any true relationships (i.e., false negatives), then sensitivity becomes more important. On the other hand, if minimizing false positives is critical to ensure the accuracy of results, then precision takes precedence.

\[
sensitivity\ (overall) = \frac{{TP}_{All}}{{TP}_{All} + FN}
\]

\[
precision\ (overall) = \frac{{TP}_{All}}{{TP}_{All} + FP}
\]

\[
sensitivity\ (parent\ offspring) = \frac{{TP}_{PO}}{{TP}_{PO} + {FN}_{PO}}
\]

\[
\ precision\ (parent\ offspring) = \frac{{TP}_{PO}}{{TP}_{PO} + {Total}_{0,1}}
\]

\({Total}_{0,1}\)= all relationships classified as zero or first-degree, including TP sibship and FP.

\[
\ sensitivity\ (sibling) = \frac{{TP}_{Sib}}{{TP}_{Sib} + {FN}_{Sib}}
\]

\[
\ precision\ (sibling) = \frac{{TP}_{Sib}}{{TP}_{Sib} + {Total}_{0,1,2}}
\]

\({Total}_{0,1,2}\ \)= all relationships classified as zero-, first-, or second-degree, including TP parent-offspring relationships and FP.

\subsection*{Modelling sensitivity and precision}\label{modelling-sensitivity-and-precision}

To explore the relationship between filtering criteria (maximum locus missingness, maximum sample missingness, and minimum minor allele frequency) and the kinship estimates (overall sensitivity and precision), we employed random forest (RF) models for each kinship estimation method within each species. We selected RF models because they handle non-linear relationships well and estimate variable importance. Utilizing the \texttt{randomForest} function from RandomForest v4.7-1.1 \parencite{liaw_classication_2002}, we generated models for sensitivity (\texttt{sensitivity \textasciitilde{} locus\_missingness + sample\_missingness + maf}) and precision (\texttt{precision \textasciitilde{} locus\_missingness + sample\_missingness + maf}) for each kinship method within each species, specifying ntrees = 5,000. Subsequently, we extracted partial dependence plots for each predictor variable from each model using the \texttt{pdp} function from pdp v0.8.1 \parencite{greenwell_pdp_2017}.

\subsection*{Algorithm performance in diverse scenarios}\label{algorithm-performance-in-diverse-scenarios}

To assess the performance of methods under adverse conditions, we examined specific subsets of individuals within the structured \emph{Acacia suaveolens} and \emph{Hakea sericea} datasets. These subsets represent more complex scenarios that pose challenges for kinship algorithms. The individuals included were either from highly admixed (individuals with the largest ancestral population proportion below 75\%), from a large genetic structure group (over 90\% ancestry from a majority structure group in the species), or a small structure group (over 90\% ancestry from a minority structure group in the species; Table S4). We plotted the kinship values within sites that were known parent-offspring, known sibling, or unknown relatives by structure group.

\subsection*{Case study: \emph{Dodonaea viscosa} seed production area}%\label{case-study-dodonaea-viscosa-seed-production-area}

Seed production areas (SPAs) are areas where plants are managed to produce diverse and viable seed that can be harvested with the purpose of targeted revegetation strategies \parencite{gibson-roy_florabank_2021}. Seeds are harvested en mass and often used across large areas that have been cleared or degraded. It is important that SPAs minimize kinship and optimize genetic diversity representativeness to avoid inbreeding and maximized adaptive potential. Utilizing reduced-representation data, we assess the relatedness among individuals across three SPAs dedicated to \emph{Dodonaea viscosa} . This analysis seeks to establish whether adjustments need to be made within SPAs to lower the likelihood of mating between closely related individuals, which can diminish the quality of the seeds produced. Such investigations provide practical insights into how kinship assessments can inform on-ground applications.

The seeds harvested from these areas primarily support restoration efforts in the region bordering the Murray River, spanning New South Wales and Victoria (Australia). In our analysis, a total of 486 individuals were examined, comprising 305 from wild populations, 53 from SPA-A, 77 from SPA-B, and 51 from SPA-C. These figures represent approximately 15\%, 80\%, and 100\% of the respective SPA populations.

Using the same packages as above, we filtered the data for a minimum reproducibility of 96\%, fixed SNPs, and one SNP per locus. For the PCA, STRUCTURE, $F_{ST}$, and phylogenetic network analyses, we further filtered the data for 20\% locus missingness and 20\% sample missingness, resulting in 2,003 loci after filtering. The phylogenetic network nexus file was produced by transforming the genotype data to Euclidean distances using the \texttt{dist} function in R and writing to a file using the \texttt{splitstree} function from the RSplitsTree v 0.1.0 package. The nexus file was reimported using the \texttt{read.nexus.networx} function from the phangorn v 2.11.1 package and plotted using the tanggle v1.8.0 package.

Based on the outcomes of our kinship method comparison, we chose to assess \textit{D. viscosa} kinship using PLINK and KING Robust implemented in SNPRelate. In PLINK analyses, the data underwent filtering for 20\% locus missingness and 10\% minimum minor allele frequency. PLINK was run twice: once on the whole dataset, and once on the separate subgroups identified throughout the investigation. KING Robust was run on data filtered for a maximum of 10\% locus missingness with no minimum minor allele frequency.

\section*{Results}\label{results}

\subsection*{Species population structure}\label{species-population-structure}

\emph{Banksia aemula} , \emph{Banksia serrata} , \emph{Hakea teretifolia} , and \emph{Acacia terminalis} exhibited no significant population structuring, with the optimal \emph{k} value being 2, a result of isolation by distance. Among these, \emph{Banksia aemula} and \emph{Banksia serrata} represent ideal wild populations characterized by high gene flow between all sites. \emph{Hakea teretifolia} has a stronger pattern of isolation by distance, as does \emph{Acacia terminalis} which was collected over a range of 2,000 km (Figures S1-S18).

\emph{Acacia suaveolens} and \emph{Hakea sericea} introduce further complexity. \emph{Acacia suaveolens} has three distinct genetic groups (\emph{k} = 3) with limited gene flow between them, with one genetic group with a notably restricted range (Figures S4-S6). \emph{Hakea sericea} is highly structured with five distinct groups (\emph{k} = 5), experiencing limited gene flow even between proximate sites (Figures S14-S15). Additionally, \emph{H. sericea} demonstrates an exceptionally high degree of inbreeding within sites due to a preference for selfing \parencite{lu-irving_capturing_2023}.

\subsection*{Overall performance of methods}\label{overall-performance-of-methods}

Goudet's beta dosage method exhibited high sensitivity (\textgreater{} 0.8) across all species and filtering scenarios but demonstrated consistently low precision (\textless{} 0.25) (Figure \ref{fig:ppv_sens}). Notably, Goudet's precision improved in the analyses of the \emph{Banksia} datasets, which are homogeneous and have a higher proportion of related individuals than the other species (Table \ref{table:species_summary}). However, precision significantly decreased when families were randomly subset to three individuals (Figure S28).

Both KING methods demonstrated high precision, approaching 1 in most scenarios (\textgreater{} 0.9; Figure \ref{fig:ppv_sens}). However, they were conservative, exhibiting the lowest overall sensitivity compared with other methods (sensitivity typically \textless{} 0.5), although KING Robust generally had higher sensitivity than KING Homo (Figure \ref{fig:ppv_sens}). Both KING methods had the highest sensitivity in analyses of \emph{H. sericea,} which is highly inbred and structured, and the lowest sensitivity for the \emph{Banksia} species and \emph{A. terminalis} , which were more outbreeding homogeneous datasets.

PC-Relate demonstrated high sensitivity (\textgreater{} 0.9) across all filtering scenarios, with varying precision. Precision was relatively high (\textgreater{} 0.7) in the analyses of the \emph{Banksia} species, while it ranged from 0.9 to 0.1 for \emph{A. terminalis} and \emph{H. teretifolia} . In the analyses of species with more complex population structure, \emph{A. suaveolens} and \emph{H. sericea} , overall precision was consistently low (\textless{} 0.25) (Figure \ref{fig:ppv_sens}).

PLINK demonstrated variable sensitivity and precision, generally maintaining high precision across various filtering scenarios. However, precision decreased in certain cases, and this was particularly noticeable in the analysis of the \emph{A. terminalis} dataset (Figures \ref{fig:pdp_ppv}, S25). Sensitivity was variable for most datasets but generally approached 1 under optimal filtering scenarios (Figure \ref{fig:ppv_sens}).

RelateAdmix showed similar patterns to PC-Relate, with high sensitivity (\textless{} 0.7) but variable precision (0--1) between species datasets. Both methods displayed low precision in analyses of structured species (Figure \ref{fig:ppv_sens}).\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/ppv-sens-plot2/ppv-sens-plot2}
\caption{{\label{fig:ppv_sens}     \internallinenumbers
    Overall precision and overall sensitivity for each kinship method across diverse filtering criteria and various species. Each point on the plot represents the performance outcome of a specific combination of filtering criteria, including maximum locus missingness (LM), maximum sample missingness (SM), and minimum minor allele frequency (MAF). The data points are colour-coded based on the respective species. The prominent red point in each panel designates the specific run chosen for closer examination in Figures \ref{fig:kin_hist_po} and \ref{fig:kin_hist_sib}. One set of filters was designated per method, chosen because they produced superior outcomes across all species for the corresponding kinship method (Figures S19-S24). Information about these selected filtering parameters is presented in the bottom panel of each column corresponding to the respective method.%
}}
\end{center}
\end{figure}

\subsection*{Performance under optimal filters}\label{performance-under-optimal-filters}

In our datasets, parent-offspring relationships can be classified as either zero-degree (due to apomixis) or first-degree relationships. Sibling relationships can be zero (clones), first (full siblings), or second degree (stepsiblings). While apomixis is possible, it is expected to be less common than sexual reproduction. Additionally, half-siblings are more likely than full siblings in wild plants with open mating. Although the true proportions of zero-, first-, or second-degree relationships within the dataset are unknown, we can qualitatively assess the distribution of kinship values using the methods.

Goudet's beta dosage consistently appears to overestimate kinship between parent and offspring, with a notably higher proportion of these relationships exceeding the zero-degree threshold (Figure \ref{fig:kin_hist_po}). The same trend is observed for sibship, where an improbable proportion is designated as zero-degree (Figure \ref{fig:kin_hist_sib}).

In contrast, KING Homo considerably underestimates both parent-offspring (Figure \ref{fig:kin_hist_po}) and sibling (Figure \ref{fig:kin_hist_sib}) relationships, with many relationships going undetected. KING Robust produces similar results but is more sensitive to all known relationships.

PC-Relate also tends to overestimate parent-offspring kinship, with a large proportion of relationships being zero-degree, although to a lesser extent than Goudet. It is also sensitive to sibship but generally detects a higher proportion of zero- and first-degree sibling relationships than PLINK, RelateAdmix, or KING. Precision is typically higher with this method than with Goudet, except in the analysis of \emph{H. sericea} .

PLINK is slightly less sensitive to parent-offspring relationships, notably less sensitive to sibship, and generally shows a more realistic distribution of kinship values for siblings in most species, with many relationships designated as second degree. Additionally, PLINK exhibits a smaller zero-degree peak in \emph{B. aemula} and \emph{B. serrata} parent-offspring data, which might be apomictic relationships not as distinct in other methods. RelateAdmix yields kinship distributions similar to those of PLINK for both parent-offspring and sibling relationships. The putative apomictic relative cluster is also apparent.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/kin-hist-parent-offspring/kin-hist-parent-offspring}
\caption{{\label{fig:kin_hist_po} \internallinenumbers
    Histograms depicting kinship results for known parent-offspring relationships. Parent-offspring sensitivity (PO. Sens), parent-offspring precision (PO. PPV), and overall precision (PPV) for each specific method and species run are given in the top left of each panel. Background shading corresponds to the theoretical kinship levels for zero-degree relatives (dark), first-degree relatives (medium), and second-degree relatives (light). Zero- and first-degree results were considered true positives for parent offspring relationships in the sensitivity and precision calculations. The results showcased herein are derived from a specific set of filtering parameters per method (Figures S19-S24). Details about these filtering parameters can be found in the bottom panel of each column, aligned with the respective method.%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/kin-hist-siblings/kin-hist-siblings}
\caption{{\label{fig:kin_hist_sib} \internallinenumbers
    Histograms depicting kinship results for known sibling relationships. Sibling sensitivity (Sib. Sens), sibling precision (Sib. PPV), and overall precision (PPV) for each specific method and species run are given in the top left of each panel. Background shading corresponds to the theoretical kinship levels for zero-degree relatives (dark), first-degree relatives (medium), and second-degree relatives (light). Zero, first, and second-degree results were considered true positive for sibling relationships the sensitivity and precision calculations. The results showcased herein are derived from a specific set of filtering parameters per method (Figures S19-S24). Details about these filtering parameters can be found in the bottom panel of each column, aligned with the respective method.%
}}
\end{center}
\end{figure}

\subsection*{Impact of filtering}\label{impact-of-filtering}

The percent variance explained by the random forest (RF) models indicates how well the model explains the data. We found that the sensitivity of the Goudet, KING, and PLINK methods was highly explainable by the data (\textgreater78\% in most cases). PC-Relate and RelateAdmix sensitivity was less explicable by the models, with 9--79\% of the variance explained. Overall, precision was typically less explicable than sensitivity by the RF models (Table S3).

Random forest models produce an increase in node purity (IncNodePurity) and increase in mean squared error (IncMSE) for each predictor variable, which are indicators of the importance to the predictive power of the model (Figure S26). For the sensitivity RF models, minor allele frequency was most important for the Goudet method, locus missingness for KING Homo and KING Robust methods, and locus missingness and minor allele frequency for the PLINK method. There were no clear trends in the RelateAdmix RF models, which had lower explanatory power for the outcomes. Similarly, trends were less obvious for the precision RF models, which explained a lower percentage of the variance.

The partial dependence plots of the RF models show the effect of the predictor variables upon the response variable. In the case of the sensitivity RF models (Figure \ref{fig:pdp_sens}), Goudet had very high sensitivity regardless of predictor variables or species, consistent with the results in Figure \ref{fig:ppv_sens}. This was similar for PC-Relate and RelateAdmix. The sensitivity of KING homo and KING robust did not have strong trends dependent upon minor allele frequency or sample missingness, but showed a decay then plateau with increasing maximum locus missingness. The sensitivity of both KING methods is highly dependent on the dataset (Figure \ref{fig:pdp_sens}). PLINK showed an increase and then plateau in sensitivity with increasing minor allele frequency for most species. PLINK also showed a linear decrease in sensitivity with increasing maximum locus missingness (Figure \ref{fig:pdp_sens}).

The partial dependence plots of the precision RF models showed that for the Goudet method, increase in minimum minor allele frequency or maximum locus missingness led to greater precision in the analyses of the \emph{Banksia} species. Precision also increased then plateaued with increase in locus missingness in PC-relate and PLINK. With PLINK, increase in minimum minor allele frequency decreased the precision somewhat. There were no clear trends in RelateAdmix (Figure \ref{fig:pdp_ppv}), consistent with the low variance explained by these models (Table S3).\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/pdp-sens/pdp-sens}
\caption{{\label{fig:pdp_sens} \internallinenumbers
    Partial dependence plots generated by random forest models illustrating the predicted change in overall sensitivity due to filter variables. Each line represents the prediction from the random forest model, colour-coded by species. The top row shows the predicted sensitivity changes with increasing minimum minor allele frequency. The middle row shows the predicted changes with an increase in maximum sample missingness. The bottom row shows the predicted changes with an increase in maximum locus missingness.%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/pdp-ppv/pdp-ppv}
\caption{{\label{fig:pdp_ppv} \internallinenumbers
    Partial dependence plots generated by random forest models illustrating the predicted change in overall precision due to filter variables. Each line represents the prediction from the random forest model, colour-coded by species. The top row shows the predicted precision changes with increasing minimum minor allele frequency. The middle row shows the predicted changes with an increase in maximum sample missingness. The bottom row shows the predicted changes with an increase in maximum locus missingness.%
}}
\end{center}
\end{figure}

\subsection*{Performance in admixed and minority genetic groups}\label{performance-in-admixed-and-minority-genetic-groups}

We looked specifically at admixed, minority, and majority genetic groups within the structured species \emph{Acacia suaveolens} and \emph{Hakea sericea} to qualitatively compare how each method performs when the assumptions are potentially violated. Although we do not know the exact relationships among untracked individuals at the same site (i.e., not known mothers or seedlings), we generally assume that these relationships are less close than known parent-offspring or sibling connections overall. For known relationships, we anticipated that more seedlings in parent-offspring pairs result from sexual reproduction (first degree) rather than apomixis (zero degree). Similarly, in sibling relationships, we expect more half-siblings (second degree) than zero or first degree.

Goudet's method produced significantly higher kinship values between known relatives and unknown individuals in small structure groups of both \emph{A. suaveolens} and \emph{H. sericea} compared with larger or admixed groups (predominantly zero degrees; Figure S27). This pattern was similarly observed in PC-relate and PLINK, where most pairwise connections in small structure groups reached the maximum values possible for the methods. Both KING methods appeared to be resilient to the impact of small structure groups, showing minimal discernible differences between results from small and large structure groups. RelateAdmix was resilient to the effect of small structure groups, yielding more comparable distributions of kinship values between small and large groups (Figure S27). However, it still tended to overestimate kinship in "unknown" relationships, consistent with the previously identified issue of low precision.

\subsection*{\texorpdfstring{Case study: \emph{Dodonaea viscosa}}{Case study: Dodonaea viscosa}}\label{case-study-dodonaea-viscosa}

The PCA and phylogenetic network of \emph{Dodonaea viscosa} indicated that there was substructure within the dataset (Figure \ref{fig:dodovisc_pop}). We manually created subclusters in the dataset based on these results. These subgroups aligned with a distinction observed in the $F_{ST}$ analysis, where between sites within subgroups $F_{ST}$ was typically \textless{} 0.4, and between sites in different subgroups $F_{ST}$ was \textgreater{} 0.4 irrespective of distance (Figure \ref{fig:dodovisc_pop}H). Notably the subgroups were not geographically isolated (Figure \ref{fig:dodovisc_pop}D). The number of ancestral populations determined using STRUCTURE and the Evanno method was 2, fewer than but broadly aligned with the manual subgroups identified (Figure \ref{fig:dodovisc_pop}G).

We compared the outcomes of PLINK and KING Robust under our recommended filtering parameters. As expected, PLINK run on the whole population appeared to overestimate kinship in the minority subgroups compared with the other methods (Figure \ref{fig:dodovisc_kin}A--C). PLINK run separately on the subgroups had a relatively strong concordance with KING Robust, which gave us confidence that it was not producing an excess of false positives.

The PLINK by subgroup analysis revealed a significant number of relatives, especially within the SPAs (Figure \ref{fig:dodovisc_kin}D). The proportion of individuals with at least one first-, second-, or third-degree relative was higher in SPAs than among wild-collected individuals (1--45\% higher; \ref{fig:dodovisc_kin}E). Additionally, in all SPAs, individuals exhibited a greater number of detectable relationships across all degrees of relationship, excluding zero degrees (first to fourth degree), when compared with their wild counterparts (Table S5).\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/DodoVisc-combined2/DodoVisc-combined2}
\caption{{\label{fig:dodovisc_pop} \internallinenumbers
    Population structure analyses of Dodonaea viscosa. (A) site collection map, (B) PCA of samples, and (C) phylogenetic network of samples coloured by seed production area (SPA). (D--F) The same plots with points coloured by subgroups. (G) STRUCTURE-estimated ancestry of samples for the optimal K=2. (H) Pairwise $F_{ST}$ and geographic distance between sites. Points are coloured blue when the pairwise comparison is within a subgroup, and red when the pairwise comparison is between subgroups.%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/DodoVisc-kin-and-relatives/DodoVisc-kin-and-relatives}
\caption{{\label{fig:dodovisc_kin} \internallinenumbers
    Outcomes of kinship analyses of the Dodonaea viscosa samples. (A) Pairwise kinship values estimated using PLINK on the whole dataset versus KING Robust. (B) Pairwise kinship values estimated using PLINK separately on each subgroup versus KING Robust. (C) Pairwise kinship values estimated PLINK on the whole dataset versus PLINK separately on each subgroup. (D) Network of relationships detected by PLINK by subgroup. Nodes are individuals, coloured by Seed Production Area (SPA). Edges are detected relationships coloured by degree. (E) The percentage of individuals in each SPA with at least one detected first-, second-, third-, or fourth-degree relative, as well as the percentage of individuals with no detected relatives (unrelated).%
}}
\end{center}
\end{figure}

\section*{Discussion}\label{discussion}

Our evaluation of six methods of kinship estimation has revealed species-specific challenges, with population structuring significantly influencing method performance. While some methods exhibited robustness in certain species (PLINK, KING Homo, KING Robust), others encountered precision issues (RelateAdmix, PC-relate, Goudet's beta dosage). Based on the results of our study, we offer guidance on the appropriate utilization of specific kinship methods with reduced-representation sequencing data to effectively address research objectives.

\subsection*{Challenges posed by species-specific characteristics}\label{challenges-posed-by-species-specific-characteristics}

The species-specific challenges observed in our study highlight the importance of understanding population dynamics and reproductive strategies when selecting kinship estimation methods. The differentiation in population structuring among species, such as the highly homogeneous \emph{Banksia} species and the more complex \emph{Hakea sericea} , influenced the performance of each method. \emph{Banksia aemula} and \emph{Banksia serrata} , with minimal population structuring, demonstrated the robustness of several methods, particularly Goudet's beta dosage and KING methods. However, the challenges posed by \emph{Hakea teretifolia} , \emph{Acacia terminalis} , \emph{Acacia suaveolens} , and \emph{Hakea sericea} , each with a unique population structure, underscore the need for methodological adaptability.

\subsection*{Performance of kinship estimation methods}\label{performance-of-kinship-estimation-methods}

We tested these kinship methods on a variety of datasets, from simple to more biologically complex, using a wide range of data filters. Our findings indicate that each method behaves differently on diverse datasets, and that their performance is affected by the chosen filtering parameters, albeit to varying extents.

Goudet's beta dosage method is intuitive, in that it provides a kinship value that is relative to the average kinship of all pairwise relationships in the dataset. It has already been shown that the accuracy of this method is greatly reduced in datasets with fewer SNPs \parencite{goudet_how_2018}, so we anticipated that it might not be the most appropriate for our reduced-representation sequencing data.

In practice, Goudet's method is very sensitive, but it appeared to overestimate kinship among known relatives (Figures \ref{fig:kin_hist_po} and \ref{fig:kin_hist_sib}) and its estimates are extremely imprecise under most data filters in most species (Figure \ref{fig:ppv_sens}). It performed best in the analyses of the \emph{Banksia} species, which are homogeneous and had a high proportion of relatives in the dataset. Precision was very poor in the cases where there was any degree of population structure or restricted gene flow, or where a large proportion of the dataset was unrelated individuals. There were also no clear trends in precision and sensitivity across all datasets in relation to filtering parameters (Figures \ref{fig:pdp_sens}, \ref{fig:pdp_ppv}, S19-S24), which makes it difficult to have confidence in this method when applying to datasets in which the relationships among individuals are unknown. This method also overinflates kinship in minority structure groups (Figure S27), which is unfortunate given that it is supposed to be able to cope with population heterogeneity \parencite{goudet_how_2018}.

PC-Relate and RelateAdmix produced similar results to Goudet. Both were highly sensitive to all known relationships across all datasets and all filter values, but their precision was very dependent upon the dataset. In some species, these methods produced consistently precise kinship estimates (\emph{B. aemula} and \emph{B. sericea} ). In other species, the results were very variable and did not have a consistent association with the filtering parameters (\emph{H. teretifolia} and \emph{A. terminalis} ). In the worst cases, where populations were highly structured, precision was always poor (\emph{A. suaveolens} and \emph{H. sericea} ; Figure \ref{fig:ppv_sens}). This is regrettable given that these methods are designed to address population structure but seem to be less adept at handling it than older methods like PLINK and KING. Neither PC-Relate nor RelateAdmix is consistently more sensitive and precise across all species. However, RelateAdmix might be more useable because it uses a 0--0.5 scale, making the outcomes more biologically meaningful, and it does not seem to overinflate known relative kinship as much as PC-Relate (Figures \ref{fig:kin_hist_po} and \ref{fig:kin_hist_sib}). RelateAdmix also provides more consistent kinship estimation between small and large structure groups (i.e., kinship is not artificially inflated in minority subpopulations; Figure S27).

Although Goudet, RelateAdmix, and PC-Relate methods performed poorly here, we note that they were being tested on much smaller datasets than they were originally intended for and have previously been tested upon. The original evaluations of these methods used datasets of 20,000 to 600,000 SNPs \parencite{moltke_relateadmix_2014, conomos_model-free_2016, ramstetter_benchmarking_2017, goudet_how_2018}. In our investigation, the application of optimal filters for sensitivity and precision in each method yielded datasets with fewer than 40,000 SNPs for PC-Relate, 25,000 for RelateAdmix, and 5,000 for Goudet (Figures S19-S24, Table S2). Given this context, it is unsurprising that these algorithms, typically known for their accurate performance, encounter challenges when applied to our reduced-representation sequencing data and as a result should not be considered as useful for smaller datasets.

KING Homo and KING Robust were both very precise across different filtering parameters, with a negligible false positive rate. They had limited sensitivity, however, detecting less than 60\% of known relationships under most scenarios. Sensitivity was highly species dependent; KING was much more sensitive in our analyses of \emph{H. sericea} , which is highly structured and has a higher degree of inbreeding than the more homogeneous species (Figure \ref{fig:ppv_sens}). Additionally, KING was more sensitive when filtered with lower maximum locus missingness (Figure \ref{fig:pdp_sens}), and KING Robust was more sensitive than KING Homo overall, indicating that the within-family correction was improving the detection of relatives. Both KING methods are also resilient to small structure groups and do not overinflate kinship in minority subpopulations (Figure S27) so would be particularly suitable for species with complex population structure or extensive inbreeding.

PLINK had the most variable outcomes of any method depending on the filtering parameters, though the association was also more predictable compared with the other kinship methods (Figure \ref{fig:ppv_sens}, Table S3). PLINK is typically less sensitive and more precise at higher maximum locus missingness and lower minimum minor allele frequency and becomes more sensitive and less precise as maximum locus missingness decreases and minimum minor allele frequency increases (Figures \ref{fig:pdp_ppv}, \ref{fig:pdp_sens}, S19-S24). This is interesting because the widely accepted view of minor allele frequency filtering in kinship estimation is that lower is better \parencite{manichaikul_robust_2010, goudet_how_2018}, which is not the case here. The sensitivity and precision of PLINK is somewhat species dependent, but not as much as KING (where sensitivity is highly species dependent; Figure \ref{fig:pdp_sens}) or Goudet, RelateAdmix, and PC-Relate (where precision is highly species dependent; Figure \ref{fig:pdp_ppv}). The trade-off between sensitivity and specificity of PLINK means that in most cases it is not 100\% sensitive, but it can remain precise even at higher sensitivities across different species' data unlike the other methods.

PLINK does have some important limitations, however. Many pairwise relationships, especially sibship, had estimated kinship of 0 (the minimum value for this method), which is not seen in the other methods (Figure \ref{fig:kin_hist_po}). PLINK also overestimated kinship in small structure groups (Figure S27), which is expected given the use of population allele frequencies and assumption of Hardy-Weinberg equilibrium in kinship calculation, with no correction for population structure. Both PLINK and RelateAdmix offer the additional advantage of distinguishing between first- and zero-degree relatives in the \emph{Banksia} datasets when assessing parent-offspring kinship (Figure \ref{fig:kin_hist_po}). This capability proves valuable for identifying instances of apomixis or asexual clones in plants.

\subsection*{\texorpdfstring{Choosing an appropriate method }{Choosing an appropriate method }}\label{choosing-an-appropriate-method}

We aimed to conduct a comprehensive and unbiased comparison of kinship methods, but practical constraints within our datasets introduce limitations. Our intention was to assess the performance of these methods in real-world scenarios with non-model species. However, working with such data required us to make assumptions about unrelated individuals and necessitated the exclusion of a substantial amount of data (e.g., kinship information within 10 km where relationships were unknown). This exclusion means that our precision estimates are likely to be inflated. Another limitation arises from the prevalence of closely related individuals in our datasets, a situation that is uncommon in most wild species. Although we explored the impact of this factor on our results, we observed no consistent change in method performance (Table S3). We uphold the validity of our method comparison, as it offers a unique evaluation of the performance of kinship estimation methods using accessible genomic data (reduced-representation sequencing) from wild species that experience restricted gene flow, population structure, and variable degrees of inbreeding, which could not be easily reproduced by simulations.

Overall, there is no single kinship estimation method or set of filtering parameters that achieves perfect detection across all the different scenarios. The outcome is partially dependent upon the dataset itself, including quality and biological underpinnings which are an inextricable component of working with multiple non-model species, and partially dependent on the filtering. Based on the results of our evaluations, we provide some recommendations for kinship estimation.

First, PLINK is the best general option due to its adept handling of both homogeneous and complex datasets, striking a balance between precision and sensitivity. However, users should be mindful of its sensitivity to filtering, and the method might overestimate kinship in minority structure groups. If a distinct minority structure group is identified, PLINK should be run separately on this subgroup. We suggest employing a minimum minor allele frequency of 10\% and a maximum locus missingness of 20--30\% to optimize the trade-off between sensitivity and precision.

Second, for cases where precision is more important than sensitivity, KING Robust is very consistently precise across all datasets. It is a particularly good choice for species that have a high degree of structure and/or inbreeding. If KING is used, we recommend using the lowest maximum locus missingness threshold that would be reasonable for the data of interest (\textless{} 10\% ideally, while retaining thousands of SNPs) and low or no minimum minor allele frequency.

Third, in most cases, we do not recommend using Goudet, PC-Relate, or RelateAdmix for estimating pairwise relationships in non-model species using reduced-representation sequencing data. These methods struggled to maintain overall precision in our analyses, particularly in scenarios involving population structure or restricted gene flow. The poor performance of these methods is likely to be due primarily to the limited number of SNPs in our reduced-representation sequencing datasets but are probably more suitable for high-density SNP datasets.

Finally, we suggest using multiple methods, such as applying PLINK to the entire population, as well as to separate subpopulations, alongside KING Homo and KING Robust. It is prudent to use different methods and assess the consistency of results. Different methods operate on distinct assumptions and are informative under varying conditions, making it valuable to compare their outcomes.

\subsection*{\texorpdfstring{Case study: \emph{Dodonaea viscosa}}{Case study: Dodonaea viscosa}}\label{case-study-dodonaea-viscosa-1}

Our investigation of \emph{Dodonaea viscosa} provides an applied example of the outcomes of our methodological evaluation. Through PCA and phylogenetic network analysis, we identified strong population structuring, with minimal gene flow between subgroups. The presence of structure meant that we needed to consider kinship overestimation in the minority structure groups.

Our PLINK analysis by subgroups identified numerous relatives within both the SPAs and wild populations of \emph{D. viscosa} . In SPA-B, two extensive families were discovered, comprising solely first- or second-degree relatives. Similarly, a sizable family was found in SPA-A, along with smaller clusters or connections among individuals across different SPAs (Figure \ref{fig:dodovisc_kin}D). These findings suggest that plants in the SPAs likely originate from seeds collected from the same plant rather than being sourced from diverse locations. These establishment practices have led to a higher proportion of relatives in the SPAs compared with the wild population (Figure \ref{fig:dodovisc_kin}E; Table S5). High relatedness in the SPAs increases the risk of inbreeding and loss of diversity, lead to lower offspring fitness, and limited adaptive potential, which contravene the objectives of restoration projects. Being able to accurately assess kinship levels enables us to provide recommendations to SPA managers on how to reconsider and update the current planting scheme to enhance the SPA's overall diversity and yields, as well as improve restoration targets and outcomes.

\subsection*{Conclusion}\label{conclusion}

The results of this study bring to mind the aphorism "all models are wrong, but some are useful." While all the kinship estimation methods investigated here have incorporated various measures to address the shortcomings of their predecessors, each comes with notable advantages and disadvantages. Describing the complexities of non-model species is challenging, and there is no single method with perfect accuracy in all situations. Nevertheless, these methods remain valuable tools for gaining insights into the complex dynamics of non-model biological systems.

\section*{Acknowledgements}\label{acknowledgements}

We acknowledge the Traditional Owners and Custodians of the lands on which this research was conducted. We pay our respects to their Elders, past, present, and emerging, and recognize their ongoing connection to Country, culture, and knowledge systems. We acknowledge that the knowledge and wisdom of Indigenous peoples have sustained and enriched these lands for millennia, and we are committed to fostering respectful relationships and engaging in meaningful collaboration with Indigenous communities in our research.

We extend our thanks to all staff and volunteers at the Botanic Gardens of Sydney who contributed to the project. Special thanks to Tom Pyne, Susan Rutherford, Alexander Carey, Katrin Quiring, Allison Bernardi, Jamily Machado da Silva, Alasdair Scott, Tom Lucchese, Joel Cohen, Lotte von-Richter, and Veronica Viler for their assistance in collecting and growing seedlings. Thanks to the managers of the seed production areas for allowing us to include their case study in this publication. Thanks to Jason Bragg for administrating and maintaining the computing facilities used to conduct our analyses.

\section*{References}\label{references}
\sloppy

\printbibliography[heading=none]

\section*{Data accessibility }\label{data-accessibility}

DArT SNP data will be made available on Dryad after an embargo period as we are still working on these data.

\section*{Funding}\label{funding}

This study was partially supported by the NSW Environmental Trust (grant 2016RD 0084) and the Restore and Renew project \parencite{rossetto_restore_2019}.

\section*{Author contributions }\label{author-contributions}

M.M. and M.R. acquired funding for collection, sequencing, and computing. P.L-I. and M.M. designed and executed the collection, germination, and sequencing preparation of the seedling datasets. M.M. and M.R. designed executed the collection of the \emph{D. viscosa} dataset. E.S.M. designed, analysed, and generated figures for the kinship comparison analyses and \emph{D. viscosa} case study. E.S.M. wrote the original manuscript, and S.Y.W.H., M.R., P.L-I., and M.M. worked on subsequent drafts.

\selectlanguage{english}
\FloatBarrier
\end{document}

