DOI: 10.5061/dryad.hx3ffbgkt
title: The X chromosome of insects likely predates the origin of Class Insecta

description: Sex chromosomes have evolved independently multiple times, but why some
 are conserved for more than 100 million years whereas others turnover
 rapidly remains an open question. Here, we examine the homology of sex
 chromosomes across nine orders of insects, plus the outgroup springtails.
 We find that the X chromosome is likely homologous across insects and
 springtails; the only exception is in the Lepidoptera, which has lost the
 X and now has a ZZ/ZW sex chromosome system. These results suggest the
 ancestral insect X chromosome has persisted for more than 450 million
 years – the oldest known sex chromosome to date. Further, we propose that
 the shrinking of gene content of the Dipteran X chromosome has allowed for
 a burst of sex-chromosome turnover that is absent from other speciose
 insect orders.
description: This analysis uses publicaly available data. It has been
 processed through a series of custom scripts (included) and bioinformatics
 software programs.
description: Data files and pipeline for: "The X chromosome of insects likely
 predates Class Insecta" In Review: Evolution, 2023
 ######################################################################################### All data files are publicly available. Table S1 from supplement with this information: Isoperla grammatica (Plecoptera): genome: GCA\_945910005.1 (link: https://www.ncbi.nlm.nih.gov/assembly/GCA\_945910005.1/) RNA: ERR9881692 (link: https://www.ncbi.nlm.nih.gov/sra/?term=ERR9881692) Schistocerca americana (Orthoptera) genome: GCF\_021461395.2 (link: https://www.ncbi.nlm.nih.gov/assembly/GCF\_021461395.2/) Bombyx mori (Lepidoptera) genome: GCF\_014905235.1 (link: https://www.ncbi.nlm.nih.gov/assembly/GCF\_014905235.1/) Lucilia cuprina (Diptera) genome: GCF\_022045245.1 (link: https://www.ncbi.nlm.nih.gov/assembly/GCF\_022045245.1/) Coccinella septempunctata (Coleoptera) genome: GCF\_907165205.1 (link: https://www.ncbi.nlm.nih.gov/assembly/GCF\_907165205.1/) Chrysoperla carnea (Neuroptera) genome: GCF\_905475395.1 (link: https://www.ncbi.nlm.nih.gov/assembly/GCF\_905475395.1/) Acyrthosiphon pisum (Hemiptera) genome: GCF\_005508785.2 (link: https://www.ncbi.nlm.nih.gov/assembly/GCF\_005508785.2/) Sinella curviseta (Collembola) genome: GCA\_004115045.3 (link: https://www.ncbi.nlm.nih.gov/assembly/GCA\_004115045.3/) RNA: SRR7948082 (link: https://www.ncbi.nlm.nih.gov/sra/?term=SRR7948082) DNA: SRR7948081 (link: https://www.ncbi.nlm.nih.gov/sra/?term=SRR7948081) Allacma fusca (Collembola) genome: GCA\_947179485.1 (link: https://www.ncbi.nlm.nih.gov/assembly/GCA\_947179485.1) RNA: SRR1653171 (link: https://www.ncbi.nlm.nih.gov/sra/?term=SRR1653171) Ischnura elegans (Odonata) genome: GCF\_921293095.1 (link: https://www.ncbi.nlm.nih.gov/assembly/GCF\_921293095.1/) ######################################################################################### ######################################################################################### Included example dataset: Damsel.gff (gff for damselfly produced through Genespace (see below) - used as input in various parts of pipeline column info: Damsel.chr - chromosome location start - gene start end - gene end Damsel - gene name strand - strand (+/-) Damsel.ord - gene order ######################################################################################### <br> ######################################################################################### Annotating genomes for I. grammatica, S. curviseta, A. fusca 1. Downloading Data from SRA module load SRA-Toolkit/2.11.2 fasterq-dump --threads 4 --split-files 2. FastQC on Data from SRA fastqc -f fastq 3. Trimming Data from SRA module load trimmomatic/0.39 java -jar trimmomatic-0.39.jar PE sample.fastq sample.fastq sample\_1.paired.fq sample\_1.unpaired.fq sample\_2.paired.fq sample\_2.unpaired.fq ILLUMINACLIP:TruSeq2-PE.fa:2:30:10:2:keepBothReads LEADING:3 TRAILING:3 MINLEN:36 4. Transcriptome Assembly module load hisat2 module load tophat module load cufflinks module load bowtie2/2.4.4 module load python/2.7 module load samtools module load hisat2 module load stringtie #Build Genome Index in HISAT2 hisat2-build genome\_index #Aligning reads to genome using HISAT2 hisat2 --phred33 -p 50 --novel-splicesite-outfile hisat2/splicesite.txt -S hisat2/accepted\_hits.sam -x genome\_index -1 sample\_1.paired.fq.gz -2 sample\_2.paired.fq.gz #SamtoBam in Samtools, sorting bamfile samtools view -bS -o hisat2/accepted\_hits.bam hisat2/accepted\_hits.sam samtools sort -o hisat2/accepted\_hits.sorted.bam hisat2/accepted\_hits.bam #Make GTF file Stringtie - will be slightly modified as input, renamed to gff for GENESPACE stringtie hisat2/accepted\_hits.sorted.bam -o transcripts.gtf 5. Using Transdecoder to extract ORFs, using blast to search for homology in uniprot databse module load ncbi-blast/2.2.31+ module load R/4.2.2 /Transdecoder5.5/TransDecoder-TransDecoder-v5.5.0/util/gtf\_genome\_to\_cdna\_fasta.pl transcripts.gtf genome.fasta &gt; transcripts.TDCL.fasta /Transdecoder5.5/TransDecoder-TransDecoder-v5.5.0/TransDecoder.LongOrfs -t transcripts.TDCL.fasta blastp -query transcripts.TDCL.fasta.transdecoder\_dir/longest\_orfs.pep -db uniprot\_sprot.fasta -max\_target\_seqs 1 -outfmt 6 -evalue 1e-5 -num\_threads 16 &gt; blastp.outfmt6 /Transdecoder5.5/TransDecoder-TransDecoder-v5.5.0/TransDecoder.Predict -t transcripts.TDCL.fasta --retain\_blastp\_hits blastp.outfmt6 6. Using SpliceFinder\_2.pl to find longest isoform (SpliceFinder\_2.pl enclosed in scripts), output will be used as input for GENESPACE and reciprocal blasts module load blat perl SpliceFinder\_2.pl transcripts.TDCL.fasta.transdecoder.pep 7. Vim commands to alter the output of the gtf (step 4) and fasta file (step 6) to make it suitable for GENESPACE input Steps to edit gtf :%s/gene\_id /gene\_id=/g :%s/transcript\_id /transcript\_id=/g :%s/"//g ## rename transcripts.gtf to transcripts.gff ## rename transcripts.TDCL.fasta.transdecoder.pep.long to transcripts.TDCL.fasta.transdecoder.pep.long.faa ######################################################################################### ######################################################################################### Identifying the X chromosome in S. curviseta 1. Aligning reads from SRA to genome module load bowtie2/2.4.4 bowtie2-build GCA\_004115045.3\_ASM411504v3\_genomic.fna Sinella --threads 2 bowtie2 -x Sinella -1 SRR7948081\_1.fastq -2 SRR7948081\_2.fastq --end-to-end --sensitive -p 8 -S SRR7948081.sam --threads 8 2. Extract unique reads grep -vw "XS:i" SRR7948081.sam &gt; SRR7948081.unique.sam 3. Compute coverage in 10k windows module load soap/coverage soap.coverage -sam -cvg -i SRR7948081.unique.sam -onlyuniq -p 8 -refsingle GCA\_004115045.3\_ASM411504v3\_genomic.fna -window SRR7948081window 10000 4. R script to identify chromosome with lowest coverage, plot coverage per chromosomes (Figure S1) IdentifyX.Sinella.R ######################################################################################### ######################################################################################### Preparing files for T. cristinae (pipeline altered since no chromosome-wide genome assembly) 1. Download files from: https://github.com/DarrenJParker/Timema\_sex\_chr\_evol\_code/tree/main/data Tce\_b3v08.fasta Tce\_b3v08.max\_arth\_b2g\_droso\_b2g.gff TCE\_1000\_chr\_info\_and\_counts.txt 2. Extract AA sequences from gff and genome Download getAnnoFastaFromJoingenes.py: https://github.com/Gaius-Augustus/Augustus/tree/master getAnnoFastaFromJoingenes.py -g Tce\_b3v08.fasta -o Timema.fasta -3 Tce\_b3v08.max\_arth\_b2g\_droso\_b2g.gff -s FILTER 3. Make GFF from sequences maintained in step 2 (lots of sequences removed because of premature stop codons) grep "&gt;" Timema.fasta.aa &gt; Timema.sequence.txt (then remove &gt; in vim) python script to make new gff: get.filtered.gff.py (included) cut -d$'\\t' -f1,2,5,6 Timema.Sequence.gff &gt; Filter.txt Filter.txt - used as input in R scripts to assess homology and compute expectations Timema.fasta.aa - used as input for reciprocal blast analysis ######################################################################################### ######################################################################################### Commands to perform reciprocal blasts (identify 1:1 orthologs) between each species and the reference (commands identical for I. elegans or S. curviseta) 1. Blast: Reference Species as Query, Focal as target module load ncbi-blast/2.2.31+ makeblastdb -in Focal.fa -dbtype prot blastp -query Reference.fa -db Focal.fa -out blast.out.txt -evalue 1e-4 -num\_threads 4 -max\_target\_seqs 1 -outfmt 6 /trinityrnaseq-Trinity-v2.8.4/util/misc/blast\_outfmt6\_group\_segments.pl blast.out.txt Reference.fa Focal.fa &gt; blast.Ref.Foc.group.txt 2. Blast: Focal Species as Query, Reference as target makeblastdb -in Reference.fa -dbtype prot blastp -query Focal.fa -db Reference.fa -out blast.out.txt -evalue 1e-4 -num\_threads 4 -max\_target\_seqs 1 -outfmt 6 /trinityrnaseq-Trinity-v2.8.4/util/misc/blast\_outfmt6\_group\_segments.pl blast.out.txt Focal.fa Reference.fa &gt; blast.Foc.Ref.group.txt 3. Obtain reciprocal best hits - requires minimum length of 50%, minimum ID 40% (get\_recip\_50.py enclosed in scripts - requires python2) ######################################################################################### ######################################################################################### R Scripts to assess homology between X chromosomes Inputs required: Focal.gff, Reference.gff, Reciprocal Best hits (recip\_matches\_min50.txt, output from python script "get\_recip\_50.py above") Note that gffs require some light editing before input: change column chr -&gt; .chr change column ID -&gt; change column ord -&gt; .ord Grasshopper as an example: Grasshopper.chr start end Grasshopper strand Grasshopper.ord 1. Using Damselfly as reference Stats.Damsel.Expectation.Allacma.r Stats.Damsel.Expectation.Aphid.r Stats.Damsel.Expectation.Fly.r Stats.Damsel.Expectation.Grass.r Stats.Damsel.Expectation.Lacewing.r Stats.Damsel.Expectation.Ladybug.r Stats.Damsel.Expectation.Silkworm.r Stats.Damsel.Expectation.Stonefly.r Stats.Damsel.Expectation.Sinella.r Stats.Damsel.Expectation.Timema.r 2. Using Sinella (Springtail) as reference Stats.Sinella.Expectation.Allacma.r Stats.Sinella.Expectation.Aphid.r Stats.Sinella.Expectation.Fly.r Stats.Sinella.Expectation.Grass.r Stats.Sinella.Expectation.Lacewing.r Stats.Sinella.Expectation.Ladybug.r Stats.Sinella.Expectation.Silkworm.r Stats.Sinella.Expectation.Stonefly.r Stats.Sinella.Expectation.Sinella.r Stats.Sinella.Expectation.Timema.r ######################################################################################### ######################################################################################### R Scripts to compute significance using Monte Carlo approach Inputs required: Focal.gff, Reference.gff, Reciprocal Best hits (recip\_matches\_min50.txt, output from python script "get\_recip\_50.py above") 1. Using Damselfly as reference Stats.Damsel.Allacma.r Stats.Damsel.Aphid.r Stats.Damsel.Fly.r Stats.Damsel.Grass.r Stats.Damsel.Lacewing.r Stats.Damsel.Ladybug.r Stats.Damsel.Silkworm.r Stats.Damsel.Stonefly.r Stats.Damsel.Sinella.r Stats.Damsel.Timema.r 2. Using Sinella (Springtail) as reference Stats.Sinella.Allacma.r Stats.Sinella.Aphid.r Stats.Sinella.Fly.r Stats.Sinella.Grass.r Stats.Sinella.Lacewing.r Stats.Sinella.Ladybug.r Stats.Sinella.Silkworm.r Stats.Sinella.Stonefly.r Stats.Sinella.Timema.r ######################################################################################### ######################################################################################### Synteny plots in GENESPACE (v0.94) 1. Make synteny plot 4 taxa (Damselfly, Grasshopper, Lacewing, Ladybug) script: FourSpecies.Synteny.Genespace.R 2. Make pairwise synteny plots for 8 taxa against Damselfly script: Pairwise.Synteny.Genespace.R ######################################################################################### ######################################################################################### Identifying shared genes on the X Note: this was run with all species with an X (all groups except Silkworm) and again excluding Diptera Output is data used for Figure 2 1. Use peptide sequence output from GENESPACE processing as input for Orthofinder module load orthofinder orthofinder -f Orthofinder/ 2. Use custom python script to extract single-copy orthologs across all species Input: Orthofinder/Orthogroups/Orthogroups\_SingleCopyOrthologues.txt Input: Orthofinder/Orthogroups/Orthogroups.tsv Script: GetSingleCopy.py (included in scripts) Output: All9SC.txt 3. Use R Script to count number of shared 1:1 orthologs at each node Input: All9SC.txt (produced above) Input: GFF files from GENESPACE Script: GetGenesOnX9taxa.R (included in scripts) #########################################################################################