A maturational shift in the frontal cortex synaptic transcriptional landscape underlies schizophrenia-relevant behavioural traits: A congenital rat model

Disruption of brain development early in life may underlie the neurobiology behind schizophrenia. We have reported more immature synaptic spines in the frontal cortex (FC) of adult Roman High-Avoidance (RHA-I) rats, a behavioural model displaying schizophrenia-like traits. Here, we performed a whole transcriptome analysis in the FC of 4 months old male RHA-I (n = 8) and its counterpart, the Roman Low-Avoidance (RLA-I) (n = 8). We identiﬁed 203 signiﬁcant genes with overrepresentation of genes involved in synaptic function. Next, we performed a gene set enrichment analysis (GSEA) for genes co-expressed during neurodevelopment. Gene networks were obtained by weighted gene co-expression network analysis (WGCNA) of a transcriptomic dataset containing human FC during lifespan (n = 269). Out of thirty-one functional gene networks, six were signiﬁcantly enriched in the RHA-I. These were differentially regulated during infancy and enriched in biological ontologies related to myelination, synaptic function, and immune response. We validated differential gene expression in a new cohort of adolescent ( < = 2 months old) and young-adult ( > = 3 months old) RHA-I and RLA-I rats. The results conﬁrmed overexpression of Gsn, Nt5cd1, Ppp1r1b, and Slc9a3r1 in young-adult RHA-I, while Cables1 , a regulator of Cdk5 phosphorylation in actin regulation and involved in synaptic plasticity and maturation, was signiﬁcantly downregulated in adolescent RHA-I. This age-related expression change was also observed for presynaptic components Snap25 and Snap29. Our results show a different maturational expression proﬁle of synaptic components in the RHA-I strain, supporting a shift in FC maturation underlying schizophrenia-like behavioural traits and adding construct validity to this strain as a neurodevelopmental model


