All animal trials were approved by the South China Agriculture University Institutional Animal Care and Use Committee (identification code: 2019076, 10 June 2019). All animal procedures were performed according to the regulations and guidelines established by this committee and international standards for animal welfare.

PBMCs samples from three ALV-J infected chickens (#5, #12, and #15) at 21 DPI were prepared as previously described (Dai et al., 2020). In particular, when compared with uninfected chickens, a CD8^high+ population increased in the ALV-J infected chicken PBMCs at 21 DPI, which formed three stable populations of CD8+ T lymphocytes, including CD8 ^high+, CD8^medium+, and CD4⁺CD8^low+ T lymphocytes. CD4⁺CD8^low+ T also named CD8^low+ T cell in this study. The CD8^high+ T lymphocyte included two phenotypes, CD8αα and CD8αβ (Dai et al., 2020). To compare and analyze the difference of these T lymphocytes populations, we subjected the sorted cell subsets to SMART-Seq2 based scRNA-seq analysis.

A fluorescence-activated cell sorting machine (FACS Aria II, Becton Dickinson, New Jersey, USA) was used to sort a single cell into each well of a 96-well PCR plate containing 2.5μL of 10× Lysis Buffer (Vazyme# N711). Pooled PBMCs were from the mixed equal amount of PBMC samples of three ALV-J infected chickens at 21 DPI. For the isolation of CD8^high+, CD8^medium+, and CD4⁺CD8^low+ (CD8^low+) populations, the pooled PBMCs from ALV-J infected chicken were stained for APC‐conjugated mouse anti‐chicken CD3⁺, FITC‐conjugated mouse anti‐chicken CD4⁺, and PE‐conjugated mouse anti‐chicken CD8α⁺ monoclonal antibodies (Southern Biotech, Birmingham, USA). For CD8^high αα⁺ and CD8^high αβ⁺ population isolation, the pooled PBMCs from ALV-J infected chicken were stained for APC‐conjugated mouse anti‐chicken CD3⁺, PE‐conjugated mouse anti‐chicken CD8α^+, and FITC‐conjugated mouse anti‐chicken CD8β⁺ monoclonal antibodies (Southern Biotech, Birmingham, USA) as described previously (Dai et al., 2020). Each population of five T lymphocyte subtypes was analyzed in four replications, and each replication sorted 100 single cells for subsequent SMART-Seq2 analysis. Furthermore, the gating strategy for each population was shown in Supplementary Figures 1A, B .

Library construction and sequencing were completed by Gene Denovo (Guangzhou, China). Total RNA from the lysed cells was briefly released, enriched by Oligo (dT) primers, and reverse transcribed into cDNA using the Discover-sc WTA Kit V2 (Vazyme). cDNA was purified with VAHTS DNA Clean Beads. Sequencing libraries were constructed according to the TruePrepTM DNA Library Prep Kit V2 (Vazyme #TD503). Lastly, library quality was assessed on an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and sequenced using Illumina Novaseq6000 by Gene Denovo Biotechnology Co. (Guangzhou, China).

Reads obtained from the sequencing machines were further filtered using fastp (Chen et al., 2018) (version 0.18.0). The parameters were as follows: 1) removal of reads containing adapters; 2) removal of reads containing more than 10% of unknown nucleotides (N); 3) removal of low-quality reads containing more than 50% of low-quality bases (Q-value ≤20). Furthermore, the ribosome RNA (rRNA) mapped reads were removed via the short reads alignment tool Bowtie2 (Langmead and Salzberg, 2012) (version 2.2.8). The remaining clean reads were mapped to the chicken genome of GRCg6a (Zerbino et al., 2018) using HISAT2.2.4 (Kim et al., 2015). FPKM (fragment per kilobase of transcript per million mapped reads) value was calculated to quantify the gene expression abundance using RSEM (Li and Dewey, 2011) software. Principal component analysis (PCA) was performed with the R package gmodels (http://www.rproject.org/) to reveal various samples’ repeatability and potential relationship.

Differential expression RNA analysis was performed using DESeq2 (Love et al., 2014) software comparing two different groups. Differentially expressed genes (DEGs) were identified with the parameter of false discovery rate (FDR) <0.05 and absolute fold change (FC) ≥2. The Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa and Goto, 2000), an important public pathway-related database, was used for pathway enrichment analysis. The analysis identified significantly enriched metabolic pathways or signal transduction pathways in DEGs compared with the whole genome background with a false discovery rate threshold (FDR) < 0.05.

