# Fast genomic analysis of aquatic bird populations from short single-end reads considering sex-related pitfalls

Pierre Faux, Jean C.P. Oliveira, Davidson Campos, Gisele P.M. Dantas, Thais

Augusta Maia, Camila Dergan, Pedro Cassemiro, Gisele Lobo Hajdu, José

Santos-Júnior, Fabrício Santos

► To cite this version:

Pierre Faux, Jean C.P. Oliveira, Davidson Campos, Gisele P.M. Dantas, Thais Augusta Maia, et al.. Fast genomic analysis of aquatic bird populations from short single-end reads considering sex-related pitfalls. Ecological Informatics, 2020, 56 (5), pp.101058. 10.1016/j.ecoinf.2020.101058. hal-04189517

HAL Id: hal-04189517

https://hal.science/hal-04189517v1

Submitted on 30 Aug 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC0 1.0 - Universal - International License
1 Fast genomic analysis of aquatic bird populations from short

2 single-end reads considering sex-related pitfalls

3

4 Running title: FAST ANALYSIS OF GBS READS OF BIRDS

5

6 Pierre Faux¹,*,a, Jean C. P. Oliveira¹, Davidson P. Campos¹, Gisele P. M. Dantas², Thais

7 Augusta Maia¹,², Camila G. Dergan¹, Pedro M. Cassemiro¹, Gisele Lobo Hajdu³, José E.

8 Santos-Júnior¹, Fabrício R. Santos¹,*

9

10 ¹Departamento de Genética, Ecologia e Evolução, Instituto de Ciências Biológicas,
11 Universidade Federal de Minas Gerais, Belo Horizonte, Minas Gerais, Brazil

12 ²Programa de Pós-Graduação em Biologia de Vertebrados, PUC Minas, Belo Horizonte,
13 Minas Gerais, Brazil

14 ³Departamento de Genética, Universidade do Estado do Rio de Janeiro, Rio de Janeiro,
15 Brazil.

16 *Address:* Room L3-244, Departamento de Genética, Ecologia e Evolução, ICB,
17 UFMG, Av. Antônio Carlos, 6627, 31270-010, Belo Horizonte, MG, Brazil

18 *ªPresent address:* 4ème aile rouge, UMR 7268, 27 Boulevard Jean Moulin, 13005
19 Marseille, France

20 *Corresponding authors:* Pierre Faux - pierrefaux@gmail.com, and Fabrício R. Santos -

21 fsantos@icb.ufmg.br
22

**Highlights**

23

* GBS may reveal sex biased genotypes that interfere with the population analyses

24

* A new pipeline is proposed to discard sequences derived from the avian W chromosome

25

26

* Discarding these W reads enhances subsequent population structure analyses

27

* Spotting these sequences allows an almost perfect determination of the samples sex

28

**Abstract**

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

Population genomics is transforming the analyses of biodiversity, particularly the studies in conservation genetics. However, genomic approaches based on short sequencing reads are often subjected to complex molecular protocols and bioinformatics procedures, which are prone to many artifactual biases. In preliminary analyses using common pipelines for handling, filtering and processing genotyping-by-sequencing (GBS) data from three aquatic bird species, we detected a bias likely associated to the sex of individuals of the genus *Sula*, which could be due to artifacts with sex chromosome derived reads. Thus, we have developed GBirdS, a new pipeline for avian population genetics from GBS data, implementing a new approach to identify molecular markers likely to be linked to the avian W chromosome. The pipeline was applied to a total of 113 samples from three avian species (*Sula leucogaster*, *Sula dactylatra* and *Mergus octosetaceus*), including 54 samples for which the molecular sexing was available. For 53 out of these 54 samples, the sex determination achieved using GBirdS agreed with molecular sexing. We have then removed from subsequent population analyses the sequences observed only in individuals determined as heterogametic (likely to be females) by our pipeline. Removing these sequences resulted in greater average fixation indexes between subpopulations of these three species. Therefore, we conclude that the
sequences removed by our pipeline had an artificial effect on interindividual differences and
population structure.

**Keywords:** Aquatic birds; Genotyping-by-Sequencing; Single Nucleotide Polymorphisms;
Bird sexual chromosomes

## 1. Introduction

Next-generation sequencing (NGS) has revolutionized the access to the structure and
variation of the genome and provided high throughput data that can be readily used for
population analyses. Among novel techniques of Reduced-Representation Sequencing (RRS,
[1]), Genotyping-By-Sequencing (GBS, [2]) has many advantages for population genetics
studies of non-model species, especially for wildlife populations and conservation genomics
([3]). As a quantitative indicator of its power, the GBS ratio between the number of markers
generated and the sequencing costs may be particularly high. The method is technically quite
straightforward to implement, and thanks to barcode tagging, it is possible to process
simultaneously (i.e. on the same sequencing lane) numerous samples at a time ([4]).

This technique, however, has a few drawbacks, mostly issues related to handling the large to
very large generated datasets. Introducing their new, computationally-effective pipeline
TASSEL-GBS, Glaubitz et al. (2014,[5]) were stating with great accuracy that "*the new*
bottleneck is the efficient bioinformatics analysis of the vast and ever-expanding sea of data".
Not only the prominent issue of low-coverage encourages development of innovative
bioinformatics pipelines (such as Stacks [6, 7] or PyRAD [8], applied to population studies of
many non-model organisms), but also the ignorance of the nature of short-length sequences
generated by "randomly" targeting the whole genome. In wildlife studies, little is known
about the karyotypic structure, the ploidy and the paralogy level between chromosomes of the
studied species. Moreover, in many cases, sequences generated by RRS are the first genome-
wide molecular information ever obtained. Among others, allele drop-out (ADO, [9]) is an example of how genotype calling (and subsequent population genetics analyses) is affected by the unavailability of other molecular knowledge about the studied species.

Sex chromosomes may also generate issues (e.g. miscalling sex-specific markers) of that type when using RRS techniques for population studies because of differences between heterogametic and homogametic sexes (e.g. female birds are ZW, and male birds are ZZ). For instances, the lack of karyotypic (and sex chromosomes) information could hamper the determination of the autosomal/non-autosomal origin of short-length sequence reads.

Moreover, in birds this issue can be serious due to a large disparity of sizes between sex chromosomes (Z and W), even considering closely related species ([10, 11]). Recent studies using RRS techniques (e.g. [12] with lizard species, [13] with rockfish species or [14] with fur seals) have resulted in new ways to identify sex-specific molecular markers. Nevertheless, in these studies, all the samples ([12, 13]) or at least a group of them ([14]) were previously sexed. The issue of determining the autosomal/non-autosomal origin of short-length reads would get even more problematic when the studied samples are not easily sexed, as some bird species displaying little sexual dimorphism. In such case, PCR-based sex determination (molecular sexing; see [15] and [16]) offers a way to distinguish males from females.

In this paper, we are dealing with three aquatic bird species expressing low (the Brazilian merganser *Mergus octosetaceus* Vieillot, 1817, and the masked booby *Sula dactylatra* Lesson, 1831) to moderate sexual dimorphism (the brown booby *S. leucogaster* Boddaert, 1783). None of these three species has a reference genome; moreover, little is known about their karyotypes (e.g. see [17] for *Sula sp.*). Population genomics can be even more intricate since only few of the *M. octosetaceus* and *S. leucogaster* samples were previously sexed (on the field and/or using molecular sexing [16]). Preliminary PCA analyses involving birds of
93 the genus *Sula* have suggested a likely bias associated to the gender, i.e., individuals of the
94 same sex had a propensity to cluster together (see results).

95 Hence to perform population analyses with RRS data accounting for potentially sex-linked
96 artifacts, we have developed *GBirdS*, a new pipeline for avian population genetics from GBS
97 reads, implementing a reverse approach described in this study. Instead of first identifying
98 molecular markers linked to sex chromosomes from ascertained samples and therefrom
99 assigning sex to samples and considering their sexes in population analyses, we first isolate,
100 from unsexed individuals, a group of linked sequences whose distribution pattern over
101 samples makes it likely to be sex-related. These likely sex-linked sequences are then
102 considered as irrelevant for downstream population analyses, which we compare to analyses
103 obtained by existing pipelines (Stacks [7] and iPyRAD – cited in [18]). Our main objective in
104 this study is to develop a new bioinformatic pipeline for RRS analyses that is able to discard
105 sequences linked to the W chromosome. The secondary objectives are to assess if the
106 subsequent inference of the (unknown) sex of the samples is correctly achieved by comparing
107 to molecular sexing and to assess how much the discarded sequences affect the estimation of
108 the population structure.

## 2. Material and Methods

110
2.1. Reduced-representation sequencing data from aquatic bird species from Brazil

111
2.1.1. Data collection

112 Three avian species were collected: two of them belong to the genus *Sula* (the brown booby *S. leucogaster*, the masked booby *S. dactylatra*; Order: Suliformes) and one belongs to the genus *Mergus* (the Brazilian merganser *M. octosetaceus*; Order: Anseriformes). The samples of *S. leucogaster* (43) and *S. dactylatra* (53) were collected from February to September 2015 in three Brazilian archipelagos of Trindade & Martin Vaz (TRI), Fernando de Noronha (NOR) and São Pedro & São Paulo (ASP), and samples of *M. octosetaceus* were collected between
118 2008 and 2016 in the states of Minas Gerais (Serra da Canastra - CAN, and Alto Paranaíba - ALP), Tocantins (Jalapão - JAL) and Goiás (Chapada dos Veadeiros - VEA); see Table 1 for more details.

121 All these 131 samples were sequenced (details in next section) and were then filtered based on rules for outliers detection among samples implemented in our pipeline (detailed in Supplementary Information SI2). Furthermore, two ASP samples of *S. leucogaster* were filtered out beforehand due to an unintended mixing of DNA samples during genomic library preparation, finally leaving 113 samples for bioinformatic analyses (Table 1).

126 **2.1.2. Generation of short single-end reads**

