Coronavirus in the Skin of the Pangolin: Implications for SARS-CoV-2 Transmission?

Siwei Deng Wenzhou-Kean University: Kean University Wenzhou Campus Xuechen Tian Wenzhou-Kean University: Kean University Wenzhou Campus Robert Belshaw Wenzhou-Kean University: Kean University Wenzhou Campus Jinfeng Zhou China Biodiversity Conservation and Green Development Foundation (CBCGDF) Siyuan Zhang China Biodiversity Conservation and Green Development Foundation (CBCGDF) Yixin Yang Wenzhou-Kean University: Kean University Wenzhou Campus Chang Huang Wenzhou-Kean University: Kean University Wenzhou Campus Weikang Chen Wenzhou-Kean University: Kean University Wenzhou Campus Hailu Qiu Wenzhou-Kean University: Kean University Wenzhou Campus Siew Woh Choo (  cwoh@wku.edu.cn ) Wenzhou-Kean University: Kean University Wenzhou Campus


Introduction
Coronavirus disease 2019 , caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection, has become a global pandemic 1 . At the beginning of the pandemic in early 2020, it was proposed that the Malayan pangolin (Manis javanica) was an intermediate host of bat-to-human transmission of SARS-CoV-2 2,3 , but subsequent sequence analysis has shown this not to be the case 4 . Nevertheless, the pangolin coronavirus (pCoV) is closely related to SARS-CoV-2 and may provide insights into its pathology.
Clinical and histopathological studies of COVID-19 patients reported some dermatologic manifestations such as petechiae (a rash and haemorrhagic dot-like areas) [5][6][7][8] . The reason for the skin eruption in COVID-19 patients is unknown, but may be due to a SARS-CoV-2 interaction with skin [5][6][7][8] . Also, it has been suggested that angiotensin-converting enzyme 2 (ACE2), which is used by SARS-CoV-2 to enter the host 9 , can be highly expressed in keratinocytes 10,11 . These results suggested that percutaneous transmission might be an overlooked route of SARS-CoV-2 infection.
In this study, we investigated a possible pCoV skin infection in a Malayan pangolin named "Dahu" seized by customs in the Guangdong province (China) before dying 12 . We previously reported that Dahu's brain (cerebrum and cerebellum) and lungs were infected by pCoV 12 and here we con rm that Dahu's skin was also infected. With RNA-seq data, we found that the most up-regulated molecular pathway was the same SARS-CoV-2-associated pathway previously identi ed in the lungs of infected humans 13 . The main differences in pathway dysregulation could be attributed to the pangolin's unique immune characteristics, namely inactivation of the interferon epsilon (IFNE) and interferon-induced with helicase C domain 1 (IFIH1, also known as MDA5) 14,15 . We also found that many endogenous retroviral (ERV) genes were differentially expressed in Dahu skin, possibly because of the pCoV infection, suggesting that they might play a role in CoV infection. These ndings indicate that pCoV was replicating within the skin tissue of Dahu.

Overview
We present an overview of the Methods, with further methodological details provided in the Supplementary Information. Brie y, Dahu skin tissue was harvested by veterinary o cers at Guangzhou Leader Animal Hospital. RNA was extracted from a lesion in the skin of the neck using Qiagen RNeasy Mini Kit and sequenced using Illumina HiSeq with a 2 × 150bp paired-end strategy. We compared our infected skin with published data from the healthy skin of another pangolin (project accession PRJNA283328; run SRR3923846). The read counts, normalised expression values and differentially expressed genes (DEGs) were generated according to the Cu inks pipeline 16 . For comparative analysis, the list of DEGs was compared with the gene list generated from human COVID-19 patients from the study of Blanco-Melo et al. (project accession PRJNA615032) 17 . Functional enrichment analyses including gene set enrichment analysis (GSEA) and over-representation analysis (ORA) were performed using clusterPro ler v3.18.0 18 . ERV elements were identi ed using TBLASTN 19 and in-house scripts.
Phylogenetic trees were constructed using MEGA-X 20 . Sequences were aligned using the MUSCLE algorithm 21 . The maximum-likelihood trees were inferred using the Tamura-Nei DNA substitution model and nodal support was estimated using 1,000 replicates.