Weighted gene co-expression network analysis (WGCNA) is a system biology method used to describe the correlation patterns among genes across multiple samples and identify modules of highly correlated genes. After filtering 9862 genes with low expression in all samples, WGCNA (v1.47) package in R (Langfelder and Horvath, 2008) was used to construct co-expression modules based on the gene expression values. Module Eigengenes were used to calculate the correlation matrices with samples to identify biologically significant modules. Intramodular connectivity of each gene was calculated, and high connectivity tended to be hub genes with potentially important functions. The networks were visualized using Cytoscape (v3.7.1). For genes in each module, KEGG pathway enrichment analysis was conducted to analyze the biological functions of each module.

Gene expression pattern analysis was used to cluster genes sharing similar expression patterns for multiple samples. The expression data of each sample were normalized to 0, log2 (v1/v0), log2 (v2/v0), and then clustered by Short Time-series Expression Miner software (STEM, version 1.3.11) to examine the expression pattern of DEGs (Ernst and Bar-Joseph, 2006). The parameters were set as follows: 1) maximum unit change in model profiles between time points was 1; 2) maximum output profiles number was 20, and similar profiles were merged; 3) a minimum ratio of FC of DEGs was no less than 2.0. The clustered profiles with p-values ≤0.05 were considered as significant profiles. Then the DEGs in each profile were subjected to KEGG pathway enrichment analysis.

A protein-protein interaction (PPI) network was first identified based on WGCNA analysis as described above. Next, the interaction network was supplemented with String v10 (Szklarczyk et al., 2015), which determined genes as nodes and interactions as lines in a network. The network file was visualized using Cytoscape (v3.7.1) (Shannon et al., 2003) software to present a core and hub gene biological interaction.

