The subjects used in this study were recruited during 2005–2011 primarily from the Mid-South states in the United States as part of the Mid-South Tobacco Case-Control (MSTCC) study. For a detailed description of these samples, please refer to our previously reported studies (Yang et al., 2015; Yang and Li, 2016). All materials related to subject recruitment were approved by the local ethics committees, and all subjects provided written informed consent.

Depressive symptoms were measured with a 20-item version of the Center for Epidemiological Studies Depression (CES-D) scale (Radloff, 1977), and the score of each subject from the baseline examination was used (Harlow et al., 1991; Kuchibhatla et al., 2012). Questionnaires assessing various other characteristics of interest were also administered, and individuals exhibiting other substance abuse (except smoking) or other psychiatric disorders (except depression) were excluded. As previously reported (Yang et al., 2015; Yang and Li, 2016; Jiang et al., 2019; Xu et al., 2020), genomic DNA was isolated from the peripheral blood of each patient, and genotyping was performed using the Illumina Infinium HumanExome BeadChip v 1.0 (Illumina Inc., San Diego, CA, United States). This exome chip includes more than 240,000 functional exonic variants that have been detected in several sequencing data sets, and it was designed to concentrate on rare variants.¹ A total of 4,817 unrelated participants with both genotypic and phenotypic data was enrolled in this study. More detailed descriptions of these participants are shown in Supplementary Table 1.

We excluded SNPs from further analysis if they: (1) were located on sex chromosomes or the mitochondrial genome; (2) had a call rate of <0.95; (3) deviated from the Hardy-Weinberg equilibrium at a p-value of <1 × 10^–6; or (4) had a minor allele frequency (MAF) of <0.01. After these quality control steps, 38,943 SNPs and 4,817 individuals remained for the analyses reported in this paper.

Single-marker association analysis was carried out with a linear regression model adjusted for age, sex, ethnicity, smoking status, and the first three MDS using PLINK (Purcell et al., 2007). Association analysis results were considered significant after Bonferroni correction for the 38,943 SNPs (i.e., p < 0.05/38,943 = 1.28 × 10^–6). Considering that only 11 SNPs were included on the Illumina Infinium HumanExome BeadChip for FZD6, we carried out SNP imputation with IMPUTE2 (Howie et al., 2012) after phasing in SHAPEIT2 (Delaneau et al., 2013) using the 1000 Genome Phase 3 data (Build 38) as a reference, which extended the number of SNPs in FZD6 from 11 to 132 after excluding those with low imputation quality (info < 0.3). Then, we conducted gene-based association analysis to examine the presence of association between FZD6 and the CES-D score using SNP-set Kernel Association Test (SKAT) software (Ionita-Laza et al., 2013). The EWAS dataset used in the current study has been submitted to public database and it is available from NCBI GEO (Accession: GSE148375).²

SWISS-MODEL³ was used to generate a three-dimensional (3D) structure for wild-type (WT) FZD6 protein (Arnold et al., 2006). A template of human Smoothened in complex with cholesterol (SMTL ID 5I7d.1), which was searched through a BlastP algorithm against the PDB database, showed 27.18% identity and was used for homology modeling of FZD6. The mutant models also were constructed and analyzed in relation to the WT by PyMOL.⁴

Changes in protein stability resulting from a single point mutation were predicted by I-Mutant 3.0⁵ based on the Support Vector Machine (SVM) algorithm. It predicted the ΔΔG value as a regression estimator and sign of the stability change, and classified mutations into three categories: neutral (−0.5 − 0.5), decrease (<−0.5), and increase (>0.5) (Schymkowitz et al., 2005). We also used Mutation Cutoff Scanning Matrix (mCSM; Pires et al., 2014b), which applies the concept of graph-based structural signatures to predict the effect of a single point mutation on protein stability. DUET, an integrated computational web server, applies mCSM to predict missense mutation effects on protein stability (Pires et al., 2014a). The ΔΔG value < 0 means a decrease in the stability; otherwise, an increase in the stability. NetSurfP 2.0⁶ is a sequence-based tool to predict protein structures, including secondary structure, structural disorder, and backbone dihedral angles for each residue of the input sequences (Klausen et al., 2019).

Human embryonic kidney 293T (HEK293T) cells were purchased from the American Type Culture Collection (Manassas, VA, United States). Cells were cultured in Dulbecco’s Modified Eagle Medium (Hyclone, Logan, UT, United States) supplemented with 10% fetal bovine serum (Gibco, Grand Island, NY, United States) and 1% penicillin/streptomycin (Gibco) in a humidified atmosphere of 5% CO₂ at 37°C.