Introduction
Schizophrenia is a complex disorder where geneenvironmental interactions pose an imprint on early neurodevelopment, priming susceptibility to risk factors later in life ( Owen et al., 2016 ). The genomic architecture behind this vulnerability is large, with this polygenicity also being highly pleiotropic as many of the common risk variants and rare copy number variations also relate to other disorders. These include autism spectrum disorders (ASD) and attention-deficit/hyperactivity disorder (ADHD) ( Gudmundsson et al., 2019 ;Reay and Cairns, 2020 ;Hall et al., 2021 ). Relevant pleiotropic loci are mostly located within genes that play important roles in neurodevelopmental processes ( Cross-Disorder Group of the Psychiatric Genomics, 2019 ) supporting that schizophrenia is a neurodevelopmental disorder related to ASD and ADHD ( Gudmundsson et al., 2019 ). In view of this, focusing on symptom phenotypes and their neural substrates will provide a better understanding of the aetiology behind schizophrenia.
Most animal models of schizophrenia have been developed to mimic specific symptoms resembling primarily positive (psychotic) symptoms ( Jones et al., 2011 ). However, when dealing with complex multifaceted, multicausal, and polygenic symptoms and traits such as in schizophrenia, animal models with a more varied presentation of behavioural traits, analogous to negative and cognitive symptoms, offer a stronger translational validity . One of these models is the inbred Roman High-Avoidance (RHA-I) rat strain. Contrary to bottom-up approaches where a given pathophysiological mechanism is assumed and linked to a specific behavioural symptom by gene-targeting approaches, the RHA-I strain has been selected according to a top-down approach where a given behavioural phenotype is recreated and its neurogenetic underpinnings investigated ( Giorgi et al., 2019 ;Fernandez-Teruel et al., 2021 ). Thus, impaired executive function has been reported repeatedly in the RHA-I strain when compared to its counterpart, the inbred Roman Low-Avoidance (RLA-I) rat strain and outbred rat stocks as external controls ( Lopez-Aumatell et al., 2009 ;Moreno et al., 2010 ;Diaz-Moran et al., 2012 ;Martinez-Membrives et al., 2015 ;Oliveras et al., 2015 ;Esnal et al., 2016 ). Underlying this behavioural phenotype, the RHA-I strain displays alterations in dopamine, serotonin, and glutamate neurotransmitter systems ( Giorgi et al., 2003 ;Klein et al., 2014 ;Oliveras et al., 2017 ;Fomsgaard et al., 2018 ), grey and white matter structures ( Rio-Alamos et al., 2017, 2018, frontal cortex (FC) activation ( Tapias-Espinosa et al., 2019 ) as well as increased susceptibility to neonatal handling and social isolation ( Rio-Alamos et al., 2017, 2018Sánchez-González et al., 2020 ). The RHA-I is thereby a model with strong face and construct validity that can confer us with important new knowledge regarding the neurobiology behind schizophrenia-like traits.
In this regard, in the FC we have recently reported differential expression of synaptic markers during neurodevelopment as well as increased thin-spine density in pyramidal neurons suggestive of a more immature cortical endophenotype in the RHA-I strain ( Elfving et al., 2019 ;Sanchez-Gonzalez et al., 2021 ). This is highly in-teresting since FC maturation is crucial for development of executive functions and cognitive skills ( Catts et al., 2013 ;Silbereis et al., 2016 ), functions that are impaired in schizophrenia ( Weickert et al., 2000 ). A failure to reach the final stage of cortical maturation in individuals at risk, resulting in retainment of a more immature cortex, has been proposed to be at the core of schizophrenia ( Catts et al., 2013 ;Gao et al., 2022 ). If true, there may be an opportunity for preventive interventional approaches in individuals genetically at risk. In this respect, it is important to determine the critical time-windows of neurodevelopment. By definition, adolescence is considered the critical period for schizophrenia onset as this is where the first symptoms usually are manifested ( Patel et al., 2021 ). According to the dual-hit hypothesis, the disruption in the neurodevelopmental trajectory, starting from early embryonic development and up to childhood, would confer the brain with increased vulnerability leading to the disease manifestation later in life ( Vasistha et al., 2020 ;Guerrin et al., 2021 ;Malwade et al., 2022 ). Furthermore, a recent single cell analysis of the prefrontal cortex in schizophrenia patient tissue revealed that enhanced synaptic plasticity in upper cortical layers in schizophrenia is derived from dysregulation of developmental transcription factors ( Batiuk et al., 2022 ).
To address developmental origin of schizophrenia in the Roman rat strains, we mapped the FC transcriptome in adult Roman rats and investigated whether there is a shift in the RHA-I compared to the RLA-I strain in functional gene networks differentially expressed during early neurodevelopment. To identify co-expressed gene networks associated with neurodevelopment, we analysed the BrainCloud © dataset that contains microarray-based dorsolateral prefrontal cortex transcriptomes of 269 individuals throughout the human lifespan ranging from foetal to 80 years old ( Colantuoni et al., 2011 ). Previous studies have revealed neurodevelopmental processes to be associated with specific co-expressed gene modules ( Kang et al., 2011 ;Li et al., 2018 ). The strength of applying co-expressed gene modules as functional units is that it allows us to compare them against specific internal and external traits ( Langfelder and Horvath, 2008 ). Here, we wanted to compare them to time windows of development and behavioural traits. We validated genes of interest across brain maturation in a new cohort of adolescent and young-adult RHA-I and RLA-I rats.

Animals
Roman male rats from the breeding colony at the Dept. Psychiatry and Forensic Medicine, Universitat Autònoma de Barcelona, were used. All animals were housed in pairs in macrolon cages (50 × 25 × 14) and kept with food and water ad libitum, maintained under a 12:12h light-dark cycle (lights on at 08:00 a.m.) and with controlled temperature (22 ± 2 °C) and humidity (50-70%). Animals were euthanised in accordance with the Spanish Royal Decree (RD 53/2013) for the protection of experimental animals and with the European Communities Council Directive (2010/63/EU). For the transcriptomic analysis, naïve adult male inbred RHA-I (n = 8) and RLA-I (n = 8) rats were sacrificed at approximately P120. A second cohort of naïve male inbred Roman rats for real-time quantitative Polymerase Chain Reaction (qPCR) were sacrificed either at P44-46 (95 -165 g) (considered adolescent RHA-I (n = 7) and RLA-I (n = 5)) or at P90 (considered young-adult RHA-I (n = 8) and RLA-I (n = 8)). Rats' weight range at P90-P120 was 320-400 g. The brains were removed, and their FC dissected by free hand. The tissue was transferred to RNase-free tubes, frozen with liquid nitrogen and then stored at -80 °C until further analysis.

RNA extraction and sample preparation
RNA extraction was performed as described in . Briefly, RNA was extracted with miRNeasy mini kit (Qiagen; cat. no. 217004) for the transcriptomics cohort. For the adolescent and young-adult cohort, RNA was extracted with RNeasy Lipid Tissue Mini Kit (Qiagen; cat. no. 74804) automated in the QIAcube (Qiagen). Subsequently, the samples were DNAse treated with a Turbo DNA-free kit (Ambion; cat. no. AM1907). RNA integrity (RIN) values and concentrations were determined using the 2100 Agilent Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA). One sample was excluded in the first cohort due to low RIN ( < 5), ending with a total of 7 RHA-I and 8 RLA-I samples for the transcriptomics analysis. RNA was kept at -80 °C until further analysis.

