Transcriptome Analysis of Duck and Chicken Brains Infected with Aquatic Bird Bornavirus-1 (ABBV-1)

Aquatic bird bornavirus 1 (ABBV-1) is a neurotropic virus that infects waterfowls, resulting in persistent infection. Experimental infection showed that both Muscovy ducks and chickens support persistent ABBV-1 infection in the central nervous system (CNS), up to 12 weeks post-infection (wpi), without the development of clinical disease. The aim of the present study was to describe the transcriptomic profiles in the brains of experimentally infected Muscovy ducks and chickens infected with ABBV-1 at 4 and 12 wpi. Transcribed RNA was sequenced by next-generation sequencing and analyzed by principal component analysis (PCA) and differential gene expression. The functional annotation of differentially expressed genes was evaluated by gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis. The PCA showed that the infected ducks sampled at both 4 and 12 wpi clustered separately from the controls, while only the samples from the chickens at 12 wpi, but not at 4 wpi, formed a separate cluster. In the ducks, more genes were differentially expressed at 4 wpi than 12 wpi, and the majority of the highly differentially expressed genes (DEG) were upregulated. On the other hand, the infected chickens had fewer DEGs at 4 wpi than at 12 wpi, and the majority of those with high numbers of DEGs were downregulated at 4 wpi and upregulated at 12 wpi. The functional annotation showed that the most enriched GO terms were immune-associated in both species; however, the terms associated with the innate immune response were predominantly enriched in the ducks, whereas the chickens had enrichment of both the innate and adaptive immune response. Immune-associated pathways were also enriched according to the KEGG pathway analysis in both species. Overall, the transcriptomic analysis of the duck and chicken brains showed that the main biological responses to ABBV-1 infection were immune-associated and corresponded with the levels of inflammation in the CNS.


