How Physical and Molecular Anthropology Interplay in the Creation of Biological Profiles of Unidentified Migrants

The skeletal sex and ancestry of unidentified human crania can be inferred both from physical and from molecular features. This paper depicts and discusses the experiences of physical and molecular anthropologists on a set of commingled crania from the largest Mediterranean shipwreck disaster on 18 April 2015, in order to facilitate identification of human crania. Twenty-one disarticulated crania that were recovered from the above-mentioned shipwreck were analyzed to estimate skeletal sex and ancestry, following a physical and a molecular pipeline. The physical analyses applied morphological and metric methods that provided posterior probabilities for the crania to be classified into a sex or ancestral group. The molecular analyses were performed on petrous bones via a shotgun sequencing approach that allowed us to determine the sex of each individual and to retrieve the complete mitochondrial genome, Y chromosome single nucleotide polymorphisms, up to 597573 SNPs across the human genome from each individual. The morphometric sex analyses showed that most crania belonged to male individuals, although some estimations remained uncertain or undetermined. Inconsistent results were obtained for ancestry estimation as well, since morphological methods classified the crania mostly as European/White, in contrast to the most numerous African forms determined by craniometric analyses. This quite agreed with molecular analyses that identified only African males. Overall, undetermined and contrasting results were obtained between disciplines, preventing the creation of reliable and sound biological profiles that could provide guidance on the sex and ancestral group of the victims. Therefore, the times may not be mature for a merger of physical and molecular anthropology. However, future investigations of this research avenue would pave the way to the possible development of novel tools, methods, and wider reference databases that could address the limitations of both disciplines.


Introduction
Recently, humanitarian emergencies related to forced migrations have intensified, and thus, 82.4 million people are estimated to be forcibly displaced worldwide [1]. Unfortunately, migration journeys and attempts at crossing borders have often resulted in mass fatalities. In Europe, this phenomenon has consistently increased since 2013, especially in the Mediterranean Sea, with countless tragedies that, most times, have remained silent and unreported [2]. Since 2014, the Mediterranean route has accounted for at least 25,000 migrant deaths [3], and therefore, it is considered to be the deadliest path from African and Eastern countries towards Europe. Within this plight, the victims' rights to identity and the relatives' rights to know should always be granted [4]. Therefore, forensic investigators have been requested to give a name to a victim, and, as part of their procedure, to first determine their biological profile that will provide insightful information to be compared against missing persons' data [5]. To do so, forensic experts must begin, particularly when working on skeletal remains, with estimations of biological sex, ancestry, age at death, and stature.
Skeletal sex and ancestry can be inferred both from anthropological evidence and from molecular evidence. Traditionally, physical anthropologists have deemed the human cranium as a reliable indicator for sex and ancestry estimation of unknown remains. On the basis of morphology and metrics, the cranium presents dimorphic traits between sexes [6,7] and it has also been demonstrated that variations of morphological cranial features exist between populations, according to the ancestral groups [8,9]. To this end, qualitative and quantitative methods have been developed that include morphological observations and metric analyses. Starting from the work of Buikstra and Ubelaker [10], Walker [7] proposed a scoring system of cranial traits that included the nuchal crest, mastoid process, supraorbital margin, glabella, and mental eminence and developed sex determination discriminant functions to determine the biological sex. This was the first study to provide a statistical framework for a subjective evaluation of cranial morphological traits for sex estimation. However, the author suggested caution in using the functions on populations other than the reference. The issue of population-specific equations was addressed by Krüger et al. [11]. The authors reported low accuracies of Walker's method when the discriminant functions were applied to a South African sample and provided modified equations. Similarly, Cappella et al. [12] developed specific regression models for an Italian contemporary cemetery population. The above-mentioned issues highlight that the morphological sex estimation method has some weak points when performed on unknown samples. Craniometric analysis represents a reliable mean for biological sex and ancestry estimation, although the lack of reference samples from which the statistical framework is inferred could complicate interpretation of the results [13]. However, in investigations such as the identification of unknown migrants, craniometric data have proven to be valuable indicators for classification into ancestral groups [5].
Since ancestry can be considered to be the genetic inheritance everyone receives from their ancestors and from population members that have occupied the same place of origin for long periods of time, molecular anthropologists could also provide valuable support to infer the biogeographical origin via DNA analysis. Recent advances in molecular technology associated with the possibility to also simultaneously genotype hundreds of markers from degraded bone samples have resulted in accurate and efficient estimation of an individual ancestry that could implement anthropological studies with additional evidence. DNA profiling has played an essential role in the positive identification of victims of several humanitarian cases, including mass burials in Guatemala [14,15] and Bosnia-Herzegovina [15], and mass disasters such as the World Trade Center, the Indonesian Tsunami [16,17], and the massacre at Fosse Ardeatine in Italy [18] (where reference DNA profiles from family members are available for comparison with those of the victims) but this is largely different from constructing a biological profile when there is no hypothesis of identity. From this perspective, DNA analyses on bones could provide further information to predict, for example, the biogeographical ancestry of a person of interest.
The fusion of physical and genetic anthropological methods, although apparently useful, may create further complications in the interpretation of a biological profile. In this article, we aim to share the experience of the application of both methods on a set of commingled crania belonging to the largest Mediterranean shipwreck disaster on 18 April 2015, in order to facilitate identification of human crania. In the context of the examination of the remains of circa 1000 victims, the present study focused on 21 disarticulated crania of these victims that were analyzed from a physical and genetic perspective in order to infer biological sex and ancestry. The final purpose, along with aiding identification, is to present the challenges that the practitioners faced when working on such a peculiar set of unknown remains and to highlight how forensic anthropologists (physical and molecular) can try to interplay in a conjoined investigation to provide more robust information and how this combination may still need further research for proper integration.