Results
The presence of pCoV subgenomic mRNA in the skin To test whether Dahu skin was infected by pCoV (in addition to its lungs and brains 12 ), we searched for viral sequences among the transcriptome data. A total of 193 reads mapped to the reference pCoV genome ( Figure S1). The distribution of these mapped reads was consistent with the corresponding locations of the subgenomic mRNAs 22 . We also observed an individual read that spanned precisely the splicing sites: this read was 150 bp in length and its 5' 47 bp mapped to the 5' region of the pCoV genome, while its 3' 103 bp mapped to the 3' region of the same genome. This read indicates that the CoV RNA has been processed in the host cell 23 . Consistently, we detected the N gene (which is used in diagnostic human SARS-CoV-2 testing) in all the tested Dahu organ using qRT-PCR 12 , including the skin.
To con rm the identity of the viral RNA in Dahu skin, we compared consensus sequences from our mapped reads with CoV from other species (Figure 1; Table S1). We observed that the pCoV genome and genes from Dahu skin were almost identical (sequence identity = 99.2-100%) to the counterparts from another pangolin, Guangdong pCoV isolate MP789 24 , con rming the presence of pCoV RNA in the Dahu skin.

Dahu skin transcriptional responses to pCoV
To explore the host skin transcriptional responses, we identi ed 3,201 DEGs (1,810 up-regulated and 1,391 down-regulated) in Dahu skin (Table S2). Consistent with the results observed in human COVID-19 patients 13 , we found terms associated with immunity processes were signi cantly enriched and upregulated, such as terms related to non-interferon cytokine production, T-cell, neutrophil, myeloid cell, and mast cell differentiation and activation ( Figure 2D). Apoptotic signalling pathways were negatively regulated, especially apoptotic processes of immune cells, such as myeloid cells and leukocytes. Also, rRNA metabolic process, RNA splicing, mRNA transcription, translational processes, and protein folding were also up-regulated. Speci cally, we found genes associated with viral transcription and viral gene expression were up-regulated. Other terms, including complements and cellular response to biotic stimulus were up-regulated, while cell cycle processes were down-regulated.
Strikingly, we found that the coronavirus disease-COVID-19 pathway was the most signi cant upregulated pathway ( Figure 2E; Table S3; Figure S2). Ribosomal proteins were up-regulated ( Figure S3A; Figure 2E; Table S3) as CoV needs ribosome frameshifting to translate and replicate 25,26 . It has been shown that the nonstructural protein 1 of the SARS-CoV-2 is a major virulence factor that interferes with host RNA translation by binding to 40S ribosomal subunit 27,28 . Li summarised the functional relationships between host ribosomal proteins and viral infection 29 , and suggested most interactions are bene cial for viral protein translation and replication. We found ribosomal proteins that are crucial for viral infection were up-regulated in Dahu skin ( Figure S3), such as RPL3 30 , RPL18 [31][32][33][34] , and RPL24 35 .
Moreover, we found RPL9 and RPL22 which help virus particle assembly and viral gene expression, were up-regulated in Dahu-skin, but not human patient lungs ( Figure S3B), probably a new strategy to promote its replication in pangolin 36,37 .
Other up-regulated pathways include pathways in cancer, TNF signalling pathway, complement and coagulation cascades, chemokine signalling pathway, immunity pathway, and viral interaction/infection related pathways ( Figure 2E; Table S3). Furthermore, the DEGs were down-regulated in pathways associated with linoleic acid and lipid metabolism ( Figure 2E; Table S3). Again, these host responses are consistent with the core pathways identi ed in the human SARS-CoV-2 infection 13 .