The FZD6 DNA fragment containing the C allele of rs61753730 was synthesized and cloned into the pGL3-promoter (Promega, Madison, WI, United States). The G allele was introduced into the rs61753730-C allele construct using a QuikChange Lightning Multi Site-Directed Mutagenesis Kit (Agilent, Santa Clara, CA, United States), and its sequence was verified by Sanger sequencing. When HEK293T cells reached 80–90% confluence, transfection was conducted with a pGL3-promoter construct and pRL-TK vector (Promega) by Lipofectamine 2000.

For the luciferase assay, cells were lysed 48 h after transfection and assayed with the Dual-Luciferase Report Assay System (Promega). For the mRNA decay assay, actinomycin D (5 μg/ml) was added 24 h after transfection, and total RNA was extracted at 0, 1, or 3 h. The quantitative real-time PCR (qRT-PCR) was then conducted to determine the relative expression of the luciferase reporter gene.

Total RNA was isolated from the cells after treatment and homogenized by an ultrasonic disruptor in TRIzol reagent (Life Technologies, Grand Island, NY, United States) according to the manufacturer’s instructions. The purity and quantity of total RNA were measured at the optical densities of 260 and 280 nm with NanoDrop 2000 (Thermo Scientific, Waltham, MA, United States). First-strand cDNAs were synthesized from 1 μg of total RNA using an iScript cDNA Synthesis Kit (Bio-Rad, Hercules, CA, United States). The qRT-PCR amplification was conducted in a volume of 10 μl containing 5 μl of 2 × Power SYBR Green PCR Master Mix (Applied Biosystems), mixed with equal amounts of sense and antisense primers (2.5 μl; final concentration 250 nM each), and 2.5 μl of diluted cDNA in a 384-well plate using 7900 HT Fast Real Time PCR system (Applied Biosystems) as described previously (Yang et al., 2016; Han et al., 2019; Chen et al., 2020; Fan et al., 2021). The PCR conditions were as follows: 50°C for 2 min, 95°C for 10 min, followed by 40 cycles of 95°C for 15 s and 60°C for 1 min. The primers for the reporter gene are forward: 5′-TGCACATATCGAGGTGGACATC-3′ and reverse: 5′-TGCCAACCGAACGGACAT-3′; both were synthesized at Sango Biotech (Shanghai) Co., Ltd. The expression of the gene of interest was normalized to the housekeeping gene glyceraldehyde 3-phosphate dehydrogenase (GAPDH) and then analyzed using a comparative C_t method (Winer et al., 1999).

All animals were purchased from GemPharmatech Co., Ltd., (Nanjing, China) and maintained in the Animal Core Facility at Zhejiang University. Mice were kept under the standard housing conditions on a 12-h light/dark cycle. Food and water were freely available. Mice were deeply anesthetized before sacrifice. All experimental procedures used in the study were approved by the Animal Use and Care Committee of the First Affiliated Hospital of Zhejiang University.

The Fzd6-knockin (Fzd6-KI) mice based on C57BL/6J were generated via the optimized CRISPR/Cas9 system at GemPharmatech Co., Ltd., with almost no off-target effects detected (Shen et al., 2014). Briefly, Cas9 mRNA, sgRNA, and donor compounds were co-microinjected into fertilized C57BL/6J mouse zygotes. The sgRNA directed to Cas9 endonuclease cleavage near the mutation site and created a double-strand break. Such a break was repaired and resulted in a Q152E point mutation, which was then inserted into exon 4 by homologous recombination and we obtained the Fzd6-KI mice.

DNA was isolated from the tail of each pup, and PCR was conducted for genotyping. The PCR products were electrophorized to screen if it is wild-type or mutant mouse. Then, Sanger sequencing was used to validate the gene sequence, especially at the mutant site. The sequence of sgRNA was 5′-CACCAAAATCCAATGTCTCT-3′, and the donor sequence was 5′-CCTGTAAC TTCTCATCCACACACAGAGCTTTCTGGGCCACAGAAGA AATCAGATGAAGTCCCACGCGACATTGGATTTTGGTGT CCAAAGCACCTTAGGACTTCCGGGGACCAAGGCTAT-3′. The primers for genotyping were forward 1: 5′- TCT GTGAATGCAGCAAAGTCATGG-3′, reverse 1: 5′- ATCCA ATGTCGCGTGGGACTTC-3′; forward 2: 5′-TCTGTGAATGC AGCAAAGTCATGG-3′, reverse 2: 5′-GTCTCTCTGGGTATC TGAATCGTC-3′. The primers for sequencing was 5′-GTC TCTCTGGGTATCTGAATCGTC-3′. After genotyping, through crossbreeding of the heterozygote mice, we obtained Fzd6-KI homozygous and WT mice, which were used for the following studies. For the identification of some representative animals with different genotypes, please see the attached Supplementary Material.

