## 2 Materials and Methods
### Mendelian Randomization  
Pair-wise associations found to be significant using cross-trait LDSC, and suggested to be part of the same clusters based on k-means clustering, were tested with two-sample MR using the TwoSampleMR R package [44; 45]. MR uses genetic variants as instrumental variables (IV) to estimate the causal effects of an exposure on an outcome [46]. To infer the causal influence of the exposure, the ratio between the SNP effect on the outcome over the SNP effect on the exposure is computed [45]. The method relies on the fact that, based on [PERSON]'s laws of inheritance and the fixed nature of germline genotypes, the alleles at a SNP are expected to be random with respect to potential confounders [45]. The results obtained from multiple SNPs associated with the exposure allow us to obtain an overall estimate of the causal effect of a potential exposure on an outcome. Two-sample MR can be performed even if the SNP-exposure effects and the SNP-outcome effects are obtained from separate studies, thus allowing to leverage pre-existing large GWAS. In order to assess the direction of the association between the selected pairs of traits, we repeated the analyses twice for each pair, considering each trait to be the exposure or the outcome, alternatively. IVs were selected based on significant association with the exposure at a genome-wide threshold (\(p<5\times 10^{-8}\)). Since it is important to ensure that IVs are independent, significant genetic variants were clumped using European data from the IEU GWAS database as recommended [44]. Default parameters in the TwoSampleMR package were used for clumping (\(r^{2}>0.001\) in the range of 10,000 Mb) and the SNP with the lowest \(p\)-value was retained. After clumping, exposure and outcome genetic data were harmonized to obtain effects and standard errors for each instrument SNP available for the exposure and outcome traits, using the Two-SampleMR package. During this step, palindromic SNPs with intermediate allele frequencies are removed by the package.  
One important assumption of MR is that genetic variants used as IVs should exert an effect on the outcome only through their effect on the exposure. A violation of this assumption is called horizontal pleiotropy (i.e., a condition in which a genetic variant exerts an effect on the outcome through different pathways) and can cause bias in the MR analysis. We checked the intercept term in MR Egger regression [47], as well as the global and the distortion tests implemented in the MR-PRESSO package [48], in order to assess whether directional horizontal pleiotropy was driving the results of the MR analysis. MR analyses were conducted with four widely used different methods implemented in the TwoSampleMR package: MR-Egger regression, weighted median estimator, inverse-varianceweighted (IVW) and simple mode. In addition, the MR analyses were conducted using the raw and outlier corrected estimates of the mendelian randomization pleiotropy residual sum and outlier (MR-PRESSO) test, implemented in the MR-PRESSO R package [48]. Based on the observation of a significant bidirectional association between MDD and insomnia, for these traits the analyses were also repeated using the different versions of the MDD and insomnia datasets which are publicly available. Specifically, we used as exposure SNPs extracted from the versions of the datasets limited to 10,000 variants selected by the authors of the two GWAS for which data for participants from 23 andMe are also included, as well as the list of significant and independent genetic variants reported by the two GWAS (n = 102 for MDD [7] and n = 554 for insomnia [33]). Aim of these analyses was to verify whether a significant bidirectional association between the two traits could still be observed when using IVs extracted from datasets including a limited number of variants genotyped in a larger number of participants.