Library preparation
Following the manufacturer instructions of DNA Library Preparation Kit (Illumina Nextera XT; cat. no. FC-131-1024), tagmentation mix was added to each cDNA sample ( ∼0.5 ng) up to a 10 μl volume, and the tagmentation reaction run at 55 °C, followed by neutralisation with NT buffer for 5 min at room temperature. Nextera Polymerase Chain Reaction (PCR) master mix (NPM), index primer i5, Fig. 1 Analysis workflow for the study. A Weight Gene Co-expression Network Analysis (WGCNA) was performed on an already existing transcriptomic dataset with human FC gene expression throughout lifespan, the BrainCloud®. Co-expressed gene modules were transformed into gene sets for performing a gene set enrichment analysis (GSEA) on a transcriptomic dataset from adult RHA-I and RLA-I cohort. Gene sets significantly enriched were validated in a new cohort of adolescent and young-adult RHA-I and RLA-I rats, by selecting representative genes and investigating gene expression across the different age groups. Pathway analysis was performed on co-expressed gene modules enriched in the WGCNA analysis and on the list of differentially expressed genes from the Roman rat cohort. Co-expressed gene modules were also correlated to the different time windows of neurodevelopment. Figure created by BioRender. and index primer i7 from the Index Kit (Illumina Nextera XT; cat. no. FC-131-1001) was used for preparing the sequencing libraries. Used indices: N701-N706, S502-S504, and S517. PCR was performed using the conditions described in ( Picelli et al., 2014 ), and so was purification of the PCR product (libraries). Libraries were stored at 4 °C overnight. All libraries were individually normalised and diluted to a final concentration of 2 nM and pooled for sequencing.

RNA sequencing
The RNA sequencing was performed on an Illumina NextSeq 500 in a single-end, unstranded format (Illumina, FC-404-2005) according to manufacturer's instructions, and the libraries were sequenced at a depth of approximately 26 million reads. For each sample, four technical replicates fastq files were obtained.

The BrainCloud® dataset
We downloaded the BrainCloud® dataset that includes 269 human samples, ranging from foetal development (negative ages) through ageing (80 years), and contains microarray mRNA data from dorsolateral prefrontal cortex Brodmann areas (DLPFC BA) 46/9 from the GEO repository (GSE30272) together with the annotation file and the demographic data. mRNA levels are normalised to a reference RNA, that comprises a pool of all samples. Subjects had no severe neuropathology nor neurological or neuropsychiatric diagnoses ( Colantuoni et al., 2011 ). In the microarray design there are multiple probes which are targeting the same gene. In this case the probes were aggregated to their mean value and as a result 17,160 genes left for the analysis. We created an additional trait for each sample, corresponding to their time-window (TW) of neurodevelopment defined as: TW1 (ranging from 0-11 months), including 16 samples; TW2 (ranging from 1-16 years) including 47 samples; TW3 (ranging from 17-29 years), including 49 samples; and TW4 (ranging from 30-78 years), including 115 samples.

Weighted gene co-expression network analysis (WGCNA)
WGCNA was performed on the BrainCloud® dataset using the WGCNA R-package version 1.70-3 ( Langfelder and Horvath, 2008 ). Following the WGCNA tutorial, we first check for outliers and unwanted variability in the BrainCloud dataset associated with batch effect, demographic traits, and tissue confounders/brain quality markers such as RIN, pH, and postmortem interval. This was done

