In search of the genetic variants of human sex ratio at birth: Was Fisher wrong about sex ratio evolution?

The human sex ratio (fraction of males) at birth is close to 0.5 at the population level, an observation commonly explained by Fisher's principle. However, past human studies yielded conflicting results regarding the existence of sex ratio-influencing mutations-a prerequisite to Fisher’s principle, raising the question of whether the nearly even population sex ratio is instead dictated by the random X/Y chromosome segregation in male meiosis. Here we show that, because a person’s offspring sex ratio (OSR) has an enormous measurement error, a gigantic sample is required to detect OSR-influencing genetic variants. Conducting a UK Biobank-based genome-wide association study that is more powerful than previous studies, we detect an OSR-associated genetic variant, which awaits verification in independent samples. Given the abysmal precision in measuring OSR, it is unsurprising that the estimated heritability of OSR is effectively zero. We further show that OSR’s estimated heritability would remain virtually zero even if OSR is as genetically variable as the highly heritable human standing height. These analyses, along with simulations of human sex ratio evolution under selection, demonstrate the compatibility of the observed genetic architecture of human OSR with Fisher’s principle and suggest the plausibility of presence of multiple human OSR-influencing genetic variants.

Methods
GWAS:

When conducting the GWAS in the UKB, we did not simply use the sibling sex ratio as the trait, because of the difficulty in accounting for different estimation errors of the sibling sex ratio for different families as a result of the variation in family size.  For example, individual A has one brother and zero sister, while individual B has four brothers and one sister.  Although A has a higher sibling sex ratio than B, B’s siblings obviously provide stronger evidence for a male-biased sibling sex ratio than A’s siblings.  To properly weigh the data by the family size, we considered the birth of each sibling as an independent event.  In the above example, we would associate A’s genotype with one male birth and associate B’s genotype with four male births and one female birth.  In GWAS, a male birth is coded as 1 and a female birth is coded as 0.  The UKB participants have a total of 873,715 full siblings, leading to an unprecedented statistical power.  In our GWAS in the UKB, we included genetic sex, year of birth, and the first ten genetic principle components as covariates.

Gene-based test:

We performed two gene-based association analyses.  First, we analyzed the UKB-based GWAS summary statistics through the R package sumFREGAT for autosomal protein-coding genes (N = 17,389).  All SNPs within the transcribed region of a gene derived from the European samples in the 1000 Genome Project were used in the test.  We implemented the optimal unified test (SKAT-O), principal component analysis-based test (PCA), and aggregated Cauchy association test (ACAT-V) in sumFREGAT.  For all three methods, weights were uniformly assigned for all alleles [beta.par = c(1, 1) in sumFREGAT] with other settings left at default values.  Variant correlation matrix files (one file per gene) were needed for the gene-based analysis, and we used the pre-calculated matrices from 1KG European samples provided by the R-package development team (http://mga.bionet.nsc.ru/sumFREGAT).  The input data were pre-processed using the R package function prep.score.files() with the reference file provided by the R-package development team (http://mga.bionet.nsc.ru/sumFREGAT).  The P values in the three tests were then combined by the omnibus aggregated Cauchy association test (ACAT-O) in sumFREGAT.

Second, we performed a gene-based burden test using rare missense variants (MAF < 1%) in the UKB whole exome sequencing data.  The burden test assumes that rare variants are functionally disruptive and therefore have the same direction of effect.  To properly weigh OSR of UKB participants by their heterogenous measurement errors, we generated a plink bed file that contained burden scores of all genes for all UKB individuals using the “--write-mask” option in REGENIE.  The annotation file that specifies the functional class of each SNP and the corresponding gene required in this step was provided in the UKB Research Analysis Platform (see https://dnanexus.gitbook.io/uk-biobank-rap/science-corner/using-regenie-to-generate-variant-masks), which included protein coding genes in autosomes, X, and Y chromosomes (N = 18,845).  We chose to include all loss-of-function and missense SNPs to calculate the burden score.  In the default setting, the burden score is calculated as the maximum number of alternative alleles across sites of a gene, being 0, 1, or 2 (see REGENIE online documentation for details, https://rgcgithub.github.io/regenie/options/).  We then used this gene-level bed file to perform association analysis on the sibling sex following the same procedure describe in the “GWAS” section.

Simulating the genetic architecture of sex ratio following that of standing height

To simulate the genetic architecture of sex ratio following that of human standing height, we obtained the hypothetical sex ratio of a participant of European ancestry in the UKB through the following four steps.  First, we computed the hypothetical sex ratio of a participant by dividing the participant’s standing height by twice the mean standing height of all UKB participants of European ancestry.  Second, we performed a multiple regression on hypothetical sex ratio; the independent variables included genetic sex, age, age squared, and the first ten genetic principal components but not SNPs.  Third, we obtained the regression residual of each participant, which is the difference between the hypothetical sex ratio computed in the first step and that predicted by the multiple regression model in the second step.  Fourth, the covariate-corrected hypothetic sex ratio was set to be the regression residual in the preceding step plus 0.5.  GWAS was subsequently performed on the covariate-corrected hypothetic sex ratio.  SNP-based heritability of the covariate-corrected hypothetical sex ratio was computed.  Based on the covariate-corrected hypothetical sex ratio, we generated the sexes of each participant’s offspring with 20 replicates.  To ensure comparability with the original GWAS data, we assumed that each participant had the same number of offspring as the number of siblings in the UKB.  We then conducted a GWAS using the simulated sexes of all offspring and estimated the SNP-based heritability of the estimated hypothetical sex ratio.

Simulations of human sex ratio evolution

We used SLiM 3 to simulate sex ratio evolution in humans.  A non-Wright-Fisher model with separate sexes and non-overlapping generations was enabled in the simulation, along with the human demographic history described by the default example code in SLiM 3 (see SLiM manual, https://messerlab.org/slim/, p. 136-142).  The diploid genome has a pair of 1000-nt chromosomes, and the recombination rate is 1×10-3 per site per generation such that one recombination per chromosome per generation is expected.  In every generation, males and females will mate randomly, and each mating will result in one offspring.  The random mating continues until the number of offspring matches the expected population size in the next generation.  To achieve the mutation-drift-selection equilibrium, the population was pre-evolved for 73,105 generations (10 times the effective population size) in every simulation.

The mutation rate varied from 1×10-6 to 1×10-2 per genome per generation.  The mean mutation size () varied from 0.00125 to 0.16.  Given , the actual size of a mutation is sampled from an exponential distribution with a mean of .  The genetic effect of the mutation is set to be paternal.  Thirty simulation replications were performed for each combination of mutation rate and mean size.

Under the directional selection scenario, we assumed that the optimal OSR changed from the default value of 0.5 to around 0.52 at 800,000 years before present.  To set the optimal OSR at around 0.52, we introduced unbalanced parental investments by reduce the future mating probability of individuals who have had daughters: future mating probability = 1 – 0.1 × number of daughters.  The optimal OSR is 0.524, which was estimated by averaging sex ratios at the last 10 time points of 10-generation intervals in all simulations where mutation rate is 0.01 and mean mutation size is 0.00125, 0.0025, or 0.005.

The heritability of sex ratio (with measurement error) was calculated by dividing the variance of genetically expected sex by the variance of observed sex.  To obtain the number of detectable variants, we used the UKB statistical power map generated earlier (Fig. 1c).  A SNP was considered detectable if its detectability exceeded 0.9.  Key statistics such as the heritability of sex ratio (with measurement error), number of detectable variants, and number of variants in each simulation replicate were calculated by averaging sex ratios at the final 10 time points where consecutive time points were separated by 10 generations.  These statistics from the 30 replicates were used to plot the mean, maximum, and minimum in Fig. 4.