(2023) 24:16
Ong et al. BMC Molecular and Cell Biology
https://doi.org/10.1186/s12860-023-00475-4

BMC Molecular andCell Biology

Open Access

RESEARCH

Using RNA‑seq to identify suitable
housekeeping genes for hypoxia studies
in human adipose‑derived stem cells
Huan Ting Ong1,2*, Cecilia M. Prêle1,2,3,4 and Rodney J. Dilley1,2,5Abstract
Background Hypoxic culture conditions have been used to study the impact of oxygen deprivation has on gene
expression in a number of disease models. However, hypoxia response elements present in the promoter regionsof some commonly used housekeeping genes, such as GAPDH and PGK1, can confound the relative gene expression analysis. Thus, there is ongoing debate as to which housekeeping gene is appropriate for studies investigating
hypoxia-induced cell responses. Specifically, there is still contradicting information for which housekeeping genes
are stable in hypoxia cultures of mesenchymal stem cells. In this study, candidate housekeeping genes curated fromthe literature were matched to RNAseq data of normoxic and hypoxic human adipose-derived stem cell cultures todetermine if gene expression was modulated by hypoxia or not. Expression levels of selected candidates were used to
calculate coefficient of variation. Then, accounting for the mean coefficient of variation, and normalised log twofold
change, genes were ranked and shortlisted, before validating with qRT-PCR. Housekeeping gene suitability were thendetermined using GeNorm, NormFinder, BestKeeper, comparativeCt , RefFinder, and the Livak method.
Results Gene expression levels of 78 candidate genes identified in the literature were analysed in the RNAseq dataset generated from hADSC cultured under Nx and Hx conditions. From the dataset, 15 candidates with coefficientof variation ≤ 0.15 were identified, where differential expression analysis results further shortlisted 8 genes with least
variation in expression levels. The top 4 housekeeping gene candidates, ALAS1, RRP1, GUSB, and POLR2B, were chosen
for qRT-PCR validation. Additionally, 18S, a ribosomal RNA commonly used as housekeeping gene but not detectedin the RNAseq method, was added to the list of housekeeping gene candidates to validate. From qRT-PCR results, 18S
and RRP1 were determined to be stably expressed in cells cultured under hypoxic conditions.
Conclusions We have demonstrated that 18S and RRP1 are suitable housekeeping genes for use in hypoxia studies
with human adipose-derived stem cell and should be used in combination. Additionally, these data shown that thecommonly used GAPDH and PGK1 are not suitable housekeeping genes for investigations into the effect of hypoxia inhuman adipose-derived stem cell.
Keywords Reference genes, qRT-PCR, RNAseq, Hypoxia, Human adipose-derived stem cells, ADSC

*Correspondence:
Huan Ting Ong
huanting.ong@uwa.edu.au
Full list of author information is available at the end of the article© The Author(s) 2023. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which
permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to theoriginal author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images orother third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line
to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory
regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of thislicence, visit http://​creat​iveco​mmons.​org/​licen​ses/​by/4.​0/. The Creative Commons Public Domain Dedication waiver (http://​creat​iveco​
mmons.​org/​publi​cdoma​in/​zero/1.​0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Ong et al. BMC Molecular and Cell Biology(2023) 24:16

Background
Housekeeping genes (HKG) are stable, constitutively
expressed reference genes which are used to normalise
the expression of genes-of-interest in a variety of molecular techniques. HKG need to be ubiquitously and stably
expressed in order to accurately quantify any treatmentrelated differential gene expression. For this purpose,
it is usual to select HKG involved in supporting cellularfunctions [1]. Typically, HKG are selected from related
literature, or most commonly determined from previous
work performed in the same laboratory, sometimes without sufficient regard to the stability of that HKG to thedifferent treatments being applied. The stability of HKG
expression should be examined for each experimental
design, and specific cell types [2], as it cannot be assumed
that a “one-size-fits-all” approach is suitable.
Hypoxic (Hx) culture conditions have been used
to model cellular physiological responses in various
ischemic disorders, such as ischemic brain injury [3, 4],tumour progression and hypoxic microenvironment [5,6], diabetes and obesity [7, 8], and toxin-induced reactive oxygen species in chronic kidney diseases [9]. The
molecular response to hypoxia is typically driven by
hypoxia response elements (HRE) present in the promoter region of genes. Functional HRE contain a target
cis-acting DNA sequence that binds hypoxia-inducible
factor -1 (HIF-1) protein complex which regulates gene
expression [10, 11]. As hypoxia is a fundamental change
to culture conditions, other genes regulated by hypoxia
can be involved in metabolic switching pathways suchas redox homeostasis and glycolysis, hence a large number of complex interactions make it difficult to select for
HKG based on the absence of HRE alone. A number ofcommonly used HKG contain HRE and are regulated byhypoxia, for example GAPDH [12–14] and PGK1 [10]. A
list of 114 genes which are induced by hypoxia is shown
in Additional file 1 [15, 16]).
While various HKG have been reported as suitable forhypoxia studies (summarised in Table 1), there is ongoing debate in the literature as to which HKG are stably
expressed in hypoxia [17–27], especially in mesenchymal stem cells [20, 28–30]. Thus, using a standard protocol for identifying multiple reference genes, we selected
RPL13A, TBP and GAPDH from the literature and performed a preliminary qRT-PCR to study gene regulation in Hx human adipose-derived stem cells (hADSC),
but found that these genes, unlike what they’ve beenreported, were not stably expressed (Additional file 2).
This represents a major setback which requires an alternative solution.
Alternative to testing potential HKG from the literature one at a time, approaches have been developed to
identify and quantify a panel of selected HKG candidates

Page 2 of 15

Table 1 Literature reported HKG suitable for hypoxia studies
Gene Name

Gene Symbol Reference

18S rRNA

18S

[17]

beta-actin

ACTB[18, 19]

ribosomal protein lateral stalk subunit
P0

RPLP0

[18]

beta-glucuronidase

GUSB[18, 20]transferrin receptor

TFRC

[18]

beta-2-microglobulin

B2M[21, 22]

ribosomal protein L13a

RPL13A[21, 23]

RNA polymerase II

RPII

[21]

28S rRNA

28S[17, 24]

hypoxanthine phosphoribosyltransferase 1

HPRT1[23, 25]

tyrosine 3/tryptophan 5-monooxygenase activation protein

YWHAZ

[20]

TATA-box binding protein

TBP[20, 22]glyceraldehyde 3-phosphate dehydrogenase

GAPDH[19, 26, 28–30]

peptidylprolyl isomerase A

PPIA

[27]

based on the gene expression values from qRT-PCR data,such as GeNorm [31] which uses a geometric averaging
method, NormFinder [32] a model-based variance estimation method, BestKeeper [33] a pair-wise comparison method, comparative Ct method [34] by relative
expression, and RefFinder [35] which uses a geometricallyaveraging weighted ranks for the listed methods above.
However, these approaches require numerous qRT-PCR
experiments with a set of at least 12 HKG candidates andinvolve a larger number of samples and reagents.
With the advancement of next generation sequencing technology, RNA sequencing (RNAseq) is an attractive tool that can be used to obtain a broad comparison
between transcriptome profiles [1]. Some recent studies
have used RNAseq datasets to select stable HKG [36–38]
customized for a unique experimental design. Lemmaet al. 2018 [37] selected HKG from RNAseq data using
normalised log2 fold change (L2FC), while Tang et al.
2017 [36] and Wang et al. 2019 [38] normalised for coefficient of variation of mRNA transcript expressions whenselecting for HKG. Due to the complex nature of gene
expression and regulation by hypoxia, RNAseq may be aneffective way to screen large transcriptome data for HKG
selection.
In this study, we aim to use RNAseq data from hADSC
to screen HKG from the literature and validate theirexpressions under hypoxic regulation, in order to identify
the most suitable HKG in Hx hADSC. RNAseq datasets
were generated from hADSC from 5 donors and cultured under Nx or Hx conditions in vitro. RNAseq andbioinformatic analysis was used to identify differences

