# **Impacts of Molecular Structure on Nucleic Acid-Protein Interactions**

Edited by

Vaclav Brazda and Richard Bowater

Printed Edition of the Special Issue Published in *International Journal of Molecular Sciences*

www.mdpi.com/journal/ijms

## **Impacts of Molecular Structure on Nucleic Acid-Protein Interactions**

## **Impacts of Molecular Structure on Nucleic Acid-Protein Interactions**

Editors **Vaclav Brazda Richard Bowater**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editors* Vaclav Brazda Biophysical Chemistry and Molecular Oncology Institute of Biophysics of the Czech Academy of Sciences Brno Czech Republic

Richard Bowater School of Biological Sciences University of East Anglia Norwich United Kingdom

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *International Journal of Molecular Sciences* (ISSN 1422-0067) (available at: www.mdpi.com/journal/ ijms/special issues/Nucleic Acid Protein).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-7445-5 (Hbk) ISBN 978-3-0365-7444-8 (PDF)**

© 2023 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

 He studied at the Faculty of Science of Masaryk University in Brno and completed a postdoctoral fellowship at the University Health Network in Toronto, Canada, in 2005-2006. He teaches at the Faculty of Chemistry of the Brno University of Technology. His research investigates nucleic acids, the relationship between the structure and function of DNA, and their interactions with proteins. He is a co-author of the broadly used bioinformatics server bioinformatics.ibp.cz.

#### **Richard Bowater**

Richard Bowater is Professor of Biochemistry and Molecular Biology Education in the School of Biological Sciences at the University of East Anglia where he teaches biochemistry and molecular biology to all levels of university students. Richard has authored many primary publications and reviews that focus on his research interests, and also delivered pedagogical commentaries about teaching biochemical concepts to diverse audiences. He became a Senior Fellow of the Higher Education Academy (SFHEA) in the UK in 2015, and is a member of the Biochemical Society, Microbiology Society, and the Royal Society of Biology.

## **Preface to "Impacts of Molecular Structure on Nucleic Acid-Protein Interactions"**

This reprint presents a collection of research findings published in the Special Issue of the *International Journal of Molecular Sciences,* titled "Impacts of Molecular Structure on Nucleic Acid-Protein Interactions". These types of interactions are indispensable for many basic biological processes and the research findings reported here have important implications across a range of human diseases. We thank all authors for their contributions to the Special Issue and hope that you will find inspiration in their research findings.

#### **Vaclav Brazda and Richard Bowater**

*Editors*

### *Editorial* **Impacts of Molecular Structure on Nucleic Acid–Protein Interactions**

**Richard P. Bowater 1,\* and Václav Brázda 2,\***

1 School of Biological Sciences, University of East Anglia, Norwich Research Park, Norwich NR4 7TJ, UK

2 Institute of Biophysics of the Czech Academy of Sciences, Královopolská 135, 612 00 Brno, Czech Republic