Introduction
Bornaviruses are enveloped viruses with a negative-sense single-stranded RNA genome, which are classified in the Bornaviridae family, Orthobornavirus genus [1]. The RNA genome of these viruses encodes at least six proteins, which are, from the 3 to 5 genomic ends, nucleoprotein (N), phosphoprotein (P), non-structural protein (X), matrix protein (M), glycoprotein (G), and RNA-dependent RNA polymerase (L) [2]. Bornaviruses are characterized by nuclear replication, the establishment of persistent infection (non-lytic cycle), and marked neurotropism in infected species, leading to chronic inflammation of the central and peripheral nervous system, resulting in neurological signs [3].
In mammalian species, Borna disease virus 1 (BoDV-1, the type-species of the Orthobornavirus genus, causes Borna disease, a neurological affliction of sheep and horses in central Europe [4]. Aquatic bird bornavirus-1 (ABBV-1), in the Orthobornavirus avisaquaticae species, has been identified in wild waterfowl in North America and Europe, and it has

Principal Component and Differential Gene-Expression Analysis
Principal component analysis (PCA) was carried out to evaluate the contribution of infection status (i.e., treatment) and week post-infection (i.e., time) to explaining the variability of the dataset. PCA plots were created using DESeq2 on the transformed geneexpression estimates using the ggplot2 R library. Venn diagrams were created using the R package, ggVenn v0.1.8.
Differential gene expression (DGE) analysis was conducted by comparing relative gene expression of infected birds with the time-matched control group using the DESeq2 R package, v1.30.0 [30]. First, we examined the total number of differentially expressed genes (tDEGs), as defined by adjusted p-values < 0.05 and absolute log 2 fold change > 0 (|log 2 FC| > 0), which is equivalent to absolute fold change > 1 (|FC| > 1). Further, we evaluated the genes that were highly differentially expressed (hDEG), as defined by adjusted p values < 0.05 and absolute log 2 fold change > 2 (|log 2 FC| > 2), which is equivalent to absolute fold change > 4 (|FC| > 4). The high threshold (|log 2 FC| > 2) was used to increase reproducibility, since replicability of gene calling is correlated with the magnitude of their differential expression [31,32]. The p-values were corrected using the Benjamini-and-Hochberg method for controlling the false-discovery rate [33].

Functional Analysis of DEG using Gene Ontology (GO) Terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Enrichment
GO terms and KEGG pathway enrichment analysis were used to investigate the functional significance of the hDEG (|log 2 FC| > 2). GO enrichment analysis was performed using the GOseq R package, v1.42.0 [34]. GO terms with Benjamini-and-Hochberg-adjusted p-values < 0.05 were considered significantly enriched. KEGG pathway enrichment anal-Viruses 2022, 14, 2211 4 of 21 ysis was performed using the g:Profiler web server version e105_eg52_p16_e84549f [35]. For chickens, the g:GOSt function from the g:Profiler web server was used to perform KEGG enrichment analysis with Gallus gallus selected as the organism and 0.05 selected as the g:SCS significance threshold. For ducks, KEGG pathway enrichment analysis was performed indirectly, due to limited pathway annotation for duck genes, as follows. First, duck hDEGs were converted to chicken ortholog identifiers using the g:Orth function on g:Profiler. Next, the g:GOSt function was used to perform KEGG pathway enrichment analysis on the chicken orthologs of duck genes, with Gallus gallus selected as the organism and 0.05 selected as the g:SCS significance threshold.

Microscopic Pathology and Magnitude of Virus Replication in the Brains of Chickens and Ducks Infected with ABBV-1 upon Intracranial Administration
The brains of Muscovy ducks and White Leghorn chickens were evaluated for histopathology and quantification of RNA virus copies by RT-qPCR. The group means are reported in Table 1, and individual bird data are reported in Supplemental Table S1.  4 3 ± 0.6 2 ± 0.6 0 3.75 ± 1.07 1 Data adapted from Iverson et al. [15]. 2 Data adapted from Iverson et al. [16]. 3 log 10 viral RNA copies/150 ng of tissue RNA, average ± standard deviation. 4 Median ± standard deviation of the semi-quantitative severity of microscopic inflammation (0-6).
No birds developed neurological signs during the course of the experiment, with the exception of a single chicken that underwent emergency euthanasia and was not sampled for the RNAseq analysis [15,16]. When present, inflammatory lesions were similar in both species and were characterized by lymphocytic cuffs in the perivascular spaces, occasionally spilling into the adjacent neuroparenchyma. No evidence of neuronal necrosis was observed. The ducks showed a high intensity of inflammation at 4 wpi (average, 3), which decreased at 12 wpi (average, 2). The chickens showed no inflammation at 4 wpi, while at 12 wpi, they had marked inflammation (average, 3.75). The magnitude of the viral replication, as expressed by RNA copy number per 150 ng of total tissue RNA, was similar in the ducks at the two time points, while it was higher in the chickens at 12 wpi compared to at 4 wpi.

Overview of Sequencing, Trimming, and Alignment Analysis of Duck-and Chicken-Brain Tissues Infected with ABBV-1 or Mock Infected
The number of raw reads for the duck samples (n = 28) fell between 66,851,188 and 136,268,072 (median, 98,259,110). After the adaptor sequences and low-quality scorecontaining bases (Phred score < 30) were trimmed from the reads using Trimmomatic [26], more than 99.90% of the raw reads remained as surviving reads in each sample. When the surviving reads were aligned to the Anas platyrhynchos platyrhynchos reference genome (CAU_duck1.0.98), between 72.64% and 80.56% (median, 76.76%) of the reads were aligned, resulting in a mean coverage per base of between 38.44 and 80.03 (median, 53.61). The number of genes covered by the reads was between 9703 and 10,080 (median, 9877) (Supplemental Table S2).
The number of raw reads for the chicken samples (n = 28) fell between 82,253,152 and 170,553,534 (median, 101,211,921). After the adaptor sequences and low-quality score-containing bases (Phred score < 30) were trimmed from the reads using Trimmomatic [26], more than 99.90% of the raw reads remained as surviving reads in each sample. When the surviving reads were aligned on the Gallus gallus reference genome (Gallus gallus.GRCg6a.98), between 91.25% and 96.60% (median, 94.57%) of the reads were aligned (Supplemental Table S3), resulting in a mean coverage per base of between 49.75 and 101.57 (median, 67.06). The number of genes covered by the reads was between 11,697 and 11,990 (median, 11,820.50) (Supplemental Table S3). In the duck cohort, as shown by the PCA plots ( Figure 1A), there was a separation between most of the control and infected ducks, with PC1 and PC2 accounting for 42% and 35% of the variance, respectively. The control birds displayed a tight grouping, while the infected birds had a more dispersed variance. While the timing of the infection could not clearly separate the birds sampled at 4 and 12 wpi, the ducks at 4 wpi grouped farthest away from the control birds along the PC1 axis. In the chicken cohort, PC1 and PC2 accounted for 43% and 28% of the total variance, respectively ( Figure 1B). There was a clear separation between the control group and the infected birds at 12 wpi along the PC1 axis; however, the infected birds at 4 wpi clustered with the controls. Considering the pathology data, the PCA analysis suggested that differences in the transcriptomic profiles between the infected and control birds paralleled the presence of inflammation, as only groups with mononuclear inflammation clustered separately from the controls. On the other hand, these data suggest that ABBV-1 infection alone, without inflammation, did not sufficiently alter the transcriptomic profiles compared to the uninfected controls, as seen in the chickens at 4 wpi.

DGE Analysis
A differential gene-expression analysis was performed separately for each species, to compare transcript-expression levels between infected and age-matched control birds. A summary of the number of differentially expressed genes for both species is reported in Table 2; furthermore, contextualized data with the results from the experimental challenge are represented in Figure 2. The ducks had a greater number of tDEGs at 4 wpi (2283) compared to 12 wpi (979) ( Table 2), with 78.4% and 79.0% of the tDEGs being upregulated at 4 wpi and 12 wpi, respectively. Of the 2283 tDEGs at 4 wpi, 1313 were known, and 970 were novel genes. Of the 979 tDEGs at 12 wpi, 479 were known genes and 500 were novel genes.

DGE Analysis
A differential gene-expression analysis was performed separately for each species, to compare transcript-expression levels between infected and age-matched control birds. A summary of the number of differentially expressed genes for both species is reported in Table 2; furthermore, contextualized data with the results from the experimental challenge are represented in Figure 2. The ducks had a greater number of tDEGs at 4 wpi (2283) compared to 12 wpi (979) ( Table 2), with 78.4% and 79.0% of the tDEGs being upregulated at 4 The number of hDEGs was 797 at 4 wpi (34.9% of tDEGs) and 228 at 12 wpi (23.5% of tDEG). Of these, 217 hDEGs were shared between 4 and 12 wpi ( Figure 3A). Nearly all (98.4%) and 100% of the hDEGs were upregulated at 4 and 12 wpi (Table 2), respectively, and all of the 217 DEGs shared between the time points were also upregulated.
The chickens had 4328 tDEGs at 4 wpi and 4718 tDEGs at 12 wpi. At 4 wpi, more genes were downregulated (58.5%), while at 12 wpi, this pattern reversed, and there were more upregulated genes (58.9%) ( Table 2). Of the 4328 tDEGs at 4 wpi, 3325 were known genes and 1003 were novel genes. Of the 4718 tDEGs at 12 wpi, 3506 were known genes and 1213 were novel genes. The number of hDEGs was drastically reduced, constituting only 3.5% (n, 153) and 21.1% (n, 997) of the tDEGs at 4 ad 12 wpi, respectively. At 4 wpi, most (94.8%) of the hDEGs were downregulated, while at 12 wpi, 98.2% of the hDEGs were upregulated. Interestingly, only 58.9% of the tDEGs were upregulated at this time point  (Table 2). Only 11 of the hDEGs were shared between 4 and 12 wpi ( Figure 3B), including both up-and downregulated genes.
The chickens had 4328 tDEGs at 4 wpi and 4718 tDEGs at 12 wpi. At 4 wpi, more genes were downregulated (58.5%), while at 12 wpi, this pattern reversed, and there were more upregulated genes (58.9%) ( Table 2). Of the 4328 tDEGs at 4 wpi, 3325 were known genes and 1003 were novel genes. Of the 4718 tDEGs at 12 wpi, 3506 were known genes and 1213 were novel genes. The number of hDEGs was drastically reduced, constituting only 3.5% (n, 153) and 21.1% (n, 997) of the tDEGs at 4 ad 12 wpi, respectively. At 4 wpi, most (94.8%) of the hDEGs were downregulated, while at 12 wpi, 98.2% of the hDEGs were upregulated. Interestingly, only 58.9% of the tDEGs were upregulated at this time point (Table 2). Only 11 of the hDEGs were shared between 4 and 12 wpi ( Figure 3B), including both up-and downregulated genes.

Figure 2.
Overview of viral RNA copies and pathology scores in the brains of ducks and chickens alongside total and highly differentially expressed genes in those tissues. In (A), data are presented for Muscovy ducks at 4 wpi and 12 wpi. In (B), data are presented for chickens at 4 wpi and 12 wpi.
The left y-axis shows both the log10 RNA copies and pathology scores. The right y-axis shows the number of differentially expressed genes. The x-axis shows the weeks post-infection. tDEG stands for total differentially expressed genes and hDEG stands for highly differentially expressed genes. The viral-RNA and pathology-score data used in this figure are adapted from [15,16]. * DEG = differentially expressed genes; ** log2FC = log2 fold change; ^ indicates the total number of tDEGs or hDEGs and the percentage of tDEGs or hDEGs that are upregulated or downregulated are in parenthesis; # indicates the total number of hDEGs and the percentage of tDEGs that are hDEGs are in parenthesis.

Figure 2.
Overview of viral RNA copies and pathology scores in the brains of ducks and chickens alongside total and highly differentially expressed genes in those tissues. In (A), data are presented for Muscovy ducks at 4 wpi and 12 wpi. In (B), data are presented for chickens at 4 wpi and 12 wpi.
The left y-axis shows both the log 10 RNA copies and pathology scores. The right y-axis shows the number of differentially expressed genes. The x-axis shows the weeks post-infection. tDEG stands for total differentially expressed genes and hDEG stands for highly differentially expressed genes. The viral-RNA and pathology-score data used in this figure are adapted from [15,16]. In summary, there were differences in the global differential gene-expression pat- the brains of ABBV-1-infected Muscovy ducks at 4 wpi were compared with those at 12 wpi. The majority of hDEGs occurred at 4 wpi, and 217 hDEGs were shared between 4 and 12 wpi. In (B), the hDEGs in the brains of ABBV-1-infected chickens at 4 wpi were compared with those at 12 wpi. The majority of hDEGs occurred at 12 wpi, and only 11 hDEGs were shared between the two time points.
In summary, there were differences in the global differential gene-expression patterns of the duck and chicken brains in response to ABBV-1 infection, mostly in association with the degree of encephalitis. For the ducks, the number of tDEGs and hDEGs was highest at 4 wpi, when inflammation was high, and decreased at 12 wpi, when lower inflammationseverity scores were registered. Regardless of the time point, most of the hDEGs were upregulated in the ducks. For the chickens, this difference was noted mainly for the hDEGs, which dramatically increased from 4 wpi, when no inflammation was detected, to 12 wpi, the time of high inflammation. Notably, the majority of the tDEGs (58.5%) and hDEGs (94.8%) were downregulated in the chickens at 4 wpi, when infection was established, but inflammation was not detected.

Highest-Ranking hDEGs (|log 2 FC| > 2) in Duck and Chicken Brains in Response to ABBV-1 Infection Are Associated with Immune Functions, Regulation of Gene Expression, and the Cell-Membrane Location
We further examined the 20 highest-ranking known and unknown (novel) hDEGs in the duck and chicken brains at 4 and 12 wpi. The known genes, along with their descriptions (as provided in the Ensembl database), are reported in Tables 3 and 4; the unknown (novel) genes and either their descriptions (as provided on the Ensembl database) or their corresponding protein names (as provided in the UniProt database) are listed in Supplemental Tables S4 and S5. Not every group produced at least 20 hDEGs. The genes of interest are mentioned below.   Table 3). Fourteen of the 20 highest-ranking upregulated hDEGs are associated with immune functions: LAG3, CXCR5, JCHAIN, CTLA4, AICDA, TRAT1, SERPING1, IFNG, XCR1, ZAP70, IL22RA2, VPREB3, PAX5, and CD79B. Four are associated with the regulation of gene expression: SPIC, HOXA4, PLAC8, and APOBEC1 (Table 3). Of the five downregulated hDEGs, three are associated with the cell-membrane region (CDH15, TMEM233, and ATP6V0D2), and two are associated with the regulation of gene expression (ZIC1 and ZIC4) ( Table 3). Considering the novel genes, there were more than 20 upregulated hDEGs and only 8 downregulated hDEGs (Supplemental Table S4). Sixteen of the twenty highest-ranking upregulated novel genes are associated with immune function, as evidenced by the generic protein names, such as 'Ig-like domain-containing protein' or 'SCY domain-containing protein', on the UniProt database (Supplemental Table S4). Eight of the downregulated genes are labeled as 'long non-coding RNA' in the Ensembl database (Supplemental Table S4).