Fig. 2
RNAseq data analysis for the RHA-I and RLA-I rats. In (A) is illustrated the volcano plot with the differentially expressed (DE) genes. In blue, genes downregulated, and in red, genes upregulated in the RHA-I compared to the RLA-I. In (B) are indicated the significant GO pathways and the significant KEGG pathways for the DE genes. by creating a sample network based on squared Euclidean distance (adjacency-function) and converting the sample traits into a colour representation. Samples with a standardised connectivity value below a threshold of -2.5 (default) were removed along with samples clustering together based on traits (unwanted variance). The final data frame used for the WGCNA contained 227 individuals. In order to satisfy the scale free topology criterion recommended by ( Langfelder and Horvath, 2008 ), we set the soft thresholding power at 4. Next, a one-step network construction and module detection were performed, considering a signed hybrid network, using biweight midcorrelation analysis ( Fig S1A). Co-expressed gene modules were identified, and for each module an eigengene was automatically defined. The module eigengene (ME) is representative of the gene expression profiles that best characterise the gene correlations within that module. Each gene was assigned to a module (or to a group of genes without enough relationship to any module), and received a module membership (MM) value, which indicates the intramodular connectivity of the gene to its specific module. For each gene module, a gene list was created by extracting the genes from the dataset with a MM > 0.6. The gene lists were used for the downstream GSEA and pathway analysis.

Gene set enrichment analysis (GSEA)
We used the R-package biomaRt (version 2.48.3) ( Durinck et al., 2009 ) for converting the gene lists generated from the BrainCloud® dataset into rat homologs gene sets. GSEA was performed using the R-package fgsea (version 1.18.0). Genes from the RHA-I/RLA-I sequencing data were ranked based on gene stats (p-value) in RHA-I against RLA-I in decreasing order, and GSEA run with adjusted p (adj.p) < 0.05, a minimum gene set size set to 5, maximum gene set size set to 600, and n-permutation set to 1000. Leading edge analysis was included in the GSEA analysis.

Real-Time qPCR analysis
Before complementary DNA (cDNA) synthesis of the new cohort of adolescent and young-adult RHA-I/RLA-I, the RNA concentration of the samples was adjusted to match the sample with the lowest concentration. RNA was reversely transcribed using random primers and Superscript IV Reverse Transcriptase (Sigma-Aldrich, St. Louis, MO, USA) following the manufacturer's instructions with an RNA  Table S4. Each real-time qPCR reaction (10 μl total volume) contained 5 μl SYBR Green (Sigma-Aldrich, St. Louis, MO, USA), 0.5 μl of 10 μM primer pair mix, 1.5 μl DEPC water, and 3 μl diluted cDNA. The thermal conditions for the PCR were 3 min at 95 °C to activate the hot-start iTaqDNA polymerase, followed by 40 cycles of 10 s denaturation at 95 °C, 30 s annealing at 60 °C, and 60 s extension at 72 °C. Each run was completed by dissociation curve analysis to confirm the amplification specificity and absence of primer dimers, 1 min 95 °C, 30 s 60 °C, and 30 s 95 °C. Primers were only included if the regression coefficient was above 0.9 and the amplification efficacy between 90% and 110%. Stability comparison of the expression of the reference genes was conducted with the NormFinder software ( http://moma.dk/normfinder-software ). The most stable combination of the reference gene expression was Actb/Hprt1. The value from each individual sample was normalised with the geometric mean of the relevant reference gene combination.