127 The GBS data was generated in Ecomol, a genomic service company from the ESALQ University (Piracicaba, SP, Brazil), following a protocol adapted from Elshire *et al*. (2011, [2]). Prior to sequencing, the integrity of the samples was evaluated by agarose gel electrophoresis (2%) to check if there was any genomic DNA degradation. Then, DNA purity and quantity were evaluated using Nanodrop 2000 and Qubit® 2.0 devices and following recommendations of the manufacturers. A total of 131 samples (see Table 1) were submitted to the GBS library construction protocol. All DNA samples were diluted to a concentration of 10 ng/µL. For enzyme digestion, 10 µL of DNA of each sample, 16 µL of water, 1 µL of *Pst*I enzyme (20 U/µL – New England Biolabs) and 3 µL of buffer were added separately. This solution was incubated at 37°C during 60 minutes and the confirmation of the digestion was performed using a 2% agarose gel.

138 The *Pst*I enzyme was chosen because it provided the best results in *in vitro* tests and used in previous studies (involving birds [19] and other animals [20, 21]). After digestion, different barcodes of varying length (5 to 9 bp) were used to tag each sample for posterior individual identification of reads. The barcode tagged DNA was then lyophilized in the vacuum centrifuge (45°C during 120 minutes) and ligated to adapters according to the GBS protocol
143

[2]: (1) addition of 1 µL of T4 Ligase enzyme (New England Biolabs), 5 µL of 10x buffer and 19 µL of water, (2) incubation to 22°C during 120 minutes, and (3) enzyme inactivation through incubation at 65°C during 30 minutes. After this, the samples were purified using the QIAquick PCR Purification Kit (Qiagen, CA) following recommendations of the manufacturer, and then amplified with primers complementary to the restrictions fragments and adapters that ligate the PCR products to the flow cell oligonucleotides. The PCR temperatures and times were as follows: 1x(72°C – 5 min; 98°C – 30s); 18x(98°C – 30s; 65°C – 30s; 72°C – 30s); 1x(72°C – 5 min). The amplified products were checked on a 2% agarose gel and purified with AMPure XP beads (Beckman-Coulter) in order to filter out small DNA fragments. The resulting GBS library was quantified by qPCR and with the Bioanalyzer (Agilent Technologies), and sequenced using the Illumina HiSeq 2000, generating ~100 bp single-end reads.

144

145

146

147

148

149

150

151

152

153

154

155

## 2.2. A pipeline for fast genotype calling from short single-end reads

156

### 2.2.1. Working hypotheses

157

Hereafter, we assume that no reference genome is available and that observed sequences may originate indifferently from nuclear (autosomes and sex chromosomes) and mitochondrial DNA. Also, we assume that sequences observed only in some samples from different sexes are likely to be sex-linked, depending (1) on the frequency of the pattern of observation of these sequences among the samples and (2) on the probability that such a pattern would be the actual sexing given a prior expected proportion of females ($p_{fem}$) in the studied population.

158

159

160

161

162

163

Furthermore, we consider that two identical restriction-site associated DNA sequences are both originated from the same locus if some conditions are fulfilled: (1) not being over-sequenced, (2) not having an extremely high or low GC content. These sequences are then compared to each other and only those showing a high degree of similarity are further considered as potential polymorphisms of the same locus. We will then call genotypes by

164

165

166

167
168

matching such sequences between them and considering those that have a high degree of
similarity as originating of the same putative loci.

169

170

In the course of our description, we will use the term “sequence” to refer to an actual
biological sequence of DNA and the term “read” to refer to the digital sequence of DNA
base-pairs produced by the sequencer.

171

172

173

### 2.2.2. Demultiplexing

174

The list of barcodes and their corresponding sequencing lane numbers and sample identifiers
are first loaded in memory. Each sequencing lane file (FASTQ format only) is then
sequentially read and reads are parsed in three parts: barcode, remnant of the cutting site and
sequence associated to the restriction site. Only reads with unambiguously identified barcode
and remnant of the cutting site are further considered; reads whose barcodes and/or remnant
of the cutting site include one or more unidentified base-pairs ("N") are trashed. However,
please note that a change of one base-pair is allowed in the barcode part. In a similar manner,
only associated sequences with all identified base-pairs are kept (first quality check; base pair
level). Once read, these sequences are stored per sample in binary format (see *Supplementary Information SI2*). Hereafter, we will refer to them as *reads* since they constitute the
meaningful part of the read.

175

176

177

178

179

180

181

182

183

184

185

### 2.2.3. Constitution of the library of reliable sequences

186

We proceed in two steps to uniquely sort the reads and quickly enumerate their numbers of
occurrences per sample. First, reads are uniquely sorted, sample after sample, into identical
sequences and the number of occurrences (hereafter referred to as the *coverage* of the
sequence) is counted down for each unique sequence. Secondly, all unique sequences of all
samples are pooled together and uniquely sorted. After sorting them, a coverage table
(sequences vs. samples) is produced and some quality checks are performed using that table,
first on sequences, then on samples (see *Supplementary Information SI2* and *SI3* for more

187

188

189

190

191

192
193 details). Sequences are discarded if either their GC-content or their maximum coverage per
194 sample is considered as outlier. Samples are discarded if their total amount of reads is lower
195 than a user-defined number (set to 500,000 by default) or if they share a low proportion of
196 identical sequences with most of other samples. The last check is achieved using a simple
197 heuristic and aims to remove samples originating from other species. The user may choose to
198 turn it off in case of interspecific studies. The remaining sequences are split in two libraries:
199 the first one gathers all sequences that were observed only once and in one sample (hereafter
200 referred to as *singleton*), and the second one gathers sequences observed in more than one
201 sample (hereafter referred to as *reliable sequences*).