Chickens
At 4 wpi, there were only five upregulated hDEGs but more than twenty downregulated hDEGs (Table 4). Of the upregulated genes, four are associated with immune functions: IRF4, IFIT5, RSAD2, and LYG2. For the downregulated genes, two are associated with responses to pathogens (VWCE and NOS1) and eight are associated with the cell-membrane location: SLC16A8, BAIAP2L2, CGN, KCNH6, SLC10A4, GPRC5C, KCNK3, and DMP1. Considering the novel genes, there were only three upregulated hDEGs, but more than twenty downregulated hDEGs (Supplemental Table S5). For the three upregulated genes, either they were labeled as 'long non-coding RNA' on the Ensembl database, or their proteins were labeled as 'uncharacterized protein' on the UniProt database. For the twenty highest-ranked downregulated hDEGs, three are associated with immune functions: ENSGALG00000053055 (B30.2/SPRY domain-containing protein), ENSGALG00000054795 (TED_complement domain-containing protein), and ENSGALG00000007740 (Wiskott-Aldrich-syndrome-protein-family member). Five are labeled as 'long non-coding RNA', according to the Ensembl database, and four encode 'uncharacterized proteins', according to the UniProt database. The rest of the novel genes could not be characterized further. At 12 wpi, there were more than 20 upregulated hDEGs, but only 10 downregulated hDEGs ( Table 4). Fifteen of the twenty highest-ranking upregulated hDEGs are associated with immune functions: IFNG, IL18RAP, JCHAIN, LAG3, CTLA4, CCL5, CXCL13L3, CRTAM, CCL1, XCL1, CCL19, POU2AF1, CXCL13L2, GZMA, and AICDA (Table 4). Of the ten downregulated genes, three are associated with the cell-membrane region: BAIAP2L2, KCNH8, and KCNK3. Considering the novel genes, there were more than 20 upregulated hDEGs, but only 8 downregulated hDEGs (Supplemental Table S5). Nineteen of the twenty highest-ranked upregulated genes are associated with immune functions and possess generic protein names, such as 'Ig-like domain-containing protein' or 'Immunoglobulin lambda like polypeptide 1 , on the UniProt database (Supplemental Table S5). Of the downregulated genes, one (ENS-GALG00000029168) encodes a 'collectrin domain-containing protein', and four are labeled as 'long non-coding RNA' on the Ensembl database. One novel gene (ENSGALG00000043906), which was downregulated at both 4 wpi and 12 wpi, encodes the protein 'chicken D-serine dehydratase'. The rest of the novel genes could not be characterized further.

