Cotton (Gossypium spp.) is a leading economic crop that is grown in more than 50 countries. The cottonseeds, once regarded as the by-product of fibre production, contain a rich supply of unsaturated fatty acids, proteins and vitamins. To date, the annual production of cottonseeds has the potential to meet the protein requirements for 550 million people globally, which shows great potential as a food resource amidst a growing food shortage (Janga et al., 2019). However, the utilization of cottonseed for food purposes is limited owing to the presence of ‘pigment glands’, which contains gossypol and its derivatives that are toxic to humans (Gao et al., 2020).

To study how pigment gland cells differentiate and to reveal the gene regulatory network in gland morphogenesis, scRNA-Seq was performed using a pair of NILs (gland cotton ‘CCRI12’ and glandless cotton ‘CCRI12gl’). The 1-week-old cotyledons were enzymatically digested, and the purified protoplasts were labelled with a 10x genomics barcode for high-throughput sequencing (Figure 1a). A total of 9186 individual cells, including 4790 cells from ‘CCRI12’ and 4396 cells from ‘CCRI12gl’, were obtained after cell filtering process (Figure S1, Table S1) and were divided into 12 clusters based on highly variable genes (Figure 1b, Figure S2).

To verify and correct cell group classifications, the expression profile of reported marker genes in 12 cell clusters of cotton cotyledons was studied. The clusters 0, 1 and 4 were identified as spongy mesophyll cells (SMC) due to the enrichment of a photosynthesis-related gene LHCB, and clusters 3 and 6 were identified as palisade mesophyll cells (PMC) due to the high expression of RBCS. The dominantly expressed GSTF9 marked cluster 5 as epidermal cells (EPC), and GSTL3 identified cluster 10 as the primordial cells (PRC) that could differentiate. In addition, cell type that specifically expressed MYB44, PXY, LTP and CYP82A3 marked the clusters 2, 7, 8 and 11 as guard cells (GC), xylem cells (XC), parenchyma cells (PAC) and phloem cells (PHC), respectively (Figures 1c,d, Table S2).

No well-known marker gene for pigment gland cells has been reported to date. GoPGF is the key factor that controls the biogenesis of pigment glands (Ma et al., 2016). However, the expression of GoPGF in different cell types has not been studied. Cluster 9 was identified in the cotyledons of gland cotton ‘CCRI12’ but not glandless cotton ‘CCRI12gl’, and GoPGF was specifically detected in the cells of cluster 9. This led to the tentative annotation of cluster 9 of the cotyledon cells as pigment gland cells (PGC; Figure 1c).

A pseudotime analysis was performed to uncover the differentiation relationships of cotyledon cell types. The study of the individual cell distribution and trajectory revealed that the PRCs originated earlier than the PGCs, suggesting that the PGCs could have differentiated from the PRCs (Figure 1e, Figure S3). In addition, four representative genes were selected to show their distribution and expression levels in PRCs and PGCs (Figure 1f).

To explore the potential regulators of gland development of cotton, the highly expressed genes in each cell cluster were studied. A total of 9325 DEGs were obtained with 1430 DEGs preferentially expressed in PGC, while the other cell clusters contained a range from 572 to 1704 (Table S3). Other than GoPGF, the previously reported GhERF105 (Wu et al., 2021), which determines the biogenesis of pigment glands in cotton leaves, was also identified as PGC-specific gene in our scRNA-Seq data. These results suggested the reliability of scRNA-Seq analyses in pigment gland cells and confirmed the accuracy of our classification of cell types.

To date, the regulators involving in pigment gland biogenesis that have been identified are all TFs, including CGF1, CGF2, GoPGF/CGF3, GaGRAS/GoSPGF and GhERF105. Therefore, this study focused on the TFs that were preferentially expressed in PGCs (Figure 1g, Table S4). qPCR revealed that most of the identified TFs were highly expressed in gland cells, while they were expressed at very low levels in the mesophyll cells (Figure S4). In addition, five candidate genes, including GoPGF, were selected for RNA in situ hybridization. These results showed that these genes have strong hybridization signals in the glandular structure (Figure 1h).