The adult (2–3-month old) male and female homozygous Fzd6-KI and WT mice were subjected to behavioral tests (n = 7–9 mice/genotype). All mice were habituated to the testing room for 1 h prior to each test. After completion of all relevant behavioral tests, all mice were sacrificed.

The chronic social defeat stress (CSDS) model of depression was carried out as described previously (Golden et al., 2011; Harris et al., 2018; Shen et al., 2019). Briefly, prior to each experiment, more than the number of 50 retired male breeder CD1 mice (known as aggressor mice) were screened for three consecutive days to determine the aggressive characteristics of each animal, which was done by placing the screener C57BL/6J mouse directly into the home cage of the CD1 mice for 3 min with the presence of an aggressor. The following two criteria were used to select the CD1 aggressor mice: (1) mouse must attack in at least two consecutive sessions during three 3-min screening sessions; and (2) the latency to initial aggression must be less than 60 s. Then, each experimental mouse (male Fzd6-KI mouse) was introduced into the home cage of a CD1 aggressor for 5 min. At the end of 5 min of social defeat, the experimental mice were moved to the other side separated from the aggressors with a perforated Plexiglas divider where animals were maintained for 24 h. After 24 h of sensory but not physical contact, the above procedure was repeated daily for 10 consecutive days. The non-defeated control animals were also pair-housed with an identical home cage setup. The two mice on each side were never allowed physical or sensory contact with CD1 mice or their home mates. Each experimental mouse was faced with a different aggressor every day. The social interaction test was conducted 24 h later after 10 days of social defeat.

Social avoidance behavior was evaluated by the social interaction (SI) test in an open-field apparatus (44 cm × 44 cm × 44 cm) consisting of two consecutive phases each for 2.5 min. During the first phase (target absent), experimental animals were placed individually in the apparatus containing an empty wire-mesh enclosure (10 cm × 6.5 cm × 25 cm) with the target CD1 aggressor absent, and their movements were monitored with Any-maze video tracking software. In the second phase (target present), a novel aggressive CD1 mouse was placed in the wire-mesh enclosure with the same metrics measured. Between the two phases, the experimental animals were returned to their home cages for an approximately 30 s interval. The time spent in the interaction zone (a 14 cm × 26 cm rectangular area projecting 8 cm around the wire-mesh enclosure) and in the corner zones (a 9 cm × 9 cm area projecting from both corner joints opposing the wire-mesh enclosure) were counted.

The experimental animals were placed individually in the center of the open field test (OFT) apparatus (44 cm × 44 cm × 44 cm) with a dim light for 10 min. A video camera was arranged directly above the apparatus to monitor the movement of each animal. Locomotor activity was measured as the total distance traveled. Anxiety-like behavior was defined by the distance traveled, time spent in the central zone (15 cm × 15 cm), and center entries during the 10 min period. The apparatus was cleaned with 75% ethanol after each experiment. Data were collected by Any-maze software (Any-maze, Stoelting, Wood Dale, IL, United States).

The forced swim test (FST) was conducted as previously described (Yang et al., 2018). Each experimental animal was placed individually in a cylinder (12 cm diameter; 25 cm height) filled with water (23–25°C) and allowed to swim for 6 min. The depth of water in the cylinder (∼16 cm) was set to prevent animals from touching the bottom with their tails or hind limbs. The behaviors were recorded with a camera from the side. The total duration of immobility (no movement except keeping the head above water) in the last 4 min was counted.

The sucrose preference test (SPT) was conducted according to published protocols (Monteggia et al., 2007; Chaudhury et al., 2013; Shen et al., 2019) with minor changes. Animals were single housed and habituated with two bottles of water for two consecutive days. After habituation, the mice were allowed free access to a two-bottle choice of 2% sucrose solution or tap water. The positions of the two bottles were switched (left to right and right to left) every 6 h to avoid side preference. The total duration of the test was 48 h. The bottles were weighed at the start and end of each testing period. Sucrose preference was calculated as sucrose consumption/(water consumption + sucrose consumption) × 100% during the test phase.