Ong et al. BMC Molecular and Cell Biology(2023) 24:16

in the transcriptome between donors and culture conditions. Specifically, the stability of candidate HKG incells cultured under Nx and Hx conditions was determined by calculating the coefficient of variation (CV)of both normalised and transformed read counts, andL2FC obtained from differential gene expression. Potential HKG candidates were identified by low CV (≤ 0.15),
ranked and were validated using qRT-PCR where relative
expression of each selected HKG was calculated using
the 2−Ct Livak method [39], GeNorm method, NormFinder method, the BestKeeper method, Ct method
and RefFinder method.

Results
Differential gene expression analysis

To identify which genes were differentially expressed by
hADSC cultured under Hx conditions, the RNAseq transcriptional profile of five primary hADSC cell cultures(hADSC-A, -B, -C, -D, -E) grown under Hx conditions
were compared to the profile of hADSC cells cultured
under standard conditions (Nx). Analysis of the transcriptome revealed a total of 2800 differentially expressedtranscripts; Fig. 1 shows the top 50 transcripts of at least|2| L2FC. DESeq2 was also used to cluster each dataset
according to the similarities in transcript count distribution (Fig. 2). Principal Component Analysis (PCA, Fig. 3)
identified distinct clusters of differentially expressed
genes in hADSC cultured in Nx and Hx conditions, aswell as a donor-to-donor variability in transcriptome
profile between the hADSC lines.
Comparison of hypoxic effects on HKG identified
from the literature

Candidate HKG were identified from the literature(ACTB, ALAS1, B2M, GAPDH, GUSB, HPRT1, PGK1,
POLR2B, PPIA, RPL13A, RPLPO, RRP1, TBP, TFRC andYWHAZ) and matched for normalised gene expression
in the RNAseq dataset generated from hADSCs cultured
under Hx and Nx conditions (Fig. 4). Normalised counts
indicated a similar distribution between Hx and Nx datasets for 7 of the candidates (ACTB, ALAS1, B2M, GUSB,POLR2B, RRP1 and TBP). In contrast, another 7 genes(GAPDH, HPRT1, PGK1, PPIA, RPL13A, RPLP0 and
YWHAZ) increased expression in Hx samples comparedto Nx hADSC, and TFRC showed reduced expression in
Hx hADSC.
Selection of HKG from the RNAseq dataset

To compare the levels of gene expression, data was firstexpressed in transcripts per million (TPM) normalised
counts and DESeq2-normalised counts for the selected
HKG. Both TPM and DESeq2-normalised counts werethen used to calculate for CV. before the list of HKG were

Page 3 of 15

ranked smallest to largest for both CV and L2FC (Additional file 3). Then, the average CV was calculated from
HKG with CV ≤ 0.15, before normalising CV and L2FC
to the smallest corresponding values of the 78 candidates
and summed. HKG were ranked from smallest to thelargest value for the sum of these four methods. Table 2
shows the top 8 (least variation) HKG candidates calculated and the top 4 listed in Table 3 for validation with
qRT-PCR, with data analysis workflow summarised in
Additional file 4.
Selection for low CV in the Hx hADSC found 8 of
78 HKG candidates screened had a CV ≤ 15%. The
top 4 HKG were validated further by qRT-PCR alongwith PGK1 and 18S mRNA expression detected in fourhADSC cell lines cultured under Nx and Hx conditions.
Table 3 lists the raw Ct values for each HKG (ALAS1,
RRP1, 18S, PGK1 and GUSB) in each hADSC (POLR2B
was not detected in all hADSC but was expressed in acontrol cell line (keratinocytes) at raw Ct = 25.0 ± 0.92).
Of the expressed HKG, GUSB showed very high Ct values in hADSC-A, -B and -D, indicating very low expression of the gene (keratinocytes showed expression at raw
Ct = 25.5 ± 0.08). Hence, both POLR2B and GUSB were
unsuitable HKG for hADSC in hypoxia studies due to the
low to no expression levels on qRT-PCR.
Ct values were validated using GeNorm (geometricaverage), NormFinder (model-based variance estimation), BestKeeper (pair-wise correlations), comparativedelta-Ct, RefFinder (geometric average of weighted ranking) and Livak (2−Ct ) methods to quantify the stabilityof HKG expression in hADSC cultured under Nx and Hx
conditions. Figure 5a-e summarises the ranking of the 5HKG, where the closer each value is to 1, the more stably
expressed each HKG were. Figure 5f shows the relative
expression (2−Ct ) of each gene using one of the HKG
as a reference gene at each time. Based on the combined
analysis of all 6 methods (Fig. 5g), 18S and RRP1 wereidentified as most stable of the 5 HKG tested in our culture conditions.
Further validation of the suitability of 18S and RRP1 as
HKG (Fig. 6) was performed. The expression of a known
hypoxia-regulated gene – VEGFA, was normalised to
18S and RRP1 and the levels of expression in cells treatedwith Hx compared to those cultured under Nx conditions. The relative expression of VEGFA was stable across
2 biological replicates of hADSC (hADSC-C and -D) cultured in hypoxia or normoxia. Figure 6 shows calculated
relative expression (2−Ct ) of VEGFA was up-regulated
6.8-fold ± 1.4 when using 18S as HKG, compared to a
12-fold ± 0.6 upregulation when using RRP1 as HKG.
Used in combination, using the geometric mean of
the two selected HKG, a 9.5-fold ± 1.0 upregulation in
VEGFA gene expression was measured.

Ong et al. BMC Molecular and Cell Biology(2023) 24:16

Page 4 of 15

Fig. 1 Heatmap for top 50 transcripts differentially expressed in hADSC under Hx and Nx. Normalised expression of the top 50 transcript variantsfrom 5 human adipose-derived stem cell lines in hypoxic and normoxic conditions

Ong et al. BMC Molecular and Cell Biology(2023) 24:16

Page 5 of 15

Fig. 2 Sample-to-sample distances between hADSC lines. The sample-to-sample distances were plotted as a heatmap to show that greatest
clustering occurred within hypoxic (Hx, purple) and normoxic (Nx, orange) treatment groups, followed by individual hADSC lines and thenreplicatesFig. 3 Principal Component Analysis (PCA) of hADSC datasets. PCA plot of hADSC lines (A-E, colour-coded) in hypoxic (Hx, triangles) and normoxic(Nx, circles) treatment groups. The first principal component (PC1) represents the maximum variance direction of data distribution at 52%, and thesecond principal component (PC2) represents the second largest variation at 24%

Discussion
Since some typical HKG were found to be regulated in
hADSC by hypoxia, this study used RNAseq datasets toscreen the transcriptome for suitable HKG candidates tovalidate differential gene expression levels, regulated byhypoxia, by qRT-PCR. Further validation was performed
with qRT-PCR and analysed using 6 established methods to calculate the variation in relative gene expression.

Ong et al. BMC Molecular and Cell Biology(2023) 24:16

Page 6 of 15

Fig. 4 Distribution of counts for normalised expression of selected housekeeping genes. Normalised counts from the differential expression
analysis by RNAseq of 5 hADSC donor lines in either Hx or Nx conditions are plotted for 15 selected housekeeping gene candidates