ABBV-1 Differentially Regulates Long Non-Coding RNA (lncRNA) Expression in Both Duck and Chicken Brains
The examination of the novel hDEGs by searching each against the Ensembl database (Ensembl release 106: April 2022) showed that many were lncRNAs. A similar examination of the known hDEGs showed no lncRNAs among them. The total number of highly differentially expressed novel lncRNAs is reported in Table 5, and the 20 highest-ranked are shown in Supplemental Table S6. In the ducks at 4 wpi, 95 lncRNAs were highly differentially expressed, of which 87 were upregulated and 8 were downregulated. At 12 wpi, 22 lncRNAs were highly differentially expressed, all of which were upregulated. In the chickens at 4 wpi, only 13 lncRNAs were highly differentially expressed, of which 2 were upregulated and 11 were downregulated. At 12 wpi, 123 lncRNAs were highly differentially expressed, of which 119 were upregulated and 4 were downregulated (Table 5). Overall, the differential expression of the lncRNAs followed the transcriptional dynamics of the other genes and the magnitude of encephalitis, with the highest upregulation in the ducks at 4 wpi and in the chickens at 12 wpi. Some of the lncRNAs showed upregulation of up to 7.33-and 7.95-log 2 fold in the ducks and chickens, respectively. The functional gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses of these novel lncRNAs using g:Profiler (version e106_eg53_p16_65fcd97) showed no enrichment of either the GO terms or the KEGG pathways.