A 1.5-kb promoter upstream of the GoPGF initiation codon was cloned to drive the expression of GUS in the cotton gland cultivar ‘Coker312’. The transgenic lines that expressed GUS were obtained and used for histochemical staining. As shown in Figure 1i, a strong and clear GUS staining was restricted to the pigment glands. To our knowledge, this study is the first to use GUS staining to demonstrate that the transcription of GoPGF is restricted to gland cells.

Virus-induced gene silencing was utilized to quickly screen the candidate genes that controlled the formation of pigment glands. The results suggested that knock down of some candidate genes, including ERF13 and MYB14, mildly reduced the gland density (Figure S5). Notably, the GH_A05G3906 could modulate the contents of gossypol without changing the number of pigment glands, which suggested a possible biosynthetic pathway of sesquiterpene metabolism that is independent of pigment gland biogenesis (Figure S5). Among all the candidates, JUNGBRUNNEN 1 (GhJUB1) is of particular interest. Knock down of GhJUB1 inhibits gland biogenesis and the accumulation of gossypol.

The GhJUB1-silenced plants (TRV:JUB1) exhibited dramatically reduced pigment glands in newly growing tissues (Figure 1j–l). In addition, GhJUB1-silenced cotton plants exhibited gossypol levels of 15% in the leaves and 18% in the stems compared with those of the control plants (Figure 1m). These results revealed that GhJUB1 regulates gland morphogenesis, which was similar to that of GoPGF. To study the relation between GhJUB1 and GoPGF, the expression of GhJUB1 was studied in the GoPGF-silenced cotton (TRV:PGF), and the results showed that the expression of GhJUB1 dramatically decreased to an undetectable level (Figure 1n), which suggests that GhJUB1 could be downstream of GoPGF to control the biogenesis of pigment glands.

The pigment gland of cotton is a highly distinctive structure, which provides an ideal system to study cell differentiation and organogenesis. Our study indicates that the initiation of cell differentiation of pigment glands is highly correlated with the specific expression of key genes. One of the major constraints in the study of glandular development of cotton is the lack of natural glandless mutants. The scRNA-Seq data that we provide is invaluable for producing novel glandless mutants, which will greatly accelerate the breeding of commercially desired cotton varieties with glandless seeds.

棉花 ( Gossypium spp.) 是一种主要的经济作物，在 50 多个国家都有种植。棉籽曾被视为纤维生产的副产品，含有丰富的不饱和脂肪酸、蛋白质和维生素。迄今为止，棉籽的年产量有可能满足全球 5.5 亿人的蛋白质需求，这在日益严重的粮食短缺中显示出作为食物资源的巨大潜力（Janga 等人，2019年 ）。然而，由于含有对人体有毒的棉酚及其衍生物的“色素腺体”，棉籽的食品用途受到限制（Gao 等人，2020年 ）。

为了研究色素腺细胞如何分化并揭示腺体形态发生中的基因调控网络，使用一对 NIL（腺棉“CCRI12”和无腺棉“CCRI12gl”）进行 scRNA-Seq。1 周龄的子叶被酶消化，纯化的原生质体用 10x 基因组学条形码标记用于高通量测序（图 1a）。经过细胞过滤过程（图 S1，表 S1）后，共获得 9186 个单个细胞，包括来自“CCRI12”的 4790 个细胞和来自“CCRI12gl”的 4396 个细胞，并根据高度可变基因分为 12 个簇（图 1b，图 S2)。

为了验证和纠正细胞群分类，研究了棉花子叶的 12 个细胞簇中报告的标记基因的表达谱。由于光合作用相关基因LHCB的富集，簇 0、1 和 4 被鉴定为海绵状叶肉细胞 (SMC)，由于RBCS 的高表达，簇 3 和 6 被鉴定为栅栏叶肉细胞 (PMC )。显性表达的GSTF9将第 5 簇标记为表皮细胞 (EPC)，而GSTL3将第 10 簇标记为可以分化的原始细胞 (PRC)。此外，特异性表达MYB44、PXY、LTP和CYP82A3的细胞类型将簇 2、7、8 和 11 分别标记为保卫细胞 (GC)、木质部细胞 (XC)、薄壁组织细胞 (PAC) 和韧皮部细胞 (PHC)（图 1c、d、表 S2）。