202 Some files are written by the end of this step. A first file ("*reads_counts_per_sample*")
203 contains the counts of reads from both libraries (singleton and reliable sequences) and the
204 yield correction factor (*YCF*) per sample. The *YCF* accounts for the discrepancy in total
205 number of reads among samples (for *N*<sub>A</sub> and *N*<sub>min</sub> the total numbers of reads respectively of
206 sample *A* and of the sample with the least reads, *YCF*<sub>*A*</sub> = *N*<sub>*A*</sub>/*N*<sub>min</sub>). Only the reliable
207 sequences are used in the further steps. Writing these files allows the user to start new
208 analyses from this step, hence avoiding the most computationally-expensive steps
209 (demultiplexing and constitution of the libraries).

210 The .us files therefore contains **all sequences that were read at least twice: either with a**
211 **coverage of one in at least two samples, or with a coverage of two or more in a single sample.**
212 Before downstream analyses, we apply coverage filters on reliable sequences: they are kept if
213 their minimum coverage per sample and/or their minimum total coverage on all samples is
214 greater than a threshold (set to 3 for both, by default). It is worth noting that these filters apply
215 on reliable sequences before they are clustered in reciprocal pairs (i.e. before genotype calling
216 – see section 2.2.5 for more details) and will thus be overtaken by coverage filters applied to
217 genotype (COVmed and COVsin – see section 2.2.6).
### 2.2.4. Isolation and removal of sequences likely to be linked to chromosome W

Say a *sexing* is a possible determination of the sex of all samples from a given population. If the population contains *N* samples, $2^N$ sexing are possible. These sexing could be coded as logical *N*-long vectors (*true* is female; *false* is male), e.g. [*T*, *T*, *T*, *F*, *T*, ..., *F*, *F*, *T*]. Here, our aim is to heuristically determine which one of these $2^N$ possible sexing is the most likely. An exhaustive search among all possible sexing would however be computationally expensive, because $2^N$ may be very high (more than 1 million possible sexing when $N = 20$). But if we assume that at least one reliable sequence would originate from the chromosome W - we actually expect that there would have much more than only one -, the number of possible sexing to investigate is greatly reduced. Under this assumption, we only have to investigate the $N_P$ patterns of observation of reliable sequences among samples ($N_P \ll 2^N$ when *N* is large). A pattern of observation of a given sequence states in which samples the sequence was observed (*true*) or not (*false*), e.g. [*F*, *T*, *F*, *F*, *T*, ..., *T*, *F*, *T*].

To do so, we first list all the patterns of observation of reliable sequences among samples. Considering then each one of them separately as the true sexing, made of $N_F$ females and $N_M$ males ($N_F + N_M = N$), we compute:

1. the occurrence $O_i$ of the pattern (i.e. the number of reads in the library of reliable sequences that follow this pattern of distribution);

2. the probability $P_i$ that the population contains the $N_F$ females and $N_M$ males of that pattern, as the binomial distribution $B(N, p_{fem})$ where $p_{fem}$ is the expected proportion of females in population (set to 0.5 by default);

For the *i*-th pattern among $N_P$, we compute the score $S_i$ as $O_i \cdot P_i$. The top-scoring pattern is considered as the true sexing and the sequences whose distribution matches with that pattern are considered as linked to chromosome W. Those sequences are eventually discarded from
242

library. Besides, a score returns the strength of the signal captured by the highest score; it is
computed as 1 - S<sub>top2</sub>/S<sub>top1</sub>, where S<sub>top1</sub> and S<sub>top2</sub> are the top-1 and top-2 scores.

243

244

2.2.5. Genotype calling

245

All remaining sequences in a library are compared to each other using a framework similar to
that of UNEAK ([22]). In a few words, the rationale behind this framework is to search for
reciprocal pairs of sequences with high similarity; such pair is further considered as
originating from the same locus. "Reciprocal pair" means here that each member of the pair
has high similarity only with the other member and not with other sequences in library. The
resulting loci are thus bi-allelic. We consider two sequences as similar if they bear few single-
nucleotide polymorphisms (SNPs) or a short-length insertion-deletion (indel). The number of
SNPs allowed and the length of indels are user-defined; by default, we respectively set them
to 2 and 9. Our implementation is based on our own comparison algorithm, which is
specifically conceived for short-length sequences: the algorithm has quadratic and linear
complexities respectively with the number of SNPs allowed and the length of indels (more
details in Supplementary Information SI2).

246

247

248

249

250

251

252

253

254

255

256

257

2.2.6. Genotypes filtering and kinship computation

258

Filters are then applied to the set of called genotypes. We discard markers (1) if they are
monomorphic, or (2) if the median coverage of one of their alleles is lower than a user-
defined threshold COVmed (set to 6 by default), or (3) if the coverage of the minor allele is
lower than a user-defined threshold COVsin (set to 10 by default). The last rule applies only
to markers whose minor allele is a singleton (the allele was observed only in one sample
among the population, in a heterozygous site); we consider that in such case the possibility
that a singleton allele would be an error is more likely than the opposite (not an error). For all
remaining markers, additional statistics are computed: sample call rate (CR), minor allele
frequency (MAF), sample minimal coverage, *p*-value from Hardy-Weinberg Equilibrium test

259

260

261

262

263

264

265

266
and, given a population map, *F*<sub>ST</sub> estimation by the third formula of Weir and Cockerham (1984, [23]). Among these genotypes, those with CR and MAF greater or equal than user-defined values (respectively 0.7 and 0.02 by default) are kept apart to compute the kinship matrix using GBS with depth adjustment (KGD). The KGD matrix was proposed by Dodds *et al.* (2015, [24]) to take uncertainty linked to low coverage into account when computing kinship coefficients.

### 2.2.7. Code availability

The main body of the pipeline (demultiplexing, sort of reads into unique sequences per sample and for the whole population, quality checks, identification of sex-linked sequences, genotype calling and filtering) is wrapped in a single program written in Fortran 95 and the final steps (computation on outputs: matrix, fixation index and principal component analysis) are wrapped in a Matlab® function. All sources and executable are freely available on GitHub (https://github.com/pierrefaux/gbirds.git).

## 2.3. Comparisons to other pipelines

In order to assess the behavior of our new pipeline, we have also processed the GBS raw data with two state-of-the-art pipelines: Stacks ([6, 7], version 1.44) and iPyRAD (cited in [18], version 0.6.27). Both of them are used with default values. Samples discarded by one of three pipelines are also discarded from other pipelines in order to compare computational requirements.

For comparison purposes, we set up four different selections of markers (labeled SI to SIV) for each of the three pipelines:

288 - SI is the largest selection: it includes all genotypes called by a given pipeline;

289 - SII gathers all markers whose COVmed and COVsin are respectively greater or equal to 6 and 10;
291

- SIII gathers all markers of SII whose CR is greater or equal to 7 and whose MAF is greater or equal to 2%;

292

293

- SIV gathers all markers of SII whose MAF is greater or equal to 2% and that are observed in all samples (CR equal to 1).

294

295

Selections SI and SII are set up to quantify the impact of filters (COVmed and COVsin, MAF and CR) on the actual number of markers returned by our pipeline (markers actually found in SIII).

296

297

298

Then, using selections SIII (which corresponds to the default filtering rules of our pipeline) and SIV, we perform comparisons between pipelines on their estimates of the population structure. In first place, we compute various correlations between pipelines using KGD matrices built on SIII (as in our pipeline): (1) on both diagonal and off-diagonal values of the KGD matrix and (2) on the two first principal components of the KGD matrix. In second place, selection SIV (SIII with no missing genotype) is used to structure each population: using Admixture ([25]), we determine the best number of sub-populations (*K*) for each pipeline and each species, and check if the samples are clustered in the same sub-populations throughout pipelines. Admixture finds out which *K* has best predictive accuracy by applying a 5-fold cross-validation procedure, for values of *K* ranging from 2 to 5. In third place, we estimate the differentiation between the sub-populations by computing the fixation index *F*<sub>ST</sub> (using the third formula of Weir and Cockerham, 1984 [23]).

299

300

301

302

303

304

305

306

307

308

309

310

Besides, the usage of computational resources of each pipeline, for each species, are further recorded: we take note of the computing time and of the memory peak. In order to avoid very long runs, the number of central processing units (CPUs) is requested before running the jobs: 16 CPUs for both Stacks and iPyRAD, as both allow parallelized work, and 1 CPU for GBirdS, as parallelization is not implemented in our pipeline. The memory peak refers to the maximum usage of random-access memory (RAM), measured as the maximum of the

311

312

313

314

315
316 resident set size (RSS) of the job (when only 1 CPU is requested – as for GBirdS), or as the
317 maximum of the sum of the RSSs of each CPU (when 16 CPUs are requested – as for Stacks
318 and iPyRAD). The computational comparisons take place on Intel® Xeon® Gold (clock
319 speed: 2.60 GHz, cache: 22 Mbytes and RAM: 192 Gbytes) computing nodes.

320 **2.4. Determination of the sex of the samples (sexing)**

321 On-field and/or molecular sexing were not previously available for all the 113 samples kept
322 after filters for the population studies. Among the 36 *S. leucogaster* samples kept, 24 of them
323 originating from São Pedro & São Paulo Archipelago (tagged “ASPxx”) were sexed both on
324 the field (based on beak measurements) and using molecular sexing (see below). Also, only
325 30 out of 31 *M. octosetaceus* samples kept after filters were sexed using molecular sexing,
326 and only 19 of these 30 samples were also sexed on the field (by cloaca eversion). Thus, in
327 order to estimate the sex of non-sexed samples, we have also aligned the library of reliable
328 sequences to ten avian Z chromosomes and computed coverage statistics from perfectly
329 matching sequences. Protocols and computations are described here below.

330 **2.4.1. Molecular sexing**

331 To perform the PCR molecular sexing, we used primers for the CHD gene (coding for
332 chromodomain helicase DNA binding protein), first described on the W chromosome (CHD-
333 W) of birds ([15]), and later found its homologue on the Z chromosome (CHD-Z) ([16]). The
334 CHD genes are highly conserved among different animal species ([26]), but the intronic
335 regions are less conserved and present different sizes ([16]). The different sizes of Z and W
336 CHD homologues makes possible the sexual identification of bird samples through
337 polymerase chain reaction (PCR), followed by visualization of the region amplified in agarose
338 gel electrophoresis, since the females (ZW) would form two bands in the gel, whereas the
339 males (ZZ) would form one single band.
First, all 54 available samples (30 *M. octosetaceus* and 24 *S. leucogaster*) had their DNA extracted following the phenol-chloroform protocol ([27]). The integrity of the samples was evaluated by agarose gel electrophoresis (2%) where it was checked if there was any genomic DNA degradation. DNA purity and concentration were evaluated on the Nanodrop 2000 and Qubit® 2.0 devices, following recommendations of the manufacturers. After these steps the samples were submitted to PCR amplifying the CHD gene with the primers CHD1F (5' - TATCGTCAGTTTCCTTTTCAGGT - 3') and CHD1R (5' - CCTTTTATTGATCCATCAAGCCT - 3') described by Lee *et al*. (2010, [28]) and the result was visualized on a 2% agarose gel electrophoresis.

### 2.4.2. Alignment of reliable sequences to avian Z chromosomes

The full sequences of the Z chromosomes of 10 avian species were downloaded from the NCBI genome repository (ftp://ftp.ncbi.nlm.nih.gov/genomes; see *Supplementary Information SI6* for accession numbers). These ten species belong from different avian orders: one (*Anas platyrhynchos*) belongs to the same order as *M. octosetaceus* (Anseriformes), four of them are Galliformes (*Coturnix japonica*, *Gallus gallus*, *Meleagris gallopavo* and *Numida meleagris*), one (*Athene cunicularia*) belongs to Strigiformes and the remaining four are Passeriformes (*Cyanistes caeruleus*, *Ficedula albicollis*, *Parus major*, *Taenipygia guttata*). The libraries of reliable sequences of the three species (hence, including reliable sequences likely to be linked to chromosome W) were aligned against the ten Z chromosomes using *bwa* ([29]). For each of the three species, we kept apart any sequence perfectly matching (whatever the read sense) with all the ten Z chromosomes.

Here we postulate that any sequence complying with such a condition may be a Z-segment conserved through avian lineages. If this would be true and if all samples were sequenced at an equal and high coverage, we would expect that coverage to be, on average, twice as high for males (ZZ) in comparison to females (ZW). This assumption may unfortunately be
affected by various effects, among others: (1) by the low coverage, i.e. the fact that a sequence was not observed in a sample does not mean that this sample does not have this Z-linked segment; (2) by the variation in total number of reads among samples, i.e. the fact that a Z-linked segment has a high coverage in a sample may be explained by a higher coverage that the average for that sample, rather than by the sex of that sample; and (3) by polymorphisms between the Z chromosomes that a male sample carries, i.e. a male carrying haplotypes Z1/Z1 is expected to have twice the coverage of a male carrying Z1/Z2.

To account for some of these effects, the number of reads of all sequences likely to be a conserved Z-segment was summed for any sample and then divided by the YCF of the sample (this factor accounts for the variation in total number of reads among samples – see section *Constitution of the library of reliable sequences*). Samples are finally ranked according to that corrected coverage. Males are expected to cluster at the top of the ranking and females at the bottom, the average values of males being twice the average value of females.

378
379

# 3. Results

380

## 3.1. Preliminary analyses suggesting sex-biased artifacts

381

We have started our study by processing the GBS data in Stacks (v.1.44) with reference

382

genomes from two closely related species (the great Cormorant – *P. carbo* – for *Sula* sp. and

383

the mallard duck – *A. platyrhynchos* – for *M. octosetaceus*), and performing preliminary

384

population analyses with a PCA clustering approach. Processed GBS data derived from *Sula*

385

*leucogaster* and *S. dactylatra* generated PCA graphs showing individuals forming groups

386

separated by islands (as expected), but also two groups within each island clustering only

387

females or only males (graph not shown). This suggested a likely bias associated to gender,

388

i.e., individuals of the same sex had a tendency to cluster together.

389

This hypothesis was further supported by the alignment analysis of GBS reads. The two

390

sequences (#1 and #2) of *S. leucogaster* samples outlined in Supplementary Information SI1

391

are aligned against the reference genome of *P. carbo*. Both sequences had high matching

392

proportion with the same locus on the reference genome (SI1a). Therefore, processing these

393

alignments with Stacks led to consider the two sequences as two alleles of the same locus. For

394

24 samples who were sexed on field, that locus was homozygous in males and heterozygous

395

in females (SI1b). As bird females are the heterogametic sex (ZW), a likely explanation

396

would be that the sequence #1 originated from the W chromosome and was incorrectly

397

associated to a paralog located on the Z chromosome (sequence #2). That hypothesis was

398

further corroborated by the sequencing depths (SI1c): sequence #2 (putatively originating

399

from Z) was 1.7x more covered in males than in females.

400

## 3.2. Output files produced by GBirdS

401

Once the original reads are uniquely sorted in sequences (per sample and for the whole

402

population), two binary files are outputted per sample; they store the libraries of singletons

403

and reliable sequences respectively with suffixes .usg and .us. The sizes of these binary files
404 are quite small (up to a few tens of Mbytes), making a new analysis from original data quite
405 easy.

406 After finishing, the pipeline returns two libraries of sequences (in FASTA format): the
407 sequences remaining after all cleaning steps and the sequences likely to be linked to W
408 chromosome. These two libraries come along with two space-delimited files of same length
409 containing, for each sequence, the coverage (number of reads) of each sample. Other outputs
410 are a detailed list of all markers, the genotypes, the KGD matrix, and the projections of
411 samples along the two first principal components of the KGD matrix. More details on the
412 contents and formats of these files are given at the end of Supplementary Information SI2.

### 3.3. Isolation of sequences likely to be linked to chromosome W

414 For each studied species, the pipeline sorted all the different patterns of observation of
415 reliable sequences among samples (considering sequences with a minimum coverage of 3). In
416 line with the assumption detailed in section 2.2.4, the pipeline assumed that the top-scoring
417 pattern shows the actual sexing. The sequences pertaining to that pattern of observation were
418 kept apart, as they were considered to be linked to chromosome W. For each species, the
419 isolated pattern was much more frequent than the average: from Table 2, it was from 360x
420 (*M. octosetaceus*) to 1035x (*S. dactylatra*) more frequent than the average. These patterns
421 were observed respectively in only 19, 23 and 13 samples for *S. leucogaster*, *S. dactylatra* and
422 *M. octosetaceus*, suggesting those samples would share sequences that originate from
423 chromosome W, hence would be females.

424 The proportion of these sequences per sample (computed as the total number of reads of the
425 sequences matching to the pattern divided by the total coverage of all reliable sequences with
426 a minimum coverage of 3) was computed and is plotted in Figure 3. The intraspecies
427 distributions have small standard deviations (average proportion ranging from 2.6E-3 to 7.7E-
428 3 and standard deviations from 3.5E-4 to 5E-4) and mirror the number of sequences matching
429

the pattern for each species (lowest number for *M. octosetaceus*, and highest number for *S.*

430

*dactylatra* – see Table 2).

431

**3.4. Comparisons to other pipelines**

432

**3.4.1. Genotype calling**

433

Comparisons in terms of “marker yield” are given in Table 3. Stacks outperforms our pipeline

434

(GBirdS) and iPyRAD in number of unfiltered markers (selection SI), even though it calls

435

only SNPs and not indels. Among species, Stacks returns from 1.15x to 1.33x more markers

436

than GBirdS.

437

The number of markers returned by GBirdS (selection SIII) is on average 4.8x lower than the

438

initial number of markers. This factor is greater than for other pipelines (~3.3x for Stacks and

439

~3.9x for iPyRAD). Another difference between pipelines is related to the markers discarded

440

from the initial selection (SI): in our pipeline, ~96% of the markers were discarded by rules

441

on COVmed and COVsin while this is the case for ~82% and ~72% of the markers with

442

Stacks and iPyRAD, respectively. The number of markers that comply to all filtering rules

443

(SIII) remains ~1.6x higher with Stacks, when compared with iPyRAD or GBirdS. However,

444

it is worth noting that there is a greater proportion of these markers that were not observed in

445

all samples with Stacks (~72% vs. ~64% with iPyRAD, and ~57% with our pipeline),

446

eventually leading to similar numbers of markers observed in all samples and complying with

447

all filtering rules (selection SIV).

448

**3.4.2. Comparisons using Admixture**

449

Prediction errors returned by Admixture (see Table 4) show that all the three pipelines agree

450

on splitting the populations of *S. leucogaster* and *S. dactylatra* in two sub-populations (*K*=2).

451

In both cases, these subpopulations match with the geographic origins. However, for *M.*

452

*octosetaceus*, Stacks finds out that the best *K* is 2 whereas both iPyRAD and GBirdS suggest

453

*K* equal to 3 as the best model. Regarding geographic origins (referring to Table 1), all

454

pipelines cluster together the *M. octosetaceus* JAL and VEA samples on the one hand and the
455 CAN samples on the other hand; ALP samples are either all joint to the CAN cluster (with
456 Stacks), either distributed between the CAN cluster and a third cluster, only made of ALP
457 samples (with iPyRAD and GBirdS).

458 Prediction errors are always lower with Stacks. In the case of *S. leucogaster*, the Admixture
459 model with *K* equal to 3 using Stacks has a prediction error lower than the lowest ones
460 obtained by iPyRAD or GBirdS.

461 **3.4.3. Comparisons using KGD matrix**

462 Looking at how first principal components of the kinship matrices correlate between pipelines
463 is another way to compare them, on their estimations of the population structure. In all the
464 three populations, samples are ordered in a highly similar manner along the first PC (in Table
465 6, correlations always greater than 98.9%) by the three pipelines. Differences between our
466 pipeline and the two other pipelines arise when ranking the samples along the second PC for
467 *S. dactylatra* (results suggest that the second PC returned by GBirdS would be the third PC
468 returned by Stacks and iPyRAD), or even further (third PC) for *S. leucogaster*.

469 Regarding kinship coefficients (see Table 6), estimates returned by the three pipelines are
470 very similar in the case of *S. leucogaster* (correlations greater than 99.4%), moderate in the
471 case of *M. octosetaceus* and low in the case of *S. dactylatra*. In the former case, correlations
472 of diagonal elements are particularly low, thus suggesting that estimation of inbreeding
473 coefficients would result in high differences between pipelines.

474 **3.4.4. Average fixation indexes**

475 Greater values of average fixation index *F*<sub>ST</sub> are always obtained with our pipeline (Table 5).
476 Estimations of *F*<sub>ST</sub> between subpopulations are nonetheless very similar from a pipeline to
477 another.

478 In terms of population structure, the most differentiated subpopulations within species are
479 those of *S. leucogaster* (NOR and ASP; see Table 1 for population details), exhibiting *F*<sub>ST</sub>
480 values ranging from 21% (from iPyRAD) to 25% (from GBirdS). Conversely, the closest

481 subpopulations are those of *S. dactylatra* (*F*<sub>ST</sub> ~2-3%). Among *M. octosetaceus*

482 subpopulations, the closest ones are CAN and ALP (*F*<sub>ST</sub>~6-9%), followed by the CAN-VEA

483 pair (*F*<sub>ST</sub> ~8-11%).

484 **3.4.5. Markers and PCs matching gender distribution**

485 The main specificity of our pipeline is to spot markers that most likely originate from

486 chromosome W and to filter them out. Conversely, some markers that perfectly match with

487 the pattern of observation of W-linked sequences by GBirdS have gone through the filters

488 with Stacks and iPyRAD (Table 7). For instance, almost 1% of the markers returned by

489 Stacks and used for the computation of KGD matrix (selection SIII) show perfect association

490 with the pattern in the case of *S. leucogaster*. These markers may thus affect the population

491 structure analyses by clustering samples according to that pattern of observation of sequences.

492 Correlations between first PCs and the pattern (see Table 7) suggest that such a clustering

493 would be found in *S. dactylatra*: the second PC correlates moderately (Stacks) to very well

494 (iPyRAD) with that pattern. That effect is clearly visible when plotting the samples against

495 the two first PCs obtained using Stacks or iPyRAD (see Figure 2). The same effect may be

496 observed, to a lower extent though, in *S. leucogaster* (the pattern is moderately reflected in the

497 third PC – see Supplementary information SI4 and SI5).

498 **3.4.6. Computational requirements**

499 As well for computing time as for memory usage, GBirdS outperforms Stacks and iPyRAD

500 (see Table 8): with 1 CPU, GBirdS runs in less than 1 hour and requires maximum 8.1

501 Gbytes. In terms of computing time, with 16 CPUs, the fastest iPyRAD run (*S. leucogaster*)

502 finishes in more than 17 hours while runs of Stacks take between 7 and 11 hours. In terms of
503 memory usage, iPyRAD requires more or less twice the memory needed by GBirdS and

504 Stacks has the highest peaks (36 to 52 Gbytes) among the three pipelines. However, this peak
505 is due to a single and short (about a few minutes) step in the Stacks pipeline (the
506

"populations" command); besides this demanding step, the memory usage of Stacks is comparable to that of GBirdS.

507

508

The computing time required by GBirdS mainly covers the first steps of the pipeline (Figure 1(a): sequencing lane parsing, quality checks and unique sort of sequences per sample).

509

510

Among that steps, the longest time was required by the lane parsing (including the first level quality check on base-pairs Phred score). That step took from 16 (*M. octosetaceus*) to 39 minutes (*S. dactylatra*) with GBirdS, whereas the comparable step was more rapidly performed by iPyRAD for all three species (maximum 15 minutes), and achieved in similar timings by Stacks.

511

512

513

514

515

### 3.5. Determination of the sex of samples

516

Molecular sexing of 30 *M. octosetaceus* and 24 *S. leucogaster* samples has resulted in identifying 13 females and 17 males in *M. octosetaceus*, and 11 females and 13 males in *S. leucogaster*. That result was in line with on-field sexing for the 43 samples sexed on the field.

517

518

519

The sexing results are detailed in Supplementary Information SI7, alongside with the observation/non-observation of sequences likely to be linked to chromosome W. For *S. leucogaster*, the isolated sequences found to be likely linked to chromosome W were all sequenced in the 11 females and none of them was sequenced in the 13 males. For *M. octosetaceus*, the isolated sequences found to be likely linked to chromosome W were all sequenced in 12 females out of 13, and none of them was sequenced in the 17 males. The mismatching sample was sexed only by molecular sexing (not on field) and the molecular sexing was confirmed after replication for that particular sample. Therefore, the pattern of observation of the sequences likely to be linked to chromosome W has a mismatch rate of 1/54 (2%) with molecular sexing (0% with on-field sexing).

520

521

522

523

524

525

526

527

528

529

The alignment of the library of reliable sequences against the ten avian Z chromosomes resulted in identifying few hundreds of sequences showing whole-length sequence match with

530
all these Z chromosomes: 206 for *S. leucogaster*, 271 for *S. dactylatra* and 299 for *M.*

*octosetaceus*. Considering the on-field/molecular sexing as the true sexing for the 54 tested

samples, the average corrected coverage of these Z-matching sequences was, respectively for

males versus females: 285.3 ± 13.0 vs. 188.6 ± 12.6 (hence, 1.51x higher for males) for *S.*

*leucogaster*, and 183.7 ± 22.1 vs. 117.2 ± 15.2 (hence, 1.57x higher for males) for *M.*

*octosetaceus*. In Figure 4, all samples are plotted along the axis of corrected coverage in Z-

matching sequences; a plain-blue circle indicates that the sequences likely to be linked to

chromosome W were all observed in that sample, conversely, an open-red circle indicates that

none of these sequences was observed in that sample. The tested samples are annotated with

their on-field/molecular sexing on the left-hand side. The two mismatching females (one in

both tested species) visually cluster with males, though they remain at the lower-edge of the

male clusters. Similarly, the *S. leucogaster* mismatching male cluster with the females, at the

upper-edge of the female clusters. Considering the observation/non-observation of sequences

likely to be linked to chromosome W as the true sexing, the average corrected coverage of the

Z-matching reads turns to (males versus females – computed using all samples): 281.6 ± 14.3

vs. 183.0 ± 22.0 (hence, 1.54x higher for males) for *S. leucogaster* and 181.8 ± 23.0 vs. 114.6

± 12.3 (hence, 1.59x higher for males) for *M. octosetaceus*.

548
549 **4. Discussion**

550 **4.1. Development of a new GBS pipeline: learnings, limitations and perspectives**

551 The proposed pipeline achieved comparable population analyses as those of state-of-the-art
552 pipelines with a much lesser usage of computing resources, making it useful for personal
553 computers (low amount of memory required, speed of analyses) and for various types of RRS
554 experiments (e.g. RAD-seq, ddRAD-seq; not tested yet). This efficiency in terms of
555 computational resources is mainly imputable to the simplicity of the implemented methods
556 (for sequences clustering, for sex identification and for genotype calling) because the most
557 time-consuming step of our pipeline is the parsing of sequencing lane, a step that is achieved
558 either in comparable times or in shorter times by the two other pipelines. Before scaling with
559 number of samples, that parsing time actually scales with the number and the length of the
560 sequencing lanes: the *M. octosetaceus* samples all come from only one sequencing lane
561 (containing 192 M of reads and parsed in 16 minutes) while the *Sula sp.* samples came from
562 three sequencing lanes (containing from 216 to 230 M of reads and parsed in minimum 38
563 minutes). In terms of peak of memory usage, our pipeline linearly scales with the total
564 number of reads: through the three species, we have observed a constant rate of ~37.8 Mbytes
565 per million of reads kept in memory. This rate is in line with expectations: at the peak of
566 usage, the two main variables allocated have the total number of kept reads as first dimension
567 and the variable in third rank (according to memory usage) have the number of unique
568 sequence as first dimension, which scales almost linearly with the total number of reads. At
569 that stage of the pipeline, the number of samples affects the variables in fourth rank of
570 memory usage and beyond. These empirical rules have yet to be confirmed but may help
571 future users to check for applicability of the pipeline on low-resources platforms.

572 The main limitation of our pipeline also relates to its simple theoretical background: we only
573 aim at sorting reads, gathering them into unique sequences, identifying patterns of distribution
of these sequences and clustering them into reciprocal pairs. This simplicity, which was
inspired by that of UNEAK [22], makes it low flexible (unable to handle overlapping reads;
does not work with paired-ends data), statistically poor (no statistical conclusions are drawn
on genotype calling) and not as efficient as Stacks and iPyRAD in terms of data usage
(singleton sequences are trashed by GBirdS while they are aligned against consensus
sequences by Stacks and iPyRAD). Therefore, our new pipeline is not meant to replace
existing, powerful, widely used pipelines such as Stacks and iPyRAD but instead, it should be
used to check for the presence of a sex-related pattern before using these pipelines. Please
note another current limitation: it was designed to work with avian species and has not yet
been tested on other animal species, although it should be applicable to all species with
differentiated sex chromosomes.

Though the sex identification method was tested on two different species (S. *dactylatra* being
kept apart as no molecular/on-field sexing was available for that species) from two different
avian orders, our algorithm has yet to be assessed, e.g. through simulation studies. The
following parameters are worth testing: (1) other species (not only avian species, as already
mentioned), including variations in size of the smaller sex chromosome; (2) various sizes of
population, because the sex identification is shape in a way that it requires a minimal number
of samples (of both sexes) to efficiently identify the sex of samples; and (3) various types of
restriction enzymes (as their distribution over the chromosomes may affect sex identification
in our method). As proposed by a reviewer, it would also be worth simulating the proportion
of short (<100 bp) reads originating from both sexual chromosomes that cluster together
based on known genomes.

Two major features should be added to the pipeline in the near future. First, the most relevant
may surely be the parallelization of the code. Lane parsing would be the main step to
parallelize since it is the most time consuming. Various solutions could be applied in GBirdS:
599 (1) reading various lane files simultaneously, (2) reading the same lane file with various file
600 units, or (3) a simultaneous translation (from characters to integers) of various parts of the
601 same read. Secondly, the possibility to include prior information about the sex of some
602 samples: weighting the score of each pattern by its concordance with the prior information
603 would be the easiest way to integrate that information.

604

605
606

## 4.2. Identification of female samples using GBirdS may be affected by the size of the W chromosome

607 In the Results section of this paper, we have outlined some elements showing that the isolated
608 sequences are, at least partially, linked to sex; hence, isolating these sequences subsequently
609 allows inferring the unknown sex of samples. On the one hand, most sexed samples agree
610 with the sexing based on isolated sequences (98% of them). On the other hand, the
611 mismatching samples cluster accordingly to that GBirdS sexing when considering their
612 coverage in sequences perfectly matching to ten avian chromosomes Z (in Figure 4, the
613 variance of corrected coverage grouped by sex are lower with GBirdS sexing than with on-
614 field/molecular sexing, except for M. octosetaceus males). However, the GBirdS sexing
615 approach remains imperfect: for one sample (M. octosetaceus) out of 54, the sexing achieved
616 by GBirdS disagrees with molecular sexing (please note that this sample was not sexed on-
617 field). Interestingly, the mismatching sample is located on the edge of the male cluster in
618 Figure 4, indicating that a misclassification using that technique is more likely for this
619 mismatching sample than for the samples located at the opposite extremes of the clusters. It is
620 also worth noting that this mismatching sample is a female (according to molecular sexing)
621 who was considered as a male by GBirdS. In other words, no GBS sequences likely to be
622 linked to W was observed in that sample. Furthermore, the number of these sequences are
623 much lower in M. octosetaceus than in *Sula sp*. (see Table 2), which does not lead to a visible
624 artifact in population structure of the Brazilian merganser (see Table 7: no clear association
625 between the 5 first PCs and the pattern distribution). Nevertheless, this reveals a limitation of
626 sexing individuals based on GBirdS using RRS data: in cases where few sequences originate
627 from the W chromosome, a female individual could be inferred as a male.

628 The number of artifactual sequences can also be implicitly related to the actual size of the W
629 chromosome in a given species. Assuming that the filtered sequences are effectively linked to
630 chromosome W, the proportion of the genome they cover (Figure 3) could be interpreted as an
631 estimation of the relative sizes of the chromosome W between each species. Though
632 inaccurate (the estimation relies on the assumption that the restriction enzyme uniformly
633 targets all genomic regions), these estimations are in line with reported lengths of W
634 chromosomes of related taxa (great cormorant for *Sula sp.*; mallard for *M. octosetaceus*, [30]).

635 **4.3. Discarding sequences likely to be linked to W chromosome enhances population**
636 **structure analyses in *Sula* sp.**

637 Our study highlights the presence of an artifactual structuration pattern, apparently biased
638 towards the gender of samples in two aquatic bird species of the *Sula* genus from Brazil:
639 some of the first principal components and the distribution of inferred sex (Table 7) are
640 moderately (*S. leucogaster*) to highly correlated (*S. dactylatra* – reflected in Figure 2.a-b).

641 Conversely, in *M. octosetaceus*, a gender-biased structuration pattern was never observed.

642 For the *Sula sp.*, removing the sequences matching the sex-related pattern enhances
643 population analyses: although Admixture analyses do not show large differences between
644 pipelines, the fixation index results (see Table 5) are always greater for the genotypes called
645 by GBirdS than for those called by Stacks or iPyRAD. We interpret this result as an evidence
646 that the signal captured is a reproducible artifact, which is not related to the geographic
647 localizations of samples, but to their sex. In worst cases, the presence of that artifact may
648 greatly affect the population study. The case of *S. dactylatra* is quite illustrative: on Figure 2b
649 (results from iPyRAD, marker selection SIII), the second PC classifies neatly the samples
650 carrying sex-related sequences versus those not carrying such sequences. A confirmation is
given in Table 7: the correlation between the distribution pattern and the second PC is close to -1. The case is very similar for the results obtained with Stacks. One can find a first explanation in the same table: 0.7% of the markers in the iPyRAD selection have their distribution perfectly matching with the sex-related pattern. An investigation on how these markers were called reveals that in most cases, those markers were shaped by grouping sequences that most likely come from both sex chromosomes, hence calling as hetero/homo-zygous the hetero/homo-gametic samples.

Removing the sequences linked to chromosome W can therefore improve downstream analyses. However, it may be viewed as incomplete: indeed, if these sequences are genuinely originating from chromosome W, then cleaning the data by discarding them may also involve discarding reads from chromosome Z. Unfortunately, spotting reads originating from these sequences in a *de novo* approach (without using any reference genome), *streamlined way cannot be achieved*.

## 5. Conclusions

In this study, we introduce *GBirdS*, a new pipeline for fast genotype calling from short-length, single-end, barcoded reads for species without reference genomes. The main purpose of this pipeline is to isolate a pattern of observation/non-observation of sequences likely to be linked to chromosome W and discard them for further genotype calling and population studies.

Compared to state-of-the-art pipelines, the main advantage of GBirdS lies in its computational performances: GBirds is faster and uses less memory because it does not use singleton sequences and rely on a simple model for pairing similar sequences.

Besides, the isolated pattern of observations of sequences likely to be linked to chromosome W can be used to infer the unknown sex of samples. Compared to molecular sexing, this "GBS-based" sexing was found to be correct for 53 samples out of 54. Therefore, we assume
675 that the ability of GBirdS to correctly infer the sex of samples is taxon specific, actually better
676 for taxa that have larger W chromosomes, or that present higher number of W-derived reads.
677 For such species, it was found that discarding the sequences likely to be linked to
678 chromosome W enhances the population structure analyses; these sequences work as an
679 artifact, masking inter-individual and inter-population genetic differences.

680 **6. Acknowledgements**

681 P.F. acknowledges the Graduate Program of Genetics of UFMG (Belo Horizonte, Brazil) and
682 Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) of Brazil for
683 funding, through the Programa Nacional de Pós Doutorado (PNPD), the postdoctoral research
684 in which this study took place. Data sets were processed in Sagarana HPC cluster, CPAD-
685 ICB-UFMG, and we thank Professor Miguel J. Ortega for technical assistance. This project
686 was funded by the grant 458581/2013-5 CNPq (Brazil) with logistic support from Marinha do
687 Brasil. Sampling permit for sea birds was issued by SisBio-ICMBio 45906-1 and ethical
688 authorization by CEUA-PUC Minas (47/2014). Samples for the Brazilian merganser are part
689 of the collection of the Brazilian National Action for conservation of ICMBio, or obtained
690 with SisBio-ICMBio permits 17051 and 12662, with ethical procedures approved by CEUA-
691 UFMG (202/2007).

692 **7. Author contributions**

693 PF developed the software, analyzed the data and wrote the first version of the manuscript;
694 JCPO, GPMD, TAM, CGD and JESJ collected the sea bird samples, PM and DPC did
695 molecular sexing of samples, GPMD and DPC accompanied the production of GBS library at
696 Ecomol/ESALQ; DPC analyzed the *Mergus* data; JCPO, JESJ and FRS drafted the
697 discussion; FRS coordinated the study, analyses and writing of the manuscript; all authors
698 read and approved the final version of the manuscript.

699

700 **8. Declaration of interests**
701 The authors declare that they have no conflict of interest.

702 **9. References**

703 1. Davey JW, Hohenlohe PA, Etter PD, Boone JQ, Catchen JM, Blaxter ML. Genome-wide genetic marker discovery and genotyping using next-generation sequencing. Nat Rev Genet. 2011;12:499–510.

706 2. Elshire RJ, Glaubitz JC, Sun Q, Poland JA, Kawamoto K, Buckler ES, et al. A Robust, Simple Genotyping-by-Sequencing (GBS) Approach for High Diversity Species. PLoS ONE. 2011;6:e19379.

709 3. Narum SR, Buerkle CA, Davey JW, Miller MR, Hohenlohe PA. Genotyping-by-sequencing in ecological and conservation genomics. Mol Ecol. 2013;22:2841–7.

711 4. Andrews KR, Good JM, Miller MR, Luikart G, Hohenlohe PA. Harnessing the power of RADseq for ecological and evolutionary genomics. Nat Rev Genet. 2016;17:81–92.

713 5. Glaubitz JC, Casstevens TM, Lu F, Harriman J, Elshire RJ, Sun Q, et al. TASSEL-GBS: A High Capacity Genotyping by Sequencing Analysis Pipeline. PLoS ONE. 2014;9:e90346.

716 6. Catchen JM, Amores A, Hohenlohe P, Cresko W, Postlethwait JH. Stacks: Building and Genotyping Loci De Novo From Short-Read Sequences. G3amp58 GenesGenomesGenetics. 2011;1:171–82.

719 7. Catchen J, Hohenlohe PA, Bassham S, Amores A, Cresko WA. Stacks: an analysis tool set for population genomics. Mol Ecol. 2013;22:3124–40.

721 8. Eaton DAR. PyRAD: assembly of de novo RADseq loci for phylogenetic analyses. Bioinformatics. 2014;30:1844–9.

723 9. Gautier M, Gharbi K, Cezard T, Foucaud J, Kerdelhué C, Pudlo P, et al. The effect of RAD allele dropout on the estimation of genetic variation within and between populations. Mol Ecol. 2013;22:3165–78.

726 10. Zhou Q, Zhang J, Bachtrog D, An N, Huang Q, Jarvis ED, et al. Complex evolutionary trajectories of sex chromosomes across bird taxa. Science. 2014;346:1246338.

728 11. Rutkowska J, Lagisz M, Nakagawa S. The long and the short of avian W chromosomes: no evidence for gradual W shortening. Biol Lett. 2012;8:636–8.

730 12. Gamble T, Zarkower D. Identification of sex-specific molecular markers using restriction site-associated DNA sequencing. Mol Ecol Resour. 2014;14:902–13.

732 13. Fowler BLS, Buonaccorsi VP. Genomic characterization of sex-identification markers in Sebastes carnatus and Sebastes chrysomelas rockfishes. Mol Ecol. 2016;25:2165–75.

734 14. Stovall WR, Taylor HR, Black M, Grosser S, Rutherford K, Gemmell NJ. Genetic sex assignment in wild populations using genotyping-by-sequencing data: A statistical threshold approach. Mol Ecol Resour. 2018;18:179–90.

737 15. Ellergren H. First gene on the avian W chromosome (CHD) provides a tag for universal sexing of non-ratite birds. Proc R Soc Lond B. 1996;263:1635–41.

739 16. Griffiths R, Korn RM. A CHD1 gene is Z chromosome linked in the chicken Gallus domesticus. Gene. 1997;197:225–9.
741 17. Malleret MM. Citogenética comparativa entre espécies de aves Suliformes: Sula dactylatra, Sula leucogaster (Sulidae) e Fregata magnificens (Fregatidae). Master Thesis. Universidade Federal de Minas Gerais; 2017.

744 18. Spriggs EL, Eaton DAR, Sweeney PW, Schlutius C, Edwards EJ, Donoghue MJ. Restriction-Site-Associated DNA Sequencing Reveals a Cryptic Viburnum Species on the North American Coastal Plain. Syst Biol. 2019;68:187-203.

747 19. Toews DPL, Brelsford A, Grossen C, Milá B, Irwin DE. Genomic variation across the Yellow-rumped Warbler species complex. The Auk. 2016;133:698-717.

749 20. White TA, Perkins SE, Heckel G, Searle JB. Adaptive evolution during an ongoing range expansion: the invasive bank vole (Myodes glareolus) in Ireland. Mol Ecol. 2013;22:2971-85.

752 21. De Donato M, Peters SO, Mitchell SE, Hussain T, Imumorin IG. Genotyping-by-Sequencing (GBS): A Novel, Efficient and Cost-Effective Genotyping Method for Cattle Using Next-Generation Sequencing. PLoS ONE. 2013;8:e62137.

755 22. Lu F, Lipka AE, Glaubitz J, Elshire R, Cherney JH, Casler MD, et al. Switchgrass Genomic Diversity, Ploidy, and Evolution: Novel Insights from a Network-Based SNP Discovery Protocol. PLoS Genet. 2013;9:e1003215.

758 23. Weir BS, Cockerham CC. Estimating F-statistics for the analysis of population structure. evolution. 1984;38:1358-70.

760 24. Dodds KG, McEwan JC, Brauning R, Anderson RM, van Stijn TC, Kristjánsson T, et al. Construction of relatedness matrices using genotyping-by-sequencing data. BMC Genomics. 2015;16. doi:10.1186/s12864-015-2252-3.

763 25. Alexander DH, Novembre J, Lange K. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 2009;19:1655-64.

765 26. Kahn NW, St. John J, Quinn TW. Chromosome-specific intron size differences in the avian CHD gene provide an efficient method for sex identification in birds. The Auk. 1998;115:1074-8.

768 27. Sambrook J, Russell DW. Molecular cloning: a laboratory manual. 3rd edition. CSHL Press; 2001.

770 28. Lee JC-I, Tsai L-C, Hwa P-Y, Chan C-L, Huang A, Chin S-C, et al. A novel strategy for avian species and gender identification using the CHD gene. Mol Cell Probes. 2010;24:27-31.

773 29. Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009;25:1754-60.

775 30. Rutkowska J, Lagisz M, Nakagawa S. The long and the short of avian W chromosomes: no evidence for gradual W shortening. Biol Lett. 2012;8:636-8.

777 31. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, Bender D, et al. PLINK: A Tool Set for Whole-Genome Association and Population-Based Linkage Analyses. Am J Hum Genet. 2007;81:559-75.

780

781
783 **Table 1.** Sample collection information and descriptive statistics of the NGS data.

<table><thead><tr><th rowspan="2">Species</th><th rowspan="2">Sampling locations<br>(code)</th><th colspan="6">Number of reads after filtering<br>samples</th></tr><tr><th>N1</th><th>N2</th><th>avg</th><th>med</th><th>min</th><th>max</th></tr></thead><tbody><tr><td><i>S. leucogaster</i></td><td>São Pedro & São Paulo<br>Archipelago (ASP),<br>Fernando de Noronha<br>Archipelago (NOR)</td><td>43</td><td>36</td><td>4,240,039</td><td>3,848,212</td><td>1,861,463</td><td>7,019,443</td></tr><tr><td><i>S. dactylatra</i></td><td>Fernando de Noronha<br>Archipelago (NOR),<br>Trindade & Martim Vaz<br>Archipelago (TRI)</td><td>53</td><td>46</td><td>4,643,619</td><td>4,206,670</td><td>2,257,164</td><td>7,813,045</td></tr><tr><td><i>M. octosetaceus</i></td><td>Serra da Canastra (CAN),<br>Jalapão (JAL), Alto<br>Paranaíba (ALP),<br>Chapada dos Veadeiros<br>(VEA)</td><td>37</td><td>31</td><td>5,794,931</td><td>5,421,581</td><td>1,043,704</td><td>12,997,656</td></tr></tbody></table>

784 N1: Number of collected samples; N2: Number of samples after filtering samples

785

786 **Table 2.** Statistics of patterns of observation of sequences among samples and isolation of a pattern likely to be linked to chromosome W, for the three species.

787

<table><thead><tr><th rowspan="2"></th><th rowspan="2"># of<br>sequences<br>considered¹</th><th rowspan="2"># of<br>observation<br>patterns</th><th rowspan="2">Avg<br>#sequences<br>per<br>pattern</th><th colspan="4">Top-scoring pattern</th></tr><tr><th>#sequences</th><th>#samples</th><th>Score</th><th>Signal<br>strength<br>score²</th></tr></thead><tbody><tr><td><i>S. leucogaster</i></td><td>431,377</td><td>199,922</td><td>2.16</td><td>1894</td><td>19</td><td>236.96</td><td>89.5%</td></tr><tr><td><i>S. dactylatra</i></td><td>485,102</td><td>265,876</td><td>1.82</td><td>1884</td><td>23</td><td>220.44</td><td>94.9%</td></tr><tr><td><i>M. octosetaceus</i></td><td>480,654</td><td>216,972</td><td>2.21</td><td>795</td><td>13</td><td>76.35</td><td>85.7%</td></tr></tbody></table>

788 ¹Were considered only reliable sequences with minimum coverage per sample equal to 3.

789 ²Computed as $1 - S_{top2}/S_{top1}$, thus bound by [0,1].

790

791
<table><caption>**Table 3.** Number of markers obtained by three different pipelines and under four different filtering rules leading to four selections (numbered SI to SIV), for the three species.</caption><thead><tr><th rowspan="2">Selection number and filters</th><th colspan="3"><i>S. leucogaster</i></th><th colspan="3"><i>S. dactylatra</i></th><th colspan="3"><i>M. octosetaceus</i></th></tr><tr><th>Stacks</th><th>iPyRAD</th><th>GBirdS</th><th>Stacks</th><th>iPyRAD</th><th>GBirdS</th><th>Stacks</th><th>iPyRAD</th><th>GBirdS</th></tr></thead><tbody><tr><td>SI: Unfiltered</td><td>64,261</td><td>44,648</td><td>55,815</td><td>86,718</td><td>66,502</td><td>72,738</td><td>47,111</td><td>33,981</td><td>43,903</td></tr><tr><td>SII: COVmed&gt;=6 &amp; COVsin&gt;=10</td><td>27,013</td><td>20,126</td><td>13,609</td><td>34,471</td><td>25,429</td><td>14,596</td><td>21,973</td><td>19,078</td><td>12,797</td></tr><tr><td>SIII: COVmed&gt;=6 &amp; COVsin&gt;=10 &amp; CR&gt;=0.7 &amp; MAF&gt;=0.02</td><td>18,865</td><td>11,248</td><td>11,481</td><td>24,389</td><td>14,844</td><td>12,360</td><td>15,932</td><td>10,693</td><td>11,898</td></tr><tr><td>SIV: COVmed&gt;=6 &amp; COVsin &gt;=10 &amp; CR&lt;=1 &amp; MAF&gt;=0.02</td><td>8,087</td><td>6,006</td><td>7,686</td><td>6,090</td><td>5,297</td><td>4,589</td><td>2,512</td><td>2,114</td><td>3,327</td></tr></tbody></table>

COVmed = allele median coverage (rules applies to both alleles); COVsin = singleton allele coverage; CR = sample call rate; MAF = minor allelic frequency

<table><caption>**Table 4.** Prediction errors for Admixture models under different number of subpopulations (K), for the three pipelines and the three species.</caption><thead><tr><th></th><th>K</th><th>Stacks</th><th>iPyRAD</th><th>Pipeline</th></tr></thead><tbody><tr><td rowspan="4"><i>S. leucogaster</i></td><td>2</td><td><strong>0,376</strong></td><td><strong>0,420</strong></td><td><strong>0,423</strong></td></tr><tr><td>3</td><td>0,392</td><td>0,439</td><td>0,448</td></tr><tr><td>4</td><td>0,523</td><td>0,573</td><td>0,565</td></tr><tr><td>5</td><td>0,544</td><td>0,737</td><td>0,777</td></tr><tr><td rowspan="4"><i>S. dactylatra</i></td><td>2</td><td><strong>0.450</strong></td><td><strong>0.467</strong></td><td><strong>0.484</strong></td></tr><tr><td>3</td><td>0.497</td><td>0.518</td><td>0.533</td></tr><tr><td>4</td><td>0.556</td><td>0.569</td><td>0.594</td></tr><tr><td>5</td><td>0.603</td><td>0.622</td><td>0.664</td></tr><tr><td rowspan="4"><i>M. octosetaceus</i></td><td>2</td><td><strong>0.562</strong></td><td>0.641</td><td>0.606</td></tr><tr><td>3</td><td>0.664</td><td><strong>0.639</strong></td><td><strong>0.589</strong></td></tr><tr><td>4</td><td>0.615</td><td>0.727</td><td>0.712</td></tr><tr><td>5</td><td>0.700</td><td>0.789</td><td>0.723</td></tr></tbody></table>

Lowest prediction errors per species and per pipeline are boldfaced.
<table><caption>**Table 5.** Fixation index *F*<sub>ST</sub>, computed between pairs of populations within each of the three species, averaged on the SIII selection of markers obtained by the three pipelines.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Pop. 1</th><th rowspan="2">Pop. 2</th><th colspan="3">Fixation index F<sub>ST</sub><br>(averaged on all SIII markers)</th></tr><tr><th>Stacks</th><th>iPyRAD</th><th>GBirdS</th></tr></thead><tbody><tr><td><em>S. leucogaster</em></td><td>ASP</td><td>NOR</td><td>22.8%</td><td>20.8%</td><td><strong>24.4%</strong></td></tr><tr><td><em>S. dactylatra</em></td><td>NOR</td><td>TRI</td><td>2.6%</td><td>2.2%</td><td><strong>2.8%</strong></td></tr><tr><td rowspan="6"><em>M. octosetaceus</em></td><td>CAN</td><td>JAL</td><td>15.6%</td><td>13.4%</td><td><strong>16.1%</strong></td></tr><tr><td>CAN</td><td>ALP</td><td>7.9%</td><td>6.7%</td><td><strong>8.1%</strong></td></tr><tr><td>CAN</td><td>VEA</td><td>10.5%</td><td>8.6%</td><td><strong>10.9%</strong></td></tr><tr><td>JAL</td><td>ALP</td><td>20.1%</td><td>18.0%</td><td><strong>20.5%</strong></td></tr><tr><td>JAL</td><td>VEA</td><td>18.1%</td><td>16.0%</td><td><strong>18.4%</strong></td></tr><tr><td>ALP</td><td>VEA</td><td>14.0%</td><td>12.2%</td><td><strong>14.1%</strong></td></tr></tbody></table>

Highest average *F*<sub>ST</sub> per comparison of populations are boldfaced.
<table><caption>**Table 6.** Comparisons between the KGD matrices obtained by three different pipelines, by correlations<fnref n="1" /> between the three first principal components (PC1, PC2 and PC3) and the diagonal (D) and off-diagonal (OD) elements of the KGD matrix, for the three species.</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="5">GBirdS</th><th colspan="5">Stacks</th></tr><tr><th>PC1</th><th>PC2</th><th>PC3</th><th>D</th><th>OD</th><th>PC1</th><th>PC2</th><th>PC3</th><th>D</th><th>OD</th></tr></thead><tbody><tr><td rowspan="10">S. leucogaster</td><td rowspan="5">Stacks</td><td>PC1</td><td style="background-color: #cccccc;">&gt;99.9%</td><td>0.4%</td><td>-0.2%</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>PC2</td><td>-0.4%</td><td style="background-color: #cccccc;">97.7%</td><td>9.0%</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>PC3</td><td>-0.1%</td><td>-15.7%</td><td>44.0%</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>D</td><td></td><td></td><td></td><td style="background-color: #cccccc;">99.6%</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>OD</td><td></td><td></td><td></td><td></td><td style="background-color: #cccccc;">99.9%</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="5">iPyRAD</td><td>PC1</td><td style="background-color: #cccccc;">&gt;99.9%</td><td>0.8%</td><td>-0.5%</td><td></td><td></td><td style="background-color: #cccccc;">&gt;99.9%</td><td>0.4%</td><td>-0.1%</td><td></td><td></td></tr><tr><td>PC2</td><td>-0.8%</td><td style="background-color: #cccccc;">97.7%</td><td>6.2%</td><td></td><td></td><td>-0.4%</td><td style="background-color: #cccccc;">99.3%</td><td>0.3%</td><td></td><td></td></tr><tr><td>PC3</td><td>0.0%</td><td>-15.6%</td><td>40.4%</td><td></td><td></td><td>0.0%</td><td>-2.3%</td><td style="background-color: #cccccc;">94.6%</td><td></td><td></td></tr><tr><td>D</td><td></td><td></td><td></td><td style="background-color: #cccccc;">99.4%</td><td></td><td></td><td></td><td></td><td style="background-color: #cccccc;">99.8%</td><td></td><td></td></tr><tr><td>OD</td><td></td><td></td><td></td><td></td><td style="background-color: #cccccc;">99.8%</td><td></td><td></td><td></td><td></td><td style="background-color: #cccccc;">99.9%</td><td></td></tr><tr><td rowspan="10">S. dactylatra</td><td rowspan="5">Stacks</td><td>PC1</td><td style="background-color: #cccccc;">99.5%</td><td>0.6%</td><td>0.2%</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>PC2</td><td>-0.3%</td><td>-24.0%</td><td>55.1%</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>PC3</td><td>5.7%</td><td style="background-color: #cccccc;">78.9%</td><td>45.6%</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>D</td><td></td><td></td><td></td><td style="background-color: #cccccc;">59.4%</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>OD</td><td></td><td></td><td></td><td></td><td style="background-color: #cccccc;">95.5%</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="5">iPyRAD</td><td>PC1</td><td style="background-color: #cccccc;">98.9%</td><td>1.6%</td><td>-1.9%</td><td></td><td></td><td style="background-color: #cccccc;">99.5%</td><td>-3.7%</td><td>2.4%</td><td></td><td></td></tr><tr><td>PC2</td><td>5.1%</td><td>-12.8%</td><td>26.7%</td><td></td><td></td><td>5.4%</td><td style="background-color: #cccccc;">89.3%</td><td>-11.1%</td><td></td><td></td></tr><tr><td>PC3</td><td>5.8%</td><td style="background-color: #cccccc;">72.3%</td><td>50.8%</td><td></td><td></td><td>3.5%</td><td>9.6%</td><td style="background-color: #cccccc;">97.1%</td><td></td><td></td></tr><tr><td>D</td><td></td><td></td><td></td><td style="background-color: #cccccc;">38.3%</td><td></td><td></td><td></td><td></td><td style="background-color: #cccccc;">62.0%</td><td></td><td></td></tr><tr><td>OD</td><td></td><td></td><td></td><td></td><td style="background-color: #cccccc;">87.9%</td><td></td><td></td><td></td><td></td><td style="background-color: #cccccc;">95.1%</td><td></td></tr><tr><td rowspan="10">M. octosetaceus</td><td rowspan="5">Stacks</td><td>PC1</td><td style="background-color: #cccccc;">99.9%</td><td>-0.4%</td><td>-1.0%</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>PC2</td><td>1.0%</td><td style="background-color: #cccccc;">99.9%</td><td>-0.4%</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>PC3</td><td>1.5%</td><td>0.4%</td><td style="background-color: #cccccc;">99.5%</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>D</td><td></td><td></td><td></td><td style="background-color: #cccccc;">90.2%</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>OD</td><td></td><td></td><td></td><td></td><td style="background-color: #cccccc;">99.6%</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="5">iPyRAD</td><td>PC1</td><td style="background-color: #cccccc;">99.9%</td><td>-2.0%</td><td>-2.3%</td><td></td><td></td><td style="background-color: #cccccc;">99.9%</td><td>-1.1%</td><td>-1.4%</td><td></td><td></td></tr><tr><td>PC2</td><td>2.6%</td><td style="background-color: #cccccc;">99.4%</td><td>2.7%</td><td></td><td></td><td>1.8%</td><td style="background-color: #cccccc;">99.7%</td><td>2.5%</td><td></td><td></td></tr><tr><td>PC3</td><td>3.5%</td><td>-4.6%</td><td style="background-color: #cccccc;">98.5%</td><td></td><td></td><td>1.8%</td><td>-4.5%</td><td style="background-color: #cccccc;">97.2%</td><td></td><td></td></tr><tr><td>D</td><td></td><td></td><td></td><td style="background-color: #cccccc;">77.1%</td><td></td><td></td><td></td><td></td><td style="background-color: #cccccc;">96.3%</td><td></td><td></td></tr><tr><td>OD</td><td></td><td></td><td></td><td></td><td style="background-color: #cccccc;">98.9%</td><td></td><td></td><td></td><td></td><td style="background-color: #cccccc;">99.4%</td><td></td></tr></tbody></table>

<footnote n="1">¹The more distant value from 0, the darker the cell background.</footnote>
<table><caption>**Table 7.** Match between filtered markers and the pattern of observation of sequences that was found likely to be linked to chromosome W by GBirdS, for the three pipelines and the three species.</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="3">S. leucogaster</th><th colspan="3">S. dactylatra</th><th colspan="3">M. octosetaceus</th></tr><tr><th>Stacks</th><th>iPyRAD</th><th>GBirdS</th><th>Stacks</th><th>iPyRAD</th><th>GBirdS</th><th>Stacks</th><th>iPyRAD</th><th>GBirdS</th></tr></thead><tbody><tr><td rowspan="2"># of markers matching with pattern</td><td>SIII</td><td>178</td><td>69</td><td>0</td><td>142</td><td>106</td><td>0</td><td>10</td><td>107</td><td>0</td></tr><tr><td>SIV</td><td>87</td><td>8</td><td>0</td><td>27</td><td>0</td><td>0</td><td>5</td><td>49</td><td>0</td></tr><tr><td rowspan="5">Correlations<fnref n="1" /> between PCs and pattern</td><td>PC1</td><td>21.0%</td><td>21.0%</td><td>-28.1%</td><td>1.1%</td><td>5.7%</td><td>1.0%</td><td>6.3%</td><td>6.5%</td><td>7.0%</td></tr><tr><td>PC2</td><td>30.8%</td><td>29.0%</td><td>-11.0%</td><td>-80.1%</td><td>-96.5%</td><td>7.0%</td><td>18.4%</td><td>21.9%</td><td>17.3%</td></tr><tr><td>PC3</td><td>59.1%</td><td>44.6%</td><td>-4.2%</td><td>12.3%</td><td>3.0%</td><td>-20.0%</td><td>13.0%</td><td>19.7%</td><td>14.5%</td></tr><tr><td>PC4</td><td>25.3%</td><td>42.4%</td><td>-32.1%</td><td>11.9%</td><td>-7.1%</td><td>9.2%</td><td>16.3%</td><td>15.5%</td><td>13.0%</td></tr><tr><td>PC5</td><td>-22.6%</td><td>4.4%</td><td>-9.4%</td><td>-25.8%</td><td>0.7%</td><td>-32.1%</td><td>-17.9%</td><td>-24.4%</td><td>-17.1%</td></tr></tbody></table>

<table><caption>**Table 8.** Usage of computational resources (computing time, memory peak) and requested number of CPUs, for each pipeline and each species.</caption><thead><tr><th></th><th></th><th>Number of CPUs</th><th>Sula leucogaster</th><th>Sula dactylatra</th><th>Mergus octosetaceus</th></tr></thead><tbody><tr><td rowspan="3">Computing time (minutes)</td><td>Stacks</td><td>16</td><td>425</td><td>645</td><td>528</td></tr><tr><td>iPyRAD</td><td>16</td><td>1027</td><td>1404</td><td>1176</td></tr><tr><td>GBirdS</td><td>1</td><td>48</td><td>53</td><td>26</td></tr><tr><td rowspan="3">Memory peak (Gbytes)</td><td>Stacks</td><td>16</td><td>36.7</td><td>43.7</td><td>51.3</td></tr><tr><td>iPyRAD</td><td>16</td><td>13.2</td><td>16.0</td><td>11.3</td></tr><tr><td>GBirdS</td><td>1</td><td>5.8</td><td>8.1</td><td>6.5</td></tr></tbody></table>

<footnote n="1">¹The more distant the value from 0, the darker the cell background.</footnote>
823

## 11.Figures

824

**Figure 1.** Graphical overview of the pipeline. (a) All reads produced by sequencing are

825

demultiplexed and checked for base-pair quality (Phred score), barcodes and cutting-site

826

remnant are removed. The reads are then uniquely sorted into unique sequences, which are

827

then checked for quality (removal of outliers for GC-content and over-coverage) and divided

828

between singletons (observed only in one individual) and shared (observed in more than one

829

individual) libraries. Meanwhile, samples are checked for quality (total number of reads,

830

contamination). (b) Patterns of observation of sequences from the shared library are uniquely

831

sorted. Their frequencies and the probability that they reflect the gender distribution are

832

computed. The library is finally split between final sequences used for genotype calling and

833

sequences likely to be linked to chromosome W. (c) All-versus-all comparisons of remaining

834

sequences in order to search for small differences (SNPs, short indels) between these

835

sequences. Genotypes are called from reciprocally-matching sequences.
(a) From the FASTQ file produced by sequencing

<figure><img src="image_2.png" /></figure>

(b) From the library of shared reads

<figure><img src="image_3.png" /></figure>

(c) From the final library of reads

<figure><img src="image_4.png" /></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Plotting the *Sula dactylatra* samples against the two first principal components, obtained by three different pipelines. The first row features the first principal components of the KGD matrix (selection SIII) while the second row features the first principal components of the genotype matrix of markers observed in all samples (selection SIV). The three columns correspond to the three pipelines: (a) Stacks, (b) iPyRAD and (c) GBirdS. For any sample, the marker shape tells if the reads likely to linked to chromosome W by GBirdS were observed (triangle) or not (circle) in that sample and the marker color indicates the geographic origin (yellow for TRI and blue for NOR).</figcaption></figure>

838

839

840

841

842
844

Figure 3. Distribution of the proportion ($p_w$) of sequences matching to the pattern found to be linked to chromosome W, for samples where the pattern was observed and for each of the three species.

845

<figure><img src="image_6.png" /></figure>

846
847

**Figure 4.** Corrected coverage in sequences matching to ten avian chromosomes Z, for all three species. Samples in which sequences likely

848

to be linked to chromosome W were observed are plotted with plain-blue circles and samples in which these sequences were not observed are

849

plotted with open-red circles. Samples sexed by on-field/molecular sexing are annotated with their sex.

<figure><img src="image_7.png" /></figure>