Primers for quantitative reverse transcription PCR (qRT-PCR) were designed using the National Center for Biotechnology Information (NCBI) Primer-BLAST program. The GAPDH gene was used as an internal control. All primers used in this study are listed in Supplementary File 1 . A total of 2×10⁵ live cells of CD8⁺ T cell or CD4⁺CD8^low+ double-positive T cell from three ALV-J infected chickens (#5, #12, and #15) at 21 DPI were respectively sorted for total RNA extraction. Furthermore, the gating strategy for each population is shown in Supplementary Figure 1C . qRT-PCR was performed on an ABI7500 Real-Time PCR system (Applied Biosystems, USA) using ChamQ Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China) according to the manufacturer’s specifications. Expression levels of the tested reference genes were determined by CT values and calculated by 2^-△△Ct (Livak and Schmittgen, 2001). Data were collected from three biological samples in each group, and each sample was performed in triplicate.

Statistical comparisons were made by GraphPad Prism 8 (GraphPad Software Inc., San Diego, CA, USA). The results were presented as mean ± SEM. The paired or unpaired t-test was used for statistical comparison. *P < 0.05, **P < 0.01, ***P < 0.001.

We used the fluorescent cell sorting based SMART-Seq2 platforms to perform full-length scRNA-seq on five T lymphocytes populations including CD4⁺CD8^low+, CD8^medium+, CD8^high+, CD8^high αα⁺, and CD8^high αβ⁺ T cells in PBMCs from ALV-J infected chicken at 21 DPI ( Figure 1A and Supplementary Figures 1A, B ). Each T cell population obtained in four replicate experiments showed excellent repeatability based on the PCA and Pearson’s correlation analysis ( Figure 1B and Supplementary Figure 2 ). Then we explored the correlation patterns among genes across the twenty samples via WGCNA analysis. Genes were clustered into 22 correlated modules ( Figure 1C and Supplementary File 2 ). Moreover, we found that CD8^high αα⁺ T cells were correlated with the green module in which genes were enriched in the immune pathway “Cytokine–cytokine receptor interaction” ( Figures 1D, E and Supplementary File 3 ), while CD4⁺CD8^low+ and CD8^medium+ T cells were related to the brown module in which genes were enriched in immune pathways “FoxO signaling pathway” and “NOD-like receptor signaling pathway” ( Figures 1F, G and Supplementary File 4 ).

To further investigate the heterogeneous populations of the five chicken T lymphocytes populations, we first performed the trend analysis of CD8^high+, CD8^high αα^+,and CD8^high αβ⁺ T cells. CD8^high+ T cell contained both CD8^high αα⁺ and CD8^high αβ⁺ T cells. The most significant profile, profile 3, displayed a slightly increased then decreased trend in the order of CD8^high+ T cell, CD8^high αα⁺ T cell, and CD8^high αβ⁺ T cell ( Figures 2A, B ). Next, 253 genes in profile 3 were subjected to KEGG enrichment analysis, in which genes were mainly enriched in the “Cytokine–cytokine receptor interaction” pathway ( Figure 2C , Supplementary File 5 ). The heat map results showed that genes enriched in the “Cytokine–cytokine receptor interaction” pathway were highly expressed in the CD8^high αα⁺ T cell population ( Figure 2D and Supplementary File 5 ).

Similarly, trend analysis for CD4⁺CD8^low+, CD8^medium+, and CD8^high+ T cells was performed. Profile 7 displayed a significantly up-regulated pattern, and profile 0 showed a down-regulated pattern in the order of CD4⁺CD8^low+, CD8^medium+ and CD8^high+ T cells ( Figures 3A – C ). In addition, genes in profile 7 were mainly enriched in the “Cytokine–cytokine receptor interaction” and “Toll-like receptor signaling pathway” and were highly expressed in the CD8^high+ T cells population as could be expected ( Figures 3D, F and Supplementary File 6 ). For example, immune-related genes including IL18RAP, CCR6, CCL4, CCL5, FAS, and IFNG were predominantly expressed in the CD8 ^high+ T cells population ( Figure 3F ). Genes in profile 0 were mainly involved in the “FoxO signaling pathway” and “TGF-beta signaling pathway” and were highly expressed in the CD4⁺CD8^low+ (CD8^low+) T cells population ( Figures 3E, F and Supplementary File 6 ).

The seven predominantly expressed genes in the CD8^low+ T cells population include CCND1, ROCK1, FOXO1, Smad7b, FOXO3, S1PR4, and ENSGALG00000037160, and we also conducted an interaction network analysis based on WGCNA analysis and STRING database. The results implied that CCND1, ROCK1, FOXO1, FOXO3, and ENSGALG00000037160 were the important hub genes ( Figure 3G , marked in red), and their regulation on chicken T cells functions had not been reported. In mammals, it has been reported that FOXO1 or FOXO3 plays an essential role in specifying the program of T cell differentiation, most importantly in the pathway leading to the development and function of regulatory T (Treg) cells (Kerdiles et al., 2010; Ouyang et al., 2012). Further, CCND1 has been reported to be negatively correlated with T cell activation (Lu et al., 2019). ROCK was also the key regulator of T-lymphocyte development, activation, and differentiation (Saoudi et al., 2014), and the RhoA-ROCK pathway is the specific pathway for Th1, Th17, and Treg responses (Liu et al., 2020). This information suggests that the function of CD4⁺CD8^low+ T cells may be involved in negatively regulating the activity of T cells, based on the high expression of CCND1, ROCK1, FOXO1, and FOXO3.

To further compare the gene expression characteristics of CD4⁺CD8^low+, CD8^medium+, CD8^high+, CD8^high αα^+,and CD8^high αβ⁺ T cell populations, a total of 23 significant genes enriched in the “Cytokine–cytokine receptor interaction” pathway were screened, and their expression levels were quantified in a heat map ( Figure 4A and Supplementary File 7 ). We also performed an interaction network analysis of these genes based on the WGCNA analysis and the STRING database ( Figure 4B ). The results showed that the hub genes including IL17A, CD40LG, IL2RA, IL2RB, CCR6, CCL20, TNFRSF25, TNFRSF11, and IL1R1 were predominantly expressed in CD8^high αα⁺ T cells ( Figure 4 ). Conversely, hub genes including CD4, IL31RA, and TNFRSF18 were highly expressed in CD4⁺CD8^low+ T cells ( Figure 4 ).

Given the limited cell numbers of CD8^high+, CD8^high αα⁺ and CD8^high αβ⁺ T cells in PBMCs, we could not obtain sufficient cell numbers to perform qRT-PCR detection. Thus, we sorted the CD4⁺CD8^low+ double-positive T cells and CD8⁺ T cells containing mixed CD8^medium+ and CD8^high+ T cell populations for qRT-PCR. The sorting gate strategy is shown in Supplementary Figure 1C . CD4⁺CD8^low+ double-positive T cells and CD8⁺ T cells from three ALV-J infected chickens were sorted respectively to evaluate the reliability of the smart-seq results, and 20 DEGs were selected to validate the relative expression levels in the CD8⁺ T cells compared to CD4⁺CD8 ^low+ T cells using qPCR. As shown in Figures 5A, B , the trends in expression of these randomly selected DEGs were consistent between smart-seq data and qRT-PCR data. Specifically, FOXO3, FOXO1, CCND1, GATA3, and IL31RA had a low expression, IL2RB, CD8A, NK-lysin, IFNG, Granzyme K, Granzyme A, and CCL5 were highly expressed in CD8⁺ T cells compared to CD4⁺CD8^low+ double-positive T cells ( Figures 5A, B and Supplementary File 8 ).