main.tex

\documentclass[12pt]{article}
\usepackage[utf8]{inputenc}
\usepackage[letterpaper,margin=1.0in]{geometry}
%% Some formatting stuff
\usepackage{authblk}
\usepackage{fancyhdr}
%\usepackage{lineno}
\usepackage{siunitx}
\usepackage{hyperref}
\usepackage{booktabs, array, longtable}
\usepackage{multirow}
\pagestyle{fancy}
\setlength{\headheight}{14.5pt} % Fix fancyhdr warning
% Custom citep command for biblatex
\newcommand{\citep}{\parencite}
\fancyhead[R]{\textbf{Cephalopod Dosage Compensation}}
% for figures
\usepackage{graphicx}
\usepackage{wrapfig}
\usepackage{float}
\usepackage{breakcites}
\usepackage{nameref}
\usepackage{amsmath}
\usepackage{xspace}


%\usepackage[figuresonly,nolists,nomarkers]{endfloat}
%\renewcommand{\processdelayedfloats}{}

\graphicspath{ {./figures/} }

\newcommand{\comment}[1]{{\color{blue} #1}}
\newcommand{\illex}{\textit{Illex illecebrosus}\xspace}

% for hyperlinks
\hypersetup{  
    colorlinks=true,
    citecolor=black,
    urlcolor=cyan,
    linkcolor=blue
    }
\urlstyle{same}

% for highlighting text
\usepackage{xcolor}
\usepackage{soul}
\usepackage{tikz}

% Colors for chromatin modifier tables
\definecolor{fcol}{HTML}{C0392B}
\definecolor{mcol}{HTML}{2471A3}
\definecolor{switchcol}{HTML}{8E44AD}
\newcommand{\isoswitch}{\textsuperscript{\textcolor{switchcol}{\textdagger}}}

% bibliography
\usepackage[backend=biber,style=authoryear]{biblatex}
\addbibresource{refs.bib}
\addbibresource{methods_refs_fixed.bib}

%\linenumbers
\renewcommand*{\bibfont}{\fontsize{10}{12}\selectfont}


\newcommand{\beginsupplement}{%
    \setcounter{table}{0}
    \renewcommand{\thetable}{S\arabic{table}}%
    \setcounter{figure}{0}
    \renewcommand{\thefigure}{S\arabic{figure}}%
    \renewcommand{\theHtable}{supp.\arabic{table}}%
    \renewcommand{\theHfigure}{supp.\arabic{figure}}%
}

\def\changemargin#1#2{\list{}{\rightmargin#2\leftmargin#1}\item[]}
\let\endchangemargin=\endlist 

\title{Dosage compensation of the Z chromosome is ancestral and conserved across coleoid cephalopods}
\author[1, 2]{Scott T. Small}
\author[1, 2]{Silas Tittes}
\author[3,4]{Thomas Desvignes}
\author[1,4]{John H. Postlethwait}
\author[1, 2]{Andrew D. Kern}
\affil[1]{\small{University of Oregon, Institute of Ecology and Evolution}}
\affil[2]{\small{University of Oregon, Department of Biology}}
\affil[3]{\small{Department of Biology, University of Alabama at Birmingham}}
\affil[4]{\small{University of Oregon, Institute of Neuroscience}}

\date{\small{\today{}}}

\begin{document}

\maketitle

\begin{abstract}
Chromosome-wide dosage compensation has evolved repeatedly in animals, but its mechanisms are
known from only a handful of model systems. Coleoid cephalopods (squids, cuttlefish, and octopuses)
share a ZZ/ZO sex determination system in which females carry a single copy of a large Z chromosome,
making dosage compensation a potential necessity and allowing a decisive test of whether compensation
is a lineage-wide feature and how it is implemented. Using a new chromosome-level genome assembly for the squid \textit{Illex illecebrosus}
and comparative transcriptomics spanning seven cephalopod species across $>$300 million years,
we show that Z-linked expression is near-equalized between sexes across coleoids, indicating that
dosage compensation originated early in coleoid evolution.
We find that the Z is unusually enriched for chromatin-modifying genes and that sex-biased isoform
usage of a small set of Z-linked chromatin regulators tracks variation in Z-linked expression,
implicating chromatin-state control in chromosome-wide equalization. Strikingly, LINE-rich regions
on the Z are enriched for genes that escape compensation, consistent with TE-associated heterochromatin
forming barriers to upregulation. We further identify cis-element motifs within gene bodies associated with
compensated Z genes and discover two sex-biased Z-linked long non-coding RNAs,
conserved across $>$300~MY of cephalopod divergence,
containing conserved RNA-binding protein motifs that suggest an additional RNA-mediated regulatory axis.
These results establish cephalopods as an independently evolved model for dosage compensation and reveal a mechanism that
reuses familiar components---chromatin regulators, cis-elements, repeats, and lncRNAs---in an
unexpected configuration.
\end{abstract}

%% Abstract option B
% \begin{abstract}
% Chromosome-wide dosage compensation has evolved in multiple animal lineages, yet its genomic
% basis is understood in only a few systems. Coleoid cephalopods share a conserved Z chromosome,
% providing an opportunity to test whether dosage compensation is widespread and to infer its
% molecular architecture in a deeply divergent clade. Here we present a chromosome-level genome
% assembly for the northern shortfin squid \textit{Illex illecebrosus} and analyze comparative
% transcriptomes from seven cephalopod species spanning five orders and >300 million years of
% divergence. Across species, Z-linked expression is near-equalized between sexes, indicating
% that Z dosage compensation originated early in coleoid evolution and predates the split
% between octopuses and decapods.
% The Z chromosome is strongly enriched for chromatin-modifying genes (16 of 188 conserved Z genes),
% and these genes are themselves compensated, consistent with a model in which Z-linked epigenetic
% factors participate in maintaining chromosome-wide regulatory state.
% Sex-biased isoform usage of four Z-linked chromatin regulators is associated with 57\% of
% variance in Z:autosome expression, with the H3K9 demethylase \textit{KDM3A} as the strongest
% predictor.
% In contrast to the mammalian X repeat hypothesis, LINE-dense regions on the Z are enriched
% for genes that escape compensation, consistent with TE-associated heterochromatin limiting
% upregulation. We identify gene-body cis-element motifs that discriminate compensated from
% escaped Z genes, functionally analogous to \textit{Drosophila} MSL Recognition Elements, and
% we discover two deeply conserved sex-biased Z-linked long non-coding RNAs (\textit{Zmast} and
% \textit{Zfest}) with conserved RNA-binding protein sites suggestive of a parallel splicing-regulatory
% axis. Together, these results establish coleoids as a tractable new model for the evolution of
% chromosome-wide dosage compensation and implicate an emergent mechanism integrating chromatin
% modifiers, cis-elements, repetitive-element domain structure, and regulatory lncRNAs.
% \end{abstract}

\section*{Introduction}
Dosage compensation is the regulatory mechanism by which gene expression between
individuals with differing numbers of sex chromosomes is equalized.
Multiple lineages have independently evolved mechanisms to address this imbalance:
in mammals, one X chromosome is transcriptionally silenced in females via the \textit{Xist} long non-coding RNA \citep{lyon1961gene, brockdorff2015dosage};
in \textit{Drosophila}, males upregulate their single X chromosome via the Male-Specific Lethal (MSL) complex \citep{lucchesi2015dosage, conrad2012dosage};
and in \textit{Caenorhabditis elegans}, hermaphrodites downregulate both X chromosomes by half \citep{meyer2005x}.

The diversity of these mechanisms reflects the fact that dosage compensation
has evolved independently in many lineages,
each time in response to the degeneration of the sex-limited chromosome
\citep{lucchesi1978dosage}.
As the Y (or W) chromosome loses functional genes,
selection favors compensatory upregulation of the remaining copy
in the heterogametic sex.
Studies of neo-sex chromosomes in \textit{Drosophila} have shown
that this process can occur rapidly and through diverse molecular paths,
with different species independently co-opting the MSL dosage compensation machinery
via distinct mutational mechanisms \citep{ellison2019convergent}.

Despite its evolutionary utility, dosage compensation is not ubiquitous or mechanistically uniform.
In birds (males ZZ, females ZW), bulk RNA-seq studies long suggested only incomplete,
gene-by-gene compensation \citep{itoh2007dosage, mank2009w},
but recent work has revealed multi-layered mechanisms---including
increased transcriptional burst frequency and elevated translational rates in females
\citep{papanicolaou2025multilayered}
and male-specific miRNA-mediated transcript degradation
\citep{fallahshahroudi2025mirna}---that
achieve near-complete protein-level balance.
In Lepidoptera, \textit{Bombyx mori} achieves compensation
through partial repression of both Z chromosomes in males \citep{rosin2022dosage},
and in livebearing fish, complete X chromosome compensation
has evolved de novo in some lineages but not others \citep{darolti2019extreme}.
The diversity of solutions---even among ZW systems---underscores
that the molecular paths to dosage balance are highly contingent
on the preexisting genetic architecture of each lineage.

Cephalopods were long thought to lack genetic sex chromosomes:
classical karyotype analyses failed to reveal heteromorphic chromosomes
in octopuses, squids, or cuttlefish,
leading to hypotheses of environmental or polygenic sex determination.
Recent genomic work has overturned this view.
In \textcite{coffing2025cephalopod}, we demonstrated that across cephalopods,
females are hemizygous for a large chromosome,
establishing a ZZ/ZO system shared across all coleoid (squids, cuttlefish, and octopuses) lineages surveyed.
This Z chromosome is conserved across coleoids,
originating over 300 million years ago in a common ancestor.
In \textit{Nautilus pompilius}, a representative of the only surviving non-coleoid cephalopod lineage,
the syntenic region is used as an X chromosome \citep{torrado2025nautilus},
indicating that the ancestral sex chromosome has been independently repurposed
in the two surviving cephalopod lineages.
These discoveries position coleoid cephalopods among the few known invertebrate
groups with deeply conserved sex chromosomes,
comparable in age to the insect X chromosome \citep{toups2023x}
and substantially older than mammalian ($\sim$170~MY) or avian ($\sim$100~MY) sex chromosomes.

The open question is then if and how cephalopods achieve dosage compensation. 
In a recent report, \textcite{papanicolaou2024z} provided the first evidence of dosage compensation
in any cephalopod, reporting M:F expression ratios of 1.04--1.15
across adult somatic tissues (brain sub-regions and arm) in two octopus species---notably
more complete than the 1.4--1.6 typical of birds \citep{mank2009w}.
They also identified two Z-linked lncRNAs with opposing sex specificities:
\textit{Zmast} (256-fold male-biased) and \textit{Zfest} (18-fold female-biased),
conserved between species separated by $\sim$2.5~MY,
whose molecular functions remain uncharacterized.
No comparable studies have been conducted outside of octopus,
leaving the prevalence and molecular basis of compensation across cephalopods unknown.

Here we present a chromosome-scale genome assembly of the northern shortfin squid,
\illex, and identify its Z chromosome via synteny and coverage analyses.
We then use transcriptomic data to test for dosage compensation across seven cephalopod species
spanning five orders and over 300 million years of divergence,
establishing that compensation of the Z chromosome is ancestral to coleoid cephalopods.
To investigate the molecular underpinnings of this compensation,
we characterize DNA methylation patterns,
identify an enrichment of chromatin-modifying genes on the conserved Z chromosome,
and show that sex-biased isoform switching of chromatin remodelers
predicts chromosome-wide compensation levels.
We find that LINE elements and CTCF insulators define structural domains on the Z,
paralleling Lyon's repeat hypothesis for the mammalian X \citep{lyon1998x},
and discover gene-body cis-elements that distinguish compensated Z genes from escapees,
analogous to the MSL Recognition Elements in \textit{Drosophila}.
Finally, we show that the sex-specific lncRNAs \textit{Zmast} and \textit{Zfest}
are conserved across all species examined.
Together, these results reveal a multi-layered dosage compensation system
in one of the oldest known sex chromosome systems in animals.


\begin{figure}[tbp]
    \centering
    \includegraphics[width=\linewidth]{figures/Figure1_update.pdf}
    \caption{Conserved synteny of seven coleoid cephalopod genomes.
    An uncalibrated ultrametric phylogeny (left) relates the seven species used in this study.
    Ribbons connect syntenic blocks between adjacent species, colored by \illex\ chromosome identity.
    The Z chromosome (red, rightmost) is conserved as a single syntenic unit across all species,
    although its name varies (Z in \illex, \textit{O.\ bimaculoides}, and \textit{Sepia esculenta};
    chromosome 46 in \textit{Sthenoteuthis oualaniensis}; 43 in \textit{D.\ pealeii}; 43 in \textit{E.\ scolopes};
    NC16 in \textit{O.\ sinensis}).}
    \label{fig:phylo}
\end{figure}

\section*{Results}
\subsection*{A chromosome-scale genome assembly for \textit{Illex illecebrosus}}
We generated a chromosome-scale genome assembly for an adult female \illex
using PacBio HiFi long reads, Illumina short reads, and Hi-C chromatin conformation capture data.
Hi-C scaffolding resolved 46 chromosomes,
consistent with the decapodiform karyotype of 2n = 92 \citep{gao1990karyological},
and the final assembly spans 4.14 Gb with 95.6\% BUSCO completeness
(mollusca\_odb12; Table~\ref{tab:assembly_stats}).
Gene annotation combining stranded total RNA-seq from both individuals
and PacBio Iso-Seq from the female (mantle tissue; see Methods),
together with \textit{ab initio} prediction, identified 23,192 protein-coding genes.
Because the RNA-seq libraries were prepared with a total RNA protocol (not poly-A selected),
non-polyadenylated transcripts including many lncRNAs are represented.
We also generated deep PacBio HiFi sequencing from a second adult male individual,
enabling direct comparison of read depth between sexes across all chromosomes.
Full assembly statistics and repeat composition are provided in the
\hyperref[sec:supp_assembly]{Supplemental Results}
(Tables~\ref{tab:assembly_stats} and \ref{tab:repeat_stats}).

\subsection*{Identification and conserved synteny of the Z chromosome in \illex}
Comparison of male and female read depth across assembled chromosomes
revealed that a single chromosome is at hemizygous coverage
in the female relative to the male, which we inferred to be the Z chromosome (Figure~\ref{fig:seq_depth}).

We next performed comparative genomic analyses with other cephalopod species for which sex chromosomes have been characterized.
Whole-genome conserved synteny analysis revealed a striking pattern of chromosomal conservation across coleoid cephalopods (Figure~\ref{fig:phylo}).
We compared the \illex assembly to published genomes of
\textit{Octopus bimaculoides}, \textit{O. sinensis}, \textit{Euprymna scolopes}, \textit{Sepia esculenta}, \textit{Doryteuthis pealeii}, and \textit{Sthenoteuthis oualaniensis}.
Chromosome-scale alignments showed that within decapods, synteny across all chromosome arms
is largely maintained, with a number of rearrangements, fusions, and fissions.
However, when comparing out to octopus, synteny breaks down on every chromosome except the Z.
The Z chromosome has maintained its integrity as a distinct chromosomal unit for over 300 million years of evolution,
making it one of the oldest known sex chromosome systems in animals.
%These syntenic relationships, combined with coverage analysis showing reduced female-to-male mapping ratios on chromosome 42,
%confirm that this chromosome represents the Z chromosome in \illex, consistent with the ZZ/ZO sex determination system described in other coleoid cephalopods.

\subsection*{Partial dosage compensation of the Z chromosome in \illex}

To assess whether \illex exhibits dosage compensation for Z-linked genes,
we analyzed sex-stratified RNA-seq data from somatic mantle tissue of adults.
We quantified gene expression levels across all chromosomes in a single male (ZZ) and female (ZO),
the same individuals used for genome assembly and coverage analysis,
using DESeq2 normalization (Figure~\ref{fig:cross_species_dc}).
We use two complementary metrics throughout:
Z:A is the ratio of mean Z-linked to mean autosomal expression level (TPM, filtered to genes with TPM~$> 1$;
a value of 1.0 indicates Z expression equal to the autosomal baseline),
and M:F is the per-gene median ratio of male to female expression
(1.0 = complete compensation; 2.0 = no compensation;
percent equalization = $(1/\text{M:F}) \times 100$,
e.g.\ M:F of 1.15 gives $(1/1.15) \times 100 = 87\%$).
Expression levels on autosomes were comparable between sexes,
and Z-linked gene expression in females was comparable to autosomal levels
despite females carrying only a single Z chromosome (Z:A = 1.59;
values above 1.0 in both sexes reflect the generally higher expression
of Z-linked genes relative to the genomic average),
indicating that the hemizygous Z undergoes transcriptional upregulation.
However, Z-linked genes showed modest but significant male bias,
with a median M:F ratio of 1.15 compared to 0.99 for autosomes
(p = 0.037, Wilcoxon rank-sum test; Figure~\ref{fig:cross_species_dc}).
These results are consistent with dosage compensation of the Z chromosome in \illex,
with females achieving $\sim$87\% of male expression levels for Z-linked genes,
though the absence of biological replicates for this species
limits the precision of this estimate.

\begin{figure}[tbp]
    \centering
    \includegraphics[width=\textwidth]{figures/illex_dosage_compensation_3panel.pdf}
    \caption{Dosage compensation of the Z chromosome in \illex.
    Left panels: gene expression levels [log2(normalized count + 1)] across all chromosomes
    in females (ZO, top) and males (ZZ, bottom).
    The Z chromosome (orange) shows expression comparable to autosomes (blue) in both sexes,
    with Z:A ratios of 1.59 in females and 1.44 in males.
    Right panel: male-to-female expression ratios for autosomal genes (n = 11,459)
    and Z-linked genes (n = 188).
    The Z chromosome shows a modest but significant male bias
    (median M:F = 1.15 vs.\ 0.99 for autosomes; p = 0.037, Wilcoxon rank-sum test),
    consistent with partial dosage compensation (n = 1 per sex).
    Variation in median expression among autosomes (e.g.\ the relatively low values for chr15) reflects differences in gene number and composition across chromosomes.
    Dosage compensation results for six additional cephalopod species, most with biological replication, are shown in Figures S1--S6.}
    \label{fig:cross_species_dc}
\end{figure}

\begin{table}[tbp]
\centering
\caption{Dosage compensation of the Z chromosome across coleoid cephalopod species.
Z~M:F and Auto~M:F are the median per-gene male-to-female expression ratios
for Z-linked and autosomal genes, respectively
(1.0 = complete compensation; 2.0 = no compensation).
The p-value is from a Wilcoxon rank-sum test comparing Z and autosomal M:F distributions.
n~(F/M) indicates the number of female and male RNA-seq libraries.}
\label{tab:cross_species_dc}
\begin{tabular}{llrrrr}
\toprule
\textbf{Species} & \textbf{Order} & \textbf{Z M:F} & \textbf{Auto M:F} & \textbf{p-value} & \textbf{n (F/M)} \\
\midrule
\textit{Illex illecebrosus} & Oegopsida & 1.15 & 0.99 & 0.037 & 1/1 \\
\textit{Sthenoteuthis oualaniensis}$^*$ & Oegopsida & 1.13 & 1.04 & 0.010 & 16/16 \\
\textit{Doryteuthis pealeii} & Myopsida & 1.10 & 1.12 & 0.29 & 10/9 \\
\textit{Euprymna scolopes} & Sepiolida & 1.08 & 1.00 & 0.003 & 3/3 \\
\textit{Euprymna berryi} & Sepiolida & 0.97 & 0.97 & 0.48 & 2/2 \\
\textit{Sepia officinalis} & Sepiida & 1.03 & 0.99 & 0.021 & 4/4 \\
\textit{Octopus bimaculoides} & Octopoda & 1.08 & 1.03 & 0.20 & 15/15 \\
\bottomrule
\end{tabular}
\vspace{0.5em}

{\small $^*$Technical replicates from a single individual per sex; lacks biological replication.}
\end{table}

\subsection*{Dosage compensation is ancestral and conserved across coleoid cephalopods}

To determine whether dosage compensation of the Z chromosome is a shared feature of coleoid cephalopods
or a lineage-specific adaptation in \illex,
we analyzed sex-stratified, publicly available RNA-seq data from
six additional species representing five coleoid orders
across both major coleoid lineages (Decapodiformes and Octopodiformes):
\textit{Sthenoteuthis oualaniensis} (Oegopsida),
\textit{Doryteuthis pealeii} (Myopsida),
\textit{Euprymna scolopes} and \textit{Euprymna berryi} (Sepiolida),
\textit{Sepia officinalis} (Sepiida),
and \textit{Octopus bimaculoides} (Octopoda).
Sample sizes per species range from one to 15 libraries per sex
(Table~\ref{tab:cross_species_dc}).
For \textit{S. oualaniensis}, the 16 libraries per sex are technical replicates
from a single individual; for \textit{O. bimaculoides}, the 15 libraries per sex
represent multiple tissues from several individuals.
All other species have independent biological replicates.
For each species, we identified Z-linked genes via synteny with the \illex Z chromosome
and compared expression levels between males (ZZ) and females (ZO).

In every species examined, Z-linked gene expression in females was upregulated
to match or approach the level observed in males,
indicating active dosage compensation of the hemizygous Z chromosome
(Figure~\ref{fig:cross_species_dc}; Figures S1--S6).
Median male-to-female expression ratios for Z-linked genes ranged from 0.97 to 1.15 across species
(Table~\ref{tab:cross_species_dc}),
all within 15\% of 1.0 and far below the 2.0 expected in the absence of compensation.
Autosomal genes showed M:F ratios centered at 1.0 in all species,
confirming that the observed Z chromosome patterns reflect sex-linked regulation
rather than systematic biases in library preparation or normalization.

In three of the seven species---\textit{D. pealeii},
\textit{E. berryi}, and \textit{O. bimaculoides}---Z chromosome
M:F ratios were indistinguishable from autosomal ratios
(p = 0.20--0.48; Table~\ref{tab:cross_species_dc}),
indicating complete dosage compensation.
\textit{E. scolopes} and \textit{S. officinalis} showed statistically significant but minimal male bias
(M:F = 1.03--1.08; p $<$ 0.02), indicating near-complete compensation.
In \illex{} (n = 1/1) and \textit{S. oualaniensis},
Z-linked genes showed modest but significant residual male bias
(M:F = 1.13--1.15; p $<$ 0.05),
consistent with partial compensation---though the \illex{} estimate
should be interpreted cautiously given the lack of biological replication.
Even in these species, the degree of male bias is far less
than that observed in birds, where Z-linked M:F ratios typically range from 1.4 to 1.6
\citep{mank2009w}.

The presence of dosage compensation across all seven species,
spanning five orders and over 300 million years of divergence,
demonstrates that compensation of the Z chromosome is ancestral to coleoid cephalopods
and has been conserved throughout their evolutionary history.
Importantly, cephalopod dosage compensation is an independent evolutionary origin of dosage compensation 
among major animal lineages, distinct from the mechanisms observed in mammals, insects, and birds,
and thus offers the opportunity to discover novel mechanisms of sex chromosome regulation.


\begin{figure}[tbp]
    \centering
    \begin{tikzpicture}
      \node[inner sep=0pt] (fig) {\includegraphics[height=\textwidth, angle=90, trim=9.5cm 4cm 4.3cm 2.1cm, clip]{figures/fig3.pdf}};
      \fill[white] ([xshift=-0.2cm,yshift=0.2cm]fig.north west) rectangle ++(2cm,-0.9cm);
    \end{tikzpicture}
    \caption{Female-specific hypermethylation of the Z chromosome in \illex.
    (A) Differential DNA methylation (Male $-$ Female) across all chromosomes.
    Autosomal chromosomes (blue) show values centered near zero, indicating equivalent methylation between sexes.
    The Z chromosome (red) shows strongly negative values, indicating higher methylation in females.
    (B) Summary comparison of differential methylation between autosomes and the Z chromosome.
    The Z chromosome is significantly hypermethylated in females relative to males
    (p = 5.37$\times$10$^{-44}$, Wilcoxon rank-sum test).
    The extreme statistical significance reflects the large number of individual genes
    in the rank-sum test (188 Z-linked vs.\ 11,459 autosomal), which provides
    high power to detect distributional shifts even when box-plot whiskers overlap.}
    \label{fig:methylation}
\end{figure}
\subsection*{DNA methylation patterns on the Z chromosome}

The conservation of dosage compensation across coleoid cephalopods
raises the question of what molecular mechanisms underlie this regulation.
Given the role of epigenetic modifications in gene regulation and dosage compensation in other systems,
we investigated DNA methylation patterns across the \illex genome.
We leveraged PacBio HiFi kinetic signatures from the male and female individual to assess CpG methylation levels.
Genome-wide, we observed a positive correlation between gene expression and methylation levels in both sexes,
with more highly expressed genes showing elevated methylation (females: r=0.396, p=2.16$\times$10$^{-9}$; males: r=0.318, p=5.10$\times$10$^{-6}$).
This pattern is consistent with gene-body methylation associated with active transcription,
as has been observed across invertebrates \citep{zemach2010genome, glastad2011dna}
including molluscs \citep{manner2021inference}.

Strikingly, when comparing methylation levels between sexes across chromosomes,
we found that the Z chromosome is significantly hypermethylated in females relative to males (Figure~\ref{fig:methylation}).
Differential methylation (Male $-$ Female) was calculated for each chromosome,
revealing that autosomal chromosomes show values centered near zero,
indicating equivalent methylation between sexes (Figure~\ref{fig:methylation}A).
In contrast, the Z chromosome exhibited strongly negative differential methylation,
demonstrating that females have substantially higher methylation on the Z than males
(p = 5.37$\times$10$^{-44}$, Wilcoxon rank-sum test; Figure~\ref{fig:methylation}B).
The Z chromosome also shows greater variance in differential methylation than individual autosomes,
likely reflecting the smaller number of Z-linked genes and the hemizygous state in females,
which eliminates allelic averaging.

This female-specific hypermethylation of the Z chromosome is notable given the positive correlation
between gene-body methylation and expression observed genome-wide.
The pattern suggests that elevated gene-body methylation on the hemizygous female Z
may contribute to its transcriptional upregulation toward diploid levels,
consistent with gene-body methylation facilitating active transcription.
If this association reflects a conserved feature of coleoid dosage compensation,
it would represent the longest-maintained epigenetic sex-chromosome signature described
in any animal---exceeding the $\sim$170~MY mammalian and $\sim$100~MY avian systems
by more than 100~MY.
Direct comparison of Z versus autosomal methylation within each sex
confirms that the effect is larger in females than males
(\hyperref[sec:supp_methylation]{Supplemental Results}; Figure~\ref{fig:meth_auto_vs_z}).
We report differential methylation as Male~$-$~Female rather than as a ratio
because the subtraction is symmetric around zero and directly interpretable
when comparing chromosomes with different baseline methylation levels.
However, the relationship between methylation and dosage compensation remains complex:
gene-body methylation correlates with accumulated transcript levels,
which reflect both transcription rate and mRNA stability,
and methylation differences alone do not fully explain the residual male bias in Z-linked expression.


\subsection*{The Z chromosome is enriched for chromatin-modifying genes}
Strong conservation of the Z chromosome across coleoid cephalopods implies
that its gene content has been subject to long-term selective pressures that
likely include the maintenance of dosage compensation.
If dosage compensation depends on chromosome-wide epigenetic regulation,
then genes encoding proteins that modify chromatin structure---histone
methyltransferases, acetyltransferases, remodelers, and related
factors (hereafter ``chromatin-modifying genes'')---might be selectively
retained on the Z to maintain \textit{cis}-regulatory control.
To test this prediction, we used reciprocal best-hit (RBH) orthology mapping
to identify genes with conserved Z-linkage
at increasing levels of phylogenetic stringency.
This gene-level orthology analysis complements the whole-chromosome
synteny shown in Figure~\ref{fig:phylo}: whereas synteny analysis reveals
that the Z chromosome is preserved as a macro-syntenic block,
the orthology approach asks which individual genes have remained on the Z
versus translocated to autosomes.
For these analyses, we included an eighth species,
\textit{Octopus vulgaris}, whose genome is available
and in which dosage compensation has been recently reported \citep{papanicolaou2024z}.
We anchored the analysis on \textit{O. bimaculoides}, for which we have the
deepest expression data (30 somatic samples; see below).
Of 639 genes annotated on the \textit{O. bimaculoides} Z chromosome,
20 are chromatin-modifying genes.
Requiring that a gene be Z-linked in all four of
\illex, \textit{O. bimaculoides}, \textit{S. oualaniensis}, and \textit{O. vulgaris}
retained 188 genes (29\% of the 639);
extending the requirement to all eight species reduced this to 132
(Table~\ref{tab:conserved_z_genes}).
The drop from 188 to 132 conserved genes when extending from four
to eight species represents 70\% retention---modest given the
$>$300~MY evolutionary distances involved.

\begin{table}[tbp]
\centering
\caption{Chromatin gene enrichment on the Z chromosome at increasing phylogenetic stringency.
As Z-linkage is required across more species, chromatin genes are preferentially retained.
OR and p-value are from one-sided Fisher's exact test comparing
conserved versus non-conserved Z genes (chromatin versus non-chromatin;
20 total chromatin genes among 639 Z-linked genes in \textit{O. bimaculoides}).}
\label{tab:z_enrichment}
\begin{tabular}{lrrrrl}
\toprule
\textbf{Stringency} & \textbf{Total genes} & \textbf{Chromatin} & \textbf{\% Chromatin} & \textbf{OR} & \textbf{p-value} \\
\midrule
4 species & 188 & 16 & 8.5\% & 10.4 & $2.9 \times 10^{-6}$ \\
5 species & 169 & 16 & 9.5\% & 12.2 & $5.6 \times 10^{-7}$ \\
6 species & 160 & 14 & 8.8\% & 7.6 & $2.2 \times 10^{-5}$ \\
7 species & 148 & 13 & 8.8\% & 6.7 & $5.9 \times 10^{-5}$ \\
8 species & 132 & 13 & 9.8\% & 7.8 & $1.5 \times 10^{-5}$ \\
\bottomrule
\end{tabular}
\end{table}

At every stringency level, the conserved Z-linked gene set
was strongly enriched for chromatin-modifying and epigenetic regulatory functions.
Of the 188 genes conserved across four species,
16 are chromatin-associated (8.5\%)---a $\sim$10-fold enrichment
over the 4 chromatin genes among 451 non-conserved Z genes
(Fisher exact $p < 10^{-6}$).
Remarkably, 13 of these 16 chromatin genes remain Z-linked in all eight species,
demonstrating that chromatin modifiers are among the most deeply conserved
elements of the Z chromosome
(Table~\ref{tab:z_enrichment}).

Among the Z-linked chromatin genes conserved across all eight species,
the most striking feature is the co-localization of regulators
from multiple independent chromatin-modifying pathways.
These include the H3K4 methyltransferase \textit{SET1B},
the MLL-complex scaffold \textit{MEN1},
the H2B ubiquitin ligase \textit{BRE1/RNF40},
the histone acetyltransferase \textit{CREBBP/CBP},
the SWR1 H2A.Z remodeler \textit{Domino/SRCAP},
the HIRA-complex chaperone \textit{UBN1},
and the heterochromatin regulator \textit{BAHD1}---all Z-linked in every species examined.
Additional deeply conserved Z-linked chromatin modifiers include
\textit{ING5} (HBO1/MOZ-associated PHD reader),
\textit{KIF4A} (chromatin condensation),
\textit{ELM2-domain} (NuRD/HDAC remodeler),
\textit{THOC6} (THO/TREX mRNA export), and \textit{PAN3} (mRNA deadenylase).
\textit{TDRD3} (Tudor-domain methyl reader), present in the 4-species core,
drops out at 6 species, suggesting lineage-specific loss from the Z.
More generally, translocation of a Z-linked gene to an autosome
would remove it from the compensation machinery that upregulates
Z-linked loci in females, potentially producing a dosage imbalance
that selects against such rearrangements.
Outright gene loss, by contrast, escapes this constraint.
Notably, the WRAD scaffold subunits shared by all COMPASS-family complexes
(\textit{WDR5}, \textit{RBBP5}, \textit{ASH2L}) are autosomal,
indicating that the Z-linked chromatin regulators
do not constitute a single complex
but rather represent convergent enrichment
of multiple independent chromatin pathways on the sex chromosome.
The involvement of multiple pathways---histone methylation, acetylation,
variant-histone deposition, and chromatin remodeling---is consistent
with the multi-layered nature of dosage compensation itself,
which likely requires coordinated regulation across several chromatin marks
rather than a single switch.

\begin{table}[tbp]
\centering
\caption{Expression of conserved Z-linked chromatin modifier genes
in \textit{O. bimaculoides} (30 somatic samples, 15F:15M).
log$_2$FC indicates male-to-female fold change; positive values indicate male bias.}
\label{tab:compass}
\vspace{0.5em}
\begin{tabular}{llrr}
\toprule
\textbf{Gene} & \textbf{Pathway} & \textbf{log$_2$FC} & \textbf{$p_\text{adj}$} \\
\midrule
\textit{THOC6} & THO/TREX & +0.52 & $9.6 \times 10^{-4}$ \\
\textit{Domino/SRCAP} & SWR1/H2A.Z & +0.59 & 0.027 \\
\textit{PAN3} & Deadenylation & +0.94 & 0.076 \\
\textit{CREBBP/CBP} & HAT & +0.61 & 0.089 \\
\textit{ELM2-domain} & NuRD/HDAC & +1.02 & 0.095 \\
\textit{MEN1} & MLL/H3K4me & +1.21 & 0.175 \\
\textit{TDRD3} & Tudor reader & $-$0.32 & 0.386 \\
\textit{SET1B} & H3K4me & +0.18 & 0.474 \\
\textit{ING5} & HBO1/MOZ & +0.15 & 0.571 \\
\textit{UBN1} & HIRA/H3.3 & +0.13 & 0.606 \\
\textit{BAHD1} & Heterochromatin & +0.23 & 0.628 \\
\textit{BRE1/RNF40} & H2Bub1 & +0.05 & 0.774 \\
\textit{KIF4A} & Condensation & +0.19 & 0.928 \\
\bottomrule
\end{tabular}
\end{table}

\subsection*{Sex-biased chromatin modifier expression in \textit{O. bimaculoides}}

To determine whether these Z-linked chromatin genes escape dosage compensation
or are themselves compensated,
we examined their expression in \textit{O. bimaculoides},
which provides the greatest statistical power
(30 somatic samples, 15F:15M).
Most conserved Z-linked chromatin genes show balanced expression
between sexes (Table~\ref{tab:compass}):
\textit{SET1B}, \textit{BRE1}, \textit{ING5}, \textit{TDRD3}, \textit{UBN1}, \textit{BAHD1}, \textit{CREBBP}, and \textit{KIF4A}
are all non-significant ($p_\text{adj} > 0.05$).
Three genes deviate from this pattern:
\textit{MEN1} trends male-biased (log$_2$FC = +1.21, $p_\text{adj}$ = 0.175),
\textit{THOC6} is significantly male-biased
(+0.52, $p_\text{adj}$ = $9.6 \times 10^{-4}$),
and \textit{ELM2-domain} trends male (+1.02, $p_\text{adj}$ = 0.095).
The compensation of most Z-linked chromatin genes is consistent with
a self-reinforcing regulatory loop:
if these genes encode the epigenetic machinery
that maintains permissive chromatin on the Z chromosome,
their own compensation would sustain expression from a single copy in females.

Genome-wide, PyDESeq2 identified 2,470 significantly sex-biased genes
in \textit{O. bimaculoides}
($p_\text{adj} < 0.05$, $|\text{log}_2\text{FC}| > 1$),
including 27 chromatin modifier isoforms (15 female-biased, 12 male-biased)
spanning SWI/SNF, KDM, CHD, HAT, and heterochromatin reader families
(\hyperref[sec:supp_chromatin]{Supplemental Results};
Tables~\ref{tab:chromatin-female} and \ref{tab:chromatin-male}).
Among these, four genes---\textit{PBRM1}, \textit{SMARCC2}, \textit{KDM3A}, and \textit{Domino}---show
significant isoforms in \emph{both} directions,
producing sex-biased isoform switching rather than simple up- or down-regulation.
In \textit{PBRM1} (Polybromo-1, PBAF complex),
the dominant female isoforms (.2/.4) and male isoform (.1) differ by 10-fold
in expression magnitude ($|\text{log}_2\text{FC}| > 9$),
indicating near-exclusive sex-specific promoter usage.
\textit{SMARCC2} (BAF170, SWI/SNF core subunit) and \textit{KDM3A} (H3K9 demethylase)
show more moderate switching ($|\text{log}_2\text{FC}|$ of 1--2),
and SUPPA2 alternative splicing analysis confirms sex-biased isoform usage
in both genes (\textit{SMARCC2}: 9 significant events; \textit{CHD5}: 6 significant events).
In both cases, females predominantly express full-length isoforms
while males produce truncated or NMD-targeted variants from alternative promoters.

To test whether sex-biased isoform usage of chromatin modifiers
directly predicts Z chromosome compensation levels,
we regressed expression-matched Z:A ratios on female isoform proportions
of the four genes with bidirectional isoform switching:
\textit{PBRM1}, \textit{SMARCC2}, \textit{KDM3A}, and \textit{Domino}.
A multiple regression model explained 56.8\% of Z:A variance
across 30 \textit{O. bimaculoides} somatic samples
($R^2 = 0.568$, $F = 8.20$, $p = 2.27 \times 10^{-4}$;
Figure~\ref{fig:za_regression}),
four-fold more than sex alone ($R^2 = 0.144$).
\textit{KDM3A}, an H3K9 demethylase, was the sole individually significant predictor
(coefficient = +17.96, $p = 0.008$, Bonferroni $p = 0.032$;
Figure~\ref{fig:kdm3a_za}):
samples with higher proportions of the female \textit{KDM3A} isoform
showed greater Z:A compensation.
This result links isoform-level regulation of a specific chromatin modifier
to chromosome-wide dosage compensation output,
consistent with a model in which sex-biased alternative splicing
of chromatin remodelers tunes Z chromosome expression.

\subsection*{LINE enrichment on the Z chromosome parallels the mammalian X}

Lyon (1998) proposed that LINE-1 elements on the mammalian X chromosome
serve as ``way stations'' for the spreading of X-inactivation from the \textit{Xist} locus,
an idea supported by the observation that the human X has twice the LINE density of autosomes
and that LINE-poor regions preferentially escape inactivation.
More generally, repeat-rich landscapes may facilitate
chromosome-scale chromatin states
by providing a homogeneous substrate
for epigenetic machinery to recognize and propagate along a chromosome
\citep[e.g.,][]{huang2022species}.
We note, however, that TE accumulation on sex chromosomes
can also arise from reduced effective population size
and recombination rates
rather than positive selection for a regulatory function;
the two explanations are not mutually exclusive.

In \textcite{coffing2025cephalopod} we discovered that LINE enrichment
is a conserved feature of cephalopod Z chromosomes:
in \textit{O. bimaculoides}, the Z harbors LINEs
at approximately twice the density of the genome-wide background
(Mann-Whitney $U$ test on 1~Mb windows, $p < 10^{-4}$),
and the same pattern holds in \textit{O. sinensis},
\textit{E. scolopes}, and \textit{S. esculenta}
(all $p < 10^{-3}$).
This led us to propose LINE density
as a signature of Z chromosomes across cephalopods.

We confirmed this enrichment in \textit{O. bimaculoides}
(Z ranks 1st of 30 chromosomes at 212,484 bp/Mb
versus an autosomal mean of 136,429 bp/Mb;
1.56-fold, $z = 4.92$ relative to the autosomal standard deviation)
and additionally found that SINEs are correspondingly depleted on Z
(42,370 bp/Mb versus an autosomal mean of 51,737 bp/Mb;
0.82-fold, $z = -1.68$),
producing a LINE-high/SINE-low signature
that parallels the mammalian X chromosome.
The dominant LINE family on Z by total content is RTE-BovB
(5.66~Mb, 45\% of Z LINE content; 1.53-fold enriched),
but the most enriched family is L1-Tx1 (2.06-fold, 2.38~Mb),
followed by Dong-R4 (1.69-fold, 0.71~Mb).

Crucially, LINE density on the Z predicts which genes escape dosage compensation.
Male-biased Z genes---those retaining the default twofold excess
expected from ZZ dosage---reside in regions with significantly
higher LINE density than female-biased genes
(mean 174.2 vs.\ 122.0 bp/kb in 10~kb flanking windows centered on each gene;
Mann-Whitney $U$ test, $p = 0.029$;
Figure~\ref{fig:line_density}).
This signal is driven by the aggregate LINE landscape rather than any single family:
RTE-BovB density alone does not differ significantly between sex-bias categories
(Mann-Whitney $p = 0.251$; Figure~\ref{fig:line_density}),
consistent with multiple LINE families contributing to the heterochromatic barrier.
The effect is nonetheless strikingly categorical:
in the highest RTE-BovB density quartile,
89\% of sex-biased genes are male-biased (8M:1F),
while lower quartiles show no consistent trend
(Figure~\ref{fig:rte_quartile}).

RTE-BovB and Dong-R4 occupy distinct spatial compartments on Z
(Spearman correlations across 500~kb sliding windows).
RTE-BovB is strongly correlated with gene density
($\rho$ = +0.68, $p = 7.4 \times 10^{-17}$),
marking the euchromatic compartment,
while Dong-R4 anticorrelates with gene density
($\rho$ = $-$0.60, $p = 1.3 \times 10^{-12}$)
and with RTE-BovB ($\rho$ = $-$0.61, $p = 4.6 \times 10^{-13}$),
defining the heterochromatic compartment.

FIMO scanning for CTCF binding motifs (JASPAR MA0139.2, $p < 10^{-4}$)
revealed that the Z chromosome also carries the highest CTCF density
of any chromosome after correcting for GC content and chromosome size
(+6.49 motifs/Mb above expectation;
rank 1/30 chromosomes, $z = 3.85$ relative to autosomal distribution).
The excess is predominantly intergenic (88\% of $\sim$460 excess sites),
consistent with insulator function.
One exceptional locus contains 35 CTCF motifs within a 44~kb intergenic region
separating the female-biased transcription factor \textit{MEIS2}
(log$_2$FC = $-$1.90, $p_\text{adj}$ = 0.043) from an adjacent male-biased gene
(log$_2$FC = $+$2.12, $p_\text{adj}$ = 0.036),
consistent with a strong insulator boundary between oppositely regulated domains.
CTCF and RTE-BovB are positively correlated along the Z
($\rho$ = +0.534, $p = 2.6 \times 10^{-18}$),
but at the local scale, CTCF motifs occupy gaps between LINE elements
rather than their immediate flanks
(observed inter-element distance 4,318 bp vs.\ expected 715~bp).

Hi-C chromatin conformation data from female \illex
confirmed this structural organization.
The Z chromosome exhibited the strongest A/B compartmentalization
of any chromosome in the genome ($z = 3.36$ for PC1 variance).
All Z-linked chromatin modifier genes reside in the A compartment (euchromatic),
while the pioneer transcription factor \textit{FoxA2}
uniquely occupies the B compartment---consistent with its role
in opening closed chromatin.
TE families define compartment boundaries:
in \illex, Penelope and RTE-RTE elements mark the B compartment
(Spearman $\rho$ with PC1 = $-$0.377 and $-$0.284, respectively),
while Satellite elements mark the A compartment ($\rho$ = +0.168).
Different TE families define compartments in different species---RTE-BovB in
\textit{O. bimaculoides} versus Penelope in \illex---but
the structural principle of TE-defined chromatin domains on Z is conserved.
Thus, while LINEs on the mammalian X are thought to facilitate the
\textit{spreading} of dosage compensation \citep{lyon1998x},
LINEs on the cephalopod Z instead demarcate regions that
\textit{escape} it---an inverse relationship that may reflect
fundamentally different compensation mechanisms acting on a shared
repeat architecture.

\begin{figure}[tbp]
    \centering
    \includegraphics[width=0.9\linewidth]{figures/line_density_sex_bias_boxplots.png}
    \caption{LINE enrichment on the \textit{O. bimaculoides} Z chromosome predicts sex-biased expression.
    Left: total LINE density in 10~kb flanking regions of Z-linked genes stratified by sex-bias category.
    Male-biased genes reside in LINE-dense regions (174.2 bp/kb)
    while female-biased genes occupy LINE-sparse regions (122.0 bp/kb; Mann-Whitney $p = 0.029$).
    Right: RTE-BovB density alone does not differ significantly among categories ($p = 0.251$),
    indicating that the signal reflects the aggregate LINE landscape rather than a single family.}
    \label{fig:line_density}
\end{figure}

\subsection*{Gene-body cis-elements distinguish compensated Z genes from escapees}

In \textit{Drosophila}, the MSL complex achieves dosage compensation
by recognizing $\sim$150 Chromatin Entry Sites (CES)
that overlap gene bodies rather than promoters.
CES contain a 21-bp GA-rich MSL Recognition Element (MRE),
which is modestly enriched on the X ($\sim$2-fold)
and whose presence rather than copy number determines MSL targeting
\citep{alekseyenko2008high}.
To search for analogous cis-elements in cephalopods,
we performed discriminative motif discovery (STREME)
comparing compensated and male-biased (``escapee'') Z genes in \textit{O. bimaculoides},
identifying 10 Z-exclusive motifs enriched in compensated gene sequences.

When scanned across gene bodies (gene start to end + 2~kb downstream) using FIMO,
these motifs showed a striking association with compensation status.
Among the 47 significantly sex-biased Z genes,
all 18 female-biased genes carried at least one motif in their gene body,
while all 9 motif-free sex-biased genes were male-biased
(Fisher exact $p = 0.007$; Figure~\ref{fig:motifs}).
The signal survived gene-size matching:
among genes below the Z median size (42~kb),
7 female-biased genes carried motifs
while all 9 motif-free genes were male-biased
(Fisher exact $p = 0.014$).

Two motifs showed individually detectable effects.
Motif A (\texttt{AGTTTTCCAAGSMGAM}; 16~nt) was present in 44 of 290 Z genes (15\%)
and showed the only individually significant sex-bias shift
(Mann-Whitney $p = 0.015$).
Motif B (\texttt{AAACATACCAGCAAGAAATAA}; 21~nt) was present in 148 Z genes (51\%)
and showed significant enrichment for female-biased genes (Fisher $p = 0.023$).
These two motifs co-occurred in 36 Z genes---significantly more
than expected by chance (Fisher OR = 5.50, $p < 0.0001$)---suggesting
they form a composite element with $\sim$41~bp spacing.
Genes carrying both motifs had the strongest female shift
(median log$_2$FC = $-$0.304 vs.\ $+$0.327 for genes with neither;
Kruskal-Wallis $p = 0.036$; Figure~\ref{fig:motif_cooccurrence}).

Critically, the sex-bias association is Z-specific.
The same motifs occur abundantly on autosomes,
but autosomal genes carrying Motif A and B together
show no sex-bias shift (median log$_2$FC = $+$0.355, Mann-Whitney $p = 0.72$
vs.\ autosomal non-carriers).
This mirrors the \textit{Drosophila} MRE system,
where MREs are present genome-wide
but the MSL complex acts exclusively on the X because
\textit{roX} lncRNAs target the complex in \textit{cis}.

TOMTOM comparison against JASPAR 2024 yielded suggestive matches:
Motif A to the Rel-homology domain family (NFAT/NF-$\kappa$B; best $E = 0.245$)
and Motif B to the Forkhead domain family (best $E = 0.218$).
The Forkhead match for Motif B is notable given that \textit{FoxA2},
a Z-linked Forkhead pioneer transcription factor capable of opening closed chromatin,
is among the conserved Z-linked genes.

To test whether LINE density modulates cis-element function,
we modeled sex-bias ($\log_2\text{FC}$) of 298 expressed Z genes
as a function of 10~kb flanking RTE-BovB density,
gene-body motif presence (Motif~A and Motif~B separately),
and their interactions.
The full model is significant ($R^2 = 0.055$, $F = 3.37$, $p = 0.006$)
and reveals a specific interaction:
RTE-BovB density drives male bias
(coefficient $= +0.039$~per~bp/kb, $p = 0.004$),
Motif~B presence promotes compensation
(coefficient $= +2.39$, $p = 0.043$),
but the RTE-BovB $\times$ Motif~B interaction is significantly negative
(coefficient $= -0.037$, $p = 0.007$),
meaning that Motif~B's compensating effect is abolished
in RTE-BovB-dense flanking chromatin.
In the highest RTE-BovB quartile, 100\% of significantly
sex-biased genes are male-biased (5M, 0F),
consistent with complete failure of motif-mediated compensation
above a density threshold.

\begin{figure}[tbp]
    \centering
    \includegraphics[width=\linewidth]{figures/genebody_motif_contingency_panels.png}
    \caption{Gene-body cis-elements predict compensation status on the Z chromosome.
    Among 47 significantly sex-biased Z genes,
    all 18 female-biased genes carry at least one Z-exclusive motif in their gene body,
    while all 9 motif-free sex-biased genes are male-biased (Fisher exact $p = 0.007$).
    The signal survives gene-size matching ($p = 0.014$ among genes $<$ 42~kb).
    These motifs parallel \textit{Drosophila} MSL Recognition Elements
    in their gene-body location, modest enrichment, and binary presence/absence effect.}
    \label{fig:motifs}
\end{figure}

\begin{figure}[tbp]
    \centering
    \includegraphics[width=0.9\linewidth]{figures/motif_cooccurrence_lfc.png}
    \caption{Co-occurrence of Motifs A and B predicts female-shifted expression on Z.
    Genes carrying both Motif A (\texttt{AGTTTTCCAAGSMGAM})
    and Motif B (\texttt{AAACATACCAGCAAGAAATAA})
    have median log$_2$FC = $-$0.304 (female-shifted),
    compared to $+$0.327 for genes with neither motif (Kruskal-Wallis $p = 0.036$).
    The two motifs co-occur more than expected by chance (Fisher OR = 5.50, $p < 0.0001$),
    suggesting a composite cis-regulatory element.
    The sex-bias association is Z-specific:
    autosomal genes carrying both motifs show no expression shift.}
    \label{fig:motif_cooccurrence}
\end{figure}

\subsection*{Z-linked long non-coding RNAs: \textit{Zmast} and \textit{Zfest}}

Among the most extreme sex-biased Z-linked transcripts in every species examined
are two long non-coding RNAs with opposing sex specificities.
The male-biased \textit{Zmast} (Z-male-specific transcript)
shows extraordinary expression asymmetry in \textit{O. bimaculoides}:
the core exon (4,685~bp) has a male-to-female ratio of 5,063:1
(Figure~\ref{fig:zmast_zfest}),
far exceeding the 2-fold difference expected from gene dosage alone.
This extreme male bias is conserved across all eight cephalopod species examined.

\paragraph{Eight-species annotation and alignment.}
We annotated \textit{Zmast} in all eight species
using a unified pipeline (see Methods).
Spliced transcript lengths range from 5,709~bp (\textit{O. bimaculoides})
to 9,576~bp (\textit{D. pealeii}),
with six of eight species in the 7--10~kb range
(Table~\ref{tab:zmast_annotation}).
\textit{O. bimaculoides} is the shortest because its 4,685~bp core exon
is highly repetitive and cannot be fully assembled from short reads alone;
\textit{S. officinalis} is limited by sequencing depth.
MAFFT L-INS-i alignment of all eight sequences
produced 14,726 alignment columns,
of which 8.6\% are gapless---consistent with the rapid sequence turnover
expected of lncRNA over $\sim$300~MY.

\begin{table}[tbp]
\centering
\small
\caption{\textit{Zmast} annotation across eight cephalopod species.}
\label{tab:zmast_annotation}
\begin{tabular}{llrrl}
\hline
Species & Spliced (bp) & Exons & Genomic span & Z chromosome \\
\hline
\textit{D. pealeii} & 9,576 & 2 & 9.6~kb & Dpe43 \\
\textit{O. vulgaris} & 9,067 & 2 & 18.0~kb & OX597829.1 \\
\textit{E. berryi} & 8,110 & 3 & 21.9~kb & Eber\_ch44 \\
\textit{S. oualaniensis} & 7,722 & 3 & 18.1~kb & HiC\_scaffold\_46 \\
\textit{E. scolopes} & 7,316 & 5 & 22.1~kb & CM044529.1\textsuperscript{a} \\
\textit{I. illecebrosus} & 7,168 & 2 & 15.3~kb & Z \\
\textit{S. officinalis} & 5,802 & 3 & 5.9~kb & OZ199189.1 \\
\textit{O. bimaculoides} & 5,709 & 6 & 169.1~kb & Z \\
\hline
\end{tabular}

\smallskip
\noindent\textsuperscript{a}\textit{Zmast} resides on a translocated Z fragment (CM044496.1/LG10) in \textit{E.\ scolopes}; the main Z chromosome is CM044529.1 (LG43).
\end{table}

\paragraph{Conserved architecture.}
Sliding-window conservation analysis (see Methods)
identified 19 conserved blocks spanning 20--1,570 alignment columns.
Nine blocks are present in all eight species,
seven in seven species, and three in six.
The alignment defines six structural domains
(Figure~\ref{fig:zmast_arch}):
a 5$^\prime$ Head (blocks B1--B2, 6/8 species),
Core A (B3--B9, all 8 species),
Core B (B10--B12, 6--8 species),
a Variable region (B13--B14),
a Poly(U) Sponge (B15, 7/8 species, absent from \textit{O. bimaculoides}),
and a 3$^\prime$ Tail (B16--B19).
Core A has the highest species occupancy
but, as described below, the lowest sequence-level conservation;
Core B comprises three small islands of 20--30~columns each,
separated from Core A by a $\sim$1,000-column gap
with 35--68\% occupancy.

\paragraph{Cross-species validation of RBP binding sites.}
FIMO scanning of each species' \textit{Zmast} against 258 metazoan RBP
position weight matrices (ATtRACT database, $p < 10^{-4}$; see Methods)
identified hundreds of binding sites per species
and suggested dense, multi-valent RBP binding across \textit{Zmast}.
However, systematic cross-species validation---mapping
all eight species' FIMO hits to alignment coordinates
and testing whether the same motif occurs at homologous positions---dramatically
revised this picture.
Of 132 motifs appearing in $\geq$4 species within the same block,
only 8 (6\%) are truly conserved
(non-low-complexity, at homologous alignment positions);
57 (43\%) are positionally conserved but driven
by species-specific homopolymers or dinucleotide repeats;
and 67 (51\%) occur at non-homologous positions
within the same block (independent hits).
In total, only $\sim$99~bp of specifically conserved sequence exists
across the entire 7,168~bp lncRNA (1.4\%).

The eight validated cross-species RBP binding sites
are concentrated in two domains.
Six reside in the 5$^\prime$ Head (blocks B1--B2):
RBM24 (6 species, identical \texttt{AGTGTGA}),
ROX8 (4 species, identical \texttt{CCATTTT}),
SNF (4 species), UAF-2 (4 species, identical \texttt{TTTTCAGG}),
MEX-5 (4 species, identical \texttt{TTAATAAT}),
and SRSF7 (5 species).
The remaining two are in the Poly(U) Sponge (block B15):
FBF-1 (5 species, identical \texttt{TGTGTATTTA})
and SYP (3 species).
Block~1 of the 5$^\prime$ Head is the deepest conservation hotspot in \textit{Zmast},
with 40\% perfectly conserved columns
and six non-homopolymer conserved sequence runs
including a 16~bp element (\texttt{TCTCTCACCCCATCCC})
shared across six species.
Notably, \textit{ROX8} is the \textit{Drosophila} ortholog of mammalian TIA-1,
which promotes U1~snRNP recruitment at the \textit{msl-2} 5$^\prime$ splice site
and is directly antagonized by Sex-lethal to repress dosage compensation in females,
and FBF-1 is a Pumilio-family translational repressor
with a central role in germline sex determination
in \textit{C. elegans}.

By contrast, Core A---despite containing all eight species---is
mostly AT-rich positional scaffolding.
Only block~B4 has genuinely conserved sequence runs;
blocks B5--B9 are spacers in which each species has evolved
its own low-complexity insertions
(poly-C in \illex, poly-T in \textit{D. pealeii},
(AT)$_n$ microsatellite in \textit{S. officinalis}, poly-A in \textit{O. vulgaris}).
Block~B14 (Variable region, 730~columns, all 8 species)
is the largest false positive:
all 18 ``shared'' FIMO motifs are at non-homologous positions.

\paragraph{Two sponge strategies.}
The dominant conserved feature of \textit{Zmast} is not specific RBP binding sites
but AT-richness (66--78\% across all blocks in all species)
and distributed poly(U) tracts.
However, poly(U) density differs markedly between lineages:
decapods average 57~bp/kb (range 52--79~bp/kb),
while octopods average 32~bp/kb
(\textit{O. bimaculoides}: 26.4~bp/kb, maximum tract only 8~nt;
every other species has tracts $\geq$14~nt).
ELAV/SXL binding to \textit{Zmast} is entirely poly(U)-mediated:
no ELAV-family position weight matrix produces
a sequence-specific hit on any \textit{Zmast} in any species.
Thus poly(U) tract density is the sole determinant of SXL sponge capacity.

Octopod \textit{Zmast} has compensated for poly(U) erosion
by enriching sequence-specific binding sites
for the RBM family (2.8$\times$ enrichment over decapods),
PCBP (1.9$\times$), and PUF (1.6$\times$).
This suggests a shift from a quantity-based sponge
(dense poly(U) tracts that bind all U-recognizing RBPs simultaneously)
to a quality-based sponge
(fewer but more diverse sequence-specific motifs).
Both lineages share CELF, PUF/FBF-1, and MEX-5 binding---the
irreducible regulatory core maintained for $\sim$300~MY
regardless of poly(U) content.

Stoichiometric analysis indicates that \textit{Zmast} poly(U) binding sites
are in 73--1,563-fold molar excess over RBP mRNAs in males,
though this captures only 2--18\% of total RBP protein---sufficient
to shift sensitive splice-site decisions past their inclusion thresholds
but not to globally deplete any single factor.
Importantly, \textit{MSL2} mRNA contains zero poly(U) tracts,
ruling out the \textit{Drosophila} \textit{Sxl}$\rightarrow$\textit{msl-2}
translational regulation mechanism.

\paragraph{\textit{Zfest}: eight-species annotation and alignment.}
The female-biased counterpart, \textit{Zfest} (Z-female-specific transcript),
is transcribed antisense to the \textit{Amnionless} (\textit{AMN}) gene
on the Z chromosome.
\textit{Zfest} expression is $\sim$18-fold higher in females than males in octopus,
while \textit{AMN} is consistently male-biased across all species examined---a
reciprocal sex-bias pair conserved over $\sim$300~MY.
We annotated \textit{Zfest} in all eight species
using a strandedness-aware scoring pipeline (see Methods).
Spliced transcript lengths range from 2,950~bp (\textit{E. scolopes})
to 19,114~bp (\illex)
(Table~\ref{tab:zfest_annotation}),
with far greater size variation than \textit{Zmast}
(6.5-fold range vs 1.7-fold).
\textit{O. vulgaris} is the only unspliced species (1 exon, 13.7~kb);
all others have 3--6 exons.
MAFFT L-INS-i alignment of all eight sequences
produced 26,275 alignment columns.
Critically, \textit{Zfest} remains Z-linked in all eight species,
including \textit{E. scolopes} where \textit{Zmast}
has translocated to an autosome (LG10)---making
\textit{Zfest} the more stably Z-linked of the two lncRNAs.

\begin{table}[tbp]
\centering
\small
\caption{\textit{Zfest} annotation across eight cephalopod species.}
\label{tab:zfest_annotation}
\begin{tabular}{llrrl}
\hline
Species & Spliced (bp) & Exons & Genomic span & Z chromosome \\
\hline
\textit{I. illecebrosus} & 19,114 & 6 & 21.7~kb & Z \\
\textit{O. vulgaris} & 13,712 & 1 & 13.7~kb & OX597829.1 \\
\textit{D. pealeii} & 12,109 & 5 & 17.6~kb & Dpe43 \\
\textit{O. bimaculoides} & 10,960 & 5 & 13.5~kb & Z \\
\textit{S. officinalis} & 8,499 & 6 & 21.9~kb & OZ199189.1 \\
\textit{S. oualaniensis} & 6,048 & 3 & 6.3~kb & HiC\_scaffold\_46 \\
\textit{E. berryi} & 3,749 & 5 & 10.8~kb & Eber\_ch44 \\
\textit{E. scolopes} & 2,950 & 3 & 5.6~kb & CM044529.1 \\
\hline
\end{tabular}
\end{table}

\paragraph{Conserved architecture and \textit{AMN} antisense overlap.}
Conservation analysis identified three conserved blocks
(Figure~\ref{fig:zfest_arch}),
compared to 19 in \textit{Zmast}---far fewer blocks
but each substantially larger.
Block~2 is the only block present in all eight species
(4,229 alignment columns);
Block~1 spans 1,829 columns in 7 species;
and Block~3 spans 175 columns in 5 species.
Together these contain $\sim$6,200 conserved columns,
comparable to \textit{Zmast} ($\sim$5,000)
despite the much greater size variation.

The defining structural feature of \textit{Zfest}
is its antisense overlap with \textit{AMN},
a single-pass transmembrane receptor required for
vitamin B12/cubilin endocytosis.
In all three species with genomic annotations
(\illex, \textit{S. oualaniensis}, \textit{O. bimaculoides}),
\textit{Zfest} overlaps the same five C-terminal \textit{AMN} coding exons
(CDS9--13, spanning the cytoplasmic tail through the extracellular domain).
The total exon-CDS overlap is remarkably consistent
across species ($\sim$730~bp across five exon pairs),
despite enormous variation in total \textit{Zfest} length.
Two of the three conserved blocks overlap \textit{AMN} coding exons:
Block~2 is antisense to CDS12--13
(transmembrane and cytoplasmic domains)
and Block~3 to CDS10--11 (extracellular stalk).
Only Block~1 is entirely \textit{AMN}-free,
providing the clearest evidence
for \textit{Zfest}-autonomous functional conservation
independent of \textit{AMN} protein-coding constraint.

\paragraph{Conserved RBP binding sites in \textit{Zfest}.}
Cross-species FIMO analysis with coordinate mapping
and low-complexity filtering (see Methods)
identified \textbf{26 truly conserved} RBP binding sites
across 20 unique RBPs---three times more
than \textit{Zmast}'s eight validated sites.
No site is conserved in all eight species;
the maximum is five (MATR3 and HOW),
reflecting a strong phylogenetic gradient:
decapods and sepiolids form a well-conserved core
(\textit{S. officinalis}: 21/26 sites,
\textit{E. scolopes}: 18,
\textit{E. berryi}: 17),
while octopods are much more divergent
(\textit{O. vulgaris}: 4,
\textit{O. bimaculoides}: 3,
\textit{D. pealeii}: 1).

The \textit{AMN}-free Block~1 carries five conserved sites:
MATR3 (5 species, identical \texttt{AATCTTG}),
a nuclear matrix protein involved in RNA retention;
PABPC3 and PABPC5 (4 species each);
RBM41 (4 species); and ROX8/TIA-1 (3 species).
Block~2 (\textit{AMN} antisense) carries 19 sites,
including the two most deeply conserved motifs:
HOW (5 species, identical 12-mer \texttt{AAATACTACCAC}),
a QKI/STAR-family muscle-specific splicing regulator;
and ETR-1 and CELF4 (4 species each),
CELF-family regulators of alternative splicing.
Germline-associated RBPs---PAPI (poly(A) polymerase interactor),
MEX-5 (P-granule partitioning),
and YBX2 (oocyte mRNA masking)---are
also concentrated in Block~2,
suggesting dual somatic and germline function.

\paragraph{\textit{Zfest} is structurally distinct from \textit{Zmast}.}
\textit{Zfest}'s RBP landscape differs qualitatively from \textit{Zmast}'s.
\textit{Zfest} is not poly(U)-enriched
and no ELAV/SXL-family position weight matrix produces
a sequence-specific hit on any \textit{Zfest} in any species---it
is not a sponge for SXL-family proteins.
MATR3 (nuclear matrix retention) is enriched 4.9-fold in \textit{Zfest}
relative to \textit{Zmast},
while ELAV/HuR, PTBP, and RBFOX are depleted ($<$0.5-fold).
\textit{Zfest} has equal or higher sequence-specific
RBP density than \textit{Zmast} (12.7/kb vs 11.9/kb in \illex),
while \textit{Zmast} has $\sim$1.5-fold higher total density
owing to its poly(U) tracts.
Where \textit{Zmast} operates as a multi-valent poly(U) sponge,
\textit{Zfest} appears to function through specific linear RBP recruitment---more
comparable to a regulatory scaffold than a molecular decoy.

Cephalopod MSL2 lacks the CXC domain
that recognizes \textit{roX} stem-loops in \textit{Drosophila},
ruling out a \textit{roX}-like chromatin-targeting mechanism for either lncRNA.
Infernal searches with \textit{roX1} and \textit{roX2} covariance models
produced no significant hits on any \textit{Zfest} sequence,
and full Rfam scanning (4,108 families) found zero structural matches,
confirming that \textit{Zfest} is not a \textit{roX}-like RNA.
\textit{Zfest} has one conserved stem-loop in Block~2
(R-scape covariation $E = 0.003$, 28 base pairs, 7-nt hairpin loop,
with a compensatory C--G to G--U change between decapods and octopods),
but this is a localized structural element
within a predominantly unstructured lncRNA.

The specific molecular targets of \textit{Zmast} and \textit{Zfest}
remain uncertain.
The poly(U) tracts and conserved RBP binding architecture of \textit{Zmast}
are consistent with a sponge model
in which male-produced \textit{Zmast} sequesters RNA-binding proteins
that would otherwise promote female-pattern gene regulation,
while \textit{Zfest}'s antisense overlap with \textit{AMN}
and nuclear-retention RBP signature (MATR3, HOW/QKI)
suggest a distinct role in regulating \textit{AMN} expression
or nuclear RNA metabolism on the female Z chromosome.
Definitive identification of the molecular functions
will require CLIP-seq or RNA pulldown experiments in cephalopod tissues.

\begin{figure}[tbp]
    \centering
    \includegraphics[width=\linewidth]{figures/zmast_architecture_figure.pdf}
    \caption{Conserved architecture of \textit{Zmast} across eight cephalopod species.
    \textbf{(A)} Alignment overview. Horizontal bars show per-species occupancy
    across 14,726 alignment columns (MAFFT L-INS-i).
    Conserved blocks (B1--B19) are colored by the number of species present
    (red = 8, orange = 7, yellow = 6).
    \textbf{(B)} Conservation score (weighted pairwise identity $\times$ occupancy,
    Gaussian-smoothed $\sigma$ = 100~columns) with conserved poly(U) hotspots
    (orange triangles; $\geq$5 consecutive columns with $\geq$4 species having U/T).
    \textbf{(C)} RBP motif density from FIMO scanning of the \illex sequence
    (ATtRACT metazoan database, $p < 10^{-4}$).
    Grey: poly(A/U)-driven hits; green: sequence-specific hits.
    Note that these densities are species-specific
    and do not represent cross-species conservation.
    \textbf{(D)} Six-domain architecture in \illex coordinates.
    Domain borders indicate conservation level:
    thick solid = genuinely conserved across species (5$^\prime$ Head, Poly(U) Sponge),
    dashed = positional scaffold (Core A, Core B, Variable),
    thin solid = mixed (3$^\prime$ Tail).
    Validated cross-species RBP binding sites are listed below each domain.}
    \label{fig:zmast_arch}
\end{figure}

\begin{figure}[tbp]
    \centering
    \includegraphics[width=\linewidth]{figures/zfest_architecture_figure.pdf}
    \caption{Conserved architecture of \textit{Zfest} across eight cephalopod species.
    \textbf{(A)} Alignment overview. Horizontal bars show per-species occupancy
    across 26,275 alignment columns (MAFFT L-INS-i).
    Three conserved blocks (B1--B3) are highlighted;
    color indicates the number of species present
    (red = 8, orange = 7, gold = 5).
    \textbf{(B)} Conservation score (weighted pairwise identity $\times$ occupancy,
    Gaussian-smoothed $\sigma$ = 100~columns) with truly conserved RBP binding sites
    marked as diamonds (colored by functional category).
    \textbf{(C)} RBP motif density from FIMO scanning of the \illex sequence
    (ATtRACt metazoan database, $p < 10^{-4}$),
    weighted by species occupancy to dampen single-species regions.
    Grey: poly(A/U)-driven hits; green: sequence-specific hits.
    Overlapping hits per RBP gene were merged.
    \textbf{(D)} Domain architecture in \illex coordinates.
    Domains are colored by \textit{Amnionless} (\textit{AMN}) antisense overlap:
    the 5$^\prime$ Head (blue) is \textit{AMN}-free,
    while subsequent domains correspond to \textit{AMN} CDS12--13
    (transmembrane/cytoplasmic),
    CDS10--11 (extracellular stalk),
    and CDS9 (extracellular domain).
    Validated conserved RBP binding sites are annotated per domain.
    \textbf{(E)} \textit{AMN} gene track showing coding exon positions
    in antisense orientation.}
    \label{fig:zfest_arch}
\end{figure}

\section*{Discussion}

Chromosome-wide dosage compensation has evolved repeatedly,
but the detailed mechanisms are known from only a handful of lineages---mammals
\citep{lyon1961gene},
\textit{Drosophila} \citep{lucchesi2015dosage, conrad2012dosage},
\textit{C. elegans} \citep{meyer2005x}---and
even in birds, the best-studied ZZ/ZW system,
the full extent of compensation has only recently been appreciated
\citep{papanicolaou2025multilayered, fallahshahroudi2025mirna}.
Here we show that cephalopods have independently evolved
dosage compensation approaching complete equalization
that has been conserved across 300~MY of coleoid divergence.
The underlying mechanism assembles
familiar components in a configuration seen nowhere else:
Z-enriched chromatin modifiers whose own compensation suggests
a self-reinforcing loop,
repetitive element architecture that opposes rather than assists upregulation,
gene-body cis-elements functionally analogous to \textit{Drosophila} MREs,
and deeply conserved sex-specific lncRNAs.

% \begin{figure}[H]
%     \centering
%     \includegraphics[width=\linewidth]{figures/dosage_compensation_model_figure.pdf}
%     \caption{Proposed model for cephalopod Z chromosome dosage compensation.
%     Z copy number sets \textit{Zmast} expression (males produce $\sim$5,000-fold more),
%     which sequesters U-rich RNA-binding proteins in males,
%     altering sex-specific splicing programs that control chromatin remodeling capacity.
%     Z-linked chromatin modifier genes maintain compensation via a self-reinforcing loop,
%     while gene-body cis-elements mark individual genes for compensation.
%     LINE elements (RTE-BovB) and CTCF insulators provide the structural architecture
%     for chromatin domain organization along the Z chromosome.}
%     \label{fig:model}
% \end{figure}

\subsection*{Dosage compensation is ancestral in coleoid cephalopods}

\textcite{papanicolaou2024z} provided the first evidence of dosage compensation
in any cephalopod, in two closely related octopus species separated by $\sim$2.5~MY.
Our results generalize this finding across seven species,
five orders, and over 300~MY of divergence
(Table~\ref{tab:cross_species_dc}).
That compensation is present in every species examined---and
at the mRNA level, more complete than in birds---demonstrates
that it is ancestral to the coleoid cephalopods,
not a lineage-specific innovation.

\subsection*{A self-reinforcing chromatin loop on the Z chromosome}

The most unexpected feature of the cephalopod Z chromosome
is the $\sim$10-fold enrichment of chromatin-modifying genes
among its conserved syntenic core (Table~\ref{tab:z_enrichment}).
These 16 genes span multiple independent pathways---histone methylation,
acetylation, ubiquitination, variant deposition, and heterochromatin regulation---yet
nearly all are themselves compensated (Table~\ref{tab:compass}).
This arrangement is consistent with a self-reinforcing feedback loop:
if these genes encode the epigenetic machinery
that maintains permissive chromatin on the Z,
their own compensation would sustain single-copy expression in females.
In this model, any erosion of compensation would diminish chromatin modifier output,
further eroding compensation---a destabilizing dynamic
that the compensated state must actively resist.
The architecture is reminiscent of autoregulatory chromatin loops
in developmental gene regulation,
but here it operates at the scale of an entire chromosome.

\subsection*{Chromatin remodelers as effectors of compensation}

The isoform-level regulation of chromatin modifiers
provides a quantitative association between epigenetic machinery
and chromosome-wide compensation output.
In \textit{O. bimaculoides}, sex-biased isoform proportions
of just four chromatin remodelers account for 57\% of Z:A variance
across 30 somatic samples---four-fold more than sex alone
(Figure~\ref{fig:za_regression}).
\textit{KDM3A}, an H3K9 demethylase, is the sole individually significant predictor,
suggesting a link between the removal of repressive chromatin marks
and the degree of Z chromosome upregulation in females.

Importantly, the specific downstream circuitry differs between species:
the \textit{SMARCC2}/\textit{CHD5} splicing circuit in \textit{O. bimaculoides}
is absent in \textit{S. oualaniensis},
where \textit{MSL2} is instead strongly female-biased.
What is conserved is not the wiring
but the upstream architecture:
Z-linked chromatin modifiers and sex-specific lncRNAs
in every species examined.

\subsection*{Lyon's repeat hypothesis, inverted}

\textcite{lyon1998x} proposed that LINEs on the mammalian X serve as
``way stations'' for the spreading of \textit{Xist}-mediated silencing---an
idea supported by the 2-fold LINE enrichment on the human X,
the tendency of LINE-poor regions to escape inactivation
\citep{bailey2000molecular},
and the active role of LINE-1 transcripts in nucleating
facultative heterochromatin during X-inactivation
\citep{chow2010line}.
On the mammalian X, LINEs and compensation are positively associated.
The cephalopod Z presents the mirror image:
the Z is LINE-enriched, yet LINE-dense regions are where genes
\textit{escape} compensation (Figure~\ref{fig:line_density}).

We propose that this inversion reflects a fundamental asymmetry
in compensation strategy.
Mammalian X-inactivation silences one of two copies;
LINEs propagate that silencing signal.
Cephalopod Z compensation upregulates a single copy---an
active chromatin-opening program that must overcome
TE-associated heterochromatin.
Host genomes constitutively repress TEs via H3K9me2/3,
and these repressive marks spread into flanking sequence
\citep{lee2017pervasive, huang2022species};
on the \textit{Drosophila} neo-X, TE-rich regions acquire
H3K9me2/3 heterochromatin that the MSL complex must overcome
\citep{zhou2013epigenome}.
In ZZ males, both copies experience this silencing equally---no
compensation is needed.
In ZO females, the compensation machinery must clear
LINE-associated heterochromatin to upregulate the single Z.
The male bias of genes in LINE-dense regions thus reflects
not passive hemizygosity but a failure of the female chromatin program
to penetrate TE-dense heterochromatin---a failure
quantified by the significant RTE-BovB $\times$ Motif~B interaction
in our regression model ($p = 0.007$),
which shows that cis-element function is abolished
in LINE-dense flanking chromatin.

\textit{KDM3A} provides a candidate enzymatic link.
As an H3K9 demethylase, it removes the exact mark deposited by TE silencing---and
is the sole significant predictor of Z:A compensation
in the isoform regression.
The resulting model is simple:
\textit{KDM3A}'s female isoform sets the demethylation rate;
flanking LINE density sets the H3K9me deposition rate.
Where demethylation wins, motifs are accessible and compensation succeeds.
Where deposition wins, motifs are buried and the gene escapes.
CTCF insulators, enriched on Z and positioned in gaps between LINE clusters,
may sharpen these domain boundaries.

This framework applies beyond cephalopods.
In \textit{Drosophila}, individual co-opted TEs nucleate compensation
by carrying MSL recognition elements to new locations
\citep{ellison2013dosage},
and in the guppy \textit{Poecilia picta},
a recent wave of piggyBac-like transposons carrying YY1 motifs
remodeled the X chromosome epigenome,
enabling male-specific hypomethylation and de novo dosage compensation
\citep{metzger2023transposon}.
In both cases, TEs serve as vehicles for cis-regulatory element dispersal
\citep{britten1971repetitive, lee2024causes}---but
the bulk TE landscape may simultaneously resist upregulation
through heterochromatin spreading.
Whether this negative association between bulk TE density and
compensation holds on the \textit{Drosophila} or \textit{Poecilia} X
remains to be tested.
Across cephalopods, the specific TE families turn over---RTE-BovB
in \textit{O. bimaculoides}, Penelope in \illex---but
the two-compartment architecture persists,
with chromatin modifiers uniformly in the euchromatic A compartment,
suggesting that TE-mediated domain organization on the Z
is under long-term selection.

\subsection*{Gene-body cis-elements as MRE analogs}

In \textit{Drosophila}, the MSL complex identifies X-linked targets
through gene-body Chromatin Entry Sites
containing a GA-rich MSL Recognition Element (MRE).
The composite cis-elements we identified on the cephalopod Z
are functionally analogous:
they reside in gene bodies, discriminate compensation status,
and---most strikingly---are present genome-wide
yet produce sex-biased expression only on the Z chromosome.
This mirrors the \textit{Drosophila} system precisely:
MREs are ubiquitous, but MSL targeting is X-specific
because \textit{roX} lncRNAs direct the complex in \textit{cis}.
If \textit{Zmast} or \textit{Zfest} serve an analogous addressing function,
the Z-restricted effect of otherwise ubiquitous motifs
is explained without invoking Z-specific sequence divergence.
That Motif~B matches a Forkhead domain consensus---and
\textit{FoxA2}, a Forkhead pioneer transcription factor,
is itself Z-linked and compensated---suggests
a candidate \textit{trans}-acting partner,
though experimental identification remains essential.

\subsection*{\textit{Zmast} and \textit{Zfest}: deeply conserved lncRNAs on the Z}

The conservation of \textit{Zmast} and \textit{Zfest}
across 300~MY of cephalopod divergence
is exceptional for non-coding RNA.
Fewer than 5\% of mammalian lncRNAs survive beyond $\sim$90~MY
\citep{hezroni2015principles},
and \textit{Xist} itself is a eutherian innovation ($\sim$160~MY).
The persistence of these lncRNAs implies strong purifying selection---but
the nature of that constraint is more nuanced
than simple sequence conservation would suggest.

For \textit{Zmast}, cross-species validation reveals
that only $\sim$99~bp of specifically conserved sequence exists
across 7,168~bp (1.4\%).
What is conserved is AT-richness, poly(U) content,
and the six-domain architecture---not individual binding sites,
87\% of which are low-complexity artifacts
or positionally non-homologous.
The conservation is at the level of biophysical properties
(multi-valent sponge capacity) rather than specific protein contacts,
analogous to binding-site turnover in promoter evolution.
The two cephalopod lineages illustrate this clearly:
decapods maintain dense poly(U) tracts (57~bp/kb)
that bind all U-recognizing RBPs simultaneously,
while octopods, having eroded poly(U) content to 32~bp/kb,
have compensated with 2--3-fold enrichment of sequence-specific sites
for RBM, PCBP, and PUF-family proteins.
Both strategies maintain the same irreducible core---CELF,
PUF/FBF-1, and MEX-5 binding---suggesting
that sponge capacity per se, not the specific mechanism,
is the target of selection.

The eight genuinely conserved RBP sites nonetheless point
to a suggestive mechanistic parallel with \textit{Drosophila}.
The \textit{Zmast} 5$^\prime$ Head contains conserved binding sites
for ROX8 (TIA-1), SNF (\textit{sans fille}), and UAF-2---the
three factors that in \textit{Drosophila} promote \textit{msl-2} intron splicing
and are antagonized by Sex-lethal to prevent dosage compensation
in females \citep{forch2001modulation, cline1999snf, merendino1999inhibition}.
All three sites coexist on a single $\sim$850~bp domain.
These predictions are entirely computational---no biochemical validation exists---but
their evolutionary conservation across 4--6 species
spanning $\sim$300~MY is the strongest available evidence
for functional significance.
The octopus lineage provides a natural test:
AT/TA repeat expansions have obliterated the ROX8, SNF, and UAF-2 sites
while preserving MEX-5, suggesting that germline-partitioning function
is under stronger constraint than the splicing-factor sites.

These observations motivate two non-exclusive functional models.
In a titration model, the massive male-specific expression
of \textit{Zmast} (5,000:1 M:F) creates a stoichiometric excess
of poly(U) binding sites (73--1,563-fold over cognate RBP mRNAs),
sequestering RNA-binding proteins away from endogenous targets in males.
Our stoichiometric analysis indicates that this excess
is sufficient to shift sensitive splice-site decisions
past their inclusion thresholds
without globally depleting any single factor---consistent
with the sex-biased isoform switching we observe
at Z-linked chromatin modifiers.
In a platform model, the 5$^\prime$ Head nucleates
a ribonucleoprotein complex---analogous to the \textit{roX} boxes
in \textit{Drosophila} or the A-repeat of \textit{Xist}---that
recognizes the gene-body cis-elements (Motifs A and B)
we identified on compensated Z genes.
Under either model, the lncRNAs provide chromosomal address
(\textit{cis}-linkage to Z),
the motifs provide gene-level specificity,
and the protein partners provide catalytic activity---explaining
why the motifs produce sex-biased effects only on Z
despite being present genome-wide.

\textit{Zfest} likely operates through a distinct mechanism.
Unlike \textit{Zmast}, it is not poly(U)-enriched
and carries no ELAV/SXL-family binding sites,
ruling out a sponge function.
Its enrichment for MATR3 (nuclear matrix retention)
and conserved antisense overlap with \textit{AMN}
instead suggest a role in nuclear RNA metabolism
or direct regulation of \textit{AMN} expression on the female Z---a
possibility underscored by the reciprocal sex-bias
of \textit{Zfest} (female) and \textit{AMN} (male)
conserved across $\sim$300~MY.

Identifying the protein partners of \textit{Zmast} and \textit{Zfest}---through
CLIP-seq, RNA pulldown, or proximity labeling---is
the single most important experiment for testing these models.
Equally important is determining which cell types express these lncRNAs
and when during development their expression initiates;
our bulk adult somatic RNA-seq cannot resolve these questions,
but single-cell and developmental time-course data
would clarify whether \textit{Zmast} and \textit{Zfest}
act globally or in specific tissue contexts.
Complementary approaches---CRISPR disruption of \textit{Zfest},
ATAC-seq to map sex-specific chromatin accessibility across the Z,
and ChIP-seq for H3K4me3 at Z-linked chromatin modifier loci---would
directly test the causal predictions of the feedback loop
and are now feasible in genetically tractable species such as \textit{E. berryi}.

\subsection*{Conclusions}

Dosage compensation in cephalopods represents an independent evolutionary solution
to a universal problem---balancing gene expression from an asymmetric sex chromosome.
That this solution converges on the same functional toolkit
used by mammals and \textit{Drosophila}---lncRNAs, chromatin modifiers,
and gene-body cis-elements---while assembling them
in a fundamentally different architecture
suggests that the design space for chromosome-wide expression regulation
is surprisingly constrained.
At the same time, the rapid turnover of downstream wiring
between cephalopod species, set against 300~MY of conserved upstream architecture,
illustrates how evolutionary stability can coexist with mechanistic flexibility.
Cephalopods thus offer not only a new model for sex chromosome biology
but a window into the general principles
by which genomes evolve chromosome-scale regulatory programs.

\section*{Materials and Methods}

\subsection*{Sample acquisition and tissue collection}

% TODO: Thomas to add specimen source details
Tissue was flash-frozen in liquid nitrogen and stored at $-80$~\textdegree C
during transport to the University of Oregon.
Mantle tissue (1~cm$^2$) was excised from one female individual (F24)
and one male individual (M71)
and delivered to the Genomics Core at the University of Oregon (GC3F) for sequencing.
All animals used in this study were handled under an approved IACUC protocol (AC-AABT8673),
including the use of deep anesthesia and procedures designed to minimize suffering.

\subsection*{Sequencing}

High-molecular-weight genomic DNA was extracted from 1~cm$^2$ mantle tissue
from two adult \illex, a female (F24) and a male (M71),
using Nanobind DNA disks following the animal tissue protocol.
DNA quantity and integrity were assessed using an Agilent Bioanalyzer.

HiFi (CCS) libraries were prepared and sequenced at the University of Oregon Genomics Core (GC3F).
Circular consensus reads were generated using PacBio ccs
and retained as HiFi reads at Q20 or higher.
The female sample (F24) was sequenced across five Sequel~II lanes and one Revio lane,
and the male sample (M71) was sequenced across two Revio lanes.

Hi-C libraries were prepared from F24 adult female brain tissue
by Phase Genomics using a multi-enzyme protocol
incorporating DpnII (GATC), DdeI (CTNAG), HinfI (GANTC), and MseI (TTAA),
and sequenced on one lane of an Illumina NovaSeq 6000 in SP mode at 150~bp paired-end.

RNA was isolated from the same F24 and M71 specimens
using the Zymo Research Direct-zol RNA Miniprep Kit (R2050).
Total RNA libraries were prepared using the Watchmaker RNA library preparation kit,
generating RF-oriented paired-end reads,
and sequenced on an Illumina NovaSeq 6000 at 150~bp paired-end.

Full-length cDNA libraries were prepared for F24
using the PacBio Kinnex protocol and sequenced on the Revio platform.
Raw reads were processed using the IsoSeq3 workflow as described below.

\subsection*{Genome assembly and scaffolding}

Genome size, heterozygosity, and repetitive content were estimated
using GenomeScope v2.0.1 \citep{Vurture2017}.
A subset of HiFi reads was used with KMC tools \citep{Kokot2017}
to generate a 21-mer frequency spectrum histogram,
which was analyzed with ploidy set to 2.

De novo genome assembly was performed using hifiasm v0.25.0 \citep{Cheng2021}
with Hi-C--assisted phasing.
Dual-haplotype self-scaffolding was enabled
and the haploid genome size prior (\texttt{--hg-size 3400m})
was provided based on GenomeScope results.
Assembly completeness was evaluated using compleasm v0.2.7 \citep{Huang2023Compleasm}
with the mollusca\_odb12 lineage.
Haplotypic duplications were identified and removed
using purge\_dups v1.2.6 \citep{Guan2020},
with read and assembly alignments generated by minimap2 v2.28 \citep{Li2018}.

Potential synthetic contamination was assessed using NCBI FCS-adaptor,
and biological contamination was assessed using FCS-GX v0.5.5.
Adaptors were removed by trimming,
and contigs outside the Mollusca lineage were removed from the assembly.

Hi-C reads were aligned to the draft assembly using bwa-mem2 v2.3,
and pairs were processed with pairtools v1.1.2 \citep{Open2C2024Pairtools}.
Scaffolding was performed using YaHS v1.2.2 \citep{Zhou2023YaHS}.
Contact maps were generated following the YaHS Juicer/JBAT workflow,
and assemblies were manually curated in Juicebox \citep{Durand2016Juicebox}.

Scaffolds were named in reference to \textit{Doryteuthis pealeii}
using BUSCO loci mapping from compleasm full\_table.txt outputs,
assigning scaffold names based on the highest number of shared BUSCO loci (minimum $>$100).

\subsection*{Genome annotation}

Adapter sequences and low-quality bases were trimmed from RNA-seq reads
using fastp v0.24.3 \citep{Chenfastp}.
Ribosomal RNA loci were identified using barrnap v0.9,
and residual rRNA-derived reads were removed by alignment
to a combined SILVA/barrnap reference
using bowtie2 v2.5.4 \citep{Langmeadbt2} in local high-sensitivity mode.

Non-rRNA reads were split into 20 files using SeqKit \citep{Shen2016SeqKit}
and aligned to the genome using HISAT2 v2.2.1 \citep{Kim2019HISAT2}
with options \texttt{--dta -q --max-intronlen 250000 --min-intronlen 20 --rna-strandness RF}
and GNU parallel \citep{Tange2018Parallel}.
Resulting BAM files were merged and sorted using samtools v1.21 \citep{Danecek2021SAMtools}.
Iso-Seq reads were processed using the IsoSeq3 pipeline (refine, cluster, collapse)
and aligned to the genome using minimap2 v2.29 \citep{Li2018}
in splice-aware high-quality mode
(\texttt{-uf -MD -ax splice:hq -G 1000000 --secondary=no}).

Gene models were generated using EviANN v2.0.5 \citep{Petersen2025EviAnn}
with short- and long-read alignments and protein homology evidence
from OrthoDB v12.0 (Lophotrochozoa clade combined with BUSCO proteins
from available cephalopod genomes).
The resulting GFF3 was standardized using AGAT v1.5.1 \citep{DainatAGATZenodo},
UTRs were added using PASA v2.5.3 \citep{Haas2003PASA},
and isoform structures were evaluated using SQANTI3 v5.5.4 \citep{Tardaguila2018SQANTI}.
The final annotation identified 23,192 protein-coding genes.

Long non-coding RNAs were identified using FEELnc v0.2.1 \citep{Wucher2017FEELnc}.
Transcript models were assembled using StringTie v3.0.0 \citep{Pertea2015StringTie}
from both Iso-Seq (\texttt{-conservative -L}) and RNA-seq alignments
and merged with the Iso-Seq GTF as guide annotation.

Functional annotation was performed using EnTAP v2.3.0 \citep{Hart2019EnTAP}
against eggNOG, NCBI nr, RefSeq Protein, Swiss-Prot, and InterProScan databases.
Protein sequences with high similarity to bacteria, fungi, viridiplantae, or archaea
were removed from the annotation.

\subsection*{Repeat, telomere, and centromere annotation}

Transposable elements were annotated using EarlGrey v7.0.2 \citep{Barber2025EarlGrey}
with a Mollusca-focused repeat library \citep{Smit2013RepeatMasker}.
Telomeric repeat arrays were identified using tidk v0.2.7 \citep{Brown2025tidk},
first in explore mode (\texttt{--minimum 5 --maximum 12})
then in search mode (\texttt{--string TTAGGG}).
Candidate centromere regions were inferred using RepeatObserver v1 \citep{Schultz2025RepeatOBserver}
by intersecting higher-order repeat intervals
with the lowest 2.5\% Hi-C read coverage across each chromosome.
Genome-wide tRNAs were annotated using tRNAscan-SE v2.0.12 \citep{Chan2021tRNAscanSE}
with default eukaryotic parameters,
retaining predictions with Infernal score $>$ 40.

\subsection*{Identification of the Z chromosome}

PacBio HiFi reads were aligned using minimap2 \citep{Li2018}
with the \texttt{map-hifi} preset.
Sex-specific coverage was summarized using mosdepth v0.3.3 \citep{Pedersen2018Mosdepth}
with options \texttt{-n --fast-mode --by 500000}.
Chromosomes were classified as candidate sex-linked
if they exhibited an approximately two-fold male-to-female depth ratio.

Synteny-based confirmation was performed using two complementary approaches.
First, GENESPACE v1.2.3 \citep{Lovell2022GENESPACE},
integrating OrthoFinder v2.5.5 \citep{Emms2019OrthoFinder}
and MCScanX \citep{Wang2012MCScanX},
was used to infer syntenic relationships among
\textit{O. bimaculoides}, \textit{O. sinensis}, \textit{E. scolopes},
\textit{S. esculenta}, \textit{D. pealeii},
and \textit{S. oualaniensis} \citep{Li2022FlyingSquid}.
Second, chromosome-level synteny was independently visualized
using ODP v0.3.3 \citep{Schultz2024ODPv033}.
The phylogenetic tree was generated using STAG \citep{Emms2019OrthoFinder}
with \textit{Nautilus pompilius} as outgroup.

\subsection*{Search for a W chromosome and pseudo-autosomal regions}

To evaluate whether W-linked sequence was present,
phased assembly outputs from hifiasm were examined for female-specific contigs.
All unplaced contigs were aligned to male-specific contigs using minimap2,
and contigs not aligning over 80\% of their length were considered candidates.
Candidate W-linked contigs were screened using sex-specific coverage differences
and by inspecting female-only presence across datasets.
Protein and nucleotide sequences of sex-linked genes reported in \textit{Nautilus}
\citep{Torrado2025NautilusSex} were used as queries
in DIAMOND similarity searches \citep{buchfink2015diamond} against the assembled genome.
Unplaced or low-confidence loci were additionally evaluated using transcript evidence
(RNA-seq and Iso-Seq assemblies) and screened for female-specific expression
indicative of W linkage.

Putative pseudo-autosomal regions were investigated by combining
sex-specific coverage and SNP density patterns along the Z chromosome.
Coverage was computed using mosdepth \citep{Pedersen2018Mosdepth},
and SNPs were called from long-read alignments
using Clair3 v1.2.0 \citep{Chang2024Clair3}.
Regions were considered PAR candidates if they showed
near-equal male and female coverage
and elevated heterozygosity consistent with ongoing recombination.

\subsection*{Transcript quantification}

Transcript sequences were extracted using gffread from final annotations
and corresponding genome FASTAs.
Transcript abundance was quantified using Salmon v1.10.3 \citep{patro2017salmon}
with a decoy-aware index constructed by concatenating
the transcriptome FASTA with the genome FASTA.
Quantification used selective alignment with
\texttt{--validateMappings --gcBias --seqBias}.
For \textit{S. oualaniensis}, a v4 annotation was generated by merging
cross-species Iso-Seq (from \illex collapsed transcripts),
RNA-seq StringTie assemblies, and LiftOn-guided gene models
using StringTie --merge, producing 53,703 transcripts.
Expression values are reported as Transcripts Per Million (TPM).

\subsection*{Illex differential expression and normalization}

Differential expression testing for \illex was performed
using DESeq2 v1.50.2 \citep{Love2014DESeq2}.
Counts were normalized and sex-biased expression was assessed
by comparing male and female expression profiles.
Effect sizes were summarized as $\log_2$(M/F) using pseudocount stabilization.

\subsection*{Dosage compensation analysis}

Z:A ratios were calculated for each sample as:
$\text{Z:A} = \bar{x}_Z / \bar{x}_A$,
where $\bar{x}_Z$ and $\bar{x}_A$ are the mean TPM values
for Z-linked and autosomal genes, respectively,
after filtering to genes with TPM $> 1$
and excluding the top 5\% of expressed genes per chromosome class,
following the method of \textcite{papanicolaou2024z}.
Z-linked genes were identified in each species
via synteny with the \illex Z chromosome.
M:F expression ratios were calculated per gene
and compared between Z and autosomal genes
using the Wilcoxon rank-sum test.

\subsection*{Differential expression analysis}

Differential expression was performed using PyDESeq2 v0.4.8 \citep{muzellec2023pydeseq2}
for each species.
For cross-species comparisons,
a uniform model ($\sim$sex) was applied to somatic samples
with tissue covariates excluded (\texttt{--no-tissue})
to ensure comparability.
Temperature-stressed samples from \textit{O. bimaculoides}
(PRJNA948369 warm/cold conditions; PRJNA1199547 warm/cold)
were excluded prior to analysis.
Multiple testing correction used the Benjamini-Hochberg procedure (FDR $<$ 0.05).

\subsection*{Cross-species orthology mapping}

Reciprocal best-hit (RBH) orthology was established between \illex (16,656 genes)
and each target species using DIAMOND BLASTP \citep{buchfink2015diamond}
(v2.1.9; \texttt{--very-sensitive},
E-value $< 10^{-10}$, $>30\%$ identity, $>50\%$ query coverage).
Z-linked gene content was compared across four species
(\illex, \textit{O. bimaculoides}, \textit{S. oualaniensis}, \textit{O. vulgaris}).
Gene functional categories were assigned using a keyword-based classifier
encompassing COMPASS/MLL, HAT, HDAC, CHD, EZH, KDM, Reader, H2Bub,
Chaperone, THO, and PAF families.

\subsection*{RepeatMasker and LINE density analysis}

Repeat annotations for \textit{O. bimaculoides}
were obtained from RepeatMasker output for assembly GCA\_046866305.1.
LINE density was calculated per chromosome as total LINE bp per Mb of chromosome length.
Per-gene LINE density was calculated in 10~kb flanking windows
centered on each expressed Z gene.
RTE-BovB and Dong-R4 densities were computed separately
and compared across 500~kb sliding windows along the Z chromosome.

\subsection*{CTCF motif scanning}

The Z chromosome and all autosomes were scanned for CTCF binding motifs
using FIMO (MEME Suite v5.5.9) with the JASPAR MA0139.2 position weight matrix
at $p < 10^{-4}$.
Motifs were classified as genic (overlapping merged gene intervals)
or intergenic.
Chromosome-level CTCF density was modeled as a function of GC content
and chromosome size (linear regression, $R^2 = 0.794$),
and Z chromosome residuals were compared to the autosomal distribution.

\subsection*{Hi-C compartment analysis}

Hi-C data from a female \illex individual were processed
at 50~kb resolution from a .mcool file.
A/B compartments were identified by eigenvector decomposition of the observed/expected
contact matrix, with the PC1 sign oriented so that positive values
correspond to the A (gene-dense, euchromatic) compartment.
Compartment strength was measured as PC1 standard deviation per chromosome,
and TE family associations with compartment identity were assessed
by Spearman correlation of TE density with PC1 scores across 50~kb bins.

\subsection*{Discriminative motif discovery}

STREME (MEME Suite v5.5.9) was used for discriminative motif discovery
with the following comparisons:
compensated vs.\ male-biased Z gene flanks (10~kb) and promoters (5~kb upstream),
and compensated vs.\ female-biased Z gene flanks and promoters.
Parameters: motif width 6--18~bp, $p < 0.05$, both strands, maximum 5 motifs per comparison.
Discovered motifs were compared to the JASPAR 2024 CORE non-redundant database
using TOMTOM ($E < 0.1$).
RepeatMasker overlap was computed for each motif site.

\subsection*{Gene-body motif scanning}

The 10 Z-exclusive motifs from STREME analysis
were scanned across gene bodies (gene start to gene end + 2~kb downstream)
of 290 expressed Z genes and 1,180 expression-matched autosomal controls
using FIMO ($p < 10^{-6}$ for presence/absence classification,
$p < 10^{-4}$ for density analysis).
Five tandem-duplicated genes at the Z terminus were excluded.
Gene-size confounding was assessed by restricting analysis
to genes below the Z median size (42~kb)
and by partial correlation of motif density with log$_2$FC controlling for gene size.

\subsection*{LINE density and motif interaction model}

To test whether flanking LINE density modulates gene-body motif efficacy,
we modeled sex bias ($\log_2\text{FC}$ from DESeq2)
of 298 expressed Z-linked genes as a multiple linear regression.
Per-gene RTE-BovB density was calculated in 10~kb flanking windows
(5~kb upstream + 5~kb downstream of gene boundaries)
from RepeatMasker annotations, expressed as bp/kb.
Gene-body motif presence was determined from FIMO scans
($p < 10^{-6}$) for Motif~A (STREME-15) and Motif~B (STREME-20) separately.
The full model was $\text{lfc} \sim \beta_0 + \beta_1 \cdot \text{BovB}
+ \beta_2 \cdot \text{MotifA} + \beta_3 \cdot \text{MotifB}
+ \beta_4 \cdot \text{BovB} \times \text{MotifA}
+ \beta_5 \cdot \text{BovB} \times \text{MotifB}$,
fitted by ordinary least squares using statsmodels v0.14.
For genes with multiple transcripts in the DESeq2 output,
the transcript with the highest baseMean was retained.
Significance of individual coefficients was assessed by $t$-test;
overall model significance was assessed by $F$-test.

\subsection*{Alternative splicing analysis}

Sex-biased alternative splicing was assessed using SUPPA2.
Transcript-level TPM values from Salmon were used to calculate
percent spliced-in (PSI) values for local alternative splicing events
(skipped exon, alternative first/last exon, alternative 5$'$/3$'$ splice site,
mutually exclusive exon, retained intron).
Differential splicing between sexes was tested using SUPPA2 diffSplice
with 1,000 bootstrap iterations ($p_\text{adj} < 0.05$, $|\Delta\text{PSI}| > 0.1$).

\subsection*{\textit{Zmast} and \textit{Zfest} identification}

Candidate loci corresponding to \textit{Zmast} and \textit{Zfest}
were identified by sequence homology searches
using reference sequences from \textcite{papanicolaou2024z}.
Similarity searches were performed using DIAMOND \citep{buchfink2015diamond}
in sensitive mode against the \illex genome assembly,
predicted proteins, and curated transcript models.
To evaluate whether candidate \textit{Zfest}-like sequences
corresponded to full-length antisense lncRNAs,
candidate hits were cross-referenced with Iso-Seq--derived transcript models
and SQANTI3 classifications \citep{Tardaguila2018SQANTI}
to confirm transcript completeness, splice-junction support,
and strand orientation relative to the \textit{Amnionless} gene.
For each locus, the genomic neighborhood was characterized
by identifying flanking genes within $\pm$5~Mb
and assessing whether conserved neighbors matched
the expected syntenic context.

\subsection*{\textit{Zmast} annotation and alignment}

\textit{Zmast} was annotated in all eight species using a three-step pipeline.
First, the syntenic region containing \textit{Zmast} was extracted
from each genome assembly by identifying the conserved flanking genes
UBL3 and CG31559 via the cross-species orthology table,
retrieving their genomic coordinates from the species annotation,
and extracting the intervening region with 10~kb padding
using \texttt{samtools faidx}.
Second, for species with available male RNA-seq data
(\textit{O. bimaculoides}, \textit{S. officinalis}, \textit{S. oualaniensis}),
reads were extracted from existing BAM files for the syntenic region,
male samples were merged, and StringTie v2.2.1 was run
in de novo mode (\texttt{-m 100 -f 0.01})
to assemble transcripts.
Female BAMs were extracted in parallel for coverage comparison.
Third, candidate exons were identified by BLAT alignment
of both the species' existing \textit{Zmast} reference
and the \illex \textit{Zmast} as a universal query.
Candidates were scored by BLAT hit quality,
male:female coverage ratio, and poly(U) tract density
($\geq$5~nt tracts per kb).
For species without BAMs (\illex, \textit{D. pealeii},
\textit{E. scolopes}, \textit{E. berryi}, \textit{O. vulgaris}),
the pipeline performed syntenic region extraction
and BLAT scoring against existing reference sequences.
Each species' canonical \textit{Zmast} sequence was frozen
at \texttt{\{species\}/reference/Zmast\_current.fa}.

Multiple sequence alignment was performed with MAFFT v7.520
using the L-INS-i algorithm
(\texttt{--localpair --maxiterate 1000}),
producing a 14,726-column alignment.
Conserved blocks were identified by sliding-window analysis
(window = 50 alignment columns, step = 10 columns)
scoring each window by a combined metric of
occupancy (fraction of species with non-gap sequence)
and pairwise identity (fraction of identical nucleotides
among non-gap species, weighted by $n_\text{species}/8$).
Windows with combined score $\geq$3.0 were retained
and merged into 19 conserved blocks.

\subsection*{\textit{Zfest} annotation and alignment}

\textit{Zfest} was annotated in all eight species
using a strandedness-aware scoring pipeline.
Candidate exons in the syntenic \textit{AMN} region
were identified by BLAT alignment of
the \illex \textit{Zfest} as a universal query
and scored by alignment quality, antisense orientation
relative to \textit{AMN},
and female:male coverage ratio where RNA-seq data were available.
Sense-strand contamination was excluded
by requiring candidate exons to be transcribed antisense to \textit{AMN}.
Each species' canonical \textit{Zfest} sequence was frozen
at \texttt{results/zfest\_multispecies.fa}.
Multiple sequence alignment was performed with MAFFT v7.520
using the L-INS-i algorithm,
producing a 26,275-column alignment.
Conserved blocks were identified using the same sliding-window approach
as for \textit{Zmast}
(Gaussian smoothing $\sigma$ = 100~columns,
combined occupancy-identity score $\geq$0.30,
minimum 20 columns,
minimum 4-species mean occupancy),
yielding three conserved blocks.

\subsection*{RBP motif scanning and cross-species validation}

RNA-binding protein motif scanning was performed with FIMO
(MEME Suite v5.5.9, \texttt{--thresh 1e-4 --text --norc})
using 258 metazoan RBP position weight matrices
from the ATtRACT database \citep{giudice2016attract}.
Each species' \textit{Zmast} and \textit{Zfest} sequences were scanned independently.
FIMO hits were classified as poly(A/U)-driven
(matched sequence $\geq$80\% A or T/U)
or sequence-specific.

Cross-species validation of RBP binding sites
used the eight-species alignment to test positional conservation.
For each species, a coordinate map was constructed
relating ungapped sequence positions to alignment columns.
FIMO hits were then projected onto alignment coordinates.
A motif was classified as ``positionally conserved''
if hits in $\geq$3 species fell within 50 alignment columns
of one another,
and as ``truly conserved'' if additionally
the matched sequences were non-low-complexity
(not homopolymer- or dinucleotide-repeat-driven).
Motifs appearing in $\geq$4 species within the same conserved block
but at non-overlapping alignment positions
were classified as ``scattered'' (independent hits).

Poly(U) tract conservation was assessed
by identifying all poly(U/T) tracts $\geq$5~nt in each species,
mapping them to alignment coordinates,
and scoring positional overlap ($\pm$10 columns) across species.
Conserved hotspots were defined as $\geq$5 consecutive alignment columns
in which $\geq$4 species have a T/U residue.

Conserved sequence runs were identified as stretches of $\geq$6
consecutive alignment columns in which $\geq$6 species
share identical nucleotides,
excluding pure homopolymer runs (e.g., AAAAAA, TTTTTT).

\subsection*{DNA methylation analysis}

HiFi reads containing base-modification tags were aligned to the final genome assembly
using pbmm2 v1.17.0.
CpG methylation probabilities were computed from aligned HiFi BAMs
using pb-CpG-tools,
which derives site-level 5mC probabilities from MM and ML auxiliary tags.
Site-level methylation outputs were generated for whole-genome visualization
and for annotation-aware summaries.
For gene- and region-level summaries,
CpG-site probabilities were aggregated within each feature
to generate a mean methylation probability per feature.
Features were classified as highly methylated
if the aggregated probability exceeded 0.80.
Differential methylation was calculated per chromosome as the mean
male minus female methylation difference across gene bodies.
The Z chromosome was compared to autosomal chromosomes
using the Wilcoxon rank-sum test.

\subsection*{Data and code availability}

% TODO: fill in accession numbers and repository URLs
The genome assembly, annotation, transcriptome FASTA, predicted proteins,
and auxiliary tracks (TE annotation, telomere calls, centromere intervals, methylation tracks)
will be deposited at NCBI under BioProject accession [TBD].
Raw sequencing data (PacBio HiFi, Hi-C, stranded total RNA-seq, and Iso-Seq)
will be deposited in NCBI SRA.
Public datasets used for comparative analyses are listed with accessions
in Table~\ref{tab:dc_samples}.
Workflow scripts and custom notebooks used for figure generation
will be deposited at [TBD].

\section*{Acknowledgements}

\printbibliography
\clearpage
\beginsupplement

\section*{Supplemental Results}

\phantomsection
\subsection*{Genome assembly and annotation of \illex}
\label{sec:supp_assembly}

We generated a high-quality chromosome-scale genome assembly for a female \illex individual
using a combination of PacBio HiFi long reads, Illumina short reads, and Hi-C chromatin conformation capture data (Fig.~\ref{fig:phylo}A).
The final assembly spanned 4.14 Gb with a contig N50 of 34.71 Mb and scaffold N50 of 75.92 Mb (Table~\ref{tab:assembly_stats}).
BUSCO analysis indicated 95.6\% completeness against the metazoan gene set,
demonstrating high assembly quality.
Hi-C scaffolding resolved 46 chromosomes, consistent with the decapodiform karyotype of 2n = 92 \citep{gao1990karyological}.
GenomeScope2 analysis of the female HiFi k-mer spectrum (k = 21)
estimated a haploid genome size of 2.45~Gb (79.8\% unique sequence),
heterozygosity of 1.06\%, and a duplication rate of 67.3\%
(Figure~\ref{fig:genomescope}).
The bimodal k-mer distribution is consistent with high heterozygosity
and extensive repetitive content typical of cephalopod genomes.

\begin{figure}[H]
    \centering
    \includegraphics[width=0.8\linewidth]{figures/genomescope.png}
    \caption{GenomeScope2 k-mer profile of the female \illex{} HiFi reads (k = 21).
    The first peak ($\sim$17$\times$) corresponds to heterozygous k-mers
    and the second peak ($\sim$33$\times$) to homozygous k-mers.
    Estimated haploid genome size: 2.45~Gb; heterozygosity: 1.06\%;
    unique sequence: 79.8\%.}
    \label{fig:genomescope}
\end{figure}

To create a comprehensive gene annotation, we generated RNA-seq data from multiple tissues
using both Illumina short-read and PacBio Iso-Seq long-read technologies.
Annotation using RNA-seq data from multiple tissues identified 23,192 protein-coding genes,
with a mean intron length of 38 kb.
BUSCO analysis of the protein set indicated 93.1\% completeness (83.2\% single-copy, 6.7\% duplicated),
confirming the high quality of the gene annotation.

Repeat annotation revealed that 58.2\% of the \illex genome consists of repetitive elements (Table~\ref{tab:repeat_stats}).
The repeat landscape is dominated by unclassified repeats (29.1\% of the genome),
followed by LINE elements (12.6\%), DNA transposons (6.3\%), and low-complexity regions (5.1\%).
Penelope elements, a clade of retrotransposons common in invertebrates, comprised 3.3\% of the genome.
LTR retrotransposons (2.5\%), SINE elements (1.7\%), and rolling-circle transposons (0.6\%)
were also identified.
In total, we annotated 1,889 distinct classifications of unclassified repeats and
over 7.8 million individual repeat instances across all categories,
reflecting the complex repetitive architecture typical of cephalopod genomes.

\phantomsection
\subsection*{Z chromosome hypermethylation in females}
\label{sec:supp_methylation}

The main text reports that the Z chromosome is hypermethylated in females
relative to males, based on differential methylation across chromosomes
(Figure~\ref{fig:methylation}).
An alternative view directly comparing CpG methylation fractions
between autosomes and the Z within each sex confirms this pattern
(Figure~\ref{fig:meth_auto_vs_z}).
In both sexes, Z-linked genes are significantly more methylated than autosomal genes
(male: $p = 4.23 \times 10^{-10}$; female: $p = 1.48 \times 10^{-20}$,
Wilcoxon rank-sum tests),
but the effect is substantially larger in females (median Z methylation 0.41 vs.\ 0.19 for autosomes)
than in males (0.31 vs.\ 0.18),
consistent with the female-specific upregulation of the hemizygous Z
being accompanied by elevated gene-body methylation.

\begin{figure}[H]
    \centering
    \includegraphics[width=0.6\linewidth]{figures/methylation_auto_vs_z.png}
    \caption{CpG methylation of autosomal versus Z-linked genes in male and female \illex.
    Z-linked genes (red) are significantly more methylated than autosomal genes (blue) in both sexes,
    but the difference is more pronounced in females ($p = 1.48 \times 10^{-20}$)
    than in males ($p = 4.23 \times 10^{-10}$),
    consistent with gene-body methylation accompanying transcriptional upregulation
    of the hemizygous female Z.}
    \label{fig:meth_auto_vs_z}
\end{figure}

\phantomsection
\subsection*{Sex-biased chromatin modifier expression in \textit{O. bimaculoides}}
\label{sec:supp_chromatin}

Genome-wide differential expression analysis in \textit{O. bimaculoides}
(30 somatic samples, 15F:15M; PyDESeq2, $p_\text{adj} < 0.05$, $|\text{log}_2\text{FC}| > 1$)
identified 2,470 genes with significant sex-biased expression:
1,454 male-biased and 1,016 female-biased.
On the Z chromosome, 63 genes were significantly sex-biased (35 male, 28 female),
with a median log$_2$FC of +0.002,
confirming effective compensation at the chromosome level.

Among chromatin modifier genes specifically,
27 isoforms showed significant sex-biased expression,
15 female-biased and 12 male-biased
(Tables~\ref{tab:chromatin-female} and \ref{tab:chromatin-male}).
These encompass multiple chromatin-modifying families:
SWI/SNF remodeling complex subunits (\textit{PBRM1}, \textit{SMARCC2}, \textit{SMARCD1}, \textit{ARID2}),
lysine demethylases (\textit{KDM5A}, \textit{KDM3A}),
chromodomain helicases (\textit{CHD5}, \textit{CHD1L}),
histone acetyltransferases (\textit{KAT7}, \textit{KAT6B}),
heterochromatin readers (\textit{CBX5}/HP1$\alpha$, \textit{CDYL}, \textit{SMCHD1}),
WRAD scaffold members (\textit{RBBP5}, \textit{KMT2E}),
and bromodomain readers (\textit{BAZ2B}, \textit{BRD7}).
Three Z-linked genes show sex bias:
the SWR1 remodeler \textit{Domino/SRCAP} exhibits isoform switching
(isoform .3 female-biased, isoform .1 male-biased),
and \textit{THOC6} is male-biased.

\input{tables/chromatin_sex_biased_tables}

\begin{figure}[H]
    \centering
    \includegraphics[width=0.9\linewidth]{figures/rte_bovb_quartile_sex_bias.png}
    \caption{RTE-BovB density quartile analysis on the Z chromosome.
    Among significantly sex-biased Z genes,
    the highest RTE-BovB density quartile (Q4) shows 89\% male-biased genes (8M:1F),
    while lower quartiles show no consistent pattern.
    The effect is categorical rather than graded:
    female-biased genes are excluded from the highest-density RTE-BovB chromatin.}
    \label{fig:rte_quartile}
\end{figure}

\phantomsection
\subsection*{CTCF insulation and spatial clustering on the Z chromosome}
\label{sec:supp_ctcf}

Although the Z chromosome carries the highest CTCF density of any chromosome
(main text),
this enrichment is diffuse rather than boundary-specific.
If CTCF insulators specifically demarcated compensation domains,
intergenic regions at transitions between compensated and sex-biased genes
should have higher CTCF density than regions between genes of the same status.
This is not observed:
median CTCF density is 27.6/Mb between pairs of compensated genes,
23.2/Mb at compensated--biased boundaries,
and 38.8/Mb between pairs of sex-biased genes
(Mann-Whitney $p = 0.128$ for compensated--biased vs.\ compensated--compensated boundaries).

Nevertheless, CTCF density modulates co-regulation of adjacent gene pairs.
Adjacent expressed Z genes with no intervening CTCF sites
show significantly correlated sex-bias
($r = +0.250$, $p = 0.015$, $n = 94$ pairs),
while pairs separated by $\geq$3 CTCF sites show no correlation
($r = -0.066$, NS; Fisher $z = 2.25$, $p = 0.025$).
However, this effect is confounded with intergenic distance
(partial correlation controlling for gap size: $r = -0.019$, $p = 0.77$),
suggesting that closely spaced genes share regulatory environments
and CTCF density co-varies with physical separation.

Consistent with this compartment-level organization,
sex-biased and compensated genes are non-randomly distributed along Z.
A Wald-Wolfowitz runs test shows significant spatial clustering
of genes with the same compensation status
(observed 75 runs vs.\ 91.2 expected; $z = -2.92$, $p = 0.004$).
Lag-1 spatial autocorrelation of log$_2$FC is significant
($r = +0.132$, $p = 0.034$)
but decays by lag-2 ($r = +0.042$, NS),
indicating that co-regulation extends only to immediate neighbors.

\begin{figure}[H]
    \centering
    \includegraphics[width=\linewidth]{figures/illex_chrom_seq_depth.png}
    \caption{PacBio HiFi read depth across \illex chromosomes in female and male individuals.
    Normalized read depth is shown for each chromosome.
    Autosomes show comparable depth between sexes,
    while the Z chromosome (chr42) shows reduced coverage in the female
    relative to the male, consistent with hemizygosity in the ZO female.}
    \label{fig:seq_depth}
\end{figure}

\begin{table}[H]
\centering
\caption{Genome assembly and annotation statistics for \textit{Illex illecebrosus}}
\label{tab:assembly_stats}
\begin{tabular}{llr}
\toprule
\textbf{Category} & \textbf{Metric} & \textbf{Value} \\
\midrule
\multirow{5}{*}{Genome size} & Total assembly length (bp) & 4,143,127,294 \\
 & Chromosome & 46 \\
 & Chromosome assembly & 3,468,808,377 \\
 & Unplaced scaffolds & 3,755 \\
 & Unplaced scaffolds (Mb) & 674,318,917 \\
\midrule
\multirow{4}{*}{Contiguity} & Contig number & 179 \\
 & Contig N50 (Mb) & 34.71 \\
 & Scaffold N50 (Mb) & 75.921 \\
 & Largest scaffold (Mb) & 119.467 \\
\midrule
\multirow{4}{*}{Completeness mollusca} & BUSCO complete (\%) & 95.57 \\
 & BUSCO single-copy (\%) & 93.01 \\
 & BUSCO duplicated (\%) & 0.48 \\
 & BUSCO missing (\%) & 2.08 \\
\midrule
\multirow{4}{*}{Accuracy} & QV (Merqury) & 36.5 \\
 & Mapping rate (\%) & 99.7 \\
 & Split-read rate (\%) & 12.54 \\
 & Depth & 37.442 \\
\midrule
Repetitive content & Repeats (\%) & 58 \\
\midrule
\multirow{2}{*}{Annotation} & Protein-coding genes & 23,192 \\
 & Mean intron length (kb) & 38 \\
\midrule
\multirow{4}{*}{Completeness (protein)} & BUSCO complete (\%) & 93.12 \\
 & BUSCO single-copy (\%) & 83.19 \\
 & BUSCO duplicated (\%) & 6.65 \\
 & BUSCO fragmented (\%) & 3.28 \\
\bottomrule
\end{tabular}
\end{table}

\begin{table}[H]
\centering
\caption{Repeat annotation statistics for the \textit{Illex illecebrosus} genome}
\label{tab:repeat_stats}
\begin{tabular}{lrrrr}
\toprule
\textbf{Classification} & \textbf{Coverage (bp)} & \textbf{Count} & \textbf{Proportion} & \textbf{Distinct Classifications} \\
\midrule
DNA & 247,981,481 & 779,103 & 0.063 & 242 \\
LINE & 437,157,312 & 1,188,842 & 0.126 & 230 \\
LTR & 84,997,764 & 171,363 & 0.025 & 81 \\
Low-complexity & 176,258,651 & 1,527,910 & 0.051 & 1,726 \\
Penelope & 113,970,328 & 311,297 & 0.033 & 121 \\
Rolling-circle & 20,935,150 & 86,259 & 0.006 & 9 \\
SINE & 59,793,525 & 164,976 & 0.017 & 29 \\
Unclassified & 1,009,867,061 & 3,351,152 & 0.291 & 1,889 \\
\midrule
\textbf{Total} & & & \textbf{0.582} & \\
\bottomrule
\end{tabular}
\end{table}

\begin{table}[H]
\centering
\caption{Z-linked chromatin and transcription factor genes conserved across cephalopod species.
These 22 genes (16 chromatin, 6 transcription factor) were identified among 188 conserved Z-linked genes
shared by four species (\textit{O.~bimaculoides}, \textit{I.~illecebrosus}, \textit{S.~oualaniensis},
\textit{O.~vulgaris}). The ``Spp'' column indicates the number of species (out of eight) in which each gene
remains Z-linked at increasing phylogenetic stringency. Thirteen chromatin genes and four transcription
factors survive at the most stringent 8-species threshold (132 conserved genes).
Chromatin genes are enriched 2.8-fold on Z ($p < 0.001$); transcription factors 33-fold ($p < 10^{-6}$).
Ortholog identifiers for all eight species are provided in Supplementary Table~S1.}
\label{tab:conserved_z_genes}
\begin{tabular}{lllll}
\toprule
\textbf{Gene} & \textbf{\textit{O. bimac} ID} & \textbf{Category} & \textbf{\textit{Illex} LOC} & \textbf{Spp} \\
\midrule
\multicolumn{5}{l}{\textit{Chromatin/epigenetic (16)}} \\
\addlinespace[2pt]
\textit{ING4/ING5}    & obimac\_0008617 & HBO1/MOZ PHD reader        & LOC\_00003818 & 8 \\
\textit{UBN2/UBN1}    & obimac\_0008620 & HIRA chaperone            & LOC\_00003702 & 8 \\
\textit{SET1B}        & obimac\_0008636 & H3K4 methyltransferase     & LOC\_00003880 & 8 \\
\textit{KIF4A}        & obimac\_0008680 & Chromatin condensation     & LOC\_00003783 & 8 \\
\textit{CREBBP/CBP}   & obimac\_0008915 & Histone acetyltransferase  & LOC\_00003925 & 8 \\
\textit{MEN1}         & obimac\_0008998 & MLL-complex scaffold       & LOC\_00003829 & 8 \\
\textit{ARID4B}       & obimac\_0009001 & ARID chromatin reader      & LOC\_00003831 & 8 \\
\textit{BAHD1}        & obimac\_0009029 & H3K27 reader               & LOC\_00003878 & 8 \\
\textit{PAN3}         & obimac\_0009030 & mRNA deadenylase           & LOC\_00003838 & 8 \\
\textit{ELM2-domain}  & obimac\_0009036 & Chromatin remodeler        & LOC\_00003928 & 8 \\
\textit{THOC6}        & obimac\_0009064 & THO/TREX export            & LOC\_00003826 & 8 \\
\textit{Domino/SRCAP} & obimac\_0009065 & SWR1 H2A.Z remodeler       & LOC\_00003883 & 8 \\
\textit{BRE1/RNF40}   & obimac\_0009144 & H2B ubiquitin ligase       & LOC\_00003828 & 8 \\
\textit{TDRD3}        & obimac\_0008806 & Tudor methyl reader        & LOC\_00003695 & 6 \\
\textit{CRAMP1}       & obimac\_0008683 & Chromatin remodeler        & LOC\_00003786 & 5 \\
\textit{FBL}          & obimac\_0009198 & rRNA methyltransferase     & LOC\_00015020 & 5 \\
\addlinespace[4pt]
\multicolumn{5}{l}{\textit{Transcription factor (6)}} \\
\addlinespace[2pt]
\textit{FOS/AP-1}     & obimac\_0008668 & bZIP TF                    & LOC\_00003775 & 8 \\
\textit{TBX3}         & obimac\_0008756 & T-box TF                   & LOC\_00003751 & 8 \\
\textit{FoxA2}        & obimac\_0009054 & Forkhead pioneer TF        & LOC\_00003746 & 8 \\
\textit{OTX2}         & obimac\_0009208 & Homeobox TF                & LOC\_00003911 & 8 \\
\textit{ZNF226}       & obimac\_0008819 & Zinc finger TF             & LOC\_00003931 & 6 \\
\textit{MEIS2}        & obimac\_0009102 & Homeobox TF                & LOC\_00003687 & 5 \\
\bottomrule
\end{tabular}
\end{table}

% Cross-species chromatin comparison table moved to Supplementary Data

\begin{figure}[H]
    \centering
    \includegraphics[width=\linewidth]{figures/soual_dosage_compensation_3panel.png}
    \caption{Dosage compensation of the Z chromosome in \textit{Sthenoteuthis oualaniensis}.
    Left panels: gene expression levels across all chromosomes in females (ZO, top) and males (ZZ, bottom).
    Z:A ratios are 0.95 in females and 0.99 in males.
    Right panel: M:F expression ratios for autosomal genes (n = 26,351) and Z-linked genes (n = 419).
    The Z chromosome shows significant male bias (median M:F = 1.13 vs.\ 1.04 for autosomes; p = 0.010).}
    \label{fig:soual_dosage}
\end{figure}

\begin{figure}[H]
    \centering
    \includegraphics[width=\linewidth]{figures/dpe_dosage_compensation_3panel.pdf}
    \caption{Dosage compensation of the Z chromosome in \textit{Doryteuthis pealeii}.
    Left panels: gene expression levels across all chromosomes in females (ZO, top) and males (ZZ, bottom).
    Z:A ratios are 0.81 in females and 0.84 in males.
    Right panel: M:F expression ratios for autosomal genes (n = 6,899) and Z-linked genes (n = 135).
    The Z chromosome shows a median M:F ratio of 1.10 vs.\ 1.12 for autosomes (p = 0.29),
    consistent with effective dosage compensation.
    10F + 9M somatic samples (PRJNA255916 + PRJNA906054); DESeq2 normalization.}
    \label{fig:dpe_dosage}
\end{figure}

\begin{figure}[H]
    \centering
    \includegraphics[width=\linewidth]{figures/esco_dosage_compensation_3panel.pdf}
    \caption{Dosage compensation of the Z chromosome in \textit{Euprymna scolopes}.
    Left panels: gene expression levels across all chromosomes in females (ZO, top) and males (ZZ, bottom).
    Z:A ratios are 0.84 in females and 0.92 in males.
    Right panel: M:F expression ratios for autosomal genes (n = 15,690) and Z-linked genes (n = 237).
    The Z chromosome shows a median M:F ratio of 1.08 vs.\ 1.00 for autosomes (p = 0.003),
    consistent with near-complete dosage compensation.
    3F + 3M somatic samples (PRJNA1156170); DESeq2 normalization.}
    \label{fig:esco_dosage}
\end{figure}

\begin{figure}[H]
    \centering
    \includegraphics[width=\linewidth]{figures/eber_dosage_compensation_3panel_PRJNA1336213.png}
    \caption{Dosage compensation of the Z chromosome in \textit{Euprymna berryi}.
    Left panels: gene expression levels across all chromosomes in females (ZO, top) and males (ZZ, bottom).
    Z:A ratios are 0.98 in females and 0.91 in males.
    Right panel: M:F expression ratios for autosomal genes (n = 17,006) and Z-linked genes (n = 250).
    The Z chromosome shows a median M:F ratio of 0.97 vs.\ 0.97 for autosomes (p = 0.48),
    consistent with complete dosage compensation.
    2F + 2M somatic samples; DESeq2 normalization.}
    \label{fig:eber_dosage}
\end{figure}

\begin{figure}[H]
    \centering
    \includegraphics[width=\linewidth]{figures/soff_dosage_compensation_3panel.png}
    \caption{Dosage compensation of the Z chromosome in \textit{Sepia officinalis}.
    Left panels: gene expression levels across all chromosomes in females (ZO, top) and males (ZZ, bottom).
    Z:A ratios are 1.42 in females and 1.39 in males.
    Right panel: M:F expression ratios for autosomal genes (n = 17,647) and Z-linked genes (n = 213).
    The Z chromosome shows a slight but significant male bias (median M:F = 1.03 vs.\ 0.99 for autosomes; p = 0.021),
    indicating near-complete dosage compensation.
    4F + 4M somatic samples; DESeq2 normalization.}
    \label{fig:soff_dosage}
\end{figure}

\begin{figure}[H]
    \centering
    \includegraphics[width=\linewidth]{figures/obimac_dosage_compensation_3panel.png}
    \caption{Dosage compensation of the Z chromosome in \textit{Octopus bimaculoides}.
    Left panels: gene expression levels across all chromosomes in females (ZO, top) and males (ZZ, bottom).
    Z:A ratios are 1.09 in females and 0.88 in males.
    Right panel: M:F expression ratios for autosomal genes (n = 15,620) and Z-linked genes (n = 294).
    The Z chromosome shows a median M:F ratio of 1.08 vs.\ 1.03 for autosomes (p = 0.20),
    consistent with near-complete dosage compensation.}
    \label{fig:obimac_dosage}
\end{figure}


\begin{figure}[H]
    \centering
    \includegraphics[width=0.9\linewidth]{figures/za_isoform_regression_scatter.png}
    \caption{Isoform proportions of chromatin modifiers predict Z:A dosage compensation ratio
    in \textit{O. bimaculoides}.
    Scatter plots of female isoform proportion versus expression-matched Z:A ratio
    for four chromatin modifier genes (\textit{PBRM1}, \textit{SMARCC2}, \textit{KDM3A}, \textit{Domino}).
    The primary model (no covariates) explains 56.8\% of Z:A variance
    ($R^2 = 0.568$, $F = 8.20$, $p = 2.27 \times 10^{-4}$),
    with \textit{KDM3A} as the sole significant predictor
    (coefficient = +17.96, $p = 0.008$, Bonferroni $p = 0.032$).
    30 somatic samples (15F:15M).}
    \label{fig:za_regression}
\end{figure}

\begin{figure}[H]
    \centering
    \includegraphics[width=0.8\linewidth]{figures/kdm3a_za_scatter_by_sex.pdf}
    \caption{\textit{KDM3A} female isoform proportion versus Z:A ratio stratified by sex
    in \textit{O. bimaculoides}.
    The positive association between \textit{KDM3A} demethylase isoform switching
    and Z chromosome compensation is evident in both sexes,
    with females (ZO) clustering at higher \textit{KDM3A} female isoform proportions
    and higher Z:A ratios.}
    \label{fig:kdm3a_za}
\end{figure}

\begin{figure}[H]
    \centering
    \includegraphics[width=0.8\linewidth]{figures/zmast_vs_zfest_scatter.png}
    \caption{\textit{Zmast} and \textit{Zfest} expression in \textit{O. bimaculoides}.
    Scatter plot of \textit{Zmast} (male-biased) versus \textit{Zfest} (female-biased) TPM
    across somatic samples.
    \textit{Zmast} shows a male-to-female ratio of 5,063:1 for the core exon,
    while \textit{Zfest} is $\sim$18-fold female-biased.
    The mutually exclusive expression pattern is consistent with
    opposing regulatory roles in Z chromosome dosage compensation.}
    \label{fig:zmast_zfest}
\end{figure}

\begin{longtable}{lllll}
\caption{RNA-seq samples used in the dosage compensation analysis.
All samples are adult somatic tissues unless otherwise noted.
Accession prefixes: SRR (NCBI SRA), CRR (CNGB Sequence Archive).}
\label{tab:dc_samples} \\
\toprule
\textbf{Species} & \textbf{Study} & \textbf{Accession} & \textbf{Sex} & \textbf{Tissue} \\
\midrule
\endfirsthead
\toprule
\textbf{Species} & \textbf{Study} & \textbf{Accession} & \textbf{Sex} & \textbf{Tissue} \\
\midrule
\endhead
\midrule
\multicolumn{5}{r}{\textit{Continued on next page}} \\
\endfoot
\bottomrule
\endlastfoot
% Illex illecebrosus (1F/1M)
\multicolumn{5}{l}{\textit{Illex illecebrosus} (1F, 1M)} \\
\addlinespace[2pt]
 & This study & F24\_full & F & Mantle \\
 & This study & M71\_full & M & Mantle \\
\addlinespace[4pt]
% Sthenoteuthis oualaniensis (16F/16M)
\multicolumn{5}{l}{\textit{Sthenoteuthis oualaniensis} (16F, 16M)} \\
\addlinespace[2pt]
 & CRA004867 & CRR317115 & M & Photophore \\
 & CRA004867 & CRR317116 & M & Tentacle \\
 & CRA004867 & CRR317117 & M & Suckers \\
 & CRA004867 & CRR317118 & M & Heart \\
 & CRA004867 & CRR317119 & M & Kidney \\
 & CRA004867 & CRR317120 & M & Photophore \\
 & CRA004867 & CRR317121 & M & Tentacle \\
 & CRA004867 & CRR317122 & M & Suckers \\
 & CRA004867 & CRR317123 & M & Heart \\
 & CRA004867 & CRR317124 & M & Kidney \\
 & CRA004867 & CRR317157 & M & Hepatopancreas \\
 & CRA004867 & CRR317159 & M & Brain \\
 & CRA004867 & CRR317160 & M & Eye \\
 & CRA004867 & CRR317165 & M & Hepatopancreas \\
 & CRA004867 & CRR317167 & M & Brain \\
 & CRA004867 & CRR317168 & M & Eye \\
 & CRA004867 & CRR317137 & F & Kidney \\
 & CRA004867 & CRR317138 & F & Hepatopancreas \\
 & CRA004867 & CRR317140 & F & Brain \\
 & CRA004867 & CRR317141 & F & Kidney \\
 & CRA004867 & CRR317142 & F & Hepatopancreas \\
 & CRA004867 & CRR317144 & F & Brain \\
 & CRA004867 & CRR317145 & F & Eye \\
 & CRA004867 & CRR317151 & F & Eye \\
 & CRA004867 & CRR317161 & F & Photophore \\
 & CRA004867 & CRR317162 & F & Tentacle \\
 & CRA004867 & CRR317163 & F & Suckers \\
 & CRA004867 & CRR317164 & F & Heart \\
 & CRA004867 & CRR317169 & F & Photophore \\
 & CRA004867 & CRR317170 & F & Tentacle \\
 & CRA004867 & CRR317171 & F & Suckers \\
 & CRA004867 & CRR317172 & F & Heart \\
\addlinespace[4pt]
% Doryteuthis pealeii (10F/9M)
\multicolumn{5}{l}{\textit{Doryteuthis pealeii} (10F, 9M)} \\
\addlinespace[2pt]
 & PRJNA255916 & SRR1522987 & M & Giant fiber lobe \\
 & PRJNA255916 & SRR1522988 & M & Optic lobe \\
 & PRJNA255916 & SRR1725163 & M & Gill \\
 & PRJNA255916 & SRR1725164 & M & Branchial heart \\
 & PRJNA255916 & SRR1725167 & F & Iridophore skin \\
 & PRJNA255916 & SRR1725169 & M & Marginal pen epithelium \\
 & PRJNA255916 & SRR1725171 & M & Ventral pen epithelium \\
 & PRJNA255916 & SRR1725172 & F & Chromatophore skin \\
 & PRJNA255916 & SRR1725213 & M & Buccal ganglia \\
 & PRJNA255916 & SRR1725235 & M & Stellate ganglion \\
 & PRJNA255916 & SRR1725236 & M & Vertical lobe \\
 & PRJNA906054 & SRR23062180 & F & Tentacle rim \\
 & PRJNA906054 & SRR23062181 & F & Skin \\
 & PRJNA906054 & SRR23062182 & F & Olfactory organ \\
 & PRJNA906054 & SRR23062183 & F & Brain \\
 & PRJNA906054 & SRR23062189 & F & Arm sucker rims \\
 & PRJNA906054 & SRR23062200 & F & Arm \\
 & PRJNA906054 & SRR23062211 & F & Arm \\
 & PRJNA906054 & SRR23062212 & F & Arm \\
\addlinespace[4pt]
% Euprymna scolopes (3F/3M)
\multicolumn{5}{l}{\textit{Euprymna scolopes} (3F, 3M)} \\
\addlinespace[2pt]
 & SRP412216 & SRR22588893 & F & Central brain \\
 & SRP412216 & SRR22588898 & F & Central brain \\
 & SRP412216 & SRR22588901 & F & Left optic lobe \\
 & SRP412216 & SRR22588910 & M & Central brain \\
 & SRP412216 & SRR22588918 & M & Central brain \\
 & SRP412216 & SRR22588927 & M & Left optic lobe \\
\addlinespace[4pt]
% Euprymna berryi (2F/2M)
\multicolumn{5}{l}{\textit{Euprymna berryi} (2F, 2M)} \\
\addlinespace[2pt]
 & PRJNA1336213 & SRR35682649 & F & Skin \\
 & PRJNA1336213 & SRR35682652 & F & Brain \\
 & PRJNA1336213 & SRR35682659 & M & Brain \\
 & PRJNA1336213 & SRR35682663 & M & Skin \\
\addlinespace[4pt]
% Sepia officinalis (4F/4M)
\multicolumn{5}{l}{\textit{Sepia officinalis} (4F, 4M)} \\
\addlinespace[2pt]
 & PRJNA242869 & SRR5204441 & F & Stellar ganglion \\
 & PRJNA242869 & SRR5204442 & M & Stellar ganglion \\
 & PRJNA242869 & SRR5204443 & F & Optic lobe \\
 & PRJNA242869 & SRR5204444 & M & Optic lobe \\
 & PRJNA242869 & SRR5204445 & F & Subesophageal mass \\
 & PRJNA242869 & SRR5204446 & M & Subesophageal mass \\
 & PRJNA242869 & SRR5204447 & F & Supraesophageal mass \\
 & PRJNA242869 & SRR5204448 & M & Supraesophageal mass \\
\addlinespace[4pt]
% Octopus bimaculoides (15F/15M)
\multicolumn{5}{l}{\textit{Octopus bimaculoides} (15F, 15M)} \\
\addlinespace[2pt]
 & PRJNA285380 & SRR2047107 & M & Posterior salivary gland \\
 & PRJNA285380 & SRR2047111 & M & Suckers \\
 & PRJNA285380 & SRR2047116 & M & Retina \\
 & PRJNA285380 & SRR2047118 & M & Optic lobe \\
 & PRJNA285380 & SRR2047120 & F & Axial nerve cord \\
 & PRJNA285380 & SRR2048495 & F & Subesophageal brain \\
 & PRJNA285380 & SRR2048496 & F & Subesophageal brain \\
 & PRJNA285380 & SRR2048497 & F & Subesophageal brain \\
 & PRJNA285380 & SRR2048498 & F & Subesophageal brain \\
 & PRJNA285380 & SRR2048521 & F & Supraesophageal brain \\
 & PRJNA285380 & SRR2048522 & F & Supraesophageal brain \\
 & PRJNA285380 & SRR2048523 & F & Supraesophageal brain \\
 & PRJNA285380 & SRR2048524 & F & Supraesophageal brain \\
 & PRJNA285380 & SRR2048525 & F & Supraesophageal brain \\
 & PRJNA1199547 & SRR32191430 & M & Brain \\
 & PRJNA1199547 & SRR32191431 & M & Brain \\
 & PRJNA1199547 & SRR32191434 & M & Arm \\
 & PRJNA1199547 & SRR32191435 & M & Arm \\
 & PRJNA1336213 & SRR35682634 & F & Brain \\
 & PRJNA1336213 & SRR35682635 & F & Brain \\
 & PRJNA1336213 & SRR35682636 & M & Skin \\
 & PRJNA1336213 & SRR35682646 & M & Sucker rim \\
 & PRJNA1336213 & SRR35682657 & M & Olfactory organ \\
 & PRJNA1336213 & SRR35682668 & M & Lip \\
 & PRJNA1336213 & SRR35682679 & M & Hectocotylus \\
 & PRJNA1336213 & SRR35682686 & F & Skin \\
 & PRJNA1336213 & SRR35682687 & F & Oviduct \\
 & PRJNA1336213 & SRR35682688 & F & Oviduct \\
 & PRJNA1336213 & SRR35682690 & M & Eye \\
 & PRJNA1336213 & SRR35682691 & M & Brain \\
\end{longtable}

\end{document}