18S and RRP1 were identified as stable HKG suitablefor use in hypoxia studies in hADSC in vitro. The suitability of HKG candidates was determined by having low
variation in the gene expression levels across test andcontrol samples. Furthermore, the raw Ct values of 18S
and RRP1 fell within the Ct range expected of high copy
number and commonly used HKG genes, which usually
range between 14 to 32 raw Ct [40] (Ct values are inverse
proportional to gene expression levels). In addition, theCt of these candidate HKG were no more than 10 raw Ct

values apart from the Ct of the genes-of-interest which
is important to reduce variation as a result of amplification errors. qRT-PCR validation identified the following
candidate HKGs suitable for hypoxia studies in hADSC,
these genes have been ranked from the most stable toleast stable: 1) 18S, 2) RRP1, 3) ALAS1, 4) PGK1 and 5)
GUSB.
Comparison between the RNAseq-screening approachand qRT-PCR (the gold standard), identified only RRP1
and PGK1 as suitable HKG for hypoxia studies. The result

Ong et al. BMC Molecular and Cell Biology(2023) 24:16

Page 7 of 15

Table 2 Housekeeping gene candidates – selected on smallest coefficient of variance in TPM and DESeq2-normalised counts
Gene

Norm

CV calculated
A

GUSB
RRP1
ALAS1
POLR2B
CHFR
GGA1
YWHAZ
NONO

B

L2FC
C

D

CV

Norm CV

Norm L2FC

E

TPM

0.19

0.37

0.11

0.04

0.27

DE

0.09

0.08

0.14

0.04

0.28

TPM

0.30

0.30

0.08

0.13

0.15

DE

0.08

0.08

0.09

0.10

0.16

TPM

0.15

0.35

0.18

0.03

0.14

DE

0.11

0.08

0.06

0.06

0.14

TPM

0.12

0.33

0.27

0.06

0.38

DE

0.13

0.05

0.12

0.05

0.38

TPM

0.24

0.15

0.23

0.13

0.12

DE

0.05

0.30

0.10

0.08

0.11

TPM

0.42

0.47

0.11

0.15

0.14

DE

0.20

0.16

0.13

0.15

0.14

TPM

0.06

0.05

0.07

0.22

0.08

DE

0.28

0.34

0.18

0.18

0.08

TPM

0.06

0.06

0.11

0.09

0.14

DE

0.23

0.36

0.19

0.13

0.15

Arbitrary
Ranking
Sum

Rank

1.1

1.0

a

2.1

1

0.02

1.2

2.2

3.4

2

0.05

0.01a

1.0a

6.4

7.4

3

-0.12

0.01

1.1

14.0

15.1

4

0.12

0.01

1.2

14.3

15.5

5

-0.15

0.02

1.8

17.9

19.7

6

0.39

0.01

1.1

47.2

48.2

7

0.46

0.02

1.2

55.2

56.4

8

a

0.01

0.01

-0.02

Coefficient of variance was calculated from TPM and DESeq2-normalised values with L2FC values extracted from DESeq2. The average CV was calculated from HKG
with CV ≤ 0.15 (bolded), before normalising CV and L2FC to the smallest corresponding values (a) of the 78 candidates and summed. HKG were ranked from smallest
to the largest value for the sum of these four methods

for GUSB, ALAS1 and POLR2B were not consistent
between these methods. 18S, an rRNA not identified in
the transcriptome data since the library preparation wasenriched for mRNA, was validated by qRT-PCR against
VEGFA, a known hypoxia-regulated gene, was one of themost stable genes tested in this study, evident also in aprevious study [17]. This represents a limitation in using
RNAseq to identify for HKG as this approach does notaccount for total RNA and may exclude potential nonmRNA HKG candidates.
One criterion for HKG candidates was to be constitutively expressed in normal cell functions [1]. 18S rRNA
is a small subunit rRNA, which together with the large
subunit rRNA forms ribosomes essential for mRNAtranslation. RRP1, encodes for the nuclear protein ribosomal RNA processing protein 1 homolog A, which is
recruited at late stages of nucleologenesis during mitosis[41]. From the validation experiments using VEGFA as aknown hypoxia-regulated gene, both 18S and RRP1 fulfilled the selection of HKG based on cellular functions,and stable replicable relative expressions with VEGFA in
2 biological hADSC replicates (hADSC-C and -D).
Both 18S and RRP1 are suitable HKG for hADSC
hypoxia studies either individually, or used in combination. 18S is expressed at a higher level (Ct ~ 20) and RRP1
has a lower expression (Ct ~ 30), which allows users theoption to select the HKG expressed at a similar level totheir gene-of-interest to not overestimate relative gene
expression analysis. Alternatively, according to the Minimum Information for Publication of Quantitative RealTime PCR Experiments (MIQE) guidelines [28, 42], itis recommended to use a combination of multiple HKG
using geometric mean.
Apart from stable HKG, the results also identified genesthat are not suitable as HKG for hADSC hypoxia studies.
ALAS1 is a mitochondrial enzyme where cellular respiration and cell response to hypoxia takes place [43]. It was
found to be significantly down-regulated ­(2−△△Ct ≤ -1.4fold change) when normalised to 18S and hence notsuitable HKG for hADSC hypoxia studies. PGK1, encoding for a glycolytic enzyme involved in glucose metabolism and a known regulator of hypoxia, was confirmed
as a Hx-regulated gene in hADSC in both RNAseq and
qRT-PCR results, acting as a positive control measure ingene regulation under hypoxia. GUSB, ranked less stable
than PGK1, encodes for beta-glucuronidase, which is an
enzyme located in the lysosomes involved in glycosaminoglycan degradation. GUSB was expressed at highCt values, indicating low gene expression, suggesting adown-regulation of GUSB and the production of betaglucuronidase as a response to hypoxia regulation andhence not an appropriate HKG for this study.
In addition to the genes validated in qRT-PCR, othercommonly used HKG were also shown to be significantlyOng et al. BMC Molecular and Cell Biology(2023) 24:16

Page 8 of 15

Table 3 qRT-PCR analysis for selected housekeeping gene candidates in hypoxic and normoxic hADSC cultures
hADSC

A

Repeats

1

2

B

1

2

C

1

2

D

1

ALAS1

RRP1

18S

PGK1

GUSB

Hx

Nx

Hx

Nx

Hx

Nx

Hx

Nx

Hx

Nx

27.2

25.0

29.1

29.2

21.9

20.3

21.5

25.4

27.8

24.9

29.7

29.6

21.7

20.1

21.5

25.4

28.0

25.1

29.9

29.6

21.5

20.2

21.7

25.6

29.6

26.3

29.8

30.5

29.5

26.1

> 35

27.5

29.6

30.4

> 35

29.9

27.3

26.2

29.5

30.5

> 35

27.3

24.9

23.9

25.6

25.7

18.7

20.7

20.1

22.0

25.0

23.9

25.6

25.7

18.6

20.6

20.1

22.0

25.0

24.0

25.8

26.1

18.8

20.5

20.1

21.9

29.2

28.4

30.8

29.8

> 35

> 35

28.9

28.1

30.5

29.7

> 35

> 35

29.1

28.0

30.8

29.7

> 35

> 35

26.5

25.9

31.3

31.5

19.5

19.4

22.4

25.7

27.1

26.3

26.3

25.8

31.1

31.2

19.7

19.3

22.4

25.7

27.1

26.3

26.4

25.7

31.1

31.1

19.6

19.3

22.4

25.9

27.2

26.2

29.7

29.1

29.9

29.5

20.9

20.9

22.2

25.5

31.8

33.7

29.5

29.2

29.7

29.3

20.2

20.6

22.4

25.7

31.2

34.2

29.6

29.0

