Transcriptome analysis reveals potential immune function-related regulatory genes/pathways of female Lubo goat submandibular glands at different developmental stages

Ethics statement

All animal experiments were approved by the Committee on the Management and Use of Laboratory Animals of Shandong Agricultural University (Permit Number: SDAUA-2018-052) and conducted in accordance with the “Guidelines for Experimental Animals” of the Ministry of Science and Technology (Beijing, China). Every effort was made to minimize pain.

Sample collection and RNA extraction

The collection of experimental animals was carried out according to the principle of the 3Rs (Reduction, Replacement, and Refinement) for animal experimentation. Goats were separated into three different groups based on their developmental stages: the 1-month-old group, the 12-month-old group and the 24-month-old group. Three healthy female goats (Lubo goats) from each group (total, nine goats) were selected as experimental animals. The animals were obtained from Shandong Fengxiang Animal Husbandry and Seed Industry Technology Co., Ltd. (Laiwu District, Jinan, Shandong Province, China). All goats were raised under the same conditions with natural lighting and free access to food. All goats were fasted from food for 24 h and fasted from liquids for 2 h, and then xylazine hydrochloride (Huamu Animal Health Products Co., Ltd., China) was injected intramuscularly to induce general anaesthesia. After slaughter and exsanguination, the submandibular glands were immediately collected surgically, placed in liquid nitrogen and then placed in a −80 °C freezer for long-term storage. To distinguish the samples, the submandibular salivary gland samples from 1-month-old goats were labelled A1-L, A2-L and A3-L; the samples from 12-month-old goats were labelled B3-L, B4-L and B5-L; and the samples from 24-month-old goats were labelled C2-L, C3-L and C5-L. Total RNA was extracted from the 9 submandibular gland samples using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). The integrity and quality of the extracted RNA was determined by measuring the absorbance at 260/280 nm (A260/A280) using an Agilent BioAnalyzer 2100 (Agilent, Shanghai, China). The short read alignment tool Bowtie2 (2.2.8) (Langmead & Salzberg, 2012) was used to map the reads to a ribosomal RNA (rRNA) database. The rRNA-mapped reads were then removed, and all other RNA reads were retained.

Library construction

Total RNA was randomly cleaved into short fragments. To synthesize first-strand cDNA, the RNA fragments were used as templates in a reaction with random hexamers, buffered dNTPs (with dUTP instead of dTTP) and RNase H. Second-strand cDNA was synthesized using DNA polymerase I. The cDNA was purified using a QIAquick PCR kit (Qiagen) and eluted with EB buffer (10 mM Tris-HCl, pH 8.0). The sticky ends were repaired, and an adaptor sequence was ligated to the A base at the 3′ end of each cDNA molecule; the second strand was subsequently degraded with uracil-N-glycosylase (UNG). The ligated products were reverse transcribed by PCR amplification, and the products were purified to generate a cDNA library. After generation, the cDNA library was sequenced using an Illumina HiSeq™ 4000 from Gene Denovo Biotechnology Co. (Guangzhou, China).

Alignment with reference genome

