Conserved exchange of paralog proteins during neuronal differentiation

Ensembl Compara paralog genes resources

Paralogs annotation for Homo sapiens (GRch38.p13) Danio rerio (GRCz11) Mus musculus (GRCm38.p6) Rattus norvegicus (Rnor_6.0), were downloaded from Ensembl (v102) via biomart (http://www.ensembl.org/biomart/martview/f04b3aa8b5c7f463e3edf9fa58d205a7). Duplicated paralog pairs (e.g., Paralog₁ | Paralog₂; Paralog₂ | Paralog₁) were removed from each dataset, so that only unique pairs (Paralog₁ | Paralog₂) were retained.

Protein complexes resources

Protein Complexes definition were taken from Ori et al (2016). Members of protein complexes were mapped by orthology in D. rerio and R. Norvegicus using the bioconductor package “biomaRt” (Durinck et al, 2009) using as reference the H. sapiens protein complexes definitions.

Publicly available data used in this study

Zebrafish embryo development data were obtained from White et al (2017) (Table S1). Human Proteome and Transcriptome data across tissues were obtained from Wang et al (2019) (Table S2). For these specific datasets, multi-mapping between reads was handled as it follows. For zebrafish embryo development data, the author used htseq-count to assign reads to its specific transcript. For this dataset reads that map to multiple genes were discarded. For human tissue atlas, Cufflinks v2.1.1 was used to assign reads to different transcript. In these cases multi-mapped reads were uniformingly split between genes. Protein identification and LFQ intensity values (Log₂) in cultured human iPSCs, NPCs and differentiated neurons, were obtained from Table S2 from Djuric et al (2017) and Frese et al (2017). Rat neuronal differentiation data published in Frese et al (2017), were downloaded from PRIDE (http://proteomecentral.proteomexchange.org/cgi/GetDataset?ID=PXD005031) and analyzed again as described below.

Early embryos (6 hpf)

Wild type embryos were removed from their chorions using 1 ml of pronase (stock 30 mg/ml) in 40 ml buffer E3 and incubated for 10–15 min with gentle shaking every 2 min in a small beaker. The supernatant was removed and the embryos were washed four to five times using buffer E3. The embryos were splitted into batches of around 250–300 per 1.5 ml tube. 1 ml of deyolking buffer (55 mM NaCl, 1.8 mM KCl, 1.25 mM NaHCO₃) was added per tube and everything passed twice through a 200 μl pipet tip. The tubes were incubated at RT in a shaker at 1,100 rpm for 5 min and afterwards spun at 300g for 30 s to remove the supernatant. The embryos were washed using 1 ml of wash buffer (110 mM NaCl, 3.5 mM KCl, 2.7 mM CaCl₂ and 10 mM Tris/HCl, pH 8.5), shaken at 1,100 rpm at RT for 2 min and spun as above to remove the supernatant. The wash step was repeated twice. The deyolked and dissociated embryos were resuspended in 400 μl wash buffer and passed through a 40 μm cell strainer to remove undissociated cells. The merged cells were washed as above and resuspended in 110 μl PBS and counted using a hemocytometer.

Late embryos (24 hpf)

After 24 hpf, the NBT dsRed positive embryos were manually sorted. Up to the addition of the deyolking buffer all steps were the same as for early embryos. After the addition of 1 ml deyolking buffer per tube, the embryos were passed 10 times through a 1,000 μl pipet tip, followed by washing twice with deyolking buffer and four times with washing buffer. For better cell dissociation the embryos were rinsed once with Accumax (Millipore) and then resuspended in 1 ml Accumax and transferred to a 15 ml tube. The embryos were incubated at RT for 5 min at the lowest speed of the vortex mixer. The embryos were dissociated by pipetting for 2 min using a 1,000 μl pipet tip, 2 min incubation on the vortex mixer and 1 min of additional pipetting. The cells were spun for 1 min at 300g at RT and washed twice using 1 ml of PBS with 0.5% BSA. 400 μl of PBS with 0.5% BSA was used per tube to resuspend the cells afterwards passed through a 40 μm cell strainer and merged. DNAse I (10 mg/ml in water; Roche) 170 U/ml and 10 mM MgCl₂ was added. Cells expressing the DsRed fluorescent protein were FAC sorted with a MoFlo cell sorter (Beckman Coulter GmbH) to obtain a highly enriched fraction for neuronal cells.

Animal management practices

All mice were maintained in specific pathogen-free conditions, with food and water available ad libitum. The animal room had a constant temperature of 21°C ± 2°C, 55% ± 15% humidity, and controlled lighting (12 h light/dark cycle). The location for animal keeping was animal house TH4 at Leibniz Institute on Aging (Fritz Lipmann Institute). Breeding was license-free and performed under §11 TierSchG. Euthanasia and organ removal were performed under the internal §4 TierSchG licences O_CK_18-20 and O-CK_21-23. Euthanasia of mice was performed in a chamber with controlled CO2 fill rate according to “Directive 2010/63/EU annex IV of the European Parliament and the Council on the protection of animals used for specific purposes.”

Mouse neuronal cell culture

Cortical neurons were isolated from wild type murine embryonic brains (E15.5) of mixed background (FVB/NJ, C57BL/6, 129/Sv) and differentiated in glia-conditioned neurobasal medium. Briefly, meninges were removed, cortices were isolated, minced and dissociated in trypsin EDTA (Invitrogen), solution for 15 min at 37°C. The supernatant was removed and the tissue was washed three times with trituration solution (10 mM Hepes, 1% penicillin/streptomycin, 10 mM L-glutamine, 1% BSA, 10% FBS, and 0.008% DNase in HBSS) and homogenized in trituration solution using fire polished glass pipettes. For the mouse in vitro neuronal differentiation data, neurons were counted and pellets containing 1 million cells (DIV0) were prepared and frozen until further use. In addition, 1 million cells were seeded on poly-L-lysine–coated 6-cm plates containing 4 ml glia-conditioned plating medium (1% penicillin/streptomycin, 1 mM sodium pyruvate, 0.5% glucose, 10 mM Hepes 1× B27 supplement, 10% FBS, and 10 mM L-glutamine in MEM). After 24 h, the plating medium was substituted by glia-conditioned neurobasal medium (10 mM Hepes, 1× B27 supplement, 5 mM L-glutamine in NBM). Neurons were collected at DIV3 and DIV10. To this end, neurons were scraped off in cold PBS and obtained cell suspensions were transferred to a microcentrifuge tube and centrifuged for 5 min at 4°C and 500g. The obtained pellets were washed with PBS twice and frozen until further use.

For preparation of glia-conditioned mediums, a primary astroglia culture was established. For this purpose, brains were isolated from 15.5 d old embryos, the meninges were removed, the cerebral hemispheres were minced and afterwards dissociated in trypsin solution for 15 min at 37°C. Finally, the tissue was homogenized by pipetting and cells were plated on a 10 cm dish containing glia medium (1% penicillin/streptomycin, 1 mM sodium pyruvate, 0.5% glucose, 10 mM Hepes, 20 mM L-glutamine, and 10% FBS in MEM) and grown to confluence. For preconditioning of neurobasal medium or plating medium, the media were added to the glia feeder cultures and collected after 24 h.

Sec23a and Sec23b knockdown in mouse neuronal differentiation

For the Sec23 paralogs knockdowns, cortical neurons were isolated from C57BL/6JRj mouse embryo (Janvier), as described above. Then, freshly isolated neurons (5 million cells per nucleofection reaction) were transfected using the 4D-Nucleofector™ X Unit and the P3 Primary Cell 4D Nucleofector X kit (Lonza), as indicated. Cells were transfected with 250 nM of siRNA and 1 μl of control pMax GFP (Nucleofector X kit; Lonza), using the CU-133 program. Immediately after transfection, cells were plated on poly-L-lysine (Sigma-Aldrich)–coated 10-cm plates containing 10 ml of glia-conditioned plating medium: 1% penicillin/streptomycin, 1 mM sodium pyruvate (Sigma-Aldrich), 0.5% glucose, 10 mM Hepes, 1× B27 supplement (Invitrogen), 10% FBS, 10 mM L-glutamine in MEM (31095-052; Invitrogen), and incubated at 37°C. After 1 d, the medium was replaced with glia-conditioned neurobasal medium: 10 mM Hepes, 1× B27 supplement, 5 mM L-glutamine in NBM (Invitrogen). After 3 d in culture, neurons were washed twice with PBS, detached using Trypsin EDTA (3–5 min, 37°C), collected in 5 ml of PBS with 2% FBS, and pelleted by centrifugation (450g, 8 min, room temperature). Pellets were resuspended in 0.3 ml PBS with 2% FBS, and GFP-positive cells were labeled with Sytox Blue Dead Cell Stain (viable staining) (Molecular Probes; Thermo Fisher Scientific) and sorted directly in 200 μl of 2× lysis buffer (200 mM Hepes pH 8.0, 100 mM DTT, and 4% SDS) using a BD FACSAria Fusion with the Software BD FACSDiva 8.0.1 and 9.0.1 (BD Biosciences), using 488 nm laser and 530/30 filter for the GFP signal and laser 405 nm and 450/50 filter for the Sytox blue.

Sample preparation and dimethyl labeling for zebrafish stem cells and neurons

Cells were lysed by addition of Rapigest (Waters) and urea to a final concentration of 0.2% and 4 M, respectively, and sonicated for 3 × 30 s to shear chromatin. Before protein digestion, samples were stored at −80°C. Samples were quickly thawed and sonicated for 1 min. DTT was added to a final concentration of 10 mM and incubated for 30 min with mixing at 800 rpm to reduce cysteines. Then 15 mM of freshly prepared iodoacetamide (IAA) was added and samples were incubated for 30 min at room temperature in the dark to alkylate cysteines. Afterwards, 1:100 (w/w) LysC (Wako Chemicals GmbH) was added for 4 h at 37°C with mixing at 800 rpm. Then urea concentration was diluted to 1.5 M with HPLC water and 1:50 (w/w) trypsin (Promega GmbH) was added for 12 h at 37°C with mixing at 700 rpm. Afterwards, the samples were acidified with 10% TFA and the cleavage of Rapigest was allowed to proceed for 30 min at 37°C. After spinning the sample for 5 min at 130,00x g at room temperature the supernatant was transferred to a new tube to proceed with peptide desalting. For desalting and cleaning-up of the digested sample, C-18 spin columns (Sep-Pak C18 Classic Cartridge; Waters) were used. A vacuum manifold was used for all washing and elution steps. First, the columns were equilibrated with 100% methanol and then washed twice with 5% (vol/vol) acetonitrile (ACN) and 0.1% (vol/vol) formic acid (FA). The sample was loaded two times and then the column was washed two times with 5% (vol/vol) ACN and 0.1% (vol/vol) FA. The undifferentiated cell samples were labeled using a “light” labeling reagent and the FACS-sorted neuronal cells were labeled using an “intermediate” labeling reagent inducing a mass shift of 28 or 32 D, respectively (Boersema et al, 2009). Formaldehyde and the D-isotopomer of formaldehyde react with primary amines of peptides (N-terminus and side chains of lysines) and generate a mass shift of 4 D. The labeling reagents consisted of 4.5 ml 50 mM sodium phosphate buffer (mixture of 100 mM NaH₂PO₄ and 100 mM Na₂HPO₄), pH 7.5, 250 μl 600 mM NaBH3CN and 250 μl 4% formaldehyde for light or 4% deuterated formaldehyde for intermediate labeling reagent, per sample. After the labeling procedure, the column was washed three times with 5% (vol/vol) ACN and 0.1% (vol/vol) FA. For elution 50% (vol/vol) ACN and 0.1% (vol/vol) FA was used. Labeled peptides from undifferentiated cells and FACS sorted neurons were pooled, dried in a vacuum concentrator, and resuspended in 20 mM ammonium formate (pH 10.0), to be ready for high pH reverse-phase peptide fractionation. To dissolve the dried samples, they were vortexed, mixed for 5 min at maximum speed in a thermomixer and sonicated for 90 s. The samples were stored at −20°C.

High pH reverse-phase peptide fractionation for dimethyl labeled samples

Offline high pH reverse-phase fractionation was performed using an Agilent 1200 Infinity HPLC System equipped with a quaternary pump, degasser, variable wavelength UV detector (set to 254 nm), peltier-cooled autosampler, and fraction collector (both set at 10°C). The column was a Gemini C18 column (3 μm, 110 Å, 100 × 1.0 mm; Phenomenex) with a Gemini C18, 4 × 2.0 mm SecurityGuard (Phenomenex) cartridge as a guard column. The solvent system consisted of 20 mM ammonium formate (pH 10.0) as mobile phase A and 100% acetonitrile as mobile phase B. The separation was accomplished at a mobile phase flow rate of 0.1 ml/min using the following linear gradient: 99% A for 2 min, from 99% A to 37.5% B in 61 min, to 85% B in a further 1 min, and held at 85% B for an additional 5 min, before returning to 99% A and re-equilibration for 18 min. Thirty two fractions were collected along with the LC separation that were subsequently pooled into 10 fractions. Pooled fractions were dried in a speed-vac and resuspended in 5% (vol/vol) ACN and 0.1% (vol/vol) FA and then stored at −80°C until LC–MS/MS analysis.

Sample preparation for in vitro differentiated mouse neurons

Frozen cell pellets of in vitro differentiated mouse neurons (∼1 million cells per sample) were thawed and resuspended in 100 μl of 1× PBS. An equivalent amount of 2× lysis buffer (200 mM Hepes pH 8.0, 100 mM DTT, 4% SDS) was added to the lysate, for a total volume of 200 μl. For neurons treated with Sec23a/b or control siRNA, cells (between 40,000 and 180,000 cells) were sorted directly into 2x lysis buffer. Samples were then sonicated in a Bioruptor Plus (Diagenode) for 10 cycles with 1 min ON and 30 s OFF with high intensity at 20°C. Samples were then boiled for 10 min at 95°C, and a second sonication cycle was performed as described above. The lysates were centrifuged at 184,07g for 1 min. Subsequently, samples were reduced using 10 mM DTT for 15 min at 45°C, and alkylated using freshly made 15 mM IAA for 30 min at room temperature in the dark. Subsequently, proteins were precipitated using acetone and digested using LysC (Wako sequencing grade) and trypsin (Promega sequencing grade), as described in Buczak et al (2020). The digested proteins were then acidified with 10% (vol/vol) trifluoroacetic acid. The eluates were dried down using a vacuum concentrator, and reconstituted samples in 5% (vol/vol) acetonitrile, 0.1% (vol/vol) formic acid. For data-independent acquisition (DIA)–based analysis (siRNA treated neurons), samples were transferred directly to an MS vial, diluted to a concentration of 1 μg/μl, and spiked with iRT kit peptides (Biognosys) before analysis by LC–MS/MS. For tandem mass tag (TMT)–based analysis (time course of in vitro differentiation), samples were further processed for TMT labeling as described below.

TMT labeling and high pH reverse-phase peptide fractionation

After desalting, peptides were dried in a vacuum concentrator and buffered using 0.1M Hepes buffer, pH 8.5 (1:1 ratio), for labeling, and then sonicated in a Bioruptor Plus for five cycles with 1 min ON and 30 s OFF with high intensity. 10–20 μg peptides were taken for each labeling reaction. TMT-10plex reagents (Thermo Fisher Scientific) labeling was performed by addition of 1 μl of the TMT reagent. After 30 min of incubation at room temperature with shaking at 600 rpm in a thermomixer (Eppendorf), a second portion of TMT reagent (1 μl) was added and incubated for another 30 min. After checking labeling efficiency, samples were pooled, desalted with Oasis HLB μElution Plate, and subjected to high pH fractionation before MS analysis. Offline high pH reverse-phase fractionation was performed using a Waters XBridge C18 column (3.5 μm, 100 × 1.0 mm; Waters) with a Gemini C18, 4 × 2.0 mm SecurityGuard (Phenomenex) cartridge as a guard column on an Agilent 1260 Infinity HPLC, as described in Buczak et al (2020). Forty-eight fractions were collected along with the LC separation, which were subsequently pooled into 16 fractions. Pooled fractions were dried in a vacuum concentrator and then stored at −80°C until LC–MS/MS analysis.

Data-dependent acquisition for dimethyl labeled samples (zebrafish neurons and stem cells)

The 10 fractions obtained by high pH fractionation were analyzed using a nanoAcquity UPLC system (Waters GmbH) connected online to a LTQ-Orbitrap Velos Pro instrument (Thermo Fisher Scientific GmbH). Peptides were separated on a BEH300 C18 (75 μm × 250 mm, 1.7 μm) nanoAcquity UPLC column (Waters GmbH) using a stepwise 145 min gradient between 3% and 85% (vol/vol) ACN in 0.1% (vol/vol) FA. Data acquisition was performed using a TOP-20 strategy where survey MS scans (m/z range 375–1,600) were acquired in the Orbitrap (R = 30,000 FWHM) and up to 20 of the most abundant ions per full scan were fragmented by collision-induced dissociation (normalized collision energy = 35, activation Q = 0.250) and analyzed in the LTQ. Ion target values were 1 × 10⁶ (or 500 ms maximum fill time) for full scans and 1 × 10⁵ (or 50 ms maximum fill time) for MS/MS scans. Charge states 1 and unknown were rejected. Dynamic exclusion was enabled with repeat count = 1, exclusion duration = 60 s, list size = 500 and mass window ±15 ppm.

Data-dependent acquisition for TMT-labeled samples (mouse in vitro differentiation)

The 16 fractions obtained by high-pH fractionation were resuspended in 10 μl reconstitution buffer (5% [vol/vol] acetonitrile, 0.1% [vol/vol] TFA in water) and 3 μl were injected. Peptides were separated using the nanoAcquity UPLC system (Waters) fitted with a trapping (nanoAcquity Symmetry C18, 5 μm, 180 μm × 20 mm) and an analytical column (nanoAcquity BEH C18, 2.5 μm, and 75 μm × 250 mm). The outlet of the analytical column was coupled directly to an Orbitrap Fusion Lumos (Thermo Fisher Scientific) using the Proxeon nanospray source. Solvent A was water, 0.1% (vol/vol) formic acid, and solvent B was acetonitrile, 0.1% (vol/vol) formic acid. The samples were loaded with a constant flow of solvent A at 5 μl/min, onto the trapping column. Trapping time was 6 min. Peptides were eluted via the analytical column at a constant flow of 0.3 μl/min, at 40°C reconstitution buffer (5% [vol/vol] acetonitrile, 0.1% [vol/vol] TFA in water), and 3.5 μl were injected. Peptides were eluted using a linear gradient from 5 to 7% in 10 min, then from 7% B to 28% B in a further 105 min and to 45% B by 120 min. The peptides were introduced into the mass spectrometer via a Pico-Tip Emitter 360 μm OD ×20 μm ID; 10 μm tip (New Objective), and a spray voltage of 2.2 kV was applied. The capillary temperature was set at 300°C. Full-scan MS spectra with mass range 375–1,500 m/z were acquired in profile mode in the Orbitrap with resolution of 60,000 FWHM using the quad isolation. The RF on the ion funnel was set to 40%. The filling time was set at a maximum of 100 ms with an AGC target of 4 × 105 ions and 1 microscan. The peptide monoisotopic precursor selection was enabled along with relaxed restrictions if too few precursors were found. The most intense ions (instrument operated for a 3 s cycle time) from the full scan MS were selected for MS2, using quadrupole isolation and a window of 1 D. HCD was performed with collision energy of 35%. A maximum fill time of 50 ms for each precursor ion was set. MS2 data were acquired with a fixed first mass of 120 m/z. The dynamic exclusion list was with a maximum retention period of 60 s and relative mass window of 10 ppm. For the MS3, the precursor selection window was set to the range 400–2,000 m/z, with an exclude width of 18 m/z (high) and 5 m/z (low). The most intense fragments from the MS2 experiment were co-isolated (using Synchronus Precursor Selection = 8) and fragmented using HCD (65%). MS3 spectra were acquired in the Orbitrap over the mass range 100–1,000 m/z and resolution set to 30,000 FWMH. The maximum injection time was set to 105 ms, and the instrument was set not to injections for all available parallelizable time.

Data-independent acquisition (Sec23a/b knockdowns)

Peptides were separated in trap/elute mode using the nanoAcquity MClass Ultra-High Performance Liquid Chromatography system (Waters; Waters Corporation) equipped with a trapping (nanoAcquity Symmetry C18, 5 μm, 180 μm × 20 mm) and an analytical column (nanoAcquity BEH C18, 1.7 μm, 75 μm × 250 mm). Solvent A was water and 0.1% formic acid, and solvent B was acetonitrile and 0.1% formic acid. 1 μl of the samples (∼1 μg on column) were loaded with a constant flow of solvent A at 5 μl/min onto the trapping column. Trapping time was 6 min. Peptides were eluted via the analytical column with a constant flow of 0.3 μl/min. During the elution, the percentage of solvent B increased in a nonlinear fashion from 0 to 40% in 120 min. Total run time was 145 min including equilibration and conditioning. The LC was coupled to an Orbitrap Exploris 480 (Thermo Fisher Scientific) using the Proxeon nanospray source. The peptides were introduced into the mass spectrometer via a Pico-Tip Emitter 360-μm outer diameter × 20-μm inner diameter, 10-μm tip (New Objective) heated at 300°C, and a spray voltage of 2.2 kV was applied. The capillary temperature was set at 300°C. The radio frequency ion funnel was set to 30%. For DIA data acquisition, full scan mass spectrometry (MS) spectra with mass range 350–1,650 m/z were acquired in profile mode in the Orbitrap with resolution of 120,000 FWHM. The default charge state was set to 3+. The filling time was set at a maximum of 60 ms with a limitation of 3 × 10⁶ ions. DIA scans were acquired with 40 mass window segments of differing widths across the MS1 mass range. Higher collisional dissociation fragmentation (stepped normalized collision energy; 25%, 27.5%, and 30%) was applied and MS/MS spectra were acquired with a resolution of 30,000 FWHM with a fixed first mass of 200 m/z after accumulation of 3 × 10⁶ ions or after filling time of 35 ms (whichever occurred first). Datas were acquired in profile mode. For data acquisition and processing of the raw data, Xcalibur 4.3 (Thermo Fisher Scientific) and Tune version 2.0 were used.

Data processing for dimethyl-labeled samples (zebrafish and rat neuronal differentiation)

Software MaxQuant (version 1.5.3.28) was used to search the MS.raw data. For D. rerio, the raw data were searched against the D. rerio UniProt database release: 2018_03, whereas for R. norvegicus, the .raw files from Frese et al (2017), were downloaded from PRIDE repository PXD005031 and searched against the UniProt R. norvegicus database release 2019_08. Both datasets were searched appending a list of common contaminants. The data were searched with the following modifications: Carbamidomethyl (C) (fixed) and Oxidation (M) and Acetyl (Protein N-term; variable). For D. rerio 2 labels, Light L (DmethLys0 and DmethNterm0) and Heavy H (DmethLys4 and DmethNterm4) were selected representing the stem cell and neurons, respectively. For the re-analysis of R. norvegicus data from Frese et al (2017), three different labels were used: Light L (DmethLys0 and DmethNterm0), Medium M, (DmethLys4 and DmethNterm4), and Heavy H (DmethLys8 and DmethNterm8). For identification, match between runs was selected with a match time window of 2 min, and an alignment time window of 20 min. The mass error tolerance for the full scan MS spectra was set at 20 ppm and for the MS/MS spectra at 0.5 D. A maximum of two missed cleavages was allowed. Identifications were filtered at 1% false discovery rate (FDR) at both peptide and protein levels using a target-decoy strategy (Elias & Gygi, 2007). From each experiment, iBAQ values (Schwanhäusser et al, 2011) and ratios between labels were extracted from the ProteinGroups.txt table. Differential expression analysis was performed using the mean of the normalized ratios between labels. The R package fdrtool (Strimmer, 2008) was used to calculate P-values and q values for the different comparisons, on the Log₂ transformed mean ratios.

Data processing for TMT10-plex data (mouse in vitro differentiation)

TMT-10plex data were processed using Proteome Discoverer v2.0 (Thermo Fisher Scientific). raw files were searched against the fasta database (UniProt M. musculus database, reviewed entry only, release 2016_11) using Mascot v2.5.1 (Matrix Science) with the following settings: Enzyme was set to trypsin, with up to one missed cleavage. MS1 mass tolerance was set to 10 ppm and MS2 to 0.5 D. Carbamidomethyl cysteine was set as a fixed modification, whereas oxidation of methionine and acetylation (N-term) were set as variable. Other modifications included the TMT-10plex modification from the quantification method used. The quantification method was set for reporter ions quantification with HCD and MS3 (mass tolerance, 20 ppm). FDR for peptide-spectrum matches (PSMs) was set to 0.01 using Percolator 13 (Brosch et al, 2009). Reporter ion intensity values for the PSMs were exported and processed with procedures written in R (v.4.0.5) and R studio server (v.1.2.5042 and 1.4.1106), as described in Heinze et al (2018). Briefly, PSMs mapping to reverse or contaminant hits, or having a Mascot score below 15, or having reporter ion intensities below 1 × 10³ in all the relevant TMT channels were discarded. TMT channels intensities from the retained PSMs were then log₂ transformed, normalized and summarized into protein group quantities by taking the median value using MSnbase (Gatto & Lilley, 2012). At least two unique peptides per protein were required for the identification and only those peptides with no missing values across all 10 channels were considered for quantification. Protein differential expression was evaluated using the limma package (Ritchie et al, 2015). Differences in protein abundances were statistically determined using the t test moderated by the empirical Bayes method. P-values were adjusted for multiple testing using the Benjamini-Hochberg method (FDR, denoted as “adj. P”) (Benjamini & Hochberg, 1995).

Data processing for DIA samples (Sec23a/b knockdowns)

DIA libraries were created by searching the DIA runs using Spectronaut Pulsar (v13) and Biognosys. The data were searched against species specific protein databases (UniProt M. musculus release 2016_01) with a list of common contaminants appended. The data were searched with the following modifications: carbamidomethyl (C) as fixed modification, and oxidation (M), acetyl (protein N-term). A maximum of two missed cleavages was allowed. The library search was set to 1% FDR at both protein and peptide levels. Libraries contained a total of 101,659 precursors, corresponding to 5,708 and 6,003 protein groups, respectively. DIA data were then uploaded and searched against this spectral library using Spectronaut Professional (v.14.10) and default settings. Relative quantification was performed in Spectronaut for each pairwise comparison using the replicate samples from each condition using default settings, except: data filtering set to Qvalue sparse, and imputation to RunWise. Differential abundance testing was performed using a paired t test between replicates. The data (candidate tables) and protein quantity data reports were then exported for further data analyses.

Data processing for human neuronal differentiation data

Protein identifications and LFQ intensity values (Log₂) in cultured iPSCs, NPCs and differentiated neurons, were obtained from the original Table S2 published in Djuric et al (2017). Differential expression analysis between the different conditions was performed on the log₂ LFQ intensity using the limma package (Ritchie et al, 2015).

Analysis of paralog pairs during development and across tissues

For the zebrafish development data (White et al, 2017), TPMs were used to calculate paralog pairs Pearson correlation coefficients. For the Human Tissue Atlas (Wang et al, 2019), Log₂(FPKM) and Log₂(IBAQ) were used to calculate correlation of paralog protein and transcript pairs. In all datasets, only genes and proteins identified in at least five time-points/tissues were considered for correlation analysis. Coefficient of variations (σ/mean protein or transcript expression along time points/tissues) was also calculated for every gene in each datasets. Genes that have at least one paralog in the genome according to Ensembl Compara were labeled as “Have Paralogs,” and used for further analysis. From all the possible paralog pairs, three categories were created. The first one indicates all the possible paralog gene pairs, the second one indicates paralog pairs residing in the same protein complexes according to definitions from Ori et al (2016), and the third one given by the exclusion between the two, indicating all other paralog pairs, namely paralog pairs that do not reside in the same complexes. For every paralog pair, the mean sequence identity was then calculated as the mean reciprocal identity retrieved from the Ensembl database. The relationship between sequence identity and co-expression between paralog pairs, was evaluated using Pearson R correlation coefficient, and visualized through a Generalized Additive Model.

Protein complex analysis during zebrafish development and across human tissues

For each datasets, proteins were annotated with the different protein complex definitions. Only protein complexes with at least five members present in each of the dataset were retained for analysis. For each of these complexes, all the possible pairwise correlations between complex members were considered, and from those the median value was used to calculate a median complex co-expression. We defined stable and variable complexes using the top and bottom 25% of the distribution, respectively. (1− median Pearson correlation) was also used to define then a measure of protein complex stoichiometric variability, as shown in Fig 1F and G. The distribution of correlations was then compared with a distribution of randomly assembled complexes of same size and complex members obtained by randomly assigning proteins/transcripts to complexes. For each dataset, the fraction of paralog pairs present was considered as the number of complex members that have at least one paralog in the genome divided by the total size of each protein complex. Finally for each subunit, we calculated expression correlation values with the other members of the same complex, taking the median of this value as a measure of co-expression for that specific subunit.

Paralog regulation during neuronal differentiation

For each datasets, differentially expressed proteins between different conditions (Log₂ Fold-Change > 0.58 and adjusted P-value, or fdr tools P-value < 0.05) were selected. Proteins were annotated as “Have Paralogs” if they had at least one paralog annotated in the genome. For each comparison, we then considered all possible paralog pairs present in the data and identified unique paralog pairs that displayed concerted regulation (same Log₂ Fold Change sign for both paralogs) or opposite regulation (different Log₂ Fold Change sign).

Complex members co-expression analysis for neuronal differentiation data

For calculating complex members stoichiometric variability, we adapted a previously established pipe-line (Gehring, 2021). For each condition and datasets, only protein complexes that had at least five quantified members were considered. Then for each subunit in each complex, the median euclidean distance of fold change between that subunit and all other complex members was calculated. The distance obtained was compared with a distribution of distances for 2,500 members from random complexes of equal size, obtained by randomly assigning proteins identified in the data to protein complexes. By comparing the two distributions we obtained a probability value for each subunit of observing lower distances with the complexes. Low P-values indicate high coexpression, denoted as stoichiometric stability, and vice versa.

eggNOG mapping

Fasta proteomes sequences used for MS protein quantification of the different dataset were annotated using emapper-2.1.4-2 (Cantalapiedra et al, 2021), based on eggNOG orthology data (Huerta-Cepas et al, 2019). Sequence searches were performed using the software MMseqs2 (Steinegger & Söding, 2017). For each proteome, eggNOG annotation was performed using default parameters, using the Vertebrate level mapping.

Conserved exchange of paralog proteins

For each dataset, protein quantification values were used to calculate paralog ratios across conditions. The log₂ paralog ratio between all possible quantified paralog pairs in each replicate was calculated for all the conditions tested. For each dataset, the significance of paralog ratio changes was assessed using the R package limma (Ritchie et al, 2015) considering replicates information. We considered only ratio changes relative to the first time point of each neuronal differentiation dataset. For comparison across species, each paralog pair was mapped to its relative eggNOG. Only paralog pairs where both entries could be mapped to a valid eggNOG were retained. After eggNOG mapping, shared eggNOG pairs between species were used to assess if specific paralog substitution were shared across different organisms, and for each specific comparison we combined the P-values using Fisher’s combined probability test from the metaRNASeq R package (https://cran.r-project.org/web/packages/metaRNASeq/index.html). Combined P-values were corrected for multiple testing using the Benjamin-Hochberg correction (Benjamini & Hochberg, 1995). Because in some cases multiple proteins can map to the same eggNOG, for each pair and condition the mean value was considered for both ratio differences and P-values. From this analysis, we considered as “conserved” only paralog gene pairs identified in all species and whose log₂ ratio changes were consistent in sign in at least five of the seven neuronal differentiation comparisons, with combined adjusted P ≤ 0.05.

Transcription factor enrichment analysis

Mouse promoter region sequences were obtained via the R package biomaRt (Durinck et al, 2009), selecting a 2,000-bp region upstream of the transcription starting site. TFs list were obtained from the JASPAR database (Fornes et al, 2020) via the TFBSTool BioconductoR package (Tan & Lenhard, 2016). Selecting a minimal score for enrichment of 80%. TFs that were not identified in a specific promoter were given an enrichment score of 0.

GO enrichment analysis

Overrepresentation analysis of GO terms was performed with the R package topGO (Alexa & Rahnenfuhrer, 2021). Fisher test was used to estimate the expected proportion for different terms and obtain a P-value indicating the enrichment score for each specific GO term. Gene set enrichment analysis was performed with the topGO R package using a Kolmogorov–Smirnov test on the cumulative ranked distributions. For both enrichments P-values were adjusted using Hommel’s correction, GO terms were considered significant if their adjusted P-values were below the value of 0.05. The R package rrvgo (https://ssayols.github.io/rrvgo/) was used to summarize and reduce redundancy of the enriched GO terms using default settings.

Figure generation

Data visualization was performed with R (v.4.0.5) and R studio server (Version 1.4.1106) using the ggplot2 package (Wickham, 2009). Figs 1A, 2A, and 4A and graphical abstract were created with https://BioRender.com.