30.0

29.2

20.3

20.4

22.4

25.7

32.3

> 35

31.6

29.0

33.6

32.5

26.4

25.2

25.5

29.7

> 35

30.7

31.8

28.9

33.8

32.6

26.3

25.1

25.5

29.8

> 35

31.0

31.8

29.0

33.5

32.3

26.5

25.3

25.5

29.4

> 35

30.8

Raw Ct values (n = 4 cell lines, technical replicates n = 3) of ALAS1, RRP1, 18S, PGK1 and GUSB in hADSC-A, -B, -C and -D under normoxic (AN1, AN2, BN1, BN2, CN1, CN2,
DN1) and hypoxic (AH1, AH2, BH1, BH2, CH1, CH2, DH1) conditions. In samples where Ct values were 35 or above, gene expression was determined as being below
the level of detection and these data were not included

regulated by hypoxia in hADSC (adj-P ≤ 0.05) transcriptome data. These included GAPDH ­(23.12 = 8.7fold), HPRT ­(20.53 = 1.4-fold), PGK1 ­(23.53 = 11.6-fold),
PPIA ­(20.78 = 1.7-fold), RPL13A ­(21.23 = 2.3-fold), RPLP0­(21.58 = 3.0-fold), YWHAZ ­(20.39 = 8.7-fold) and TFRC­(2–3.23 = -9.4-fold). The literature indicated that HRE wereidentified in the promoter regions of PGK1 [10], GAPDH[13] and TFRC [44], but not HPRT, PPIA, RPL13A, RPLP0
or YWHAZ, suggesting an indirect hypoxia response tothese genes. Hence, this study provides added information to exclude genes that were unsuitable as HKG for
Hx hADSC, in the absence of HRE in corresponding promoter regions.
In comparison to the only other study in the literaturecomparing HKG in hADSC cultured under hypoxic conditions [20], this study identified different findings for 18S
and GUSB being suitable HKG, but also similar findings
in identifying unsuitable HKG such as GAPDH. While
it is accepted that GAPDH is regulated by hypoxia andnot suitable to be used as a HKG, the difference in findings for 18S and GUSB might be contributed to hADSC
being cultured in slightly different growth conditions in
Fink et al., such as high versus low passages, and hypoxia

cultures over 1 and 2 weeks. Our culture conditions havebeen optimised for the experimental question, and hence
reinforcing the need to validate suitable HKG according
to experimental design.
Limitations of study

RNA sequencing was performed in two separate batches
and the differential expression analysis did not correct
for batch effects. However, the selection of HKG using
RNAseq data for CV calculation was independently
performed between the hADSC lines. As the second
sequencing experiment (for hADSC-A, -B, -C and -D)
was sequenced more in-depth (25 to 57 million reads)
than the first (for hADSC-E at 5.9 to 8.7 million reads),
more confidence was given if CV ≤ 0.15 was more consistent in hADSC-A, -B, -C and -D when selecting for
HKG after the calculations.
This study described a straightforward RNAseq-based
screening approach to identify potential HKG, with theassumption that normalised count reads represent gene
expressions. However, the TPM normalisation method
may underestimate genes that are expressed at very
low levels. For example, hADSC-C mean expression

Ong et al. BMC Molecular and Cell Biology(2023) 24:16

Page 9 of 15

Fig. 5 Comparative study of published methods for quantitating variance in normalised Ct values of housekeeping genes in hADSC under Nx andHx culture conditions. Raw Ct values of the 5 HKG in hADSC under Nx and Hx conditions were normalised and ranked using A) GeNorm (geometricmean), B) NormFinder (model-based variance estimation), C) BestKeeper (pair-wise correlations) D) comparative delta-Ct, E) RefFinder (geometricmean of each ranking), F) relative expression calculated with ­2−△△Ct (Livak and Schmittgen 2001) using all HKG candidates as reference genes in a
step-wise manner to compare between Hx and Nx, and G) The final ranking of all HKG candidates in each hADSC

Ong et al. BMC Molecular and Cell Biology(2023) 24:16

Page 10 of 15

Materials and methods
Cell culture maintenanceFig. 6 Validation of 18S and RRP1 housekeeping genes by qRT-PCR
analysis for hypoxic stimulation of VEGFA expression. Relative
expression of VEGFA in hypoxic versus normoxic treated hADSC(2−Ct ). Relative expression was calculated using the geometricmean of both 18S and RRP1 reference genes (N = 3 replicates,
experiment repeated twice. Data are represented as mean ± SD)

for CDKN1A was 163.5 TPM ± 19.1 TPM at CV 11.7%,
but for CHFR was 0.6 ± 0.1 TPM at CV = 23.1%. In this
scenario, CDKN1A would be defined as “stable expression” and included for subsequent analysis, but CHFR,
where there was only 0.1 TPM standard deviation ofcount reads, would be excluded. This may be resolved if
the selection for potential HKG is performed for genes
expressed at similar levels as specific gene-of-interests ofthe study.
Additionally, this study may be limited by specificity of probe/primer designs in qRT-PCR, as POLR2B
was expressed in both RNAseq normalised counts, butwhen used for qRT-PCR validation, were found to not be
expressed at detectable levels. Further evaluation of qRTPCR primers for POLR2B isoforms in hADSC may be
warranted in future studies.

Conclusions
In conclusion, this study found RNAseq data able toidentify candidate housekeeping genes for use in hypoxia
studies of human adipose-derived mesenchymal stem
cells. After validation by qRT-PCR, 18S rRNA and RRP1
were identified as suitable HKG candidates and are recommended as stable housekeeping genes for future
hypoxia-induced gene expression studies in hADSC. The
differential expression of several commonly used housekeeping genes, provided further evidence to eliminate
them as suitable for HKG in future hypoxia studies ofhuman adipose-derived stem cells.

hADSC from 5 donors (A, B, C, D and E) were purchasedfrom Lonza (Lonza, Walkersville, Maryland) or ZenBio(ZenBio, Research Triangle Park, North Carolina), where
cells were previously characterised (general information
on hADSC donors and cell characterisation are provided
in Additional file 5). hADSC were cultured in Dulbecco’s
Modified Eagle’s medium (DMEM) containing 1 g/L
D-glucose (Gibco, Carlsbad, California), 2 mM GlutaMAX™ (hADSC-E in 4 mM L-glutamine) and 110 mg/L
Sodium Pyruvate, supplemented with 10% foetal bovine
serum (FBS, Bovogen, Keilor East, Australia) and 1% Penicillin–Streptomycin (P/S; Gibco, Carlsbad, California).
Spent media were replaced every 2–3 days. Cells wereincubated in standard culture conditions of 37 °C, with
5% carbon dioxide (­CO2) and 21% atmospheric oxygen­(O2) in a humidified cell culture incubator unless otherwise specified, and passaged at approximately 80% confluence. hADSC were used between passage 5 to 8.
A human eardrum keratinocyte cell line characterised previously [45] was cultured in DMEM containing
4.5 g/L D-glucose (Gibco, Carlsbad, California), 4 mM
L-glutamine) and 110 mg/L Sodium Pyruvate, supplemented with 10% FBS and 1% P/S. Spent media were
replaced every 3–4 days. Cells were incubated in standard culture conditions of 37 °C, with 5% ­CO2 and 21%
atmospheric oxygen (­O2) in a humidified cell culture
incubator and passaged at approximately 90% confluence.
Keratinocytes of passage 36 were used.
Hypoxia conditioning