**\*** Correspondence<EMAIL_ADDRESS>(R.P.B<EMAIL_ADDRESS>(V.B.)

**Keywords:** cruciform; DNA base sequence; DNA structure; G-quadruplex; nucleic acid–protein interactions; Z-DNA

Interactions between nucleic acids and proteins are some of the most important interactions in biology because they are the cornerstones for fundamental biological processes, such as replication, transcription, and recombination. Nucleic acids can adopt a wide range of structural conformations and this structural flexibility plays critical roles in their interactions with proteins [1]. This Special Issue of the International Journal of Molecular Sciences reports on diverse representatives of such interactions (Figure 1) across a wide range of biological systems.

RNA molecules can adopt great structural diversity due to the range of intramolecular interactions that can be formed within the single-stranded molecules [2]. DNA molecules are typically presented as double-stranded, right-handed B-form helices as their canonical structures, and this maximizes the thermodynamic stability of the molecule [3]. However, a significant body of research emphasizes that alternative, non-canonical DNA structures can exist, including double-stranded, left-handed Z-DNA, but also multi-strand structures such as G-quadruplexes (G4s), intercalated-motifs (i-motifs), triplexes, and cruciform structures. These non-B-DNA structures are usually characterized by the occurrence of single-stranded regions (loops) and/or sites of disrupted base pair stacking (junctions between continuous B-form DNA and the alternative structure) [3].

Variations in the structures of nucleic acids offer different binding sites for proteins, including interactions that focus on a range of sequence- and structure-specific nucleic acid targets. The structures of nucleic acids influence different aspects of biological activity, including physiological and pathological functions [4,5], themes that are addressed in this Special Issue. The collection of articles involve biophysical, biochemical, molecular biological and bioinformatics approaches and cover different biological systems, but there are some common themes among them: several refer to computational biology or bioinformatics approaches, or highlight additional information in DNA sequences [6–9]; several articles refer to different types of non-B DNA structures [10–12], with specific interests in quadruplexes [13–16]; several articles address different approaches and outcomes from proteins binding to DNA structures [11–13,17].

The wealth of DNA sequence information provided by genome-sequencing projects has brought new insights into the primary sequences of genomes and also about possible sequence-dependent local secondary structures [3,18]. Advances provided by such genome sequences are exemplified by the Human Genome Project, with complete telomere-totelomere sequences being finalised in 2022 [19,20]. As highlighted, this Special Issue includes several articles that report on computational biology or bioinformatics studies of DNA sequences [6–9], identifying unexpected and additional information within them (Figure 1A). In a mini review, Bartas et al. summarize current knowledge about the amino acid composition of various nucleic-acid-binding proteins, highlighting differences across proteins that bind in a sequence-specific manner compared to those that recognize local

**Citation:** Bowater, R.P.; Brázda, V. Impacts of Molecular Structure on Nucleic Acid–Protein Interactions. *Int. J. Mol. Sci.* **2023**, *24*, 407. https:// doi.org/10.3390/ijms24010407

Received: 16 December 2022 Accepted: 21 December 2022 Published: 26 December 2022

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

non-B-DNA structures and those that recognize both types of properties of nucleic acids [6]. Bioinformatic studies of repetitive DNA sequences in *Drosophila melanogaster* polytene chromosomes show that chromatin structure plays a crucial role in the regulation of gene activity [7]. Recent advances reported by Choi et al. demonstrate that nucleic acids may provide useful tools for building complex logic circuits [8]. Finally, for this grouping of articles, Víglaský explores the organization of genetic information in nucleic acids using a novel orthogonal representation, which proves to be useful in predicting the likelihood of particular regions of nucleic acids to form non-canonical motifs [9]. proteins that bind in a sequence-specific manner compared to those that recognize local non-B-DNA structures and those that recognize both types of properties of nucleic acids [6]. Bioinformatic studies of repetitive DNA sequences in *Drosophila melanogaster* polytene chromosomes show that chromatin structure plays a crucial role in the regulation of gene activity [7]. Recent advances reported by Choi et al. demonstrate that nucleic acids may provide useful tools for building complex logic circuits [8]. Finally, for this grouping of articles, Víglaský explores the organization of genetic information in nucleic acids using a novel orthogonal representation, which proves to be useful in predicting the likelihood of particular regions of nucleic acids to form non-canonical motifs [9].

acid composition of various nucleic-acid-binding proteins, highlighting differences across

*Int. J. Mol. Sci.* **2023**, *24*, x FOR PEER REVIEW 2 of 4

**Figure 1.** Nucleic acids provide a wide array of sequences and structures that are useful for biological processes. (**A**) Genome sequences store large amounts of information that can be accessed for biological processes and structure prediction by various computational algorithms or for building complex logic circuits. (**B**) Nucleic acids can adopt a range of structures, including those indicated. From left to right: (i) double-stranded, right-handed B-DNA; (ii) double-stranded, left-handed Z-DNA; (iii) two intramolecular hairpins can come together to form a cruciform; (iv) G-quadruplex, formed from four strands that can be parts of one molecule (as shown) or from different molecules; (v) a triplex can be formed when three strands come together, which can be parts of one molecule (as shown) or from different molecules. (**C**) The variety of structures of nucleic acids offer opportunities to be recognised by other molecules, such as proteins. The range of structures shown here may be recognised by different proteins, as indicated by the different colours. From the earliest days of genome sequence analysis, it was recognized that natural **Figure 1.** Nucleic acids provide a wide array of sequences and structures that are useful for biological processes. (**A**) Genome sequences store large amounts of information that can be accessed for biological processes and structure prediction by various computational algorithms or for building complex logic circuits. (**B**) Nucleic acids can adopt a range of structures, including those indicated. From left to right: (i) double-stranded, right-handed B-DNA; (ii) double-stranded, left-handed Z-DNA; (iii) two intramolecular hairpins can come together to form a cruciform; (iv) G-quadruplex, formed from four strands that can be parts of one molecule (as shown) or from different molecules; (v) a triplex can be formed when three strands come together, which can be parts of one molecule (as shown) or from different molecules. (**C**) The variety of structures of nucleic acids offer opportunities to be recognised by other molecules, such as proteins. The range of structures shown here may be recognised by different proteins, as indicated by the different colours.

DNA molecules contain a wide array of repeating sequences [3]. These types of sequences are particularly prone to adoption of non-canonical DNA structures, such as G4s, triplexes, and cruciforms (Figure 1B), which are all explored in this Special Issue [10–16]. Zhao and Usdin review the range of structures that can form in specific trinucleotide From the earliest days of genome sequence analysis, it was recognized that natural DNA molecules contain a wide array of repeating sequences [3]. These types of sequences are particularly prone to adoption of non-canonical DNA structures, such as G4s, triplexes, and cruciforms (Figure 1B), which are all explored in this Special Issue [10–16]. Zhao and Usdin review the range of structures that can form in specific trinucleotide repeats, highlighting how their expansion in length is important in the pathology of fragile X-related disorders in humans [10]. Left-handed Z-helices can form in both DNAs and RNAs with appropriate sequences, and searching databases containing protein structures identified

novel proteins predicted to bind them [11]. A different type of repetitive DNA sequence, inverted repeats, can adopt cruciform structures and many proteins have now been validated to bind to them [12]. A series of articles provide insights about G4s. Bezzi et al. suggest that putative G4s found in the SARS-CoV-2 RNA genome and the cellular proteins likely to interact with them may constitute interesting targets for antiviral drugs [13]. Putative G4s in viruses are explored further in a study that reveals a positive correlation between their frequencies in double-stranded DNA viruses and their hosts from archaea, bacteria, and eukaryotes, indicating that their close coevolution leads to reciprocal mimicking of genome organization [14]. The potential of compounds to target G4s was explored for Rhodamine 6G, which was shown to have high selectivity for G4s with parallel topology [15]. A bioinformatic study combined with circular dichroism measurements identified a stable G4 that is evolutionarily conserved amongst plants sensu lato (in Archaeplastida), and this may form an additional layer of regulatory networks [16].

The wide array of structures that can be adopted by nucleic acids offer different opportunities for proteins (and other molecules) to bind to, leading to different types of outcomes [4,5]. Some proteins recognise sequence-specific targets, but an increasing number are being shown to interact with non-canonical structural aspects of nucleic acids (Figure 1C). We have already referred to some articles in this Special Issue that describe such interactions [11–13]. Another study took advantage of available datasets and discovered new correlations between specific amino acid deviations in p53 proteins, showing a direct association between specific amino acid residues in the protein and changes in p53 functionality, and further highlighting the importance of p53 protein in processes that influence lifespan and aging [17].

To summarize, this Special Issue of the International Journal of Molecular Sciences reports on representatives of interactions between nucleic acids and proteins, with an emphasis on understanding how the structure of the nucleic acid influences such interactions. It is important to characterize these molecular complexes because many are essential requirements for the viability of cellular life due to their involvement in fundamental aspects of nucleic acid metabolism. It is now clear that the structural flexibility of the nucleic acids plays critical roles in their interactions with proteins, with important implications across a range of human diseases, including cancer and some infectious diseases. A deeper understanding of these molecular interactions will require the use of complementary methods and techniques [1]. As is described in this Special Issue, biophysical, biochemical, molecular biological and bioinformatics approaches will deliver useful advances across a wide range of biological systems.

**Author Contributions:** Writing—original draft preparation, review and editing, R.P.B. and V.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** Václav Brázda is supported by the Czech Science Foundation, grant number 22-21903S.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** No new data were created or analyzed in this study. Data sharing is not applicable to this article.

**Conflicts of Interest:** Both authors declare no conflict of interest.

#### **References**

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

## *Review* **Amino Acid Composition in Various Types of Nucleic Acid-Binding Proteins**

**Martin Bartas , Jiˇrí Cerve ˇn, Simona Guziurov ˇ á, Kristyna Slychko and Petr Peˇcinka \***

Department of Biology and Ecology, Institute of Environmental Technologies, Faculty of Science, University of Ostrava, 710 00 Ostrava, Czech Republic<EMAIL_ADDRESS>(M.B<EMAIL_ADDRESS>(J.C.); ˇ<EMAIL_ADDRESS>(S.G<EMAIL_ADDRESS>(K.S.)

**\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** Nucleic acid-binding proteins are traditionally divided into two categories: With the ability to bind DNA or RNA. In the light of new knowledge, such categorizing should be overcome because a large proportion of proteins can bind both DNA and RNA. Another even more important features of nucleic acid-binding proteins are so-called sequence or structure specificities. Proteins able to bind nucleic acids in a sequence-specific manner usually contain one or more of the well-defined structural motifs (zinc-fingers, leucine zipper, helix-turn-helix, or helix-loop-helix). In contrast, many proteins do not recognize nucleic acid sequence but rather local DNA or RNA structures (G-quadruplexes, i-motifs, triplexes, cruciforms, left-handed DNA/RNA form, and others). Finally, there are also proteins recognizing both sequence and local structural properties of nucleic acids (e.g., famous tumor suppressor p. 53). In this mini-review, we aim to summarize current knowledge about the amino acid composition of various types of nucleic acid-binding proteins with a special focus on significant enrichment and/or depletion in each category.

**Keywords:** DNA; RNA; protein binding; G-quadruplex; triplex; i-motif; Z-DNA; Z-RNA; cruciform; amino acid composition

#### **1. Introduction**

Interactions between proteins and nucleic acids (DNA and RNA) are central to all aspects of maintaining and accessing genetic information. Nucleic acid-binding proteins are mostly composed of at least one DNA or RNA-binding domain where the interfacing with amino acids takes place in a specific or nonspecific manner [1]. Identification of nucleic acidbinding proteins is one of the most important tasks in molecular biology. Currently, nucleic acid-binding proteins can be identified and further characterized by several experimental techniques, including pull-down assays [2,3], yeast one-hybrid system [4,5], electrophoretic mobility shift assays [6,7], chromatin immunoprecipitation [8,9], and by other specialized techniques [10,11]. However, it is time-consuming and expensive to identify nucleic acid-binding proteins by experimental approaches [12]. With the easy availability of a large amount of protein sequence data, there is a rapid development of computational approaches and prediction tools that can rapidly and reliably identify nucleic acid-binding proteins [13,14]. Several such tools model nucleic acid-binding abilities based on protein amino acid composition [15,16]. There is a growing interest in so-called noncanonical nucleic acid structures and proteins that preferentially bind them [17–24]. Noncanonical nucleic acid structures are DNA and RNA structures different from their basic form, i.e., double-stranded right-handed DNA or single-stranded RNA, and are often formed by simple nucleotide repeats [25–28]. Physiologically, they are represented mainly by Gquadruplexes [29], i-motifs [30], triplexes [31], R-loops [32], slipped hairpins [33], DNA cruciforms [34], RNA hairpins [35], and Z-DNA [36]. These DNA/RNA structures have important biological functions [37–42] and contribute to many human diseases [43–46]. It became more and more evident, that proteins preferentially interacting with these

**Citation:** Bartas, M.; Cerve ˇn, J.; ˇ Guziurová, S.; Slychko, K.; Peˇcinka, P. Amino Acid Composition in Various Types of Nucleic Acid-Binding Proteins. *Int. J. Mol. Sci.* **2021**, *22*, 922. https://doi.org/10.3390/ ijms22020922

Received: 29 December 2020 Accepted: 16 January 2021 Published: 18 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

structures share distinct amino acid features/fingerprints [47,48]. This mini-review aims to focus on the amino acid composition of various types of DNA and RNA-binding proteins and to compare the amino acid composition of proteins that prefer binding to different noncanonical forms of nucleic acids.

#### **2. Amino Acid Composition of Nucleic Acid-Binding Proteins**

According to the Gene Ontology (GO) knowledgebase, there are 5037 nucleic acidbinding proteins (filtering GO:0003676 term by "protein") with experimental evidence in *Homo sapiens* [49–51]. Of this number, 2572 are annotated as RNA-binding and 2439 as DNAbinding proteins (some proteins have both functions). 1768 human proteins are known to bind DNA in a sequence-specific manner. It would be interesting to quantify the overall amount of proteins binding nucleic acids in a structure-specific manner. Unfortunately, there is no such category yet. We strongly suggest revisions in this manner. Inspiration can be found in the following review papers/databases focused on specific properties of proteins binding to G-quadruplexes [19,52–54], cruciforms [55], and Z-DNA/Z-RNA [56].

#### *2.1. History*

Amino acid composition of some nucleic acid-binding proteins was intensively studied at the beginning of the 70s, when Koichi Iwai et al. determined that "calf-thymus histones comprise five main types which differ in amino acid composition and electrophoretic mobility: A glycine-rich, arginine-rich histone (also known as f2al or IV); a glutamic-acid-rich, arginine-rich histone (fe or III); a leucine-rich, intermediate type histone (f2a2 or IIb1); a serine-rich, slightly lysine-rich histone (f2b or IIb2); and an alanine-rich, very lysine-rich histone (f1 or I)" [57], by using specialized chromatographic technique followed by polyacrylamide gel electrophoresis. In 1975, from the comparison of 68 representative proteins and frequencies of 61 codons of the genetic code, it was found that the average amounts of lysine, aspartic acid, glutamic acid, and alanine are above the levels anticipated from the genetic code, and arginine, serine, leucine, cysteine, proline, and histidine are below such levels [58]. There are a couple of examples from the 90s and 2000s when amino acid substitution in nucleic acid-binding protein abolished its function, e.g., an arginine to lysine substitution in the bZIP (Basic Leucine Zipper) domain of an opaque-2 mutant in maize abolished specific DNA-binding [59], missense mutations (Met175Arg and Ser191Asn) abolishing DNA-binding of the osteoblast-specific transcription factor OSF2/CBFA1 in human patients with cleidocranial dysplasia [60], or impaired RNA-binding of fragile X mental retardation protein upon missense mutation IIe-304→Asn in one of its KH domain [61]. Recent advantages in sequencing and bioinformatic methods allow us to directly compare the amino acid composition of thousands of (not only) human nucleic acid-binding proteins [62,63]. One of the most popular programs for this purpose is, e.g., composition profiler [64], which is a web-based tool for semi-automatic discovery of enrichment or depletion of amino acids, either individually or grouped by their physicochemical or structural properties [64]. Scientists often find themselves in the situation when they only have a sequence of new "hypothetical" protein, derived mainly from transcriptome sequencing, and want to deduce its function [65]. In case that no meaningful alignment to protein with known function is available, there is still a way to get some useful information using only primary amino acid sequence and its composition. In 2003, Cai and Lin used a protein's amino acid composition and support vector machine (SVM) prediction to decide if protein belongs to one of three classes—rRNA-, RNA-, or DNA-binding [66]. Currently, there are also user-friendly web-based prediction tools called DNAbinder and PseDNA-Pro, which can predict if the submitted protein sequence has DNA-binding ability [12,67].

#### *2.2. Methods to Inspect the Amino Acid Composition of Proteins*

Several approaches are used to inspect the amino acid composition of nucleic acidbinding proteins. Basically, we can divide the methods into in vitro and in silico. In vitro approaches are necessary to obtain a sequence of the protein of interest. Although the development of large-scale genomic sequencing has greatly simplified the procedure of determining the primary structures of proteins, the genomic sequences of many organisms are still unknown, and also modifications such as post-translational events (citrullination, deamidation, polyglutamylation, . . . ) may prevent proper determination of the protein sequence [68]. Then, the complete characterization of the primary protein structure often requires a mass spectrometry method with minimal assistance from genomic data, i.e., de novo protein sequencing [68,69]. In silico approaches are based mostly on previous knowledge about primary protein sequence. There is currently a plentitude of bioinformatics tools designed for that purpose, see, e.g., [64,70–73].

#### *2.3. Amino Acid Composition of Nucleic Acid-Binding Proteins*

Nucleic acid-binding proteins are traditionally divided into two categories. The first category comprises proteins with the ability to bind DNA, and the second category comprises proteins that bind to RNA. This division is quite outdated, mainly because, from the historical perspective, proteins that bind RNA were typically considered as functionally distinct from proteins that bind DNA and studied independently. Interestingly, current gene ontology analyses reveal that DNA-binding is potentially a major function of the mRNA-binding proteins [74]. Nonetheless, several studies inspecting amino acid composition of DNA and/or RNA-binding proteins were published [75,76] and find that particular amino acid residues are generally enriched or depleted within these protein categories (see Table 1).

Another, even more important division of nucleic acid-binding proteins is based on a so-called sequence or structure-specific type of binding. Proteins able to bind nucleic acids in a sequence-specific manner usually contain one or more of the well-defined structural motifs. One of such motifs, zinc-finger, binds DNA (or RNA) through specific interaction with nucleotides and sugar-phosphate backbone. Tandem repeating of slightly different zinc-finger motifs in protein then allows to recognizing its consensus nucleic acid-binding sequence specifically. Cysteine and histidine amino acid residues are crucially important to coordinate Zn2+ binding in the largest and best-characterized subgroup of zinc-finger binding proteins named the Cys2His2 fold subgroup [77,78]. Other well-defined sequencespecific motifs—leucine zipper, helix-turn-helix, or helix-loop-helix—are listed in Table 1, together with their common signatures of amino acid residues.

**Table 1.** Types of nucleic acid-binding proteins. This table summarizes the main categories of nucleic acid-binding proteins. There are two points of view. At first, we can simply divide these proteins into DNA and RNA-binding ones (and a relatively small category of proteins that are able to bind both DNA and RNA). Secondly (and more importantly), we can distinguish proteins that specifically bind known sequence motifs (sequence-specific DNA/RNA-binding) and proteins, which specifically bind local DNA/RNA structures. Besides, keep in mind that this table is very simplified, and categories are divided to be reader-friendly. In fact, many of the DNA/RNA-binding proteins combine sequence and structure-specific binding mechanisms.

**Table 1.** *Cont.*

In contrast, many proteins do not recognize nucleic acid sequence but rather local DNA or RNA structures (G-quadruplexes, i-motifs, triplexes, cruciforms, left-handed DNA/RNA form, and others) [19,30,41,55,90]. Finally, there are also proteins recognizing both sequence and local structural properties of nucleic acids (e.g., famous tumor suppressor p53 [91], Myc-associated zinc finger protein (MAZ) [92,93], and many RNA-binding proteins [94])—these proteins usually contain sequence-specific binding domain(s) together with domain(s)/region(s) with preference to noncanonical nucleic acid structures [95–97]. In 2016, Wang et al. analyzed the abundance of intrinsic disorder in the DNA- and RNAbinding proteins in over 1000 species from Eukaryota, Bacteria, and Archaea domains of life [98]. They have revealed a very interesting phenomenon that DNA-binding proteins had significantly increased disorder content and were significantly enriched in disordered domains in Eukaryotes but not in Archaea and Bacteria. The RNA-binding proteins were significantly enriched in the disordered domains in Bacteria, Archaea, and Eukaryota, while the overall abundance of disorder in these proteins was significantly increased in Bacteria, Archaea, animals, and fungi [98]. Disordered domains or regions are also extensively present in chromatin-binding proteins [99,100]. Interestingly, some disordered proteins or regions show very high structural specificity to the different types of noncanonical nucleic acids. For instance, human protein SRSF1 (Serine/arginine-rich splicing factor 1) contains several intrinsically disordered regions [101], which are compositionally enriched in glycine (14.11% of overall amino acid residues) and arginine (17.39% of overall amino acid residues) content. It was previously shown that SRSF1 has a high affinity to RNA G-quadruplex structure [102]. Subsequent analyses have shown that the dataset of 77 G-quadruplex binding proteins is significantly globally enriched in arginine, glycine, aspartic acid, asparagine, and valine, and depleted in cysteine and other amino acid residues [47] (Figure 1). Finally, the common amino acid motif in the form of RGRGRGRGGGSGGSGGRGRG was derived, and most of the currently known G-quadruplex binding proteins contain at least some modification of it [47]. Using this motif, a new dataset of G-quadruplex binding proteins was predicted from the set of all human DNA/RNA-binding proteins [47], and some of them were independently experimentally validated (e.g., CIRBP, which is a cold-inducible RNA-binding protein in the study by Huang and colleagues [103]). A similar study focused on an amino acid composition of cruciform binding proteins was also published, and the significant enrichment for lysine and serine amino acid residues has been revealed [48]

(Figure 1). Unpublished results also indicate distinct amino acid profiles in Z-DNA/RNA and triplex binding proteins, both significantly enriched in aspartic acid and isoleucine and depleted in cysteine residues [89] (Figure 1). In future studies, it would be interesting to specifically analyze local amino acid composition (only in the nucleic acid interaction sites and their close neighborhood) in these proteins. Unfortunately for the vast majority of them, the knowledge about exact DNA/RNA binding site(s) is still missing.

**Figure 1.** Significantly enriched and depleted amino acid residues (one letter aa code) in the dataset of G-quadruplex binding proteins (top), cruciform binding proteins, triplex binding proteins, and Z-DNA-binding proteins. Using Bonferroni correction, only values lower than 0.0025 were taken as significant (*p* <0.0025; *p* <0.0010; *p* <0.0001). The size of arrows indicates the significance of enrichment/depletion on scale (highest, moderate, lowest). Figure compiled using data from [47,48,89]. Created with BioRender.com.

As was shown above, proteins that preferentially recognize noncanonical nucleic acid structures often have a distinct amino acid composition with particular significant global enrichment and/or depletion of different amino acid residues. Noncanonical structures and proteins preferentially binding them often play a critical role in physiological molecular processes [32,104,105], but also in the progression of human diseases, such as various cancer types and neurodegenerative diseases, reviewed in [55,106,107]. Knowledge about the amino acid composition of various proteins binding noncanonical nucleic acids can be utilized as an additional clue/fingerprint in discovering novel noncanonical nucleic acid-binding protein candidates and therapeutically utilized [108–110].

The scheme below depicts sequence and structure-specific nucleic acid-binding phenomena in a nutshell (Figure 2).

Almost every year, multiple novel noncanonical nucleic acid-binding proteins are identified. This year was, for instance, found that Guanine Nucleotide-Binding Protein-Like 1 (GNL1) binds RNA G-quadruplex structures in genes associated with Parkinson's disease [111], or that Small Nuclear Ribonucleoprotein Polypeptide A (SNRPA) directly binds to the BAG-1 mRNA through the G-quadruplex which can modulate BAG-1 expression level [112] (anti-apoptotic BAG-1 protein is known to be overexpressed in colorectal cancers [113]). Prediction of proteins that preferentially bind noncanonical DNA/RNA

structures, therefore, should be a logical first step towards rapid identification of novel therapeutic targets for future treatment of severe human diseases.

**Figure 2.** Types of nucleic acid-binding. The nucleic acid-binding mechanism can be basically divided into two main categories—sequence and structure-specific binding. (Left) Sequence-specific binding proteins recognize the variety of known DNA/RNA sequences via specific interaction with well-characterized protein motifs (zinc-fingers, helix-loop-helix, leucine zipper, helix-turn-helix, etc.). (Right) Structure-specific binding proteins recognize specific local structure(s) of nucleic acids, e.g., G-quadruplexes, i-motifs, cruciforms, triplexes, Z-DNA, and many others. In fact, it is a very common phenomenon that protein with the sequence-specific binding also prefers local DNA/RNA structure in its binding site or within the near neighborhood (e.g., p53), which is indicated by vertical black dashed line and arrows. Created with BioRender.com.

#### **3. Closing Remarks**

The global or local amino acid composition of nucleic acid-binding proteins is often overlooked and an unjustly underestimated parameter. Mainly statistically significant enrichment or depletion of particular amino acid residues may serve as a promising tool to predict novel proteins with a similar function, as it was confirmed e.g., for G-quadruplex binding proteins.

**Author Contributions:** Conceptualization, M.B. and P.P.; methodology, M.B.; resources, M.B., S.G., and K.S.; writing—original draft preparation, M.B., S.G, K.S., and J.C.; writing—review and editing, ˇ J.C and P.P.; visualization, M.B.; supervision, P.P. All authors have read and agreed to the published ˇ version of the manuscript.

**Funding:** This research was funded by University of Ostrava, SGS01/PˇrF/2020, and by National Agency for Agricultural Research (NAZV) of Czech Republic grant no. QK1810391 "Utilization of genomic and transcriptomic approaches to create genetic resources and breeding materials of poppy with specific traits".

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** No new data were created or analyzed in this study. Data sharing is not applicable to this article.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Abbreviations**

#### **References**

## *Article* **Chromatin Structure and "DNA Sequence View": The Role of Satellite DNA in Ectopic Pairing of the** *Drosophila* **X Polytene Chromosome**

**Aleksandr V. Zhuravlev 1,\* , Gennadii A. Zakharov 1,2, Ekaterina V. Anufrieva 3 , Anna V. Medvedeva 1 , Ekaterina A. Nikitina 1,3 and Elena V. Savvateeva-Popova 1**

**Abstract:** Chromatin 3D structure plays a crucial role in regulation of gene activity. Previous studies have envisioned spatial contact formations between chromatin domains with different epigenetic properties, protein compositions and transcription activity. This leaves specific DNA sequences that affect chromosome interactions. The *Drosophila melanogaster* polytene chromosomes are involved in non-allelic ectopic pairing. The mutant strain *agnts3*, a *Drosophila* model for Williams–Beuren syndrome, has an increased frequency of ectopic contacts (FEC) compared to the wild-type strain *Canton-S* (*CS*). Ectopic pairing can be mediated by some specific DNA sequences. In this study, using our Homology Segment Analysis software, we estimated the correlation between FEC and frequency of short matching DNA fragments (FMF) for all sections of the X chromosome of *Drosophila CS* and *agnts3* strains. With fragment lengths of 50 nucleotides (nt), *CS* showed a specific FEC–FMF correlation for 20% of the sections involved in ectopic contacts. The correlation was unspecific in *agnts3*, which may indicate the alternative epigenetic mechanisms affecting FEC in the mutant strain. Most of the fragments that specifically contributed to FMF were related to 1.688 or 372-bp middle repeats. Thus, middle repetitive DNA may serve as an organizer of ectopic pairing.

**Keywords:** *Drosophila*; polytene chromosomes; *Canton-S*; *agnostic*; ectopic pairing; 1.688 repeats; 372-bp repeats

#### **1. Introduction**

Spatial organization of the cell nucleus is an important factor defining the regulation of gene activity, as well as the processes of DNA replication, recombination and reparation. During interphase, chromosomes occupy separate territories in the nucleus, being radially arranged: gene-rich chromosome territories are localized toward the interior, while genepoor territories are close to periphery [1,2]. In human cells, regions of increased gene expression (ridges) are clustered in spatially distinct gene-enriched domains characterized by irregular forms and low chromatin condensation. These domains are predominantly located toward the nuclear interior. On the contrary, antiridges are relatively gene poor, condensed, transcriptionally inactive and localize closer to the cell envelope. Mechanisms behind the formation and maintenance of such 3D structures are still unclear, possibly due to the interaction between some unknown DNA sequences with nuclear matrix, specific proteins and/or non-coding RNAs [3]. Studying such mechanisms is necessary for understanding the process of gene regulation at the system level. As gene juxtaposition in nuclei facilitates specific chromosomal translocations, 3D chromatin structures can also predict the genetic rearrangements leading to carcinogenesis [4].

**Citation:** Zhuravlev, A.V.; Zakharov, G.A.; Anufrieva, E.V.; Medvedeva, A.V.; Nikitina, E.A.; Savvateeva-Popova, E.V. Chromatin Structure and "DNA Sequence View": The Role of Satellite DNA in Ectopic Pairing of the *Drosophila* X Polytene Chromosome. *Int. J. Mol. Sci.* **2021**, *22*, 8713. https://doi.org/10.3390/ ijms22168713

Academic Editors: Richard Bowater and Vaclav Brazda

Received: 30 July 2021 Accepted: 11 August 2021 Published: 13 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Chromosome territories are not rigidly fixed spatial units and show a significant percentage of intermingling, mostly at their borders, that is influenced by gene transcription [5]. Some genes are able to change their location in nuclei, being brought together with the help of actin and myosin motor proteins [6]. Transcription mainly occurs within the nuclear areas enriched in RNA polymerase II (RNAPII), known as transcription factories [7]. The expression level for given genes depends on their proximity to such a factory. The constitutively active genes may nucleate the factory, whereas the others relocalize to it upon their induction [4,7].

In addition to diploid cells, some organisms, such as *Diptera* species, also have polytene cells where chromatids are not segregated after multiple duplications. The giant polytene chromosomes of *Drosophila melanogaster* 3rd instar larvae are characterized by specific banding patterns, which can be revealed by electron and light microscopy. The densely packed thick "black" bands with high DNA content are transcriptionally repressed chromatin areas with low gene density. The largest among them are the intercalary heterochromatin (IH) bands, being late-replicating, under-replicated genomic areas prone to chromosomal breaks, constrictions and non-allelic ectopic contacts formation. *Drosophila* polytene chromosomes harbor about 250 IH sites. The "grey" bands are partially decondensed and more transcriptionally active compared to the IH. Interbands are the most active and the least condensed genomic areas, characterized by the "open" chromatin structure. Bands are united into cytological sections, such as 1A, 1B, 1C . . . and up to 20F for the X chromosome [8–10]. Each type of band contains specific proteins and is enriched for specific genetic elements. Interbands mostly contain the promoters of the constantly active housekeeping genes, being associated with open chromatin proteins such as CHRIZ/Chromator. The "grey" bands contain multiple active genes but lack CHRIZ, being enriched with RNApol II. The IH bands are composed of tissue-specific genes, being associated with SUUR, D1, lamin B and histone H1 proteins [11]. Notably, such structures are not unique for polytene nuclei, as the chromatin folding and protein composition are conserved in different fly tissues, being closely related to the morphology of the polytene chromosomes. At the same time, the ability to form distant contacts in polytene chromosomes is restricted by their lack of flexibility [12]. Replication timing is similar between polytene and diploid *Drosophila* cells. The late replicating black bands in various tissues correspond to silent chromatin types, with borders enriched for SUUR, lamin and H3K27me3 [12]. This makes polytene nuclei a convenient model to study 3D nuclear organization.

As previously shown by Horchstrasser et al. [13], polytene chromosomes occupy specific spatial domains similar to chromosome territories in diploid nuclei. Chromosomes extend across the nucleus in a Rabl orientation (i.e., their centromeres group near one pole of the nucleus and telomeres near the opposite pole). Chromosomes are coiled in a right-handed fashion, with their 2L and 2R arms being mostly next to each other, as well as the 3L and 3R arms. The loci enriched in IH and ectopic contacts are oriented toward the envelope. However, there are no stable intrachromosomal interactions beyond the distance of two cytological divisions. Thus, chromosomal configuration varies significantly, even though the contacts between the distant loci through the ectopic fibers were not addressed in that study.

Chromosome conformation capture technologies make it possible to investigate nuclear 3D organization with a resolution of the order of one to tens kb. The chromosome conformation capture (3C) method is used to estimate the average frequency of a contact between the two known chromosomal loci in a cell population. Chromosome conformation capture-on-chip (4C) technology allows spatial contacts of a selected genomic site to be assayed with all unknown distant sites. In the chromosome conformation capture carbon copy (5C) method, a massive analysis of contacts between specific loci across the entire genome is performed. The Hi-C method reveals spatial contacts at the level of the whole genome [14]. The complex net of such physical contacts uncovers the chromosome topology in detail, assaying interactions between genes and their regulatory elements.

Using Hi-C technology, the spatial structure of the *Drosophila* genome was studied for both diploid and polytene nuclei with a resolution of 15 kb [15]. Polytene bands nearly correspond to topologically associated domains (TADs) with a mean size of 195 kb. They are conserved for polytene and diploid nuclei. TADs are persistent throughout fly development, and are formed by the axial condensation of the chromatin fiber. The putative role of a TAD is DNA compaction rather than regulation of gene activity. Stable interactions between the different TADs were not observed, with is consistent with the variable longrange chromosomal conformation in the nucleus [13]. In the other study on *Drosophila* embryonic nuclei, the long-range interactions between Polycomb-repressed domains were found. The hierarchically organized domains were associated with active and repressive epigenetic modifications of chromatin [16]. Ectopic pairing was not considered in Hi-C studies, as its frequency was low and such long-range contacts were beyond the limit of the resolution. Thus, light and electron microscopy remain the most appropriate methods by which to study the ectopic pairing.

Ectopic contact is morphologically observed as an intimate association of IH bands or as an unstructured fiber connecting two IH bands. As the paired chromosome sections seem to be covalently linked due to the chromatids recombination/reparation, ectopic pairing is not disrupted upon squashing with acetic fixation [12,17]. The method of squashed preparations was used to estimate frequency of ectopic contacts (FEC) in several *Drosophila* strains. FEC can be expressed as the total number of ectopic contacts between a given section pair.

*agnts3* is a *D. melanogaster* mutant with a dysfunction of LIM kinase 1 (LIMK1), the main regulator of actin polymerization in nervous cells. This mutant strain shows multiple cognitive impairments, being the model object for Williams–Beuren syndrome [18]. a/trich *agnostic* locus (X:11AB) is predisposed to mutations. Its length varies, probably due to spontaneous unequal recombination [19]. Impairment of LIMK1 and actin dynamics affects the spatial organization of chromatin [20]. For *agnts3*, as well as the wild-type strains *Canton-S (CS)*, *Berlin* and *Oregon-R*, multiple polymorphisms were found in the *limk1* gene and flanking sequences. In *agnts3*, there is also mobile S-element insertion downstream *limk1* [21,22]. The *agnts3* FE*C*s are significantly higher compared to the wildtype strains [22]. Though profiles of ectopic pairing have shown a significant inter-strain variation, the pairing often occurs at the same loci. Thus, DNA sequence itself may define which loci can form a contact, whereas the epigenetic factors and/or activity of specific genes, such as *limk1*, affect FEC values by changing chromatin properties and tendency to pair.

There are several models of ectopic contact formation with the IH bands: (1) Pairing of "sticky ends" of the short repeated sequences within the areas of the DNA breaks that occur due to under-replication. (2) Pairing between the extended homologous DNA sequences. Generally, ectopic pairing occurs between the areas of chromosomes that do not show a significant homology, though it has been observed for some bands. (3) DNA branch migration upon replication mistakes, presumably due to the restricted homology between the associating sequences. (4) Pairing mediated by specific heterochromatin-associated proteins [9]. At least in the first three cases, presence of identical DNA sequences is crucial for pairing according to the principle of complementarity. In other words, high FECs should correlate with high frequencies of matching short DNA fragments (FMF) for the contacting regions.

The method of squashed preparations is imprecise, as it permits the localization of areas of contacts with a resolution of tens to hundreds of thousands bp. Hence, it gives no information about the specific DNA sequences involved in ectopic pairing. The bioinformatics approach helps to handle this problem. To estimate correlations between Drosophila FEC and FMF, we designed software called Homology Segment Analysis [23]. The current version of the software performs the following:

Steps 1–3 can be performed for sections within chromosome parts of different sizes.

Using a previous version of the software, we showed a positive correlation between *Berlin*/*agnts3* FEC and FMF for identical short (30–50 nt) DNA fragments. In that research, we specifically focused on the X:11AB region and its contacts with the other sections of the X chromosome. Most of the fragments found to putatively make impact into ectopic pairing were similar to the middle repetitive DNA 372-bp sequence and the 1.688 g/mL satellite DNA family [24]. The distribution and properties of 372-bp indicate its possible role in *Drosophila* dosage compensation and primary sex determination [25]. The 1.688 satellite DNA is abundant in *Drosophila* genome (2%), being localized both in heterochromatin and in euchromatin domains (1860 and 168 copies, respectively), mainly on the X chromosome. This satellite family includes 360/359-bp, 353-bp and 257-bp subfamilies [26]. The 1.688 satellite DNA is known to produce small RNAs that participate in the localization of male-specific lethal complex (MSL) on the X chromosome, increasing male survival [27]. The above point to a striking connection between the ectopic pairing and non-coding RNAdependent processes of dosage compensation. 359-bp also produces a long non-coding RNA that interacts with centromeres of all major chromosomes, participating in their mitotic segregation [28].

In this study, we calculated FEC-FMF correlations for all pairs of the X chromosome sections (1A–20F) in *CS* and *agnts3* strains at fragment lengths of 10–60 nt. Each section begins with an IH band, hence all of them are theoretically able to participate in ectopic pairing. The effect of the proximity of chromosome sections on FEC-FMF correlation value was estimated by analyzing the average correlation values within the chromosome parts of different sizes. For fragments specifically making the contribution to FMF, the biological nature was determined using NCBI Blast software. Most of them showed a high percentage of identity with 1.688 and 372-bp repeats, as well as related genes. These repeated sequences were either concentrated in genomic regions predisposed to ectopic pairing or governing pairing themselves, via DNA-DNA complementary binding or indirectly with the help of some unknown protein or RNA factors.

#### **2. Results**

#### *2.1. FEC-FMF Correlations for the Whole X Chromosome*

The average values of statistically significant FEC-FMF Spearman rho correlations (R) were calculated for all the pairs of sections of the X chromosome (Figure 1). R specific (RSP) varied within 0.18–0.3, corresponding to a rather weak positive FEC-FMF correlation [29]. RSP was slightly higher for *agnts3* compared to *CS*, especially at a fragment length (L) of 45–50 nt (L45-50), possibly due to a larger FEC in *agnts3*. Repeats exclusion did not significantly affect RSP, except for its small decrease at L25 in *CS*. For both *CS* and *agnts3* strains, there were no significant RSP differences from those calculated at L50.

**Figure 1.** R values at different fragment lengths. *X* axis: L (nt). *Y* axis: R (conventional units). Difference: # from *CS*, ˆ from the case with excluded repeats, \* from the case with unspecific correlation, shading – difference from R calculated at L50 (two-sided Mann–Whitney U-test; *p* < 0.05). Standard error of mean is shown. Here and below: SP—specific, SP\_NR—specific with repeats exclusion, UN—unspecific, UN\_NR—unspecific with repeats exclusion. Sampling number: for specific correlations, *n* = 6–20 (*CS*), 5–28 (*agnts3*); for unspecific correlations, *n* = 568–2235 (*CS*), 667–2987 (*agnts3*). Total number of R estimations: for specific correlations, *n* = 90 (*CS*), 95 (*agnts3*); for unspecific correlations, *n* = 10,620 (*CS*), 11,210 (*agnts3*).

To check the correlation specificity, R unspecific (RUN) was calculated as the average value of statistically significant rho correlations for all the different section pairs A and B. Such a shuffle of sections permits estimation of the "false" correlation between the inappropriate FEC and FMF values. RUN was nearly the same as RSP; thus, if rho correlation appears to be statistically significant by chance, its average value does not differ from that of the "true" correlation. RUN was higher for *agnts3* compared to *CS* at L10–40; however, there were no interstrain difference at a larger L. For both strains, RUN grew with the length of fragments, reaching maximum values at L50–55. The exclusion of short repeats reduced RUN values for small Ls (10–20). Thus, microsatellites seem to have a significant impact on unspecific correlations. At a larger L (45–60), repeats exclusion did not significantly affect RUN. For both strains, it was mostly impossible to distinguish between RSP and RUN. Thus, the average rho value rather weakly reflected the probability of ectopic pairing.

The picture was different for the proportions (P) of statistically significant correlations (their share of all correlations). P specific (PSP) showed a non-linear variation along with L growth (Figure 2). For *agnts3*, the first PSP maximum (about 0.3) was observable at L15 (i.e., 30% of all sections forming the ectopic contacts had a significant FEC-FMF correlation); then, PSP dropped to 0.12 at L30, returning to about 0.15 at L35. Finally, PSP dropped to 0.05 at L50, being nearly equal to the probability of finding an FEC-FMF correlation by

chance (*p* < 0.05). This corresponds to the virtually complete absence of section-to-section matching for fragments longer than 60 nt (FMF was zero for the most sections). For *CS*, the whole picture was the same, but PSP changes with L were not significant. Most importantly, we observed a striking interstrain PSP difference at L45–50: PSP value remained rather high for *CS* but not for *agnts3*. The fragment length of 50 nt has been previously shown as optimal to detect FEC-FMF correlations for the other *Drosophila* wild-type strain, *Berlin* [24]. Thus, at least in some strains, high FMF (L50) can serve as a predictor of ectopic pairing with a probability of about 20%.

**Figure 2.** P values at different fragment lengths. *X* axis: L (nt). *Y* axis: P (conventional units). Difference: # from *CS*, ˆ from the case with excluded repeats, \* from the case with unspecific correlation, shadingdifference from P calculated at L50 (Chi-square test; *p* < 0.05). Standard error of sample proportion is shown.

P unspecific (PUN) values showed a decrease along with the L increase, down to 0.05–0.07, which was close the theoretically expected *p* value of 0.05. For *CS*, PSP is significantly higher than PUN at L25–60. For *agnts3*, there were no significant differences between PSP and PUN, except at L25. Hence, *agnts3* strain clearly demonstrated less FEC-FMF correlation specificity compared to the wild-type strain. Repeats exclusion significantly decreased PUN. This also proves that short simple repeats abundant in *Drosophila* genome make a significant impact into unspecific FMF-FEC correlations. The exclusion of microsatellites decreased the interstrain PSP difference, though PSP remained significantly higher for *CS* at L30. Putatively, short repeats are parts of some longer DNA fragments that specifically have an impact on FEC-FMF correlation, therefore repeat filtration can simultaneously remove these fragments from samplings used to calculate FMF.

The observable P-L dependence seems to be the summary effect of the following trends: (1) Both PSP and PUN are high for short fragments that are randomly distributed within genome. (2) Both PSP and PUN are low for long fragments, but PUN drops more rapidly along with L increase. Exclusion of short repeats lowers the probability of FEC-FMF correlation, so we do not recommend it be performed for large Ls. High FMF calculated at L50 seems to indicate the higher probability of specific ectopic pairing for *CS* X chromosome sections. However, the share of statistically significant correlations remains rather low (about 20%). Thus, ectopic pairing is mostly governed by some other factors not taken into account in the FMF calculated according to the above scheme.

#### *2.2. FEC-FMF Correlation for Chromosomal Regions of the Different Length*

Ectopic contacts are mainly formed between sections that are not too far from each other (i.e., are separated by not more than 10–20 sections) [30]. To check the influence of the intersection distance (D) on R and P values, we performed FEC-FMF calculations for sections within chromosome zones of specific length D (for a more detailed description, see Section 4.2, Stage 9).

The average R values for different D are shown in Figure 3. In *CS*, RSP–RUN difference (RDIFF) was about 0.05–0.1 for small fragments (L10) and small D (10–25). However, in most cases it was insignificant because of the high R variance. For larger D, RDIFF was small but statistically significant, save for the highest D values obtained for very small samplings. The highest RDIFF could be observed at L30, with its maximum at D40–45. At L50, there was no significant *CS* RDIFF, save for a few D values. The picture was similar for *agnts3*, but its RDIFF was smaller at L30 and even became negative at L50 (D < 50). The interstrain difference was small in most cases, but at L30 *CS* had apparently lower RUN compared to the mutant strain. Generally, the highest RSP and RDIFF could be observed for rather small Ds (up to 40), confirming the hypothesis that ectopic contacts mainly occur between spatially close bands.

The exclusion of microsatellites affected the *CS* RDIFF, which approached zero at L30 and increased at L50 with most D values (Supplementary Materials, Figure S1). RDIFF (L50) remained maximal for small D. Hence, the highest FEC-FMF positive correlation can be observed within the short areas of the X chromosome. In *agnts3*, the area of D with positive RDIFF shrinked at L10 and increased at L50, probably due to the increase in correlation specificity after repeats filtration. However, the range of D with the negative RDIFF increased at L30. Hence, *agnts3* shows pronounced negative FEC-FMF correlations for relatively small Ds.

For the *CS* P values, we can see an obvious increase in PSP–PUN difference (PDIFF) along with increases in D and L values (Figure 4). At L50, occurrence of specific correlations became significantly higher for almost the whole range of D. The highest values of both PSP and PUN could be observed in *agnts3* at L10, but PDIFF was relatively small. There was also an increase in *agnts3* PDIFF along with D growth at L30, but it completely vanished at L50 and even became negative for the mean D values. At L50, the interstrain PSP difference was maximum and PUN was equal for both strains, being close to the theoretically predicted value of 0.05. This confirms results obtained for the X chromosome without division onto D zones. Thus, the most specific FEC-FMF correlation for the whole X chromosome can be observed in *CS* at L50.

The exclusion of microsatellites did not generally change the picture for *CS* PDIFF, the value and significance of which increased along with the growth of D and L (Figure S2). PUN was about 0.05, similar to that for the whole X chromosome (see Figure 2). PSP (L50) was somewhat lower compared to the case without repeats exclusion, probably because of simultaneous filtration of long fragments specifically making impact into FMF. For *agnts3* , PDIFF was mostly negative at L30. This may indicate some inverse relationship between FEC and FMF for that strain. At L50, the mutant PDIFF was negative for D < 30 and positive for D > 60, though in both cases its values were small. Thus, in contrast to *CS*, *agnts3* does not show a stable increase in PDIFF along with D and L growth.

**Figure 3.** The average R values for the chromosomal zones of different intersection distances (D). *X* axis: D—zone length (in sections). *Y* axis: R (conventional units). Standard error of mean is shown. The areas of D with statistical differences are shown by straight lines above the diagram: green—SP vs. UN, orange—SP (*CS*) vs. SP (*agnts3*), grey—UN (*CS*) vs. UN (*agnts3*) (two-sided Mann–Whitney U-test; *p* < 0.05, *n* = 120–D).

**Figure 4.** The average P values for the chromosomal zones of different intersection distances (D). *X* axis: D—zone length (in sections). *Y* axis: P (conventional units). Standard error of mean is shown. The areas of D with statistical differences are shown by straight lines above the diagram: green—SP vs. UN, orange—SP (*CS*) vs. SP (*agnts3*), grey—UN (*CS*) vs. UN (*agnts3*) (two-sided Mann–Whitney U-test; *p* < 0.05, *n* = 120–D).

#### *2.3. Sections Prone to Ectopic Pairing*

Some chromosomal areas are known to be prone to ectopic pairing. The total FEC number (FECTOT) is a number of ectopic contacts between a given section and all the other sections of the X chromosome. There was a strong positive FECTOT correlation for *CS* and *agnts3*: rho = 0.757 (*p* < 0.001, *n* = 119). After replacing the exact FEC with the indices of presence (1) or absence (0) of contacts, the interstrain correlation decreased but remained significant: rho = 0.540 (*p* < 0.001). Thus, in both strains, ectopic contacts are mainly formed by the same sections, and sections with a high FEC are usually localized at the same position (see also [30]).

According to [31], there are five chromatin features (F) increasing the probability of the *Drosophila* X chromosome ectopic pairing, such as *Dm225* and *Dm234b* genes hybridization sites, ectopic conjugation, weak points, late replication and giant palindromes. For each F, we assigned the value of 1 or 0 to sections depending on whether they had or did not have a specific F. Then we calculated the F index (Ind) as the sum of values for each section (FSUM) divided by 10. As expected, there was a positive correlation between FECTOT and FSUM: rho (*CS*) = 0.357 (*p* < 0.001); rho (*agnts3*) = 0.252 (*p* < 0.01). Replacing the exact values with indices of presence (1) or absence (0) made correlation insignificant.

At the same time, we did not observe any correlation between FECTOT and the total FMF number for each section (FMFTOT) calculated at L30 or L50. FEC-FMF correlation was specifically observed only for a set of the X chromosome sections (Figure 5). Three of them (10B, 11D and 18B) coincided for both *Drosophila* strains. Only the minor part of sections (26.1% for *CS*, and 18.2% for *agnts3*) had chromatin features predisposing them to ectopic pairing. Considering a rather moderate FECTOT-FSUM correlation, we state that the indicated specific chromatin properties can facilitate ectopic pairing, but are not necessary for that. For most sections, ectopic pairing seems to be governed by some other mechanisms possibly related to specific DNA sequences within the interacting areas.

**Figure 5.** X chromosome sections showing statistically significant FEC-FMF correlations. Ind—the index of chromatin features (F). Y value: RCORR (for *CS* and *agnts3*) and Ind values (C.U.).

#### *2.4. The Biological Nature of Sequences Making Impact into Ectopic Pairing*

To reveal the molecular nature of the fragments most contributing to FMF, we analyzed their sequences using NCBI Blast. The first fifty L30 and L50 fragments with the maximum number of occurrences (NOs) were assayed, both for all sections (set I) and for a set of sections showing *CS*- or *agnts3*-specific FEC-FMF correlations (sets II and III, respectively). To simplify the analysis, we did not consider fragments with NOs ≤ 10 and considered only one of the fragments with the equal NO values.

BLAST analysis made it possible to divide all matching fragments into six classes: (1) Microsatellites, such as an, tn, (at)n, (aat)n, (gata)n, (agata)n, (tcccag)n and so on, the maximum motif length being of six. (2) Fragments showing a high percentage identity (90–100%) with 1.688 and 372-bp repeats, as well as with some genes such as *c11.1*, *kl22 Drak* intron and related. (3) Fragments showing a high percentage identity with sequences of *c11.1* and *kl22 Drak* but not with 1.688 or 372-bp repeats. (4) Fragments of transposon HB1. (5) Fragments of long non-coding RNA genes and retrotransposon roo-900. (6) Other; mostly the sequences with an unknown molecular nature, or those for which BLAST returned no result. Up to 31% of fragments (with NOs < 50) remained non-annotated. The results are present in Figure 6. The full list of fragments, along with their NO, molecular nature and BLAST sequence identity, is shown in Table S1.

**Figure 6.** Short DNA fragments that contribute to FMF. Difference: # from I, non-annotated sequences included when calculating the proportions; \$ from I, non-annotated sequences excluded (Chi-square test of two proportions; *p* < 0.05, *n* > 200).

In set I, microsatellites constituted the major class of fragments. The above is not surprising, as the program searches fragments one by one, with a step of one nucleotide, generating a large number of identical short sequences from an extended area of repeats. At the same time, the simplest repeats, such as (at)n and (ta)n, prevailed at L30, but not at L50. The microsatellite (tcccag)n, appeared to have a maximum NO at L50 and a third-rank NO at L30, indicating its abundance in *Drosophila* genome. The NO was also high for fragments belonging to class 3, while class 2 was nearly absent at L50. Transposon parts constituted at least 3% of L50 fragments.

In set II, containing fragments of sections with *CS*-specific FEC-FMF correlations, NOs for microsatellites significantly decreased, especially at L50. At the same time, NOs for 1.688/372-bp-related sequences increased many times (up to 72% at L50). PSP was at its maximum for *CS* at L50 (see Figures 2 and 4). This clearly indicates that microsatellites mostly contribute to unspecific FMF, and 1.688/372-bp-related sequences mostly contribute to FEC-FMF specific correlation. The latter is in agreement with our previous data obtained for *Berlin* [24]. (at)n was the only microsatellite that seemed to have an impact on *CS*specific correlation. The fact that (at)n and complementary sequences constituted 66% of microsatellites and 31% of all *CS*-specific fragments at L30 explains why the NO was high for *CS*-specific sets of selected microsatellite fragments. Sequences of class 3 also contributed to FEC-FMF correlation at L30.

On the contrary, for *agnts3*-specific fragments (set III) we observed an increase in microsatellite NOs, which reached their maximum (at)n at L30 and (gata)n at L50. For classes 2 and 3, the NOs insignificantly grew at L30, although for class 3 the NO dropped at L50. In *agnts3*, PSP was only slightly higher than PUN at L30, with no difference at L50 (see Figures 2 and 4). The positions of sections with specific FEC-FMF correlations were different for *CS* and *agnts3* (Figure 5). Taken together, our data show that 1.688/372-bprelated sequences are associated with ectopic pairing in *CS*, with either a minor or no association for *agnts3* .

#### **3. Discussion**

Ectopic pairing is a long-range interaction that occurs with a relatively low frequency between the IH bands of the *Drosophila* polytene chromosomes. Though ectopic contacts can be easily observed using light microscopy, their molecular nature and functional role in the nuclear 3D organization remain obscure. The *agnts3* strain was shown to have a significantly higher FEC compared to the wild-type fly strains. At the same time, some IH bands are more prone to form ectopic contacts compared to the others. In *Drosophila* reciprocal hybrids of *Berlin* and *agnts3*, specific bands of the X chromosome demonstrate either matroclinic or patroclinic inheritance of high FEC values [24]. Thus, ectopic pairing obviously has both genetic and epigenetic basis.

As molecular processes of ectopic pairing are not studied in detail, we cannot say a priori how long a zone of local DNA pairing should be to initiate a contact formation. Neither can we say whether g/c rich sequences will interact preferably over a/t rich due to the higher binding energy or if the tendency will be reversed, as a/t regions are often nucleosome free, easy to melt, frequent in genome, and therefore may be more prone to misalign. In our study, we have made a simple assumption that all cases of local DNA match will increase the probability of the ectopic contact formation. The optimal length of the areas should be the tread off the specificity of matching and the probability to find a relatively long specific area in genome. We have revealed the optimal fragment length to be about 50 nt. The ability to consider cases of incomplete matching would let us work with longer fragments, but the current version of Homology Segment Analysis software does not permit us to work with local mismatches. FMF values do not depend on the relative orientation of fragments (e.g., they are the same for the sequential parts of long (at)n repeats and several short non-overlapping (at)n repeats). The fragments of all chromosome sections were considered, though only the IH bands participate in ectopic pairing. We believe that the above has a rather small effect on FMF, as the IH bands are densely packed, containing much more DNA than interbands. FEC was zero for most of the section pairs, and the range of FEC variation was rather narrow: in most cases, the FEC equaled 1 or 2 in its Spearman rank correlation. In addition to DNA sequence, epigenetic factors greatly influence FEC values.

As a result, we obtained a rather moderate average FEC-FMF Spearman correlation rho value (R of 0.3), as well as the proportion of section pairs for which the correlation was significant (P of 0.2). Nevertheless, for *CS*, P appeared to be higher than that of unspecific correlations, both for the whole X chromosome and for almost the entire range of the lengths of its parts. Moreover, we found that not all DNA sequences, but mainly those related to 1.688 or 372-bp repeats, contribute to FEC-FMF correlations. This cannot be explained by the predominance of such sequences in the *Drosophila* X chromosome, as, according to our data, the most common matching repeats were (at)n at L30 and (tcccag)n at L50 (Table S1). The other microsatellites also contributed to FMF, but only (at)n seemed to make an impact on the *CS*-specific FEC-FMF correlation. In *agnts3*, FEC-FMF correlations are generally unspecific, despite the higher FEC for this strain. Similarly, the X—X:11AB FEC-FMF correlation was lower in *agnts3* compared to *Berlin* [24]. At L30, specific correlations still could be observed for about 10% of bands participating in *agnts3* ectopic pairing. Microsatellites mostly contributed to *agnts3*-specific sets of matching fragments.

Satellites are multi-copy tandem DNA repeats classified according to the length of their monomers: microsatellites (1–10 bp), minisatellites (10–100 bp) and satellites (>100 bp). *Drosophila* satellites are involved in multiple cell processes, such as dosage compensation, heterochromatin establishment, gene activity regulation, maintenance of genomic architecture, chromosome segregation and development [32]. Satellite DNA is one of the most abundant and fast evolving components of genome, being the major part of the constitutive heterochromatin in eukaryotes. There are at least 14 families of highly repeated *D. melanogaster* satellite DNA [33]. The most abundant satellite repeats of the X chromosome are aatat, aagag, 359/372/260-bp, and IGS [34].

The 372-bp satellite is a conservative a/t rich *Drosophila* repeat concentrated on the euchromatin of the X chromosome, located mainly between cytogenetic regions 4 and 15, with about 300–400 copies per haploid genome. Its distribution and sequence features suggest its participation in the primary fly sex determination and dosage compensation [25]. This satellite DNA is homologous to the 1.688 g/mL class of satellites predominantly localized to the centromere heterochromatin of the X chromosome [34,35]. The 1.688- 3F satellite expresses a siRNA-generating hairpin dsRNA that increases males survival, regulating the male-specific lethal complex (MSL) positioning on the X chromosome. It is possible that siRNA affects dosage compensation by modifying chromatin at the 1.688 repeat [27]. It should be noted that the *agnostic* locus in *agnts3* does not show dosage compensation [36,37]. Along with a lack of FEC-FMF correlation, this may indicate some deregulation in 1.688 or 327-bp activity in *agnts3*. The 1.688 sequence shows a significant intraspecific nucleotide divergence: 10% for heterochromatin and 27% for euchromatin. Hence, the chromatin structure seems to influence the rate of 1.688 evolution [26]. In our study, the identity between the matching fragments with strain-specific FEC-FRF correlations and 1.688 satellites from 3C and 10Ep bands [35] was about 100% at L30 and slightly below 100% at L50.

Polytene chromosomes normally have low flexibility, extending across the nucleus in a Rabl orientation, with no stable interaction between loci 1-2 cytological divisions apart [13]. The dysfunction of the two interband-associated proteins, Chromator and JIL-1 kinase, leads to dramatic impairment of polytene chromosome morphology (i.e., band misalignment, curling and numerous ectopic contact formations) [38]. As the pattern of ectopic contacts is determined early in development, it should rather reflect the morphology of diploid embryonic nuclei with long-range interactions (e.g., between Polycomb-repressed domains) [16]. Thus, the functional role of ectopic pairing in chromatin spatial organization is questionable. At the same time, it seems to represent some aspects of nuclear 3D structures typical of the early stages of fly development.

The molecular processes that govern ectopic pairing are still poorly understood. In the *Drosophila* interphase nucleus, blocks of heterochromatin tend to associate with each other. The association does not require similar sequences such as (aagag)n satellite to be located within the contacting areas. Presumably, their associations are mediated by proteins that recognize the general features of heterochromatin, such as specific histone modifications, repetitiveness, late replication and low activity [39]. Nevertheless, the restricted homology at the IH areas containing DNA breaks may be important for ectopic contact formations. Probably, it results from the joining and ligation of truncated DNA ends between the IH bands. The level of suppressor of under-replication *(SuUR)* gene expression positively correlates with FEC. *SuUR* overexpression enhances the IH under-replication and ectopic pairing only before the third instar larval stage [17].

In accordance with the above, FEC as a phenotypic trait is determined at the embryonic stage: the high temperature (37 ◦C) applied at that stage increases FEC in *CS* without an effect in *agnts3*. This seems to be connected to strain-specific properties of heterochromatin that begins to form at this stage [22]. High FEC values in *agnts3* may reflect the increase in strain-specific recombination/reparation activity, as well as the activity of some chromatin proteins, such as the Polycomb group, HP1, Chromator and JIL-1. It is interesting to note that both JIL and *agnts3* LIMK1 genes bring the insertion of a mobile S-element that may

theoretically cause some interaction between these genes [21,40]. *agnts3* is also characterized by significant changes in miRNA expression profile compared to the wild-type strains *CS* and *Berlin* [21,24]. The mutant-specific heterochromatin properties may affect the wide profile expressions of genes, including non-coding RNA genes.

On the contrary, miRNAs may affect the expression of heterochromatin proteins and the tendency of the IH bands to form ectopic contacts. The targets of some of miRNAs participating in the development of human neurodegenerative disorders are Swi/Snf-like chromatin remodeling complex, REST factor that recruits histone deacetylases (HDAC), and SIRT1, a NAD+ -dependent HDAC involved in heterochromatin formation. miR-124 and miR-34c negatively regulate HDAC1/2 and SIRT1, respectively [41]. *Drosophila* miR-124 and miR-34 are decreased in *agnts3* compared to *CS* [21]. Hence, the high percentage of heterochromatin in *agnts3* may be due to the increase in HDAC activity.

*agnts3* is shown to impair the activity of LIMK1, the main regulator of actin remodeling [18]. The filamentous and globular actin differently affect chromatin conformation, as well the activity of HDAC [42]. Actin is widely involved in the regulation of genetic apparatus, including transcription machinery [20]. Active forgetting also depends on an LIMK1-dependent signaling cascade [43]. This reveals a possible connection between the molecular processes of heterochromatin formation in the early embryogenesis of fruit flies and the *Drosophila* ability to learn, memorize and forget.

Heterochromatin proteins such as the Polycomb group and HP1 may influence ectopic pairing by bringing the IH bands closer together in space or linking them to the envelope [17]. HP1 promotes the formation of chromosome loops, facilitating the coalescence of dispersed middle repeats such as micropia retrotransposon and non-coding RNA gene αγ [44]. H3K9me3 modification of 1.688 satellite creates a binding site for HP1 [45,46]. Hence, FEC-FMF correlation can be theoretically explained by HP1-dependent juxtaposition of the 1.688-containing under-replicated IP bands, followed by ligation of the double-stranded DNA ends. A similar role of 1.688-3C in ectopic pairing of the spermatocyte X chromosome was proposed in [34,47]. Previously, 1.688 was proposed to influence chromatin architecture by interacting with proteins of the nuclear matrix, such as Topoisomerase II and satellite binding protein. The 1.688 satellite may also participate in long-range interactions regulated by siRNA-dependent chromatin modifications [48]. Though 372-bp repeats are localized to euchromatin, they may also play a role in bringing together chromosome sections that form ectopic contacts.

In summary, our computational data confirm the hypothesis that *Drosophila* satellite DNA such as 1.688 and related sequences, can participate in long-range interactions between the IH bands. The lower FEC-FMF correlation for *agnts3* relative to *CS* may indicate less specificity of ectopic pairing in the mutant strain, similar to Chromator/JIL-1 mutants. Some proteins or non-coding RNAs, possibly produced by 1.688-like repeats, can mediate the interaction. Alternatively, they can affect heterochromatin properties and/or DNA replication within IH bands, making them more prone to pairing. Further studies are necessary to prove this experimentally.

#### **4. Materials and Methods**

#### *4.1. FEC Matrices*

FEC matrices for the X chromosome sections of *CS* and *agnts3* strains were taken from [30]. Files containing FE*Cs* for all X chromosome sections of both strains were supplied with Homology Segment Analysis software. Each matrix was built on the data obtained by orcein staining of squashed preparations of *Drosophila* 3rd instar larvae. For each strain, 30 larvae were taken, and about 20 nuclei were examined for each larva. FE*Cs* were calculated as the total number of ectopic contacts between the given section pair. The example of an ectopic contact is shown in Figure S3.

#### *4.2. Bioinformatics Analysis*

FEC-FMF correlations were estimated using Homology Segment Analysis software (Zhuravlev A.V., Zakharov G.A.; Pavlov Institute of Physiology, Saint Petersburg, Russia). The program is written on Python3. The calculations in the paper were performed using the latest version of the program, freely available at [23]. For both strains, we used the *D. melanogaster* X chromosome sequence genome assembly Release 6 (dm6) [49]. As there are no full genomic sequences for these strains, and S-elements are relatively short and can be found at various positions in different *Drosophila* strains, we do not consider it justified to include *limk1* S-element in *agnts3* genome sequence to compute its FMF. All interstrain differences are expressed here as FEC differences. The borders of the X sections were chosen according to Flybase data (www.flybase.org; accessed on 20 February 2021): 1A–20F, except 20B missing some nucleotides, totally 119 sections.

To install and run all the scripts, see the Readme file in the main directory. The computational algorithm is as follows:

For all B except A: FMF(A) is a set of FMF(A-B); FEC(A) is a set of FEC(A-B). FEC(A-B) is a value obtained from FEC matrix for a given *Drosophila* strain, being the number of contacts between sections A and B.

	- a. For each section A: calculation of FEC(A)-FMF(A) Spearman rank correlation coefficient (rho; *p* < 0.05, *n* = 119).
	- b. Calculation of the average rho value (R) and proportion of statistically significant FEC-FMF cases (P). P is calculated as follows: P = *n* (*p* < 0.05)/total *n* of estimations for which Spearman correlation data were obtained.
	- a. For each D: a sequential selection of the X chromosome zones (Z) of D length, with a step of 1 section (e.g., for D = 30, 90 different Z are selected, starting from 1 (Z1, or 1A–5F) and up to 90 (Z90, or 15F–20F)). The section notations A–F are equivalent to 1–6, so Z90 of D30 is also denoted as 156–206.
	- b. For each Z(D): specific and unspecific R and P calculation, as described in Stages 5–8, taking into account only the sections within Z. Currently, analysis of unspecific correlations is time consuming, taking up about 3 h for each L. So some cases (e.g., with specific L values or excluded repeats) may be omitted to speed up the processes. For each D, R(D) and P(D) values constitute samplings for further statistical analysis. The sampling size *n* is equal to the number of Z(D): *n* = 120–D.

The scheme of the Stages 1–9. is also given in Figure S4.

#### *4.3. Statistical Analysis*

All analyses were performed using scripts included in the Homology Segment Analysis software package.

#### *4.4. BLAST Analysis of Fragments Contributing to FMF*

#### **5. Conclusions**

Using our Homologous Segment Analysis software, we have shown a specific positive correlation between FEC and FMF for about 20% of the *CS* X chromosome sections involved in ectopic pairing. Most of the 50 nt fragments specifically contributing to FMF appeared to be related to 372-bp or 1.688 middle repeats. Thus, our bioinformatics approach lets us to handle the problem caused by the low resolution of the method of squashed preparations, which does not give information about the specific sequences involved in ectopic pairing. Using the experimental data on chromatin properties obtained by Hi-C and other 3Crelated methods with significantly higher resolutions can substantially increase correlation values and validity. Moreover, Homology Segment Analysis can be easily applied to search correlations between FMF and every feature associated with pairs of genomic regions both in *Drosophila* and in other species. For example, it can be used to find DNA motifs involved in contact formations, as well as the binding of proteins or RNA that mediate such interactions and, thereby, define nuclear 3D organization.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/ijms22168713/s1.

**Author Contributions:** Conceptualization, A.V.Z.; Methodology, A.V.Z.; Software, A.V.Z. and G.A.Z.; Validation, A.V.Z.; Investigation, A.V.Z. and E.V.A.; Resources, E.V.S.-P.; Data Curation, A.V.Z. and E.V.A.; Writing—Original Draft Preparation, A.V.Z., E.V.A., A.V.M., G.A.Z. and E.V.S.-P.; Writing— Review & Editing, A.V.Z., E.A.N., A.V.M., G.A.Z. and E.V.S.-P.; Visualization, A.V.Z. and G.A.Z.; Supervision, E.V.S.-P. and E.A.N.; Project Administration, E.V.S.-P.; Funding Acquisition, E.V.S.-P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the State Program of the Russian Federation 47 GP "Scientific and Technological Development of the Russian Federation" (2019-2030) (63.1) and the Russian Foundation for Basic Research (N 20-015-00300 A).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The FEC data used in this paper are freely available at: https://drive. google.com/drive/folders/1i\_SOIR3cxFN1951akkAMGcOgXwd71NG5?usp=sharing (accessed on 12 August 2021).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

### *Article* **Cellular Computational Logic Using Toehold Switches**

**Seungdo Choi, Geonhu Lee and Jongmin Kim \***

Department of Life Sciences, Pohang University of Science and Technology, 77 Cheongam-ro, Pohang 37673, Gyeongbuk, Korea<EMAIL_ADDRESS>(S.C<EMAIL_ADDRESS>(G.L.) **\*** Correspondence<EMAIL_ADDRESS>Tel.: +82-54-279-2322

**Abstract:** The development of computational logic that carries programmable and predictable features is one of the key requirements for next-generation synthetic biological devices. Despite considerable progress, the construction of synthetic biological arithmetic logic units presents numerous challenges. In this paper, utilizing the unique advantages of RNA molecules in building complex logic circuits in the cellular environment, we demonstrate the RNA-only bitwise logical operation of XOR gates and basic arithmetic operations, including a half adder, a half subtractor, and a Feynman gate, in *Escherichia coli*. Specifically, de-novo-designed riboregulators, known as toehold switches, were concatenated to enhance the functionality of an OR gate, and a previously utilized antisense RNA strategy was further optimized to construct orthogonal NIMPLY gates. These optimized synthetic logic gates were able to be seamlessly integrated to achieve final arithmetic operations on small molecule inputs in cells. Toehold-switch-based ribocomputing devices may provide a fundamental basis for synthetic RNA-based arithmetic logic units or higher-order systems in cells.

**Keywords:** toehold switch; arithmetic operation; RNA–RNA interaction; molecular computing; reversible computing

#### **1. Introduction**

Synthetic biology aims to create technologies for designing and building biological systems with programmable and predictable dynamics [1]. Since the demonstration of synthetic biological circuits in living cells over two decades ago [2,3], considerable progress has been made towards more sophisticated artificial cellular functions, such as feedback oscillation [4], combinatorial logic computation [5–7], and information storage [8,9]. In principle, synthetic circuits can be constructed using any biological molecule as a backbone. Natural and engineered protein regulators can provide the framework to implement logic circuits and computing devices [10,11], including de-novo-designed components [12]. Still, the construction of synthetic biological circuits presents numerous challenges, including the lack of composability [13], limited modularity [14], unpredictable cross-reactivity [15], cellular resource usage [16,17], and unexpected idiosyncratic behavior in real-world applications [18]. Therefore, novel approaches for synthetic biological circuits, including the development of readily characterized, standardized, and modular components are required to overcome the innate difficulties in managing and programming cellular behavior towards large, complex synthetic systems [19].

Nucleic-acid-based genetic devices have made remarkable progress in molecular computing and may provide the required platform for scalable synthetic biological systems. Complex logic circuits and advanced computing systems have been implemented using toehold-mediated strand displacement, including a bistable circuit [20], a square-root circuit [21], neural networks for memory [22] and pattern recognition [23], and an arithmetic logic unit [24]. Furthermore, these molecular computing systems are amenable to computational design and analysis [25,26]. While DNA strand-displacement circuits have been demonstrated in live cells [27] and for live-cell imaging [28], DNA logic gates are generally not suitable for in vivo applications due to the challenges in generating single-stranded

**Citation:** Choi, S.; Lee, G.; Kim, J. Cellular Computational Logic Using Toehold Switches. *Int. J. Mol. Sci.* **2022**, *23*, 4265. https://doi.org/ 10.3390/ijms23084265

Academic Editors: Vaclav Brazda and Richard Bowater

Received: 8 March 2022 Accepted: 10 April 2022 Published: 12 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

DNAs [29] and maintaining stable gate configurations [30] in cells. In comparison, RNA provides several advantageous characteristics for synthetic biological computing devices. RNA strands can be easily programmed to interact in a designed fashion due to the singlestranded nature of RNA. The co-transcriptional folding of RNA allows the formation of stable secondary structures [31], suitable for natural and synthetic riboswitches [32]. Furthermore, RNA signals can be easily modulated in a tunable manner [33] and can also be amplified with several RNA and protein counterparts [34–38]. Therefore, RNA has been exploited as a platform to engineer gene expression programs that operate robustly in vivo [39–44].

Building on the success of synthetic RNA regulatory parts and inspired by natural RNA regulators [45], several de-novo-designed RNA regulators have been utilized for synthetic biological devices with a large library of well-characterized parts [40,41]. As an example, toehold switches control gene expression in trans via well-established Watson-Crick base pairing of switch and trigger RNA molecules (Figure 1a) [40]. Unlike conventional riboregulators [39], the toehold switches remove nearly all the sequence constraints, exhibit a wide dynamic range, and show excellent programmability with a large library of orthogonal parts. The versatility of toehold switches for synthetic genetic circuit construction is exemplified by the recent developments in cellular logic computation [46], translational repressing riboregulators [47], incoherent feed-forward loop circuits [48], synthetic transcription terminators [49], the protein quality control system [50], modulators of riboswitch circuits [51], and regulators of mammalian cells [52]. Beyond cellular circuits, the toehold switches find use in other platforms, such as cell-free systems [53–57] and paper-based diagnostic devices [58–60] for broader applications.

**Figure 1.** De-novo-designed toehold switch and toehold-switch-based NIMPLY gate. (**a**) Scheme of toehold switch operation. The toehold switch has repressed the translation state through the secondary structure sequestering the RBS and start codon. Linear-linear interaction between the toehold switch and the trigger RNA exposes the RBS and start codon with the strand displacement process; therefore, the translation of the downstream gene is resumed. (**b**) Scheme of NIMPLY gate operation. Antisense RNA has extended overhang sequences at both ends and can inhibit the trigger RNA through sequence displacement or complementary binding. Thus, the toehold switch is reverted to the initial OFF state.

In particular, toehold switches may form the basis for constructing an arithmetic logic unit (ALU) in vivo. A generalized toehold switch architecture, termed ribocomputing devices, concatenated several toehold switch sensor domains and utilized the self-assembly of RNA species to compute multi-input AND/OR/NOT operations [46]. The design

flexibility of toehold switches, if effectively utilized, can lead to the streamlined design and construction of a basic form of synthetic biological ALU. Previous work has demonstrated basic ALUs, including a half adder and a half subtractor, in bacteria and mammalian cells [5,61,62]. These binary calculators can perform bitwise calculations across two 1-bit input signals and serve as building blocks for higher-level systems. A half adder takes two bits of information and generates two output bits: one for the sum and one for the carry. The sum bit can be calculated via an XOR gate, and the carry bit can be calculated via an AND gate. Thus, it is straightforward to calculate the carry bit, as previously shown, but a toehold-switch-based XOR gate needs to be engineered. Unlike previous XOR gate implementation, for example, in vitro [63], in prokaryotes [61,62,64,65], and in eukaryotes [66–68], a toehold-switch-based XOR gate is an RNA-only synthetic logic device. Building on previous work in which an antisense RNA that titrates a cognate trigger RNA molecule can be used for a NIMPLY gate operation (Figure 1b), a toehold-switch-based XOR gate with a compact architecture can be obtained by concatenating two orthogonal switches in an OR-gate fashion. A half subtractor can be analogously constructed with an XOR gate and a NIMPLY gate. In addition, a Feynman gate, one of the reversible logic gates that map input and output signals in a one-to-one manner [69], can be obtained using an XOR gate and a BUFFER gate. In summary, we present the binary operation of an XOR gate, cellular arithmetic operations with a half adder and a half subtractor, and a Feynman gate in *E. coli* using a de-novo-designed toehold switch and antisense RNAs. Synthetic RNA-based ALUs could lay the foundation for making sophisticated molecular devices with neural-network-like capabilities for biomedical applications.

#### **2. Results**

*2.1. XOR Gate*

#### 2.1.1. Design of XOR Gate with Toehold Switches

The NIMPLY gate often used in synthetic biology and genetic circuits [61,70] was previously demonstrated using toehold switches [46]. A NIMPLY B is equivalent to A AND (NOT B), and an XOR gate can be constructed using two NIMPLY gates connected via an OR gate as follows: A XOR B = (A NIMPLY B) OR (B NIMPLY A). Thus, we sought to first demonstrate two orthogonal NIMPLY gates. The mechanism for NIMPLY gates is analogous to previous work where the switch RNA is activated by the trigger RNA (A), and the antisense RNA (B) deactivates the trigger RNA via direct hybridization or toehold-mediated strand displacement to separate the trigger RNA bound to the switch RNA. The extended overhang sequences at both ends provide the thermodynamic driving force to shift the equilibrium towards trigger and antisense RNA binding rather than trigger and switch RNA binding.

To implement a NIMPLY gate in *E. coli*, the three circuit components—switch RNA, trigger RNA, and antisense RNA—should be selected from the existing library with appropriate modifications (Figure 2a). We selected a pair of second-generation toehold switches with large dynamic range and strong orthogonality. These two switches are connected with a 9-nt linker sequence to create an OR gate, as previously demonstrated [46]. The overhang sequences of both trigger RNAs and antisense RNAs were designed via the RNA secondary structure prediction software NUPACK [71–75] (Table S1). Fifteen nucleotide overhangs were attached to both ends of the trigger and antisense RNAs, and a single nucleotide bulge was inserted between the overhang region and the switch binding domain to prevent the formation of long double-stranded RNA that could be targeted for degradation by RNase III [76,77]. The design candidates were analyzed for ensemble defect [78], overhang accessibility, and crosstalk in silico to select the best designs to be tested in experiments.

**Figure 2.** Design of combined NIMPLY gates and optimization strategies. (**a**) Scheme of combined NIMPLY gates composed of two orthogonal NIMPLY gates connected by a 9-nt linker. Trigger and antisense RNAs control the translation states of NIMPLY gates. (**b**) Performance of initial design for combined NIMPLY gates. T1 and T2 denote trigger RNAs that activate NIMPLY gates 1 and 2, respectively, and A1 and A2 denote antisense RNAs that annihilate trigger RNAs 1 and 2, respectively. (**c**,**d**) Effect of the location of extended overhangs on trigger and antisense RNAs. TR and AR indicate trigger RNA and antisense RNA, respectively. Absence of trigger (TR = 0) indicates that only the NIMPLY gate RNA is present. (**e**) Effect of the presence of bulge on antisense RNA. TR and AR indicate trigger RNA and antisense RNA, respectively. Full match means that no bulge was introduced in the antisense RNA design. (**f**) Flow cytometry GFP fluorescence histograms for the NIMPLY complex with full-match antisense RNAs. T1, T2, A1, and A2 indicate trigger and antisense RNAs for switches 1 and 2 as above, and d represents a decoy RNA that does not interact with the switch RNA. T7 RNA polymerase was induced by 1 mM IPTG in *E. coli* BL21 DE3 strain. GFP fluorescence was measured on the flow cytometry (error bars indicate ± SD from three biological replicates). Cellular autofluorescence was subtracted in all cases. Autofluorescence level was measured from cells not bearing a GFP-expressing plasmid.

The NIMPLY gates were tested in *E. coli* BL21 DE3 strain with the switch, trigger, and antisense RNAs expressed from separate low, medium, and high copy plasmids, respectively. The RNA components were under the control of a T7 promoter, and genomically encoded T7 RNA polymerase was induced by Isopropyl β-d-1-thiogalactopyranoside (IPTG). GFP was used to characterize the switch output performance via flow cytometry. First, the OR gate functionality was verified in the absence of antisense RNAs, where the GFP fluorescence was increased at least 100-fold in the presence of either trigger RNAs (Figure S1). Next, the switch, trigger, and antisense RNAs were expressed together in cells to evaluate the function of the NIMPLY gates. The design with the least expected intramolecular and intermolecular structures showed the best performance among the candidates (Figure 2b). The performance was evaluated by dividing the ON state, with a cognate trigger RNA and a non-cognate antisense RNA, by the repressed state, with cognate trigger and antisense RNA pairs. Consequently, we observed increases of 11.2-fold and 43.3-fold for T1-A1 and T2-A2, respectively. Other design candidates with expected secondary structures within trigger RNAs showed relatively poor functionality (Figure S2). Since the design variants on previous NIMPLY gate designs [46] were not extensively characterized, we aimed to further explore and optimize the design choices to enhance the functionality of the NIMPLY gates and hence the performance of the synthetic XOR gate.

#### 2.1.2. Optimization Strategies for Toehold-Switch-Based XOR Gate

For the design variants of the trigger and antisense RNAs within the NIMPLY gate, we mainly adjusted the location of the extended overhangs and the presence of bulge. First, we investigated whether the location of the overhang could affect the functionality of the trigger or antisense RNAs. The GFP fluorescence output for trigger RNA 2 with overhangs showed a stronger reduction than the trigger RNA 2 without the overhang sequences (Figure S3). We hypothesized that the close proximity of RBS within switch 2 and the 50 overhang of trigger RNA 2 affect the access of RBS through steric hindrance. Therefore, trigger RNA variants with only a 50 overhang or a 30 overhang were constructed, and the impact of the overhang location on the switch performance was investigated. Trigger RNAs with only a 50 extended overhang showed weak repression by antisense RNA for switch 1 and weak activation for switch 2, indicating that the 50 extended overhang could reduce performance. On the other hand, trigger RNAs with only a 30 extended overhang showed improved performance for both switches compared to the trigger RNAs with both 50 and 3 0 overhang domains (Figure 2c). The antisense RNAs were similarly modified to test the impact of overhang domains: antisense RNAs with only a 5' overhang showed improved fold repression, while antisense RNAs with only a 3' overhang showed weak repression comparable to the antisense RNAs without the extended overhangs (Figure 2d).

Other design candidates were analyzed for the impact of overhang locations on the trigger and antisense RNAs, and a similar trend was observed (Figure S4). Although the 50 extended overhang can be considered disposable, simply removing the existing 50 overhang caused crosstalk in some cases because it was not considered during the design phase (Figure S5). Fortunately, the apparent crosstalk was negligible when the 3' overhang trigger RNAs were paired with antisense RNAs with both overhang sequences (Figure S6). Additionally, an expanded hairpin loop of the switch RNA was explored to help reduce the potential steric hindrance of trigger RNA on the RBS. The increased hairpin size in switch RNA increased the ON level expression but also generally increased the OFF-state leakage (Figure S7). Together, we observed the impact of extended overhang locations on both trigger and antisense RNAs and trade-offs in the switch RNA hairpin loop size on the performance of NIMPLY gates.

To further enhance the functionality of the NIMPLY gates, we investigated the effect of bulges within the trigger–antisense RNA complex on the repression efficiency. The antisense RNA presumably works in one of two ways: (1) dissociating the trigger from the switch or (2) capturing the trigger freely floating in the cytoplasm [46]. Single nucleotide bulges located between the overhang and switch binding domain can act as an energy barrier to the strand displacement pathway that removes the trigger from the switch [79,80]; in that case, the direct hybridization mechanism would be predominant. In order to increase

the repression efficiency of antisense RNA, an antisense RNA with extended overhangs but without bulges was designed and tested. The repression efficiency was enhanced nearly 10-fold on trigger 1 (Figure 2e), and the combined NIMPLY gates with optimization exhibited 48.5-fold and 65.6-fold improvements, respectively (Figure 2f). Therefore, we successfully constructed two orthogonal NIMPLY gates with large dynamic ranges using optimization strategies on switch, trigger, and antisense RNA designs. These may serve as useful strategies for other toehold-switch-based logic circuit designs and potentially for other RNA regulatory devices as well.

#### *2.2. In Vivo Characterization of XOR Gate*

Encouraged by the optimized NIMPLY gates, we then aimed to construct an XOR gate with two chemical inducers as inputs: IPTG and anhydrotetracycline (aTc). An XOR gate provides a true output with an odd number of true inputs (Figure 3a). In the case of the chemically inducible XOR gate that we aimed to construct, the GFP output should be high when either IPTG or aTc is present, but not both. To achieve this, Lac and Tet operators were strategically placed downstream of T7 promoters that drive the expression of trigger and antisense RNAs, such that the trigger RNA of one NIMPLY gate and the antisense RNA of the other NIMPLY gate are simultaneously induced by the same chemical inducer for both inducers (Figure 3b). Specifically, an IPTG induction promotes the expression of trigger RNA 1 and antisense RNA 2, such that the output GFP expression is high. The process works similarly for aTc induction. However, the simultaneous treatment of both inducers results in the expression of both trigger RNAs as well as both antisense RNAs, such that the GFP expression is inefficient. While the Lac and Tet operator sequences are also expressed upon the expression of the trigger and antisense RNAs, the expected secondary structure changes on the core signaling parts of the trigger and antisense RNAs were not noticeable (Figure S8).

At the molecular level, the optimized NIMPLY gates previously characterized were used along with the overhang deletion and bulge deletion strategies. The switch RNA that combines two orthogonal switches in an OR-gate fashion is expressed from a low copy plasmid. To facilitate strong repression by the antisense RNAs, both the trigger RNAs were expressed from a medium copy plasmid, and both the antisense RNAs were expressed from a high copy plasmid. The performance of the XOR gate was evaluated in *E. coli* BL21 AI strain, where genomically encoded T7 RNA polymerase was induced by arabinose. A number of basic molecular interactions were verified for the XOR gate: the crosstalk between the switch and antisense RNAs was negligible (Figure S9); both trigger RNAs, despite the attached additional operator sequences, could turn on the switch RNAs (Figure S10); the antisense RNAs could annihilate the cognate trigger RNA activities, as expected, with little crosstalk (Figure S11). When all the components were put together and the chemical inducers were used, the XOR gate functioned as expected, with a high ON state for either IPTG or aTc input but with a low OFF state for no inducer or both inducer cases (Figure 3c). An XOR gate using trigger RNAs with both 50 and 30 overhangs was also shown to be functional, albeit with a reduced ON state for trigger 2 (Figure S12). Furthermore, the GFP output pattern changed sharply as the concentration of inducers was adjusted, indicating that the XOR gate showed a switch-link function suitable for digital circuits (Figure 3d). When incorporated within larger logic circuits, this digital logic ensures an all-or-none response to a variety of inputs and provides a robust output signal regardless of input perturbations [81], thus conveying information with less noise and high accuracy for decision-making processes [82,83].

**Figure 3.** Toehold-switch-based XOR gate. (**a**) Schematics of XOR gate configuration. IPTG and aTc were used as input signals A and B for the XOR gate, and the output signal was visualized through GFP fluorescence. The truth table of the XOR gate indicated the ON and OFF states of the XOR gate depending on the combination of inducer molecules. (**b**) Genetic blueprint of trigger and antisense cassettes and schematics of the XOR gate. Lac operator was placed upstream of T1 and A2, and Tet operator was placed in front of T2 and A1. The optimized extended overhangs were used. (**c**) Performance of toehold-switch-based XOR gate. T7 RNA polymerase was induced with the pretreatment of 1% (*w*/*w*) arabinose in *E. coli* BL21 AI strain. XOR gate components were induced by 0.5 mM IPTG and 100 ng/mL aTc. GFP fluorescence was measured via flow cytometry. Cellular autofluorescence was subtracted in all cases. Autofluorescence level was measured from cells not bearing a GFP-expressing plasmid. Statistical analysis was performed to compare each state of the XOR gate. (Two-tailed Student's *t*-test; \*\*\*\* *p* < 0.0001; Error bars indicate ± SD from three biological replicates) (**d**) Heat map plot of XOR gate output. The color scale was ranged between the minimum and maximum values of the XOR gate output. IPTG and aTc were treated in gradient concentration as described in the table. Each point of the heat map indicates the median value of three replicates.

#### 2.2.1. Cellular Arithmetic Operation of a Half Adder and a Half Subtractor

Building on the RNA devices that were modularized and rigorously characterized earlier, the logical complexity of synthetic RNA circuits can be further increased. As a test case, we focused on basic binary calculators: the half adder and the half subtractor. A half adder takes two input bits and generates two output bits that require an XOR gate for SUM output and an AND gate for CARRY output (Figure 4a). A half subtractor can be

analogously constructed, where an XOR gate computes the DIFFERENCE output and a NIMPLY gate computes the BORROW output. Fortunately, a high-performance AND gate was available from the toehold switch library, and another orthogonal NIMPLY gate was constructed with available toehold switches after NUPACK analysis. The functionality of the AND gate and the NIMPLY gate were verified in isolation (Figure S13). Then, these new genetic elements were incorporated into expression cassettes in the same plasmid backbones as before. The XOR gate with GFP output was used to compute the SUM and DIFFERENCE output bits, and the newly introduced AND gate and NIMPLY gate with mCherry output were used to compute the CARRY and BORROW output bits in the half adder and half subtractor, respectively (Figure 4b). To investigate the performance of the binary calculators at the single-cell level, we characterized the system by flow cytometry. For all input combinations, the half adder and half subtractor showed correct ON and OFF states with statistically significant differences (Figures 4c and S14). Nevertheless, as the genetic complexity and the number of heterologous expression cassettes increased, a concomitant decrease in circuit performance was observed when compared to the single XOR gate. Therefore, we checked all combinations of RNA–RNA interactions with NUPACK <IP_ADDRESS> and confirmed that on-target MFE structures were maintained, albeit with some unintended crosstalk interactions (Table S2). Further improvements in circuit elements, as well as contexts such as promoter arrangements and spacer sequences, may allow for the successful implementation of even more complex circuits such as a full adder.

**Figure 4.** Binary operation of half adder and half subtractor. (**a**) Schematic of half adder and half

subtractor configurations. IPTG and aTc were used as input signals, and the output signals were visualized through GFP and mCherry fluorescence. GFP was assigned to the XOR gate output, and mCherry was assigned to AND or NIMPLY gate output. The truth table of the binary calculators indicated the ON and OFF states of each binary calculator depending on the combination of inducer molecules. Diff. denotes the Difference bit of the half subtractor. (**b**) Schematic of toehold-switchbased half adder and half subtractor. Trigger and antisense RNAs under the same inducer control are shown in boxes. (**c**) Flow cytometry GFP and mCherry fluorescence histograms for the half adder and the half subtractor. The presence of IPTG and aTc was displayed within each panel in brackets. T7 RNA polymerase was induced with the pretreatment of 1% (*w*/*w*) arabinose in *E. coli* BL21 AI strain. RNAs of the half adder and the half subtractor were induced by 1 mM IPTG and 200 ng/mL aTc. Statistical analysis was performed for comparing each state of the binary calculators. (Two-tailed Student's *t*-test; \*\* *p* < 0.01; \*\*\* *p* < 0.001; \*\*\*\* *p* < 0.0001).

#### 2.2.2. Cellular Reversible Logic Operation of Feynman Gate

Reversible computing, a nonconventional form of computing with one-to-one mapping of inputs and outputs, may be useful for biomolecular diagnostic and sensing applications. One of the reversible computing devices, a Feynman gate, can also be constructed using a method similar to other binary calculators, using an XOR gate and a BUFFER gate. Due to the unique output patterns, it is also called as a controlled NOT gate when the output signal (Q) changes from BUFFER gate to NOT gate in response to the input signal (Figure 5a). At the molecular level, we executed the same set of sequences as those for the half subtractor, except that a trigger 3 without extended overhangs was used. The RNA-based Feynman gate possessed both a functional switching ability and the capacity for information storage (Figure 5b and Figure S15). The circuit acted as a BUFFER gate for input B in the absence of input A, whereas the circuit functioned as a NOT gate for input B in the presence of input A. Furthermore, information about the input combinations was preserved in the Feynman gate because of the one-to-one manner of input to output mapping. Overall, we demonstrated that toehold-switch-based ribocomputing designs may prove useful for reversible computing in cells.

**Figure 5.** Demonstration of toehold-switch-based Feynman gate. (**a**) Schematic of the Feynman gate. Truth table of the Feynman gate indicated the ON and OFF states for each of the binary calculators depending on the combination of inducer molecules. Transition denotes the functional transition of a Buffer gate to a NOT gate. IPTG and aTc were used as input signals. Trigger and antisense RNAs under the same inducer control are shown in boxes. (**b**) Performance of Feynman gate. T7 RNA polymerase was induced with the pretreatment of 1% (*w*/*w*) arabinose in *E. coli* BL21 AI strain. RNAs of Feynman gate were induced by 1 mM IPTG and 200 ng/mL aTc. GFP and mCherry fluorescence were measured on flow cytometry. Cellular autofluorescence was subtracted in all cases. Autofluorescence level was measured from cells not bearing a GFP- or mCherry-expressing plasmid. Statistical analysis was performed for comparing each state of the Feynman gate. (Two-tailed Student's *t*-test; \*\*\* *p* < 0.001; \*\*\*\* *p* < 0.0001; Error bars indicate ± SD from three biological replicates).

#### **3. Discussion**

In this study, we present the binary operation of XOR gates, cellular arithmetic operation including a half adder and a half subtractor, and a Feynman gate in single-cell *E. coli* using de-novo-designed toehold switches and antisense RNAs. A systematic approach was taken where the basic building block, a NIMPLY gate, was optimized and then subsequently used for an XOR gate, which in turn could be used for basic ALUs. While the NIMPLY gate design was previously demonstrated [46], it was simply used as a proxy for a NOT gate, and further optimization of its performance was limited. Hence, we investigated several design candidates using a number of parameters, including the ensemble defect, the overhang domain accessibility, and the cross-reactivity in silico. The design candidates with accessible overhang domains generally showed better performance (Figures 2b and S2). Crucially, several optimization strategies can improve the performance of NIMPLY gate designs. A thorough analysis of the impact of extended overhang sequences revealed that a 30 extended overhang on the trigger RNA and a 50 extended overhang on the antisense RNA improved performance compared to having overhangs on the other location or on both sides. We reasoned that the negative effect of a 50 extended overhang on the trigger RNA might be due to the potential interference on the ribosome binding to the RBS of the switch. One piece of evidence in support of this hypothesis is that the trigger with a 50 extended overhang showed improved functionality for toehold switches with increased loop length. Although the mechanistic reasoning on the impact of overhang locations on the antisense RNA is unclear, there still may be physical interference during the initiation stage of trigger and antisense RNA interactions. Recognizing that the single nucleotide bulges located between the overhang and switch binding domains can act as an energy barrier to strand displacement [79,80], we tested the antisense RNA with no bulges and observed improved repression efficiency. These optimization strategies laid the foundation for constructing more complex systems building on the NIMPLY gate designs.

Notably, the antisense RNA designs can be extended to other related synthetic RNA regulators. As an example, a recently reported 3-way junction (3WJ) repressor [47] can be analogously regulated using the antisense RNA design for trigger RNAs (Figure S16). The output characteristics can be considered as an implementation of an IMPLY gate (Figure S17). If applied to the previously reported NOR gate constructed using the 3WJ repressor, an XNOR gate could be constructed similar to the toehold-switch-based XOR gate reported here (Figure S17). Recognizing that NAND gate outputs are distinct from those of XOR gates in the no input case, a NAND gate can be constructed from the current XOR gate by changing the inducible promoters of the trigger RNAs to constitutive promoters (Figure S17). Another important class of de-novo-designed RNA regulator, the small transcription activating RNA (STAR) [41,42], was also subject to antisense RNA-based regulation (Figure S18). Both the T181- and AD1-based STAR designs were successfully regulated using antisense RNA that targets the toehold-binding domain and several bases within the stem-binding domain of the STAR trigger RNAs. Together, these findings indicate that the antisense RNA regulators can be adapted in a straightforward manner to other synthetic RNA regulators and can potentially be used to scale up the complexity of synthetic RNA-based regulatory circuits.

The successful demonstration of a synthetic XOR gate can be seen as a benchmark for systematic synthetic gene circuit construction. Previously, several lines of work demonstrated RNA-based XOR gates, including sRNA [62], miRNA [67], and gRNA [68] that encompass bacterial cells as well as mammalian cell lines. Still, the repression mechanism within the XOR gates relied on protein regulators such as phage-encoded λ repressor protein (CI) [62]. Thus, our demonstration of an RNA-only XOR gate provided a distinct design approach for synthetic XOR gates with performance comparable to the previous work in bacterial cells [62]. More importantly, these RNA-only logic gates can be seamlessly combined for basic ALUs, a half adder and a half subtractor, with performance rivaling previous work [61,62]. Despite thorough in silico analysis and screening for optimized system composition, the performance of basic ALUs showed fold changes less than those of indi-

vidual gates. There are a number of potential pitfalls in the circuit construction, including the leaky expression of promoters, the unexpected interaction between components, and the cellular burden on synthetic RNA expression. Fortunately, these shortcomings can be mitigated with alternative tightly controlled promoters, such as AraBAD or rhaBAD [84–86] or novel synthetic promoters [87], and by the division of load to different cell populations with multicellular networks [88–90].

Herein, we provided a framework for constructing several synthetic RNA-only basic ALUs with de-novo-designed toehold switches at the single-cell level. This design paradigm offers excellent programmability with simple structural specifications rather than sequence constraints. First, the concatenation of switch RNAs can effectively reduce the encoding space of genetic programs and, therefore, enable the operation of complex systems in *E. coli*. Second, the ALUs can be designed with almost no sequence constraints with in silico screening and optimization. Third, a large library of orthogonal toehold switches provides the required parts for building complex systems. Lastly, the system inherits the general advantages of RNA-based operations, including a fast response time, reduced resource usage, and multiplexing [46,91,92]. Recent developments on degradation-tunable RNAs in combination with toehold switches may provide further design flexibility [93]. Notably, a variety of ALU circuits using DNA strand displacement reactions [24] showcases the power of nucleic-acid-based molecular computations. Still, the demonstrations of ALUs in living cells are limited in complexity and scope. The toehold-switch-based ribocomputing circuits could open a new avenue to exploring the rich design space of synthetic RNA-based ALUs, building up to higher-order systems such as a full adder and a full subtractor, ultimately leading to neural-network-like functions in cells.

#### **4. Materials and Methods**

#### *4.1. E. coli Strains and Plasmid Construction*

The following *E. coli* strains were used in this study: BL21 DE3 (Invitrogen; F− *omp*T *hsd*SB (rB − mB −) *gal dcm),* BL21 AI (Invitrogen; F− *omp*T *hsd*SB (rB − mB −) *gal dcm ara*B::T7RNAP-*tet*A), and DH5α (Invitrogen; *endA1 recA1 gyrA96 thi-1 glnV44 relA1 hsdR17*(rK − mK + ) λ −).

The backbones for the plasmids used in this research were taken from the commercial vectors pET15b, pCDFDuet, pCOLADuet, and pACYCDuet (EMD Millipore). The switch RNA of the NIMPLY complex was constructed using ACTS Type II N3 and ACTS Type II N7 from previous research [46] and was constructed in pACYCDuet. All the trigger RNAs and trigger cassettes were constructed in pCDFDuet. All the antisense RNAs and antisense cassettes were constructed in pET15b. The switch RNAs of the AND gate and the NIMPLY gate were constructed in pCOLADuet. All constructs were cloned via blunt end ligation [94], Gibson Assembly [95], circular polymerase extension cloning (CPEC) [96], and/or round-the-horn site-directed mutagenesis [97]. The plasmid architecture and specific part sequences are listed in Tables S3–S11. Plasmids were constructed in *E. coli* DH5α and purified using the EZ-PureTM plasmid Prep Kit. Ver. 2 (Enzynomics). Plasmid sequences were confirmed by DNA sequencing after every cloning step. Plasmids were transformed through chemical transformation [98].

#### *4.2. Cell Culture and Induction Condition*

For in vivo experiments, *E. coli* BL21 DE3 and AI strains were used; they contain chromosomally integrated T7 RNA polymerase under the control of IPTG-inducible lacUV5 promoter and arabinose-inducible PBAD promoter, respectively. For the in vivo characterization of the NIMPLY complex in Figure 2, chemically transformed *E. coli* BL21 DE3 cells were cultured on LB agar plates (1.5% agar) with appropriate antibiotics: pACYC-Duet (25 µg/mL Chloramphenicol), pCDFDuet (50 µg/mL Spectinomycin), and pET15b (100 µg/mL Ampicillin). Single colonies were grown overnight (~16 h) in 96-well plates with shaking at 800 rpm, 37 ◦C. Overnight cultures were diluted 1/100-fold into fresh media and returned to shaking (800 rpm, 37 ◦C). After 80 min, cell cultures were induced

with 1 mM IPTG (Promega) and returned to the shaker (800 rpm, 37 ◦C) until fluorescence measurement after 3 h 30 min. For the experiments on the toehold-switch-based XOR gates, a half adder, a half subtractor, and a Feynman gate, chemically transformed *E. coli* BL21 AI cells were cultured on LB agar plates (BD biosciences; 1.5% agar) with appropriate antibiotics. All antibiotics were purchased from Gold biotechnology: pACYCDuet (25 µg/mL Chloramphenicol), pCOLADuet (50 µg/mL Kanamycin), pCDFDuet (50 µg/mL Spectinomycin), and pET15b (100 µg/mL Ampicillin). Single colonies were grown overnight (~16 h) in 96-well plates with shaking at 800 rpm, 37 ◦C. Overnight cultures were diluted 1/100-fold into fresh media and returned to shaking (800 rpm, 37 ◦C). After 80 min, cell cultures were induced with 1% (*w*/*w*) arabinose (Gold biotechnology) to produce T7 RNA polymerase for 1 h with shaking (800 rpm, 37 ◦C). Then the cell cultures were induced with 0.5 mM IPTG, 100 ng/mL aTc (Takara) for the XOR gates and 1 mM IPTG, 200 ng/mL aTc for a half adder, a half subtractor, and a Feynman gate, and returned to the shaker (800 rpm, 37 ◦C) until fluorescence measurement after 3 h 30 min.

#### *4.3. Microplate Reader Analysis*

For the experiment on the XOR gate with a gradient concentration of chemical inducers (Figure 3d), 200 µL of cell cultures were added per well on a 96-well Black Plate 33,396 (SPL) after 1 mM IPTG induction. GFP fluorescence (excitation: 479 nm, emission: 520 nm), mCherry fluorescence (excitation: 587 nm, emission: 610 nm), and OD600 were measured in a Synergy H1 microplate reader (Biotek) running Gen5 3.08 software. GFP and mCherry fluorescence levels were normalized as follows: the fluorescence of LB blank was subtracted for background normalization, and a measured fluorescence value was divided by its OD600. The number of biological replicates was three for in vivo experiments.

#### *4.4. Fluorescence Measurements Using Flow Cytometry*