Endothelial Cell Heterogeneity and Microglia Regulons Revealed by A Pig Cell Landscape at Single-Cell Level
Endothelial Cell Heterogeneity and Microglia Regulons Revealed by A Pig Cell Landscape at Single-Cell Level
Endothelial Cell Heterogeneity and Microglia Regulons Revealed by A Pig Cell Landscape at Single-Cell Level
https://doi.org/10.1038/s41467-022-31388-z OPEN
Wendi Wu3,11, Yuting Yuan12, Mingyi Pu3,11, Haoyu Wang3,4, Aiping Wu13,14, Lin Xie8, Ping Liu8, Fang Chen8,
Jacqueline Herold2, Joanna Kalucka 2,6,15, Max Karlsson 16, Xiuqing Zhang3,4, Rikke Bek Helmig17,
Linn Fagerberg 16, Cecilia Lindskog 18, Fredrik Pontén 18, Mathias Uhlen 7,16, Lars Bolund1,2,
Niels Jessen 6, Hui Jiang8, Xun Xu 3, Huanming Yang3,19, Peter Carmeliet 2,9,10,20, Jan Mulder 7,
Dongsheng Chen3,13,14 ✉, Lin Lin 2,6 ✉ & Yonglun Luo 1,2,3,6,19 ✉
Pigs are valuable large animal models for biomedical and genetic research, but insights into
the tissue- and cell-type-specific transcriptome and heterogeneity remain limited. By lever-
aging single-cell RNA sequencing, we generate a multiple-organ single-cell transcriptomic
map containing over 200,000 pig cells from 20 tissues/organs. We comprehensively
characterize the heterogeneity of cells in tissues and identify 234 cell clusters, representing
58 major cell types. In-depth integrative analysis of endothelial cells reveals a high degree of
heterogeneity. We identify several functionally distinct endothelial cell phenotypes, including
an endothelial to mesenchymal transition subtype in adipose tissues. Intercellular commu-
nication analysis predicts tissue- and cell type-specific crosstalk between endothelial cells
and other cell types through the VEGF, PDGF, TGF-β, and BMP pathways. Regulon analysis of
single-cell transcriptome of microglia in pig and 12 other species further identifies MEF2C as
an evolutionally conserved regulon in the microglia. Our work describes the landscape of
single-cell transcriptomes within diverse pig organs and identifies the heterogeneity of
endothelial cells and evolutionally conserved regulon in microglia.
T
he domestic pig (Sus scrofa domesticus) is an important including nine brain regions, retina, kidney, heart, spleen, liver,
large animal for modeling both monogenic and complex and lung, were analyzed by single-nuclei RNA sequencing
human diseases such as cystic fibrosis1, atherosclerosis2, (snRNA-seq). Among the nine different brain regions, four
Huntington’s diseases3, and diabetes4–6. Additionally, it is regions (area postrema, cerebellum, subfornical organ, and the
explored and regarded as the most promising source for vascular organ of lamina terminalis (OVoLT)) were generated by
xenotransplantation7–9 largely due to the postulation of their high this study, and the other five brain regions (frontal lobe, hypo-
resemblance to humans in organ size, structure, anatomy, thalamus, occipital lobe, parietal lobe, and temporal lobe) were
genetics, and physiological functions10,11. Recently, clinical from a study we had previously published using snRNA-seq28.
transplantation of genetically modified pig kidneys and heart has Both scRNA-seq and snRNA-seq techniques were used for single-
successfully been achieved with patients12. Despite all these great cell transcriptome analysis, which had both technical strengths
promises and progress, several species-specific cellular and and limitations31. The scRNA-seq is performed using freshly
molecular differences between pigs and humans exist. For example, isolated single cells, thus capturing all transcripts in the cells, but
pluripotency progression and metabolic transition were found to be limited by the sample processing procedures. Single-cell suspen-
different using single-cell RNA-sequencing (scRNA-seq)13. To gain sions must be prepared from fresh tissues for scRNA-seq. For
better insights into the biomedical similarity, as well as differences snRNA-seq, tissues can be snap-frozen after sampling and used
between pig and human, and to advance the applications of pigs in for nuclei extraction, thus less limited by timing. We selected both
biomedical research, a comprehensive and body-wide investigation scRNA-seq and snRNA-seq for reasons of practicality and
of the domestic pig is needed at single-cell level. resource availability. To compare the two methods, four pig tis-
High-throughput scRNA-seq technology has greatly expanded sues (liver, retina, lung, spleen) were analyzed with both scRNA-
the ability to better understand the cell composition, interactions, seq and snRNA-seq. In total, the pig single-cell transcriptome
heterogeneity, and functions in highly organized and multi- atlas includes twenty pig tissues (Fig. 1a). All samples were dis-
cellular mammalian organs under physiological or pathological sociated into single cells or single-nuclei, followed by high-
conditions14,15. Pioneering work has been completed for most throughput scRNA-seq or/and snRNA-seq library generation,
model animals and humans16–22. Since the completion of the first deep sequencing, and data analyses (Fig. 1b). Complete infor-
porcine reference genome, technological breakthroughs in gene mation of all batches of samples was summarized in Supple-
editing and cloning allow precise manipulation of the pig genome mentary Data 1. After filtering low-quality cells
in living animals23–25 and applications of pigs in biomedical (nFeatures_RNA < 200 or % mitochondria transcripts > 30%) and
research have been greatly expanded. Previous scRNA-seq studies doublets (nFeatures_RNA > 5000), high quality single-cell tran-
in pigs were restricted to a few tissues or a single organ26–28. scriptome data were obtained from 222,526 cells, of which
Recently, the pig BodyMap atlas with the conventional bulk RNA- 133,492 and 89,034 were obtained with scRNA-seq and snRNA-
seq method has been reported29. Using similar bulk RNA-seq based seq, respectively.
transcriptome profiling, we have systematically profiled the protein- We first compared single-cell gene expression profiling of
coding transcriptome in 98 pig tissues (www.rnaatlas.org)30. tissues obtained with scRNA-seq and snRNA-seq. Quality control
In this study, we report the generation of a single-cell tran- (QC) comparison of scRNA-seq and snRNA-seq demonstrated
scriptome atlas of 222,526 cells across twenty tissues of the that scRNA-seq captures more transcripts per cell and obtains
domestic pig by using scRNA-seq and single-nuclei RNA- higher percentage of mitochondria- and ribosome genes. In
sequencing (snRNA-seq). The mapping results are available contrast, snRNA-seq captures a higher percentage of protein-
through the pig single-cell atlas database for comparative analyses coding transcripts and transcription factor-encoding gene
and data exploration (https://dreamapp.biomed.au.dk/pigatlas/). transcripts, largely due to the absence of mitochondria and
In total, 58 cell types were identified, which contribute to tissue- ribosome RNAs. Despite that, when evaluating cell cycle gene
specific and shared functions between the tissues. We identified expression, the scores of S and G2M stages are similar in tissues
tissue-specific cell types, as well as common cell types shared analyzed with scRNA-seq and snRNA-seq (Supplementary
across different tissues. Commonly shared cell types also exhibit Data 2). Next, cells from the spleen, liver, lung, and retina
tissue-specific expression patterns and functions. One such cell captured by snRNA-seq and scRNA-seq were clustered and
type is vascular endothelial cells (ECs). Further analysis of ECs visualized by t-SNE plots (Supplementary Fig. S1). Cell clusters
probed rare ECs types supporting the notion of endothelial to were annotated based on the expression of canonical cell-type
mesenchymal transition (EndMT), which we further validated by markers (Supplementary Fig. S1a–d). The number captured cells
scRNA-seq of cultured ECs and induced EndMT by transforming for each cell type in the four shared tissues analyzed by scRNA-
growth factor-beta 2 (TGF-β2) treatment. These single-cell seq and snRNA-seq are highly variable (Supplementary
transcriptome data allow us to gain insights into the similarities Fig. S1e–h), which can be explained by sampling and tissue
and differences in biomedical and cellular functions between pigs processing biases, as the scRNA-seq and snRNA-seq experiments
and humans. We also performed a pan-species regulon com- were carried out in two different laboratories. Despite that, our
parison analysis covering thirteen different species and identified results showed that common cell types captured by scRNA-seq
MEF2C as the most conserved regulon for microglia evolution and snRNA-seq exhibit similar gene expression profiles and cell-
over 300 million years. type-specific markers. Most importantly, correlation analysis
between shared cell types further showed that the transcriptome
of the same cell types identified scRNA-seq and snRNA-seq are
Results well correlated (Supplementary Fig. S1i–l). Thus, we combined
Single-cell and single-nuclei RNA sequencing of four pig tis- the scRNA-seq and snRNA-seq data, batch-corrected, and
sues highly correlates in common cell types. To construct the performed an integrative analysis with Seurat32.
first body-wide single-cell transcriptome atlas of pigs, 20 pig
tissues/organs were analyzed using single-cell and/or single-nuclei
sequencing. Nine pig tissues, including visceral and subcutaneous A pig single-cell transcriptome atlas. To identify the cell types
adipose tissues, spleen, intestine, liver, lung, peripheral blood based on their single-cell transcription profile, we first performed
mononuclear cells (PBMCs), whole brain, and retina, were ana- graph-based clustering and visualized all the cells from these
lyzed by single-cell RNA sequencing (scRNA-seq). Fifteen tissues, twenty tissues using t-distributed stochastic neighbor embedding
a b
Lung c,n Frontal lobe n (1) Single-cell/nucleus (2) Single-cell/nucleus (3) Library construction
c
Liver c,n PBMC Hypothalamus n preparation barcoding
Occipital lobe n
Parietal lobe n GEMs
c
Brain Temporal lobe n
Intestine c
Sus scrofa Area postrema n
domesticus (4) Sequencing
Cerebellumn
Subfornical organn (5) Bioinformatic analysis
c,n
Spleen c,n
OVoLT n
Retina
c
Subcutaneous adipose
n Adipose
Heart Visceral adipose
c
Cell cluster EC heterogeneity Cell communication Evolutionary analysis
Kidney n
c 1 14
3
d
8
7
Spleen Lung
2
20 21791(9.79%) 25725(11.56%)
12 1 Adipose_S
2 Adipose_V
1
16 3 Area postrema Adipose_V
11
6 6 11 4 Brain 16417(7.38%) Liver
20 18 5 Cerebellum 32888(14.78%)
13 9
1 6 Frontal lobe
2 5 2 4 7 Heart
Subfornical Organ
12 3
19 5 Brain 1527(0.69%)
10 8 Hypothalamus
13 9 Intestine
15794(7.10%) Temporal Lobe
12 16 18
10 3145(1.41%)
18 10 Kidney
Heart
4 20 7 11 Liver
11
3220(1.45%)
5 6 18
12 Lung Adipose_S Kidney
11
11
13 Occipital lobe 12814(5.76%) 4255(1.91%)
1
12
8 14 OVoLT Area Postrema
2 13 10 15 6 5109(2.30%)
14
15 Parietal lobe Retina
3 Parietal Lobe
16 PBMC 12164(5.47%) 6162(2.77%)
9 8 20 18 19 17 Retina
Intestine Occipital Lobe
11 18 Spleen
10518(4.73%) 6829(3.07%)
9
tSNE2
e 21 46 9
38
f
55 38 19
56 41 log(n)+1
7 34 49 57
12 56
37 39 6 36 18
14 33 21 10
38
42 2 4 46 35 56
28 11 11
19 5 11 1010
48 21 47
34 31 9 9 9
37 47
56
55 25 55 56 27 7
18
38 3
26 2 58 47 6
22 15 20 32 19
7 21 33 51 11 47 5
20 33 10 56 18
38 57 46 56
Intersection size
23 4 4
33 18 12 45 1
49
30 24 56 51 9 23 29
47 29 38 28 28 3 3 3
26 33 23 49
1 42 47 33
4 56 21 2
16 58 45
55 3250 42 24 47 38 2 2 2 2 2 2 2
25 20 47
51 47 30
46 38 31
24 55 19 52
21 44
1
5 33 46 24
30 56 31 27 54
Oligodendrocyte progenitor cells
38 40
Distal convoluted tubule cells
38 45
CD8+ Naive T cells
Mesenchymal cells
Excitatory neurons
49
Inhibitory neurons
Oligodendrocytes
Beige adipocytes
Megakaryocytes
56 31 21
Endothelial cells
Cardiomyocytes
Memory B cells
Cholangiocytes
55 25
Erythroid cells
Dendritic cells
Macrophages
53
Purkinje cells
Immune cells
Granule cells
Naive B cells
Ciliated cells
Plasma cells
Cycling cells
Hepatocytes
Kupffer cells
Bipolar cells
Enterocytes
Neutrophils
Fibroblasts
Monocytes
Astrocytes
Clara cells
54
Podocytes
Muller glia
Mast cells
Basophils
AT1&AT2
Microglia
NK cells
Others
B cells
T cells
tSNE1
AT1
AT2
0
1 Astrocytes 21 Endothelial cells 41 Naive B cells
Adipose_S
2 AT1 22 Enterocytes 42 Neural progenitor cells
Adipose_V
3 AT1&AT2 23 Erythroid cells 43 Neutrophils
Area postrema
4 AT2 24 Excitatory neurons 44 Newly formed oligodendrocytes
Brain
5 B cells 25 Fibroblasts 45 NK cells
Cerebellum
6 Basophils 26 Granule cells 46 Oligodendrocyte progenitor cells
Frontal lobe
Heart
7 Beige adipocytes 27 Hepatic stellate cells 47 Oligodendrocytes
Hypothalamus
8 Bipolar cells 28 Hepatocytes 48 Others
Intestine
9 Cardiomyocytes 29 Immune cells 49 Plasma cells
Kidney
10 CD4+ Naive T cells 30 Inhibitory neurons 50 Podocytes
Liver
11 CD8+ Cytotoxic T cells 31 Kupffer cells 51 Proximal tubule cells
Lung
12 CD8+ Naive T cells 32 Loop of Henle’s cells 52 Purkinje cells
Occipital lobe
13 Cholangiocytes 33 Macrophages 53 Retinal ganglion cells
OVoLT
14 Cilated cells 34 Mast cells 54 Rod photoreceptor cells
Parietal lobe
15 Clara cells 35 Megakaryocytes 55 Smooth muscle cells
PBMC
16 Collecting duct cells 36 Memory B cells 56 T cells
Retina
17 Cone photoreceptor cells 37 Mesenchymal cells 57 Transient amplifying cells
Spleen
18 Cycling cells 38 Microglia 58 Ureteric epithelial cells
Subfornical organ
19 Dendritic cells 39 Monocytes Temporal lobe
20 Distal convoluted tubule cells 40 Muller glia
100 0
Celltype proportions
(t-SNE) at tissue levels (Fig. 1c). On average, the number of corroborating the general transcriptome regulation of tissue
single-cell/nuclei transcriptomes obtained from each tissue ranges origin33. We next performed cell-type clustering and annotated
from 1527 cells (subfornical organ) to 32,888 cells (liver), each cluster of cells according to the expression of canonical cell-
representing 0.69% and 14.78% of the total cells, respectively type markers for each tissue separately (Supplementary Data 3).
(Fig. 1d). Most cells from each tissue are clustered separately thus In total, we identified 234 cell clusters corresponding to 58 major
Fig. 1 A single-cell transcriptome atlas of 20 pig tissues. a Schematic diagram of organs/tissues. Superscripts “c” and “n” represent the tissue analyzed
by scRNA-seq and snRNA-seq respectively. OVoLT, vascular organ of lamina terminalis. b Schematic diagram of cDNA libraries generation and
downstream bioinformatics analyses. The scRNA-seq and snRNA-seq were constructed independently, followed by high-throughput sequencing, and
downstream bioinformatic analyses. c t-SNE visualization of all single cells in the 20 tissues. Cells are color-coded according to the tissue origin. d Pie chart
showing the number of cells and proportion of cells from each tissue after filtering low-quality cells and doublets. e t-SNE visualization of all annotated
major cell types from the 20 tissues. Cells are color-coded according to cell types. f Bar graph and intersect plots showing the presence of the 58 cell types
across the 20 tissues. Source data are provided as a Source Data file. Schematic diagrams in a and b were created with BioRender.com.
cell types with significantly enriched markers (Fig. 1e, Supple- negative regulation of neuron differentiation are enriched in bipolar
mentary Fig. S2, Supplementary Data 4). The cell types in each cells. In addition, pathways such as axonogenesis, gliogenesis, and
organ were visualized by using t-SNE, which revealed a great negative regulation of neurogenesis, regulation of neuron projection
diversity of cell types within each tissue ranging from six cell development are enriched in Müller glia. Retinal ganglion cells share
types (subfornical organ) to sixteen cell types (lung) per tissue/ many enriched pathways with Müller glia. One specific pathway
organ (Fig. 1f, Supplementary Fig. S2, Supplementary Data 3). enriched in retinal ganglion cells (RGC) is neuron recognition,
With the transcriptional profiles from these large amounts of which is in good agreement with RGC’s function in pattern
cells, we comprehensively characterized the cell types in all tis- recognition and visual processing40,41. For rod photoreceptor cells,
sues. For example, beige adipocytes (CIDEA+, TBX1+) were only pathways such as photoreceptor cell differentiation, photoreceptor
identified in adipose tissues. Cardiomyocytes (ACTN2+, cell cilium, photoreceptor outer segment, detection of light stimulus,
MYH7+, FHL2+, and TNNT2) were obtained from the heart. and phototransduction are enriched (Fig. 2c, Supplementary Data 5).
Enterocytes (CHP2+, FABP1+, and FABP6+) are only detected Subsequently, we validated the canonical cell markers in the retina
in the intestine. Cholangiocytes (MMP7+, SPP1+, and ONE- cell types by protein staining. Rhodopsin RHO, which is essential for
CUT1+), hepatic stellate cells (ACAT2+, COL1A1+, RBP1+, and vision42, was highly expressed in the rod photoreceptor cells.
COLEC11+), hepatocytes (GHR+, HAMP+, HSD11B1+, and Arrestin 3 (ARR3), which is a non-visual arrestin and binds
RPP25L+), and Kupffer cells (CD163+and VSIG4+) were iden- phosphorylated G protein-coupled receptors, was highly expressed
tified in the liver. In the lung, the AT1 (AGER+, AQP5+, in the cone photoreceptor cells. The G Protein Subunit Gamma 13
CLIC5+, and SCGB1A1+) and AT2 (SFTPB+, SFTPD+, and (GNG13) was detected strongly expressed in bipolar cells, in line
ABCA3+) were identified. Immune cell types such as T cells, B with the previous findings by scRNA-seq of retina bipolar
cells, and macrophages are the major cell types found in PBMC neurons43. CRX is a cone-rod homeobox gene expressed in cone-
and the spleen. Additionally, several types of neuron cells and rod photoreceptor cells. Lastly, CDHR1 and RBP3 were identified in
microglia were identified in brain regions (Fig. 1f, Supplementary the photoreceptor cells (Fig. 2d, e).
Data 3). Moreover, several cell types commonly shared between In the pig kidney, collecting duct cells were characterized by
tissues were identified, such as ECs across nineteen different the specific expression of AQP3, GATA2, and AQP2. Distal
tissues (except PBMC); microglia mainly across all brain regions convoluted tubule cells were indicated by the specific expression
and retina; immune cells such as T cells, B cells, and NK cells of TMEM52B. Proximal tubule cells were demonstrated with the
across multiple different tissues (Fig. 1f). These results demon- specific expression of CUBN, LRP2, SLC13A3, and SLC34A1. In
strated that the main tissue-specific cell types were identified in addition, we identified podocytes (NPHS1+, NPHS2+, WT1+,
each tissue and some cell types were distributed across tissues. It and CLIC5+) and Loop of Henle cells (SLC12A1+). In addition
is consistent with the understanding that the cells derived from to these cells, we identified a few other cell types in the kidney.
the three germ layers are widely distributed within the human Briefly, fibroblasts had a high expression of CALD1, COL6A1, and
body34. To further facilitate the sharing and utility of the resource DCN1, and the T cells demonstrated enrichment of CD3E, CD3D,
generated by this study, we constructed a single-cell tran- CD2, and CD3G. The ECs showed an increased expression of
scriptome atlas database (https://dreamapp.biomed.au.dk/ PECAM1 and NRP1 (Fig. 2f, g, Supplementary Data 6). We next
pigatlas/). performed GO enrichment analysis of the selected cell types. Our
results show that pathways i.e., metal ion, sodium ion
transmembrane transport, cellular drug response were mainly
Validation of intra-tissue cell heterogeneity in the retina and enriched in collecting duct cells. The positive regulation of
kidney. To further explore and validate the intra-tissue cell het- sodium ion transmembrane transport, regulation of transmem-
erogeneity, we selected two tissues, retina and kidney, which have brane transport, and regulation of membrane potential were
been largely used as models for ophthalmology and kidney enriched in distal convoluted tubule cells. In addition, the cell-
diseases35–39. matrix adhesion, regulation of angiogenesis, ficolin-1-rich granule
The porcine retina contains several retina-specific cell types, such lumen, and myofibril were enriched in collecting duct cells, ECs,
as bipolar cells, cone photoreceptor cells, Müller glia, retinal and podocytes. Glomerular epithelial cell differentiation, renal
ganglion cells, and rod photoreceptor cells (Fig. 2a). The main cell filtration cell differentiation, nephron development, and glomer-
types shared in multiple tissues are microglia and T cells. Bipolar ulus development were specifically enriched in podocytes.
cells express high levels of TRPM1, PCP2, and GNG13. Cone Moreover, our single-cell analysis suggested that loop of Henle
photoreceptor cells exhibit a high level of ARR3. Microglia were cells plays an important role in potassium ion homeostasis,
indicated by the specific expression of C1QA, C1QB, CSF1R, and potassium ion import, chloride ion homeostasis, and metanephric
CD68. Müller glia demonstrated a specific expression of RLBP1 and nephron tubule development (Fig. 2h, Supplementary Data 6), in
CA2, and retinal ganglion cells were annotated with high expression line with its functions in maintaining iron and water
of NEFL, THY1, and NRN1. Rod photoreceptor cells had a high homeostasis44. To validate the canonical markers expressed in
expression of PDE6A, CNGA1, and SAG (Fig. 2b, Supplementary the kidney, we further analyzed them with immunohistochem-
Data 3, 4). Gene ontology (GO) analysis further showed that istry. Our results showed that NPHS2, a podocyte-specific
axonogenesis, axonal/axonal-dendritic transport, negative regulation marker, is specifically expressed in podocytes. SLC12A1 was
of neurogenesis, regulation of neuron projection development, and specifically expressed in the loop of Henle cells, while GATA2,
C
R
on
a b c
od
e
et
ph
ph
in
p.adjust Count
Retina
En
ot
ot
al
or
CSQA RA
or
ga
do
PE 3 2B
FL F 1
0.04 4
M
A
SA O 1
1
3
ec
1
ec
CDF1R
CNE6A
C1 A−D
CABP1
Bi
ARMK
RHGA
GNPM
CA G1
PCKC
CD3G
PRTO
VWCA
RL QB
PD N 1
CD3D
C D3 E
PR
TR P2
ng
th
C1 6 8
TH FL
NRY1
HLT1
R
CDG
PT 2
8
M
NE2
ep
po
ep
0.03
NE
el
lio
ul
ia
to
la
to
ic
12
le
T
n
0.02
ro
lc
rc
rc
rc
Bipolar cells
rg
ce
ce
gl
16
el
el
el
el
0.01
lia
lls
lls
ia
ls
ls
ls
ls
T cells axonogenesis
axonal transport
Cone photoreceptor cells axo−dendritic transport
anterograde axonal transport
photoreceptor cell cilium
9+0 non−motile cilium
Retinal ganglion cells photoreceptor outer segment
endothelium development
endothelial cell differentiation
endothelial cell development
Microglia vasculogenesis
endothelial cell migration
Muller glia blood vessel endothelial cell migration
macrophage differentiation
macrophage activation
microglial cell activation
glial cell activation
negative regulation of neurogenesis
Rod photoreceptor cells negative regulation of nervous system development
Endothelial cells negative regulation of cell development
gliogenesis
tSNE2
d e
RHO ARR3 CDHR1 RHO (Rod photoreceptor cells) ARR3 (Cone photoreceptor cells) GNG13 (Bipolar cells)
Cone&Rod photoreceptor cells
Rod photoreceptor cells
Cone photoreceptor cells
2.0 5
30
4
1.5
20
1.0
3 50 μm 50 μm 50 μm
2
10
0.5 1
0 0.0 0
f g h
D
is
ta
lc
U
Kidney
on
re
Pr
Lo
p.adjust Count
C
3 B 1 3 8 1
te
21 52 72 1
vo
1
1 A1 IN 3 2AS1 S2 5 N
ox
3A 2A 4A C
ol
M C2
op
3 2 2 BM 8M M N M1
r ic
P TA P L E N E E C G CA P R T2 LD L6 N IL IT C1 H H T1 IC B P2 C1 C2 C3 PR 3D 3E 3G 2 CAFD
lu
le
En
i
Q A Q VA LD 4
m
M N E R D S A O C M L L P P L U L L L T D D D D P
A G A P TM C TM TM E E P N K B C C D E S S N N W C C LR S S S P C C C C E W
of
te
ct
0.04
ep
Fibroblasts
al
do
in
d
8
H
ith
0.03
tu
g
tu
th
Fi
en
Po
T cells
el
bu
du
bu
12
el
br
le
d
0.02
ia
Endothelial cells
ia
ob
l
ct
oc
le
's
T
e
lc
16
lc
0.01
ce
la
ce
yt
ce
ce
ce
el
el
s
es
lls
lls
lls
lls
lls
ls
ts
ls
regulation of peptidase activity
response to metal ion
sodium ion homeostasis
Podocytes cellular response to drug
sodium ion transmembrane transport
positive regulation of sodium ion transmembrane transport
positive regulation of sodium ion transport
regulation of membrane potential
regulation of sodium ion transmembrane transporter activity
endothelial cell proliferation
Proximal tubule cells regulation of endothelial cell proliferation
cell−matrix adhesion
Collecting duct cells endothelium development
regulation of angiogenesis
ficolin−1−rich granule lumen
ficolin−1−rich granule
myofibril
contractile fiber
myofilament
potassium ion homeostasis
Ureteric epithelial cells potassium ion import across plasma membrane
potassium ion import
inorganic cation import across plasma membrane
cellular potassium ion homeostasis
chloride ion homeostasis
Loop of Henle's cells metanephric nephron tubule development
monovalent inorganic anion homeostasis
tSNE2
i j
NPHS2 SLC12A1 PAX2 NPHS2(Podocytes) SLC12A1(Loop of Henle’s cells) GATA2 (Collecting duct cells)
Podocytes
50 4
40
40
3
50 μm 50 μm 50 μm
30
30 2
20 20
10 10 1
0 0 0
10.0 2.0 50 μm 50 μm
7.5 1.5
5.0 1.0
2.5 0.5
0.0 0.0
AQP2, and PAX2 were specifically expressed in collecting duct which line the vascular systems and play important roles in e.g.
cells in the porcine kidney (Fig. 2i, j). regulating immune responses, regulation of blood fluidity, car-
diovascular homeostasis, maintenance of vascular functions, have
been extensively studied by scRNA-seq in human and
Endothelial cell heterogeneity. Single-cell transcriptome analysis mouse45–49, though little is known about the single-cell tran-
also enables us to probe rare cell types. Endothelial cells (ECs), scriptome and heterogeneity of ECs in pigs. We therefore focused
Fig. 2 Validation of cell heterogeneity in retina and kidney. a t-SNE visualization of major cell types in the retina. b Heatmap of marker gene expression in
the cell types captured in the retina. c GO term enrichment analysis on marker genes in the major cell types in the retina. The hypergeometric test was
used for GO term analysis, and p values were adjusted by Benjamini & Hochberg. d t-SNE visualization of the expression of RHO, ARR3, CDHR1, GNG13, CRX,
and RBP3 in the retina. e Representative IHC staining of RHO, ARR3, GNG13, CRX, CDHR1, and RBP3 in pig retina (n = 3). f t-SNE visualization of major cell
types in the kidney. g Heatmap of marker gene expression in major cell types in the kidney. h GO term enrichment analysis on marker genes of major cell
types in the kidney. The hypergeometric test was used for GO term analysis, and p values were adjusted by Benjamini & Hochberg. i t-SNE visualization of
NPHS2, SLC12A1, PAX2, GATA2, and AQP2 expression in the kidney. j Representative IHC staining of NPHS2, SLC12A1, GATA2, AQP2, and PAX2 in pig
kidney (n = 3). Source data are provided as a Source Data file.
on the ECs which were captured in 19 tissues in our datasets. To trajectory analysis which demonstrates that the EndMT cells
characterize the heterogeneity of ECs in pig organs, we extracted undergo a dynamic developmental transition from ECs to
ECs from all 19 tissues expressing the canonical EC marker mesenchymal cells (Fig. 3f). During the EndMT process,
PECAM1, and excluded epithelial cells (EPCAM+), immune cells expression of EndMT inducing genes (TGFB2 and SNAI1) were
(PTPRC+), fibroblasts (COL1A1+), and pericytes (PDGFRB+) high in the early EndMT process, followed by gradually
which frequently co-express PECAM1 (Fig. 3a, Supplementary decreasing expression of EC-specific genes (e.g., PECAM1,
Fig. S3a, Supplementary Data 7). In total, we obtained 9520 ECs VWF, ICAM1, and CDH5), and the increased expression of
from 19 tissues and 56% of all ECs were derived from the adipose mesenchymal cell-specific genes (e.g., ACTA2, TAGLN, CD44,
tissues (Supplementary Data 7). This is expected, as we applied a VIM, and CNN1) (Fig. 3g, h). Early findings suggested that TGF-
protocol optimized for enriching ECs from adipose tissues (see β2 is a key regulator of the EndMT process55,58–61. Our
methods). After batch correction, all ECs were separated into pseudotime analysis of EndMT cells also suggested that the
21 subtypes with significantly enriched genes (Fig. 3b, c, Sup- expression of the TGF-β signaling pathway (TGFB2) a key
plementary Fig. S3b, Supplementary Data 8). Although the inducing factor driving the EndMT process (Fig. 3h).
number of analyzed ECs from each tissue are largely different, the To verify some of the markers for EndMT cells at the protein
distribution of these 21 EC clusters was quite similar between the and histological level, we analyzed the co-expression of ECs
tissues (Fig. 3c). We performed gene ontology analysis of these 21 markers (PECAM1 and VWF) and mesenchymal cell markers
EC clusters based on enriched markers (Supplementary Data 8). (ACTA2 and TAGLN) by antibody-based immunofluorescence
Two EC clusters with distinct functions were highlighted here. staining. Our results showed that EndMT cells can be clearly
The EC cluster (c12) upregulates defense marker genes (i.e., identified in adipose tissues with co-expression of VWF and
complement: C1QA, C1QB, C1QC; cathepsins: CTSS, CTSD, TAGLN (Fig. 3i), as well as PECAM1 and ACTA2 (Fig. 3j), in a
CTSB, CTSZ; cystatins: CST3, CSTB (Supplementary Fig. S3c, small fraction of ECs. Corroborating the scRNA-seq results, only
Supplementary Data 8)), suggesting that it is an EC phenotype a very small fraction of the adipose tissue ECs is EndMT. We also
with immune active features. This observation is consistent the validated another ECs-specific marker FABP4 found by our
scavenging and immune-modulating EC phenotype reported single-cell analysis of pig adipose ECs. The fatty acid-binding
previously50–53. In adult tissues, most ECs are quiescent. We protein 4 (FABP4) is a lipid transport protein, which is expressed
identified a small EC cluster (c17) significantly up-regulating cell in adipocytes and capillary ECs62. To investigate if FABP4 is
proliferation genes (CENPF, CENPE, TOP2A, TPX2) (Supple- expressed in capillary ECs across pig tissues, we analyzed pig
mentary Data 8), which we define as proliferating ECs. The liver, heart, kidney, spleen, duodenum, jejunum, cerebellum, and
presence of proliferating ECs in the analyzed pig tissues agrees brain cortex by immunohistochemistry. The results show that
with the age of the pigs used for this study (Supplementary FABP4 is strongly expressed in capillary ECs of all included pig
Data 1). tissues (Supplementary Fig. S3e) and is highly expressed in the
As EC heterogeneity can be caused by tissue types and vascular subcutaneous adipocytes and the adipose ECs (Fig. 3k). Collec-
bed, we therefore focused on the ECs from adipose tissues and tively, our results support the previous model of EndMT (Fig. 3l),
further investigated the EC heterogeneity within tissues. We next which is involved in important mesenchyme-associated physio-
performed cell clustering based on the adipose tissue-derived ECs logical and pathological processes54.
and annotated each ECs subtype based on the expression of ECs
markers and pseudotime trajectory analysis. The ECs from
adipose tissues are mainly composed of blood ECs (arterial, Validation of EndMT in cultured ECs. Cultured ECs are com-
capillary, and vein) and lymphatic ECs (Fig. 3d, e), which also mon models for studying the EndMT process, of which differ-
seems to express higher level of mesenchymal genes (Supple- entiation of ECs into mesenchymal cells can occur spontaneously
mentary Fig. S3a). We also identified three functionally distinct or through TGF-β induction63. Previously, we also identified the
phenotypes: proliferating EC, immune active EC and an EC EndMT ECs phenotype in cultured human lung tumor ECs46. To
phenotype with co-expression of the mesenchymal cell markers further validate this EndMT ECs phenotype, we isolated and
ACTA2 and TAGLN (Fig. 3d, Supplementary Data 9). Both cultured primary ECs from pig lung and aorta and analyzed them
immune active and proliferating ECs have been described above by scRNA-seq (Fig. 4a). In total, after filtering low-quality cells
and in previous scRNA-seq studies of mouse and human and doublets, 5698 and 882 cells were obtained from the cultured
ECs45–47. Here we focus on the small fraction of mesenchymal- ECs of pig lung and aorta respectively, which were further clus-
like endothelial phenotype. Previous studies have suggested the tered into 5 clusters (Fig. 4b–f) based on the expression of cell-
existence of such a specific endothelial phenotype, which is type-specific markers (Fig. 4c–e) and cell cycle analysis (Fig. 4f).
undergoing the processing of endothelial-to-mesenchymal transi- In cultured ECs from both pig lung and aorta, we identified three
tion (EndMT)54–57. The number of expressed genes per cell is clusters of proliferating ECs (G1, S and G2M phase) and an
similar between EndMT cells and other EC phenotypes intermediate EC phenotype. Because the ECs isolation protocol is
(Supplementary Fig. S3d), confirming that it is a cluster of based on enzymatic perfusion of the blood vessels (see methods),
transcriptomically distinct EC subtypes which is not caused by we also identified a cluster of fibroblasts expressing high levels of
doublets of an EC and a mesenchymal cell. To further validate the COL1A1 and COL1A2 (Fig. 4b–e). In the cultured aorta ECs, we
mesenchymal-like EC phenotype, we performed pseudotime identified a cluster of cells that highly express mesenchymal cell
a b c
Immune active EC
Tissue
Adipose-S 4000 3000 2000 1000 0 0 25 50 75 100
Adipose-V 973 Adipose−S
Brain
Proliferating EC
4365 Adipose−V
Heart
1060 Brain
Intestine
Kidney 172 Heart
tSNE2
149 Retina
4 133 Spleen
3
2
tSNE1
1
0
d Lymphatic ECs e f
EndMT cells 5 10
1
Large vein ECs
2
0 5
Large artery ECs
Capillary-venous ECs
Artery ECs Vein ECs
−5 0 2
3
Proliferating ECs
Immune active ECs
−10 −5 0 5 −10 0 10
tSNE1
Large artery ECs Capillary ECs (1) Capillary-venous ECs Large vein ECs Immune active ECs
Artery ECs Capillary ECs (2) Vein ECs Lymphatic ECs EndMT cells Proliferating cells
g h
TGFB2 PECAM1 ACTA2 NOTCH3
5.0 100
5.0 5.0
3.0
3.0 3.0
10
1.0 1.0 1.0
1
0.5 0.5 0.5
0 10 20 30 40 0 10 20 30 40 0 10 20 30 40 0 10 20 30 40
SNAI1 VWF TAGLN FABP4
30
1.000 5.0 30
3.0 10
0.100 10
3 3
0.010 1.0
1 1
0.001 0.5
0 10 20 30 40 0 10 20 30 40 0 10 20 30 40 0 10 20 30 40
ZEB2
5.0
3.0
l
1.0
0.5
0 10 20 30 40
i j k
Adipose tissue Adipose tissue Adipose tissue
VWF TAGLN PECAM1 ACTA2 Control
markers (TAGLN and ACTA2) but not fibroblast markers To further validate that EndMT can be induced by TGF-β2, we
(COL1A1 and COL1A2), which we defined as EndMT-like cells treated cultured pig aorta endothelial cells (PAECs) with TGF-β2
(Fig. 4d, e). The cultured lung and aorta ECs used for scRNA-seq (2 ng/mL) and measured the expression mesenchymal cell marker
were maintained in normal ECs growth medium. This EndMT- ACTA2 and ECs marker CD31 by fluorescence-conjugated
like cells might be induced by the presence of TGF-β2 in serum or antibody staining and flow cytometry analysis (Supplementary
contaminating mesenchymal cells from the isolation procedure. Fig. S4a). Increased ACTA2 expression was already detected in
Fig. 3 Single-cell transcriptome analysis of porcine ECs in pig tissues. a Heatmap of gene expression in the selected PECAM+/EPCAM−/PTPRC−/
COL1A1−/PDGFRB− ECs from 19 tissues. b t-SNE visualization of 21 EC clusters, which are colored according to EC cluster numbers. c Bar graph shows
numbers and percentage of ECs in tissues. ECs were colored according to clusters. d EC subtypes in adipose tissues. Cells are color-coded according
to the subtypes of ECs. e EC subtypes trajectory analysis in adipose tissues (adipose-S and adipose-V) using monocle 2 and cells on the tree are
colored by EC subtypes. f EndMT EC trajectory analysis using monocle 2 and cells on the tree are colored by EndMT states. g Pseudotime trajectory
analysis of marker genes in different states of EndMT cells. h Representative pseudotime trajectory of marker genes (TGFB2, SNAI1, PECAM1,
VWF, ZEB2, ACTA2, TAGLN, NOTCH3, FABP4) in different states of EndMT ECs. i Representative Immunofluorescence staining images of VWF and
TAGLN in adipose tissues. Arrow (red) indicates ECs expressing both VWF and TAGLN. Arrow (green) indicates ECs only expressing VWF (n = 3).
j Representative Immunofluorescence staining images of PECAM1 and ACTA2 in adipose tissue. Arrow (red) indicates ECs expressing both PECAM1
and ACTA2. Arrow (green) indicates ECs only expressing PECAM1 (n = 3). k IHC of FABP4 in ECs and adipose tissues. Arrows (green) indicate ECs.
Arrows (red) indicate adipocytes. Control was stained with an antibody against a gene not expressed in adipose tissues (n = 3). l An integrated model
of EndMT process. Source data are provided as a Source Data file. Schematic diagrams in I were created with BioRender.com.
cultured PAECs 2 days after TGF-β2 treatment (Supplementary podocytes, and proximal tubule cells through the VEGF pathway
Fig. S4). After culturing PAECs in medium with TGF-β2 for five (Fig. 5a, Supplementary Fig. S5a). The analysis suggests that the
days, the expression of ACTA2 was significantly increased in VEGF signaling pathway can be used by most renal cells for
PAECs by 6 folds compared to untreated controls (Fig. 4g, intercellular communication. Unlike the VEGF signaling path-
P = 1.44E-7, t-test). The CD31 expression was slightly decreased way, ECs only communicate with a few other cell types through
in the TGF-β2 treated PAECs, but not significant (P = 0.300, t- the PDGF pathway. For example, pig ECs communicate with
test). Our results also showed that TGF-β2 significantly inhibits immune cells (Kupffer-, T-, NK-, and B cells) in the liver, with all
the proliferation of ECs (Supplementary Fig. S4b, P = 1.33E−4, t- cell types in the heart, and with only podocytes in the kidney.
test). Consistent with the cultured PAECs, similar TGF-β2 Human ECs communicate with hepatocytes and Kupffer cells in
induced EndMT results were also obtained in cultured human the liver, with fibroblasts in the heart, and with podocytes and
umbilical vein endothelial cells (HUVECs) (Fig. 4h). Together, proximal tubule cells in the kidney through the PDGF pathway
our results demonstrated that EndMT is driven by the TGF-β2 (Fig. 5b, Supplementary Fig. S5a). Similarly, analysis of TGF-β
signaling pathway and is conversed between pig and human. and BMP pathways-mediated cell–cell communications show
similar cell-type-specific and tissue-specific preferences, as well as
Communications between ECs and other cell types. Intercellular some divergence between pig and human (Supplementary
communication between ECs and other parenchymal cells in the Fig. S5b). The presence of TGF-β-based cell–cell communication
tissue plays a vital role in the structure and function of maintaining between ECs and other tissue types further suggests a potential
normal tissue growth, development, and homeostasis64. We focused cellular mechanism in inducing EndMT in tissues. However, it
on the ECs and analyzed the communications between ECs and should be noted that the inferred cell–cell communications based
other cell types using CellChat, a tool for quantitatively inferring and on single-cell RNA sequencing could be affected by e.g.,
analyzing intercellular communications based on the differential posttranscription or posttranslational protein modifications,
expression of ligand and receptor gene pairs65. For this analysis, we completeness and biases in the genome annotation between pig
used single-cell transcriptome data from the liver, kidney, and heart, and human. We demonstrated that the pig single-cell RNA atlas
which are the three most studied pig organs in biomedical research66. provides a valuable resource for inferring cell–cell communica-
We performed a comparison with corresponding scRNA-seq data- tions within pig tissues or between pig and human tissue.
sets of the human liver, kidney, and heart (see methods) to gain a
better understanding of the cell–cell interactions in these three MEF2C is a conserved regulon in microglia evolution. Single-
organs between pig and human. We identified six cell types in the cell transcriptomic analysis not only provides good insights into
liver (hepatocytes, ECs, Kupffer cells, B cells, T/NK cells, and ery- the cellular heterogeneity and functional diversity of structural
throid cells), six cell types in the kidney (epithelial cells, podocytes, cell types across tissues, but it is a good way to uncover the
proximal tubule cells, collecting duct cells, ECs, and distal convoluted similarities and divergences of cell types across species. In this
tubule cells), and five cell types in the heart (ECs, fibroblasts, car- study, we utilized our single-cell pig atlas to analyze the cross-
diomyocytes, lymphoid cells, and myeloid cells) isolated from both tissue ECs heterogeneity and ECs conversion in adipose tissues.
pig and human tissues (Supplementary Data 10, 11). We were also interested in the cross-species cell types, which are
We next investigated signaling interactions based on ligand- mainly focused on microglia in the brain across 13 species. It
receptor pairs and performed a global communication pattern enabled us to better understand cell-type evolutions. The pig
recognition analysis to identify the key signal communications single-cell transcriptome atlas includes nine brain regions
involved in ECs. Our results suggested that ECs use the vascular (Fig. 6a). Most cell types in the brain were clustered based on the
endothelial growth factor (VEGF)67, PDGF68, TGF-β69, and regions (Fig. 6b, Supplementary Fig. S6a). We also validated a few
BMP70 signaling pathways as major communicating pathways cell types by immunohistochemistry, confirming that scRNA-seq
with other cell types in the liver, kidney, and heart in both human is a robust method for the classification of different cell types
and pig. However, single-cell-based cell communication analysis (Fig. 6c, Supplementary Fig. S6b), including SLC1A6 which is a
provided several unique insights regarding tissue-specific and marker for Purkinje cells71 expressed in pig cerebellum; CALB2
cell-type-specific communications between ECs and other cell which is an inhibitory interneuron marker; PAX6 which is a
types. The VEGF signaling is the most studied pathway in ECs highly specific marker for granule cells in cerebellum. By protein
activation, proliferation, and angiogenesis. In the liver and heart, staining, we also detected SLC17A7 in the excitatory synapse in
our results showed that there are similar sender-receiver the pig hippocampus. Lastly, protein staining of SST is consistent
communications between ECs and other cell types. In the kidney, with scRNA-seq suggesting that SST is a specific marker for
the VEGF signaling pathway is used by most renal cells for inhibitory interneurons in the pig cortex (Fig. 6c).
cell–cell communications. Particularly, similar between pigs and We focused on microglia which is important for brain
humans, there is a higher level of communication between ECs, homeostasis and involved in several brain disorders i.e.
Average expression
Lung (perfused) Aorta (near to lung) ICAM2
RPL3 1
0
RPL7A
ECs (G1) −1
Features
TOP2A
Dissociation PCNA
Proliferating ECs (S)
MCM6
Percent expressed
Digestion ACTA2
Intermediate ECs TAGLN
25
50
tSNE2
75
S100A4 100
Culture
Proliferating ECs (G2M) COL1A1
COL1A2
scRNA-seq
1)
S)
s
tSNE1
2M
EC
st
G
s(
la
s(
G
EC
ob
te
EC
s(
ECs (G1) Proliferating ECs (S) Fibroblasts
ia
br
EC
g
id
tin
Fi
rm
g
ra
tin
te
ife
Intermediate ECs Proliferating ECs (G2M)
In
ra
ol
ife
Pr
ol
Pr
d Cultured aorta ECs e f Percentage (%)
PECAM1 0 25 50 75 100
Average expression
VCAM1 ECs(G1)
Proliferating ECs (S/G2M) ICAM2 1
Intermidiate ECs
Lung
RPL7A 0
Percent expressed
MCM6 25
Stage
Intermediate ECs ACTA2 50 G1
Percentage (%)
TAGLN 0 25 50 75 100
S
75
S100A4 G2M
tSNE2
ECs(G1)
COL1A1 100
Intermediate ECs
Aorta
COL1A2
EndMT-like cells
Proliferating ECs(S/G2M)
s
lls
s
1)
2M
EC
st
ce
G
la
G
s(
EndMT-like cells
ob
te
tSNE1
e
S/
EC
lik
ia
br
s(
ed
T-
Fi
EC
dM
Fibroblasts
rm
g
te
En
ra
ife
g h
Induced EndMT in cultured PAEC by TGF-β2 at day 5 Induced EndMT in cultured HUVEC by TGF-β2 at day 5
unstain 4 103
0 20 40 60 80 100
control P = 1.44E-7
control P = 3.91E-6
TGF-β2
MFI (-BG)
MFI (-BG)
Count (%)
TGF-β2
Count (%)
3 103
4 104
2 103
2 104
1 103
2 3 4 5 6 7 0 2 3 4 5 6 7 0
10 10 10 10 10 10 10 10 10 10 10 10
Fluorescent intensity Control TGF-β2 Fluorescent intensity Control TGF-β2
CD31 CD31 CD31 CD31
unstain unstain
0 20 40 60 80 100
0 20 40 60 80 100
MFI (-BG)
MFI (-BG)
4 104 4 105
2 104 2 105
0 2 4 6 2 4 6
10 10 10 10 0 10 10 10 0
Fluorescent intensity Control TGF-β2 Fluorescent intensity Control TGF-β2
Fig. 4 Validation of EndMT ECs in cultured ECs from pig lung and aorta by scRNA-seq. a The flowchart of ECs isolation and culture from pig lung and
aorta. b tSNE visualization of cell-type annotations from cultured lung ECs. c Dot-plot of selected canonical marker genes for cell-type annotations from
cultured lung ECs. d tSNE visualization of cell-type annotations from cultured aorta ECs e Dot-plot of selected canonical marker genes for cell-type
annotations from cultured aorta ECs. f Cell cycle analysis of annotated cell types in cultured lung and aorta ECs. g Expression of ACTA2 and CD31 in
cultured PAEC treated with TGFb2 for five days (n = 3, two-sided t-test). Values are presented as mean±SD. h Expression of ACTA2 and CD31 in cultured
HUVEC treated with TGFb2 for 5 days (n = 3, two-sided t test). Values are presented as mean ± SD. Source data are provided as a Source Data file.
VEGF Kupffer cells Endothelial cells Endothelial cells Proximal tubule cells Podocytes
Fibroblasts
Myeloid cells
Distal convoluted tubule cells
Erythroid cells Lymphoid cells Endothelial cells
T/NK cells
1
Importance
1 Sender
Importance
1
Importance
Sender Sender
Receiver
Receiver Receiver
Mediator
Mediator Mediator
Influencer
Influencer Influencer
0 0 0
lls
ls
lls
lls
ls
le d
te
el
lls
bu te
lls
ce
el
s
lls
ce
ce
lls
s
ls
ls
ls
lls
es
te
lc
cy
st
lc
tu lu
ce
ce
el
el
el
ce
ce
ce
yt
cy
ct
le
la
vo
ia
do
lc
lc
rc
ia
oc
du
d
bu
ob
K
el
yo
B
el
id
on
oi
ia
ia
Po
oi
fe
th
at
ith
lo
tu
br
om
ph
el
el
hr
lc
pf
T/
do
ep
ye
in
Fi
th
th
Ep
al
yt
Ku
ta
di
ct
En
H
do
do
im
Er
Ly
is
ar
le
D
En
En
ox
ol
C
C
Pr
Endothelial cells Endothelial cells Proximal tubule cells Podocytes
Kupffer cells
Fibroblasts
Importance
Sender 1
Importance
Sender Sender
Receiver Receiver Receiver
Mediator Mediator Mediator
Influencer Influencer Influencer
0 0 0
lls
ls
lls
lls
ls
ls
s
s
lls
lls
le d
s
lls
ls
lls
es
ls
lls
te
te
el
bu te
st
ce
el
el
ce
ce
el
el
ce
ce
ce
ce
yt
lc
cy
ce
cy
la
lc
lc
tu olu
lc
rc
oc
ct
le
ob
ia
yo
d
K
do
d
ia
ia
id
ia
B
v
fe
oi
du
bu
oi
N
el
at
el
el
lo
on
br
om
el
Po
ph
pf
T/
hr
th
ep
ith
th
ye
tu
Fi
th
lc
Ku
yt
do
di
in
do
Ep
do
H
al
ta
Er
ar
Ly
ct
En
im
En
En
is
le
C
D
ox
ol
C
Pr
b
PDGF Kupffer cells Endothelial cells
Endothelial cells
Proximal tubule cells Podocytes
Fibroblasts
1
Importance
Importance
Sender Sender Sender
Receiver Receiver Receiver
Mediator Mediator Mediator
Influencer Influencer Influencer
0 0
0
lls
ls
s
ls
lls
lls
lls
ls
lls
s
ls
s
lls
lls
le d
es
ls
lls
te
st
te
el
el
bu ute
ce
el
el
ce
ce
el
ce
ce
ce
ce
yt
cy
ce
la
lc
lc
cy
lc
lc
rc
oc
l
ob
ct
le
yo
vo
K
ia
id
d
do
ia
ia
ia
B
fe
oi
oi
du
N
bu
at
el
lo
el
br
el
om
on
el
Po
ph
pf
T/
hr
th
tu
ep
ye
th
ith
tu
Fi
th
lc
Ku
yt
di
do
do
in
do
M
H
Ep
al
Er
ta
ar
Ly
ct
En
En
im
En
is
C
le
D
ox
ol
C
Pr
Kupffer cells Endothelial cells Endothelial cells Proximal tubule cells Podocytes
Fibroblasts
1
Importance
1
Importance
1
Importance
s
ls
le d
s
lls
s
lls
ls
lls
ls
lls
te
el
es
te
bu te
ce
el
st
el
lls
ce
ce
el
el
ce
ce
ce
lc
cy
ce
cy
lc
tu lu
la
lc
yt
ce
lc
rc
ct
le
vo
ia
ob
do
oc
ia
yo
d
d
ia
id
K
ia
du
bu
fe
oi
el
B
oi
el
on
el
lo
N
br
Po
at
om
el
ph
pf
th
hr
ith
tu
T/
th
ye
g
ep
Fi
th
lc
Ku
do
yt
in
di
m
do
Ep
al
do
ta
H
Er
ct
ar
Ly
En
im
En
is
En
le
C
D
ox
ol
C
Pr
Fig. 5 Comparison of cell communication and signaling pathway between pig and human. a Comparison of cell–cell communication of VEGF signaling
pathway in liver, heart, and kidney. b Comparison of cell–cell communication of PDGF signaling pathway in liver, heart, and kidney. Source data are
provided as a Source Data file.
amyotrophic lateral sclerosis72. Microglia are the primary resident evolutionary module in human brain diseases73. To characterize
immune cells of the brain, which play a critical role in many the conservation and divergence of microglia in pigs, we analyzed
physiological and pathological brain processes. The conserved the single-cell microglia datasets from 13 different species
and divergent microglia gene program in pan-species provides spanning more than 300 million years of evolution. The species
important implications for investigating the microglia contain bearded dragon lizard, turtle, hedgehog, mink, alpaca,
Fig. 6 Validation of brain cell types and comparison of microglia regulome across species. a Schematic diagram visualizing the 9 different brain regions
analyzed by this study. b t-SNE visualization of cells in the 9 brain regions. Cells were color-coded according to brain regions. Cell types were shown in
extended Fig. S6. The hypergeometric test was used for GO term analysis, and p values were adjusted by Benjamini & Hochberg. c Representative IHC of
SLC1A6, CALB2, PAX6, SLC17A7, and SST by antibody staining in the cortex or cerebellum. Arrows indicate corresponding cell types (n = 3).
d Phylogenetic tree based on the NCBI taxonomy of animals used in this study (generated via http://www.timetree.org//). e Violin plots visualizing the
expression of canonical microglia markers in the 13 species. f Conserved genetic regulatory networks in microglia within each indicated species. Light blue
nodes represent regulators and green nodes represent corresponding target genes. g Feature plots visualizing the clustering of microglia single-cell
transcriptome and expression of four TFs among the 13 species. Source data are provided as a Source Data file. Species icons in d, f were created with
BioRender.com.
sheep, pig, chinchilla, mouse, hamster, blind mole rat, macaque, of pathways such as adaptive immune response, regulation, or
and human (Fig. 6d, Supplementary Data 12). The microglia- humoral immune response in the immune system, and the central
specific marker genes C1QB and C1QC were expressed in most nervous system neuron differentiation and development for
species except for Chinchilla, Hamster, and Hedgehog, while the neuron-related pathways (Supplementary Fig. S6c). Collectively,
markers of CSF1R, APOE, and P2RY12 are covertly expressed in these results indicate that the combinations of multiple TFs
all 13 species (Fig. 6e), suggesting conservation and divergence in regulate microglia development and maintain the functional
the expression of microglia markers across species. states of microglia. Furthermore, we investigated the expression
Transcription factors (TFs) have been demonstrated as the level of conserved TFs of MEF2C, SPl1, IRF8, and ZFP36L1 across
important regulators of gene expression and with the ability to species. The results showed these conserved TFs are highly
shape different phenotypes of microglia74. We, therefore, expressed in microglia in these 13 species (Fig. 6g). Collectively,
performed single-cell genetic regulatory network (GRNs) infer- we provide a valuable resource of the conserved and divergent
ring and clustering analysis to assess TFs underlying differential GRNs program for microglia evolution across species, with
gene expression in microglia across species. In total, 1590 important implications for future development of the microglia
conserved TF-target pairs were identified, which were observed functional studies in the brain.
in at least five of thirteen species. Interestingly, two pairs of TF-
target were conserved in ten different species, which are
MEF2C_P2RY12 and MEF2C_ZFP36L1. The target P2RY12 is a Discussion
microglia-specific gene, while the target ZFP36L1 is reported as a Pigs are important large animal models for studying complex
pivotal regulator for microglial fate specification75,76. For human diseases, as well as promising alternative organ donors for
TF_target pairs covering nine species, four pairs of conserved humans due to their high similarities with humans: physiology,
TF_target were identified including MEF2C_C1QB, MEF2C_- anatomy, genetics, metabolism, and organ size4. Single-cell
FAU, MEF2C_RGS10, and MEF2C_SERBP1. The target C1QB is transcriptomic profiles of multiple organs in the mouse and
a microglia-specific marker. The target RGS10 is reported as a key human body reveals the cellular compositions and heterogeneity
regulator of proinflammatory cytokine produced in microglia for of inter-and intra- organs and offers the opportunity to overall
neuroprotective factors77. SERBP1 plays an important role in organ development, physiology, and plasticity22,34,80,81. This
proinflammatory TLR4 signaling78. However, the FAU gene study generated a multiple-organ single-cell transcriptomic atlas
encodes a ubiquitin-like protein fused to the ribosomal protein of pigs covering twenty tissues. The landscape profiles of pig
S30. These results demonstrated that MEF2C is a core TF organs depicted the transcriptomic cellular heterogeneity in each
regulating multiple key genes related to microglia functions. tissue and expand the functions of cross-tissue cell types and rare
Furthermore, the TF-target pairs in eight species share the unique cell-type identification. With a focus on endothelial cells, we
TF of MEF2C, and most of the conserved TFs in at least five identified the subpopulations of cross-tissue ECs, such as blood
species also mainly contain MEF2C. These targets of MEF2C ECs, lymphatic ECs, and several subtypes of functionally distinct
include most microglia and immune-related genes such as ECs. These cell types were also reported in the human single-cell
PTPRC, CSFER1, APOE, AIF1, CD14, and CTSS et. These results atlas. We are also able to identify functionally specific ECs, such
demonstrated that MEF2C is a conserved TF in microglia as the immune active ECs that express typical ECs markers
evolution. It is consistent with the functional TF of MEF2C (PECAM1+) and immunomodulating genes (CD68+, C1QA+,
reported regulating the fundamental functions of microglia79. The C1QB+, and AIF1+). These findings of immune active ECs are in
other importantly conserved TFs were identified for microglia good agreement with our previous scRNA-seq based ECs tax-
functions such as SPI1 and IRF8 (Supplementary Data 13). onomy in mouse and human tissues45–47 and the well-study
Subsequently, we separately analyzed the conserved TF-target immune-modulating functions of ECs82–84. In this study, the ECs
pairs shared by five of thirteen different species. For example, in were only enriched from adipose tissues through a modified
humans, the top ten of these conserved TFs were shared with four isolation protocol. The degree of ECs heterogeneity was not
other random species are ARID1, ATF4, BCLAF1, CCDC88A, revealed as detailed as the mouse ECs atlas45. To systematically
CNBP, CSDE1, EGR1, ELF1, FLI1, and FOS. In the macaque, the characterize and compare the ECs heterogeneity across different
top ten conserved TFs are ELF1, FLI1, IKZF1, IRF8, KLF6, LITAF, tissues and organs by scRNA-seq, future studies are needed to be
MEF2A, MEF2C, NFE2L2, and SFPQ. In pigs, the top ten performed using enriched the ECs by MACS and/or FACS.
conserved TFs are MEF2C, IRF8, MEF2A, SON, SPl1, ATF4, Angiogenesis is matured through ECs and non-ECs to form the
ELF1, FLI1, FOS, and HMGB1. However, for reptiles such as vascular channels. The non-endothelial microcirculation is called
lizards, the top ten conserved TFs are MEF2C, ATF4, IRF8, vascular mimicry85. In physiologic and pathological angiogenesis,
MEF2A, MEF2C, ZEB2, SPl1, ATF4, FOS, and HMGB1. These macrophages are thought to play a supportive role to promote
conserved TFs in each species contain the shared conserved TF- vascular mimicry outgrowth through cytokine secretion and
target pairs and species-specific TF-target pairs (Fig. 6f, Supple- remodeling of the extracellular matrix (ECM)86,87. In tumors,
mentary Data 13). Additionally, the rank top ten conserved TFs macrophages are believed to reprogram the pathological angio-
(MEF2C, ZFP36L1, FOS, MEF2A, IRF8, SON, HMGB1, ZEB2, genesis to serve as the major source of angiogenic factors88.
SPI1, and ATF4) in humans were indicated with the enrichment PECAM1 signaling participates in the regulation of leukocyte
detachment, T cell activation, and angiogenesis89. Therefore, the We also highlight a few limitations of the current study. First,
high expression of macrophage marker genes and PECAM1+ cell tissues used by this study are taken from animals of adult ages.
population may play an important role in vascular mimicry Thus, the development and age effects on gene expression and
during angiogenesis. EndMT ECs were identified in the pig adi- cell-type compositions in tissues cannot be addressed by the data
pose tissues and cultured ECs with co-expression of ECs markers generated by this study. When conducting the comparison of cell-
of PECAM1 and VWF, and the mesenchymal cell markers of type-specific gene expression between pig and human, effects of
ACTA2 and TAGLN. In particular, the EndMT cells underwent ages, gender, and physiological conditions cannot be addressed by
the dynamic transition stages with the decreased expression of the this study. Second, the scRNA-seq and snRNA-seq experiments
ECs markers PECAM1 and VWF while the expression of the were conducted by two laboratories, using two different sets of
mesenchymal cell markers ACTA2 and TAGLN increased. This tissues, and with technologies provided by two companies. Thus,
phenotype is consistent with the stages of EndMT: (1). EndMT is the difference in the number of cell types and fraction of each cell
initiated by autocrine and/or paracrine inflammatory signals such type captured by the two methods could be caused by the many
as TGF-β2, which was also validated by our TGF-β2 induction steps during sample processing. Since the aim of this study was
experiment of cultured PAECs and HUVECs, or response to not focused on comparing the scRNA-seq and snRNA-seq
vascular injury; (2). Transitioning ECs acquire a migratory phe- technologies, we carried out comprehensive batch correction
notype, invade under the vascular basement membrane, and and normalization to use both scRNA-seq and snRNA-seq
begin to express mesenchymal markers, such as ACTA2; (3). Cells datasets to construct the first pig single-cell transcriptome atlas.
that have undergone EndMT have lost their endothelial pheno- Most importantly, the transcriptome of sample cell types cap-
type. These EndMT-derived cells contribute to the local tured by scRNA-seq and snRNA-seq is highly correlated, sug-
mesenchymal lineage population and are likely to produce var- gesting that both methods can capture the transcriptome of cells
ious growth factors, such as TGF-β290. EndMT is an important with high fidelity. Since the number and composition of cell types
developmental process, participating in tumor formation, inva- could be affected by the single-cell transcriptomics analysis
sion, and metastasis, and has been extensively participating in methods (scRNA-seq and snRNA-seq), comparison of cell-type
several diseases by causing morphology changes and pathological abundance and composition between tissues or conditions should
processes. The mesenchymal cells derived from EndMT might be be carried out using datasets generated with the same method to
differentiated into various mesenchymal cell types for tissue avoid method-induced biases. Lastly, we provide the first insight
engineering and subsequent transplantation into the patient91. into EC heterogeneity in pig tissues. To fully characterize the EC
The most potent pathway of VEGF played the core role in the heterogeneity and identify functionally distinct EC phenotypes in
intercellular interactions between ECs and tissue-specific cell pig organs, an EC-focused single-cell transcriptome analysis
types, suggesting active interactions and communications should be carried in future study.
between ECs and other cell types in organs92. Angiogenesis In summary, we constructed the first single-cell transcriptomic
occurs in both physiologic and pathological conditions and atlas of pig organs (https://dreamapp.biomed.au.dk/pigatlas/). We
interplays with other cell types93,94. In adult tissues, most ECs are identified the tissue-specific cell types and cross-tissue cell types
quiescent, but these ECs are metabolically active and actively such as ECs, immune cells, and microglia. Moreover, some rare
involved in the regulation of several important cellular processes cell types were identified in our data, such as EndMT ECs, sug-
such as immune modulations95. In addition, the growth of gesting the rare cell types could be identified by single-cell tech-
pathological angiogenesis in human diseases such as cancers niques. The regulatory mechanism analysis of microglia across
highlights that the targeting this process should help to reduce species provides insights into the conserved TFs regulatory
both morbidity and mortality from carcinomas96. Hence, anti- module for microglia evolution. Together, our study offers an
angiogenic therapy is a novel approach for the treatment of important resource for a better understanding of pig biology,
cancers, diabetic retinopathy, and other angiogenesis-dependent xenotransplantation, evolution, development, and regenerative
diseases97. medicine research.
Pigs are an excellent model for studying genetic and somatic
evolution23. In this study, we demonstrated that the evolution of Methods
microglia, a cell type that exhibits increasing interest and is Ethical statement. The study was approved by the Institutional Review Board on
important in neurological disorders98. Approximately 1500 reg- the Ethics Committee of BGI (Approval letter reference number BGI-NO.BGI-
ulatory sequence-specific DNA-binding factors (transcription IRB18135-T1). All experimental procedures were conducted following the national
and institutional guidelines of using the experimental animals for research. All the
factors, TFs) are encoded in the human genome, which is upre- applicable institutional and national guidelines for the care and welfare of animals
gulated in a tissue-specific and cell-type-specific manner. Changes have been strictly followed for the tissue sampling procedures.
in gene expression between species could be due to changes in the
TFs and/or changes in the instructions within the regulatory Tissue dissociation and sample preparation for scRNA-seq library generation.
regions of specific genes. TF expression patterns and binding Porcine tissues for scRNA-seq library generation were collected from a local
activities could advance the understanding of how tissue specifi- slaughterhouse (Sus scrofa domesticus, three-way hybrid of Landrace, Large White
city and conserved regulatory functions across species99–101. and Duroc, age 6 months, Hårby Slagteren IvS). The fresh tissues were collected,
immediately placed on ice, and processed within 30 min. Each tissue was dis-
Here, we investigated the genetic regulatory networks of microglia sociated and digested independently. To ensure efficient digestion and cell viability,
across thirteen different species, spanning more than 300 million after tissue dissociation (as described for each organ below), cell suspensions were
years of evolution. The conserved TFs across different species filtered by cell strainers after debris removal and red blood cell lysis. The cell
demonstrated the regulatory mechanism of target genes for the viability (Hoechst 33342 (Invitrogen, Cat#H3570) and Propidium Iodide (Thermo
Fisher Scientific, Cat#P3566)) from each tissue were evaluated by flow cytometry
fundamental functions of microglia. In particular, the conserved (FC) analysis on a NovoCyte Quanteon analyzer (Acea bioscience, Inc, US) pro-
TFs of MEF2C were detected in ten species, and its regulatory vided by the FACS Core Facility, AU.
genes are important for microglia functions in the brain such as
P2RY12, ZFP36L1, RGS10, and SERBP1. These transcriptional Liver. The fresh liver was collected and sampled with five different anatomical
regulatory genes play an important role in microglia biology, such regions: left lateral lobe (LLL), left medial lobe (LML), right medial lobe (RML),
right lateral lobe (RLL), and quadrate lobe (QL). For each region, 1 g tubes were
as microglia fate specification, a key regulator of proinflammatory punched and washed twice with cold PBS. The samples were mixed and carefully
cytokine produced in microglia, suggesting that microglia per- dissected into small pieces, and then transferred into a 50 mL tube with 20 mL
form overall similar functions during species evolution73. digestion medium containing 0.5 mg/mL collagenase type II (Gibco,
Cat#17101015), 1.25 mg/mL protease (Sigma, Cat#P5147-100MG), 7.5 μg/mL Cultured ECs isolation from pig lung. One leaf of pig lung was used for perfusion
DNase I (Sigma-Aldrich, Cat#D4527-10KU) in cold HBSS. The tissue digestion with PBS from the big vessel on top. After perfusion, the lung was squeezed and
was performed at 37 °C for 15 min with gently shaking once every 5 min. The extracted the PBS. Fifteem milliliters of the digestion buffer containing 0.1% col-
reaction was stopped in 20 mL cold MACS buffer containing 0.25% BSA (Sigma- lagenase II, 0.25% collagenase IV, 75 μL DNAse I in KnockOutTM DMEM medium
Aldrich, Cat#10735096001), 2 mM EDTA, in PBS, and then filtered with 100 μm was injected into the sealed pig lung. The opening site was closed by sealing clip
cell strainer (Sigma-Aldrich, Cat#CLS431752-50EA), and the cells were stained for and left in the bag into a water bath at 37 °C for 30 min. Next, the lung was taken
the FC analysis. out and the digestion buffer has flowed away. Twenty milliliters of PBS was used to
wash the lung thrice, and the medium was collected for filtering with 100 μm cell
Spleen. The fresh spleen was collected and sampled with 1 g pieces for both ana- strainer afterward. The medium was transferred into 10 mL Falcon tube and
tomical sides including the parietal and visceral side and washed twice in cold PBS. centrifuged at 300 × g for 5 min. The cells were washed with 10 mL PBS and
The samples were mixed and carefully dissected into small pieces and transferred centrifuged for culturing into 5 mL PAEC medium (DMEM with 10% fetal calf
into a 50 mL tube with 10 mL digestion buffer containing 20 mg/mL collagenase IV serum, 1 × MEM-NEAA, 1 mM sodium pyruvate, 1 × glutamax, 1 × penicillin/
(Gibco, Cat#17104019), 1 U/mL Dispase II (Gibco, Cat#17105041), 7.5 μg/mL streptomycin) in CO2 incubator.
DNase I in 10 mL HBSS. The tissue digestion was performed at 37 °C for 15 min
with gently shaking once every 5 min. The reaction was stopped in 20 mL cold
MACS buffer containing 0.25% BSA, 2 mM EDTA, in PBS, and filtered with 100 Cultured ECs isolation from pig aorta. The aorta near to lung was separated and
μm and 40 μm cell strainer, and the cells were stained for the FC analysis. washed with PBS. Five milliliters of digestion buffer was injected into the sealed
aorta and left the aorta into a water bath at 37 °C for 30 min. Next, the aorta was
taken out, and the digestion buffer was flowed out. The medium was transferred
Retina. The porcine eyes were collected and dissected by surgical scissors, and the
into a 15 mL Falcon tube and centrifuged at 300 × g for 5 min. The collected cells
retina was transferred to cold HBSS by using forceps. Subsequently, the retina was
were washed with 5 mL PBS and centrifuged for culturing within 5 mL PAEC
dissected into small pieces and put in 10 mL digestion buffer with 2 mg/mL col-
medium.
lagenase I (Gibco, Cat#17018029), 3.75 μg/mL DNase I into 10 mL HBSS. The
tissue was digested at 37 °C for 20 min with gently shaking once every 2 min. The
digestion reaction was stopped in 20 mL MACS buffer with 0.25% BSA, 2 mM Single-cell library construction and sequencing. Library generation: Single-cell
EDTA in PBS and filtered with 40 μm cell strainer, and the cells were stained for RNA-sequencing (scRNA-seq) libraries were prepared following the manu-
the FC analysis. facturer’s instructions of GemCode Single Cell 3′ Gel Bead and Library kit (v3
Chemistry) from 10× Genomics, Inc. (Pleasanton, CA). Briefly, scRNA-seq library
Brain. The fresh porcine brain was collected including 7 different regions: neo- was generated by the cells from the different tissues: lung, liver, intestine, spleen,
cortex, cerebellar cortex, caudate nucleus, thalamus, hippocampus, hypothalamus, adipose, brain, and retina. After library construction, the library conversion was
and pons with 0.5 g pieces for each region. The samples from different regions were performed using the MGIEasy Universal DNA Library Preparation reagent kit
mixed and dissected into small pieces. The digestion was performed at 37 °C for (BGI, Shenzhen, China) for compatibility, followed by sequencing on a DNBSEQ-
30 mins in 20 mL digestion buffer containing 10 mg/mL collagenase IV, 15 μg/mL T7 platform (MGI).
DNase I, followed by gently shaking once every 5 min. Subsequently, the digestion
reaction was stopped by using MACS buffer and filtered with a 40 μm cell strainer.
After centrifuge and resuspension, the cells were stained for the FC analysis. Sample collection for nucleus extraction and snRNA-seq. The approaches for
sample collection and nuclei extraction for snRNA-seq libraries generation in this
Lung. Seven different regions: left apical lobe, left middle lobe, left main lobe, right study were followed in our previous study28. Briefly, the tissues used in this study:
apical lobe, right middle lobe, accessory lobe, and right main lobe were collected heart, kidney, spleen, liver, lung, retina, and brain regions (area_postrema, cere-
from the fresh porcine lung. From each region, a 0.5 g piece was collected and bellum, subfornical organ, and OVoLT) were carefully dissected from the healthy
washed twice in cold HBSS. The samples were dissected into small pieces and put domestic pigs (Sus scrofa domesticus, three-way hybrid of Landrace, Large White
into 50 mL Falcon tube with 20 mL digestion medium containing 1 mg/mL col- and Duroc, age 3 months) with strict compliance to the ethical guidelines. The
lagenase type II, 2.5 mg/mL collagenase type IV, 7.5 μg/mL DNase I. The samples collected tissues were washed with cold PBS, and immediately frozen in liquid
were digested at 37 °C for 30 min with gently shaking once every 5 min and nitrogen and stored in a −80 °C freezer before use. For the nuclei extraction
stopped digestion in MACS buffer. The samples were diluted into cold HBSS and process, the tissues were thawed and cut into small pieces, then transferred into a
then filtered with 100 μm and 40 μm cell strainers respectively. Subsequently, the homogenization buffer containing 20 mM Tris pH 8.0, 500 mM sucrose, 0.1% NP-
debris was removed, and the cells were stained for the FC analysis. 40, 0.2 U/mL RNase inhibitor, 1% BSA, and 0.1 mM DTT. The tissue pieces were
grinded for 15 times with tight pestles and filtered with 40 µm strainer. The
Visceral adipose tissue and subcutaneous adipose tissue. For each of the adipose samples were centrifuged at 500 × g for 10 min at 4 °C to carefully discard the
tissue depot, 10 g of each tissue was weighed and placed in the 25 mL digestion supernatant. The pellets (nuclei) were resuspended in PBS containing 1% BSA and
media containing: FBS-free KnockOut™ DMEM medium (Gibco, Cat#10829018) 20 U/µL RNase Inhibitor for later snRNA-seq library construction (MGI).
supplemented with 1% (v/v) Penicillin/Streptomycin (Thermo Fisher Scientific,
Cat#15140122), 2 × Antibiotic-Antimycotics (Thermo Fisher Scientific, Gibco
#1524-062), 1 mM Sodium Pyruvate (Thermo Fisher Scientific, Cat#11360070), Single-nuclei library construction and sequencing. The mRNA capture was
MEM Non-Essential Amino Acids Solution (MEM-NEAA) (Thermo Fisher Sci- operated on a DNBelab C4 device (MGI). cDNA amplification and libraries con-
entific, Cat#11140035), 2 mM L-Glutamine (Thermo Fisher Scientific, Gibco struction were generated using the MGI DNBelab C series reagent kit (MGI)
#25030-024), 0.2% collagenase type I, 0.25 U/mL Dispase II and 7.5 μg/mL DNase following the manufacturer’s instructions. All the libraries were sequenced on the
I. Samples were incubated at 37 °C in a water bath for at least 45 min, shaken and DNBSEQ-T7 platform.
mixed by pipetting every 5 min. At the end of the incubation time, a PBS-based
wash buffer containing 0.5% BSA, 2 mM EDTA in PBS was added and the cell
suspension was filtered through a 100 μm cell strainer, and centrifuged at 300 g for Pre-processing and quality control of scRNA-seq and snRNA-seq data. Cell
7 min. The washing step was repeated 2 times more and the cell suspension was Ranger 3.0.2 (10x Genomics) was used to process the raw sequencing data of
always re-filtered through 100 μm cell strainer and transferred to a new canonical scRNA-seq. The sequencing data from snRNA-seq was filtered, and the gene
tube to remove the access of undigested clumps and fat. expression matrix was obtained using DNBelab C Series scRNA analysis software
(MGI). The reference genome was downloaded from the Ensemble assembly:
Sscrofa11.1. Cells were only retained if the number of detected genes were greater
Intestine. The fresh intestine was collected and washed twice in cold PBS. 3 g of than 200 and less than 5000 and the percentage of detected mitochondrial tran-
intestine pieces were picked and dissected into small pieces. The digestion reaction scripts from MT genes (ATP6, ATP8, COX1, COX2, COX3, CYTB, ND2, ND3,
was at 30 min with gently shaking once every 5 min in 20 mL digestion medium ND4, ND4L, ND5, ND6) was less than 30%. The Pig ND1 gene was not included in
containing 5 mg/mL collagenase type I, 2.5 mg/mL collagenase type IV, 15 μg/mL MT-based filtering due to high sequence variant in pigs.
DNase I. After digestion, the tissues were diluted into cold HBSS and then filtered
with 100 μm and 40 μm cell strainers respectively. The debris was removed, and the
cells were stained for FC analysis. Identification of cell clusters. After filtering, unsupervised clustering was per-
formed using Seurat v332. Datasets from different sequencing libraries underwent
PBMC. The approaches for porcine PBMC isolation were based on the human normalizing and scaling. Variable genes were determined using Seurat’s “Find-
PBMC isolation protocol in our previous study102. In brief, 10 mL of pig blood VariableGenes” function with default parameters (selection.method = “vst”, nfea-
sample in a citrated vial was gently inverting to mix well. The PBMC were isolated tures = 2000). Clusters were identified via the “FindClusters” function (0.8 <
followed by density gradient centrifugation with Ficoll®-Paque Premium medium resolution < 1.5) implemented in Seurat using principal components with a P value
(Cytiva, Cat#17-5442-02). The cells were resuspended after red blood cells lysis and < 0.01 and subsequently visualized using the “RunTSNE” and “RunUMAP”
stained for the FC analysis. functions (reduction = “pca”).
Identification of differentially expressed genes (DEGs) across clusters. absolute alcohol (VWR chemicals) and xylene (Histolab). Next, paraffin (Histolab)
“FindAllMarkers” function implemented in Seurat v3 was used to identify DEGs immersion was performed using an automated Tissue Processing Center TPC 15
across clusters with the options “min.pct = 0.25, logfc.threshold = 0.25”. Multiple Duo (MEDITE) machine. After manual embedding into separate paraffin tissue
test correction for P value was performed using the Bonferroni method, and 0.05 blocks, one representative section (4um) was taken from each tissue block using a
was set as a threshold to define significance. Furthermore, cell-type identities were microtome (Microm HM 355 S, Thermo Fisher Scientific). A microm STS Section-
assigned using canonical cell-type markers. Transfer-System (waterfall) was used for section transfer into a warm water bath
(38 °C) stretching before placed on SuperFrost PlusTM slides (VWR). All slides
were dried at room temperature for 24 h followed by 50 °C overnight (LAMB
Gene ontology (GO) enrichment analysis. Gene Ontology (GO) analysis was
Windsor Incubator E18.31, Histolab).
using in the clusterProfiler 4.0 package103. The GO terms of selected genes were
enriched in the database “org.Hs.eg.db” using “enrichGO” function because of the
lack of study in pigs. Benjamini–Hochberg (BH) method was used for the multiple Immunohistochemical staining. Antibodies produced within the HPA project (see
test adjustment. Supplementary Data 14 for antibody information) with high reliability (based on
human antibody validation112,113), and over 80% sequence homology to pig
Integration of multiple datasets from different sequencing platforms. Datasets orthologs was used for immunostaining. All antibodies are published on the HPA
derived from the 10X (scRNA-seq) and the DNBelab C4 (snRNA-seq) platforms portal (www.proteinatlas.org) with more details about antibody reliability and
were integrated using the Seurat R package (version 3.2.2) after cell filtering. In tissue distribution in humans Formalin-fixed paraffin-embedded (FFPE) pig tis-
detailed, data was log1p-normalized with the SCT normalizeData method using the sues, previously validated and prepared30, was utilized for validating that scRNA-
“SCTransform” function, and subsequently scaled by the Pearson Residuals with a seq is a robust method for the classification of different cell types. The staining
scale factor of 10,000 as default using the “ScaleData” function. The top 3000 highly protocol follows our previous study by Karlsson et al, as briefly described below.
variable features were selected using the “SelectIntegrationFeatures” function, fol- Full size tissue sections as well as sections from tissue micro array (TMA) were
lowed by finding the integration anchors using the “FindIntegrationAnchors” represented by 1 mm punches moved from the donor block to an empty recipient
function, performing the integration of the data using the “IntegrateData” function. paraffin block. Deparaffinization and rehydration were performed by Autostainer
Following integration, principal component analysis was performed using the XL (ST5010, Leica biosystems) followed by heat-induced epitope retrieval in pH6.1
“RunPCA” function with default parameters, then both t-SNE (t-Distributed Sto- citrate buffer (DAKO, diluted 1:10 with deionized water) and pressure boiler
chastic Neighbor Embedding) and UMAP (Uniform Manifold Approximation and (decloaking chamber, Biocare Medical). Autostainer 480 (ThermoFisher Scientific)
Projection) dimensionality reduction methods were conducted based on the top 20 was used for automated IHC staining with UltraVision™ Quanto Detection System
principal components (PCs) using the “RunTSNE” and “RunUMAP” functions, HRP DAB-kit from Thermo Fisher Scientific including Ultra V Block, HRP
respectively. Moreover, unsupervised clusters were identified by setting the top 20 Polymer, Primary Antibody Enhancer, DAB Quanto Substrate, DAB Quanto
PCs and a clustering resolution of 1.0 using “FindNeighbors” and “FindClusters” Chromogen, and primary antibodies were diluted in Antibody Diluent OP Quanto.
functions. The Pearson correlation coefficients of cell types were calculated using After a 5 min block (Ultra V Block), primary antibodies were incubated for 30 min,
the average expression of top 3000 highly variable features and visualized using followed by the secondary HRP Polymer (ThermoFisher Scientific) and the final
pheatmap R package (1.0.12). step of 5 min DAB Quanto incubation. All steps were separated by double wash
(Tris Buffer and Tween 20). Slides were placed in water and moved to the Auto-
stainer XL (ST5010, Leica biosystems) for counterstaining (hematoxyline), dehy-
Pseudotime trajectory analysis. PECAM1 + PTPRC-EPCAM-COL1A1-PDGFRB- dration, and cover glass mounting. Image digitalization was performed with
HBBlow ECs was used to subset the ECs from adipose tissues. Then, the “Inte- Scanscope AT2 (Aperio) using a 20× objective.
grateData” function with cca methods in Seurat package were used to correct the
batch effect of ECs in adipose-V and adipose-S. The integrated ECs of both
adipose-V and adipose-S were clustered by “FindClusters” with resolution 1.0. Each Immunofluorescence staining. Two slides with adipose tissue samples were
cluster was annotated based on canonical subtype markers and cluster specific dewaxed using Histoclear. One slide of adipose tissue was then incubated with
DEGs. Both all ECs and EndMT subtype were used to perform the subsequent ACTA2 monoclonal mouse antibody (1:100 dilution) and PECAM1 polyclonal
pseudotime trajectory analysis. Monocle 2104 package was used to discover the cell rabbit antibody (1:50 dilution) overnight at 4 °C. The other slide was incubated
state transitions ECs. Genes expressed in less than 10 cells were filtered out. DEGs with VWF monoclonal mouse antibody (1:7,000 dilution) and TAGLN polyclonal
were computed by function “differentialGeneTest” in monocle2. Genes with qvalue rabbit antibody (1:100 dilution). The slides were then washed in TBS-Tween and
less than 0.01 were regarded as DEGs and sorted by qvalue using “setOrder- blocked with TNB buffer. To visualize the proteins, an HRP-conjugated swine anti-
ingFilter” function. The pseudotime trajectory was constructed by “DDRTree” rabbit antibody (diluted 1:200 in TNB buffer) was added and the slides were
algorithm with default parameters. The dynamical expression changes of selected incubated at room temperature for 30 min. The slides were washed and then
marker genes by pseudotime were visualized by “plot_genes_in_pseudotime” and incubated with TSA-FITC substrate for 15 min in the dark. After this step, the
“plot_pseudotime_heatmap” function. slides were placed in 0.1% NaN3 in PBS to inactivate the secondary antibody. They
were then incubated in the dark with an HRP-conjugated donkey anti-mouse
antibody (diluted 1:200 in TNB buffer), and after washing incubated with a TSA-
Cross-species comparison of intercellular communications between pig and Cy3.5 substrate.
human. Intercellular communication analysis was conducted using CellChat
(v0.0.1) R package with default parameters65. Pig and human liver, kidney, and
heart datasets were analyzed separately. The human liver105, kidney106, and TGF-β2 treatment-induced EndMT. Primary pig aortic endothelial cells (PAECs)
heart107 datasets were collected from previous reports. Intercellular communica- were cultured in Dulbecco’s modified Eagles’s medium (DMEM) with 4.5 g/L D-
tions analysis was performed based on cell types shared by the separate liver, glucose, L-Glutamine, and pyruvate (Gibco, #41966052), supplemented with 10%
kidney, and heart between pig and human. Cell–cell communication network was fetal bovine serum (FBS) (Sigma, #F7524), 1× non-essentials AA (MERCK,
visualized using the “netVisual_aggregate” function, centrality score was computed #M7145), and 1× penicillin/streptomycin. The PAECs were cultured on 0.1%
and visualized using the “netAnalysis_signalingRole_network” function, relative gelatin (MERCK, #D8537) coated culturing flasks (Nunc, #156499) in a 5% CO2
contribution of each ligand-receptor pair was visualized using the “netAnaly- humidified incubator at 37 °C.
sis_contribution” function. Primary human umbilical vein endothelial cells (HUVECs) were cultured on
0.1% gelatin-coated culturing flasks in M199 medium (Gibco, #22350-029)
supplemented with 2mM L-Glutamine (Gibco, #35050-061), Endothelial Cell
TF-target interaction inference. TFs gene list was downloaded from the Growth Supplement (ECGS)/ Heparin (PromoCell, #C-30120), and 20% FBS, in a
animalTFDB3.0108. Only genes expressed in more than 5% of corresponding cell 5% CO2 humidified incubator at 37 °C. The M199 medium with supplements was
types were subjected to GENIE3109 (v1.8.0) to infer putative regulatory circuits replaced three times per week. HUVECs were passaged every 14 days by washing
from expression data using tree-based ensemble methods. TF-target pairs with the cells twice with PBS, trypsinized, and centrifuged at 250 × g for 5 min.
weight value more than 0.01 were retained for primary regulatory network con- For TGF-β2 treatment, 100,000 PAECs/well were seeded in 12-well plates
struction. The total frequency of each TF-target pair present in 13 species was (Thermo Scientific, #150628) supplement with 2 ng/mL TGF-β2(Merck, T2815) in
counted to evaluate its conservation level. TF-target pairs present in at least triplicates. Control cells were cultured in normal PAECs medium without TGF-β2.
5 species were considered as ‘conserved regulomes’. Conserved regulomes were On day 2, the PAECs were trypsinized and seeded (100,000 PAECs/well) in 12-well
visualized using the igraph110 (v1.2.6) R package. plates followed by changing the medium on day 4. For HUVECs, which grow much
slower compared to PAECs, cells were cultured in control (0 ng/mL, TGF-β2) and
Construction of the PCA database. The PCA database was generated with TGF-β2 medium (2 ng/mL) for 5 days without passaging. All cells were harvested
ShinyCell111 with default settings and modified to include the introduction and for analysis of ACTA2 and CD31 expression by flow cytometry 5 days after TGF-
user guide pages. β2 treatment.
Tissue processing. All peripheral pig tissues were stored in 70% ethanol at 4 °C. Flow cytometry analysis of ACTA2 and CD31 expression. Cells (PAECs and
The pig brain tissues were stored in PBS buffer at 4 °C and changed into 70% HUVECs) were washed in PBS twice, dissociated with trypsin (Gibco, #25300054),
ethanol one week prior to paraffin embedding. We first dehydrated the tissues with and centrifuged at 400 g for 4 min. The cell pellets were washed with PBS + 5% FBS
twice and resuspended in 240 µL PBS + 5% FBS. PAECs were stained with ACTA2 9. Cooper, D. K., Ekser, B., Ramsoondar, J., Phelps, C. & Ayares, D. The role of
(AF488 Human Alpha-Smooth muscle actin, R&D, Cat#IC1420G, 1:100 dilution), genetically engineered pigs in xenotransplantation research. J. Pathol. 238,
CD31 (Porcine CD31/PECAM1 AF700, R&D, #FAB33871N-100UG, 1:100 dilu- 288–299 (2016).
tion), and CD45 (BV421 Mouse Anti-Human CD45, BD Bioscience, #563879, 10. Bustad, L. K. & McClellan, R. O. Swine in biomedical research. Science 152,
1:100 dilution). HUVECs were stained with CD31 (FITC Mouse Anti-Human 1526–1530 (1966).
CD31 (BD Bioscience, #555445), 1:100 dilution) and ACTA2 (AF488 Human 11. Lelovas, P. P., Kostomitsopoulos, N. G. & Xanthos, T. T. A comparative
alpha-Smooth muscle actin, R&D, Cat#IC1420G, 1:100 dilution) separately. Cells anatomic and physiologic overview of the porcine heart. J. Am. Assoc. Lab
were incubated on ice in the dark for 30 min. Before FC analysis, the stained cells Anim. Sci. 53, 432–438 (2014).
were washed twice with pre-cooled PBS + 5% FBS. At the final wash, the cells were 12. Porrett, P. M. et al. First clinical-grade porcine kidney xenotransplant using a
resuspended in 200 µL PBS + 5% FBS. Prior to this study, we compared the pig and human decedent model. Am. J. Transplant. (2022).
human ACTA2 amino acid sequences, which are identical. Antibody validation of 13. Liu, T. et al. Cross-species single-cell transcriptomic analysis reveals pre-
anti-ACTA2 was also validated by staining human mesenchymal stem cells. Flow gastrulation developmental differences among pigs, monkeys, and humans.
cytometry was performed with the NovoCyte Quanteon analyzer (Acea bioscience, Cell Discov. 7, 8 (2021).
Inc, US) provided by the FACS Core Facility, Aarhus University. 14. Gladka, M. M. et al. Single-cell sequencing of the healthy and diseased heart
reveals cytoskeleton-associated protein 4 as a new modulator of fibroblasts
Statistics and reproducibility. Statistical significance of differential expression activation. Circulation 138, 166–180 (2018).
gene was performed with multiple test correction for P value (Bonferroni). 15. Chen, G., Ning, B. & Shi, T. Single-cell RNA-seq technologies and related
Benjamini–Hochberg (BH) method was used for the multiple test adjustment for computational data analysis. Front. Genet. 10, 317 (2019).
the gene ontology enrichment analysis. Unpaired, two-sided t-test was used for 16. Packer, J. S. et al. A lineage-resolved molecular atlas of C. elegans
comparison of CD31 and ACTA2 expression in the induced EndMT experiments embryogenesis at single-cell resolution. Science 365, eaax1971 (2019).
in cultured ECs. A P value less than 0.05 was considered statistically significant. 17. Plass, M. et al. Cell type atlas and lineage tree of a whole complex animal by
Unless stated elsewhere, all immunohistochemistry, immunofluorescence staining single-cell transcriptomics. Science 360, eaaq1723 (2018).
and FACS experiments was performed with at least three experimental replicates. 18. Davie, K. et al. A single-cell transcriptome atlas of the aging Drosophila brain.
No statistical method was used to predetermine sample size. No data were excluded Cell 174, 982–998 e920 (2018).
from the analyses. Filtering criteria for the low-quality cells and potential doublets 19. Wagner, D. E. et al. Single-cell mapping of gene expression landscapes and
are provided in the method above and the source code. Randomization is not lineage in the zebrafish embryo. Science 360, 981–987 (2018).
related to this study. The investigators were not blinded to allocation during the 20. Tabula Muris, C. et al. Single-cell transcriptomics of 20 mouse organs creates a
experiments and the outcome assessment. Tabula Muris. Nature 562, 367–372 (2018).
21. Wang, S. et al. Single-cell transcriptomic atlas of primate ovarian aging. Cell
Reporting summary. Further information on research design is available in the Nature 180, 585–600 e519 (2020).
Research Reporting Summary linked to this article. 22. Cao, J. et al. A human cell atlas of fetal gene expression. Science 370, eaba7721
(2020).
23. Groenen, M. A. et al. Analyses of pig genomes provide insight into porcine
demography and evolution. Nature 491, 393–398 (2012).
Data availability 24. Jinek, M. et al. A programmable dual-RNA-guided DNA endonuclease in
The single-cell and single-nuclei RNA-sequencing data generated in this study have been adaptive bacterial immunity. Science 337, 816–821 (2012).
deposited in the CNGB Sequence Archive (CNSA) of China National GeneBank 25. Du, Y. et al. Piglets born from handmade cloning, an innovative cloning
DataBase (CNGBdb) under accession code “CNP0002165”. The single-cell and single- method without micromanipulation. Theriogenology 68, 1104–1110 (2007).
nuclei RNA-sequencing data generated in this study have also been deposited in the gene 26. Zhang, L. et al. A high-resolution cell atlas of the domestic pig lung and an
expression omnibus database (GEO) under accession codes “GSE196055” and online platform for exploring lung single-cell data. J Genet Genomics 48,
“GSE193975 [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi]”. All matrix data can be 411–425 (2021).
downloaded from the PCA database (https://dreamapp.biomed.au.dk/pigatlas/). All other 27. Herrera-Uribe, J. et al. Reference Transcriptomes of porcine peripheral
relevant data supporting the key findings of this study are available within the article and immune cells created through bulk and single-cell RNA sequencing. Front
its Supplementary Information files or from the corresponding author upon reasonable Genet. 12, 689406 (2021).
request. Source data are provided with this paper. 28. Zhu, J. et al. Single-cell atlas of domestic pig cerebral cortex and
hypothalamus. Sci. Bull. 66, 1448–1461 (2021).
29. Jin, L. et al. A pig BodyMap transcriptome reveals diverse tissue physiologies
and evolutionary dynamics of transcription. Nat. Commun. 12, 3715 (2021).
Code availability 30. Karlsson, M. et al. Genome-wide annotation of protein-coding genes in pig.
All codes for processing of the single-cell RNA-sequencing data have been deposited to
BMC Biol. 20, 25 (2022).
GitHub114 and available through this URL: https://github.com/Dingpw/PigAtlas.
31. Liu, F. et al. Systematic comparative analysis of single-nucleotide variant
detection methods from single-cell RNA sequencing data. Genome Biol. 20,
Received: 24 November 2021; Accepted: 16 June 2022; 242 (2019).
Published online: 24 June 2022 32. Satija, R., Farrell, J. A., Gennert, D., Schier, A. F. & Regev, A. Spatial
reconstruction of single-cell gene expression data. Nat. Biotechnol. 33,
495–502 (2015).
33. Sonawane, A. R. et al. Understanding tissue-specific gene regulation. Cell Rep.
References
21, 1077–1088 (2017).
1. Rogers, C. S. et al. Disruption of the CFTR gene produces a model of cystic
34. Han, X. et al. Construction of a human cell landscape at single-cell level.
fibrosis in newborn pigs. Science 321, 1837–1841 (2008).
2. Al-Mashhadi, R. H. et al. Familial hypercholesterolemia and atherosclerosis in Nature 581, 303–309 (2020).
35. Middleton, S. Porcine ophthalmology. Vet. Clin. N. Am. Food Anim. Pract. 26,
cloned minipigs created by DNA transposition of a human PCSK9 gain-of-
557–572 (2010).
function mutant. Sci. Transl. Med. 5, 166ra161 (2013).
36. Choi, K. E. et al. An experimental pig model with outer retinal degeneration
3. Yan, S. et al. A huntingtin knockin pig model recapitulates features of selective
neurodegeneration in Huntington’s disease. Cell 173, 989–1002 e1013 (2018). induced by temporary intravitreal loading of N-methyl-N-nitrosourea during
vitrectomy. Sci. Rep. 11, 258 (2021).
4. Prather, R. S., Lorson, M., Ross, J. W., Whyte, J. J. & Walters, E. Genetically
37. Chade, A. R., Williams, M. L., Engel, J., Guise, E. & Harvey, T. W. A
engineered pig models for human diseases. Annu. Rev. Anim. Biosci. 1,
translational model of chronic kidney disease in swine. Am. J. Physiol. Ren.
203–219 (2013).
Physiol. 315, F364–F373 (2018).
5. Walters, E. M. & Prather, R. S. Advancing swine models for human health and
38. Giraud, S. et al. Contribution of large pig for renal ischemia-reperfusion and
diseases. Mo. Med. 110, 212–215 (2013).
transplantation studies: the preclinical model. J. Biomed. Biotechnol. 2011,
6. Kleinwort, K. J. H. et al. Retinopathy with central oedema in an INS (C94Y)
532127 (2011).
transgenic pig model of long-term diabetes. Diabetologia 60, 1541–1549
39. Misra, S. et al. The porcine remnant kidney model of chronic renal
(2017).
7. Wolf, E., Kemter, E., Klymiuk, N. & Reichart, B. Genetically modified pigs as insufficiency. J. Surg. Res. 135, 370–379 (2006).
40. Schwartz, G. & Berry, M. J. 2nd Sophisticated temporal pattern recognition in
donors of cells, tissues, and organs for xenotransplantation. Anim. Front. 9,
retinal ganglion cells. J. Neurophysiol. 99, 1787–1798 (2008).
13–20 (2019).
41. Vlasiuk, A. & Asari, H. Feedback from retinal ganglion cells to the inner
8. Sykes, M. & Sachs, D.H. Transplanting organs from pigs to humans. Sci.
Immunol. 4, eaau6298 (2019). retina. PLoS ONE 16, e0254611 (2021).
42. Kandori, H. Biophysics of rhodopsins and optogenetics. Biophys. Rev. 12, 73. Geirsdottir, L. et al. Cross-species single-cell analysis reveals divergence of the
355–361 (2020). primate microglia program. Cell 181, 746 (2020).
43. Shekhar, K. et al. Comprehensive classification of retinal bipolar neurons by 74. Holtman, I. R., Skola, D. & Glass, C. K. Transcriptional control of microglia
single-cell transcriptomics. Cell 166, 1308–1323 e1330 (2016). phenotypes in health and disease. J. Clin. Investig. 127, 3220–3229 (2017).
44. Pannabecker, T. L. Structure and function of the thin limbs of the loop of 75. Walker, D. G. et al. Patterns of expression of purinergic receptor P2RY12, a
Henle. Compr. Physiol. 2, 2063–2086 (2012). putative marker for non-activated microglia, in aged and Alzheimer’s disease
45. Kalucka, J. et al. Single-cell transcriptome atlas of murine endothelial cells. Cell brains. Int. J. Mol. Sci. 21, 678 (2020).
180, 764–779 e720 (2020). 76. Weng, Q. et al. Single-cell transcriptomics uncovers glial progenitor diversity
46. Goveia, J. et al. An integrated gene expression landscape profiling approach to and cell fate determinants during development and gliomagenesis. Cell Stem
identify lung tumor endothelial cell heterogeneity and angiogenic candidates. Cell 24, 707–723 e708 (2019).
Cancer Cell 37, 21–36 e13 (2020). 77. Lee, J. K. et al. Regulator of G-protein signaling 10 promotes dopaminergic
47. Rohlenova, K. et al. Single-cell RNA sequencing maps endothelial metabolic neuron survival via regulation of the microglial inflammatory response. J.
plasticity in pathological angiogenesis. Cell Metab. 31, 862–877 e814 (2020). Neurosci. 28, 8517–8528 (2008).
48. Teuwen, L. A. et al. Tumor vessel co-option probed by single-cell analysis. Cell 78. Oishi, Y. et al. SREBP1 contributes to resolution of pro-inflammatory
Rep. 35, 109253 (2021). TLR4 signaling by reprogramming fatty acid metabolism. Cell Metab. 25,
49. Dumas, S. J. et al. Single-cell RNA sequencing reveals renal endothelium 412–427 (2017).
heterogeneity and metabolic adaptation to water deprivation. J. Am. Soc. 79. Deczkowska, A. et al. Mef2C restrains microglial inflammatory response and is
Nephrol. 31, 118–138 (2020). lost in brain ageing in an IFN-I-dependent manner. Nat. Commun. 8, 717 (2017).
50. Gottfried, E. et al. Expression of CD68 in non-myeloid cell types. Scand. J. 80. Han, X. et al. Mapping the mouse cell atlas by microwell-seq. Cell 172,
Immunol. 67, 453–463 (2008). 1091–1107 e1017 (2018).
51. Bulla, R. et al. C1q acts in the tumour microenvironment as a cancer- 81. Tabula Muris, C. A single-cell transcriptomic atlas characterizes ageing tissues
promoting factor independently of complement activation. Nat. Commun. 7, in the mouse. Nature 583, 590–595 (2020).
10346 (2016). 82. Al-Soudi, A., Kaaij, M. H. & Tas, S. W. Endothelial cells: from innocent
52. Tian, Y., Jain, S., Kelemen, S. E. & Autieri, M. V. AIF-1 expression regulates bystanders to active participants in immune responses. Autoimmun. Rev. 16,
endothelial cell activation, signal transduction, and vasculogenesis. Am. J. 951–962 (2017).
Physiol. Cell Physiol. 296, C256–C266 (2009). 83. Danese, S., Dejana, E. & Fiocchi, C. Immune regulation by microvascular
53. He, H. et al. Endothelial cells provide an instructive niche for the endothelial cells: directing innate and adaptive immunity, coagulation, and
differentiation and functional polarization of M2-like macrophages. Blood inflammation. J. Immunol. 178, 6017–6022 (2007).
120, 3152–3162 (2012). 84. Mai, J., Virtue, A., Shen, J., Wang, H. & Yang, X. F. An evolving new
54. Piera-Velazquez, S. & Jimenez, S. A. Endothelial to mesenchymal transition: paradigm: endothelial cells–conditional innate immune cells. J. Hematol.
role in physiology and in the pathogenesis of human diseases. Physiol. Rev. 99, Oncol. 6, 61 (2013).
1281–1324 (2019). 85. Dunleavey, J. M. & Dudley, A. C. Vascular Mimicry: concepts and
55. Ma, J., Sanchez-Duffhues, G., Goumans, M. J. & Ten Dijke, P. TGF-beta- implications for anti-angiogenic therapy. Curr. Angiogenes 1, 133–138 (2012).
induced endothelial to mesenchymal transition in disease and tissue 86. Barnett, F. H. et al. Macrophages form functional vascular mimicry channels
engineering. Front Cell Dev. Biol. 8, 260 (2020). in vivo. Sci. Rep. 6, 36659 (2016).
56. Zeisberg, E. M. et al. Endothelial-to-mesenchymal transition contributes to 87. Debels, H. et al. Macrophages play a key role in angiogenesis and adipogenesis
cardiac fibrosis. Nat. Med 13, 952–961 (2007). in a mouse tissue engineering model. Tissue Eng. Part A 19, 2615–2625
57. Dejana, E., Hirschi, K. K. & Simons, M. The molecular basis of endothelial cell (2013).
plasticity. Nat. Commun. 8, 14361 (2017). 88. Riabov, V. et al. Role of tumor associated macrophages in tumor angiogenesis
58. Pardali, E., Sanchez-Duffhues, G., Gomez-Puerto, M. C. & Ten Dijke, P. TGF- and lymphangiogenesis. Front Physiol. 5, 75 (2014).
beta-induced endothelial-mesenchymal transition in fibrotic diseases. Int. J. 89. Liu, L. & Shi, G. P. CD31: beyond a marker for endothelial cells. Cardiovasc.
Mol. Sci. 18, 2157 (2017). Res. 94, 3–5 (2012).
59. Bischoff, J. Endothelial-to-mesenchymal transition. Circ. Res. 124, 1163–1165 90. Potenta, S., Zeisberg, E. & Kalluri, R. The role of endothelial-to-mesenchymal
(2019). transition in cancer progression. Br. J. Cancer 99, 1375–1379 (2008).
60. Cooley, B. C. et al. TGF-beta signaling mediates endothelial-to-mesenchymal 91. Man, S., Sanchez Duffhues, G., Ten Dijke, P. & Baker, D. The therapeutic
transition (EndMT) during vein graft remodeling. Sci. Transl. Med. 6, potential of targeting the endothelial-to-mesenchymal transition. Angiogenesis
227ra234 (2014). 22, 3–13 (2019).
61. Chen, P. Y. et al. Endothelial-to-mesenchymal transition drives atherosclerosis 92. Carmeliet, P. Angiogenesis in health and disease. Nat. Med. 9, 653–660 (2003).
progression. J. Clin. Investig. 125, 4514–4528 (2015). 93. Chiaverina, G. et al. Dynamic interplay between pericytes and endothelial cells
62. Fuseya, T. et al. Ectopic fatty acid-binding protein 4 expression in the vascular during sprouting angiogenesis. Cells 8, 1109 (2019).
endothelium is involved in neointima formation after vascular injury. J. Am. 94. Ramirez-Pedraza, M. & Fernandez, M. Interplay between macrophages and
Heart Assoc. 6, e006377 (2017). angiogenesis: a double-edged sword in liver disease. Front. Immunol. 10, 2882
63. Zaniboni, A. et al. In vitro differentiation of porcine aortic vascular precursor (2019).
cells to endothelial and vascular smooth muscle cells. Am. J. Physiol. Cell 95. Falkenberg, K. D., Rohlenova, K., Luo, Y. & Carmeliet, P. The metabolic
Physiol. 309, C320–C331 (2015). engine of endothelial cells. Nat. Metab. 1, 937–946 (2019).
64. Red-Horse, K., Crawford, Y., Shojaei, F. & Ferrara, N. Endothelium- 96. Nishida, N., Yano, H., Nishida, T., Kamura, T. & Kojiro, M. Angiogenesis in
microenvironment interactions in the developing embryo and in the adult. cancer. Vasc. Health Risk Manag 2, 213–219 (2006).
Dev. Cell 12, 181–194 (2007). 97. Fallah, A. et al. Therapeutic targeting of angiogenesis molecular pathways in
65. Jin, S. et al. Inference and analysis of cell-cell communication using CellChat. angiogenesis-dependent diseases. Biomed. Pharmacother. 110, 775–785
Nat. Commun. 12, 1088 (2021). (2019).
66. Swindle, M. M., Makin, A., Herron, A. J., Clubb, F. J. Jr. & Frazier, K. S. Swine 98. Hammond, B. P., Manek, R., Kerr, B. J., Macauley, M. S. & Plemel, J. R.
as models in biomedical research and toxicology testing. Vet. Pathol. 49, Regulation of microglia population dynamics throughout development,
344–356 (2012). health, and disease. Glia 69, 2771–2797 (2021).
67. Abhinand, C. S., Raju, R., Soumya, S. J., Arya, P. S. & Sudhakaran, P. R. 99. Vaquerizas, J. M., Kummerfeld, S. K., Teichmann, S. A. & Luscombe, N. M. A
VEGF-A/VEGFR2 signaling network in endothelial cells relevant to census of human transcription factors: function, expression and evolution.
angiogenesis. J. Cell Commun. Signal 10, 347–354 (2016). Nat. Rev. Genet. 10, 252–263 (2009).
68. Erber, R. et al. Combined inhibition of VEGF and PDGF signaling enforces 100. Hemberg, M. & Kreiman, G. Conservation of transcription factor binding
tumor vessel regression by interfering with pericyte-mediated endothelial cell events predicts gene expression across species. Nucleic Acids Res. 39,
survival mechanisms. FASEB J. 18, 338–340 (2004). 7092–7102 (2011).
69. Goumans, M. J., Liu, Z. & ten Dijke, P. TGF-beta signaling in vascular biology 101. Nitta, K. R. et al. Conservation of transcription factor binding specificities
and dysfunction. Cell Res. 19, 116–127 (2009). across 600 million years of bilateria evolution. Elife 4, e04837 (2015).
70. Dyer, L. A., Pi, X. & Patterson, C. The role of BMPs in endothelial cell 102. Wang, F. et al. A single-cell approach to engineer CD8+ T cells targeting
function and dysfunction. Trends Endocrinol. Metab. 25, 472–480 (2014). cytomegalovirus. Cell Mol. Immunol. 18, 1326–1328 (2021).
71. Gold, D. A. et al. RORalpha coordinates reciprocal signaling in cerebellar 103. Wu, T. et al. clusterProfiler 4.0: a universal enrichment tool for interpreting
development through sonic hedgehog and calcium-dependent pathways. omics data. Innovation 2, 100141 (2021).
Neuron 40, 1119–1131 (2003). 104. Trapnell, C. et al. The dynamics and regulators of cell fate decisions are
72. Li, Q. & Barres, B. A. Microglia and macrophages in brain homeostasis and revealed by pseudotemporal ordering of single cells. Nat. Biotechnol. 32,
disease. Nat. Rev. Immunol. 18, 225–242 (2018). 381–386 (2014).
105. MacParland, S. A. et al. Single cell RNA sequencing of human liver reveals Professor Fred Dubee for providing critical comments and language editing for the
distinct intrahepatic macrophage populations. Nat. Commun. 9, 4383 (2018). manuscript.
106. Lake, B. B. et al. A single-nucleus RNA-sequencing pipeline to decipher the
molecular anatomy and pathophysiology of human kidneys. Nat. Commun.
10, 2832 (2019).
Author contributions
Conceptualization, Y.L., D.C., L.L., M.U., L.B., and H.Y.; Methodology, F.W., P.D., X.L.,
107. Litvinukova, M. et al. Cells of the adult human heart. Nature 588, 466–472
X.D., C.B.B, and E.S.; Investigation, F.W. (major), P.D. (major), X.L., X.D., C.B.B, E.S.,
(2020).
S.B., J.Z., L.Z., L.P.D.R., Li.L, Y.W., W.Z., Z.L., J.H., R.L., Q.Q., Y.J., W.W., Y.Y., M.P.,
108. Hu, H. et al. AnimalTFDB 3.0: a comprehensive resource for annotation and
H.W., A.W., Jac.H., J.K., and M.K.; Sequencing, L.X., P.L., F.C., and H.J.; Writing –
prediction of animal transcription factors. Nucleic Acids Res. 47, D33–D38
Original Draft, F.W. and Y.L.; Writing – review & editing, all authors; Funding acqui-
(2019).
sition, M.U., L.B., H.Y., X.X., D.C., L.L., and Y.L.; Resources, X.Z., R.B.H., L.F., C.L., F.P.,
109. Huynh-Thu, V. A., Irrthum, A., Wehenkel, L. & Geurts, P. Inferring
M.U., L.B., N.J., X.X., H.Y., P.C., J.M., D.C., L.L., and Y.L.; Supervision, X.Z., R.B.H., L.F.,
regulatory networks from expression data using tree-based methods. PLoS
C.L., F.P., M.U., L.B., N.J., X.X., H.Y., P.C., J.M., D.C., L.L., and Y.L.
ONE 5, e12776 (2010).
110. Mora, A. & Donaldson, I. M. iRefR: an R package to manipulate the iRefIndex
consolidated protein interaction database. BMC Bioinform. 12, 455 (2011). Competing interests
111. Ouyang, J. F., Kamaraj, U. S., Cao, E. Y. & Rackham, O. J. L. ShinyCell: simple The authors declare no competing interests.
and sharable visualisation of single-cell gene expression data. Bioinformatics
37, 3374–3376 (2021).
112. Edfors, F. et al. Enhanced validation of antibodies for research applications. Additional information
Nat. Commun. 9, 4130 (2018). Supplementary information The online version contains supplementary material
113. Sivertsson, A. et al. Enhanced validation of antibodies enables the discovery of available at https://doi.org/10.1038/s41467-022-31388-z.
missing proteins. J. Proteome Res. 19, 4766–4781 (2020).
114. Wang, F., Ding, P. & Luo, Y. Endothelial cell heterogeneity and microglia Correspondence and requests for materials should be addressed to Dongsheng Chen,
regulons revealed by a pig cell landscape at single-cell level. Zenodo. https:// Lin Lin or Yonglun Luo.
doi.org/10.5281/zenodo.6607859 (2022).
Peer review information Nature Communications thanks Nicole Valenzuela,
Christopher Tuggle, Lance Daharsh and the other, anonymous, reviewer(s) for their
Acknowledgements contribution to the peer review of this work. Peer reviewer reports are available
We thank the FACS core team from the Department of Biomedicine at Aarhus Uni-
versity for the help of FACS experiments and data analysis. We thank Trine Skov Reprints and permission information is available at http://www.nature.com/reprints
Petersen and Xi Xiang for their technical assistance. This project was partially supported
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in
by the Sanming Project of Medicine in Shenzhen (SZSM201612074, to L.B. and Y.L),
published maps and institutional affiliations.
Guangdong Provincial Key Laboratory of Genome Read and Write (No.
2017B030301011 to X.Xu.), Guangdong Provincial Academician Workstation of BGI
Synthetic Genomics (No. 2017B090904014 to X.Xu.), CAMS Innovation Fund for
Medical Sciences (CIFMS) (22021-I2M-1-061 to D.C.), the DFF Sapere Aude Starting Open Access This article is licensed under a Creative Commons
grant (8048-00072A to L.L.), and the Novo Nordisk Foundation (NNF21OC0071718 to Attribution 4.0 International License, which permits use, sharing,
L.L.). J.K. is supported by AIAS-CO-FUND II: GA: MSCA:754513, Lundbeckfonden: adaptation, distribution and reproduction in any medium or format, as long as you give
R307-2018-3667, Carlsberg Fonden: CF19-0687, Kræftens Bekæmpelse: R302-A17296, appropriate credit to the original author(s) and the source, provide a link to the Creative
A.P. Møller Fonden: 20-L-0317, Riisfort Fonden and Steno Diabetes Center Aarhus Commons license, and indicate if changes were made. The images or other third party
(SDCA). P.C. is supported by Grants from Methusalem funding (Flemish government), material in this article are included in the article’s Creative Commons license, unless
the Fund for Scientific Research-Flanders (FWO-Vlaanderen), the European Research indicated otherwise in a credit line to the material. If material is not included in the
Council ERC Advanced Research Grant EU- ERC74307 and the NNF Laureate Research article’s Creative Commons license and your intended use is not permitted by statutory
Grant from Novo Nordisk Foundation (Denmark). Y.L. is supported by the Independent regulation or exceeds the permitted use, you will need to obtain permission directly from
Research Fund Denmark (9041-00317B), European Union’s Horizon 2020 research and the copyright holder. To view a copy of this license, visit http://creativecommons.org/
innovation program under grant agreement No 899417, and the Novo Nordisk Foun- licenses/by/4.0/.
dation (NNF21OC0068988; NNF21OC0071031). We thank Hårby Slagteren IvS for
providing the study materials. We thank the China National GeneBank for the support of
executing the project under the framework of Genome Read and Write. We thank © The Author(s) 2022, corrected publication 2022
1
Lars Bolund Institute of Regenerative Medicine, Qingdao-Europe Advanced Institute for Life Sciences, BGI-Qingdao, BGI-Shenzhen,
Qingdao, China. 2Department of Biomedicine, Aarhus University, Aarhus, Denmark. 3BGI-Shenzhen, Shenzhen, China. 4College of Life Sciences,
University of Chinese Academy of Sciences, Beijing, China. 5Department of Biology, University of Copenhagen, Copenhagen, Denmark. 6Steno
Diabetes Center Aarhus, Aarhus University Hospital, Aarhus, Denmark. 7Department of Neuroscience, Karolinska Institutet, Stockholm, Sweden.
8
MGI, BGI-Shenzhen, Shenzhen, China. 9Laboratory of Angiogenesis and Vascular Metabolism, Center for Cancer Biology, VIB, Leuven, Belgium.
10
Department of Oncology, Leuven Cancer Institute, KU Leuven, Leuven, Belgium. 11College of Basic Medicine, Qingdao University,
Qingdao, China. 12School of Basic Medical Sciences, Binzhou Medical University, Yantai, China. 13Institute of Systems Medicine, Chinese
Academy of Medical Sciences, Peking Union Medical College, Beijing, China. 14Suzhou Institute of Systems Medicine, Suzhou, China. 15Aarhus
University of Advanced Studies (AIAS), Aarhus University, Aarhus, Denmark. 16Department of Protein Science, Science for Life Laboratory, KTH-
Royal Institute of Technology, Stockholm, Sweden. 17Department of Obstetrics and Gynecology, Aarhus University Hospital, Aarhus, Denmark.
18
Department of Immunology, Genetics and Pathology, Uppsala University, Uppsala, Sweden. 19IBMC-BGI Center, the Cancer Hospital of the
University of Chinese Academy of Sciences (Zhejiang Cancer Hospital), Institute of Basic Medicine and Cancer (IBMC), Chinese Academy of
Sciences, Hangzhou, China. 20Center for Biotechnology, Khalifa University of Science and Technology, Abu Dhabi, United Arab Emirates. 21These
authors contributed equally: Fei Wang, Peiwen Ding, Xue Liang, Xiangning Ding, Camilla Blunk Brandt. ✉email: cds@ism.pumc.edu.cn;
lin.lin@biomed.au.dk; alun@biomed.au.dk