Gene Ontology (GO) Analysis of hDEGs Demonstrates Enrichment of GO Terms Broadly
Associated with Immune Functions 3.6.1. Ducks The enrichment analysis of the hDEGs in the duck brains at 4 wpi showed that there was a total of 59 significantly enriched GO terms: 38 in 'biological process', 9 in 'cellular component', and 12 in 'molecular function'. The 20 highest-ranked enriched GO terms are shown in Figure 4A. Twelve of these are broadly associated with immune functions, antipathogen defense, and transmission signaling ('biological process' class); the terms with the highest number of hits included 'signal transduction' (n, 77), 'immune response' (n, 53), and 'G protein-coupled receptor signaling pathway' (n, 49). The 'cellular component' class included the two most enriched GO terms: 'membrane' (n, 267) and 'integral component of membrane' (n, 247). The 'molecular function' class included enriched GO terms associated with cell signaling (Figure 4A), likely as part of the activated immune response.   At 12 wpi, there were only 14 significantly enriched GO terms: 11 in 'biological process', 1 in 'cellular component', and 2 in 'molecular function' ( Figure 4A). Overall, the terms were less enriched compared to at 4 wpi, and tended to have a lower statistical significance. These GO terms are broadly associated with immune functions, anti-pathogen defense, or cell signaling. The membrane-and extracellular-space-associated GO terms, markedly enriched at 4 wpi, were no longer so at 12 wpi.

Chickens
The enrichment analysis of the hDEGs in the chicken brains at 4 wpi showed no significantly enriched GO terms, while at 12 wpi, there were a total of 107 significantly enriched GO terms: 82 in 'biological process', 8 in 'cellular component', and 17 in 'molecular function'. The 20 highest-ranked enriched GO terms are shown in Figure 4B

Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Analysis of hDEG Demonstrates Enrichment of Pathways Associated with Immune Functions or Viral Infections
The KEGG pathway analysis of the hDEGs for the chicken brains at 4 wpi did not show any enrichment. The KEGG pathway analysis of the hDEGs for the ducks at 4 wpi, the ducks at 12 wpi, and the chickens at 12 wpi showed the enrichment of seven common pathways ( Table 6) between all three groups. These were 'cytokine-cytokine receptor interaction', 'toll-like receptor signaling pathway', 'herpes simplex virus 1 infection', 'cell adhesion molecules', 'influenza A', 'NOD-like receptor signaling pathway', and 'phagosome'. Among these pathways, the highest significance was observed in the ducks at 4 wpi and the chickens at 12 wpi, probably as a reflection of the inflammatory reactions in the brains. Five additional pathways were not common to all the groups and were only enriched in some (Table 6). Overall, except for the 'cell adhesion molecules' and 'neuroactive ligandreceptor interaction' pathways, all of the significantly enriched KEGG pathways are directly associated with immune functions or responses to virus infections.

