INTRODUCTION
The World Health Organization (WHO) reports that since 1975, the number of obese people in the world has grown nearly three-fold, and currently, nearly 20 million adults are overweight or obese [
1]. Obesity not only affects one’s physical appearance but also causes various diseases [
2] that burden the economy and the medical and health infrastructure. At present, popular methods of weight loss include dieting, targeted drug therapy, and surgical liposuction. Weight loss based on long-term fasting has significant effects, but it is not healthy [
3]. Short-term starvation stress might be beneficial for cardiometabolic health by decreasing insulin resistance, blood pressure, and oxidative stress [
4]. Long-term starvation stress can lead to a deficit of effective nutrients, increase levels of corticosterone, the stress hormone, and inhibit cardiovascular function [
5].
The intestine is the main site of the digestion, absorption, and metabolism of nutrients in the body. The host intestinal mucosa directly or indirectly participates in metabolism, and its function is closely related to metabolic diseases such as obesity, diabetes, and allergies; moreover, it plays an active role in immune regulation [
6]. Studies have reported that starvation can cause changes in the intestinal structure, damage the intestinal barrier, decrease the villus height-to-crypt depth (V/C) ratio, and dysregulate the intestinal flora, which eventually leads to metabolic disorders [
7].
Long non-coding RNAs (lncRNAs) are a class of non-coding RNAs longer than 200 nucleotides that are widely present in the nucleus and cytoplasm [
8]. LncRNAs participate in the regulation of various biological processes (BPs), such as gene expression, transcriptional activation, cell cycle regulation, and disease occurrence [
9]. Recent studies have surmised that an increasing number of lncRNAs play important roles in maintaining intestinal epithelial homeostasis, promoting the proliferation and apoptosis of intestinal epithelial cells, and regulating the intestinal barrier [
10–
12]. For example, lncRNA H19 promotes the expression of miR-675-5p and miR-675-3p, which downregulate the levels of tight junction proteins 1 (TJP1, also called ZO-1) and E-cadherin, destroying the intestinal epithelial barrier [
10]. In mice with starvation-induced mucosal atrophy, lncRNA uc.173 promotes the proliferation of intestinal epithelial cells by reducing the expression of miR-195 [
11]. Transcriptome sequencing performed on ulcerative colitis tissue revealed lncRNA BC012900 overexpression. Furthermore, the overexpression of this lncRNA promotes intestinal epithelial cell apoptosis [
12].
Owing to the availability of genetic information, easy genetic modification, and low cost of rearing, mice are the preferred
in vivo models for studying human diseases [
5,
13]. However, significant differences exist in metabolic and physiological characteristics between humans and mice, which prevents researchers from applying findings from murine studies to humans for metabolism-related disease prevention and intervention strategies [
14]. In contrast, pigs are very similar to humans in terms of metabolic characteristics and organ development; therefore, pigs are an ideal animal model for human energy metabolism research [
15]. Thus, in this study, we subjected weaned piglets to different periods of starvation and used their ileum tissue for RNA sequencing analysis. In this study, we aimed to identify lncRNAs that have important functions in starvation stress, analyze their functions, and discover potential molecular targets for alleviating starvation stress, to provide a theoretical reference for subsequent in-depth research.
MATERIALS AND METHODS
Ethics statement
This research strictly adhered to the principles of animal use prescribed by the China Laboratory Animal Science Association. The study was approved by the Animal Ethics Committee of Shanxi Agricultural University (Taigu, China). The approval number for the Ethics Committee agreement was SXAU-EAW-2018P002005. The animals were humanely sacrificed as necessary to ameliorate suffering.
Sample collection
Yorkshire piglets were bred at the animal experiment station of Shanxi Agricultural University in accordance with the National Research Council (NRC) [
16]. The test animals were from three litters of half-sibling piglets of the same birth age. On the weaning day, three boars with similar body weights (25 d-old, 5.96±0.15 kg) were selected from each litter, and nine weaned Yorkshire piglets were divided into three groups according to the block design principle. The piglets were randomly divided into a long-term starvation stress group (72 h, abbreviated as IHT-72), short-term starvation stress group (48 h, abbreviated as IHT-48), and a normal group (control group, abbreviated as ICT). The piglets were sacrificed, and ileum samples, which connected to the ileocecal ligament, were collected. The samples were washed with phosphate-buffered saline and immediately frozen in liquid nitrogen until RNA extraction.
Morphological observation of ileum villi
The ileum samples (ICT, IHT-48, and IHT-72) were analyzed following hematoxylin and eosin (H&E) staining (Bosterbio, Wuhan, Hubei, China). The tissue samples were fixed in 4% paraformaldehyde for 24 h and processed using routine histological methods. Subsequently, 7 μm-thick sections were cut using a Leica RM2265 (Leica, Wetzlar, Germany) and stained. Three images were acquired (Olympus, Tokyo, Japan) for each section, and three sections were selected for each piglet. The lengths of intestinal villi were calculated and counted using SPSS ver.22.0 (IBM Corp., Armonk, NY, USA).
RNA extraction, library preparation, and high-throughput sequencing
Total RNA was extracted from the nine libraries using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. RNA integrity was detected by performing 1% agarose gel electrophoresis, and the purity and concentration were determined using an ND-2000 nucleic acid-protein analyzer (Nanodrop Technologies, Winooski, VT, USA). The ratio of absorbances at 260 nm and 280 nm was determined, with values ranging from 1.8 to 2.0, indicating an RNA concentration of approximately 1,000 ng/μL. After total RNA was extracted, the rRNAs were removed to retain mRNAs and ncRNAs. The cDNA library was generated according to the instructions of the TruSeqRapid Duo cBot Sample Loading kit (Illumina, San Diego, CA, USA) and sequenced using an Illumina HiSeqTM 4000 by Gene Denovo Biotechnology Co. (Guangzhou, China).
Transcriptome assembly
To obtain high-quality clean reads, reads were filtered using fastp [
17] (version 0.18.0). Data were cleaned by removing reads containing adapters, those containing more than 10% poly(N), and low-quality reads (containing more than 50% low-quality [Q-value≤20] bases) from the raw data. The short-read sequence alignment tool Bowtie2 [
18] (version 2.2.8) was used to map the read sequence to the ribosomal RNA (rRNA) database. An index of the reference genome (Sscrofa11.1 [GCF_000003025.6]) was built, and paired-end clean reads were mapped to the reference genome using HISAT2 (version 2.1.0). The rRNA database uses the Nucleotide Sequence Database. Transcripts were reconstructed using Stringtie software (version 1.3.4), which together with HISAT2 allowed us to identify new genes and new splice variants of known genes.
Identification and annotation of novel lncRNAs
The lncRNAs were screened according to a previously described process (
SI Appendix, Figure S1). From the remaining transcripts that overlapped (>1 bp) with pig protein-coding genes, transcripts <200 bp and single-exon transcripts were removed. The coding potential calculator (CPC) (version 0.9-r2) and Coding-Non-Coding-Index (CNCI) (version 2) tools were used to assess the coding potential of the remaining transcripts, and transcripts with a CPC score >0 and CNCI score >0 were removed. The intersection of both non-protein-coding potential results was considered lncRNAs.
Differential expression analysis of lncRNAs
Differential expression of RNA and lncRNAs among different groups was analyzed using DESeq2 software. In addition, different samples in the same group were compared with edgeR. Transcripts with a false discovery rate (FDR) less than 0.05 and an absolute fold change ≥2 were considered differentially expressed. Differentially expressed coding RNAs were subjected to enrichment analysis based on gene ontology (GO) functions and Kyoto encyclopedia of genes and genomes (KEGG) pathway analysis.
LncRNA–mRNA association analysis
lncRNAs regulate mRNA expression levels via cis and trans interactions. Target genes within 10 kb of differentially expressed lncRNAs (DElncRNA) were selected to explore cis-regulation, and DElncRNA–mRNA target gene pairs with a Pearson correlation coefficient greater than 0.99 were selected to explore co-expression regulation. DElncRNA–mRNAs were visualized using Cytoscape v3.8.2.
Gene expression detection by real time quantitative polymerase chain reaction
The extracted total RNA was diluted to 500 ng/μL and reverse-transcribed according to the instructions of the PrimeScript RT reagent kit with gDNA Eraser (TaKaRa Bio, Dalian, China). real time quantitative polymerase chain reaction (RT-qPCR) was performed on a CFX96 Real-Time PCR detection system (Bio-Rad, Hercules, CA, USA) using a SYBR Green RT-qPCR kit (TaKaRa Bio, China). 18S rRNA was used as an internal reference gene, and the sequences of primers used are listed in
Supplementary Table S1. Three technical replicates were performed for each sample, along with a non-template negative control and a negative control without reverse transcriptase. The RT-qPCR cycling parameters were as follows: pre-denaturation at 95°C for 30 s, followed by 40 cycles of denaturation at 95°C for 30 s and annealing/extension at 60°C for 30 s. Melting curve analysis comprised a denaturation step at 95°C for 30 s, 60°C for 1 min, and 95°C for 30 s.
Statistical analyses
RT-qPCR results were used to calculate the relative expression levels of the target genes using the 2−ΔΔCt method. Statistical analyses were performed using SPSS ver.22.0 (IBM Corp., USA). Differences in gene expression among the different groups were identified using analysis of variance using Duncan’s test for significance. Pearson correlation analysis was used to determine the relationship between lncRNA expression and mRNA expression and to verify RNA-seq and RT-qPCR results. Graphpad version 8.0 (GraphPad Software, San Diego, CA, USA) was used to draw histograms. All tests were two-tailed, and a p-value <0.05 was considered a significant difference.
Data availability
The sequencing data obtained have been deposited in the Sequence Read Archive with the accession number PRJNA 770271.
DISCUSSION
Despite the continuous increase in food production in modern times, starvation continues in various parts of the world either due to poverty or geological disasters. In contrast, several individuals starve voluntarily; for example, more people are inclined to adopt improper ways to lose weight by fasting for long durations [
3]. The pig was used in this experiment as an ideal model to study human metabolic diseases [
15]. In this study, we analyzed the transcriptional changes in piglets subjected to starvation stress via high-throughput sequencing and revealed the damage caused by starvation stress to the body from the perspective of lncRNA expression.
Genome-wide studies have revealed that lncRNAs play an important role in pigs [
22]. However, few studies have been conducted on the expression of starvation stress-related lncRNA expression in pigs. In this study, we established a piglet model of starvation stress and found that the intestine is one of the key tissues damaged by starvation stress. Illumina HiSeqTM 4000 high-throughput sequencing was used to construct differential expression profiles of starvation stress-related lncRNAs in the piglet ileum. By testing the relationships between the samples, we found that at the transcriptome level, negligible differences were observed between the control group and the group starved for 48 h. However, compared with the ICT and IHT-48 groups, notable differences were observed in the group starved for 72 h. In our study, 72 h starvation inhibited glycolysis, gluconeogenesis, the TCA cycle, fructose and mannose metabolism, pyruvate metabolism, and vitamin digestion and absorption. This result suggests that long-term starvation leads to an insufficient energy supply, exhaustion of glycogen storage, and nutrient deficits. Furthermore, long-term starvation affects the body’s immunity and increases the risk of disease. During 72 h of starvation stress, mRNAs related to CAMs were overexpressed. The CAMs play an important role in the occurrence of diseases and are closely related to inflammation [
23]. Further, 72 h starvation could induce Th1/Th2 cell differentiation and lead to inflammatory bowel disease and malaria. However, it might reduce the onset of diabetes, which could be related to the insufficient energy supply. In summary, prolonged starvation can cause irreversible damage to the body and at the same time provides a theoretical basis for the “golden 72 hours” associated with earthquake rescue.
A co-expression sub-network containing 26 DElncRNAs and 72 DEmRNAs was reconstructed to reveal the mechanisms underlying starvation stress and its effect. LncRNA MSTRG.19894.13 was determined to be an important lncRNA in the constructed network.
C4BPA, dehydrogenase/reductase 4 (
DHRS4),
SDSL,
UGDH, and
XYLB were identified as the target genes of this lncRNA.
C4BPA controls the classic complement activation pathway and interacts with RelA, a member of the nuclear factor kappa-B (NF-κB) family [
24]. The expression of
C4BPA is regulated by stress [
24].
DHRS4 is a peroxisomal member of the short-chain dehydrogenase/reductase superfamily. Unlike that in its human counterpart, a mutation at Thr177 with the corresponding residue Asn stabilizes DHRS4 in pigs and prevents its cold-induced inactivation [
25]. SDSL, also known as SDH2, exhibits low serine dehydratase and threonine dehydratase activities [
26]. UGDH is a cytosolic enzyme that catalyzes the oxidation of UDP-glucose to UDP-glucuronic acid, thereby participating in the biosynthesis of glycosaminoglycans, such as hyaluronic acid, chondroitin sulfate, and heparan sulfate [
27]. In patients with lung adenocarcinoma, the presence of UGDH in the nucleus is correlated with poorly differentiated cells and larger tumors and could indicate overall reduced survival [
28].
XYLB encodes a xylulose kinase-like protein, and its forced overexpression is an effective strategy to improve xylose utilization and P(3HB) production in
Burkholderia sacchari [
29]. Our results showed that MSTRG. 19894.13 regulates the expression of
C4BPA,
DHRS4,
SDSL,
UGDH, and
XYLB, and therefore, it might play an important role in regulating functions such as complement activation, enzyme activity, biosynthesis, and energy utilization.
The co-expressed target genes of MSTRG.16726.3 were determined to be aldo-keto reductase family 1 member A1 (
AKR1A1),
C3,
CFH,
GST01, MYB proto-oncogene (
MYB), and polymeric immunoglobulin receptor (
PIGR). Studies have shown that AKR1A1 can synthesize ascorbic acid in rodents [
30] and is involved in the metabolism of gamma-hydroxybutyric acid in human liver cancer-derived HepG2 cells [
30]. In contrast, C3 plays a central role in activation of the complement system and can be used as a biomarker for insulin resistance and cardiometabolic diseases [
31]. CFH plays an important role in regulating complement activation and limits the effect of innate defense mechanisms against microbial infection [
32]. GSTO1 is involved in the metabolism of xenobiotics and carcinogens, such as with cutaneous malignant melanoma [
33]. Mutations in and the overexpression of
MYB were first discovered in leukemia cells and were recently discovered in solid cancers [
34].
MYB plays an important role in hematopoietic functions and tumorigenesis [
34]. The circ-XPO4 of the small extracellular vesicles present in milk promotes the expression of
PIGR in intestinal cells by inhibiting miR-221-5p, thereby increasing the level of intestinal SIgA [
35]. In conclusion, MSTRG.16726.3 might play an important role in complement activation, the metabolism of exogenous and carcinogenic substances, resistance against pathogenic microorganisms, and improvements in intestinal immunity.
The following nine genes are targeted by MSTRG.12176.1: C-C motif chemokine ligand 25 (
CCL25),
CDX2,
GALE, hepatocyte nuclear factor 4 gamma (
HNF4G),
PDXK,
PLCB3,
RAB27B, solute carrier family 2 member 5 (
SLC2A5), and solute carrier family 34 member 3 (
SLC34A3). CCL25, a member of the chemokine CC subfamily, has also been recognized as a thymus-expressed chemokine [
36]. CCL25 and its receptor C-C motif chemokine receptor 9 (CCR9) are overexpressed in cancers such as leukemia, ovarian cancer, breast cancer, prostate cancer, liver cancer, lung cancer, and melanoma [
36].
CDX2 is a major regulator of intestine-specific genes involved in cell growth and differentiation [
37]. Decreased expression of
CDX2 is associated with mucinous tumors, lymph node involvement, and high-grade tumors [
37].
GALE encodes UDP-galactose-4-epimerase, and mutations in
GALE result in epimerase-deficient galactosemia, also referred to as galactosemia type 3, a disease characterized by early onset cataracts, liver damage, deafness, and mental retardation [
38].
HNF4G is overexpressed in colorectal cancer (CRC) and promotes CRC cell proliferation via the PI3K/AKT pathway by targeting G protein subunit gamma 12 (GNG12) and protein tyrosine kinase 2 (PTK2) [
39].
PDXK encodes a pyridoxal kinase, which converts inactive B vitamins to the active cofactor pyridoxal 5′-phosphate (PLP) [
40]. Hereditary polyneuropathy with optic atrophy can occur due to a
PDXK variant leading to impaired vitamin B6 metabolism [
40].
PLCB3 is involved in innate immunity, and studies have shown that silencing
PLCB3 enhances the inflammatory signaling cascade of toll-like receptors [
41].
RAB27B is an RAS oncogenic family member gene that regulates extracellular vesicle production in cells infected with Kaposi’s sarcoma-associated herpesvirus to promote cell survival and persistent infection [
42]. The fructose transporter encoded by
SLC2A5 is required for intestinal fructose absorption [
43].
SLC2A5 expression is induced in the intestine and skeletal muscle of patients with type 2 diabetes and in certain cancers dependent on fructose metabolism [
43]. In renal brush border cells,
SLC34A3 is involved in transporting intracellular phosphate via sodium co-transport and contributes to the maintenance of inorganic phosphate concentrations in the kidneys [
44]. Mutations in
SLC34A3 can cause hereditary hypophosphatemic rickets with hypercalciuria [
44]. Therefore, by acting on multiple target genes, MSTRG.12176.1 can play an important role in the starvation stress in piglets and affect energy metabolism, nutrient absorption, and disease occurrence and development.
In this study, 11,792 lncRNAs and 63,682 mRNAs were identified through transcriptome sequencing of ileum tissue obtained from pigs subjected to starvation stress. These included 2,500 novel lncRNAs, 509 DElncRNAs, and 3,319 DEmRNAs. The target genes of DElncRNAs were determined to be mainly involved in metabolic pathways, cellular processes, immune system processes, digestive systems, and transport activities, among others. A co-expression network induced by starvation stress in pigs was constructed, and the key lncRNAs were analyzed. However, the specific mechanism by which these lncRNAs regulate metabolism is unknown. Further studies are required to elucidate the molecular mechanisms linking lncRNAs to the regulation of starvation stress in pigs.