All experimental data were analyzed with GraphPad Prism 8.0 (San Diego, CA, United States). All data are presented as means ± standard error of the mean (SEM). We used a two-tailed unpaired t-test for two-group comparisons and one-way ANOVA followed by Bonferroni multiple testing for more than two-group comparisons. Two-way ANOVA was applied to compare data from the behavioral tests after CSDS with two factors (Genotype × Stress). A p-value of <0.05 was considered statistically significant.

The characteristics of the study population are presented in Supplementary Table 1. A total of 4,817 subjects with an average age of 42.8 years (± 13.5; SD) were included. Of them, 3,279 (68.1%) were self-reported African Americans, and 45.2% were males. In line with the population-based design of the study, the mean depressive symptom scores were 6.3 (± 9.5; SD) for the 20-item version of the CES-D.

Table 1 lists the association results of the top 11 SNPs with p < 1 × 10^–5, among which two loci showed significant association with the CES-D score after Bonferroni correction: rs61753730 (Chr. 8: 104336788; MAF = 0.012; p = 8.46 × 10^–9) in FZD6; and rs35024632 (Chr. 7: 38433726; MAF = 0.025; p = 3.08 × 10^–8) in the amphiphysin (AMPH) gene. Figure 1A shows the Manhattan plot of the association results by chromosome location, and two SNPs, located in chromosomes 7 and 8, showed signal above the commonly defined genome-wide significance threshold (p < 5 × 10^–8) (Sullivan et al., 2009). Of these two SNPs, both are missense variants resulting in amino acid changes in their encoded proteins (rs61753730 from glutamine to glutamate; rs35024632 from lysine to threonine). Figure 1B depicts the quantile-quantile plot of observed vs. expected –log₁₀ (p- value) for all 38,943 SNPs. The distribution of p- values was similar to what would be expected by chance, with a genomic control λ value of 1.000, indicating no detectable inflation of test statistics (Devlin and Roeder, 1999).

To characterize the most significantly associated SNP rs61753730 in FZD6, we performed gene-based association analysis for 11 genotyped and 121 imputed SNPs (info > 0.3) located in the gene as well as in a 2-kb region both upstream and downstream, which showed a significant association of FZD6 with the CES-D score (p = 0.0078). Of these SNPs, 102 were rare and 30 were common variants.

To further investigate the biological role of rs61753730 on FZD6 in depression, we carried out in silico analysis and found that the different alleles of rs61753730 lead to an altered protein configuration (Figure 2A). Although no obvious difference was observed in the physiochemical characteristics of the protein and disulfide bonds, several distinctions were discovered in the secondary structure and transmembrane helices (Supplementary Table 2). On the basis of the prediction by NetSurfP-2.0, the secondary structure was quite different at position 152 of FZD6 (Figure 2B), which was a coil in the WT (Gln, Q) but a helix in the mutant FZD6 (Glu, E). We further predicted the stability of the protein after mutation by I-Mutant 3.0 and found that the ΔΔG value was −0.18 kcal/mol, which means the mutation decreased the stability of the FZD6 protein because of the amino acid change from Gln (rs61753730-C) to Glu (rs61753730-G) at position 152. These findings were further supported by another structure-based prediction tool mCSM with an online server DUET (ΔΔG = −0.039 kcal/mol). Together, these in silico results indicated that the different allele of rs61753730 changed the structure of FZD6.

In determining the allele-specific function of rs61753730 on FZD6 experimentally, HEK293T cells transfected with the rs61753730-G allele construct showed a 20% decrease in luciferase activity (p = 0.014) compared with the cells transfected with the rs61753730-C allele (Figure 2C). Further, when the cells were treated with actinomycin D, the mRNA quantity of the luciferase reporter gene was significantly reduced in both C and G allele-transfected cells, whereas the construct carrying the rs61753730-C allele decayed significantly more slowly than that containing the G allele, with a p-value of 0.0349 and 0.0267 after 1 and 3 h, respectively, of actinomycin D treatment (Figure 2D).

By analyzing the expression quantitative trait loci (eQTL) dataset from GTEx Portal, we found that rs61753730 showed eQTL regulation with allele-specific significance on FZD6 mRNA expression, especially in the cerebellum (p = 1.3 × 10^–5), a brain region implicated in depression (Depping et al., 2018). These results demonstrated that the missense mutation rs61753730 contributed to the effect of FZD6 on depression.