Bioinformatics and statistical analysis
The snakemake workflow management system ( Köster et al., 2021 ) was used to combine tools for the analysis of preprocessing the reads as well as for quantification and normalisation. In order to Modules are represented by the colours assigned through the co-expression network construction and module detection. In red, the modules' hubs, which is the most connected gene in the module. (C) Module Eigengenes (MEs) for each significant module plotted for the different TWs. All six modules were significantly either up-or down-regulated in TW1, corresponding to infancy, equivalent to the direction of enrichment observed in the RHA-I in (A). Mean ±SEM, * * * * p < 0.0001, * * p < 0.01. inspect the quality of the raw data, several statistics were assessed with FastQC which assured excellent quality of the data ( http: //www.bioinformatics.babraham.ac.uk/projects/fastqc/ ). To ensure that poor quality data are not included in the analysis we paid extra attention to the sequence (per base and per sequence) quality scores and the over-represented sequences which can indicate potential adapter contamination. All the data were aligned to the reference rat genome sequence (Rattus norvegicus v.Rnor_6.0.101) that was retrieved from the Ensembl website ( http://www.ensembl.org ). The alignment was performed using the STAR software (v. 2.7.2) ( Dobin et al., 2013 ) and the data of the same samples that were sequenced on different lanes were concatenated in one file. The quality of the alignment was assessed with Alfred (v. 0.1.17) ( Rausch et al., 2019 ). Next, the gene expression was quantified using the htseq-count ( Anders et al., 2015 ) by counting the number of reads overlapping the ensemble gene annotations and resulted in a count table with the number of reads aligned in the gene features. We sequenced between 5.6 and 8.9 million reads per sample (avg 7.33 million reads) and reads mapping was between 79.1% and 85.2% uniquely mapped reads of the total reads aligned.
The DESeq2 R-package (version 1.26.0) was used for differential expression (DE) analysis. The data was pre-filtered to keep rows that had at least one read count. Alpha was set to 0.05 for FDR p-value adjustment. Ensembl IDs were converted to gene symbol, gene name and Entrez ID with help of the genome wide annotation R-package org.Rn.eg.db (version 3.13.01). The R-package cluster-Profiler (version 4.0.4) was used for gene ontology (GO) (enrichGO function) and Kyoto Encyclopaedia Genes and Genomes (KEGG) (en-richKEGG function) pathway analysis of the significant genes. For identifying significantly enriched GO terms, we set an adj.p < 0.05 (Benjamini-Hochberg (BH)) and a relaxed false discovery rate (FDR) q-value < 0.2. For identifying KEGG terms significantly enriched, we set a p-value < 0.01 (BH) and an FDR q-value < 0.2. For comparing gene modules across TWs, the mean ME of the samples grouped in each TW was calculated. We tested for normality using Shapiro-Wilk test. If the data were not normally distributed, Kruskal-Wallis and Dunn tests were used. If normally distributed, data were tested for homogeneity of variance using Levene test. If the data did not fulfill the criteria for homogeneity of variance, Welch and Games-Howell tests were used, but if the data fulfilled the criteria, ANOVA and Tukey tests with p-value < 0.05 were used. For the RT-qPCR data analysis, mixed-effects model with the Geisser-Greenhouse correction and genes of interest as matched values, followed by BH with adj.p < 0.05, was applied for comparing gene expression between RHA-I and RLA-I in both adolescent and young-adult cohort. Grubb's test was run before analysis for outlier detection. Graphs and statistical analysis were done in R and GraphPad prism. Table 2 Co-expressed gene modules extracted from the WCNA analysis on the BrainCloud® dataset. Each module has been assigned a biological function based on GO analysis.

Gene modules co-expressed during human FC development are related to synaptic function
Upon analysis of WGCNA, 37 co-expressed gene modules were identified ( Figure S1A) containing between 22 and 1,721 genes (Table S1). Each module was automatically assigned a colour with grey containing genes not assigned to a functional module (8,377 genes (49%)). A heatmap was generated showing ME/traits correlations (Fig. S1B). For all modules, no correlation between MEs and array batch, sex, race, post-mortem interval, and pH were seen. A moderate correlation ( < 0.53) was seen for some ME in relation to brain bank source (BBS). As expected, MEs correlated with age. When correlating MEs with the trait TW, the infancy time-window (TW1), and childhood and early adolescence (TW2) showed the highest correlations (Fig. S1C).
For each gene module, a cut-off at MM > 0.6 was set, and a gene list was generated for performing pathway analysis. In total, 31 modules were significantly enriched for biological processes GO terms and/or KEGG pathways. Based on this we assigned each module with a biological function ( Table 2 ;  full list Table S2). Sienna3 module was not enriched for GO terms, and therefore assumed to be an artefact and omitted from further analyses.

Gene sets enriched in the RHA-I are differentially regulated during infancy
The gene lists extracted from the co-expressed gene modules were translated to rat gene homologs and used for running a GSEA on the RHA-I and RLA-I full transcriptome data set (Table S3). Six of the gene sets were significantly enriched (adj.p < 0.05) ( Table S3). The gene sets corresponded to module Lightyellow (mitochondrial respiration), Darkgrey (neuroinflammation), Brown (myelination), Lightcyan (ion sequestration), Blue (synapse formation) and Green (synapse function) in the BrainCloud® dataset ( Fig. 3 B). The module Steelblue was also significant but due to the low number of genes ( < 15) we decided not to include this gene set in the further analysis to avoid a false positive finding. Hub genes for each of the significant modules are illustrated in Fig. 3 B. When plotting mean MEs for the gene modules according to their TWs, we observed that the Lightyellow (mitochondrial respiration), Darkgrey (neuroinflammation), Brown (myelination) and Lightcyan (ion sequestration) modules are significantly downregulated and Blue (synapse formation) and Green (synapse function) upregulated in the BrainCloud dataset for TW1, corresponding to infancy ( Fig.3 C, Table 3 ).