Comparative analysis of transcription in Dahu skin and human lungs
We further investigated the similarities and differences of host responses between Dahu skin and the human lungs by comparing DEG lists between our study and an external human lung study 17 . Our comparative analysis revealed 366 DEGs shared between species (i.e., 'common DEGs' below), 2,835 Dahu skin-speci c DEGs, and 1,527 human lung-speci c DEGs. As anticipated, the common DEGs were enriched in coronavirus diseases-COVID-19 pathway, followed by MAPK signalling pathway, apoptosis, Ctype lectin receptor signalling pathway and Kaposi sarcoma-associated herpesvirus infection. These ndings are consistent with the pCoV infection in Dahu skin ( Figure 3A).
Interestingly, interferon-speci c responses were exclusively up-regulated in human lungs compared to Dahu skin. We found interferon responses (e.g., interferon-gamma production) signi cantly enriched and up-regulated in human lung-speci c DEGs ( Figure 3C-D; Table S4), but not in Dahu skin-speci c DEGs (Table S5). We also downloaded genes related to interferon signalling pathways in Reactome pathway database 38,39 and found 131 genes were expressed in Dahu skin and/or human lungs ( Figure 3E; Table  S6). Among them, 102 genes were expressed in human lungs, while 88 genes were expressed in Dahu skin ( Figure 3E). We found 42 of them (37 up-regulated and 5 down-regulated) were differentially expressed in human lungs, but none of them was differentially expressed in Dahu skin. It has been reported that pangolins have unique immune characteristics 14,15 : both IFNE and IFIH1 (a cytoplasmic RNA sensor that helps initiates the innate immune response to viral infection) genes have been pseudogenised. IFNE plays an important antiviral role in epithelial cells [40][41][42][43][44][45][46] . Therefore, IFNE-de cient pangolin skin may be susceptible to CoV infection. Therefore, it indicates pangolin might adopt alternative strategies to ght against CoV infection over evolutionary time.
In our comparative analysis, we found three enriched pathways in Dahu skin-speci c DEGs ( Figure 3B), in which malaria and Staphylococcus aureus infection pathways were up-regulated in Dahu skin-speci c DEGs, while arachidonic acid (AA) metabolism pathways were down-regulated. These pathways may indicate unique pangolin host skin responses. Malaria pathway is commonly up-regulated after SARS-CoV-2 infection 13 , and anti-malarial drugs have shown effects on inhibiting SARS-CoV-2 replication 47 . In addition, we discovered that malaria pathway was also up-regulated in human lungs. As it is one of the three Dahu skin-speci c pathways, it indicates malaria pathway was further up-regulated in Dahu skin. Consistent with other up-regulated in ammatory responses, we observed Staphylococcus aureus infection pathway was speci cally up-regulated in Dahu skin.
In addition, the AA metabolism pathways were down-regulated in Dahu skin-speci c DEGs ( Figure 3B), indicating that these pathways were suppressed by pCoV infection. It is known that AA pathways have inhibitory effects on coronavirus replication, suggesting that lipid metabolism could be a druggable target of coronavirus-infected patients 48 . Therefore, pCoV might also suppress AA metabolism pathway in pangolin skin to bene ts its replication.

Responses of endogenous retrovirus (ERV) gene expression in the pCoV skin infection
Then, we investigated the relationships between host ERV gene expression and pCoV infection. We identi ed 6,076 ERV genes by screening 3,162 known viral proteins from Swiss-Prot across the pangolin genome (Table S7). We found 466 genes expressed in infected and/or non-infected pangolin skin, in which 348 of them (81 up-regulated versus 267 down-regulated) were differentially expressed ( Figure 4A; Figure S4; Table S8), suggesting that the exogenous pCoV might suppress the expression of ERV genes after infection to bene t its replication. We found most of the ERV DEGs were env (43%), pol (31%), and gag (16%) ( Figure 4B; Table S8), while the compositions in the genome were 28%, 38%, and 20%, respectively ( Figure S5A; Table S8). The most abundant group of ERV DEGs were most closely related to mouse intracisternal A-particle (IAP) viruses ( Figure 4D; Table S8). Mouse and Hamster IAP viruses make up 19% of the ERV DEGs, while only 11% of the IAP sequences are found in the pangolin genome ( Figure  S5C; Table S7), which may suggest the importance of IAPs in CoV infection and their possible interaction with the virus.

