What is the justification of Chargaff's second parity rule?

Question

First parity rule

The first rule holds that a double-stranded DNA molecule globally has percentage base pair equality: %A = %T and %G = %C. The rigorous validation of the rule constitutes the basis of Watson-Crick pairs in the DNA double helix model.

Second parity rule

The second rule holds that both %A ~ %T and %G ~ %C are valid for each of the two DNA strands. This describes only a global feature of the base composition in a single DNA strand.

[source]

It makes sense that in the context of dsDNA, that A = T and C = G, but I don't see an obvious reason why in a single strand of DNA, A ~ T and C ~ G.

Answer 1

This is an excellent question, which was the cause of some confusion in the field over the years:

"The validity of the second Chargaff rule was unexpected. Obviously it should be regarded as a global rule, i.e. applicable to large sections of chromosomes. Nonetheless, not being derived from a compelling principle, such as the one underlying the first rule, it remains a mystery. This is even more so, when one studies extended versions of Chargaff’s second rule." (paper)

There are a few things to note here:

1. Chargaff's rules were purely observational; they pre-existed the double helix (see your same Wikipedia ref), and thus pre-existed the double helix model (which was in fact famously a rendering of the first rule in a 3-D model, as Watson wrote about in his book "The Double Helix"). There is thus no logical or mechanistic justification for the rules; they were simply what Chargaff et al. observed for dsDNA and ssDNA. Perhaps I have misinterpreted, but I felt the question implied a belief that there must be some logical requirement for this.
2. Note that the first rule uses "=" and the second rule uses "~" ("approximately"). That is, if we believe that all bases in a genome are paired, then it follows as a mechanistic/logical requirement that Watson strand and Crick strand must match up compositionally. This is obviously not the case with ssDNA, as I think you indicate in the Q. However, it is nonetheless true that genomes have overall compositional biases; they tend to have consistent (G+C)/(A+T) ratios across the molecule (on average). More recently, it has become clear that this is true not only for single nucleotides, but also for $k$-mers, or short DNA sequences of length $k$ (paper). Indeed, these patterns have been used to match up different pieces from the same genomes (paper). These biases exist on both strands, therefore each strand of a genome is subject to these biases. (The reasons why these biases exist are often argued over, but may have something to do with mutational biases). Therefore, both ssDNA strands of a dsDNA genome have the same compositional biases, and they must complement each other's similar biases, so even within a strand the statistical expectation is Chargaff's second rule. Chargaff, again, probably didn't know many of these things when he came up with the rules. This explanation for Chargaff's second rule has been explicitly discussed in the literature, in the first paper cited.

Notably, the second rule only holds reliably for long sequences, because in short sequences sampling error may lead to a failure to match the biases of forward and reverse strands. This is more or less what everyone means when they call it a "global" rule.

Answer 2

In addition to the answer by Maximilian Press, the papers

• [1] Intrastrand Parity Rules of DNA Base Composition and Usage Biases of Synonymous Codons
• [2] Properties of a General Model of DNA Evolution Under No-Strand-Bias Conditions

show that under a no-strand-bias assumption, the base frequencies of a single strand at equilibrium, denoted as $A^*,T^*,G^*,C^*$, are such that $A^*=T^*$ and $G^*=C^*$.

I would try to explain it briefly here, but please refer to the original papers for a complete and better explanation, as some of my explanation is simplified.

---

First, we denote the base frequencies (in time $t$) of a single strand $A(t),T(t),G(t),C(t)$.
(each frequency is in $[0,1]$, and $A(t)+T(t)+G(t)+C(t)=1$)

Next, we denote the substitution rates in a single strand $r_1, ..., r_{12}$, as follows:

(diagram copied from [1])

E.g.: $$\frac{dA}{dt} = -\left(r_{1}+r_{4}+r_{6}\right)\cdot A\left(t\right)+r_{2}\cdot T\left(t\right)+r_{3}\cdot G\left(t\right)+r_{5}\cdot C\left(t\right)$$