Cables1 is downregulated in young RHA-I
Leading-edge analysis was performed after the GSEA analysis (Table S3). From the differentially expressed gene sets, we selected six genes from top 10 genes for further realtime qPCR validation in the new cohort of adolescent and young-adult RHA-I/RLA-I rats for validating gene expression changes during specific stages of development. These were Cables1, Gsn, Nt5dc1, Slc9a3r1, and Il12rb2. Ppp1r1b was also included for its specific interest in schizophrenia. In the young-adult RHA-I and RLA-I rats, there was a significant effect of strain (F(1,13) = 9.675, p < 0.01) on Gsn (p < 0.05, q < 0.05), Nt5cd1 (p < 0.01, q < 0.01), Slc9a3r1 (p < 0.0001, q < 0.001), and Ppp1r1b (p < 0.05, q < 0.05) gene levels showing a significant upregulation in the RHA-I rats compared to the RLA-I rats ( Fig 4 A). In the adolescent group, we identified a significant gene x strain interacion (F(5,68) = 9.360; p < 0.0001) for Cables1 , that was downregulated in the RHA-I rats compared to the RLA-I rats (p < 0.001, q < 0.05) ( Fig. 4 A).

Presynaptic components are divergently expressed in the FC of adolescent and young-adult RHA-I rats
Based on this observation of a divergent expression of Ca-bles1 in the adolescent RHA-I strain, we decided to add to the analysis a subset of synaptic-related gene targets that we previously reported are differentially regulated in the FC of adult RHA-I ( Elfving et al., 2019 ). We grouped them in the following categories: neurotransmitter receptors, postand presynaptic components, and neurotrophic signaling ( Fig. 4 B). In the receptor group, there was a significant interaction effect for gene, strain, and age (F(6,53) = 5.067, p < 0.01) for Grm2 which showed lower expression in the RHA-I strain, both in the adolescent (p < 0.01, q < 0.01) and young-adult group (p < 0.001, q < 0.001) when compared to the RLA-I strain. Grin2b and Drd1 on the contrary, were both nominally significantly upregulated in the RHA-I strain when compared to the young-adult RLA-I (p < 0.05, q = 0.2). For the postsynaptic targets, there was also this general increased gene expression in the young-adult RHA-I. We observed a significant interaction effect for gene, strain, and age (F(9,80) = 2.425, p < 0.05) for Nrg1 , which was significantly increased in young-adult RHA-I compared to youngadult RLA-I (p < 0.0001, q < 0.0001), but also between adolescent RHA-I and young-adult RHA-I (p < 0.05, q < 0.05). Both Homer1 (p < 0.05, q = 0.2) and Homer3 (p < 0.05, q = 0.1) were nominally significantly increased in young-adult RHA-I when compared to respectively adolescent RHA-I or young-adult RLA-I. For the presynaptic targets, there was a general increase in gene expression in the young-adult RHA-I and a decrease in the adolescent RHA-I. For Snap25 expression, there was a significant interaction effect for strain and age (F(3,27) = 8.807, p < 0.001), with the adolescent RHA-I having a lower expression (p < 0.05, q < 0.05) and the youngadult RHA-I a higher expression (p < 0.01, q < 0.01) compared to their counterparts. The same pattern was observed for Mean ±SEM. * p < 0.05, * * p < 0.01, * * * p < 0.001, * * * * p < 0.0001, with adj.p; #p < 0.05, without adj.p.
Snap29 , another SNARE component, and for synaptophysin, Syp. Lastly, when looking at neurotrophic signalling, here too there was a general increased gene expression in the young-adult RHA-I with a significant gene x strain x age effect (F(86,54) = 3.284, p < 0.01), with higher Bdnf expression in the RHA-I strain compared to the RLA-I during adulthood (p < 0.001, q < 0.05) while its receptor, Trkb , showing a steeper increase in expression from adolescent to adulthood in the RHA-I strain (p < 0.01, q < 0.02).