Discussion
In this study, pCoV RNA was found in the Dahu skin. This Cov was almost identical to pCoV MP789 from the Guangdong pangolin. Dahu and the Guangdong pangolin MP789 were seized by Guangdong antismuggling bureaus after being smuggled into China 12 and kept together at Guangdong wildlife rescue centre. Therefore, it is possible that both pCoVs originated from the same source 12 .
Our functional enrichment analyses are generally consistent with the results observed in human patients with SARS-CoV-2 infection 13 . After being infected, the immune responses were activated via immune cells (lymphocytes, leukocytes, myeloid cells, neutrophils, and mast cells), stimulating non-interferon cytokine production and preventing immune cell apoptosis 49,50 . It has been shown that the cytokine storm is also an early feature after SARS-CoV-2 infection 51 , which can be cross-validated by up-regulating rRNA metabolic process, RNA splicing, mRNA transcription, translational processes, and protein folding. Also, the virus exploits host cell machinery to proliferate in the host 52 , which up-regulates genes associated with viral transcription and viral gene expression. Complement activation and hypercoagulation are frequently accompanied by SARS-CoV-2 infection 53 . Here, we also observed the activation of complement pathways in Dahu skin. Furthermore, cell cycle processes were suppressed in pangolin skin. It has been shown that CoV can arrest cell cycle to boost viral replication e ciency 54 through mechanisms such as regulating through cyclin-CDK complex 55,56 , p53-dependent pathway 57,58 , N protein of coronavirus 59-61 , and directly interacting host cell cycle proteins 62,63 .
At the pathway level, consistent with SARS-CoV-2 infection, our analysis showed that COVID-19 pathway, immunity and in ammation (except for IFN-related pathways), cell proliferation, and coagulation pathways were the most signi cant up-regulated pathways in the Dahu skin 13 . In the CoV-infected pangolin skin, the interferon-speci c pathways were not enriched, and the expressions of many interferon pathway-related genes were undetectable or/and not signi cantly differentially expressed. In other placental mammals, IFNE provides the rst line of defence in the skin and mucosa-protected organs (e.g., lungs) and senses viral infection by triggering downstream genes involved in the interferon pathways [40][41][42][43][44][45][46] . Therefore, the IFNE-mediated pathways, including interferon-stimulated gene responses, are unlikely to be activated or up-regulated in the naturally IFNE-and IFIH1-de cient pangolins upon CoV infection 14 .
In healthy pangolin skin, many ERV genes were expressed, indicating their biological signi cance. We believe that these ERV genes are bene cial to the host, such as to boost host's immunity [64][65][66][67] , which is especially important in IFNE-de cient pangolins. Our data showed that most of the ERV genes were down-regulated after being infected, leading to our speculation that pCoV might directly or indirectly suppress the ERV genes to bene t its proliferation.
pCoV might infect the skin cells via the known receptor ACE2. However, the expression of ACE2 in Dahu skin was extremely low or undetectable. It is consistent with low interferon level because ACE2 is stimulated by interferon 68 . Also, low RNA level does not mean ACE2 is not present. We used bulk RNA-seq in this study, but ACE2 is only expressed in the skin keratinocytes, and might only be detectable using single-cell technology 10,11 . In addition, the ACE2 RNA expression is low even in the healthy human lungs, where the protein is thought to be the main receptor 69 . Another possibility is that pCoV might infect the skin cells through alternative receptors such as DPP4, which is expressed in infected Dahu skin. Studies have shown that DPP4 is a potential binding target of SARS-CoV-2, leading to suitable therapeutic possibilities 70,71 . Gianotti, Zerbi, and Dodiuk-Gad 5 have demonstrated clinical features of SARS-CoV-2 infected patient skin, and they hypothesised that the patient skin eruption might be due to the viral interaction with the skin. We cannot rule out the possibility that the skin of these patients was indeed infected by SARS-CoV-2 and novel mechanisms may exist to assist CoV in infecting the skin.
A limitation of this study is that our observations are only based on one CoV-infected pangolin. We are also unsure whether our observations are caused by the unique pseudogenised IFNE and IFIH1, although this might resemble the response of people who have weakened immunity. A possible cause of our observations is contamination by pCoV-infected blood. However, a recent study showed the enrichment of representative pathways such as porphyrin and chlorophyll metabolism, oxidative phosphorylation, and tryptophan metabolism in human SARS-CoV-2 infected blood transcriptome 72 , but we did not observe these pathways enriched in the Dahu skin transcriptome. We therefore suggest that the presence of pCoV and the systematic host's responses to CoV infection observed in Dahu skin is unlikely to be due to blood contamination.

Conclusion
We report the presence of pCoV subgenomic (spliced) mRNAs, which only occur in infected cells, and transcriptomic hallmarks of host response to CoV infection in Dahu skin. We suggest it is the strongest evidence so far for CoV infection in the skin of a placental mammal. Also, the pathway dysregulation was consistent with CoV infection of an organism lacking IFNE and IFIH1, and we observed the downregulation of ERV genes. This study emphasises the possibility of SARS-CoV-2 transmission via skin in humans.
"Animal ethics approval: This work was approved (reference number: GF(2019)BASE08) by the Biology and Science Ethics Committee of the China Biodiversity Conservation and Green Development Foundation (CBCGDF)." Figure 1 Phylogenetic trees of the coronavirus RNA and gene sequences from Dahu skin. The phylogenetic tree was generated using the maximum likelihood method. Bootstrap numbers were generated in 1,000

Figures
replicates. Nodes with bootstrap support values of 70 or greater are indicated.