To obtain high-quality clean reads, all raw reads were further filtered by removing the reads containing adapters, reads that were less than 50 nt, reads containing poly N, and low-quality reads from raw data. Then, the alignment tool Bowtie2 (2.2.8) was used to map the high-quality clean reads to an rRNA database (mismatch number: 0), and the mapped rRNA reads were removed (Langmead & Salzberg, 2012). The comparison software TopHat2 (2.1.1) was used to align the remaining reads to the goat reference genome (http://www.ncbi.nlm.nih.gov/genome/?term=goat) (Kim et al., 2013). The TopHat2-based alignment results were used to reconstruct the transcripts using Cufflinks. Based on the positions of the assembled transcripts on the reference genome and the screening criteria of a transcript length ≥200 bp and a number of exons ≥2, known transcripts and new transcripts were obtained.

Screening of differentially expressed genes (DEGs)

Transcript abundances were quantified by RSEM software (Li & Dewey, 2011). FPKM (Fragments Per Kilobase per Million) was used as the standard for screening differentially expressed transcripts. EdgeR v3.26.8 was used for differential mRNA expression analysis among the 3 age groups (Robinson, McCarthy & Smyth, 2010). To adjust the P-values, the Benjamini–Hochberg method for controlling the FDR was used (Guo & Bhaskara Rao, 2008). Genes with FPKM >1, false discovery rate (FDR) < 0.05 and a |log2(fold change)| > 2 were considered significantly differentially expressed (Trapnell et al., 2010). Differentially expressed coding RNAs were then subjected to enrichment analysis of GO functions and KEGG pathways (Kanehisa & Goto, 2000).

Based on the expression levels of the transcripts, the relationships between samples and transcripts were hierarchically clustered, and heat maps were used to present the clustering results. The expression scale of 9 groups of samples was Z-Score standardized using the R language scale function, and the heat map was constructed using ggplot two packages in R (version 3.2).

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses

GO and KEGG pathway enrichment analyses were performed using the R package clusterProfiler. The DEGs were mapped to the GO database (http://www.geneontology.org/), and the obtained genes for each term were counted to obtain a list of transcripts enriched for specific GO terms.

The DEGs were mapped to the KEGG database (https://www.genome.jp/kegg/) for KEGG pathway analysis using the same method as that used for the GO analysis. In the GO and KEGG pathway analyses, the FDR after p value calibration was used as the screening standard, and pathways with FDR values ≤0.05 were considered significantly enriched pathways.

Protein–protein interaction (PPI) network analysis and gene immune function network construction

Based on information from STRING v.10.0, a PPI network was constructed for the protein-coding genes (Szklarczyk et al., 2015). Then, credible interactions (with combined scores ≥ 0.4) were accepted for further network analysis. Cytoscape (version 3.5.1) was also used (Shannon et al., 2003). The two resulting networks were merged to obtain the final interaction network. The degree-sorted layout was used for interaction degree analysis (in which all nodes were sorted by their interaction degree values).

The gene immune function network was constructed using the R Package ‘enrichplot’ (https://github.com/YuLab-SMU/enrichplot) with 18 genes and eight GO terms selected from the PPI network and GO/KEGG analyses.

Verification using real-time quantitative PCR (RT-qPCR)

A total of 18 DEGs were randomly screened from the transcriptome sequencing results for RT-qPCR analysis. The RT-qPCR primers used for validation were designed according to the sequences of selected DEGs using Primer Premier 5.0 (Premier, Canada), which are provided in Table 1. cDNA was synthesized by using the Primerscript RT reagent Kit (TaKaRa) with 1 µg of total RNA. RT-qPCR was performed using SYBR Premix Ex Taq II (TaKaRa, DRR081A) and conducted using a Stratagene Mx3000 real-time system (USA). A reaction volume of 15 µL with 0.3 µL of forwarding primer (0.4 mmol each), 0.3 µL of reverse primer (0.4 mmol each), 7.5 µL of 2X SYBR Premix Ex Taq II reaction buffer, 0.3 mL of 50X ROX Reference Dye II, 1.5 µL of cDNA (10 pg/mL, 1 mg/mL), and 5.1 µL of RNasE-free dH₂O was used. The reaction mixtures were incubated in a 96-well plate in an Mx3000p SYBR Green real-time quantitative PCR analyser (Stratagene, USA) at 95 °C for 30 s, followed by 40 cycles of 95 °C for 5 s and 60 °C for 30 s. All reactions were performed in triplicate. ACTB was used as a reference gene. After amplification, the 2^−ΔΔCT method was used to calculate gene expression levels, the significance of gene expression between different groups was tested by Tukey’s HSD, and gene expression histograms between different groups were performed by GraphPad software.

Summary of sequencing data

We constructed nine submandibular gland cDNA libraries for Lubo goats of three ages. A total of 862,995,872 raw reads were obtained from sequencing. After screening, 853,786,462 clean reads were obtained, accounting for more than 98% of the total number of reads (Table S1). The alignment rates for the nine libraries with the goat genome were above 80%; specifically, the mapping rates for A1-L, A2-L, A3-L, B3-L, B4-L, B5-L, C2-L, C3-L and C5-L were 91.02%, 93.13%, 91.66%, 88.54%, 87.30%, 88.00%, 83.76%, 84.49% and 87.57%, respectively (Table S2). All data have been uploaded to the Gene Expression Omnibus (GEO) database under the accession number GSE144368.

The numbers of known and new mRNAs that were identified and predicted in each library are provided in Table S6. As shown in Table 2, a total of 42,686 mRNAs were detected in all libraries: 41,914 in the 1-month-old goat library (group A), 39,070 in the 12-month-old goat library (group B), and 37,923 in the 24-month-old goat library (group C). Overall, 35,451 known mRNAs were detected in all libraries, accounting for 83.05% of the total number of mRNAs. The numbers of known mRNAs in group A, group B and group C were 32,875 (77.02% of the total), 30,577 (71.63% of the total), and 29,634 (69.42% of the total), respectively. A total of 10,053 new mRNAs were found in all libraries, including 9,039 in group A, 8,493 in group B, and 8,289 in group C (Table 2). The C5-L library had the lowest numbers of known mRNAs, new reference mRNAs, and total mRNAs (23,946, 6,080, and 30,026, respectively). The A3-L library had the most known mRNAs, new mRNAs, and total mRNAs (28,581, 7,549, and 36,130, respectively).

Identification and analysis of DEGs

The DEG analysis revealed 2,706 and 2,525 differentially expressed mRNAs between group A and group B and between group A and group C, respectively, but only 52 differentially expressed mRNAs between group B and group C. These results may indicate that mRNA expression levels in the submandibular salivary glands of 12-month-old goats and 24-month-old goats are highly similar. Numerous mRNAs were downregulated in group B and group C compared with group A (Fig. 1A). The clustering diagram of the differential expression patterns between the groups (Fig. 1B) shows that there were many DEGs between group A and group B and between group A and group C, while there were no significant DEGs between group B and group C.

GO and KEGG pathway analyses

GO and KEGG enrichment analyses were performed on DEGs to further confirm their involvement in the regulation of biological processes. The GO analysis results are shown in Figs. 2 and 3. In the DEG enrichment analysis between group A and group B, 3,673 GO terms were obtained, among which 427 GO terms were significantly enriched, including 46, 16 and 365 GO terms in the cellular component, molecular function and biological process categories, respectively (FDR < 0.05) (Fig. 2A; Table S7). Between group A and group C, 3,569 GO terms were obtained, among which 362 GO terms were significantly enriched, including 41, 21, and 300 GO terms in the cellular component, molecular function and biological process categories, respectively (FDR < 0.05) (Fig. 2B; Table S8). Between group B and group C, 516 GO terms were obtained, but there were no significantly enriched annotations (FDR < 0.05) (Fig. 2C). Notably, the top four significantly enriched terms for the DEGs between group A and group B and the DEGs between group A and group C were the same; these terms were the immune system process (GO:00023761, FDR = 6.26E−38 and 1.44E−27, respectively), leukocyte activation (GO:0045321, FDR = 5.44E−33 and 2.78E−25, respectively), lymphocyte activation (GO:0046649, FDR = 9.28E−32 and 5.51E−24, respectively) and cell activation (GO:0001775, FDR = 2.57E−31 and 1.16E−23, respectively) terms. At the same time, most of the significantly enriched GO terms of group A vs. B and A vs. C in Biological Process were immune-related terms (Fig. 3).

KEGG pathway analysis helped to further elucidate the biological functions of the transcripts (Fig. 4). Analysis of the differentially expressed mRNAs among goats at 3 developmental stages revealed 769 DEGs between group A and group B, of which 456 DEGs were significantly enriched in 47 KEGG pathways (FDR < 0.05). Among all significantly enriched pathways, the systemic lupus erythaematosus (FDR = 1.90E−33) and T cell receptor (TCR) signalling pathways (FDR = 1.14E−17) were the most highly enriched (Fig. 4A). A total of 683 DEGs were identified between group A and group C, of which 390 DEGs were significantly enriched in 31 KEGG pathways (FDR < 0.05) (Fig. 4B). The pathways with the highest degrees of enrichment were again the systemic lupus erythaematosus (FDR = 3.97E−27) and TCR signalling pathways (FDR = 3.09E−19). A total of 11 DEGs were identified between group B and group C that participated in 23 signalling pathways (Fig. 4C). However, only one gene was enriched per pathway, and no pathways reached a level of significant enrichment.

Notably, more than half of all significantly enriched KEGG pathways between groups A and B and between groups A and C were related to immunity. Among the 47 significantly enriched KEGG pathways between groups A and B, 35 pathways belonged to the immune system and immune disease classes, and 12 other pathways belonged to the cancer and infectious disease classes. Similarly, among the 31 significantly enriched KEGG pathways between groups A and C, 20 pathways belonged to the immune system and immune disease classes, and 11 other pathways belonged to the cancer and infectious disease classes. These results reveal that there are transcription level differences in the submandibular glands between 1-month-old kids and adult goats. The DEGs between groups A and B and between groups A and C were all enriched in 10 identical immune system pathways. However, the Q values were different, and the upregulated genes were not exactly the same in both comparisons. The immunoglobulins involved in these 10 pathways were upregulated, while almost all the other genes were downregulated. The Q values and upregulation levels of the DEGs in the 10 immune system pathways between groups A and B and between groups A and C are shown in Table 3.

Among the DEGs, genes encoding 8 antibacterial peptides, including the BAC7.5 protein precursor, MAP28, BAC5, MAP34-B, CATHL1A, LOC102171106, LOC102169231 and LOC102170003, were relatively highly expressed in the transcriptome of group A (FPKM values: 162.88, 50.45, 46.60, 43.32, 22.67, 35.32333, 29.97333, and 123.5833, respectively). However, there was almost no expression of the above 8 peptide-encoding mRNAs in group B (FPKM values: all 0.001) and in group C (FPKM values: 0.001, 0.001, 0.026667, 0, 0.001, 0.001, 0.001, and 0.001, respectively). The above results are shown in Fig. 5A.

Screening of key regulatory genes

To further identify the key regulatory genes that caused the differential enrichment of salivary gland immune-related pathways between the different age groups, a PPI network analysis was performed. According to the KEGG pathway enrichment results, a total of 27 significantly enriched pathways overlapped between the two comparisons (A vs. B and A vs. C); these pathways contained 290 genes that overlapped between both comparisons.

All 290 genes were used for PPI network construction. A total of 2,319 interactions among 202 genes were obtained (Fig. 6). All genes in the network were ranked by interaction degree, and the top 10 were selected as potential key regulatory genes: PTPRC, CD28, SELL, LCP2, MYC, LCK, ZAP70, ITGB2, SYK and CCR7. We then performed KEGG enrichment analysis for the 10 key regulatory genes, and two signalling pathways were determined to be significantly enriched: the TCR signalling pathway (FDR = 5.0E−05) and the NF- κβ signalling pathway (FDR = 2.2E−02). Since these two pathways were also the top enriched pathways in the DEG enrichment analysis, we identified them as key pathways. Gene immune function network analysis was performed, and 8 GO terms were determined to be significantly enriched: T cell receptor signalling pathway (FDR = 1.65E−08), neutrophil mediated immunity (FDR =7.37E−05), B cell mediated immunity (FDR = 1.84E−07), regulation of alpha-beta T cell activation (FDR = 3.46E−06), positive regulation of T cell proliferation (FDR = 3.65E−06), regulation of leukocyte differentiation (FDR = 4.12E−06), positive regulation of antigen receptor-mediated signalling pathway (FDR = 5.92E−06), positive regulation of lymphocyte proliferation (FDR = 1.18E−05) (Fig. 6B).

Validation of DEGs using RT-qPCR

To verify the accuracy of the sequencing and analysis results, 18 differentially expressed mRNAs were selected for RT-qPCR validation (Fig. 7). The RT-qPCR validation results for most selected genes were consistent with the RNA-seq results; only XM_018055547.1 of the DNA (cytosinE-5)-methyltransferase 3A (DNMT3a) gene did not show differential expression in the RT-qPCR analysis. In addition, the forkhead box protein M1 (FOXM1) and interleukin 1 receptor type 2 (IL1R2) genes, which were not expressed in the salivary glands of the goats in group B and group C based on the sequencing data, were detected by RT-qPCR in these two groups.