迄今为止，尚未报道色素腺细胞的众所周知的标记基因。GoPGF 是控制色素腺体生物发生的关键因素 (Ma et al .,  2016 )。然而， GoPGF在不同细胞类型中的表达尚未得到研究。第 9 簇在腺棉“CCRI12”而非无腺棉“CCRI12gl”的子叶中被识别，并且在第 9 簇的细胞中特异性检测到GoPGF 。这导致将子叶细胞的第 9 簇初步注释为色素腺细胞（PGC；图 1c）。

进行伪时间分析以揭示子叶细胞类型的分化关系。对单个细胞分布和轨迹的研究表明，PRC 的起源早于 PGC，这表明 PGC 可能已经从 PRC 中分化出来（图 1e，图 S3）。此外，还选择了四个代表性基因来显示它们在 PRC 和 PGC 中的分布和表达水平（图 1f）。

为了探索棉花腺体发育的潜在调节因子，研究了每个细胞簇中的高表达基因。总共获得了 9325 个 DEG，其中 1430 个 DEG 优先在 PGC 中表达，而其他细胞簇的范围为 572 到 1704（表 S3）。除了GoPGF之外，先前报道的GhERF105（Wu等人，  2021 年）决定了棉花叶片中色素腺体的生物发生，在我们的 scRNA-Seq 数据中也被确定为 PGC 特异性基因。这些结果表明色素腺细胞中 scRNA-Seq 分析的可靠性，并证实了我们对细胞类型分类的准确性。

迄今为止，已确定的参与色素腺生物发生的调节因子均为TF，包括CGF1、CGF2、GoPGF/CGF3、GaGRAS/GoSPGF和GhERF105。因此，本研究重点关注优先在 PGC 中表达的转录因子（图 1g，表 S4）。qPCR 显示大多数已鉴定的转录因子在腺体细胞中高表达，而它们在叶肉细胞中的表达水平非常低（图 S4）。此外，还选择了包括GoPGF在内的五个候选基因用于 RNA原位杂交。这些结果表明这些基因在腺体结构中具有强杂交信号（图 1h）。

克隆了GoPGF起始密码子上游的 1.5-kb 启动子，以驱动GUS在棉腺栽培品种“Coker312”中的表达。获得表达GUS的转基因品系并用于组织化学染色。如图 1i 所示，强烈而清晰的 GUS 染色仅限于色素腺体。据我们所知，这项研究首次使用 GUS 染色来证明GoPGF的转录仅限于腺体细胞。

利用病毒诱导的基因沉默来快速筛选控制色素腺体形成的候选基因。结果表明，敲低一些候选基因（包括ERF13和MYB14）会轻度降低腺体密度（图 S5）。值得注意的是，GH_A05G3906可以在不改变色素腺体数量的情况下调节棉酚的含量，这表明倍半萜代谢的生物合成途径可能与色素腺体生物合成无关（图 S5）。在所有候选者中，JUNGBRUNNEN 1 ( GhJUB1 ) 特别令人感兴趣。敲低GhJUB1会抑制腺体生物发生和棉酚的积累。

GhJUB1沉默的植物 (TRV:JUB1) 在新生长的组织中表现出显着减少的色素腺体 (图1j-l)。此外，与对照植物相比，GhJUB1沉默的棉花植物在叶子中表现出 15% 的棉酚水平，在茎中表现出 18% 的棉酚水平（图 1m）。这些结果表明GhJUB1调节腺体形态发生，这与 GoPGF 相似。为了研究GhJUB1与GoPGF之间的关系，研究了GhJUB1在GoPGF沉默棉花（TRV:PGF）中的表达，结果表明GhJUB1的表达显着降低到检测不到的水平（图 1n），这表明 GhJUB1 可能在 GoPGF 的下游以控制色素腺体的生物发生。

棉花的色素腺是一种极具特色的结构，为研究细胞分化和器官发生提供了理想的系统。我们的研究表明，色素腺细胞分化的启动与关键基因的特异性表达高度相关。棉花腺体发育研究的主要限制之一是缺乏天然的无腺体突变体。我们提供的 scRNA-Seq 数据对于产生新型无腺体突变体非常宝贵，这将大大加速具有商业需求的无腺体种子棉花品种的育种。