Materials and Methods
The sample included 21 disarticulated crania that were recovered from the shipwreck of 18 April 2015. The study complies with Police Mortuary Rules (DPR 09. 10.1990 No. 285, art. 43) and the Regio Decreto (08. 31.193331. No. 1592 and is encompassed within the identification project of the victims enshrined in an MOU between the University of Milan and the Italian government. The crania were brushed with water and soap to remove dirt and soft tissue on the external surface and within the skull. Morphological and metric analyses were performed prior to sampling for DNA typing. Dental development was investigated through radiographic images performed with a portable dental X-ray unit (intraoral sensor (EZ Sensor 1.5, Vatech, Hwaseong, Republic of Korea) and X-ray tube (model Rextar, Poskom Co. Ltd., Goyang, Republic of Korea)). There were 14 crania with completed dental development (the apices of the third molar completely closed). Seven crania presented evidence of on-going root growth of third molars which were evaluated according to Mincer et al. [19]. Since the reference populations of the method are hardly representative of the ancestry of the sample, only the stage of development is provided (Table 1). For all crania, the petrous bone was chosen as the preferred anatomical element for DNA analysis [20][21][22] and was sampled with a Stryker hacksaw.
Metric data were acquired using a microscribe digitizer, and the 3Skull software [25] that obtains the coordinates x, y, and z of each landmark and calculates the inter-landmark distances of the measurements defined by Howells [26]. The crania were classified against West and East African, Somalian, Zulu, Hainan Chinese, Portuguese, and Euro-American population samples. The West African group of population samples are described in detail in Spradley and Jantz [27]. The Somalian group of population samples were likely battlefield casualties of the Italian invasion of British Somiland in 1940. The Zulu and Hainan group of population samples are from Howells [26]. The African group of population samples all date to the 19th or early 20th century. As such, the population samples are not entirely suitable as groups to compare with the modern crania represented by the shipwreck victims that may have experienced secular change. The analysis of the 21 crania was a two-step process. The first step was to estimate the population affinity of each cranium based on its morphometric similarity to a reference sample.
Step two was a search for a substructure that may exist among the crania, based on their distances from one another.
Step one used the Mahalanobis generalized distance of each of the 21 crania from the reference samples. The principal component scores for the reference samples were derived from within the covariance matrix. Then, the eigenvalues and eigenvectors were used to obtain principal component scores for each of the 21 crania used in the present study. They were used to compute the Mahalanobis distance of each cranium from the reference samples. A cranium's population affinity was assigned to the group to which it had the lowest Mahalanobis distance.
Step two used the Mahalanobis distance of each cranium from every other cranium, yielding a 21 × 21 distance matrix. The principal coordinates from this matrix were extracted using the procedure described by Gower [28]. The principal coordinates were calculated in such a way as to maximize variation among crania variances. It is important to stress that the principal coordinates are chosen without regard to the population affinity described above. The only relationship to the reference sample is that the reference sample within the covariance matrix is used to obtain the principal component scores. This only requires that the covariance matrices from which the 21 crania were drawn do not differ significantly from those of the reference samples, which is not an unreasonable assumption. After obtaining the principal coordinates, post hoc tests of whether they sort populations according to affinities identified in Step one were performed. The population affinity of each cranium was attached to its principal coordinate score. The scores could then be displayed graphically and analyzed statistically. Table 1. Summary of the sex and ancestry estimation of the 21 crania. M/F, a cranium with a higher number of probabilities that classify it as male; F/M, more probabilities classify the cranium as female. In brackets is the highest posterior probability for the cranium to be assigned to an ancestral or sex group. Skeletal sex estimation was based on the five morphological traits suggested by Walker [7]: nuchal crest, mastoid process, supraorbital margin, glabella, and mental eminence. The five traits were scored for each cranium, except for the mental eminence that could not be evaluated in all cases, since 10 crania did not present the mandible. The skeletal sex was estimated by calculating the equations that were related to the populations indicated by the ancestry estimation. When a cranium was classified as White, Krüger's equations [11] for White South African were applied, and vice versa when a cranium was classified as Black. In addition, the original equations by Walker [7] were applied, although the method does not provide specific equations for Black and White individuals. A cranium was considered to be male (M) or female (F) when all the traits considered by the equations agreed. When a cranium presented mixed traits, it was not possible to express a definitive judgment. This meant that, in some instances, the highest probability when considering one combination of traits (e.g., glabella and mastoid process) classified the cranium as female and as male when considering another couple of traits (e.g., nuchal crest and mastoid process). As such, the sex estimation is indicated as M/F, if the cranium has a higher number of probabilities that classify it as male, or as F/M if more probabilities classify the cranium as female. The metric data acquired for ancestry estimation were also used for craniometric sex estimation. The Somalian crania were from males only, but all other samples included both sexes. Sexes were assigned from cranial morphology, except for Zulu, which were samples of a known sex. Zulu crania were used to develop sex discriminant functions using three highly dimorphic dimensions, glabellar projection (GLS), bizygomatic breadth (ZYB), and mastoid height (MDH).

Sample Preparation, DNA Extraction, Whole-Genome Library Preparation, and Sequencing
Bone samples were processed in the molecular anthropology unit of the University of Florence, a state-of-the-art facility dedicated to the analysis of degraded DNA samples. To remove potential contaminants, the outer layer of petrous bones was brushed using a dentist drill with disposable tips as suggested by Pinhasi et al. [21] and irradiated with UV light (λ = 254 nm) for 45 min in a Biolink DNA Crosslinker (Biometra, Goettingen, Germany). A minimally invasive approach was followed to recover approximately 50 mg of bone powder from petrous bones, as described by Sirak et al. [29]. DNA was extracted using silica-based protocol [30] and eluted in TET buffer (10 nM Tris, 1 mM EDTA, 0.05% Tween-20) twice for a final volume of 100 µL Subsequently, a double-stranded and dual-indexed Illumina DNA library was prepared from 25 µL of each extract following previously published protocols [31,32]. For all samples, deaminated cytosines resulting from DNA damage were partially removed using an uracil-DNA-glycosylase treatment (UDG-half) as proposed by Rohland et al. [33]. Then, a unique combination of two indices per library was used for barcoding, the libraries were pooled in an equimolar amount and sequenced in 200 bp paired-end mode on a NovaSeq 6000 instrument (Illumina, San Diego, CA, USA) for a depth of~14 million reads for samples, except for 58-A1, 62, 131, and 149-1 for which deep sequencing was performed. Negative controls were used in all the experimental steps to monitor the absence of contaminants in the reagent and the environment.

Postsequencing Data Processing
Raw data were demultiplexed, allowing for a maximum of one mismatch in each index, and processed through the EAGER pipeline (v1.92.59) [34] keeping paired-end reads separate. Adapters were trimmed using Clip&Merge v1.7.4 with the parameter "perform only adapter clipping" [34] and reads shorter than 30 bp were discarded. Filtered reads were subsequently mapped against the human reference genome (Hg19) using BWA v. 0.7.10 [35] with the stringency parameter set to 0.01 and a mapping quality filter of 30. The duplicate removal was performed using DeDup (v0.12.2) [34], which considered both start and end coordinates of the reads to identify and discard identical sequences.
Before genotype calling, reads mapped onto the human genome were authenticated by deamination and fragmentation pattern analysis using mapDamage2.0 [36]. The presence of degraded DNA typical features was checked as a short average length of DNA sequences and an increased proportion of miscoding lesions at the terminal position of the molecules. Simultaneously, extraction and library preparation blanks were also analyzed for the presence of potential background contamination during the laboratory activity.

Sex Determination
The molecular sex was estimated from each individual by comparing the number of alignments to the Y chromosome and the total number of alignments to the X and Y chromosomes in the libraries prepared [37].

Genotype Calling and Principal Component Analysis for Ancestry Inference
Genotyping was carried out on the BAM files. A pseudo-haploid genotype was reconstructed for each individual with pileupCaller (Phred-scaled base quality score ≥30) by performing a calling of 597,573 alleles from the 1240 K + human origin dataset. Subsequently, genotyped data were merged with the human origin panel and the 1240 K dataset (https://reich.hms.harvard.edu accessed on 16 November 2022) selecting only 7566 modern individuals (African = 1195, European = 1237, Asian = 4563, American = 564, and Ocean = 7). The software smartpca from the EIGENSOFT package (v16000) [38] was used to compute a world PCA using all individuals of the human origin dataset.

Haplogroup Assignment of Uniparental Markers
To perform the mtDNA analysis, reads were processed with EAGER and alignment to the mitochondrial reference genome (rCRS) using CircularMapper, a tool that considers the circularity of the mtDNA [34]. Reads with mapping quality ≥30 were selected and extracted from the BAM files using SAMtools v1.7 [39], and the consensus sequence was called using mpileup and vcfutils.pl in the SAMtools package. The mitochondrial haplogroups were assigned according to the rCRS using MITOMASTER [40,41]. The D-loop variants (16024-576 bp) of each mitochondrial genome obtained from the previous analysis with MITOMASTER were uploaded into the EMPOP platform (v4/R13) (http://empop.online/ accessed on 19 December 2022) to calculate the frequency of the genetic profile within the reference database and to assess its worldwide distribution through heatmaps [42].
The Y chromosome haplogroup was also determined using the Yleaf 2.2 software [43], a Phython-based tool for effective analysis and interpretation of Y chromosome NGS data. This easy-to-use and publicly available software performs NRY single nucleotide polymorphism calling and subsequent haplogroup inference, based on ISOGG markers (https://isogg.org/tree/ accessed on 13 December 2022), starting from the BAM file. Reads with mapping quality ≥30 were selected.

Physical Data Skeletal Ancestry and Sex Estimation
Regarding morphoscopic ancestry, 14 crania were classified as White and 7 crania were classified as Black, according to the OSSA scores [23]. The decision support system HefneR [24] provided the following outcomes based on the highest probability: 10 crania were classified as European, 8 crania were classified as African, and 3 crania were classified as Asian. In detail, 12 crania were classified into a single population group with a probability equal to or above 0.8; in six cases, the highest probability was lower than 0.8 (between 0.72 and 0.50); in three crania, the highest probabilities to be classified in one group were equal to or lower than 0.50, with lower distributions in other classification groups. In five cases, inconsistent results between the two morphoscopic methods were obtained: in these instances, the same cranium was classified into two different population groups according to the OSSA score and HefneR. For example, the cranium Sample 095 was classified as White by the OSSA score, but its highest posterior probability (pp = 0.8) classified it as African by using HefneR.
An OSSA score identifies two major categories, i.e., Whites and Blacks, based on the score assigned to six of the traits suggested by Hefner [9]. Although this binary classification provides indications on the ancestral group of a crania from which one can start the ancestry analysis, it does not suggest further details. In contrast, the decision support system HefneR outputs posterior probabilities according to the affinity of all the traits by Hefner [9] for four classification groups (European, African, Native American, and Asian), thus expanding the two-way classification of the OSSA score and including two additional population groups. Although the two methods are both based on the same approach, that is, the semi-quantitative assessment of almost the same morphoscopic traits, some disagreements between the methods were observed. For example, in some instances, the OSSA score and HefneR classified the same cranium into two different population groups (i.e., African/Black and European/White). With respect to these crania, it is noted that the probabilities provided by HefneR were distributed over more than one group, with the highest value equal to or below 0.65, except the case of cranium Sample 095 whose probability to be classified as African was 0.80. Again, the results can only be partially compared between methods because of the different classification groups: three cases were classified as Asian according to HefneR with probabilities lower than 0.8; however, this population group is not considered by the OSSA score.
Based on the craniometric data, the population affinities are as follows: Somali 10, West African 7, Zulu 1, Portuguese 1, Euro-American 1, and Hainan 1. The posterior probabilities are generally low, only eight posterior probabilities are ≥0.7, five of which belong to Somali. Lower posterior probabilities are consequently distributed over several groups. Such a pattern indicates that low posterior probability crania do not fit well into any of the reference samples, which is to be expected given the likelihood that migrants come from populations other than those represented by the reference groups. Nevertheless, a strong African pattern emerges, 18 of 21 crania classify as either West African, Zulu, or Somali. The principal coordinate (PC) plot is shown in Figure 1. PC1 sorts those classified as Somali (PC1 scores <0) and from those classified as West Africa, including Zulu (PC1 scores >0). The ellipses, fitted by eye, include all or most of those with West African affinity vs. those with Somali affinity. The Somali ellipse includes two West Africans, the Euro-American, and the Portuguese. A t-test on the PC1 scores of Somali vs. West Africa and Zulu yields (t = 2.91, df = 16, p = 0.009). PC2 shows less dispersion and serves mainly to separate the cranium with Hainan affinity from the others.  The results are strongly suggestive of a substructure among these 21 crania. A course mesh can use the available samples to identify the substructure. Some groups that may have been on the ship are not represented, for example, Eritreans, and they predate the migrants by many decades. Future work may be able to more accurately identify a sub- The results are strongly suggestive of a substructure among these 21 crania. A course mesh can use the available samples to identify the substructure. Some groups that may have been on the ship are not represented, for example, Eritreans, and they predate the migrants by many decades. Future work may be able to more accurately identify a substructure.
For the morphological sex estimation, the results can be regrouped as follows: 1. Determined: When the morphoscopic methods agreed on the classification group (e.g., European and White), the sex classification followed the ancestry estimation, and therefore, the appropriate equations by Krüger et al. [11] for South African Black or South African White were applied. Eight crania belonged to this group: five crania displayed the highest probability to be classified as male for all the traits considered, whereas three crania were classified as female.

2.
Uncertain: When the morphoscopic ancestry estimation provided consistent results between methods, but the cranium presented mixed traits and therefore mixed probabilities according to the equations for sex estimation. For eight crania, a definitive judgment could not be expressed, as the coexistence of feminine and masculine traits produced mixed probabilities from the equations: four crania were F/M and four were M/F.

3.
Undetermined: When the morphoscopic ancestry estimation classified the cranium into different population groups, it was not possible to choose the appropriate recalibrated equation by Krüger et al. [11]. Therefore, the result was based on Walker [7], since the method pools White and Black individuals. Five crania belonged to this group; the sex estimation based on Walker [7] classified two crania as female and three crania as male.
Cranial traits for sex estimation present inter-population variation [8]; thus, the published methods may demonstrate poor accuracy or even be unreliable when applied to samples other than those from which they were developed [11,12,44]. In the present study, the sex estimation methods calculate discriminant functions using equations that are based on combined traits (e.g., nuchal crest and mental eminence or mastoid process and glabella); thus, the coexistence of masculine and feminine traits in the same cranium can lead to mixed results, where the cranium is classified as female according to some traits and as male according to other traits. Apart from the cases where the ancestry was indeterminate due to inconsistent results between morphoscopic methods, some crania presented mixed traits (i.e., more feminine traits combined with more masculine traits), and therefore, the discriminant functions could not produce a unanimous result. Here, in nine cases that presented such a mixture, a conclusive estimation could not be expressed; therefore, the cranium was classified as M/F or F/M according to the higher number of equations that output a higher probability for the cranium to be male or female, respectively.
Following the craniometric analysis, the sex discriminant classified 7 of the 21 crania as female, although two had approximately equal posterior probabilities (i.e., cranium 104-1 and cranium 125), and therefore no firm conclusion could be drawn. The remaining five crania had relatively high posterior probabilities for female. Based on sex established by DNA, these are clearly misclassifications. What this tells us is that cranial morphology of some migrants is more gracile than the Zulu from which the sex discriminants were derived. Table 1 summarizes the results of the physical analyses.
A human cranium may present mixed traits that prevent an exact classification into single groups, both for ancestry and sex. As observed in this study, some cases could not be classified as male or female because they presented mixed probabilities according to different traits, and the sex/ancestral group into which a cranium was classified differed according to the method. The size of the sample was limited to 21 crania that were already genotyped for identification purposes; inter-individual variability of such fluctuating features for sex and ancestry estimation within this pool of crania may be therefore clarified by increasing the sample. When applying statistical methods for ancestry and sex estimation, anthropologists must consider that the reference population data may not be appropriate and representative of the sample of the study [5]. Another factor influencing the anal-yses may be related to secular change [6,13] potentially experienced by this population, which would benefit from further investigation. As a result, indefinite ancestry and sex estimations represent a possible scenario where the strong dependence of the methods on the original reference sample is the main element that may hamper the creation of a biological profile.

Molecular Data
Sex Estimation and Skeletal Ancestry DNA was extracted from the petrous portion of the temporal bone of 21 individuals, double-stranded DNA libraries were created, and a shotgun approach was used to determine the sex of each individual and to retrieve the complete mitochondrial genome, Y chromosome single nucleotide polymorphisms (SNPs), up to 597,573 SNPs across the human genome from each individual. In Table S1 (Supplementary Materials) the shotgun results of each sample are summarized.
Subsequently, sex estimation was determined by considering the ratio of sequences aligning to the X and Y chromosomes and the molecular analysis highlighted that all crania belonged to individuals of male sex. Regarding the ancestry estimation, 597,573 SNPs were selected across the human genome. The SNPs called for ancestry estimation of each individual are summarized in Table 2. As can be observed in Table 2, the number of SNPs varies among different samples, ranging from 30,842 out of 597,573 SNPs for Sample 028 to 580,917 out of 597,573 SNPs for Sample 149. The principal component analysis plot (Figure 2) shows that our unknown samples form a genetic cluster that overlaps with the African individuals in the database. Therefore, the PCA analysis revealed relatedness/genetic similarity among all unknown individuals with the 1195 African individuals in the database. In an attempt to assess ancestry at the geographical area level, PCA was performed for the single African population after dividing it into three different areas (North Africa, East Africa, and West Africa). Individuals from Central and South Africa were not considered since the numbers of individuals in our datasets were low. It can be observed in Figure 3 that the African geographical areas form three different clusters and most of the unknown individuals (black dots in Figure 3) fall into the East Africa cluster, except for Samples 105-2, 154-1154-2, and 178-3 that fall into the West Africa cluster, as evaluated by comparison of the coordinates of database individuals and unknown samples. In an attempt to assess ancestry at the geographical area level, PCA was performed for the single African population after dividing it into three different areas (North Africa, East Africa, and West Africa). Individuals from Central and South Africa were not considered since the numbers of individuals in our datasets were low. It can be observed in Figure 3 that the African geographical areas form three different clusters and most of the unknown individuals (black dots in Figure 3) fall into the East Africa cluster, except for Samples 105-2, 154-1154-2, and 178-3 that fall into the West Africa cluster, as evaluated by comparison of the coordinates of database individuals and unknown samples.
In addition, to investigate the potential influence of sex biases in the biogeographical ancestry and to estimate when the common patrilineal and matrilineal ancestors lived, Y chromosome and mtDNA haplogroup were determined using, respectively, Yleaf and Mitomaster.
The Y chromosome results in terms of reads that mapped on the Y chromosome, number of SNPs called by Yleaf, QC scores, and haplogroup assigned by the software are reported in Table 3. An automatic prediction of the haplogroup with QC score greater than 0.81 was obtained for most samples. No prediction was shown for four samples (058 A1, 099, 100-1 and 104-5), probably because the QC score fell below 0.75. Therefore, manual interpretations of these samples' specific output files were required, as suggested by the Yleaf manual (https://github.com/genid/Yleaf accessed on 22 December 2022).
As can be observed in Table 3, most of the male individuals analyzed belonged to macrohaplogroups E (approximately 66.7% of individuals) and J (19%), except for Sample 003 and Samples 099 and 104-5 that showed, respectively, macrohaplogroups B (4.8%) and A (9.5%). As presented in different papers (for example [45][46][47][48]), haplogroups A, B, and E are overwhelmingly the most common in the African population. Whereas haplogroup J shows high frequency in the Arabian Peninsula, southern Mesopotamia, and the southern Levant as proposed by Sahakyan et al. [49]; however, the same study also highlighted a moderate frequency of this haplogroup in North and East Africa. In addition, Wood et al. [46] found that haplogroup J was concentrated in about 20% of the Afro-Asiatic studied. In addition, to investigate the potential influence of sex biases in the biogeographical ancestry and to estimate when the common patrilineal and matrilineal ancestors lived, Y chromosome and mtDNA haplogroup were determined using, respectively, Yleaf and Mitomaster.
The Y chromosome results in terms of reads that mapped on the Y chromosome, number of SNPs called by Yleaf, QC scores, and haplogroup assigned by the software are reported in Table 3. An automatic prediction of the haplogroup with QC score greater than 0.81 was obtained for most samples. No prediction was shown for four samples (058 A1, 099, 100-1 and 104-5), probably because the QC score fell below 0.75. Therefore, manual interpretations of these samples' specific output files were required, as suggested by the Yleaf manual (https://github.com/genid/Yleaf accessed on 22 Dicember 2022). Table 3. ID sample, number of reads mapped on the Y chromosome, number of SNPs used for the Yleaf analysis of each individual, and QC score, i.e., an overall quality score of the predicted haplogroup. The haplogroups that were assigned manually are shown in red.    Table 4. Regarding the mitochondrial results, the number of reads that map on the mitochondrial genome, the number of bases covered, and the haplogroup assigned by MITOMASTER are summarized in Table 5. As can be observed in Table 5, most of the individuals (90.47%) belong to haplogroup L, in particular to haplogroups L0 (5.26%), L2 (36.84%), L3 (52.63%), and L5 (5.26%), except for Samples 058 A1 and 104-5 that belong to, respectively, haplogroups T and U. The L haplogroups are African specific, indicating the African origin of our unknown individuals [50]. As proposed by Maier et al. [51] in their recent paper, the L0 subclade shows high frequency in South Africa and moderate frequency in Central and East Africa with a maximum distribution percentage (approximately 35%) of L0a (Sample 100-2) in Congo. Otherwise, the L2 subclade is distributed in most of Sub-Saharan Africa with a maximum distribution percentage (55%) of L2a (Samples 028, 131, 137, 146, 154-1, and 178-3) in Congo and with a distribution percentage of L2c (Sample 095) less than 19% in Mali [51]. The L3 subclade, instead, is distributed in most of Africa with low frequency in Botswana, South Africa and Algeria, Tunisia, and Libya. Figure S6 in the Supplementary Materials by Maier et al. [51] shows that the frequency of haplogroup L3b (Samples 062, 100-1, 178-1, and 178-2) has the greatest distribution percentage (approximately 18%) in Madagascar with high values also in Mali, Niger, and Chad. On the contrary, haplogroups L3e (Samples 003 and 105-2) and L3f (Sample 154-2) show, respectively, the greatest distribution percentage (48%) in the inner part of Western Sahara with high values also in Morocco and the greatest distribution percentage (11%) in Sudan. Haplogroups L3i (Samples 099 and 125) and L3x (Sample 149-1) show, respectively, the greatest distribution percentage (5%) in Somalia, Ethiopia, and Uganda and the maximum distribution percentage (7%) in Ethiopia. Finally, according to Maier et al. [51], haplogroup L5 (Sample 104-1) is distributed in most of Central Africa with the highest frequency in Congo (9.3%).

ID Sample
A different discussion must be made for Sample 058-A1 that belongs to haplogroup T1a. Haplogroup T makes up almost 10% of the mitochondrial genomes in Europe and approximately 8% in the Near East. Haplogroup T1 is distributed, albeit at varying frequencies, from Northwestern Africa throughout Europe, the Caucasus, and the Near East, into Western India, and across Central Asia into Siberia. The distribution of subclade T1a, despite being widespread, is nonhomogeneous and reaches values of approximately 4% and 8% also in Egypt and Algeria/Tunisia, respectively [52]. The high frequency of subclade T1a in Tunisia (38%) was also confirmed by Halim's study [53]. Haplogroup U2 is mainly distributed in Western Eurasia [54] and, although rare, haplogroup U2d has also been found in Ethiopia in both the Amhara region and Tigray [55], presumably as a result of the subsequent penetration from Western Asia and the Caucasus due to the absence of haplogroup U2d in North Africa.
As with the Y haplogroups, the worldwide distribution of mitochondrial haplogroups was assessed using EMPOP (Table 6).

Comparison of Physical and Molecular Data
In such a complex scenario, it is obvious that the interaction of physical and molecular anthropologists is of paramount importance to magnify efforts towards reliable and sound biological profiles that may provide guidance on the sex and ancestral groups of victims, in order to direct the search for relatives. However, the times may not yet be mature for such a merger.
This work presented a multi-method study that investigated and discussed the possible difficulties that anthropologists confront when studying human crania of unknown sex and origin, within the context of a migrant population.
The physical morphological analyses indicate that the sample is mainly composed of males, of White/European ancestry, although in several instances a conclusive sex and ancestry estimation was not possible due to opposing or inconsistent results between methods. Craniometrics confirm that the crania are mainly male, although the population affinity is more frequently African (Somali and West Africa). This is in agreement with the molecular analyses that show that the crania exclusively belong to male, African individuals. Based on the identity cards that were found associated with some of the victims of the 18 April 2015 shipwreck (but not necessarily associated with these crania), it seems likely that this ship carried people from Sub-Saharan countries. Clearly, no documents were found with the disarticulated crania; therefore, this is a completely unknown sample with respect to sex and ancestry. This hinders the evaluation of the accuracy of the estimations. However, this work did not aim to compare the accuracy rates of anthropological methods against DNA typing, but aimed to stress the elemental interdisciplinary teamwork of physical and molecular anthropologists and the difficulties that may emerge when dealing with biological profiles especially of migrant populations. Skeletal sex and ancestry estimation results based only on morphological traits of the cranium have been found to be 92% accurate compared to the DNA results [56]. Accuracy rates on skeletal ancestry (versus ancestry determined by documents) have been at 90% [57]. Here, it is not possible to assess the reliability of the observations because there is no positive or negative control sample, given the completely unknown origin of the sample, which represents both the major shortcoming and the purpose of this study. The contrasting physical results with respect to ancestry somehow also affect the sex estimation of these crania, since this is also usually inferred according to the ancestral group of the individual [58]. These crania were, therefore, considered to be "undetermined" both for ancestry and sex following Thomas et al. [57]. It should be pointed out that seven crania did not show complete skeletal and/or dental development. This may divert a conclusive sex and ancestry estimation by physical means given that related characteristics may be undeveloped [59], thus further reducing the reliability of the observations. Therefore, molecular analyses may prove to be essential to infer at least the genetic sex in individuals that present such an ambiguous set of traits and prevent a robust estimation [60].
Novel specific formulae that are more accurate for a given population can be suggested and tested for their validity, proving their reliability when applied to the appropriate sample [12]. However, the ancestry and sex prediction parameters, and especially the methods that strive for a quantification of the results, strongly depend upon the reference population on which they are developed, thus limiting the application to other populations [61]. This study raised the following question: How can anthropologists tackle the issue of sex and ancestry estimation when dealing with unknown specimens such as those of this study? With respect to sex, although this study lacks a positive control that could derive from the personal identification of these victims, the molecular analyses can be considered to be more reliable due to the stability of the genetic code and its determination. Regarding ancestry, this work shows that anthropology still finds itself mired in an endless conundrum: Which reference population is appropriate to estimate the ancestry of individuals? As morphology classifies crania into macro groups, craniometric analyses seem to provide more exhaustive information. Therefore, decision support systems and methods enhanced by statistical frameworks can provide valuable insights for physical analyses, but only by guiding an-thropologists through the interpretation of evidence and decision making rather than by identifying binding groups [8]. Regarding the molecular analysis, the 597,573 SNP callings in association with the PCA model were assessed to infer the biogeographical ancestry of unknown individuals. The PCA analysis (also known as exploratory analysis), which is one of the most used statistical clustering methods for the inference of the biogeographical ancestry of a person, allowed us to create four separate clusters corresponding to African, American, Asian, and European individuals and unequivocally assign all the 21 crania to the African population through a visual, intuitive, and easy-to-interpret approach. The autosomal data were also confirmed by the analysis of uniparental markers (mtDNA and Y chromosome) and the assignment of the African haplogroup to all 21 analyzed crania. In addition, athe PCA was also performed at an intra-continental level considering only the African population and dividing it into three different areas (North, West, and East Africa). The obtained results highlight that most of the unknown individuals fall into the East Africa cluster, except for Samples 105-2, 154-1154-2, and 178-3 that fall into West Africa. However, since PCA modeling is not a reliable discrimination/classification model [62] and the individuals from Central and South Africa were missing from the database used, the ancestry inference, especially at an intra-continental level, may not be accurate. The difficulty of ancestry inference at an intra-continental level was also confirmed by the results of uniparental markers. In summary, the molecular analysis clearly inferred the continent of origin of the 21 crania, but it was not possible to define a single geographical area of origin for each cranium. In these cases, particular attention should be paid to the database. Since the analysis of ancestry inference was performed by comparing the sample genotype with one or more known reference population groups, well-characterized databases with high-quality genotyping results of well-defined reference populations are critical.
Despite this, as shown in Figures 4 and 5 where the comparison between anthropological (physical and molecular) results relating to sex and ancestry are summarized, the molecular analysis support the physical data in both sex and ancestry inference, confirming the presence of 21 male crania of African origin. As proposed by Alladio et al. [62] and further developed by Pilli et al. [63], to infer intra-continental biogeographical ancestry in a forensic context, we suggest adopting a well-characterized and well-defined database and a robust classification method such as a multivariate statistical and machine learning approach associated with selected markers (AIMSNPs ancestry informative markers).
Again, our results cannot not be supported by antemortem data, but this was not the scope of the project. On a final note, physical and molecular analyses can produce contrasting results where an anthropologist is not able to confidently determine sex and ancestry, thus preventing the elaboration of the biological profiles and the identification activities. Although a strict reliance on reference groups still represents a drawback when analyzing completely unknown specimens, an approach that intertwines physical and molecular anthropology may start breaking through this everlasting issue with further comparative studies. As such, the expectation of the authors is that investigation of this research avenue would pave the way to the development of novel tools, methods, and reference databases that may help to overcome the limitations of both disciplines. (2) South African Black [11]; (3) Walker [7]; (4) craniometrics; (5) molecular method [37].

Conclusions
This paper presents the anthropological and genetic analyses for sex and ancestry estimation of 21 disarticulated crania that belong to unidentified migrants who perished in the Mediterranean shipwreck of 18 April 2015. This study highlights the drawbacks that challenge physical and molecular anthropologists when presented with such a peculiar set of remains. In some cases, physical and molecular analyses produced differing results (for sex estimation 8 out of 21 crania showed discrepant results, while for ancestry estimation 3 or 12 out of 21 crania showed discrepant results, according to the method) highlighting that the methods were still strictly linked to reference populations. Within the context of a migrant population, the application of a multimethod approach that combines molecular and physical evidence may provide valuable insights for elaborating the biological profiles of unidentified crania. However, according to the authors' experiences, a method to merge the two disciplines and to reach a univocal and conclusive result is still to be fully explored.

Data Availability Statement:
The data presented in this study are available upon request.