Discussion
We evaluated the transcriptomic profiles of aquatic bird bornavirus -1 (ABBV-1) in the infected brains of the Muscovy ducks and chickens at 4 and 12 wpi. We opted for these two time points as 4 wpi was the earliest time point at which all the birds were infected and had a relatively high magnitude of viral RNA copies in their brain tissues. The last time point was included to model chronic infection, as seen in the field. In choosing 4 wpi as our first time point in the analysis, we decreased our ability to assess transcriptomic changes that were not superimposed with inflammation (see below). Nonetheless, the use of brains with established infection was considered the most stringent inclusion criterion for this study. Moreover, age-matched controls were included to account for transcriptomic changes associated with the fast growth of poultry during their first few weeks of age.
For the entire data set, the number of aligned (mapped) reads from the samples was between 53,167,653 (lowest) and 157,453,838 (highest). The number of aligned reads per sample in the data set was well above the recommended number, approximately 36 million, for the accurate determination of the abundance in 80% of the genes [36]. The percentage of aligned reads for the chickens was between 91% and 97%, and for the ducks, it was between 70% to 81%. The reduced alignment for the Muscovy-duck samples was probably due to the use of the Anas platyrhynchos platyrhynchos (Pekin duck) reference genome instead of the Cairina moschata reference genome. The latter was not used because it is relatively recent and poorly annotated compared to the Pekin-duck reference genome, thus limiting downstream analyses. Nonetheless, a percentage of aligned reads between 70% and 90% is acceptable for further analysis [37]. The use of 7 biological replicates (i.e., 7 brain samples) per experimental unit (group/time point) is in line with the recommend minimum of 6 replicates per condition to identify differentially expressed genes [38].
Changes in the transcriptomic profiles in response to ABBV-1 infection appeared earlier in the duck brains than in those of the chickens, as shown by the PCA and DGE analyses. In the PCA analysis, the infected duck-brain tissues were clearly separated from control-brain tissues both at 4 and 12 wpi, whereas the infected chicken-brain tissues were only separated from the controls at 12 wpi. While the chicken brains had a greater number of tDEGs than the ducks at 4 wpi (4328 vs. 2283), when considering hDEG, only 153 genes were upregulated in the chicken brains (3.5% of tDEG), compared to 797 in those of the ducks (34.9% of the tDEGs). Furthermore, by 12 wpi, the transcriptome response of the ducks decreased to 979 tDEGs (228 hDEGs), while in chickens, it increased to 4718 tDEGs (997 hDEGs). This pattern correlated with the development of inflammation in the brains. The ducks had a higher brain-pathology score at 4 wpi (3.0 score) than at 12 wpi (2.0 score), while in the chickens, inflammation did not develop at 4 wpi (0 score) but did at 12 wpi (3.75 score) [15,16]. The functional evaluation of the upregulated genes further indicated that the transcriptomic changes in the brains of the ducks at both time points and in those of the chickens at 12 wpi were driven by inflammation, as most (≥70%) of the 20 most upregulated hDEGs were associated with immune functions. The reduced ABBV-1 viral load in the chicken brains could explain the lower degree of inflammation and the less-significant transcriptomic changes. In fact, at 4 wpi, the duck brains had an average of 10 6.87 viral RNA copies per 150 ng of total RNA, while the chicken brains had approximately 10 4.97 [15,16]. By 12 wpi, the average viral RNA copy in the duck brains remained relatively steady at 10 7.15 , while it increased to 10 6.56 in those of the chickens, where it was approximately 39 times higher than at 4 wpi. The lower ABBV-1 replication in the chickens relative to the ducks was not surprising, since the same was also true in vitro, where ABBV-1 was propagated less efficiently in chicken compared to duck cell lines [17,39]. Moreover, the lower concentration of inoculum in the chickens might have played a role.
The chickens at 4 wpi had a much greater proportion of downregulated tDEGs (58.5%) and hDEGs (94.8%) compared to the chickens at 12 wpi (41.1% of tDEGs and 1.8% of hDEGs) and the ducks at both 4 wpi (21.6% of tDEGs and 1.6% of hDEGs) and 12 wpi (21.0% of tDEGs and 0% of hDEGs). Numerous (8/20) downregulated hDEGs in the chicken brains at 4 wpi were associated with cell-membrane regulation, which suggests that ABBV-1 infection alone, without superimposed inflammation, may alter cell-membrane homeostasis. However, this differential gene expression was not sufficient to separate the infected chickens at 4 wpi from the controls in the PCA analysis and did not significantly enrich the GO terms or KEGG pathways. Furthermore, in the other groups, some of the downregulated hDEGs were associated with the cell membrane or cell-membrane/extracellular-matrix region, although this proportion varied (3/5 for the ducks at 4 wpi, 0 for the ducks at 12 wpi, and 3/10 for the chickens at 12 wpi). This may indicate a separate contribution of ABBV-1 infection to transcriptomic changes, distinct from inflammation. In a microarray experiment using brains from rats infected with BoDV-1 [18], cell-membrane/extracellular-matrixregion-associated genes, such as Ca-transporter ATPase3, Na/K transporting ATPase beta 1 subunit, Ras-guanine nucleotide release/exchange factor p140, and vacuolar ATP synthase 16-kDa proteolipid subunit (ATP6V0C (ATPase H+ transporting V0 subunit C)), were downregulated [18], supporting our observations. For example, the ATP6V0D2 (ATPase H+ transporting V0 subunit d2) gene was highly downregulated in ducks at 4 wpi. The other downregulated genes in our study included CHRNA2 (Cholinergic receptor nicotinic alpha 2 subunit), KCNH6 (potassium voltage-gated channel subfamily H member 6), SLC10A4 (solute carrier family 10 member 4), CBLN1 (cerebellin 1 precursor), NOS1 (nitric oxide synthase 1), and SHISA6 (shisa family member 6), which are associated with neurotransmission and synaptic physiology, suggesting the downregulation of these neuronal activities by ABBV-1. Furthermore, the BoDV-1 infection of rat neurons was shown to impair neuronal remodeling and synaptic plasticity [22,40].
When inflammation was present (ducks at 4 and 12 wpi, and chickens at 12 wpi), there was a marked upregulation of immune-check-point genes in both species; specifically, LAG3 (lymphocyte activating 3) and CTLA4 (cytotoxic T-lymphocyte-associated protein 4). LAG3 is a membrane receptor protein expressed on T cells, which binds to the MHC class II receptor, preventing its binding with the T-cell receptor (TCR) and CD4 [41]. This leads to the suppression of TCR signaling and, consequently, the suppression of CD4 + and CD8 + T cell activation. CTLA4 is a membrane receptor protein that is constitutively expressed in regulatory T cells [42]. It can competitively bind to CD80/86 on antigen-presenting cells (APCs), preventing its interaction with the CD28 on T cells, leading to the inhibition of CD4 + T-cell activation [42]. This suggests that ABBV-1-infected brain tissues may develop a mechanism to control neuroinflammation and prevent the development of functional T cells, despite their inability to prevent the accumulation of lymphocytes in the neuroparenchyma. This could explain the lack of neurological signs despite the morphological evidence of encephalitis, as well as the lack of neuronal necrosis. By contrast, the neuroinflammation caused by other orthobornaviruses, such as BoDV-1, variegated squirrel bornavirus 1 (VSBV-1), and parrot bornavirus-2 (PaBV-2), can result in disease development driven by T-cell-dominated immunopathogenesis, often resulting in histological evidence of neuronal degeneration [2,14,[43][44][45]. In one genome-wide RNAseq study, the BoDV-1 infection of newborn rats led to neurological signs in the animal, and the transcriptomic analysis of their brains showed the upregulation of inflammatory cytokines [23]. Additionally, KEGG pathways similar to those enriched in duck and chicken brains infected with ABBV-1, such as 'cytokine-cytokine receptor interaction', 'toll-like receptor signaling pathway', 'NOD-like receptor signaling pathway', and 'intestinal immune network for IgA production', were enriched in the rat brains [23]. However, very few studies have evaluated the transcriptomic profiles of viral infection in avian brains. In ducks infected with highly pathogenic H5N1, the development of neurological signs was correlated with the upregulation of immuneassociated genes, such as PRRs, IFNs, and cytokines, with no upregulation of regulatory genes, such as LAG3 or CTLA4 [46]. Future studies should evaluate the expression of immune-checkpoint genes or the relative amounts of regulatory T cells in the perivascular lymphocytic cuffs of the nervous tissue, as possible drivers of asymptomatic disease in orthobornavirus infection.
Numerous novel lncRNAs were highly differentially expressed in the brains of the ducks and chickens infected with ABBV-1; however, their functions are unknown, as these genes did not yield any significantly enriched GO terms or KEGG pathways. This is likely to have been a consequence of the only partial annotation of lncRNAs in the respective reference genomes of ducks (CAU_duck1.0) and chickens (GRCg6a). Although the function of lncRNAs in cells is broad, with roles in many diverse biological processes [47,48], most of the lncRNAs in our experimental cohort appeared to correlate with the severity of inflammation. In fact, the differential regulation of the lncRNAs was higher in the ducks at 4 wpi and in the chickens at 12 wpi, and lower in the ducks at 12 wpi and in the chickens at 4 wpi. This suggests that the differential expression of the lncRNAs occurred in the inflammatory cells rather than the neuroparenchyma, which is consistent with the established role of lncRNAs in regulating the inflammatory response in mammals [49][50][51]. Nonetheless, some of the lncRNAs were differentially expressed either in the absence of inflammation (chicken 4 wpi), or when only mild inflammation was present (ducks 12 wpi), suggesting differential expression in virus-infected neurons or glial cells. A genome-wide microarray study of primary cultures of mouse cortical neurons showed the differential regulation of 3528 and 2661 lncRNAs upon infection with BoDV-1 Strain V and Hu-H1, respectively [24]. A functional analysis showed that these lncRNAs were involved in numerous processes, including metabolic and biological regulation, cell adhesion, endocytosis, cancer, and viral infections [24]. The relationship between lncRNAs and bornavirus infection is an interesting area of future investigation.
The GO analysis of both the ducks and the chickens produced many enriched GO terms, of which most (≥60%) of the 20 most enriched were broadly associated with immune functions or anti-pathogen defense, while the remainder were associated with the membrane or extracellular location. The enrichment of these GO terms seemed to correlate with the severity of the inflammation and dovetailed with the number of DEGs. The ducks at 4 wpi had more enriched GO terms (59) than the ducks at 12 wpi (14), and in the chickens, the GO terms were only enriched at 12 wpi, as no inflammation was observed at 4 wpi. The GO terms in the cellular-component class were broadly associated with the cell membrane and extracellular space and included some of the most enriched terms in the whole analysis (i.e., 'membrane', 'integral component of membrane', and 'external side of the plasma membrane'). This may have been due, at least in part, to the links between some of the immune-associated genes and multiple GO terms. For example, the upregulated CTLA4 in the duck brains at 4 wpi is linked to five GO terms that are both related to immune response and the plasma membrane (i.e., 'immune response', 'integral component of membrane', 'membrane', 'B cell receptor signaling pathway', and 'external side of plasma membrane'), and the IFNG gene is linked to 10 GO terms. Several membrane-associated genes were downregulated in the brains of the chickens at 4 wpi; however, these did not produce significantly enriched GO terms, and it is unclear whether ABBV-1 alone caused alterations in the membrane homeostasis. When the inflammation in the duck brains at 12 wpi decreased, the enrichment of the terms associated with the membrane and extracellular space was drastically reduced, despite the fact that the birds did not have high levels of viral RNA.
One difference between the GO profiles of the ducks and the chickens was the enrichment of the terms associated with adaptive immunity. GO terms such as 'immunoglobulin production', 'positive regulation of B cell activation', 'B cell receptor signaling pathway', 'immunoglobulin complex, circulating', 'antigen binding', and 'immunoglobulin receptor binding' were enriched in the chickens, whereas the ducks had 'MHC class II protein complex' as the only GO term clearly associated with adaptive immunity. The severity of the inflammation in the chickens at 12 wpi was similar to or greater than that observed in the ducks at 4 wpi, and in both species, the amount of Pax-5-positive cells (B lymphocytes) in the perivascular cuffs was low (<~10-15%) [15,16]. Therefore, the reason why the chickens had a greater proportion of GO terms associated with B-cell activation is unknown, but it may be related to species-specific differences in the immune response to viral infection. The presence of GO terms associated with adaptive immunity agrees with the evidence of seroconversion against ABBV-1 N protein in both species [15,16].
The KEGG analysis showed the enrichment of 12 pathways in ducks at 4 wpi, 9 in ducks at 12 wpi, and 10 in chickens at 12 wpi. Seven common KEGG pathways were enriched between all the groups. As with the GO analysis, the KEGG pathway analysis showed the enrichment of pathways mainly related to immune functions or viral infections, with the highest significance in the groups with the most inflammation (the ducks at 4 wpi and the chickens at 12 wpi). Multiple pattern-recognition receptor (PRR) pathways were enriched, including the 'toll-like receptor-signaling pathway', 'NOD-like receptor-signaling pathway', 'cytosolic DNA-sensing pathway', and 'RIG-I-like receptor-signaling pathway'. The 'toll-like receptor-signaling pathway', 'NOD-like receptor-signaling pathway', and 'cytosolic DNA-sensing pathway' were also enriched in the brains of the rats after infection with BoDV-1 [23], corroborating our findings. Some of the DEGs enriched by ABBV-1 infection in these pathways included numerous PRRs genes, such as TLR1A (binds diand triacylated lipopeptides), TLR3 (binds dsRNA), TLR4 (binds bacterial lipopolysaccharide), TLR7 (binds ssRNA), DHX58 (also known as (aka) LGP2 and binds to dsRNA), IFIH1 (encodes MDA5 and regulated by DHX58), and TMEM173 (aka STING1 and binds DNA) [52][53][54][55][56][57]. While TLR4 is commonly known to be activated by bacterial lipopolysaccharide [54], it can be activated by multiple RNA viruses in different families [58] and by BoDV-1 [23]. TMEM173 is a cytosolic DNA sensor that binds DNA but can also be activated by RNA viruses through multiple proposed mechanisms [57]. Therefore, our demonstration of TLR4 and TMEM173 upregulation in the context of ABBV-1 is not without precedent.
The activation of PRRs leads to signaling cascades of multiple pathways, involving several proteins. Facilitator genes, such DDX60, TRIM25, and TRIM14 were differentially expressed in response to ABBV-1 infection and have important roles in either the activation or propagation of the signaling cascades [59][60][61][62]. Signal transductions eventually result in, but are not limited to, the activation of IRF-7, which then induces the expression of IFNα/β [63]. IFN-α/β initiates the JAK-STAT signaling-cascade pathway by binding to IFNα/β receptors (IFNAR) on the cell membrane [64]. In our study, ABBV-1 infection resulted in the differential expression of IRF-7, IFN-A, STAT1, and STAT2. Similarly, the BoDV-1 infection of rat brains led to the upregulation of IRF-7 and STAT1 transcripts [20] and KEGG enrichment of the JAK-STAT signaling pathway [23]. The IFN JAK-STAT pathway induces the expression of a group of proteins collectively called interferon stimulated genes (ISG), some of which possess antiviral activities [65]. The ISGs induced by ABBV-1 in both Muscovy ducks and chickens included OASL, RSAD2 (aka viperin), and MX1, all of which possess antiviral activities [65]. Similarly, Mx1 transcripts, OASL2 transcripts, and Mx1 proteins were also previously found to be upregulated in mice embryos and rat brains infected with BoDV-1 [20,66].
The 'intestinal immune network for IgA production' KEGG pathway was enriched in both chickens at 12 wpi (p-value adj. 6.503 × 10 −22 ) and ducks at 4 wpi (p-value adj. 4.814 × 10 −8 ), corresponding with the time points of severe inflammation. However, the statistical significance level of this enrichment was much greater in the chickens at 12 wpi than in the ducks at 4 wpi, which correlates with the enriched GO terms in the chickens associated with B-cell development and humoral immunity. Therefore, IgA production against ABBV-1 proteins was probably stimulated in the chickens and ducks. The BoDV-1 infection of rats also led to the enrichment of the 'intestinal immune network for IgA production' KEGG pathway, along with the 'antigen processing and presentation' KEGG pathway [23]. This further demonstrated the capacity of bornavirus to induce IgA in hosts.
The 'necroptosis' (or inflammatory cell death) KEGG pathway was enriched in the brains of infected ducks at 4 wpi, when the inflammation was severe. The enrichment of this pathway was further supported by the enrichment of both the 'influenza A' and 'herpes simplex virus 1 infection' KEGG pathways, as these were two viruses known to modulate necroptosis [67]. Necroptosis possibly occurred in the lymphocytes, as suggested by the histological evidence of the perivascular lymphocytic cuffs, and the lack of histological evidence of damage in neurons or glial cells. Necroptosis could have limited the excessive accumulation of lymphocytes and may explain the lower degree of inflammation in the brains of ducks at 12 wpi, compared to at 4 wpi. Since the immunopathogenesis of bornavirus infection is driven by primed cytotoxic CD8+ T lymphocytes against infected cells [44,[68][69][70][71], the necroptosis of cytotoxic lymphocytes could explain the lack of clinical signs seen in our experimental trial. This explanation, however, remains speculative. To corroborate this mechanism, the expression of necroptosis markers, such as mixed lineage kinase domain-like (MLKL) and receptor interacting protein kinase 1 (RIPK 1) and 3 (RIPK 3), in lymphocytes should be evaluated by immunohistochemistry or in situ hybridization in sections of ABBV-1-infected brains [72].

Conclusions
This study shows that the transcriptomic changes in the brains of both Muscovy ducks and chickens infected with ABBV-1 were driven by the development of non-heterophilic inflammation, as evidenced by histopathology. The majority of highly differentially expressed genes were associated with immune response and cytokine signaling. Additionally, numerous highly differentially expressed novel genes were determined to be lncRNAs. In the groups with high magnitudes of inflammation (the ducks at 4 wpi and the chickens at 12 wpi), the chickens had greater upregulation of the pathways associated with adaptive immune response and B-cell maturation than the ducks.
When infection was established, but inflammation was not developed, as in the case of the infected chickens at 4 wpi, the transcriptomic changes were dominated by downregulated genes; however, no functional pathways were significantly enriched. This suggests that the transcriptomic changes caused by ABBV-1 infection alone are of low magnitude, supporting the notion that bornaviruses are not cytolytic and can establish chronic infection without severe cellular derangement.
Separating infection from inflammation may be of no use in the context of animal studies. While in vivo studies are needed to model infection in the host, in vitro or ex vivo studies may be better suited to assessing changes in the global transcriptomic profiles caused by ABBV-1, without the confounding effect of superimposed inflammation.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/v14102211/s1. Table S1: Summary of individual-bird ABBV-1 RNA copy-number and histology-score data in brains of ducks and chickens; Table S2: Summary of sequencing analysis in infected and control Muscovy ducks; Table S3: Summary of sequencing analysis in infected and control chickens; Table S4: The 20 most highly up-and downregulated novel differentially expressed genes (DEG) in Muscovy ducks; Table S5: The 20 most highly up-and downregulated novel differentially expressed genes (DEG) in chickens; Table S6: The 20 most highly up-and downregulated novel long non-coding RNAs (lncRNAs) in ducks and chickens.