Discussion
This is the first study performing a whole transcriptome analysis of the FC in the RHA-I and RLA-I rat strains. Among the DEGs were genes previously reported to be differentially regulated in the FC of RHA-I rats. These genes included G rm2 ( Fomsgaard et al., 2018 ;Elfving et al., 2019 ), Vamp2 ( Elfving et al., 2019 ) , Syp ( Elfving et al., 2019 ), Pvalb ( Tapias-Espinosa et al., 2022 ), Drd1 ( Elfving et al., 2019 ) and Ppp1r1b ( Guitart-Masip et al., 2006 ). When performing a pathway analysis for the DEGs there was a sig-nificant GO enrichment in synaptic components, supporting our previous observations of synaptic alterations in the RHA-I strain ( Elfving et al., 2019 ;Sanchez-Gonzalez et al., 2021 ). We also identified enrichment of KEGG pathways related to cocaine and alcohol addiction, fitting well with the RHA-I rats being more prone to substance abuse and relapse ( Giorgi et al., 2007( Giorgi et al., , 2019Fattore et al., 2009 ). The Roman rats used for the transcriptomics analysis were 4 months old corresponding to adulthood in human years ( Semple et al., 2013 ;Sengupta, 2013 ). The FC is one of the last brain regions to mature, with appropriate development being crucial for high-order cognitive abilities ( Kroeze et al., 2018 ;Chini and Hanganu-Opatz, 2021 ). Therefore, a shift in maturation is proposed to be at the core of schizophrenia ( Catts et al., 2013 ;Skene et al., 2017 ). Here, we show that in adult RHA-I rats there is a discrepancy, when compared to their counterpart the RLA-I rat strain, in the transcriptomic landscape of gene networks co-expressed during infancy.
This suggests an alteration in the RHA-I FC gene expression profile of genes highly relevant to brain maturation. Infancy is the period of neurodevelopment where most expression changes occur and present the highest number of DEGs when compared to other periods of neurodevelopment ( Kroeze et al., 2018 ). Correspondingly, we observed in the BrainCloud © WGCNA analysis a high number of the coexpressed gene modules correlating with the trait infancy. The co-expressed modules we obtained in our study overlapped with the modules reported by another study applying WGCNA for mapping expression trajectories of human FC neurodevelopment ( Radulescu et al., 2020 ). It has previously been reported that a high degree of co-regulated transcriptional networks differentially regulated during late foetal and infancy periods are related to neural cell proliferation and migration, dendrite and synapse development, and myelination ( Kang et al., 2011 ). What was interesting is that many of these gene networks contained singlenucleotide polymorphisms (SNPs) related to schizophrenia risk ( Kang et al., 2011 ). Indeed, foetal and early infancy are the periods where genes harbouring risk for schizophrenia are most dynamically expressed ( Gulsuner et al., 2013 ;Clifton et al., 2019 ). In our study, the modules that were significantly enriched in the adult RHA-I compared to the RLA-I were related to synapse regulation, but also to mitochondria function, immune regulation, myelination, and ion sequestration, which are all processes that have been associated with brain development and schizophrenia ( Maschietto et al., 2015 ;Schwede et al., 2018 ;Park et al., 2020 ;Fišar, 2023 ). More specifically, one of the hub genes, NTSR2, that codes for the neurotensin receptor 2, has been linked to sensorimotor gating and locomotor activity ( Feifel et al., 2010 ;Oliveros et al., 2010 ), which are functions repeatedly shown to be disrupted in the RHA-I strain ( Oliveras et al., 2015. Also of interest, recent single-cell analysis of schizophrenia prefrontal cortex has revealed an upregulation of synaptic pathways with concomitant downregulation of metabolic pathways, and further analysis of synaptic pathways demonstrated developmental origin of dysregulated transcriptional factor network ( Batiuk et al., 2022 ).
When looking at gene expression of selected genes, we found significant upregulation of the gelsolin gene ( Gsn ) in young-adult RHA-I compared to RLA-I rats. This gene is involved in dendritic spine stabilisation through directly regulating filamentous actin assembly and disassembly ( Pontrello and Ethell, 2009 ). The increased gene expression in the young-adult RHA-I rats could be a compensatory effect to the higher number of immature dendritic spines in this strain . Other genes upregulated in the young-adult RHA-I strain were Nt5dc1, which coding variants have been associated with ADHD ( Zayats et al., 2016 ;Gudmundsson et al., 2019 ) and Slc9a3r1, also known as NHERF1 , which is a PDZ scaffold that interacts with the mGlu2/3 receptor ( Ritter-Makinson et al., 2017 ). Overexpression of Slc9a3r1 could well be a consequence of the decreased Grm2 expression in the RHA-I strain ( Klein et al., 2014 ;Fomsgaard et al., 2018 ) and an attempt to regulate glutamate homeostasis, as Slc9a3r1 also is a scaffold protein for the astrocytic glutamate transporter ( Lee et al., 2007 ). Ppp1r1b , that was also upregulated in the young-adult RHA-I rats, codes for the cAMP-regulated phosphoprotein of molecular weight 32 kDa (DARPP-32), a protein involved in the regulation of dopaminergic signalling in striatum and FC, and linked to schizophrenia ( Kunii et al., 2019 ). Ppp1r1b expression is associated with positive or negative decision outcomes ( Frank et al., 2009 ), suggesting an association between the strain effect on Ppp1r1b expression and the different outcomes of reinforcement-based and conflict-related instrumental learning tasks between the RHA-I and the RLA-I strains ( Cuenya et al., 2012 ;Corda et al., 2018 ;Fernández-Teruel et al., 2021 ;Bellés et al., 2023 ). As DARP-32 colocalises with D1R and D2R in cortex ( Rajput et al., 2009 ), the increased Ppp1r1b expression may be linked to the increased expression for Drd1 in the young-adult RHA-I strain.
We corroborated earlier observations of synaptic changes in the adult RHA-I ( Elfving et al., 2019 ). Moreover, we confirmed that the decreased Grm2 expression in the RHA-I strain due to the presence of the stop codon mutation ( Fomsgaard et al., 2018 ) is already manifested during the adolescent period of neurodevelopment. This is relevant, as it is hypothesised that dysfunction in the glutamatergic input during cortical development may be a convergence point for the aberrant maturation and network dysfunction associated with schizophrenia Gao et al., 2022 ).
Surprisingly, we observed a decreased expression of Cables 1 and the presynaptic components Snap25 and Snap29 in the adolescent RHA-I rats when compared to the adolescent RLA-I. Cables1 is an activator of the cyclin-dependent kinase 5 (Cdk5) ( Zukerberg et al., 2000 ), which is an important modulator of synaptic plasticity through regulation of NRB2 degradation ( Hawasli et al., 2007 ) as well as a stabiliser of dendritic spines into maturational states ( Huang et al., 2017 ). Lower expression of Cables1 may therefore indicate impaired synaptic maturation in the adolescent RHA-I strain. Cdk5 is also involved in recruitment and assembly of presynaptic proteins during synaptogenesis ( Easley-Neal et al., 2013 ), explaining the lower expression of Snap25 and Snap29 in the adolescent RHA-I strain.
Interestingly, schizophrenia has been associated with dysregulation of Cdk5 ( Engmann et al., 2011 ). Our results suggest that an aberrant synaptic maturation during the adolescent period leads to permanent neurodevelopmental anomalies in the RHA-I rats as they retain a more immature synaptic phenotype. This would translate into the cognitive dysfunctions characterising the adult RHA-I strain ( Fernández-Teruel et al., 2021 ). Indeed, in a recent study, a deviated maturational trajectory for sensorimotor gating and startle habituation was reported for the RHA-I rat strain . Rats were tested during prepubescence, adolescence, and adulthood, and while the RLA-I strain showed an improvement in prepulse inhibition (PPI) in adulthood, this was not so for the RHA-I, that retained the same PPI as seen during prepubescence and adolescence, resulting in a lower PPI when tested in adulthood .
In conclusion, this study contributes further into adding construct validity to the RHA-I rat strain as a genetic neurodevelopmental model of schizophrenia-relevant features. This is important as this model provides us with a tool for understanding the maturational trajectory leading to the neurobiological underpinnings behind the behavioural and cognitive disturbances associated with schizophrenia. It can help us test and develop early therapeutic interventions, as shown by a beneficial effect of neonatal handling on social interaction and volumetric brain changes in adult RHA-I rats ( Rio-Alamos et al., 2018 ; Sampedro-Viana et al., 2021 ), as well as an attenuating effect of oxytocin administration on the PPI deficits of the RHA-I strain . The next steps should aim at elucidating how we can intervene in the developmental trajectory of synaptic maturation and at what time window of neurodevelopment this will be most effective and relevant.

Role of funding sources
This study was supported by a scholarship from "Laege Sofus Carl Emil Friis og Hustru Olga Doris Friis'", and partially supported by grants PID2020-114697GB-I00 and 2017SGR-1586 (A.F-T). I.O. is supported by a "Juan de la Cierva" postdoctoral fellowship (FJC2018-038808-I; MCIN/AEI). The funding bodies had no further role in the study design; in the collection, analysis and interpretation of data; in the writing of the report; and in the decision to submit the paper for publication.

Declaration of Competing Interests
We have no conflicts of interest to declare.