hADSC were cultured to 80–90% confluence in 75 c­ m2
tissue culture flasks, rinsed with phosphate-bufferedsaline (PBS), pH 7.2 (Gibco, Carlsbad, California), andcultured in serum-free DMEM for at least 8 h before
replacing with fresh serum-free medium for a further
48 h. hADSC at this step were either incubated in standard (Nx; 37 °C, 5% C
­ O2) or Hx conditions in the GENbox Jar (BioMérieux, Craponne, France) produced byan anaerobic atmosphere generator sachet (BioMérieux,
Craponne, France) at 37 °C. A hypoxia indicator strip(BioMérieux, Craponne, France) was included in thesealed chamber to confirm gas composition (< 0.1% ­O2
and 15% ­CO2). At the end of the 48-h incubation, total
RNA was extracted.
RNA extraction

Total RNA from hADSC-E was extracted using theFavorPrep™ Tissue total RNA Mini Kit (Favorgen Biotech Corp, Ping Tung, Taiwan) according to the manufacturer’s protocol, before eluted RNA was cleaned up

Ong et al. BMC Molecular and Cell Biology(2023) 24:16

of genomic DNA (gDNA) using the RNeasy Plus Microkit (Qiagen, Hilden, Germany). Total RNA from hADSCA, -B, -C, -D, and keratinocytes were extracted using the
RNeasy Plus Micro kit according to the manufacturer’s
instructions. Briefly, cells were lysed and gDNA wasremoved via the gDNA Eliminator spin column. Total
RNA was then bound, collected and washed on the RNeasy MinElute spin column before being eluted at approximately 14 μL per sample. All RNA samples were stored
at -80 °C until required for RNA sequencing and cDNA
synthesis.
RNA samples were submitted to the Australia Genome
Research Facility (AGRF, Melbourne, Australia) for quality control, library preparation and RNA sequencing.
Briefly, RNA samples were quantified using a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, California) where samples of at least an RNA Integrity Number(RIN) of ≥ 8.0 and an A260/280 ratio of 1.8 – 2.0 were
accepted. RNA quantification results are provided in
Additional file 6.
Library preparation and RNA sequencing

For all samples, library preparation was performed using
the Illumina TruSeq® Stranded mRNA kit according to
the manufacturer’s protocol (Illumina, San Diego, California). Briefly, mRNA was purified and fragmented for
double-stranded cDNA synthesized using specific methods to improve strand specificity. Adapters were ligatedto DNA fragments to prepare for hybridization onto
a flow cell during RNA sequencing, and enriched byamplification.
For hADSC-E, samples were processed earlier in a separate batch to the rest. RNA sequencing was performedon the Illumina HiSeq 2500 platform High Output Modeusing the HiSeq Control Software (HCS) v2.2.68 and RealTime Analysis (RTA) v1.18.66.3 running on the Illumina
sequencing computer for single-end reads at 50 bp. The
Illumina bcl2fastq 2.18.0.12 pipeline was used next togenerate the sequence data.
For hADSC-A, -B, -C and -D, RNA sequencing was
performed on the Illumina NovaSeq 2500 platform
High Output Mode using the NovaSeq Control Software(NCS) v 1.6.0 and Real Time Analysis (RTA) v 3.4.4 performing real-time based calling on the NovaSeq instrument computer for single-end reads at 100 bp. The
Illumina bcl2fastq v 2.20.0.422 pipeline was used next
to generate the sequence data. All RNAseq files were
uploaded onto ArrayExpress accession E-MTAB-9979.

Page 11 of 15

Bioinformatic analysis

RNAseq raw fastq files were checked using the FastQC (v0.11.9) tool [46] to evaluate the quality of the sequenced
samples. Reads were trimmed with Trimmomatic (v0.36) [47] with parameters recommended in the Trimmomatic tool manual for single-end reads. Trimmed
reads were then aligned to the GrCh38 human reference genome assembly obtained from Ensembl using the
HISAT2 (v 2.1.0) tool [48], with the default parameters inthe HISAT2 tool manual. Mapped reads were processedusing SamBamba (v 0.6.6) [49]. For hADSC-A, -B, -C
and -D, lane 1 files were merged with their corresponding repeated lane 2 files. Multimappers were filtered out
using a mapping quality threshold of ≥ 5. Sorted and filtered Bam files were used as input for the featureCounts
function in subread [50] to count each transcript to the
GrCh38.90 reference genome. Sequencing details andprocessing steps are attached in Additional file 7. Countfiles (number of reads per transcript counted in data files)
were parsed to include the corresponding Ensembl Gene
IDs, before analysing with DESeq2 (v 1.26.0) [51] via R
studio.
Differential expression analyses were performed using
DESeq2 to compare differential gene expressions fromhADSC cultured in Hx against Nx. Low count reads(< 10) were removed and DESeq2 was used with themulti-factor design for “conditions” (Hx vs Nx) and “celltypes” (hADSC-A, -B, -C, -D and -E). Count data weretransformed using regularised logarithmic transformation (rlog function) in DESeq2 to compare between datasets in the clustering analyses. Differential expression
analysis was performed on raw count data and results
were shrunk using the Apeglm function [52] for better
estimation of distribution and exported as L2FC. Transcripts were considered differentially expressed (DE) if
L2FC ≥ 2 at adjusted-P ≤ 0.05.
In addition to the L2FC data, raw count reads from all
5 hADSC were also normalised to transcripts per million (TPM) by normalising for both gene lengths andsequence depths.
All analysis performed in R for this study can be foundin Additional file 4.
Ranking method

From the literature, a list of 78 previously known HKG
were curated. L2FC and DESeq2-normalised values wereextracted from the differential expression analysis dataof the 5 hADSC, and TPM were calculated using the formula below, for the 78 HKG.

Transcripts per million (TPM) = Raw counts ÷Gene length÷ Sequencing depth
1000Ong et al. BMC Molecular and Cell Biology(2023) 24:16

Page 12 of 15

Table 4 Description of TaqMan Assay probes used in qRT-PCR validation of housekeeping genes
Gene

Label

Amplicon (bp)Assay Design

Catalogue No

Reference

ALAS1

FAM

77

Spanning Exon 9—10

Hs00167441_m1

NM_000688.5

GUSB

FAM

96

Spanning Exon 8—9

Hs00939627_m1

NM_000181.3

POLR2B

FAM

58

Spanning Exon 6—7

Hs00946293_m1

NM_000938.2

RRP1

FAM

75

Spanning Exon 4—5

Hs01002953_m1

NM_003683.5

VEGFA

FAM

59

Spanning Exon 3—4

Hs00900055_m1

NM_001025366.2

18S

VIC

187

Within single exon

Hs99999901_s1

X03205.1

PGK1

VIC

75

Spanning Exon 4—5

Hs99999906_m1

NM_000291.3

Information of TaqMan Assay probes used, and the NCBI Reference Sequences for the corresponding genes. FAM: 5(6)-carboxy-fluorescein; VIC: 2’-chloro-7’phenyl-1,4dichloro-6-carboxy-fluorescein

The CV of each HKG between Hx and Nx was calculated with TPM or DE data using the following formula
below, for all 5 hADSC (­ CVABCDE) as well as individually­(CVA, ­C VB, ­C VC, ­C VD, ­C VE):
Co efficient of Variance (CV ) =Standard DeviationHx and Nx(TPMorDE)
MeanHx and Nx (TPMorDE)

An average of CV ( CV ) was first calculated using onlyCV ≤ 0.15. For example, ALAS1 had 7 CV values ≤ 0.15,
the average CV was calculated by averaging the 7 CV values only. Next, the top 8 genes were selected to normalise for both CV and L2FC by dividing each value with thesmallest corresponding values for an arbitrary score. For
example, from the averaged CV, ALAS1 was the smallest
at 0.013 so average CV of all other genes were normalised
to 0.013, and of the L2FC, GUSB was the smallest (closest
to 0) at 0.01 so L2FC of all other genes were normalised
to 0.01. These two arbitrary scores were then summed
and ranked from the smallest to the largest where the top
4 HKG were selected for validation via qRT-PCR. Additional file 3 contains all 78 HKG data, calculations, andranking.
cDNA synthesis