Similarly, we denote the substitution rates in the other strand (in the double helix) $r'_1, ..., r'_{12}$.

If we assume no-strand-biases$^{(*)}$ (i.e., the substitution rates are independent of the strand: $r_1=r'_1, ..., r_{12}=r'_{12}$), then using Chargaff's first parity rule we can deduce: $$r_1=r_2,r_3=r_9,r_4=r_{10},r_5=r_7,r_6=r_8,r_{11}=r_{12}$$ (because $r_7=r'_5$, etc.)

---

$(*)$ If I understand correctly, there is another implicit assumption here - each substitution rate in time $t$ is independent of the base frequencies. E.g., regardless of whether the frequency of $A$ bases in one strand is different from the frequency of $A$ bases in the other strand, we assume that $r_1=r'_1$.

---

So we can (under a no-strand-bias assumption) denote the substitution rates in a single strand $a, ..., f$, as follows:

(diagram copied from [1])

I.e.: $$ \left(\begin{matrix}\frac{dA}{dt}\\ \frac{dT}{dt}\\ \frac{dG}{dt}\\ \frac{dC}{dt} \end{matrix}\right)=\left(\begin{matrix}\left(-a-e-c\right)\cdot A\left(t\right)+a\cdot T\left(t\right)+b\cdot G\left(t\right)+d\cdot C\left(t\right)\\ a\cdot A\left(t\right)+\left(-a-e-c\right)\cdot T\left(t\right)+d\cdot G\left(t\right)+b\cdot C\left(t\right)\\ c\cdot A\left(t\right)+e\cdot T\left(t\right)+\left(-b-d-f\right)\cdot G\left(t\right)+f\cdot C\left(t\right)\\ e\cdot A\left(t\right)+c\cdot T\left(t\right)+f\cdot G\left(t\right)+\left(-b-d-f\right)\cdot C\left(t\right) \end{matrix}\right) $$ Or equivalently: $$ \left(\begin{matrix}-a-e-c & a & b & d\\ a & -a-e-c & d & b\\ c & e & -b-d-f & f\\ e & c & f & -b-d-f \end{matrix}\right)\left(\begin{matrix}A\left(t\right)\\ T\left(t\right)\\ G\left(t\right)\\ C\left(t\right) \end{matrix}\right)=\left(\begin{matrix}\frac{dA}{dt}\\ \frac{dT}{dt}\\ \frac{dG}{dt}\\ \frac{dC}{dt} \end{matrix}\right) $$

At equilibrium: $$ \left(\begin{matrix}\frac{dA}{dt}\\ \frac{dT}{dt}\\ \frac{dG}{dt}\\ \frac{dC}{dt} \end{matrix}\right)=\left(\begin{matrix}0\\ 0\\ 0\\ 0 \end{matrix}\right) $$

So we have to solve: $$ \left(\begin{matrix}-a-e-c & a & b & d\\ a & -a-e-c & d & b\\ c & e & -b-d-f & f\\ e & c & f & -b-d-f \end{matrix}\right)\left(\begin{matrix}A\left(t\right)\\ T\left(t\right)\\ G\left(t\right)\\ C\left(t\right) \end{matrix}\right)=\left(\begin{matrix}0\\ 0\\ 0\\ 0 \end{matrix}\right) $$

It turns out that every vector proportional to the following solves the equation: $$ \left(\begin{matrix}d+b\\ d+b\\ e+c\\ e+c \end{matrix}\right) $$

The relevant solution to us is a solution such that each frequency is in $[0,1]$ and $A(t)+T(t)+G(t)+C(t)=1$: $$ \left(\begin{matrix}\frac{d+b}{2\left(d+b+e+c\right)}\\ \frac{d+b}{2\left(d+b+e+c\right)}\\ \frac{e+c}{2\left(d+b+e+c\right)}\\ \frac{e+c}{2\left(d+b+e+c\right)} \end{matrix}\right) $$ I.e., $A^*=T^*$ and $G^*=C^*$.