dataesr/scientific-paragraphs-categorization huggingface dataset:
org/2013/11/04/oligotyping-best-practices/). We also removed the SNVs identi ed as outliers by anvi'o (cov_outlier_in_contig = 1, cov_outlier_in_split = 1). We used the "anvi-script-snvs-to-interactive" program to generate a visualizable SNVs database from the intra-MAG ltered variability table using the "anvi-interactive" program. We then focused our SNVs analyses on a set of SSG (Shared Single-copy Genes, that is genes that occur only once in each genome of the genome set) as de ned by the pangenomic analyses of our MAGs and Wolbachia reference genomes described in the section below. We added the SSG assignation in the visualizable SNVs database using the "anvi-import-misc-data" anvi'o program. Finally, we edited the SNVs gure using the "anvi-interactive" program.
Pangenomics We performed pangenomic analysis for i) all the Wolbachia MAGs, ii) each Wolbachia MAGs and 3 selected Wolbachia reference genomes: wPipPel isolated from Culex quinquefasciatus (NCBI Accession ID NC_010981. 1), wPipMol isolated from Culex molestus (NCBI Accession ID NZ_CTEH00000000. 1) and wPipJHB isolated from Culex quinquefasciatus (NCBI Accession ID NZ_ABZA00000000. 1). To do that, we selected and downloaded the fasta les of the three selected Wolbachia reference genomes and reformatted them using the anvi'o "anvi-script-reformat-fasta" program. We then generated a new contigs database from the reformatted fasta les using the "anvi-gen-contigs-database" program. We identi ed HMM hits using the "anvi-run-hmms" program to assign functions but also taxonomy to genes using the "anvi-run-ncbi-cogs" and "anvi-run-scg-taxonomy" programs, respectively.
We created an external genome database including all these Wolbachia reference contigs databases. We created an internal genome database including all the Wolbachia MAGs contigs database (M11, O11, O03, O07, O12) stored in the pro le and contigs databases generated during the snakemake work ow with references mode. We then generated a genome storage database from both external and internal genome databases using the "anvi-gen-genomes-storage" program. We ran the "anvi-pan-genome" program (using '--use-ncbi-blast' and '--mcl-in ation 10' ags and the "genome-name" ag) to identify gene clusters (GC) and compute the pangenome for i) all the Wolbachia MAGs to study the gene similarity inter-MAGs, ii) each Wolbachia MAG and all the 3 Wolbachia reference genomes to study the similarity between MAGs and Wolbachia reference genomes.
To note, only highly similar genes are added in a GC during the anvi'o pangenomic work ow with almost no chance for two highly similar genes to end up in distinct GC (which was essential for our study when we compared SNVs occurring in SSGCs). We nally used the "anvi-display-pan" program to display the pangenomes and visualize the distribution of gene clusters across genomes from our different genome sets. And we identi ed SSGs using a custom R script among the GC identi ed during the pangenomic analyses performed on the two genome sets detailed above.
Metapangenomics We used a custom R script (based on https://merenlab. org/data/wolbachia-plasmid/#recoveringcoverage-values-for-gene-clusters-of-the-wolbachia-pangenome-in-c-pipiens-metagenomes) to recover coverage values of M11 and O11 MAGs and compute the mean of M11 and O11 metagenome coverage values by gene cluster. We used the available results of blastn that identi ed wPipPel genes that match to WO prophage regions ("WO_in_wPip_best_hit. txt" le from https://merenlab. org/data/wolbachia-plasmid/#identifying-genesthat-correspond-to-wo-prophages detailed in  [32] ) together with the same custom R script to detect "phage like" gene clusters. We imported coverage values and WO prophage assignation to the pangenome database using the "anvi-import-misc-data" anvi'o program to build the metapangenome.
We nally ran and edited the metapangenome using the "anvi-display-pangenome" program. MLST, wsp, cidA and cidB hits We used blastn to identify MLST (MultiLocus Sequence Typing -gatB, coxA, hcpA, ftsZ, fbpA) and wsp genes from the PubMSLT [48] (Public databases for molecular typing and microbial genome diversity) in our Wolbachia MAGs. Similarly, we identi ed best hits using blastn for cidA (NCBI Accession ID from MF444963 to MF444981 for the 18 cidA variants) and cidB (NCBI Accession ID from MF444982 to MF444996 for the 14 cidB variants) genes from Bonneau et al. 2018  [49]. We nally used a custom R script to tag SNVs occurring in these hits and added them to the SNVs database using the "anvi-importmisc-data" anvi'o program.
Results Reconstruction of Wolbachia metagenome-assembled genomes (MAGs) in one midgut and four ovaries of Culex pipiens individuals Our quality ltering of raw reads sequenced from midgut and ovaries resulted in 94,024,472 and 75,040,983 million paired-end reads on average, respectively (Supplementary Table  1 ). Individual sample metagenomic assembly generated 135k and 214k contigs > 1kb (mean numbers) with 60% and 47% of recruited ltered reads for each separated organ (Supplementary Table  1 ). To estimate the extent of eukaryotic contamination in our metagenomes, we used phyloFlash  [34]  that annotates short reads based on the SILVA rRNA databases  [35]. Data suggested that the vast majority of our reads (over 99% for each sample) originated from the host genome, especially in midgut metagenomes (Supplementary Fig.
These observations suggested the occurrence of synonymous mutations without any amino acid change in this cid gene. We did not identify SNVs in single-copy cid genes in O03, O07 and O12 (Supplementary Table  13 , sheets 3 & 5 and Supplementary Figs. 5  & 6 ). To note, we reported some SNVs in multi-copy cid genes in O07 (Supplementary Note 3). In addition, we also did not detect any SNV occurring in genes annotated to MLST or wsp genes (from PubMLST) in our full dataset (Supplementary Table  12 , Supplementary Table  13 , Supplementary Figs. 2 3 4 5 6 ), in agreement with previous studies. SNVs found outside of SSG are reported in Supplementary Note 3 (also see Supplementary Table  13 , Supplementary Fig. 4 ). b. Wolbachia MAGs M11, O03, O07, and O12. From inner to outer and for each SNV, the layers correspond to i) the departure from the consensus (max value set to 1), ii) the Shannon entropy (max value set to 1), iii) the coverage (max value set to 2500X for O11, O03 and 800X for M11, O07, O12), iv) the position in codon (from 1st to 3rd if the SNV is detected in a coding gene and 0th if it is not), v) the competing nucleotides (the two most abundant nucleotides occurring in a SNV), vi) the presence / absence of SNV in a coding gene call, vii) the presence in a GC shared by the MAG and at least another one, viii) the presence in a gene shared by the MAG and at least 1 Wolbachia reference genomes.
Acknowledgments We are thankful to  Valentin Loux ,  Cédric Midoux ,  Véronique Martin  for their support on the MIGALE highcomputer during the data analysis,  Corinne Cruaud  from the Genoscope for her assistance with sequencing, and to  Blandine Trouche  for helpful discussions about SNVs analyses. We also thank  EID Mediterranée  and  Ignace Rakotoarivony  for mosquito sampling and dissections, respectively.
Data availability The raw sequencing data for the shotgun metagenomes of midgut and ovary samples are available in the European Nucleotide Archive via accession number PRJEB56379 and PRJEB26028, respectively. In addition, we made available the merged anvi'o pro les for the midgut and ovary metagenomes (10. 5281/zenodo. 7183277), the FASTA les for the ve Wolbachia Metagenome-Assembled Genomes (10. 5281/zenodo. 7183304) and the anvi'o merged pro le databases for the Wolbachia MAGs used for the SNVs and metapangenomic analyses (10. 5281/zenodo. 7183324). Availability of data and data generated or analyzed during this study are included in this published article (and its supplementary information les).
Negative controls (PCR mix only) and positive controls (RNA from viral strains or DNA from bacterial strains) were included in each run. Positive results were defined as the amplification with a cycle threshold (CT) value ≤35. A threshold value of 35 was used in each experimental run and we calculated the RFU cut-off value recommend by CFX Manager Software Version 3. 1 (Bio-Rad) in order to verify the positive cases. Results were considered to be positive when the cycle threshold value of real-time PCR was greater than the cut-off value. All quantitative real-time PCR tests were performed using a C1000 TouchÔ Thermal Cycle (Bio-Rad, Hercules, CA, USA).
Genome assemblies The male PacBio reads were assembled using 2 different genome assemblers: Flye (v. 2. 7. 1,  Kolmogorov et al. 2019)  and Miniasm (0. 3-r179, minimap2 2. 18-r1028-dirty was used for mapping and the consensus was generated using Racon v1. 4. 22)  (Li 2016; Vaser et al. 2017). The Flye assembly was polished using male A. sinica short genomic reads (trimmed with the Trimmomatic package,  Bolger et al. 2014) , and the Miniasm assembly was polished using the same male short reads using the wtpoa-cns tool from wtdbg2 (version 2. 5, Ruan and Li 2020). The 2 assemblies were then merged using quickmerge (version 0. 3,  Chakraborty et al. 2016)  with the Miniasm assembly as the query and the Flye assembly as the reference. The resulting assembly was purged based on the male pacbio read depth to remove duplicates and contig overlaps using the purge_dups program (version 1.
2. 5,  Guan et al. 2020). To scaffold the assembly into pseudo-chromosomes, PCR duplicates were first removed from the Hi-C data using the clumpify. sh script from the BBMap package  (Bushnell 2014) , and the Hi-C reads were then mapped to the genome assembly using the Arima mapping pipeline with MAPQ 5 (Arima Genomics 2021) and then scaffolded using the YaHS tool (pre-release of version 1. 1,  Zhou et al. 2022). The contact maps were visualized and manually edited on Juicebox (version 1. 11. 08, Robinson et al. 2018)  to generate the final chromosome-level assembly. The female A. sinica genome was assembled from female PacBio reads using Flye (version 2. 7. 1), and it was not polished to avoid collapsing the Z and the W scaffolds. To identify putative W scaffolds, short genomic reads from 2 A.
sinica males and 2 females were mapped to the female assembly using Bowtie2  (Langmead and Salzberg 2012). We then counted how many male and female reads mapped to each scaffold, after filtering for alignments without mismatches (by selecting only mapped reads with the CIGAR string "NM:i:0"). The female-specific k-mers inferred to obtain W-specific transcripts (section Identification of candidate W-genes with k-mer analysis below) were similarly mapped to each scaffold with Bowtie2 and counted. Scaffolds which had more than 5 female-specific k-mers, and more perfect matches in females than in males [male/(male þ female) 0. 3] were considered candidates W-derived scaffolds. The A. sp. Kazakhstan genome was assembled from 2 male short-read libraries with  Megahit (v1. 1. 4, Li et al. 2015)  and then scaffolded using SOAPdenovo-fusion (SOAPdenovo2 version 2.
Kazakhstan genomes to the new A. sinica assembly We aligned the A. sinica published transcriptome  (Huylmans et al. 2021)  to both the A. franciscana and to the A. sp. Kazakhstan genomic scaffolds using blat  (Standalone BLAT v. 36x2, Kent 2002). For each transcript, we kept only the mapping location with the highest score in each genome (using the customized script 1-besthitblat. pl). When multiple transcripts overlapped by more than 20 bp on the genome, only the transcript with the highest mapping score was kept (2-redremov_blat_v2. pl). We then used the location of the transcripts on the A. sinica genome to infer the location of the A. franciscana and A. sp. Kazakhstan scaffolds based on the transcripts they carried (AssignScaffoldLocation. pl). This script uses a majority rule to assign each scaffold to a chromosome, and then the mean location of genes on that scaffold to infer its final coordinate on the chromosome.
All scripts are available on our git page. F ST between male and female populations RNA-seq reads from 10 pooled A. sinica males and 10 pooled A. sinica females (from  Huylmans et al. 2021) , sampled from head, thorax, and gonads, were trimmed with Trimmomatic  (Bolger et al. 2014 ) and pooled by sex, and mapped separately to the male A. sinica reference genome using STAR  (Dobin et al. 2013 ) with default parameters. The resulting alignments with MAPQ score lower than 20 were filtered out and the remaining alignments were sorted using samtools view and sort functions  (Li et al. 2009). Next, a pileup file of male and female alignments was produced using the samtools-mpileup function. Finally, we used scripts from the Popoolation2 package  (Kofler et al. 2011)  to calculate F ST. The mpileup file was reformatted with the Popoolation2 mpileu-p2sync.
pl script, and the resulting synchronized file was used as an input for fst-sliding. pl script. F ST between male and female populations was calculated for windows of 1,000 nucleotides, using the fst-sliding. pl script with following options -suppress-noninformative -min-count 3 -min-coverage 10 -max-coverage 200min-covered-fraction 0. 5 -window-size 1,000 -step-size 1,000pool-size 10. We applied the same pipeline to estimate male:female F ST using head and gonad RNA-seq samples obtained from 10 males and 10 females of A. franciscana (from  Huylmans et al. 2019). The resulting F ST values were plotted based on the inferred location of the genomic scaffolds along the A. franciscana chromosomes (section Mapping of the A. franciscana and A. sp. Kazakhstan genomes to the new A. sinica assembly).
Identification of candidate W-genes with k-mer analysis We used a k-mer-based subtraction approach  (Elkrewi et al. 2021)  based on the tools included in the BBMap package  (Bushnell 2014 ) on male and female genomic and RNA-seq data from A. franciscana and A. sinica. The pipeline was applied to each species separately. In A. sinica, 2 male and 2 female DNA libraries and 2 whole body RNA-seq replicates for each sex were used (Supplementary Tables  1  and 2 ). In A. franciscana, the analysis was performed using 1 male and 1 female DNA libraries and pools of 2 RNA-seq replicates of heads and gonads for each sex, along with 1 whole body male and female RNA-seq libraries (SRR14598203 and SRR14598204). First, the shared 31-mers between the female DNA and RNA libraries were identified, and then any k-mers matching male libraries were removed.
Female RNA-seq reads containing these female-specific k-mers [with minimum k-mer fraction of 0. 6 (mkf ¼ 0. 6)] were retrieved and assembled using Trinity  (Grabherr et al. 2011) , and the perl script from the Trinity package (get_longest_isoform_seq_per_trinity_gene. pl) was used to keep only the longest isoform. The male and female genomic libraries were mapped to the assembled transcripts using Bowtie2  (Langmead and Salzberg 2012, p. 2), and candidates with a sum of female perfect matches 8 and a ratio of sum-of-females/ sum-of-males 2 were removed. The final set consisted of 402 transcripts in A. franciscana and 319 in A. sinica.
Transcriptome assemblies and expression analysis The A. sinica male transcriptome was assembled from 2 replicates of male whole body RNA-seq data  (Huylmans et al. 2021)  using Trinity  (Grabherr et al. 2011 ) in 2 different modes: denovo and genome-guided. The 2 assemblies were concatenated and then filtered using the tr2aacds. pl script from EvidentialGene  (Gilbert 2019). The transcriptome was annotated with the Pannzer annotation server (To ¨ro ¨nen et al. 2018), and mapped to the A. sinica genome using the same procedure as described in the section Mapping of W candidates to the A. sinica genome. For the expression analysis, only the first isoform was kept for each gene, and only transcripts longer than 500 bp were used in the analysis. The RNA-seq reads from the A. sinica heads, gonads, and thoraces of males and females  (Huylmans et al.
2021)  were mapped to the curated transcriptome and gene expression levels (in Transcripts per million, TPM) were obtained using Kallisto (version 0. 46. 2, Bray et al. 2016). Quantile normalization was done using NormalyzerDE  (Willforss et al. 2019). Two different A. franciscana de novo transcriptome assemblies were made using Trinity. The first using pooled RNA-seq reads from male heads and testes [2 replicates each  (Huylmans et al. 2019) ], and the second using the published whole-body male RNA-seq library (SRR14598203, Jo, Lee, Choi, Kim, Lee, and Park 2021). The 2 assemblies were concatenated and then filtered using the tr2aacds. pl script from EvidentialGene, and mapped to the A. sinica genome using the same procedure as described in section Mapping of W candidates to the A. sinica genome.
Phylogenetic trees The W candidates of A. sinica and A. franciscana were mapped reciprocally to each other using pblat (v. 36x2 with default parameters, Wang and Kong 2019), and reciprocal best hits were considered shared candidates. The W candidates of the 2 species were further mapped to their respective uncollapsed male transcriptome assemblies (see previous section) with pblat (Wang and Kong 2019) with a translated query and database, and a minimum match score of 50. The transcripts with the highest mapping score to the W candidates were used as the putative Z homologs in their respective species. The Branchinecta lindahli transcriptome  (Schwentner et al. 2018 ) was downloaded from the Crustacean Phylogeny dataset on Harvard Dataverse (https://doi. org/10. 7910/DVN/SM7DIU). B. lindahli homologs of shared W-candidates were obtained by mapping the putative Z homologs of both species to the B.
lindahli transcriptome using pblat (-minScore ¼ 50 -t ¼ dnax -q ¼ dnax) and retrieving the transcript with highest alignments score (using the customized script 2-besthitblat. pl). A transcript was considered a homolog if it mapped to at least one of the putative Z homologs of the 2 species, and when the 2 Z homologs mapped to different outgroup sequences, both outgroup sequences were retrieved and used to make 2 different alignments. The shared W candidates of A. sinica and A. franciscana, their Z homologs, and the outgroup sequences were aligned using MAFFT (version v7. 487, with the options "mafft -adjustdirection INPUT > OUTPUT,"  Katoh et al. 2002). The resulting alignments were fed to phylogeny. fr (Dereeper et al. 2008) , where the alignment was curated using Gblocks  (Talavera and Castresana 2007) , and the phylogenetic tree was constructed using PhyML  (Guindon et al.
e. had the same F2 father and A. sp. Kazakhstan mothers) as the asexual females, but were otherwise selected randomly. Illumina short-read sequencing was then performed as described in the section DNA short-and long-read sequencing. Analysis of backcross between the Aibi Lake rare male and A. sp. Kazakhstan females We sequenced 5 asexual females and 10 putatively sexual females from the F2 generation. This resulted in an average of 101 million reads per asexual female and 50 million reads per putatively sexual female. We first used SEQTK v1. 2 (https://github. com/lh3/seqtk) to randomly select a subset of reads from each asexual sample to match the number of reads of the smallest sample (to avoid biasing allele estimates toward high-coverage individuals). We removed adaptors and trimmed reads using Trimmomatic v0.
39  (Bolger et al. 2014). We then aligned the resulting paired-end reads to the genome using Bowtie2 v2. 4. 4  (Langmead and Salzberg 2012). SAM files were converted to BAM files and sorted in Samtools v. 1. 13 (Li et al. 2009). For our pooled analyses, we merged BAM files into a pooled asexual BAM file and a pooled putatively sexual BAM, and created a mpileup file in Samtools v. 1. 13. We then used Popoolation2  (Kofler et al. 2011)  to call F ST for both individual SNPs and in 1-kb windows. We used F ST computed for 1-kb windows to visualize F ST across the genome in a Manhattan plot in the R package qqman  (Turner 2018). We computed rolling medians in sliding windows of 101 consecutive SNPs on each linkage group using the rollmedian function from the package zoo  (Zeileis and Grothendieck 2005)  in R v.
14  (Li et al. 2009) , and used to call SNPs with the mpileup function of BCFtools v1. 14  (Li 2011 ). The VCF file was filtered with VCFtools v0. 1. 17 (Danecek et al. 2011)  for minimum and maximum depth (4 and 50), minimum quality score (30), and minimum minor allele frequency (0. 1). Only SNPs for which the 2 A. sp. Kazakhstan had a genotype of 0/0, and the 2 A. parthenogenetica individuals 1/1, were kept for further analyses. We computed F ST between the F2 asexual and control females using the function -weir-fst-pop in VCFtools for 10-kb windows. We then inferred ancestry of each genomic scaffold in every sample (i. e. whether they were homozygous for the A. sp. Kazakhstan haplotype, or carried a copy of the A. parthenogenetica haplotype as well) using the customized script Chromopaint.
In total, only 17 scaffolds with an assigned location on the A. sinica genome show ancestry patterns consistent with an asexuality locus (i. e. they show evidence of A. parthenogenetica ancestry in all asexual females, but not in all control females). Eleven are on the Z chromosome (vs 1 expected, P ¼ 1. 3eÀ20 with a chi-square test) and correspond to the region of high F ST , providing further support for a role of the Z chromosome in the transmission of asexuality. None of the other minor peaks of F ST are in regions with ancestry patterns consistent with asexuality loci (Supplementary Fig. 14 ), although chromosome 16 contains 3 such loci (vs 0. 9 expected, P < 0. 01 with a chi-square test). The annotation of genes located in the Z-linked candidate locus did not yield any obvious candidates (the annotation for all transcripts is provided as a Supplementary Dataset).
Data availability All genomic reads generated for this study are available at the NCBI short reads archive under Bioproject number PRJNA848277. The pipelines used to analyze the data are at https://git. ist. ac. at/ bvicoso/zsexasex2021, and important processed data files such as the new A. sinica genome assembly are provided at https://doi. org/10. 15479/AT:ISTA:11653. Supplemental material is available at GENETICS online.
Smartphone app Assessments were conducted using the MeMind app, a smartphone/web-based platform extensively described elsewhere  (Arenas-Castañeda et al. , 2020; Barrigón, Berrouiguet et al. , 2017; Berrouiguet et al. , 2015 Berrouiguet et al. , , 2017;; Jung et al. , 2019; Lopez-Castroman et al. , 2015) , which features a SMART-SCREEN protocol. This protocol is an innovative tool for mental health screening that has been designed to explore and understand the distribution of suicidal behaviours in the general population  (Arenas-Castañeda et al. , 2020). It detects and refer people at risk to current mental health facilities in the State of Yucatan, as it provides feedback about the assessment to improve help-seeking behaviour. Of note, as MeMind databases do not contain personal identification data on any individual, all data were collected anonymously using encrypted security protocols.
The use of structural equations enabled us to identify relationships and even formulate hypotheses concerning the directionality of the possible causal relationships between depression, anxiety, well-being, and suicidal behaviour. The PLS-SEM presented here consists of exploratory models aimed at depicting the overarching causal structure underlying these latent variables. These preliminary models will serve to construct and, ultimately, refine a proper prediction model for screening purposes. PLS-SEM holds great potential in this concern over traditional covariance-based SEM  (Hair et al. , 2016). Moreover, our data were assessed through validated built-in questionnaires on a smartphone other factors like academic stress, substance use, high-risk sports/hobbies and certain contextual factors not assessed in this study could also have influenced this finding.