Complementary DNA was synthesised using the ­RT2
Omniscript kit (Qiagen, Hilden, Germany) according
to manufacturer’s instructions. Briefly, 1 μg of template
RNA was reverse transcribed in accordance with the
manufacturer’s protocol for each sample, together withthe corresponding no reverse transcriptase controls(NRT) and no template controls (NTC). For each reaction, 1 μM oligo(dT)15 primers (Promega Corporation,
Fitchburg, Wisconsin) were used, where reactions wereincubated at 37 °C for 60 min.
qRT‑PCR

qRT-PCR was performed using the QuantiFast Probe
PCR kit (Qiagen, Hilden, Germany) on cDNA generated from hADSC-A, -B, -C and -D RNA samples. Gene
expressions were performed using FAM-labelled and

VIC-labelled TaqMan Assay probes (Thermo Fisher
Scientific, Waltham, Massachusetts) as listed below inTable 4. 18S was not included in our RNAseq since it’s anrRNA, VEGFA was to validate the known regulation by
hypoxia. Briefly, 2.5 μL of cDNA, 5 μL of 2 × PCR Mastermix, 0.5 μL of TaqMan Assay probes and RNase-free­H2O were combined to a total of 10 μL per reaction perwell in a white 96-well full-skirt plate, along with appropriate controls for NRT and NTC. qRT-PCR was performed on the CFX Real-Time PCR Detection System(Bio-Rad Laboratories, Hercules, California) under thefollowing conditions: 95 °C for 3 min, and 40 cycles ofdenaturation at 95 °C for 3 s and annealing at 60 °C for
30 s. Threshold cycles (Ct) were determined using the
exponential growth phase and the baseline signal fromfluorescence versus cycle number plots.
Validation of HKG stability

Using RefFinder [35], raw Ct values for each HKG of all 4
hADSC lines were entered according to the instructions
to calculate for geometric average (GeNorm), modelbased variance estimation (NormFinder), pair-wisecorrelations (BestKeeper), Ct method, as well as thegeometric mean of the weighted ranks from the 4 methods to be expressed as RefFinder rank. Then, relative gene
expression levels were determined using 2−Ct method[39] for each hADSC to compare Hx (test) against Nx(control). The overall rank was averaged for each HKG.
Relative gene expression levels (n = 3) using 2−Ct wasalso performed for VEGFA against 18S and RRP1 in 2
biological replicates of hADSC in Hx (test) against Nx(control).
Abbreviations
qRT-PCR	Real-time quantitative polymerase chain reaction
cDNA	Complementary DNA
POLR2B	RNA polymerase II subunit B
ALAS1	5’-Aminolevulinate synthase 1
GUSB	Glucuronidase beta
RRP1	Ribosomal RNA processing 1
18S	18S Ribosomal RNA

Ong et al. BMC Molecular and Cell Biology(2023) 24:16

PGK1	Phosphoglycerate Kinase 1
VEGFA	Vascular endothelial growth factor A
Ct	Threshold cycles
NRT	No reverse transcriptase control
NTC	No template control
HKG	Housekeeping gene
TPM	Transcript per million
Hx	Hypoxic
Nx	Normoxic
L2FC	Log 2 Fold Change
CV	Coefficient of Variance
hADSC	Human adipose-derived stem cells
ACTB	Beta-actin
RPLP0	Ribosomal protein lateral stalk subunit P0
TFRC	Transferrin receptor
B2M	Beta-2-microglobulin
RPL13A	Ribosomal protein L13a
RPII	RNA polymerase II
28S	28S rRNA
HPRT1	Hypoxanthine phosphoribosyltransferase 1
YWHAZ	Tyrosine 3/tryptophan 5-monooxygenase activation protein
TBP	TATA-box binding protein
GAPDH	Glyceraldehyde 3-phosphate dehydrogenase
PPIA	Peptidylprolyl isomerase A
mRNA	Messenger RNA
HRE	Hypoxia response elements
HIF-1	Hypoxia-inducible factor -1

Page 13 of 15

Availability of data and materials
The datasets generated during this study are available in the ArrayExpress
repository with the accession code E-MTAB-9979. For correspondence andmaterial requests, please address them to HTO.

Declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not Applicable.
Competing interests
The authors declare no competing interests.
Author details
1
Ear Science Institute Australia, Nedlands, Western Australia, Australia. 2 Ear Sciences Centre, The University of Western Australia, Nedlands, Western Australia,
Australia. 3 Institute for Respiratory Health, The University of Western Australia,
Nedlands, Western Australia, Australia. 4 Discipline of Medical, Molecularand Forensic Sciences, Murdoch University, Murdoch, Western Australia,
Australia. 5 Centre for Cell Therapy and Regenerative Medicine, The Universityof Western Australia, Nedlands, Western Australia, Australia.
Received: 30 November 2022 Accepted: 15 March 2023

Supplementary Information
The online version contains supplementary material available at https://​doi.​
org/​10.​1186/​s12860-​023-​00475-4.
Additional file 1. List of 114 hypoxia-inducible genes, gene symbol, andgene ID.
Additional file 2. Preliminary qRT-PCR data analysis of RPL13A, TBP,GAPDH, B2M, VEGFA, using RNA extracted of hADSC cultured in normal
and hypoxic conditions.
Additional file 3. Ranking of normalised counts of RNAseq data of
Nx-hADSC and Hx-hADSC curated for selected HKG candidates andcalculated for CV.
Additional file 4. R markdown script of differential expression analysis ofRNAseq data using DESeq2.
Additional file 5. Summary of hADSC donor information and cell characterisation used in this study.
Additional file 6. Quality control and qualification results of RNA
extracted from hADSC.
Additional file 7. RNAseq processing efficiencies represented as percentage counts of original sequenced reads
Acknowledgements
The authors would like to acknowledge Australian Genome Research Facility
for performing the sequencing.
Authors’ contributions
The project was conceptualised by HTO and RJD. HTO and RJD contributed tothe design of the methodology. CMP contributed to the experiment designfor PCR validation. HTO performed all experiments and data analyses, andwas a major contributor in writing the manuscript. CMP and RJD substantiallyrevised the manuscript. All authors contributed to the interpretation of data,read and approved the final manuscript.
Funding
This work was supported by the Ear Science Institute Australia and the Gift of
Hearing Appeal. HTO was supported by an Australian Government ResearchTraining Program (RTP) Scholarship, and the Sarich and Stokes Research
Scholarship.

