Functional differences between TSHR alleles associate with variation in spawning season in Atlantic herring
AnimalsTissue samples for expression analysis were collected on August 24 and September 9 2016 from a population of spring-spawning herring kept in captivity at University of Bergen; the rearing of captive herring was approved by the Norwegian national animal ethics committee (Forsøksdyrutvalget FOTS ID-5072). The tissue samples used for ATAC-sequencing were collected at Hästskär on June 19 2019 from wild spring-spawning Atlantic herring in the Baltic Sea, which do not require ethical permission.Genome scan and genetic diversity at the TSHR locusWe used a 2 × 2 contingency X2 test to estimate the extent of SNP allele frequency differences between seven spring- and seven autumn-spawning herring populations from the Northeast Atlantic (Supplementary Table 2), and thus, identify the major genomic regions associated with seasonal reproduction. The SNP allele frequencies were generated in a previous study16 and were derived from Pool-seq data. For the X2 test, we formed two groups, spring and autumn spawners, and summed the reads supporting the reference and the alternative alleles for the pools included in each group.To characterize genetic diversity at the TSHR locus, we calculated nucleotide diversity (π) and Tajima’s D for the same seven spring- and seven autumn-spawning Atlantic herring populations used in the genome scan (n = ~41–100 per pool) (Supplementary Table 2). The whole-genome re-sequence data of these pools were previously reported by Han et al.16 (for details of the pools used here see Supplementary Table 2). Unbiased nucleotide diversity π and Tajima’s D were calculated for each pool using the program PoPoolation 1.2.236, which accounts for the truncated allele frequency spectrum of pooled data. In brief, a pileup file of chromosome 15 was generated from each pool BAM file using samtools v.1.1037,38. Indels and 5 bp around indels were removed to exclude spurious SNPs due to misalignments around indels. To minimize biases in the π and Tajima’s D calculations, which are sensitive to sequencing errors and coverage fluctuations39, the coverage of each pileup file was subsampled without replacement to a uniform value based on a per-pool coverage distribution (the target coverage corresponded to the 5th percentile of the coverage frequency distribution, which was used as the minimum coverage allowed for an SNP to be included in the analysis) (Supplementary Table 2). We also calculated the diversity parameters but skipping the coverage subsampling step and obtained very similar results with both approaches (Supplementary Fig. 5), thus, we decided to keep working with the subsampled datasets as coverage subsampling is recommended by the software developers36. To exclude spurious SNPs associated with repetitive sequences and copy number variants, we applied a maximum coverage filter equivalent to the per-pool 99th percentile of the coverage frequency distribution (Supplementary Table 2). A minimum base quality of 20, a minimum mapping quality of 20, and a minor allele count of 2 were required to retain high quality SNPs for further analysis. Both, π and Tajima’s D statistics were calculated using a sliding window approach with a window size of 10 kb and a step size of 2 kb (the selected window-step combination offered a good genomic resolution while reducing the noise from single SNPs, after testing windows of 5, 10, 20, 40, 50, and 100 kb for non-overlapping and overlapping windows with a step size equivalent to 20% of the window size). Only windows with a coverage fraction ≥ 0.5 were included in the computations. In addition, we estimated the effective allele frequency difference, or delta allele frequency (dAF), between spring and autumn spawning groups at the TSHR locus using the formula dAF = abs(mean(spring pools)−mean(autumn pools)). For each of the diversity parameters, we assessed whether the mean differences between sets of SNPs within chr 15: 8.85–8.95 Mb (215 SNPs) and outside (214 635 SNPs) the TSHR region were statistically significant among spring- and autumn-spawning groups using a Wilcoxon test. Data postprocessing, statistical tests, and plotting were performed in the R environment40 (for specific parameters used in PoPoolation, see the associated code to this publication).Identification of the 5.2 kb structural variantSequences spanning the entire TSHR gene plus 10 kb upstream and downstream from two reference assemblies, ASM96633v115 and Ch_v2.0.218, were aligned using BLAST41 and the output was subsequently processed with a custom R script40. Repeats were annotated by CENSOR21 for the region harboring the 5.2 kb structural variant. To validate the structural variant, long-range PCR was performed with genomic DNA from two spring- and autumn-spawning Atlantic herring in a 20 μL reaction containing 0.8 mM dNTPs, 0.3 μM each of the forward and reverse 5.2kb-confirm primers (Supplementary Table 1) and 0.75 U PrimeSTAR GXL DNA Polymerase (TaKaRa) following the program: 95 °C for 3 min, 35 cycles of 98 °C for 10 s, 58 °C for 20 s and 68 °C for 2 min 30 s, and a final extension of 10 min at 68 °C.ATAC-seq analysisBSH and brain without BSH were dissected from two spring-spawning herring caught in the Baltic Sea and transported to the lab on dry ice, then kept at –80 °C before nuclei isolation. ATAC-seq libraries were constructed according to the Omni-ATAC protocol with minor modifications42. Briefly, tissue was homogenized in 2 ml homogenization buffer (5 mM CaCl2, 3 mM Mg(Ac)2, 10 mM Tris-HCl (pH = 7.8), 0.017 mM PMSF, 0.17 mM ß-mercaptoethanol, 320 mM Sucrose, 0.1 mM EDTA and 0.1% NP-40) with a Dounce homogenizer on ice. 400 μl suspension was transferred to a 2 ml tube for the density gradient centrifugation with different concentrations of Iodixanol solution. After centrifugation, a 200 μl fraction containing the nuclei band was collected, stained with Trypan blue and counted with a Countess II Automated Cell Counter (Thermo Fisher Scientific). An aliquot of 100,000–200,000 nuclei was used as input in a 50 μl transposition reaction containing 2 X TD buffer and 100 nM assembled Tn5 transposase for a 30-min incubation at 37 °C. Tagmented DNA was purified with a Zymo clean kit (Zymo Research). Purified DNA was used for an initial pre-amplification for 5 cycles, and the additional amplification cycle was determined by qPCR based on the “R vs Cycle Number” plot43. Amplified libraries were purified with a Zymo clean kit again, and library concentrations and qualities were evaluated using the 2200 TapeStation System (Agilent Technologies).ATAC-seq was performed with a MiniSeq High Output Kit (150 cycles) on a MiniSeq instrument (Illumina) and 7–9 million reads were generated for each ATAC-seq library. Quality control, trimming, mapping, and peak calling of the sequenced reads were conducted following the ENCODE ATAC-seq pipeline (https://www.encodeproject.org/atac-seq/). The trimmed reads were aligned to the Atlantic herring reference genome (Ch_v2.0.2)18 with Bowtie244 and the mapping rate was 85–95%. Duplicate reads, reads with low mapping quality and those aligned to the mitochondria genome were removed. The remaining reads (4–5 million) were subjected to peak calling by MACS245, where 22–32 K peaks were called. Sequenced library qualities were further evaluated by calculating the TSS enrichment score and checking the library complexity with the Non-Redundant Fraction (NRF) and PCR Bottlenecking Coefficients (PBC1 and 2). Finally, conserved peaks between two biological replicates were identified by evaluating the irreproducible discovery rate (IDR).Genotyping of six differentiated variants and haplotype analysisAll six genetic variants, including the 5.2 kb structural variant, two non-coding SNPs, two missense SNPs and the copy number variant of C-terminal 22aa repeat, were genotyped in 45 spring-, 67 autumn-spawning Atlantic herring and 13 Pacific herring. TaqMan Custom SNP assays were performed to genotype the four SNPs in 5 μl reactions with a template of 20 ng genomic DNA (ThermoFisher Scientific). Copy number of the C-terminal 22aa repeat was determined by the PCR product size generated with geno22aa primers. Genotyping of the 5.2 kb structural variant was performed in a PCR reaction containing two forward primers (geno5.2kb-1F and geno5.2kb-2F) and one reverse primer (geno5.2kb-R), which generated PCR products with different sizes between spring and autumn spawners. All the primers used for genotyping are listed in Supplementary Table 1.Tissue expression profiles by quantitative PCRTotal RNA was prepared from gonad, heart, spleen, kidney, gills, intestine, hypothalamus and saccus vasculosus (BSH), and brain without BSH (brain) of six adult spring-spawning Atlantic herring using RNeasy Mini Kit (Qiagen). RNA was then reverse transcribed into cDNA with a High-Capacity cDNA Reverse Transcription Kit (ThermoFisher Scientific). TaqMan Gene Expression assay (ThermoFisher Scientific) containing 0.3 μM primers and 0.25 μM TaqMan probe (Integrated DNA Technologies) was performed to compare the relative expression levels of TSHR among different tissues. qPCR with SYBR Green chemistry was used for TSHB and DIO2 in a 10 μl reaction of SYBR Green PCR Master Mix (ThermoFisher Scientific) and 0.3 μM primers, with a program composed of an initial denaturation for 10 min at 95 °C followed by 40 cycles of 95 °C for 15 s and 60 °C for 1 min. Ct values were first normalized to the housekeeping gene ACTIN, then the average expression for each gene in the gonad was assumed to be 1 for the subsequent calculation of the relative expression in other tissues.Plasmid constructsThe coding sequence for the herring single-chain TSH (scTSH) was designed following a strategy previously used for mammalian gonadotropins46 that contained an in-frame fusion of the cDNA sequences (5′–3′) of herring TSH beta subunit (NCBI: XM_012836756.1) and alpha subunit (NCBI: XM_012822755.1) linked by six histidines and then the C-terminal peptide of the hCG beta subunit. The designed sequence should generate a protein with a size of 30.6 kDa. Both scTSH and spring herring TSHR cDNA sequences were synthesized in vitro and cloned in the expression vector pcDNA3.1 by Genscript (Leiden, Netherlands). pcDNA3.1 plasmid expressing human TSHR was kindly provided by Drs. Gilbert Vassart and Sabine Costagliola (Université libre de Bruxelles, Belgium). Then, the spring herring TSHR and human TSHR plasmids were used as templates for site-directed mutagenesis to generate constructs coding for different mutant herring or human TSHRs. Plasmids for the dual-luciferase assay, including pGL4.29[luc2P/CRE/Hygro] containing cAMP response elements (CREs) to drive the transcription of luciferase gene luc2P and pRL-TK monitoring the transfection efficiency, were purchased from Promega. Five ng of each plasmid was used to transform the XL1-Blue competent cells (Agilent), plasmid DNA was subsequently extracted from 200 ml overnight culture of a single transformant clone using an EndoFree plasmid Maxi Kit (Qiagen).Cell cultureChinese hamster ovary (CHO) (ATCC CCL-61) and human embryonic kidney 293 (HEK293) (ATCC CRL-1573) cells were maintained in DMEM supplemented with 5% (CHO) or 10% FBS (HEK293), 100 U/ml penicillin, 100 μg/ml streptomycin and 292 μg/ml l-Glutamine (ThermoFisher Scientific) at 37 °C with 5% CO2. Epithelioma Papulosum Cyprini (EPC) cells (ATCC CRL-2872) were cultured in EMEM (Sigma) supplemented with 10% FBS, 100 U/ml penicillin, 100 μg/ml streptomycin, 292 μg/ml l-Glutamine and 1 mM Sodium Pyruvate (ThermoFisher Scientific) at 26 °C with 5% CO2.Production of recombinant herring scTSHCHO cells were transfected with the scTSH expression plasmid using Lipofectamine 3000 (Invitrogen), stable clones were subsequently selected with 500 μg/ml G418 (Invitrogen) and screened for producing scTSH by western blot using a polyclonal antisera against the sea bass alpha subunit47. A positive clone was expanded in 225 cm2 cell culture flasks (Corning) in culture medium containing 5% FBS until confluence, then the cells were maintained in serum-free DMEM for hormone production for 7 days at 25 °C48. After 7 days, culture medium containing scTSH or without (negative control) was centrifuged at 15000 x g for 15 min and concentrated by ultrafiltration using Centricon Plus-70 / Ultracel PL-30 (Merck Millipore Ltd.). Then, western blotting was performed to confirm TSH production. Concentrated medium containing herring scTSH was denatured at 94 °C for 5 min in 0.1% SDS and 50 mM 2-mercaptoethanol, and then treated with 2.5 units of peptide-N-glycosidase F (Roche Diagnostics) at 37 °C for 2 h in 20 mM sodium phosphate with 0.5% Nonidet P-40, pH 7.5. All samples were run in 12% SDS-PAGE in the reducing condition and transferred to a PVDF membrane (Immobilon P; Millipore Corp.), then blocked overnight with 5% skimmed milk at 4 °C. After blocking, the membrane was incubated with polyclonal antisera against the sea bass alpha subunit (dilution 1:2000) for 90 min at room temperature, washed, and then further incubated with 1:25000 goat anti-rabbit immunoglobulin G (IgG) horseradish peroxidase conjugate (Bio-Rad Laboratories) for 60 min at room temperature. Immunodetection was performed by chemiluminescence with a Pierce ECL Plus Western Blotting Substrate kit (ThermoFisher Scientific).Cell surface expressionA Rhotag (MNGTEGPNFYVPFSNKTGVVYEE) was inserted at the N-terminus of herring TSHR for flow cytometry analysis of receptor cell surface expression. Anti-Rhotag polyclonal antibody was kindly provided by Drs. Gilbert Vassart and Sabine Costagliola (Université Libre de Bruxelles, Brussels, Belgium). PBS containing 1% BSA and 0.05% sodium azide was prepared as the flow cytometry (FCM) buffer for the washing and antibody incubation steps. 2.2 × 106 EPC cells were seeded in a 100 mm poly-d-Lysine-treated petri dish the day before transfection. Each dish was transfected with 10 μg TSHR or empty pcDNA3.1 expression plasmid using 20 μl jetPRIME transfection reagent in 500 μl jetPRIME transfection buffer (Polyplus transfection). Cells were harvested 24 h after transfection, then washed once in cold PBS and fixed in 2% PFA for 10 min at room temperature. After fixation, cells were washed three times with FCM buffer, then incubated with anti-Rhotag antibody or FCM buffer (negative control) for 1 h at room temperature. Cells were washed again with FCM buffer three times and stained with Alexa Fluor 488-labeled chicken anti-mouse IgG (H + L) antibody (1:200 dilution, ThermoFisher Scientific) or FCM buffer (negative control) for 45 min in the dark. After the fluorescent staining, cells were washed three times and resuspended in FCM buffer before analysis on a CytoFLEX instrument (Beckman Coulter). A minimum of 100,000 events was recorded for each sample, fluorescence intensities of negative control and cells transfected with empty pcDNA3.1 plasmid were used as the background for gating strategy. Cell surface expression was represented by the mean fluorescence intensity of the positively stained cell population.Dual-luciferase reporter assayEPC or HEK293 cells were plated in a 48-well plate at a density of 1 × 105 cells/well the day before transfection. A total of 250 ng plasmid mixture containing pGL4.29[luc2P/CRE/Hygro], TSHR expression plasmid (or empty pcDNA3.1) and pRL-TK with the ratio of 20:5:1 was prepared to transfect each well of cells using jetPRIME transfection reagent (Polyplus). Medium was replaced by fresh medium containing 10% FBS (TSH-induced condition) or serum-free medium (constitutive activity condition) 4 h after transfection. On day three, cells were treated with serum-free medium containing different dilutions of the concentrated scTSH medium for 4 h (TSH-induced condition) or directly subjected to the luminescence measurement without TSH induction (constitutive activity condition). Luminescence was measured using a Dual-Luciferase Reporter assay (Promega) on an Infinite M200 Microplate Reader (Tecan Group Ltd., Switzerland), and luciferase activity was represented as the ratio of firefly (pGL4.29[luc2P/CRE/Hygro]) to Renilla (pRL-TK) luminescence.5′-RACE to identify the herring DIO2 TSSTotal RNA was prepared from brain of a spring-spawning Atlantic herring using the RNeasy Mini Kit (Qiagen). Six μg of the isolated RNA was used for 5′-RACE with a FirstChoiceTM RLM-RACE Kit (ThermoFisher Scientific). One μl cDNA or Outer RACE PCR product was used as PCR template in a 20 μL reaction containing 0.8 mM dNTPs, 0.3 μM of each forward and reverse primer (Supplementary Table 1) and 0.75 U PrimeSTAR GXL DNA Polymerase (TaKaRa). Amplification was carried out with an initial denaturation of 3 min at 95 °C, followed by 35 cycles of 98 °C for 10 s, 58 °C for 20 s and 68 °C for 40 s, and a final extension of 10 min at 68 °C. The final 5′ RACE product was sequenced at Eurofins Genomics (Ebersberg, Germany).Sequence conservation analysisGenomic sequences covering the TSHR locus were extracted from Ensembl Genome Browser for Atlantic herring and 11 other fish species, including Amazon molly (Poecilia formosa), denticle herring (Denticeps clupeoides), goldfish (Carassius auratus), guppy (Poecilia reticulata), Neolamprologus brichardi, Japanese medaka (Oryzias latipes), northern pike (Esox lucius), orange clownfish (Amphiprion percula), spotted gar (Lepisosteus oculatus), three-spined stickleback (Gasterosteus aculeatus) and spotted green pufferfish (Tetraodon nigroviridis). The extracted sequences were firstly aligned using progressiveCactus49,50, and a subsequent alignment was generated using the hal2maf program from halTools51 with Atlantic herring assembly (Ch_v2.0.2)18 as the co-ordinate backbone. This alignment was used for the downstream phastCons score calculation by running phyloFit24 and phastCons25 from the PHAST package with default parameters. Peaks were called by grouping signals with a minimum phastCons score of 0.2 within 500 bp region.Structure modeling of human and herring TSHRsIn order to explore the possible interactions of the variant residues with other receptor interacting proteins and to study intramolecular interactions, we built a structural homology model for the herring TSHR (herrTSHR) complexed with herring TSH and Gs-protein. The TSHR hinge region that harbors the Q370H substitution and the C-terminus containing the 22aa repeat were excluded from the homology model due to the lack of structural templates for these regions. The homology model was constructed by using the following structural templates of evolutionarily related class A GPCRs: (i) the leucine-rich repeat domain (LRRD) complexed with hormone was modeled based on the solved FSHR LRRD – FSH complex structure (Protein Data Bank (PDB) ID: 4AY9)52,53, this part of model included herring TSHR Cys33 – Asn296 and fragments of the hinge region Gln297 – Thr312 and Ser393 – Ile421; (ii) the available structural complex of β2-adrenoreceptor with Gs-protein (PDB ID: 3SN6)54 was used as the template to model the seven-transmembrane helix domain (7TMD) of herring TSHR in the active conformation; (iii) the extracellular loop 2 (ECL2) was built by using the ECL2 of μ-opioid receptor (PDB ID: 6DDE)55. To prepare the template for herring TSHR modeling, the fused T4-lysozyme and bound ligand of β2-adrenoreceptor were deleted, the ECL1 and ECL3 loops were adjusted manually to the loop length of herring TSHR. Due to the lack of third intracellular loop (ICL3) in the β2-adrenoreceptor structure, amino acid residues of herring TSHR ICL3 were manually added to the template. Since herring TSHR does not have the TMH5 proline, which is highly conserved among all class A GPCRs and responsible for the helical kinks and bulges within this region56, we assumed a rather regular (stretched) helix conformation for the herring TSHR TMH5 and therefore replaced the kinked β2-adrenoreceptor TMH5 template with a regular α-helix. Moreover, the ECL2 template was substituted with μOR ECL2 structure because of its higher sequence similarity with herring TSHR in this region. Finally, amino acid residues of this chimeric 7TMD template and FSHR N-terminus were mutated to the corresponding spring herring TSHR residues and sequence of the heterodimeric FSH ligand was substituted by the herring TSH. All homology models were generated by using SYBYL-X 2.0 (Certara, NJ, US). The 7TMD structure was then fused with FSHR N-terminus at position 421. The assembled complex was subsequently optimized by the energy minimization under constrained backbone atoms (the AMBER F99 force field was used), followed by a 2 ns molecular dynamics simulation (MD) of the side chains. The entire TSHR complex was energetically minimized without any constraints until converging at a termination gradient of 0.05 kcal/mol*Å. Next, for autumn herrTSHR modeling, the spring TSHR sequence was substituted with autumn TSHR sequence. For humTSHR, the spring TSHR sequence was substituted with human TSHR, and the herring TSH ligand was replaced by the bovine TSH sequence. Both complex models were energetically minimized until converging at a termination gradient of 0.05 kcal/mol*Å.To investigate the microenvironment around the L471M mutation at TMH2 position 2.51, local short MD’s of 4 ns on Met4712.51 (spring herrTSHR), Leu471 (autumn herrTSHR) or Phe461 (humTSHR) and its surrounding amino acids were performed. During MD simulations, backbone atoms of the entire complexes as well as all side chains, except residues at positions 1.47, 1.51, 1.54, 2.48, 2.52, and 2.55 that form the hydrophobic patch around position 2.51, were constrained.Statistics and reproducibilityResults were presented as the mean + SD (standard deviation) calculated from at least four biological replicates for each experiment, and at least two independent experiments were conducted for each assay. Unpaired two-tailed Student’s t test was performed to calculate the P-values and means were judged as statistically significant when P ≤ 0.05.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More