As shown in Supplementary Figures 1A,B, the glutamine residue at position 152 (rs61753730-C) was highly conserved from human to mouse and rat, suggesting it has functional importance, and the mouse model could be used to investigate the biological role of rs61753730 in depressive symptoms. Supplementary Figure 1C shows the strategy used to generate the rs61753730 mutant Fzd6-KI mice with the CRISPR/Cas9 system. Since rs61753730 is located in the fourth exon of FZD6, the CRISPR/Cas9 system targeted to it to obtain the Fzd6-KI mice carrying G allele. The sequence was validated by Sanger sequencing with the C allele in WT mice and the G allele in Fzd6-KI mice (Supplementary Figure 1D).

To characterize the role of rs61753730 in depression in the mouse model, the 10-day CSDS-induced depression model, a well-validated and widely used model (Golden et al., 2011), was employed (Figures 3A,B). After the 10-day CSDS, the social interaction test was performed (Figures 3C–E). When there was no target, the two-way ANOVA showed no significant interaction for Stress × Genotype (F_(1,13) = 0.04783, p = 0.8303), and both the Genotype factor (F_(1,13) = 0.3138, p = 0.5849) and Stress factor (F_(1,13) = 2.643, p = 0.1280) also showed no significant difference, respectively, which means the baseline is similar for the four groups. When the target present, although the results showed no significant interactions for Stress × Genotype (F_(1,13) = 2.212, p = 0.1608) and the Genotype factor (F_(1,13) = 2.212, p = 0.1608), the stress factor showed a close to significant difference (F_(1,13) = 4.353, p = 0.0572). As shown in Figure 3E, in the presence of target mice, the WT mice who subjected to CSDS spent less time in the interaction zone than the WT mice without CSDS (p < 0.05), and the Fzd6-KI mice spent less time in the interaction zone compared with the WT mice (p < 0.05). No significant difference was observed between Fzd6-KI mice whether subjected to CSDS or not (p > 0.05). It suggests that the Fzd6-KI mice is not susceptible to stress, however, this type of mice still can present the similar depressive-like symptoms similar with WT-CSDS mice. Thus, more depressive-like and anxiety-like behaviors need to be conducted to verify this phenomenon.

To identify the depressive-like behaviors in Fzd6-KI mice, we carried out FST and SPT experiments and used the measurement of immobility duration and sucrose preference ratio to evaluate anhedonia. In the FST, two-way ANOVA showed significance for Stress factor (F_(1,13) = 9.415, p = 0.0090), but not for Genotype factor (F_(1,13) = 1.631, p = 0.2239) and for their interactions (F_(1,13) = 3.312, p = 0.0919). Compared with the WT control mice, the immobility duration increased significantly (p < 0.05; Figure 4A) after the WT mice were subjected to CSDS. The Fzd6-KI mice exposed to CSDS also showed significantly greater immobility compared with the WT control mice (p < 0.05), whereas no significant difference was observed between the Fzd6-KI mice with and without CSDS exposure. In the SPT, two-way ANOVA showed that significant difference for Genotype factor (F_(1,12) = 7.452, p = 0.0183), but not for the Stress factor (F_(1,12) = 3.088, p = 0.1044) and their interactions (F_(1,12) = 0.08952, p = 0.7699). Although no significant difference in the sucrose preference ratio was detected between the WT mice exposed and not exposed to CSDS (Figure 4B), the results showed the same trend; and Fzd6-KI mice exposed to CSDS displayed significantly less sucrose preference than the WT control mice (p < 0.05).

As anxiety-like behaviors correlated highly with depression, we measured these behaviors in Fzd6-KI mice with the OFT, which was measured by the number of times the mice entered or spent time in the central area. Two-way ANOVA showed significant differences for Stress × Genotype interactions for total distance traveled (F_(1,13) = 8.376, p = 0.0125) and center entries (F_(1,13) = 6.317, p = 0.0259), but not for time spent in the center zone (F_(1,13) = 3.289, p = 0.0929). As shown in Figure 4C, the locomotor activities with total distance traveled (Figure 4D), center entries (Figure 4E), and time spent in the center zone (Figure 4F) were all greatly decreased after the WT mice were subjected to CSDS (p < 0.01), and the Fzd6-KI mice with or without CSDS exposure showed a comparable reduction. In Fzd6-KI mice exposed to CSDS, it also presented a significant decrease in locomotor activities compared with the WT control mice (p < 0.05), but there was no significant difference in Fzd6-KI mice with or without CSDS exposure. Taken together, these behavioral results strongly indicated that the Fzd6-KI mice presented depressive symptoms like those in WT mice subjected to CSDS, suggesting an important role of rs61753730 in depression.