References
1. Eisenberg E, Levanon EY. Human housekeeping genes, revisited. Trends
Genet. 2013;29(10):569–74. https://​doi.​org/​10.​1016/j.​tig.​2013.​05.​010.
2. Caradec J, Sirab N, Keumeugni C, Moutereau S, Chimingqi M, Matar C,
et al. ‘Desperate house genes’: the dramatic example of hypoxia. Br J
Cancer. 2010;102:1037–43. https://​doi.​org/​10.​1038/​sj.​bjc.​66055​73.
3. Marti HH, Bernaudin M, Petit E, Bauer C. Neuroprotection and angiogenesis: dual role of erythropoietin in brain ischemia. Physiology News Physiol
Sci. 2000;15:225–9. https://​doi.​org/​10.​1152/​physi​ology​online.​2000.​15.5.​
225.
4. Schmidt-Kastner R, Haberkamp M, Schmitz C, Hankeln T, Burmester T.
Neuroglobin mRNA expression after transient global brain ischemia andprolonged hypoxia in cell culture. Brain Res. 2006;1103(1):173–80. https://​doi.​org/​10.​1016/j.​brain​res.​2006.​05.​047.
5. Maxwell PH, Dachs GU, Gleadle JM, Nicholls LG, Harris AL, Stratford IJ,
et al. Hypoxia-inducible factor-1 modulates gene expression in solid
tumors and influences both angiogenesis and tumor growth. Proc Natl
Acad Sci U S A. 1997;94(15):8104–9. https://​doi.​org/​10.​1073/​pnas.​94.​15.​8104.
6. Freese KE, Kokai L, Edwards RP, Philips BJ, Sheikh MA, Kelley J, et al.
Adipose-derived stem cells and their role in human cancer development,
growth, progression, and metastasis: a systematic review. Cancer Res.
2015;75(7):1161–8. https://​doi.​org/​10.​1158/​0008-​5472.​can-​14-​2744.
7. Yin J, Gao Z, He Q, Zhou D, Guo Z, Ye J. Role of hypoxia in obesity-induced
disorders of glucose and lipid metabolism in adipose tissue. Am J Physiol
Endocrinol Metab. 2009;296(2):E333–42. https://​doi.​org/​10.​1152/​ajpen​do.​90760.​2008.
8. Hodson L, Humphreys SM, Karpe F, Frayn KN. Metabolic signatures ofhuman adipose tissue hypoxia in obesity. Diabetes. 2013;62:1417–25.
https://​doi.​org/​10.​2337/​db12-​1032.
9. Khanh VC, Ohneda K, Kato T, Yamashita T, Sato F, Tachi K, et al. Uremic
toxins affect the imbalance of redox state and overexpression of prolylhydroxylase 2 in human adipose tissue-derived mesenchymal stem cells
involved in wound healing. Stem Cells Dev. 2017;26(13):948–63. https://​doi.​org/​10.​1089/​scd.​2016.​0326.
10. Semenza GL. Hypoxia-inducible factors in physiology and medicine. Cell.
2012;148(3):399–408. https://​doi.​org/​10.​1016/j.​cell.​2012.​01.​021.
11. Semenza GL, Agani F, Booth G, Forsythe J, Iyer N, Jiang BH, et al. Structural and functional analysis of hypoxia-inducible factor 1. Kidney Intl.
1997;51(2):553–5. https://​doi.​org/​10.​1038/​ki.​1997.​77.

Ong et al. BMC Molecular and Cell Biology(2023) 24:16

12. Lu S, Gu X, Hoestje S, Epner DE. Identification of an additional hypoxia
responsive element in the glyceraldehyde-3-phosphate dehydrogenase
gene promoter. Biochimica et Biophysica Acta (BBA) Gene Struct Express.
2002;1574(2):152–6. https://​doi.​org/​10.​1016/​s0167-​4781(01)​00359-1.
13. Higashimura Y, Nakajima Y, Yamaji R, Harada N, Shibasaki F, Nakano Y, et al.
Up-regulation of glyceraldehyde-3-phosphate dehydrogenase gene
expression by HIF-1 activity depending on Sp1 in hypoxic breast cancer
cells. Arch Biochem Biophys. 2011;509(1):1–8. https://​doi.​org/​10.​1016/j.​
abb.​2011.​02.​011.
14. Graven KK, Troxler RF, Kornfeld H, Panchenko MV, Farber HW. Regulationof endothelial cell glyceraldehyde-3-phosphate dehydrogenase expression by hypoxia. J Biol Chem. 1994;269(39):24446–53. https://​doi.​org/​10.​1016/​S0021-​9258(19)​51104-8.
15. Wenger RH, Stiehl DP, Camenisch G. Integration of oxygen signaling at
the consensus HRE. Science STKE. 2005;2005(306):re12. https://​doi.​org/​10.​1126/​stke.​30620​05re12.
16. Dengler VL, Galbraith MD, Espinosa JM. Transcriptional regulation by
hypoxia inducible factors. Crit Rev Biochem Mol Biol. 2014;49(1):1–15.
https://​doi.​org/​10.​3109/​10409​238.​2013.​838205.
17. Yang Y, Fan W, Zhu L, Zhao T, Ma L, Wu Y, et al. Effects of hypoxia onmRNA expression of housekeeping genes in rat brain tissue and primarycultured neural cells. Front Med China. 2008;2:239–43. https://​doi.​org/​10.​
1007/​s11684-​008-​0045-7.
18. Bakhashab S, Lary S, Ahmed F, Schulten HJ, Bashir A, Ahmed FW, et al.
Reference genes for expression studies in hypoxia and hyperglycemia models in human umbilical vein endothelial cells. G3 (Bethesda).
2014;4(11):2159–65. https://​doi.​org/​10.​1534/​g3.​114.​013102.
19. Moein S, Javanmard S, Abedi M, Izadpanahi M, Gheisari Y. Identificationof appropriate housekeeping genes for gene expression analysis in longterm hypoxia-treated kidney cells. Adv Biomed Res. 2017;6:15. https://​doi.​
org/​10.​4103/​2277-​9175.​200790.
20. Fink T, Lund P, Pilgaard L, Rasmussen JG, Duroux M, Zachar V. Instabilityof standard PCR reference genes in adipose-derived stem cells during
propagation, differentiation and hypoxic exposure. BMC Mol Biol.
2008;9:98. https://​doi.​org/​10.​1186/​1471-​2199-9-​98.
21. Foldager C, Munir S, Ulrik-Vinther M, Søballe K, Bünger C, Lind M. Validation of suitable house keeping genes for hypoxia-cultured humanchondrocytes. BMC Mol Biol. 2009;10:94. https://​doi.​org/​10.​1186/​
1471-​2199-​10-​94.
22. Baddela VS, Baufeld A, Yenuganti VR, Vanselow J, Singh D. Suitable
housekeeping genes for normalization of transcript abundance analysis
by real-time RT-PCR in cultured bovine granulosa cells during hypoxia
and differential cell plating density. Reprod Biol Endocrinol. 2014;12:118.
https://​doi.​org/​10.​1186/​1477-​7827-​12-​118.
23. Yao L, Chen X, Tian Y, Lu H, Zhang P, Shi Q, et al. Selection of housekeeping genes for normalization of RT-PCR in hypoxic neural stem cells of
rat in vitro. Mol Biol Rep. 2012;39(1):569–76. https://​doi.​org/​10.​1007/​
s11033-​011-​0772-8.
24. Zhong H, Simons JW. Direct comparison of GAPDH, β-actin, cyclophilin,and 28S rRNA as internal standards for quantifying RNA levels underhypoxia. Biochem Biophys Res Commun. 1999;259(3):523–6. https://​doi.​
org/​10.​1006/​bbrc.​1999.​0815.
25. de Kok JB, Roelofs RW, Giesendorf BA, Pennings JL, Waas ET, Feuth T, et al.
Normalization of gene expression measurements in tumor tissues: comparison of 13 endogenous control genes. Lab Invest. 2005;85(1):154–9.
https://​doi.​org/​10.​1038/​labin​vest.​37002​08.
26. Said HM, Hagemann C, Stojic J, Schoemig B, Vince GH, Flentje M, et al.
GAPDH is not regulated in human glioblastoma under hypoxic conditions. BMC Mol Biol. 2007;8:55. https://​doi.​org/​10.​1186/​1471-​2199-8-​55.
27. Bielecka ZF, Malinowska A, Brodaczewska KK, Klemba A, Kieda C,
Krasowski P, et al. Hypoxic 3D in vitro culture models reveal distinct resistance processes to TKIs in renal cancer cells. Cell Biosci. 2017;7:71. https://​doi.​org/​10.​1186/​s13578-​017-​0197-8.
28. Grayson WL, Zhao F, Bunnell B, Ma T. Hypoxia enhances proliferation andtissue formation of human mesenchymal stem cells. Biochem Biophys Res
Commun. 2007;358(3):948–53. https://​doi.​org/​10.​1016/j.​bbrc.​2007.​05.​054.
29. Ohnishi S, Yasuda T, Kitamura S, Nagaya N. Effect of hypoxia on gene
expression of bone marrow-derived mesenchymal stem cells and mononuclear cells. Stem Cells. 2007;25:1166–77. https://​doi.​org/​10.​1634/​stemc​
ells.​2006-​0347.

Page 14 of 15

30. Isik B, Thaler R, Goksu BB, Conley SM, Al-Khafaji H, Mohan A, et al. Hypoxic
preconditioning induces epigenetic changes and modifies swine
mesenchymal stem cell angiogenesis and senescence in experimental
atherosclerotic renal artery stenosis. Stem Cell Res Therap. 2021;12:240.
https://​doi.​org/​10.​1186/​s13287-​021-​02310-z.
31. Vandesompele J, de Preter K, Pattyn F, Poppe B, van Roy N, de Paepe A,
et al. Accurate normalization of real-time quantitative RT-PCR data by
geometric averaging of multiple internal control genes. Genome Biol.
2002;3:research0034.1.
32. Andersen CL, Jensen JL, Ørntoft TF. Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation
approach to identify genes suited for normalization, applied to bladderand colon cancer data sets. Cancer Res. 2004;64(15):5245–50. https://​doi.​
org/​10.​1158/​0008-​5472.​can-​04-​0496.
33. Pfaffl MW, Tichopad A, Prgomet C, Neuvians TP. Determination of stable
housekeeping genes, differentially regulated target genes and sample
integrity: BestKeeper – Excel-based tool using pair-wise correlations.
Biotechnol Lett. 2004;26(6):509–15. https://​doi.​org/​10.​1023/B:​BILE.​00000​
19559.​84305.​47.
34. Silver N, Best S, Jiang J, Thein SL. Selection of housekeeping genes forgene expression studies in human reticulocytes using real-time PCR. BMC
Mol Biol. 2006;7:33. https://​doi.​org/​10.​1186/​1471-​2199-7-​33.
35. Xie F, Xiao P, Chen D, Xu L, Zhang B. miRDeepFinder: a miRNA analysis
tool for deep sequencing of plant small RNAs. Plant Mol Biol. 2012;80:75–
84. https://​doi.​org/​10.​1007/​s11103-​012-​9885-2.
36. Tang X, Zhang N, Si H, Calderón-Urrea A. Selection and validation ofreference genes for RT-qPCR analysis in potato under abiotic stress. Plant
Methods. 2017;13:85. https://​doi.​org/​10.​1186/​s13007-​017-​0238-7.
37. Lemma S, Avnet S, Meade M, Chano T, Baldini N. Validation of suitable
housekeeping genes for the normalization of mRNA expression for studying tumor acidosis. Int J Mol Sci. 2018;19(10):2930. https://​doi.​org/​10.​3390/​ijms1​91029​30.
38. Wang Z, Lyu Z, Pan L, Zeng G, Randhawa P. Defining housekeeping
genes suitable for RNA-seq analysis of the human allograft kidney biopsy
tissue. BMC Med Genomics. 2019;12(1):86. https://​doi.​org/​10.​1186/​
s12920-​019-​0538-z.
39. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using
real-time quantitative PCR and the 2(-delta delta C(T)) method. Methods.
2001;25(4):402–8. https://​doi.​org/​10.​1006/​meth.​2001.​1262.
40. Valente V, Teixeira SA, Neder L, Okamoto OK, Oba-Shinjo SM, Marie SK,
et al. Selection of suitable housekeeping genes for expression analysis
in glioblastoma using quantitative RT-PCR. BMC Mol Biol. 2009;10:17.
https://​doi.​org/​10.​1186/​1471-​2199-​10-​17.
41. Savino TM, Bastos R, Jansen E, Hernandez-Verdun D. The nucleolar
antigen Nop52, the human homologue of the yeast ribosomal RNA
processing RRP1, is recruited at late stages of nucleologenesis. J Cell Sci.
1999;112(Pt 12):1889–900.
42. Bustin SA, Benes V, Garson JA, Hellemans J, Huggett J, Kubista M, et al. The
MIQE guidelines: Minimum information for publication of quantitative
real-time PCR experiments. Clin Chem. 2009;55(4):611–22. https://​doi.​
org/​10.​1373/​clinc​hem.​2008.​112797.
43. Solaini G, Baracca A, Lenaz G, Sgarbi G. Hypoxia and mitochondrial oxidative metabolism. Biochim Biophys Acta. 2010;1797(6–7):1171–7. https://​doi.​org/​10.​1016/j.​bbabio.​2010.​02.​011.
44. Lok CN, Ponka P. Identification of a hypoxia response element in thetransferrin receptor gene. J Biol Chem. 1999;274(34):24147–52. https://​doi.​org/​10.​1074/​jbc.​274.​34.​24147.
45. Redmond SL, Levin B, Heel KA, Atlas MD, Marano RJ. Phenotypic andgenotypic profile of human tympanic membrane derived cultured cells. J
Mol Histol. 2011;42(1):15–25. https://​doi.​org/​10.​1007/​s10735-​010-​9303-5.
46. Andrews S, Krueger F, Segonds-Pichon A, Biggins L, Krueger C, Wingett S.
FastQC: A quality control tool for high throughput sequence data. Babraham Institute. 2012. http://​www.​bioin​forma​tics.​babra​ham.​ac.​uk/​proje​cts/​
fastqc. Accessed 11 July 2019.
47. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for
Illumina sequence data. Bioinformatics. 2014;30(15):2114–20. https://​doi.​
org/​10.​1093/​bioin​forma​tics/​btu170.
48. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low
memory requirements. Nat Methods. 2015;12(4):357–60. https://​doi.​org/​
10.​1038/​nmeth.​3317.

Ong et al. BMC Molecular and Cell Biology(2023) 24:16

Page 15 of 15

49. Tarasov A, Vilella AJ, Cuppen E, Nijman IJ, Prins P. Sambamba: fast processing of NGS alignment formats. Bioinformatics. 2015;31(12):2032–4.
https://​doi.​org/​10.​1093/​bioin​forma​tics/​btv098.
50. Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose
program for assigning sequence reads to genomic features. Bioinformatics. 2014;30(7):923–30. https://​doi.​org/​10.​1093/​bioin​forma​tics/​btt656.
51. Love MI, Huber W, Anders S. Moderated estimation of fold change anddispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.
https://​doi.​org/​10.​1186/​s13059-​014-​0550-8.
52. Zhu A, Ibrahim JG, Love MI. Heavy-tailed prior distributions for sequence
count data: removing the noise and preserving large differences. Bioinformatics. 2019;35(12):2084–92. https://​doi.​org/​10.​1093/​bioin​forma​tics/​
bty895.

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Ready to submit your research ? Choose BMC and benefit from:• fast, convenient online submission• thorough peer review by experienced researchers in your field• rapid publication on acceptance• support for research data, including large and complex data types• gold Open Access which fosters wider collaboration and increased citations• maximum visibility for your research: over 100M website views per year
At BMC, research is always in progress.
Learn more biomedcentral.com/submissions